114 #include "time_slice.h" 116 #include "graphique.h" 117 #include "utilitaires.h" 121 const Tbl& monitor_scalar(
const Scalar& uu, Tbl& resu) ;
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,
140 int ngraph0_mon = 70 ;
142 int nz_bound = nz - 2 ;
152 for (
int lz=1; lz<nz; lz++) {
157 for (
int it=0; it<
depth; it++) {
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 )
176 hijtt_old.
annule(0, nz_bound) ;
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 )
186 hatatt_old.
annule(0, nz_bound) ;
206 double Rmax = map.
val_r(nz-2, 1., 0., 0.) ;
207 double ray_des = 1.25 * Rmax ;
211 double an = 23./12. ;
212 double anm1 = -4./3. ;
213 double anm2 = 5./12. ;
216 int i_minus_one = -1 ;
220 double *l_val_A =
new double(1./Rmax) ;
221 double *l_der_A =
new double(1.) ;
228 Tbl* tmp_Acc =
new Tbl(np2, nt) ;
229 Tbl& Acc = *tmp_Acc ;
235 double* l_val_B =
new double(1./Rmax) ;
236 double* l_der_B =
new double(1.) ;
239 par_B.
add_int(i_minus_one, 2) ;
243 Tbl* tmp_Bcc =
new Tbl(np2, nt) ;
244 Tbl& Bcc = *tmp_Bcc ;
252 Tbl xijm1_b(np2, nt) ;
259 Tbl xijm1_a(np2, nt) ;
269 Tbl* cf_b_hh =
new Tbl(10) ;
271 cf_b_hh->
set(0) = 11*Rmax + 12*pdt ;
272 cf_b_hh->
set(1) = 6*Rmax*pdt ;
273 cf_b_hh->
set(2) = 0 ;
274 cf_b_hh->
set(3) = 0 ;
275 cf_b_hh->
set(4) = 11*Rmax*Rmax + 18*Rmax*pdt ;
276 cf_b_hh->
set(5) = 6*Rmax*Rmax*pdt ;
277 cf_b_hh->
set(6) = 0 ;
278 cf_b_hh->
set(7) = 0 ;
279 cf_b_hh->
set(8) = 0 ;
280 cf_b_hh->
set(9) = 0 ;
282 Tbl* kib_hh =
new Tbl(np2, nt) ;
283 Tbl& khib_hh = *kib_hh ;
286 Tbl* mb_hh =
new Tbl(np2, nt) ;
287 Tbl& mub_hh = *mb_hh ;
293 Tbl xij_mu_hh(np2, nt) ;
295 initialize_outgoing_BC(nz_bound, mu_hh_evol[
jtime] , mu_a_evol[
jtime], xij_mu_hh) ;
296 Tbl xijm1_mu_hh(np2, nt) ;
298 initialize_outgoing_BC(nz_bound, mu_hh_evol[
jtime-1] , mu_a_evol[
jtime-1],
301 Tbl xij_ki_hh(np2, nt) ;
303 initialize_outgoing_BC(nz_bound, khi_hh_evol[
jtime] , khi_a_evol[
jtime], xij_ki_hh) ;
304 Tbl xijm1_ki_hh(np2, nt) ;
306 initialize_outgoing_BC(nz_bound, khi_hh_evol[
jtime-1] , khi_a_evol[
jtime-1],
310 par_bc_hata.
add_int(nz_bound) ;
311 Tbl* cf_b_hata =
new Tbl(10) ;
313 cf_b_hata->
set(0) = 11*Rmax + 12*pdt ;
314 cf_b_hata->
set(1) = 6*Rmax*pdt ;
315 cf_b_hata->
set(2) = 0 ;
316 cf_b_hata->
set(3) = 0 ;
317 cf_b_hata->
set(4) = 11*Rmax*Rmax + 18*Rmax*pdt ;
318 cf_b_hata->
set(5) = 6*Rmax*Rmax*pdt ;
319 cf_b_hata->
set(6) = 0 ;
320 cf_b_hata->
set(7) = 0 ;
321 cf_b_hata->
set(8) = 0 ;
322 cf_b_hata->
set(9) = 0 ;
324 Tbl* kib_hata =
new Tbl(np2, nt) ;
325 Tbl& khib_hata = *kib_hata ;
328 Tbl* mb_hata =
new Tbl(np2, nt) ;
329 Tbl& mub_hata = *mb_hata ;
335 Tbl xij_mu_a(np2, nt) ;
337 initialize_outgoing_BC(nz_bound, mu_a_evol[
jtime] ,
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) ;
344 Tbl xij_ki_a(np2, nt) ;
346 initialize_outgoing_BC(nz_bound, khi_a_evol[
jtime] ,
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) ;
358 Vector beta_new(map, CON, triad) ;
382 Tbl select_scalar(6) ;
387 const Vector& hat_S = ( mom_euler == 0x0 ? zero_vec : *mom_euler ) ;
395 for (
int jt = 0; jt < nb_time_steps; jt++) {
400 if (jt%check_mod == 0) {
402 "==============================================================\n" 405 <<
", Log of central lapse: " <<
log(
nn().val_grid_point(0,0,0,0)) << endl
406 <<
"==============================================================\n" ;
411 if (jt > 0)
des_evol(m_adm,
"ADM mass",
"Variation of ADM mass",
412 ngraph0_mon, graph_device) ;
415 nn_monitor.
update(monitor_scalar(
nn(), select_scalar),
418 psi_monitor.
update(monitor_scalar(
psi(), select_scalar),
421 A_h_monitor.
update(monitor_scalar(
A_hh(), select_scalar),
424 B_h_monitor.
update(monitor_scalar(
B_hh(), select_scalar),
427 trh_monitor.
update(monitor_scalar(
trh(), select_scalar),
446 int jt_graph = jt / check_mod ;
449 double max_error = tham(0,0) ;
450 for (
int l=1; l<nz-1; l++) {
451 double xx = fabs(tham(0,l)) ;
452 if (xx > max_error) max_error = xx ;
455 if (jt > 0)
des_evol(test_ham_constr,
"Absolute error",
456 "Check of Hamiltonian constraint",
457 ngraph0_mon+1, graph_device) ;
460 max_error = tmom(0,0) ;
461 for (
int l=1; l<nz-1; l++) {
462 double xx = fabs(tmom(0,l)) ;
463 if (xx > max_error) max_error = xx ;
466 if (jt > 0)
des_evol(test_mom_constr_r,
"Absolute error",
467 "Check of momentum constraint (r comp.)", ngraph0_mon+2,
470 max_error = tmom(1,0) ;
471 for (
int l=1; l<nz-1; l++) {
472 double xx = fabs(tmom(1,l)) ;
473 if (xx > max_error) max_error = xx ;
476 if (jt > 0)
des_evol(test_mom_constr_t,
"Absolute error",
477 "Check of momentum constraint (\\gh comp.)", ngraph0_mon+3,
480 max_error = tmom(2,0) ;
481 for (
int l=1; l<nz-1; l++) {
482 double xx = fabs(tmom(2,l)) ;
483 if (xx > max_error) max_error = xx ;
486 if (jt > 0)
des_evol(test_mom_constr_p,
"Absolute error",
487 "Check of momentum constraint (\\gf comp.)", ngraph0_mon+4,
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 ;
500 evol_check.
set(i) = max_error ;
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") ;
548 Scalar bc_A = -2.*A_hata_new ;
550 evolve_outgoing_BC(pdt, nz_bound,
A_hh_evol[
jtime], bc_A, xij_a, xijm1_a,
552 A_hh_new.match_tau(par_A, &par_mat_A_hh) ;
554 Scalar bc_B = -2.*B_hata_new ;
556 evolve_outgoing_BC(pdt, nz_bound,
B_hh_evol[
jtime], bc_B, xij_b, xijm1_b,
558 B_hh_new.match_tau(par_B, &par_mat_B_hh) ;
563 + 2*mu_hh_evol[
jtime-2]) / (6*pdt) ;
569 tmp = mu_hh_evol[
jtime] ;
575 evolve_outgoing_BC(pdt, nz_bound, tmp, sbcmu, xij_mu_hh, xijm1_mu_hh,
580 + 2*khi_hh_evol[
jtime-2]) / (6*pdt) ;
581 if (sbckhi.
get_etat() == ETATZERO) {
586 tmp = khi_hh_evol[
jtime] ;
592 evolve_outgoing_BC(pdt, nz_bound, tmp, sbckhi, xij_ki_hh, xijm1_ki_hh,
596 sbcmu = (18*mu_a_evol[
jtime] - 9*mu_a_evol[
jtime-1]
597 + 2*mu_a_evol[
jtime-2]) / (6*pdt) ;
603 tmp = mu_a_evol[
jtime] ;
609 evolve_outgoing_BC(pdt, nz_bound, tmp, sbcmu, xij_mu_a, xijm1_mu_a,
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) {
620 tmp = khi_a_evol[
jtime] ;
626 evolve_outgoing_BC(pdt, nz_bound, tmp, sbckhi, xij_ki_a, xijm1_ki_a,
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++)
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++)
662 tmp = hij_tt( 1, 1 ) ;
666 tmp = hata_tt( 1, 1 ) ;
679 for (
int k = 0; k < niter_elliptic; k++) {
682 psi_new = relax * psi_new + (1.-relax) *
psi() ;
690 for (
int k = 0; k < niter_elliptic; k++) {
692 npsi_new =
solve_npsi( ener_euler, s_euler ) ;
693 npsi_new = relax * npsi_new + (1.-relax) *
npsi() ;
699 for (
int k = 0; k < niter_elliptic; k++) {
702 beta_new = relax * beta_new + (1.-relax) *
beta() ;
724 des_meridian(
hh()(2,3), 0., ray_des,
"h\\u\\gh\\gf\\d", ngraph0+13,
726 des_meridian(
hh()(3,3), 0., ray_des,
"h\\u\\gf\\gf\\d", ngraph0+14,
747 const Tbl& monitor_scalar(
const Scalar& uu,
Tbl& resu) {
762 int nr = mg.
get_nr(nzm1) ;
763 int nt = mg.
get_nt(nzm1) ;
764 int np = mg.
get_np(nzm1) ;
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
void add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Cmp log(const Cmp &)
Neperian logarithm.
void mult_cost()
The basis is transformed as with a multiplication.
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
virtual double adm_mass() const
Returns the ADM mass at (geometrical units) the current step.
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
void ylm_i()
Inverse of ylm()
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
virtual void del_deriv() const
Deletes all the derived quantities.
Tbl central_value(const Tensor &aa, const char *comment=0x0, ostream &ost=cout)
Central value of each component of a tensor.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
void save(const char *rootname) const
Saves in a binary file.
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.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Base class for coordinate mappings.
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
int jtime
Time step index of the latest slice.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
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.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
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.
Tensor field of valence 1.
void clean_all()
Deletes all the objects stored as modifiables, i.e.
Vectorial bases (triads) with respect to which the tensorial components are defined.
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
void set_dzpuis(int)
Modifies the dzpuis flag.
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void annule_hard()
Sets the Scalar to zero in a hard way.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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 .
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
void des_evol(const Evolution< double > &uu, const char *nomy, const char *title, int ngraph, const char *device, bool closeit, bool show_time, const char *nomx)
Plots the variation of some quantity against time.
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.
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $$ (Hamiltonian constraint).
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Evolution_std< Scalar > source_A_hata_evol
The potential A of the source of equation for .
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. ...
Evolution_std< Scalar > source_A_hh_evol
The A potential of the source of equation for .
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...
virtual const Vector & vec_X(int method_poisson=6) const
Vector representing the longitudinal part of .
int get_nzone() const
Returns the number of domains.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
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 .
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 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 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 > source_B_hh_evol
The potential of the source of equation for .
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Transverse symmetric tensors of rank 2.
Time evolution with partial storage (*** under development ***).
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ) ...
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Evolution_std< double > the_time
Time label of each slice.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
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...
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
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 ...
Bases of the spectral expansions.
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Evolution_std< Scalar > B_hh_evol
The potential of .
virtual const Scalar & A_hh() const
Returns the potential A of .
int get_taille() const
Gives the total size (ie dim.taille)
void mult_x()
The basis is transformed as with a multiplication by .
void arrete(int a=0)
Setting a stop point in a code.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Valeur & set_spectral_va()
Returns va (read/write version)
void save(const char *filename) const
Saves *this in a formatted file.
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
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< Vector > beta_evol
Values at successive time steps of the shift vector .
const Map & get_mp() const
Returns the mapping.
int depth
Number of stored time slices.
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
void annule_hard()
Sets the Tbl to zero in a hard way.
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
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)...
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.
Class intended to describe valence-2 symmetric tensors.
Evolution_std< Scalar > source_B_hata_evol
The potential of the source of equation for .
Transverse and traceless symmetric tensors of rank 2.
virtual const Scalar & trh() const
Computes the trace h, with respect to the flat metric ff , of .
Time evolution with full storage (*** under development ***).
Evolution_std< Scalar > A_hh_evol
The A potential of .
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )