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