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
28char 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
90namespace Lorene {
91void 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:
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
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}
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
void downdate(int j)
Suppresses a stored value.
Definition evolution.C:244
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:530
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual const Scalar & nn() const
Lapse function N 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: .
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition time_slice.h:542
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
virtual const Scalar & A_hata() const
Returns the potential A of .
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual const Scalar & B_hata() const
Returns the potential of .
virtual void del_deriv() const
Deletes all the derived quantities.
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition time_slice.h:522
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: .
const Scalar & psi4() const
Factor at the current time step (jtime ).
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
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition time_slice.h:547
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition time_slice.h:507
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition time_slice.h:552
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition time_slice.h:224
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
virtual const Vector & beta() const
shift vector at the current time step (jtime )
int depth
Number of stored time slices.
Definition time_slice.h:179
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition time_slice.h:216
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
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
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:219
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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 .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
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.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.