LORENE
tslice_conf_init.C
1 /*
2  * Method of class Time_slice_conf to compute valid initial data
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_conf_init.C,v 1.14 2016/12/05 16:18:19 j_novak Exp $
32  * $Log: tslice_conf_init.C,v $
33  * Revision 1.14 2016/12/05 16:18:19 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.13 2014/10/13 08:53:48 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.12 2014/10/06 15:13:22 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.11 2010/10/20 07:58:09 j_novak
43  * Better implementation of the explicit time-integration. Not fully-tested yet.
44  *
45  * Revision 1.10 2008/12/04 18:22:49 j_novak
46  * Enhancement of the dzpuis treatment + various bug fixes.
47  *
48  * Revision 1.9 2008/12/02 15:02:22 j_novak
49  * Implementation of the new constrained formalism, following Cordero et al. 2009
50  * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
51  *
52  * Revision 1.8 2004/05/17 19:53:13 e_gourgoulhon
53  * Added arguments graph_device and method_poisson_vect.
54  *
55  * Revision 1.7 2004/05/12 15:24:20 e_gourgoulhon
56  * Reorganized the #include 's, taking into account that
57  * time_slice.h contains now an #include "metric.h".
58  *
59  * Revision 1.6 2004/05/10 09:12:01 e_gourgoulhon
60  * Added a call to del_deriv() at the end.
61  *
62  * Revision 1.5 2004/05/03 14:48:48 e_gourgoulhon
63  * Treatment of special cases nn_jp1.etat = ETATUN and psi_jp1.etat = ETATUN.
64  *
65  * Revision 1.4 2004/04/29 17:10:36 e_gourgoulhon
66  * Added argument pdt and update of depth slices at the end,
67  * taking into account the known time derivatives.
68  *
69  * Revision 1.3 2004/04/08 16:45:11 e_gourgoulhon
70  * Use of new methods set_*.
71  *
72  * Revision 1.2 2004/04/07 07:58:21 e_gourgoulhon
73  * Constructor as Minkowski slice: added call to std_spectral_base()
74  * after setting the lapse to 1.
75  *
76  * Revision 1.1 2004/04/05 21:25:37 e_gourgoulhon
77  * First version.
78  *
79  *
80  * $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_conf_init.C,v 1.14 2016/12/05 16:18:19 j_novak Exp $
81  *
82  */
83 
84 // C headers
85 #include <cassert>
86 
87 // Lorene headers
88 #include "time_slice.h"
89 #include "unites.h"
90 #include "graphique.h"
91 #include "utilitaires.h"
92 
93 namespace Lorene {
95  const Scalar& trk_in, const Scalar& trk_point,
96  double pdt, double precis, int method_poisson_vect,
97  const char* graph_device, const Scalar* p_ener_dens,
98  const Vector* p_mom_dens, const Scalar* p_trace_stress) {
99 
100  using namespace Unites ;
101 
102  // Verifications
103  // -------------
104  double tr_uu = max(maxabs(uu.trace(tgam()), "trace tgam_{ij} u^{ij}")) ;
105  if (tr_uu > 1.e-7) {
106  cerr <<
107  "Time_slice_conf::initial_data_cts : the trace of u^{ij} with respect\n"
108  << " to the conformal metric tgam_{ij} is not zero !\n"
109  << " error = " << tr_uu << endl ;
110  abort() ;
111  }
112 
113  assert(trk_in.check_dzpuis(2)) ;
114  assert(trk_point.check_dzpuis(4)) ;
115 
116  // Initialisations
117  // ---------------
118  double ttime = the_time[jtime] ;
119 
120  trk_evol.update(trk_in, jtime, ttime) ;
121 
122  // Reset of quantities depending on K:
123  k_dd_evol.downdate(jtime) ;
124  k_uu_evol.downdate(jtime) ;
125 
126  set_hata(psi4()*psi()*psi()* uu / (2.* nn()) ) ;
127 
128  const Map& map = uu.get_mp() ;
129  const Base_vect& triad = *(uu.get_triad()) ;
130 
131  // For graphical outputs:
132  int ngraph0 = 10 ; // index of the first graphic device to be used
133  int nz = map.get_mg()->get_nzone() ;
134  double ray_des = 1.25 * map.val_r(nz-2, 1., 0., 0.) ; // outermost radius
135  // for plots
136 
137  Scalar ener_dens(map) ;
138  if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
139  else ener_dens.set_etat_zero() ;
140 
141  Vector mom_dens(map, CON, triad) ;
142  if (p_mom_dens != 0x0) mom_dens = *(p_mom_dens) ;
143  else mom_dens.set_etat_zero() ;
144 
145  Scalar trace_stress(map) ;
146  if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
147  else trace_stress.set_etat_zero() ;
148 
149  Scalar tmp(map) ;
150  Scalar source_psi(map) ;
151  Scalar source_nn(map) ;
152  Vector source_beta(map, CON, triad) ;
153 
154  // Iteration
155  // ---------
156  int imax = 100 ;
157  for (int i=0; i<imax; i++) {
158 
159  //===============================================
160  // Computations of sources
161  //===============================================
162 
163  const Vector& dpsi = psi().derive_cov(ff) ; // D_i Psi
164  const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
165  const Vector& dnn = nn().derive_cov(ff) ; // D_i N
166 
167  Sym_tensor taa = aa().up_down(tgam()) ;
168  Scalar aa_quad = contract(taa, 0, 1, aa(), 0, 1) ;
169 
170  // Source for Psi
171  // --------------
172  tmp = 0.125* psi() * tgam().ricci_scal()
173  - contract(hh(), 0, 1, dpsi.derive_cov(ff), 0, 1 ) ;
174  tmp.inc_dzpuis() ; // dzpuis : 3 -> 4
175 
176  tmp -= contract(hdirac(), 0, dpsi, 0) ;
177 
178  source_psi = tmp - psi()*psi4()* ( 0.5*qpig* ener_dens
179  + 0.125* aa_quad
180  - 8.33333333333333e-2* trk()*trk() ) ;
181 
182  // Source for N
183  // ------------
184 
185  source_nn = psi4()*( nn()*( qpig* (ener_dens + trace_stress) + aa_quad
186  - 0.3333333333333333* trk()*trk() )
187  - trk_point )
188  - 2.* contract(dln_psi, 0, nn().derive_con(tgam()), 0)
189  - contract(hdirac(), 0, dnn, 0) ;
190 
191  tmp = psi4()* contract(beta(), 0, trk().derive_cov(ff), 0)
192  - contract( hh(), 0, 1, dnn.derive_cov(ff), 0, 1 ) ;
193 
194  tmp.inc_dzpuis() ; // dzpuis: 3 -> 4
195 
196  source_nn += tmp ;
197 
198 
199  // Source for beta
200  // ---------------
201 
202  source_beta = 2.* contract(aa(), 1,
203  dnn - 6.*nn() * dln_psi, 0) ;
204 
205  source_beta += 2.* nn() * ( 2.*qpig* psi4() * mom_dens
206  + 0.66666666666666666* trk().derive_con(tgam())
207  - contract(tgam().connect().get_delta(), 1, 2,
208  aa(), 0, 1) ) ;
209 
210  Vector vtmp = contract(hh(), 0, 1,
211  beta().derive_cov(ff).derive_cov(ff), 1, 2)
212  + 0.3333333333333333*
213  contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0)
214  - hdirac().derive_lie(beta())
215  + uu.divergence(ff) ;
216  vtmp.inc_dzpuis() ; // dzpuis: 3 -> 4
217 
218  source_beta -= vtmp ;
219 
220  source_beta += 0.66666666666666666* beta().divergence(ff) * hdirac() ;
221 
222 
223  //=============================================
224  // Resolution of elliptic equations
225  //=============================================
226 
227  // Resolution of the Poisson equation for Psi
228  // ------------------------------------------
229 
230  Scalar psi_jp1 = source_psi.poisson() + 1. ;
231 
232  if (psi_jp1.get_etat() == ETATUN) psi_jp1.std_spectral_base() ;
233 
234  // Test:
235  maxabs(psi_jp1.laplacian() - source_psi,
236  "Absolute error in the resolution of the equation for Psi") ;
237 
238  des_meridian(psi_jp1, 0., ray_des, "Psi", ngraph0, graph_device) ;
239 
240  // Resolution of the Poisson equation for the lapse
241  // ------------------------------------------------
242 
243  Scalar nn_jp1 = source_nn.poisson() + 1. ;
244 
245  if (nn_jp1.get_etat() == ETATUN) nn_jp1.std_spectral_base() ;
246 
247  // Test:
248  maxabs(nn_jp1.laplacian() - source_nn,
249  "Absolute error in the resolution of the equation for N") ;
250 
251  des_meridian(nn_jp1, 0., ray_des, "N", ngraph0+1, graph_device) ;
252 
253  // Resolution of the vector Poisson equation for the shift
254  //---------------------------------------------------------
255 
256  Vector beta_jp1 = source_beta.poisson(0.3333333333333333, ff,
257  method_poisson_vect) ;
258 
259  des_meridian(beta_jp1(1), 0., ray_des, "\\gb\\ur\\d", ngraph0+2,
260  graph_device) ;
261  des_meridian(beta_jp1(2), 0., ray_des, "\\gb\\u\\gh\\d", ngraph0+3,
262  graph_device) ;
263  des_meridian(beta_jp1(3), 0., ray_des, "\\gb\\u\\gf\\d", ngraph0+4,
264  graph_device) ;
265 
266  // Test:
267  Vector test_beta = (beta_jp1.derive_con(ff)).divergence(ff)
268  + 0.3333333333333333 * (beta_jp1.divergence(ff)).derive_con(ff) ;
269  test_beta.inc_dzpuis() ;
270  maxabs(test_beta - source_beta,
271  "Absolute error in the resolution for beta") ;
272 
273  //===========================================
274  // Convergence control
275  //===========================================
276 
277  double diff_psi = max( diffrel(psi(), psi_jp1) ) ;
278  double diff_nn = max( diffrel(nn(), nn_jp1) ) ;
279  double diff_beta = max( diffrel(beta(), beta_jp1) ) ;
280 
281  cout << "step = " << i << " : diff_psi = " << diff_psi
282  << " diff_nn = " << diff_nn
283  << " diff_beta = " << diff_beta << endl ;
284  if ( (diff_psi < precis) && (diff_nn < precis) && (diff_beta < precis) )
285  break ;
286 
287  //=============================================
288  // Updates for next step
289  //=============================================
290 
291  set_psi_del_npsi(psi_jp1) ;
292 
293  n_evol.update(nn_jp1, jtime, ttime) ;
294 
295  beta_evol.update(beta_jp1, jtime, ttime) ;
296 
297  // New value of A^{ij}:
298  Sym_tensor aa_jp1 = ( beta().ope_killing_conf(tgam()) + uu )
299  / (2.* nn()) ;
300 
301  set_hata( aa_jp1 / (psi4()*psi()*psi()) ) ;
302 
303  }
304 
305  //==================================================================
306  // End of iteration
307  //===================================================================
308 
309  npsi_evol.update( n_evol[jtime]*psi_evol[jtime], jtime, ttime ) ;
310  A_hata() ;
311  B_hata() ;
312 
313  // Push forward in time to enable the computation of time derivatives
314  // ------------------------------------------------------------------
315 
316  double ttime1 = ttime ;
317  int jtime1 = jtime ;
318  for (int j=1; j < depth; j++) {
319  jtime1++ ;
320  ttime1 += pdt ;
321  psi_evol.update(psi_evol[jtime], jtime1, ttime1) ;
322  npsi_evol.update(npsi_evol[jtime], jtime1, ttime1) ;
323  n_evol.update(n_evol[jtime], jtime1, ttime1) ;
324  beta_evol.update(beta_evol[jtime], jtime1, ttime1) ;
325  hh_evol.update(hh_evol[jtime], jtime1, ttime1) ;
326  hata_evol.update(hata_evol[jtime], jtime1, ttime1) ;
327  A_hata_evol.update(A_hata_evol[jtime], jtime1, ttime1) ;
328  B_hata_evol.update(B_hata_evol[jtime], jtime1, ttime1) ;
329  trk_evol.update(trk_evol[jtime], jtime1, ttime1) ;
330  the_time.update(ttime1, jtime1, ttime1) ;
331  }
332  jtime += depth - 1 ;
333 
334  // Taking into account the time derivative of h^{ij} and K :
335  // ---------------------------------------------------------
336  Sym_tensor uu0 = uu ;
337  uu0.dec_dzpuis(2) ; // dzpuis: 2 --> 0
338 
339  for (int j=1; j < depth; j++) {
340  hh_evol.update(hh_evol[jtime] - j*pdt* uu0,
341  jtime-j, the_time[jtime-j]) ;
342 
343  trk_evol.update(trk_evol[jtime] - j*pdt* trk_point,
344  jtime-j, the_time[jtime-j]) ;
345 
346  }
347 
348  // Reset of derived quantities (at the new time step jtime)
349  // ---------------------------
350  del_deriv() ;
351 
352 }
353 }
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_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition: time_slice.h:525
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
Definition: time_slice.h:520
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 ).
virtual void del_deriv() const
Deletes all the derived quantities.
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
virtual void initial_data_cts(const Sym_tensor &uu, const Scalar &trk_in, const Scalar &trk_point, double pdt, double precis=1.e-12, int method_poisson_vect=6, const char *graph_device=0x0, const Scalar *ener_dens=0x0, const Vector *mom_dens=0x0, const Scalar *trace_stress=0x0)
Computes valid initial data by solving the constraint equations in the conformal thin-sandwich approa...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
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
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition: time_slice.h:550
Base class for coordinate mappings.
Definition: map.h:688
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
int jtime
Time step index of the latest slice.
Definition: time_slice.h:193
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
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition: time_slice.h:555
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:817
Evolution_std< Sym_tensor > k_uu_evol
Values at successive time steps of the contravariant components of the extrinsic curvature tensor ...
Definition: time_slice.h:216
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: ...
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
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...
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:353
virtual const Scalar & A_hata() const
Returns the potential A of .
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Evolution_std< Sym_tensor > k_dd_evol
Values at successive time steps of the covariant components of the extrinsic curvature tensor ...
Definition: time_slice.h:211
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
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
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:196
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
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and ...
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
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: vector.C:398
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition: time_slice.h:227
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition: time_slice.h:219
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
virtual const Scalar & B_hata() const
Returns the potential of .
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition: time_slice.h:545
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:879
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition: time_slice.h:222
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
int depth
Number of stored time slices.
Definition: time_slice.h:182
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition: time_slice.h:533
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )