LORENE
tslice_dirac_max_solve.C
1 /*
2  * Methods of class Tslice_dirac_max for solving Einstein equations
3  *
4  * (see file time_slice.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Eric Gourgoulhon & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: tslice_dirac_max_solve.C,v 1.19 2016/12/05 16:18:19 j_novak Exp $
32  * $Log: tslice_dirac_max_solve.C,v $
33  * Revision 1.19 2016/12/05 16:18:19 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.18 2014/10/13 08:53:48 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.17 2014/10/06 15:13:22 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.16 2010/10/20 07:58:10 j_novak
43  * Better implementation of the explicit time-integration. Not fully-tested yet.
44  *
45  * Revision 1.15 2008/12/02 15:02:22 j_novak
46  * Implementation of the new constrained formalism, following Cordero et al. 2009
47  * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
48  *
49  * Revision 1.14 2004/07/08 12:29:01 j_novak
50  * use of new method Tensor::annule_extern_cn
51  *
52  * Revision 1.13 2004/06/30 08:02:40 j_novak
53  * Added filtering in l of khi_new and mu_new. ki_source is forced to go to
54  * zero at least as r^2.
55  *
56  * Revision 1.12 2004/06/17 07:07:11 e_gourgoulhon
57  * Method solve_hij:
58  * -- replaced the attenuation of khi_source with tempo by a call to
59  * the new method Tensor::annule_extern_c2
60  * -- suppressed filtre_r on khi_source and khi_new
61  * -- added graphs of W^i and LW^{ij} (transverse decomp of S^ij).
62  *
63  * Revision 1.11 2004/06/15 09:43:36 j_novak
64  * Attenuation of the source for khi in the last shell (temporary?).
65  *
66  * Revision 1.10 2004/06/14 20:47:31 e_gourgoulhon
67  * Added argument method_poisson to method solve_hij.
68  *
69  * Revision 1.9 2004/06/03 10:02:45 j_novak
70  * Some filtering is done on source_khi and khi_new.
71  *
72  * Revision 1.8 2004/05/24 21:00:44 e_gourgoulhon
73  * Method solve_hij: khi and mu.smooth_decay(2,2) --> smooth_decay(2,1) ;
74  * added exponential_decay(khi) and exponential_decay(mu) after the
75  * call to smooth_decay. Method exponential_decay is provisory defined
76  * in this file.
77  *
78  * Revision 1.7 2004/05/17 19:56:25 e_gourgoulhon
79  * -- Method solve_beta: added argument method
80  * -- Method solve_hij: added argument graph_device
81  *
82  * Revision 1.6 2004/05/12 15:24:20 e_gourgoulhon
83  * Reorganized the #include 's, taking into account that
84  * time_slice.h contains now an #include "metric.h".
85  *
86  * Revision 1.5 2004/05/05 14:47:05 e_gourgoulhon
87  * Modified text and graphical outputs.
88  *
89  * Revision 1.4 2004/05/03 15:06:27 e_gourgoulhon
90  * Added matter source in solve_hij.
91  *
92  * Revision 1.3 2004/05/03 14:50:00 e_gourgoulhon
93  * Finished the implementation of method solve_hij.
94  *
95  * Revision 1.2 2004/04/30 14:36:15 j_novak
96  * Added the method Tslice_dirac_max::solve_hij(...)
97  * NOT READY YET!!!
98  *
99  * Revision 1.1 2004/04/30 10:52:14 e_gourgoulhon
100  * First version.
101  *
102  *
103  * $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_solve.C,v 1.19 2016/12/05 16:18:19 j_novak Exp $
104  *
105  */
106 
107 // C headers
108 #include <cassert>
109 
110 // Lorene headers
111 #include "time_slice.h"
112 #include "unites.h"
113 #include "graphique.h"
114 #include "proto.h"
115 
116  //----------------------------//
117  // Equation for N\Psi //
118  //----------------------------//
119 
120 namespace Lorene {
122  const Scalar* p_trace_stress) const {
123 
124  using namespace Unites ;
125 
126  const Map& map = npsi().get_mp() ;
127  Scalar ener_dens(map) ;
128  if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
129  else ener_dens.set_etat_zero() ;
130 
131  Scalar trace_stress(map) ;
132  if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
133  else trace_stress.set_etat_zero() ;
134 
135  // Source for N\Psi
136  // ----------------
137 
138  const Vector& dnpsi = npsi().derive_cov(ff) ; // D_i N\Psi
139  const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
140  const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
141  // D_k {\tilde \gamma}_{ij}
142  Sym_tensor taa = aa().up_down(tgam()) ;
143 
144  Scalar aa_quad = contract(taa, 0, 1, aa(), 0, 1) ;
145  Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
146  - 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
147 
148  Scalar source_npsi = - contract( hh(), 0, 1, dnpsi.derive_cov(ff), 0, 1 ) ;
149  source_npsi.inc_dzpuis() ;
150  source_npsi += npsi()*( 0.5*qpig*(ener_dens + 2.*trace_stress )/( psi()*psi() )
151  + 0.875*aa_quad/( psi4()*psi4() ) + 0.125*tildeR ) ;
152 
153  // Resolution of the Poisson equation for N\Psi
154  // --------------------------------------------
155 
156  Scalar npsi_new = source_npsi.poisson() + 1. ;
157 
158  if (npsi_new.get_etat() == ETATUN) npsi_new.std_spectral_base() ;
159 
160 #ifndef NDEBUG
161  // Test:
162  maxabs(npsi_new.laplacian() - source_npsi,
163  "Absolute error in the resolution of the equation for N") ;
164 #endif
165  return npsi_new ;
166 
167 }
168 
169 
170  //--------------------------//
171  // Equation for \Psi //
172  //--------------------------//
173 
174 
175 Scalar Tslice_dirac_max::solve_psi(const Scalar* p_ener_dens) const {
176 
177  using namespace Unites ;
178 
179  const Map& map = psi().get_mp() ;
180  Scalar ener_dens(map) ;
181  if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
182  else ener_dens.set_etat_zero() ;
183 
184  // Source for \Psi
185  // ---------------
186 
187  const Vector& dpsi = psi().derive_cov(ff) ; // D_i Psi
188  const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
189  const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
190  // D_k {\tilde \gamma}_{ij}
191 
192  Sym_tensor taa = hata().up_down(tgam()) ;
193 
194  Scalar aa_quad = contract(taa, 0, 1, hata(), 0, 1) ;
195 
196  Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
197  - 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
198 
199  Scalar source_psi = -contract( hh(), 0, 1, dpsi.derive_cov(ff), 0, 1 ) ;
200  source_psi.inc_dzpuis() ;
201  source_psi -= 0.5*qpig*ener_dens/psi()
202  + 0.125*( aa_quad*pow(psi(), -7) - tildeR*psi() ) ;
203 
204  // Resolution of the Poisson equation for Psi
205  // ------------------------------------------
206 
207  Scalar psi_new = source_psi.poisson() + 1. ;
208 
209  if (psi_new.get_etat() == ETATUN) psi_new.std_spectral_base() ;
210 
211 #ifndef NDEBUG
212  // Test:
213  maxabs(psi_new.laplacian() - source_psi,
214  "Absolute error in the resolution of the equation for Psi") ;
215 #endif
216 
217  return psi_new ;
218 }
219 
220 
221 
222  //--------------------------//
223  // Equation for beta //
224  //--------------------------//
225 
226 
228  const {
229 
230  // Source for beta
231  // ---------------
232  Vector source_beta =
233  - contract(hh(), 0, 1, beta().derive_cov(ff).derive_cov(ff), 1, 2)
234  - 0.3333333333333333*contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0) ;
235 
236  Sym_tensor sym_tmp = 2*nn()*aa() ;
237  source_beta += sym_tmp.divergence(ff) ;
238  source_beta.inc_dzpuis() ; // dzpuis: 3 -> 4
239 
240  // Resolution of the vector Poisson equation
241  //------------------------------------------
242 
243  Vector beta_new = source_beta.poisson(0.3333333333333333, ff, method) ;
244 
245 #ifndef NDEBUG
246  // Test:
247  Vector test_beta = (beta_new.derive_con(ff)).divergence(ff)
248  + 0.3333333333333333 * (beta_new.divergence(ff)).derive_con(ff) ;
249  test_beta.inc_dzpuis() ;
250  maxabs(test_beta - source_beta,
251  "Absolute error in the resolution of the equation for beta") ;
252 #endif
253  return beta_new ;
254 
255 }
256 
257 
258 }
virtual const Vector & beta() const
shift vector at the current time step (jtime )
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
const Scalar & psi4() const
Factor at the current time step (jtime ).
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
Base class for coordinate mappings.
Definition: map.h:688
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Tensor field of valence 1.
Definition: vector.h:188
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $$ (Hamiltonian constraint).
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
virtual Scalar solve_npsi(const Scalar *ener_dens=0x0, const Scalar *trace_stress=0x0) const
Solves the elliptic equation for (maximal slicing condition + Hamiltonian constraint) ...
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:825
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1023
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )