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