LORENE
tslice_dirac_max_evolve.C
1/*
2 * Method of class Tslice_dirac_max for time evolution
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_evolve.C,v 1.23 2016/12/05 16:18:19 j_novak Exp $
32 * $Log: tslice_dirac_max_evolve.C,v $
33 * Revision 1.23 2016/12/05 16:18:19 j_novak
34 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35 *
36 * Revision 1.22 2015/08/10 15:32:27 j_novak
37 * Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
38 *
39 * Revision 1.21 2014/10/13 08:53:48 j_novak
40 * Lorene classes and functions now belong to the namespace Lorene.
41 *
42 * Revision 1.20 2013/01/24 12:55:18 j_novak
43 * Corrected the declaration of variables for boundary conditions.
44 *
45 * Revision 1.19 2012/02/06 12:59:07 j_novak
46 * Correction of some errors.
47 *
48 * Revision 1.18 2011/07/22 13:21:02 j_novak
49 * Corrected an error on BC treatment.
50 *
51 * Revision 1.17 2010/10/20 07:58:10 j_novak
52 * Better implementation of the explicit time-integration. Not fully-tested yet.
53 *
54 * Revision 1.16 2008/12/04 18:22:49 j_novak
55 * Enhancement of the dzpuis treatment + various bug fixes.
56 *
57 * Revision 1.15 2008/12/02 15:02:22 j_novak
58 * Implementation of the new constrained formalism, following Cordero et al. 2009
59 * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
60 *
61 * Revision 1.13 2004/06/14 20:51:37 e_gourgoulhon
62 * Method solve_hij has now argument method_poisson.
63 * Its value is set **provisory** to 1 (instead of method_poisson_vect !).
64 *
65 * Revision 1.12 2004/05/31 09:09:59 e_gourgoulhon
66 * Added monitoring of khi and mu.
67 * Added writing of whole configuration in file (via Time_slice::save).
68 *
69 * Revision 1.11 2004/05/24 20:58:05 e_gourgoulhon
70 * Added graphical output of khi, mu and trh.
71 *
72 * Revision 1.10 2004/05/20 20:32:01 e_gourgoulhon
73 * Added arguments check_mod and save_mod.
74 * Argument graph_device passed to des_evol.
75 *
76 * Revision 1.9 2004/05/17 19:55:10 e_gourgoulhon
77 * Added arguments method_poisson_vect, nopause and graph_device
78 *
79 * Revision 1.8 2004/05/13 21:35:30 e_gourgoulhon
80 * Added monitoring of various quantities (as Evolution_full<Tbl>).
81 * Added function monitor_scalar.
82 *
83 * Revision 1.7 2004/05/12 15:24:20 e_gourgoulhon
84 * Reorganized the #include 's, taking into account that
85 * time_slice.h contains now an #include "metric.h".
86 *
87 * Revision 1.6 2004/05/11 20:15:10 e_gourgoulhon
88 * Added Evolution_full's for ADM mass and checks of the constraint,
89 * as well as the corresponding plots and write to files.
90 *
91 * Revision 1.5 2004/05/10 09:19:27 e_gourgoulhon
92 * Added a call to del_deriv() after set_khi_mu.
93 *
94 * Revision 1.4 2004/05/09 20:59:06 e_gourgoulhon
95 * Change of the time scheme: first solve d'Alembert equations,
96 * then psuh forward in time and solve the elliptic equation
97 * on the new slice.
98 *
99 * Revision 1.3 2004/05/06 15:26:29 e_gourgoulhon
100 * No longer necessary to initialize khi and mu.
101 *
102 * Revision 1.2 2004/05/05 14:39:32 e_gourgoulhon
103 * Added graphical outputs.
104 *
105 * Revision 1.1 2004/05/03 14:49:10 e_gourgoulhon
106 * First version
107 *
108 *
109 * $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_evolve.C,v 1.23 2016/12/05 16:18:19 j_novak Exp $
110 *
111 */
112
113// Lorene headers
114#include "time_slice.h"
115#include "param.h"
116#include "graphique.h"
117#include "utilitaires.h"
118#include "proto.h"
119
120namespace Lorene {
121const Tbl& monitor_scalar(const Scalar& uu, Tbl& resu) ;
122
123void Tslice_dirac_max::evolve(double pdt, int nb_time_steps,
124 int niter_elliptic, double relax,
125 int check_mod, int save_mod,
126 int method_poisson_vect, int nopause,
127 const char* graph_device, bool verbose,
128 const Scalar* ener_euler,
129 const Vector* mom_euler, const Scalar* s_euler,
130 const Sym_tensor* strain_euler) {
131
132 // Intermediate quantities
133 // -----------------------
134 const Map& map = nn().get_mp() ;
135 const Base_vect& triad = *(beta().get_triad()) ;
136 assert( triad == map.get_bvect_spher() ) ;
137
138 // For graphical outputs:
139 int ngraph0 = 20 ; // index of the first graphic device to be used
140 int ngraph0_mon = 70 ; // for monitoring global quantities
141 int nz = map.get_mg()->get_nzone() ;
142 int nz_bound = nz - 2 ;
143 int nt = map.get_mg()->get_nt(0) ;
144 int np = map.get_mg()->get_np(0) ;
145 int np2 = np+2 ;
146 Scalar tmp(map) ; tmp.std_spectral_base() ;
147 Base_val base_ref = tmp.get_spectral_base() ;
148 Base_val base_pseudo = base_ref ;
149 base_pseudo.mult_cost() ;
150 base_pseudo.mult_x() ;
151#ifndef NDEBUG
152 for (int lz=1; lz<nz; lz++) {
153 assert( map.get_mg()->get_np(lz) == np ) ;
154 assert( map.get_mg()->get_nt(lz) == nt ) ;
155 }
156 assert (depth > 2) ; //## for the moment, a more flexible test should be put
157 for (int it=0; it<depth; it++) {
158 assert ( hh_evol.is_known( jtime - it ) ) ;
159 assert ( hata_evol.is_known( jtime - it ) ) ;
160 assert ( A_hata_evol.is_known( jtime - it ) ) ;
161 assert ( B_hata_evol.is_known( jtime - it ) ) ;
162 assert ( A_hh_evol.is_known( jtime - it ) ) ;
163 assert ( B_hh_evol.is_known( jtime - it ) ) ;
164 }
165#endif
166
167 // Initialization of the TT fields
168 //--------------------------------
169 Sym_tensor_tt hij_tt( map, triad, ff ) ;
170 hij_tt.set_A_tildeB( A_hh_evol[jtime], B_hh_evol[jtime] ) ;
171 Sym_tensor_tt hijtt_old = hij_tt ;
172 for (int i=1; i<=3; i++)
173 for (int j=i; j<=3; j++)
174 if ( hijtt_old(i,j).get_etat() == ETATZERO )
175 hijtt_old.set( i, j ).annule_hard() ;
176 hijtt_old.annule(0, nz_bound) ;
177
178 Sym_tensor_tt hata_tt( map, triad, ff ) ;
179 hata_tt.set_A_tildeB( A_hata_evol[jtime], B_hata_evol[jtime] ) ;
180 hata_tt.inc_dzpuis(2) ;
181 Sym_tensor_tt hatatt_old = hata_tt ;
182 for (int i=1; i<=3; i++)
183 for (int j=i; j<=3; j++)
184 if ( hatatt_old(i,j).get_etat() == ETATZERO )
185 hatatt_old.set( i, j ).annule_hard() ;
186 hatatt_old.annule(0, nz_bound) ;
187
188 // Declaration / initialization of mu and khi for hh and hata
189 //-----------------------------------------------------------
190 Evolution_std<Scalar> khi_hh_evol(depth) ;
191 Evolution_std<Scalar> mu_hh_evol(depth) ;
192 Evolution_std<Scalar> khi_a_evol(depth) ;
193 Evolution_std<Scalar> mu_a_evol(depth) ;
194 Sym_tensor_trans Tij(map, map.get_bvect_spher(), ff) ;
195 for (int j=jtime-depth+1; j<=jtime; j++) {
196 tmp = hij_tt(1,1) ;
197 tmp.mult_r() ; tmp.mult_r() ;
198 khi_hh_evol.update( tmp, j, the_time[j] ) ;
199 mu_hh_evol.update( hij_tt.mu(), j, the_time[j] ) ;
200 tmp = hata_tt(1,1) ;
201 tmp.mult_r() ; tmp.mult_r() ;
202 khi_a_evol.update( tmp, j, the_time[j] ) ;
203 mu_a_evol.update( hata_tt.mu(), j, the_time[j] ) ;
204 }
205
206 double Rmax = map.val_r(nz-2, 1., 0., 0.) ; // outermost radius
207 double ray_des = 1.25 * Rmax ; // for plots
208
209 // Parameters for the evolution equations
210 //---------------------------------------
211 double an = 23./12. ;
212 double anm1 = -4./3. ;
213 double anm2 = 5./12. ;
214
215 int i_zero = 0 ;
216 int i_minus_one = -1 ;
217 int i_two = 2 ;
218
219 Param par_A ;
220 double *l_val_A = new double(1./Rmax) ;
221 double *l_der_A = new double(1.) ;
222 par_A.add_int(nz_bound, 0) ;
223 par_A.add_int(i_two, 1) ; //matching of function and derivative
224 par_A.add_int(i_zero, 2) ;// no shift in l
225 par_A.add_int(i_two, 3) ; // only for l>=2
226 par_A.add_double_mod(*l_val_A, 0) ;
227 par_A.add_double_mod(*l_der_A, 1) ;
228 Tbl* tmp_Acc = new Tbl(np2, nt) ;
229 Tbl& Acc = *tmp_Acc ;
230 Acc.annule_hard() ;
231 par_A.add_tbl_mod(Acc) ;
232 Param par_mat_A_hh ;
233
234 Param par_B ;
235 double* l_val_B = new double(1./Rmax) ;
236 double* l_der_B = new double(1.) ;
237 par_B.add_int(nz_bound, 0) ;
238 par_B.add_int(i_two, 1) ; //matching of function and derivative
239 par_B.add_int(i_minus_one, 2) ;// shift in l for tilde{B}
240 par_B.add_int(i_two, 3) ; // only for l>=2
241 par_B.add_double_mod(*l_val_B, 0) ;
242 par_B.add_double_mod(*l_der_B, 1) ;
243 Tbl* tmp_Bcc = new Tbl(np2, nt) ;
244 Tbl& Bcc = *tmp_Bcc ;
245 Bcc.annule_hard() ;
246 par_B.add_tbl_mod(Bcc) ;
247 Param par_mat_B_hh ;
248
249 Tbl xij_b(np2, nt) ;
250 xij_b.set_etat_qcq() ;
251 initialize_outgoing_BC(nz_bound, B_hh_evol[jtime] , B_hata_evol[jtime], xij_b) ;
252 Tbl xijm1_b(np2, nt) ;
253 xijm1_b.set_etat_qcq() ;
254 initialize_outgoing_BC(nz_bound, B_hh_evol[jtime-1] ,
255 B_hata_evol[jtime-1], xijm1_b) ;
256 Tbl xij_a(np2, nt) ;
257 xij_a.set_etat_qcq() ;
258 initialize_outgoing_BC(nz_bound, A_hh_evol[jtime] , A_hata_evol[jtime], xij_a) ;
259 Tbl xijm1_a(np2, nt) ;
260 xijm1_a.set_etat_qcq() ;
261 initialize_outgoing_BC(nz_bound, A_hh_evol[jtime-1] ,
262 A_hata_evol[jtime-1], xijm1_a) ;
263
264 // Parameters for the Dirac systems
265 //---------------------------------
266
267 Param par_bc_hh ;
268 par_bc_hh.add_int(nz_bound) ;
269 Tbl* cf_b_hh = new Tbl(10) ;
270 cf_b_hh->annule_hard() ;
271 cf_b_hh->set(0) = 11*Rmax + 12*pdt ; // mu
272 cf_b_hh->set(1) = 6*Rmax*pdt ; // d mu / dr
273 cf_b_hh->set(2) = 0 ; // X
274 cf_b_hh->set(3) = 0 ; // d X / dr
275 cf_b_hh->set(4) = 11*Rmax*Rmax + 18*Rmax*pdt ; // h^rr
276 cf_b_hh->set(5) = 6*Rmax*Rmax*pdt ; // d h^rr / dr
277 cf_b_hh->set(6) = 0 ; //eta
278 cf_b_hh->set(7) = 0 ; //d eta / dr
279 cf_b_hh->set(8) = 0 ; //W
280 cf_b_hh->set(9) = 0 ; //d W / dr
281 par_bc_hh.add_tbl_mod(*cf_b_hh, 0) ;
282 Tbl* kib_hh = new Tbl(np2, nt) ;
283 Tbl& khib_hh = *kib_hh ;
284 khib_hh.annule_hard() ;
285 par_bc_hh.add_tbl_mod(khib_hh,1) ;
286 Tbl* mb_hh = new Tbl(np2, nt) ;
287 Tbl& mub_hh = *mb_hh ;
288 mub_hh.annule_hard() ;
289 par_bc_hh.add_tbl_mod(mub_hh, 2) ;
290
291 Param par_mat_hh ;
292
293 Tbl xij_mu_hh(np2, nt) ;
294 xij_mu_hh.set_etat_qcq() ;
295 initialize_outgoing_BC(nz_bound, mu_hh_evol[jtime] , mu_a_evol[jtime], xij_mu_hh) ;
296 Tbl xijm1_mu_hh(np2, nt) ;
297 xijm1_mu_hh.set_etat_qcq() ;
298 initialize_outgoing_BC(nz_bound, mu_hh_evol[jtime-1] , mu_a_evol[jtime-1],
299 xijm1_mu_hh) ;
300
301 Tbl xij_ki_hh(np2, nt) ;
302 xij_ki_hh.set_etat_qcq() ;
303 initialize_outgoing_BC(nz_bound, khi_hh_evol[jtime] , khi_a_evol[jtime], xij_ki_hh) ;
304 Tbl xijm1_ki_hh(np2, nt) ;
305 xijm1_ki_hh.set_etat_qcq() ;
306 initialize_outgoing_BC(nz_bound, khi_hh_evol[jtime-1] , khi_a_evol[jtime-1],
307 xijm1_ki_hh) ;
308
309 Param par_bc_hata ;
310 par_bc_hata.add_int(nz_bound) ;
311 Tbl* cf_b_hata = new Tbl(10) ;
312 cf_b_hata->annule_hard() ;
313 cf_b_hata->set(0) = 11*Rmax + 12*pdt ; // mu
314 cf_b_hata->set(1) = 6*Rmax*pdt ; // d mu / dr
315 cf_b_hata->set(2) = 0 ; // X
316 cf_b_hata->set(3) = 0 ; // d X / dr
317 cf_b_hata->set(4) = 11*Rmax*Rmax + 18*Rmax*pdt ; // h^rr
318 cf_b_hata->set(5) = 6*Rmax*Rmax*pdt ; // d h^rr / dr
319 cf_b_hata->set(6) = 0 ; //eta
320 cf_b_hata->set(7) = 0 ; //d eta / dr
321 cf_b_hata->set(8) = 0 ; //W
322 cf_b_hata->set(9) = 0 ; //d W / dr
323 par_bc_hata.add_tbl_mod(*cf_b_hata, 0) ;
324 Tbl* kib_hata = new Tbl(np2, nt) ;
325 Tbl& khib_hata = *kib_hata ;
326 khib_hata.annule_hard() ;
327 par_bc_hata.add_tbl_mod(khib_hata,1) ;
328 Tbl* mb_hata = new Tbl(np2, nt) ;
329 Tbl& mub_hata = *mb_hata ;
330 mub_hata.annule_hard() ;
331 par_bc_hata.add_tbl_mod(mub_hata, 2) ;
332
333 Param par_mat_hata ;
334
335 Tbl xij_mu_a(np2, nt) ;
336 xij_mu_a.set_etat_qcq() ;
337 initialize_outgoing_BC(nz_bound, mu_a_evol[jtime] ,
338 mu_a_evol.time_derive(jtime, 2), xij_mu_a) ;
339 Tbl xijm1_mu_a(np2, nt) ;
340 xijm1_mu_a.set_etat_qcq() ;
341 tmp = ( mu_a_evol[jtime] - mu_a_evol[jtime-2] ) / (2.*pdt) ;
342 initialize_outgoing_BC(nz_bound, mu_a_evol[jtime-1] , tmp, xijm1_mu_a) ;
343
344 Tbl xij_ki_a(np2, nt) ;
345 xij_ki_a.set_etat_qcq() ;
346 initialize_outgoing_BC(nz_bound, khi_a_evol[jtime] ,
347 khi_a_evol.time_derive(jtime, 2), xij_ki_a) ;
348 Tbl xijm1_ki_a(np2, nt) ;
349 xijm1_ki_a.set_etat_qcq() ;
350 tmp = ( khi_a_evol[jtime] - khi_a_evol[jtime-2] ) / (2.*pdt) ;
351 initialize_outgoing_BC(nz_bound, khi_a_evol[jtime-1] , tmp, xijm1_ki_a) ;
352
353 // Quantities at new time-step
354 //----------------------------
355 Scalar n_new(map) ;
356 Scalar psi_new(map) ;
357 Scalar npsi_new(map) ;
358 Vector beta_new(map, CON, triad) ;
359 Scalar A_hh_new(map) ;
360 Scalar B_hh_new(map) ;
361 Scalar A_hata_new(map) ;
362 Scalar B_hata_new(map) ;
363
364 // Successive values of various quantities:
365 // ---------------------------------------
366 Evolution_full<double> m_adm(adm_mass(), jtime, the_time[jtime]) ;
367 Evolution_full<double> test_ham_constr ;
368 Evolution_full<double> test_mom_constr_r ;
369 Evolution_full<double> test_mom_constr_t ;
370 Evolution_full<double> test_mom_constr_p ;
371 Evolution_full<Tbl> nn_monitor ;
372 Evolution_full<Tbl> psi_monitor ;
373 Evolution_full<Tbl> A_h_monitor ;
374 Evolution_full<Tbl> B_h_monitor ;
375 Evolution_full<Tbl> trh_monitor ;
376 Evolution_full<Tbl> beta_monitor_maxabs ;
377 Evolution_full<Tbl> hh_monitor_central ;
378 Evolution_full<Tbl> hh_monitor_maxabs ;
379 Evolution_full<Tbl> hata_monitor_central ;
380 Evolution_full<Tbl> hata_monitor_maxabs ;
381 Evolution_full<Tbl> check_evol ;
382 Tbl select_scalar(6) ;
383 Tbl select_tens(6) ;
384
385 Vector zero_vec( map, CON, map.get_bvect_spher() ) ;
386 zero_vec.set_etat_zero() ;
387 const Vector& hat_S = ( mom_euler == 0x0 ? zero_vec : *mom_euler ) ;
388 Scalar lapB(map) ;
389 Scalar lapBm1 = source_B_hata_evol[jtime-1] ;
390 Scalar lapBm2 = source_B_hata_evol[jtime-2] ;
391
392 // Evolution loop
393 // --------------
394
395 for (int jt = 0; jt < nb_time_steps; jt++) {
396
397 double ttime = the_time[jtime] ;
398 k_dd() ;
399
400 if (jt%check_mod == 0) {
401 cout <<
402 "==============================================================\n"
403 << " step: " << jtime << " time = " << the_time[jtime] << endl
404 << " ADM mass : " << adm_mass()
405 << ", Log of central lapse: " << log(nn().val_grid_point(0,0,0,0)) << endl
406 << "==============================================================\n" ;
407
408 // Monitoring
409 // ----------
410 m_adm.update(adm_mass(), jtime, the_time[jtime]) ;
411 if (jt > 0) des_evol(m_adm, "ADM mass", "Variation of ADM mass",
412 ngraph0_mon, graph_device) ;
413
414
415 nn_monitor.update(monitor_scalar(nn(), select_scalar),
416 jtime, the_time[jtime]) ;
417
418 psi_monitor.update(monitor_scalar(psi(), select_scalar),
419 jtime, the_time[jtime]) ;
420
421 A_h_monitor.update(monitor_scalar(A_hh(), select_scalar),
422 jtime, the_time[jtime]) ;
423
424 B_h_monitor.update(monitor_scalar(B_hh(), select_scalar),
425 jtime, the_time[jtime]) ;
426
427 trh_monitor.update(monitor_scalar(trh(), select_scalar),
428 jtime, the_time[jtime]) ;
429
430 beta_monitor_maxabs.update(maxabs_all_domains(beta(), -1, 0x0, cout, verbose),
431 jtime, the_time[jtime]) ;
432
433 hh_monitor_central.update(central_value(hh()),
434 jtime, the_time[jtime]) ;
435
436 hh_monitor_maxabs.update(maxabs_all_domains(hh(), -1, 0x0, cout, verbose),
437 jtime, the_time[jtime]) ;
438
439 hata_monitor_central.update(central_value(hata()),
440 jtime, the_time[jtime]) ;
441
442 hata_monitor_maxabs.update(maxabs_all_domains(hata(), -1, 0x0, cout, verbose),
443 jtime, the_time[jtime]) ;
444
445
446 int jt_graph = jt / check_mod ;
447
448 Tbl tham = check_hamiltonian_constraint(0x0, cout, verbose) ;
449 double max_error = tham(0,0) ;
450 for (int l=1; l<nz-1; l++) { // all domains but the last one
451 double xx = fabs(tham(0,l)) ;
452 if (xx > max_error) max_error = xx ;
453 }
454 test_ham_constr.update(max_error, jt_graph, the_time[jtime]) ;
455 if (jt > 0) des_evol(test_ham_constr, "Absolute error",
456 "Check of Hamiltonian constraint",
457 ngraph0_mon+1, graph_device) ;
458
459 Tbl tmom = check_momentum_constraint(0x0, cout, verbose) ;
460 max_error = tmom(0,0) ;
461 for (int l=1; l<nz-1; l++) { // all domains but the last one
462 double xx = fabs(tmom(0,l)) ;
463 if (xx > max_error) max_error = xx ;
464 }
465 test_mom_constr_r.update(max_error, jt_graph, the_time[jtime]) ;
466 if (jt > 0) des_evol(test_mom_constr_r, "Absolute error",
467 "Check of momentum constraint (r comp.)", ngraph0_mon+2,
468 graph_device) ;
469
470 max_error = tmom(1,0) ;
471 for (int l=1; l<nz-1; l++) { // all domains but the last one
472 double xx = fabs(tmom(1,l)) ;
473 if (xx > max_error) max_error = xx ;
474 }
475 test_mom_constr_t.update(max_error, jt_graph, the_time[jtime]) ;
476 if (jt > 0) des_evol(test_mom_constr_t, "Absolute error",
477 "Check of momentum constraint (\\gh comp.)", ngraph0_mon+3,
478 graph_device) ;
479
480 max_error = tmom(2,0) ;
481 for (int l=1; l<nz-1; l++) { // all domains but the last one
482 double xx = fabs(tmom(2,l)) ;
483 if (xx > max_error) max_error = xx ;
484 }
485 test_mom_constr_p.update(max_error, jt_graph, the_time[jtime]) ;
486 if (jt > 0) des_evol(test_mom_constr_p, "Absolute error",
487 "Check of momentum constraint (\\gf comp.)", ngraph0_mon+4,
488 graph_device) ;
489
490 if (jt>2) {
491 Tbl tevol = check_dynamical_equations(0x0, 0x0, cout, verbose) ;
492 Tbl evol_check(6) ; evol_check.set_etat_qcq() ;
493 for (int i=1; i<=3; i++)
494 for(int j=1; j<=i; j++) {
495 max_error = tevol(i, j, 0) ;
496 for (int l=1; l<nz-1; l++) {
497 double xx = fabs(tevol(i,j,l)) ;
498 if (xx > max_error) max_error = xx ;
499 }
500 evol_check.set(i) = max_error ;
501 }
502 check_evol.update(evol_check, jtime, the_time[jtime]) ;
503 }
504 }
505
506 if (jt%save_mod == 0) {
507 m_adm.save("adm_mass.d") ;
508 nn_monitor.save("nn_monitor.d") ;
509 psi_monitor.save("psi_monitor.d") ;
510 A_h_monitor.save("potA_monitor.d") ;
511 B_h_monitor.save("potB_monitor.d") ;
512 trh_monitor.save("trh_monitor.d") ;
513 beta_monitor_maxabs.save("beta_monitor_maxabs.d") ;
514 hh_monitor_central.save("hh_monitor_central.d") ;
515 hh_monitor_maxabs.save("hh_monitor_maxabs.d") ;
516 hata_monitor_central.save("hata_monitor_central.d") ;
517 hata_monitor_maxabs.save("hata_monitor_maxabs.d") ;
518 test_ham_constr.save("test_ham_constr.d") ;
519 test_mom_constr_r.save("test_mom_constr_r.d") ;
520 test_mom_constr_t.save("test_mom_constr_t.d") ;
521 test_mom_constr_p.save("test_mom_constr_p.d") ;
522 check_evol.save("evol_equations.d") ;
523
524 save("sigma") ;
525
526 }
527
528
529 // Resolution of hyperbolic equations
530 // ----------------------------------
531 compute_sources(strain_euler) ;
532
533 A_hata_new = A_hata_evol[jtime]
534 + pdt*( an*source_A_hata_evol[jtime] + anm1*source_A_hata_evol[jtime-1]
535 + anm2*source_A_hata_evol[jtime-2] ) ;
536 B_hata_new = B_hata_evol[jtime]
537 + pdt*( an*source_B_hata_evol[jtime] + anm1*source_B_hata_evol[jtime-1]
538 + anm2*source_B_hata_evol[jtime-2] ) ;
539
540 A_hh_new = A_hh_evol[jtime]
541 + pdt*( an*source_A_hh_evol[jtime] + anm1*source_A_hh_evol[jtime-1]
542 + anm2*source_A_hh_evol[jtime-2] ) ;
543
544 B_hh_new = B_hh_evol[jtime]
545 + pdt*( an*source_B_hh_evol[jtime] + anm1*source_B_hh_evol[jtime-1]
546 + anm2*source_B_hh_evol[jtime-2] ) ;
547
548 Scalar bc_A = -2.*A_hata_new ;
549 bc_A.set_spectral_va().ylm() ;
550 evolve_outgoing_BC(pdt, nz_bound, A_hh_evol[jtime], bc_A, xij_a, xijm1_a,
551 Acc, 0) ;
552 A_hh_new.match_tau(par_A, &par_mat_A_hh) ;
553
554 Scalar bc_B = -2.*B_hata_new ;
555 bc_B.set_spectral_va().ylm() ;
556 evolve_outgoing_BC(pdt, nz_bound, B_hh_evol[jtime], bc_B, xij_b, xijm1_b,
557 Bcc, -1) ;
558 B_hh_new.match_tau(par_B, &par_mat_B_hh) ;
559
560 // Boundary conditions for hh and hata
561 //------------------------------------
562 Scalar sbcmu = (18*mu_hh_evol[jtime] - 9*mu_hh_evol[jtime-1]
563 + 2*mu_hh_evol[jtime-2]) / (6*pdt) ;
564 if (sbcmu.get_etat() == ETATZERO) {
565 sbcmu.annule_hard() ;
566 sbcmu.set_spectral_base(base_pseudo) ;
567 }
568 sbcmu.set_spectral_va().ylm() ;
569 tmp = mu_hh_evol[jtime] ;
570 if (tmp.get_etat() == ETATZERO) {
571 tmp.annule_hard() ;
572 tmp.set_spectral_base(base_pseudo) ;
573 }
574 tmp.set_spectral_va().ylm() ;
575 evolve_outgoing_BC(pdt, nz_bound, tmp, sbcmu, xij_mu_hh, xijm1_mu_hh,
576 mub_hh, 0) ;
577 mub_hh *= 6*pdt ;
578
579 Scalar sbckhi = (18*khi_hh_evol[jtime] - 9*khi_hh_evol[jtime-1]
580 + 2*khi_hh_evol[jtime-2]) / (6*pdt) ;
581 if (sbckhi.get_etat() == ETATZERO) {
582 sbckhi.annule_hard() ;
583 sbckhi.set_spectral_base(base_ref) ;
584 }
585 sbckhi.set_spectral_va().ylm() ;
586 tmp = khi_hh_evol[jtime] ;
587 if (tmp.get_etat() == ETATZERO) {
588 tmp.annule_hard() ;
589 tmp.set_spectral_base(base_ref) ;
590 }
591 tmp.set_spectral_va().ylm() ;
592 evolve_outgoing_BC(pdt, nz_bound, tmp, sbckhi, xij_ki_hh, xijm1_ki_hh,
593 khib_hh, 0) ;
594 khib_hh *= 6*pdt ;
595
596 sbcmu = (18*mu_a_evol[jtime] - 9*mu_a_evol[jtime-1]
597 + 2*mu_a_evol[jtime-2]) / (6*pdt) ;
598 if (sbcmu.get_etat() == ETATZERO) {
599 sbcmu.annule_hard() ;
600 sbcmu.set_spectral_base(base_pseudo) ;
601 }
602 sbcmu.set_spectral_va().ylm() ;
603 tmp = mu_a_evol[jtime] ;
604 if (tmp.get_etat() == ETATZERO) {
605 tmp.annule_hard() ;
606 tmp.set_spectral_base(base_pseudo) ;
607 }
608 tmp.set_spectral_va().ylm() ;
609 evolve_outgoing_BC(pdt, nz_bound, tmp, sbcmu, xij_mu_a, xijm1_mu_a,
610 mub_hata, 0) ;
611 mub_hata *= 6*pdt ;
612
613 sbckhi = (18*khi_a_evol[jtime] - 9*khi_a_evol[jtime-1]
614 + 2*khi_a_evol[jtime-2]) / (6*pdt) ;
615 if (sbckhi.get_etat() == ETATZERO) {
616 sbckhi.annule_hard() ;
617 sbckhi.set_spectral_base(base_ref) ;
618 }
619 sbckhi.set_spectral_va().ylm() ;
620 tmp = khi_a_evol[jtime] ;
621 if (tmp.get_etat() == ETATZERO) {
622 tmp.annule_hard() ;
623 tmp.set_spectral_base(base_ref) ;
624 }
625 tmp.set_spectral_va().ylm() ;
626 evolve_outgoing_BC(pdt, nz_bound, tmp, sbckhi, xij_ki_a, xijm1_ki_a,
627 khib_hata, 0) ;
628 khib_hata *= 6*pdt ;
629
630 // Advance in time
631 // ---------------
632
633 jtime++ ;
634 ttime += pdt ;
635 the_time.update(ttime, jtime, ttime) ;
636
637 // Setting As and Bs for h^{ij} and \hat{A}^{ij}
638 set_AB_hata(A_hata_new, B_hata_new) ;
639 set_AB_hh(A_hh_new, B_hh_new) ;
640
641 hij_tt.set_A_tildeB( A_hh_new, B_hh_new, &par_bc_hh, &par_mat_hh ) ;
642 for (int i=1; i<=3; i++)
643 for (int j=i; j<=3; j++)
644 for (int l=nz_bound+1; l<nz; l++)
645 hij_tt.set(i,j).set_domain(l) = hijtt_old(i,j).domain(l) ;
646 hata_tt.set_A_tildeB( A_hata_new, B_hata_new, &par_bc_hata, &par_mat_hata ) ;
647 for (int i=1; i<=3; i++)
648 for (int j=i; j<=3; j++) {
649 for (int l=nz_bound+1; l<nz; l++)
650 hata_tt.set(i,j).set_domain(l) = hatatt_old(i,j).domain(l) ;
651 hata_tt.set(i,j).set_dzpuis(2) ;
652 }
653
654 // Computation of h^{ij} at new time-step
655 hh_det_one(hij_tt, &par_mat_hh) ;
656
657 // Reset of derived quantities
658 del_deriv() ;
659
660 // Update of khi's and mu's
661 //-------------------------
662 tmp = hij_tt( 1, 1 ) ;
663 tmp.mult_r() ; tmp.mult_r() ;
664 khi_hh_evol.update( tmp, jtime, the_time[jtime] ) ;
665 mu_hh_evol.update( hij_tt.mu(), jtime, the_time[jtime] ) ;
666 tmp = hata_tt( 1, 1 ) ;
667 tmp.mult_r() ; tmp.mult_r() ;
668 khi_a_evol.update( tmp, jtime, the_time[jtime] ) ;
669 mu_a_evol.update( hata_tt.mu(), jtime, the_time[jtime] ) ;
670
671 // Resolution of elliptic equations
672 // --------------------------------
673 psi_evol.update(psi_evol[jtime-1], jtime, ttime) ;
674
675 // \hat{A}^{ij} is computed at the new time-step
676 compute_X_from_momentum_constraint(hat_S, hata_tt, niter_elliptic) ;
677
678 // Iteration on the conformal factor
679 for (int k = 0; k < niter_elliptic; k++) {
680
681 psi_new = solve_psi(ener_euler) ;
682 psi_new = relax * psi_new + (1.-relax) * psi() ;
683 set_psi_del_npsi(psi_new) ;
684 }
685
686 set_npsi_del_n(npsi_evol[jtime-1]) ;
687
688 // Iteration on N*Psi ## play with the number of iterations...
689 npsi_evol.update(psi_evol[jtime-1], jtime, ttime) ;
690 for (int k = 0; k < niter_elliptic; k++) {
691
692 npsi_new = solve_npsi( ener_euler, s_euler ) ;
693 npsi_new = relax * npsi_new + (1.-relax) * npsi() ;
694 set_npsi_del_n(npsi_new) ;
695 }
696
697 // Iteration on beta ## play with the number of iterations...
698 beta_evol.update(beta_evol[jtime-1], jtime, ttime) ;
699 for (int k = 0; k < niter_elliptic; k++) {
700
701 beta_new = solve_beta(method_poisson_vect) ;
702 beta_new = relax * beta_new + (1.-relax) * beta() ;
703 beta_evol.update(beta_new, jtime, ttime) ;
704 }
705
706 des_meridian(vec_X()(1), 0., ray_des, "\\gb\\ur\\d", ngraph0+6,
707 graph_device) ;
708 des_meridian(vec_X()(2), 0., ray_des, "\\gb\\u\\gh\\d", ngraph0+7,
709 graph_device) ;
710 des_meridian(vec_X()(3), 0., ray_des, "\\gb\\u\\gf\\d", ngraph0+8,
711 graph_device) ;
712 tmp = A_hh() ;
713 tmp.set_spectral_va().ylm_i() ;
714 des_meridian(tmp, 0., ray_des, "A\\dh", ngraph0+9,
715 graph_device) ;
716 tmp = B_hh_new;
717 tmp.set_spectral_va().ylm_i() ;
718 des_meridian(tmp, 0., ray_des, "B\\dh", ngraph0+10,
719 graph_device) ;
720 des_meridian(trh(), 0., ray_des, "tr h", ngraph0+11,
721 graph_device) ;
722 des_meridian(hh()(1,1), 0., ray_des, "h\\urr\\d", ngraph0+12,
723 graph_device) ;
724 des_meridian(hh()(2,3), 0., ray_des, "h\\u\\gh\\gf\\d", ngraph0+13,
725 graph_device) ;
726 des_meridian(hh()(3,3), 0., ray_des, "h\\u\\gf\\gf\\d", ngraph0+14,
727 graph_device) ;
728
729 arrete(nopause) ;
730 }
731
732 par_A.clean_all() ;
733 par_B.clean_all() ;
734 par_mat_A_hh.clean_all() ;
735 par_mat_B_hh.clean_all() ;
736
737 par_bc_hh.clean_all() ;
738 par_mat_hh.clean_all() ;
739
740 par_bc_hata.clean_all() ;
741 par_mat_hata.clean_all() ;
742}
743
744
745//***************************************************************************
746
747const Tbl& monitor_scalar(const Scalar& uu, Tbl& resu) {
748
749 assert( resu.get_ndim() == 1) ;
750 assert( resu.get_taille() >= 6) ;
751
752 resu.set_etat_qcq() ;
753
754 resu.set(0) = uu.val_grid_point(0,0,0,0) ;
755 resu.set(1) = max(max(uu)) ;
756 resu.set(2) = min(min(uu)) ;
757
758 const Mg3d& mg = *(uu.get_mp().get_mg()) ;
759
760 int nz = mg.get_nzone() ;
761 int nzm1 = nz - 1 ;
762 int nr = mg.get_nr(nzm1) ;
763 int nt = mg.get_nt(nzm1) ;
764 int np = mg.get_np(nzm1) ;
765
766 resu.set(3) = uu.val_grid_point(nzm1, 0, 0, nr-1) ;
767 resu.set(4) = uu.val_grid_point(nzm1, 0, nt-1, nr-1) ;
768 resu.set(5) = uu.val_grid_point(nzm1, np/2, nt-1, nr-1) ;
769
770 return resu ;
771}
772}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:325
Lorene::Base_val::mult_x
void mult_x()
The basis is transformed as with a multiplication by .
Definition base_val_manip.C:160
Lorene::Base_val::mult_cost
void mult_cost()
The basis is transformed as with a multiplication.
Definition base_val_manip.C:337
Lorene::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Lorene::Evolution_full
Time evolution with full storage (*** under development ***).
Definition evolution.h:270
Lorene::Evolution_full::update
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
Definition evolution_full.C:176
Lorene::Evolution_std
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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::Evolution::save
void save(const char *filename) const
Saves *this in a formatted file.
Definition evolution.C:589
Lorene::Evolution::time_derive
TyT time_derive(int j, int n=2) const
Computes the time derivative at time step j by means of a n-th order scheme, from the values at steps...
Definition evolution.C:504
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:479
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:474
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:469
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::clean_all
void clean_all()
Deletes all the objects stored as modifiables, i.e.
Definition param.C:177
Lorene::Param::add_double_mod
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition param.C:456
Lorene::Param::add_tbl_mod
void add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
Definition param.C:594
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:249
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::std_spectral_base
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
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:643
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:621
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:610
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:386
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:560
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1328
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:814
Lorene::Scalar::mult_r
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:211
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:803
Lorene::Sym_tensor_trans
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:611
Lorene::Sym_tensor_tt
Transverse and traceless symmetric tensors of rank 2.
Definition sym_tensor.h:933
Lorene::Sym_tensor_tt::set_A_tildeB
void set_A_tildeB(const Scalar &a_in, const Scalar &tb_in, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A and .
Definition sym_tensor_tt_etamu.C:490
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:226
Lorene::Sym_tensor::mu
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
Definition sym_tensor_aux.C:154
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim).
Definition tbl.h:400
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Time_slice_conf::npsi
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Definition time_slice_conf.C:735
Lorene::Time_slice_conf::hh_evol
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:533
Lorene::Time_slice_conf::k_dd
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ).
Definition time_slice_conf.C:633
Lorene::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime ).
Definition time_slice_conf.C:594
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:407
Lorene::Time_slice_conf::hata
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Definition time_slice_conf.C:775
Lorene::Time_slice_conf::hata_evol
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition time_slice.h:545
Lorene::Time_slice_conf::vec_X
virtual const Vector & vec_X(int method_poisson=6) const
Vector representing the longitudinal part of .
Definition time_slice_conf.C:828
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:696
Lorene::Time_slice_conf::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition time_slice_conf.C:350
Lorene::Time_slice_conf::npsi_evol
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition time_slice.h:525
Lorene::Time_slice_conf::set_npsi_del_n
virtual void set_npsi_del_n(const Scalar &npsi_in)
Sets the factor at the current time step (jtime ) and deletes the value of N.
Definition time_slice_conf.C:481
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:520
Lorene::Time_slice_conf::set_AB_hata
virtual void set_AB_hata(const Scalar &A_in, const Scalar &B_in)
Sets the potentials A and of the TT part (see the documentation of Sym_tensor for details).
Definition time_slice_conf.C:894
Lorene::Time_slice_conf::A_hata_evol
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition time_slice.h:550
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:510
Lorene::Time_slice_conf::compute_X_from_momentum_constraint
void compute_X_from_momentum_constraint(const Vector &hat_S, const Sym_tensor_tt &hata_tt, int iter_max=200, double precis=1.e-12, double relax=0.8, int methode_poisson=6)
Computes the vector from the conformally-rescaled momentum , using the momentum constraint.
Definition time_slice_conf.C:840
Lorene::Time_slice_conf::B_hata_evol
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition time_slice.h:555
Lorene::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition time_slice.h:193
Lorene::Time_slice::save
void save(const char *rootname) const
Saves in a binary file.
Definition time_slice.C:464
Lorene::Time_slice::check_hamiltonian_constraint
Tbl check_hamiltonian_constraint(const Scalar *energy_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the hamiltonian constraint is verified.
Definition tslice_check_einstein.C:82
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition time_slice_access.C:90
Lorene::Time_slice::depth
int depth
Number of stored time slices.
Definition time_slice.h:182
Lorene::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:196
Lorene::Time_slice::check_dynamical_equations
Tbl check_dynamical_equations(const Sym_tensor *strain_tensor=0x0, const Scalar *energy_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the dynamical equations are verified.
Definition tslice_check_einstein.C:142
Lorene::Time_slice::beta_evol
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:222
Lorene::Time_slice::check_momentum_constraint
Tbl check_momentum_constraint(const Vector *momentum_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the momentum constraints are verified.
Definition tslice_check_einstein.C:112
Lorene::Tslice_dirac_max::trh
virtual const Scalar & trh() const
Computes the trace h, with respect to the flat metric ff , of .
Definition tslice_dirac_max.C:550
Lorene::Tslice_dirac_max::adm_mass
virtual double adm_mass() const
Returns the ADM mass at (geometrical units) the current step.
Definition tslice_adm_mass.C:146
Lorene::Tslice_dirac_max::hh
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition tslice_dirac_max.C:477
Lorene::Tslice_dirac_max::solve_psi
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $\Psi$ (Hamiltonian constraint).
Definition tslice_dirac_max_solve.C:175
Lorene::Tslice_dirac_max::source_B_hh_evol
Evolution_std< Scalar > source_B_hh_evol
The potential of the source of equation for .
Definition time_slice.h:998
Lorene::Tslice_dirac_max::source_A_hh_evol
Evolution_std< Scalar > source_A_hh_evol
The A potential of the source of equation for .
Definition time_slice.h:992
Lorene::Tslice_dirac_max::evolve
void evolve(double pdt, int nb_time_steps, int niter_elliptic, double relax_elliptic, int check_mod, int save_mod, int method_poisson_vect=6, int nopause=1, const char *graph_device=0x0, bool verbose=true, const Scalar *ener_euler=0x0, const Vector *mom_euler=0x0, const Scalar *s_euler=0x0, const Sym_tensor *strain_euler=0x0)
Time evolution by resolution of Einstein equations.
Definition tslice_dirac_max_evolve.C:123
Lorene::Tslice_dirac_max::B_hh
virtual const Scalar & B_hh() const
Returns the potential of .
Definition tslice_dirac_max.C:539
Lorene::Tslice_dirac_max::A_hh_evol
Evolution_std< Scalar > A_hh_evol
The A potential of .
Definition time_slice.h:980
Lorene::Tslice_dirac_max::hh_det_one
void hh_det_one(int j, Param *par_bc=0x0, Param *par_mat=0x0) const
Computes from the values of A and and using the condition , which fixes the trace of .
Definition tslice_dirac_max_setAB.C:112
Lorene::Tslice_dirac_max::source_B_hata_evol
Evolution_std< Scalar > source_B_hata_evol
The potential of the source of equation for .
Definition time_slice.h:1010
Lorene::Tslice_dirac_max::compute_sources
void compute_sources(const Sym_tensor *strain_tensor=0x0) const
Computes the sources source_A_XXX_evol and source_B_XXX_evol , for the solution of the evolution equa...
Definition tslice_dirac_max_setAB.C:222
Lorene::Tslice_dirac_max::solve_npsi
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).
Definition tslice_dirac_max_solve.C:121
Lorene::Tslice_dirac_max::A_hh
virtual const Scalar & A_hh() const
Returns the potential A of .
Definition tslice_dirac_max.C:528
Lorene::Tslice_dirac_max::solve_beta
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Definition tslice_dirac_max_solve.C:227
Lorene::Tslice_dirac_max::set_AB_hh
virtual void set_AB_hh(const Scalar &A_in, const Scalar &B_in)
Sets the potentials A and of the TT part of (see the documentation of Sym_tensor for details).
Definition tslice_dirac_max_setAB.C:87
Lorene::Tslice_dirac_max::B_hh_evol
Evolution_std< Scalar > B_hh_evol
The potential of .
Definition time_slice.h:986
Lorene::Tslice_dirac_max::source_A_hata_evol
Evolution_std< Scalar > source_A_hata_evol
The potential A of the source of equation for .
Definition time_slice.h:1004
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:134
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:461
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:438
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:299
Lorene::arrete
void arrete(int a=0)
Setting a stop point in a code.
Definition arrete.C:64
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:874
Lorene::Tensor::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:680
Lorene::Tensor::inc_dzpuis
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
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:879
Lorene::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:506
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:663
Lorene::maxabs_all_domains
Tbl maxabs_all_domains(const Tensor &aa, int l_excluded=-1, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maximum of the absolute value of each component of a tensor over all the domains.
Definition tensor_calculus_ext.C:791
Lorene::central_value
Tbl central_value(const Tensor &aa, const char *comment=0x0, ostream &ost=cout)
Central value of each component of a tensor.
Definition tensor_calculus_ext.C:652
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:777
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142