43 #include "utilitaires.h" 50 #include "param_elliptic.h" 52 #include "sym_tensor.h" 56 #include "excision_surf.h" 73 expa(sph.theta_plus()),
74 dt_expa(sph.theta_plus())
91 conf_fact(exc_in.conf_fact),
96 delta_t(exc_in.delta_t),
97 no_of_steps(exc_in.no_of_steps),
99 dt_expa(exc_in.dt_expa)
196 cout <<
"===========================================================================================" << endl;
197 cout <<
"Starting the routine that computes parameters for the hyperbolic evolution of the expansion" << endl;
198 cout <<
"Those parameters should be set once and for all: do not call that routine twice" <<endl;
199 cout <<
"===========================================================================================" << endl;
206 if ((
max(expa_0))(0)<=1.e-5)
209 cout <<
"=============================================================================" << endl;
210 cout <<
" WARNING: Routine failure..." << endl;
211 cout <<
"the horizon expansion is already pretty close to zero..."<< endl;
212 cout <<
"We advise to switch from hyperbolic evolution to parabolic evolution" << endl;
213 cout <<
"This routine will not do anything relevant on parameters alpha, beta, gamma" << endl;
214 cout <<
"=============================================================================" << endl;
215 double is_expansion_adapted = 1.;
216 assert (is_expansion_adapted ==2.);
233 cout <<
"beta: " << beta << endl;
236 gamma = -beta*beta* (1./8.);
238 cout <<
"gamma: " << gamma << endl;
243 double nb_l = (*tfz).set(0).get_dim(1);
244 double nb_m = (*tfz).set(0).get_dim(2);
246 cout <<
"nb_l: " << nb_l << endl;
247 cout <<
"nb_m: " << nb_m << endl;
250 alpha = 2.*fabs(beta - (*tfz).val_in_bound_jk(0,1,0));
253 for (
int ii=0; ii <nb_m; ii++){
254 alpha_aux = 2.*fabs(beta - (*tfz).val_in_bound_jk(0,1,ii));
255 if (alpha_aux >=alpha){
260 cout <<
"alpha: " << alpha << endl;
262 for (
int ii=2; ii <nb_l; ii++)
263 for (
int jj=0; jj <nb_m; jj++){
264 alpha_aux = 2.*(beta - (*tfz).val_in_bound_jk(0,ii,jj))/(
double(ii*(ii+1)));
265 if (alpha_aux >= alpha){
266 cout <<
"higher Ylm harmonics are predominant in the expansion variation on the horizon!" << endl;
267 cout <<
"changing the values of the parameters accordingly..." << endl;
305 expa_3.std_spectral_base();
313 for (
int f= 0; f<nz; f++)
314 for (
int k=0; k<np; k++)
315 for (
int j=0; j<nt; j++) {
316 for (
int l=0; l<nr; l++) {
318 expa_3.set_grid_point(f,k,j,l) = exppa.val_grid_point(0,k,j,0);
325 expa_3.annule_domain(0);
326 expa_3.annule_domain(nz - 1);
328 expa_3.std_spectral_base();
329 expa_3.set_spectral_va().ylm();
335 Metric gamconf(gamconfcov);
342 bound_psi = 0.25*bound_psi;
348 bound_psi = bound_psi + bound_psi2;
375 theta_plus3.std_spectral_base();
380 expa_fin3.std_spectral_base();
388 for (
int f= 0; f<nz; f++)
389 for (
int k=0; k<np; k++)
390 for (
int j=0; j<nt; j++) {
391 for (
int l=0; l<nr; l++) {
393 theta_plus3.set_grid_point(f,k,j,l) = thetaplus.
val_grid_point(0,k,j,0);
394 expa_fin3.set_grid_point(f,k,j,l) = expa_fin.
val_grid_point(0,k,j,0);
399 theta_plus3.annule_domain(0);
400 theta_plus3.annule_domain(nz - 1);
401 expa_fin3.annule_domain(0);
402 expa_fin3.annule_domain(nz - 1);
407 Scalar ff =
lapse*(c_psi_lap*theta_plus3.lapang() + c_psi_fin*(theta_plus3 - expa_fin3));
408 ff.std_spectral_base();
418 Metric gamconf(gamconfcov);
421 theta_int += psi_int*tilde_s.
divergence(gamconf);
423 theta_int = theta_int/
pow(psi_int,3);
428 Scalar ff_int =
lapse*(c_psi_lap*theta_int.lapang() + c_psi_fin*(theta_int - expa_fin3));
463 theta_plus3.std_spectral_base();
468 expa_fin3.std_spectral_base();
476 for (
int f= 0; f<nz; f++)
477 for (
int k=0; k<np; k++)
478 for (
int j=0; j<nt; j++) {
479 for (
int l=0; l<nr; l++) {
481 theta_plus3.set_grid_point(f,k,j,l) = thetaplus.
val_grid_point(0,k,j,0);
482 expa_fin3.set_grid_point(f,k,j,l) = expa_fin.
val_grid_point(0,k,j,0);
487 theta_plus3.annule_domain(0);
488 theta_plus3.annule_domain(nz - 1);
489 expa_fin3.annule_domain(0);
490 expa_fin3.annule_domain(nz - 1);
499 Scalar ff =
lapse*(c_theta_lap*theta_plus3.lapang() + c_theta_fin*(theta_plus3 - expa_fin3));
500 ff.std_spectral_base();
509 Scalar ff_int =
lapse*(c_theta_lap*theta_int.
lapang() + c_theta_fin*(theta_int - expa_fin3));
516 Scalar bound_theta = theta_plus3 + k_2;
548 ff.std_spectral_base() ;
571 k_2.annule_domain(nz-1);
590 const Scalar& Excision_surf::get_BC_lapse_1(
double value)
const{
595 bound_lapse.set_spectral_va().ylm();
609 ff.std_spectral_base();
610 ff = ff -
lapse*c_lapse_fin*lapse_fin;
611 ff.std_spectral_base();
623 ff_int = ff_int -
lapse*c_lapse_fin*lapse_fin;
658 Scalar Matter = 8.*M_PI*(bb*Ee);
661 Matter.annule_domain(nz2-1);
677 for (
int k=0; k<np; k++)
678 for (
int j=0; j<nt; j++) {
696 Scalar source = 2.*
contract(Ll.
up_down(qab),0,
sph.
derive_cov2d(bb2),0) + bb2*(
contract(
sph.
derive_cov2d(Ll.
up_down(qab)),0,1)+
contract(Ll,0,Ll.
up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) + 0.25*Tplus*Tplus + 0.5*
contract(Shear.
up_down(qab),0,1,Shear,0,1)) + Matter2 + Tplus*Kappa2 + dttheta;
700 source = source + lapb;
705 source = source - source_add;
710 qab_con.
set(1,1) = 1. ;
711 qab_con.
set(1,2) = 0. ;
712 qab_con.
set(1,3) = 0. ;
724 Scalar lap_rem =
contract(qab_con, 0,1,
contract(Delta2,0,dN,0),0,1)*sqrt_q_h2 -
sph.
get_hsurf()*
sph.
get_hsurf()*
contract(hab,0,1,ddN,0,1);
726 source = sqrt_q_h2*source + lap_rem;
730 if (sph_sym ==
true){
731 Scalar lapang_s_par = (
contract(
sph.
derive_cov2d(Ll.
up_down(qab)),0,1)+
contract(Ll,0,Ll.
up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) - 0.25*Tplus*Tplus - 0.5*
contract(Shear.
up_down(qab),0,1,Shear,0,1))*sqrt_q_h2;
740 cout <<
"non_spherical case has not been treated yet!" << endl;}
752 for (
int f= 0; f<nz2; f++)
753 for (
int k=0; k<np2; k++)
754 for (
int j=0; j<nt2; j++) {
755 for (
int l=0; l<nr2; l++) {
757 bound_N3.set_grid_point(f,k,j,l) = bound_N.
val_grid_point(0,k,j,0);
763 bound_N3.annule_domain(nz2 - 1);
767 bound_N3.std_spectral_base();
768 bound_N3.set_spectral_va().ylm();
793 Scalar Matter = 8.*M_PI*(bb*Ee);
796 Matter.annule_domain(nz2-1);
810 for (
int k=0; k<np; k++)
811 for (
int j=0; j<nt; j++) {
829 Scalar source = 2.*
contract(Ll.
up_down(qab),0,
sph.
derive_cov2d(nn2 -bb2),0) + (nn2- bb2)*(
contract(
sph.
derive_cov2d(Ll.
up_down(qab)),0,1)+
contract(Ll,0,Ll.
up_down(qab),0)- 0.5*(Ricci2 + Tplus*Tminus)) -(bb2 +nn2)*(0.25*Tplus*Tplus + 0.5*
contract(Shear.
up_down(qab),0,1,Shear,0,1)) - Matter2 - Tplus*Kappa2 ;
833 source = source + lapNmb;
837 cout <<
"mean difference between old and new expansion" << endl;
855 const Vector& Excision_surf::get_BC_shift_1(
double Omega)
const{
873 V_par.
set(3) = V_phi;
879 Vector bound_shift = bound*sss + V_par;
902 ff.std_spectral_base();
919 Scalar bound_b_min_N = b_min_N + k_2;
966 ricci3.std_spectral_base();
974 for (
int f= 0; f<nz; f++)
975 for (
int k=0; k<np; k++)
976 for (
int j=0; j<nt; j++) {
977 for (
int l=0; l<nr; l++) {
986 ricci3.annule_domain(nz - 1);
994 ricci_t = 0.5*ricci3*ricci_t;
1054 Vector V_par_int = V_par + k_1V;
1110 Vector bound_V = V_par + k_2V;
1114 Vector bound_shift = bb2*sss + bound_V;
1144 pure_source = dtpsi - pure_source;
1151 Scalar source = pure_source/factor;
1195 ricci3.std_spectral_base();
1203 for (
int f= 0; f<nz; f++)
1204 for (
int k=0; k<np; k++)
1205 for (
int j=0; j<nt; j++) {
1206 for (
int l=0; l<nr; l++) {
1215 ricci3.annule_domain(nz - 1);
1223 ricci_t = 0.5*ricci3*ricci_t;
1283 Vector V_par_int = V_par + k_1V;
1339 Vector bound_V = V_par + k_2V;
1343 Vector bound_shift = bb2*sss + bound_V;
1379 Scalar Matter = 8.*M_PI*(bb*Ee);
1382 Matter.annule_domain(nz-1);
1398 for (
int k=0; k<np; k++)
1399 for (
int j=0; j<nt; j++) {
1402 Matter2.
set_grid_point(0,k,j,0) = Matter.val_grid_point(1,k,j,0);
1419 Scalar source = 2.*
contract(Ll.
up_down(qab),0,
sph.
derive_cov2d(nn2),0) + nn2*(
contract(
sph.
derive_cov2d(Ll.
up_down(qab)),0,1)+
contract(Ll,0,Ll.
up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) - 0.25*Tplus*Tplus - 0.5*
contract(Shear.
up_down(qab),0,1,Shear,0,1)) - Matter2 - Tplus*Kappa2 - dttheta;
1423 source = source + lapN;
1428 source = source - source_add;
1433 qab_con.
set(1,1) = 1. ;
1434 qab_con.
set(1,2) = 0. ;
1435 qab_con.
set(1,3) = 0. ;
1447 Scalar lap_rem =
contract(qab_con, 0,1,
contract(Delta2,0,db,0),0,1)*sqrt_q_h2 -
sph.
get_hsurf()*
sph.
get_hsurf()*
contract(hab,0,1,ddb,0,1);
1449 source = sqrt_q_h2*source + lap_rem;
1452 Scalar lapang_s_par = (
contract(
sph.
derive_cov2d(Ll.
up_down(qab)),0,1)+
contract(Ll,0,Ll.
up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) + 0.25*Tplus*Tplus + 0.5*
contract(Shear.
up_down(qab),0,1,Shear,0,1))*sqrt_q_h2;
1455 if (sph_sym ==
true){
1456 bound_b = source/lapang_par;
1462 cout <<
"non spherical case has not been treated yet!" << endl;
1475 for (
int f= 0; f<nz; f++)
1476 for (
int k=0; k<np2; k++)
1477 for (
int j=0; j<nt2; j++) {
1478 for (
int l=0; l<nr2; l++) {
1480 bound_b3.set_grid_point(f,k,j,l) = bound_b.
val_grid_point(0,k,j,0);
1486 bound_b3.annule_domain(nz - 1);
1490 bound_b3.std_spectral_base();
1491 bound_b3.set_spectral_va().ylm();
1536 ricci3.std_spectral_base();
1538 for (
int f= 0; f<nz; f++)
1539 for (
int k=0; k<np2; k++)
1540 for (
int j=0; j<nt2; j++) {
1541 for (
int l=0; l<nr2; l++) {
1550 ricci3.annule_domain(nz - 1);
1558 ricci_t = 0.5*ricci3*ricci_t;
1618 Vector V_par_int = V_par + k_1V;
1674 Vector bound_V = V_par + k_2V;
1678 Vector bound_shift = bbb2*sss + bound_V;
1705 const Scalar& Excision_surf::get_BC_Npsi_1(
double value)
const{
1723 ff.std_spectral_base();
1724 ff = ff -
lapse*c_npsi_fin*npsi_fin;
1725 ff.std_spectral_base();
1735 Scalar ff_int =
lapse*(c_npsi_lap*npsi_int.
lapang() + c_npsi_fin*(npsi_int - npsi_fin));
1742 Scalar bound_npsi = npsi + k_2;
1763 boundN.std_spectral_base();
1810 bound_npsi = bound_npsi - Kappa;
1834 bound_npsi += bound_npsi2;
1853 cout <<
"c'est pas fait!" << endl ;
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
const Scalar & theta_plus() const
Computes the outgoing null expansion .
const Scalar & get_BC_Npsi_5(double Kappa) const
Source for a Neumann BC on (N*Psi1), fixing a constant surface gravity in space and time...
Scalar * p_get_BC_lapse_4
Source of Dirichlet condtion on , based on einstein equations (conservation of isotropic gauge) ...
Metric for tensor calculation.
const Scalar & get_BC_conf_fact_2(double c_psi_lap, double c_psi_fin, Scalar &expa_fin) const
Source for the Dirichlet BC on the conformal factor, based on a parabolic driver for the conformal fa...
double delta_t
The time step for evolution in parabolic drivers.
Sym_tensor Kij
The 3-d extrinsic curvature on the slice.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Scalar * p_get_BC_Npsi_2
Source of Dirichlet boundary condition on .
Scalar * p_get_BC_Npsi_4
Source of Dirichlet boundary condition on .
const Scalar & derive_t_expa(Scalar &Ee, Vector &Jj, Sym_tensor &Sij) const
Forms the prospective time derivative for the expansion using projected Einstein equations. Does NOT modify the member dt_expa: do it by hand!
Scalar expa
The 2d expansion, directly evolved from the initial excision with Einstein Equations.
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Vector * p_get_BC_shift_4
Source of Dirichlet BC for the shift vector , partly from projection of Einstein Equations.
Cmp sqrt(const Cmp &)
Square root.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
const Vector & get_BC_shift_3(Scalar &dtpsi, double c_V_lap, double epsilon) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; Radial part is dealt with using ...
Scalar * p_get_BC_conf_fact_2
Source of Neumann boundary condition on ,.
const Scalar & get_BC_conf_fact_3(double c_theta_lap, double c_theta_fin, Scalar &expa_fin) const
Source for the Neumann BC on the conformal factor, based on a parabolic driver for the expansio...
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Scalar & get_BC_Npsi_4(double Kappa) const
Source for a Dirichlet BC on (N*Psi1), fixing a constant surface gravity in space and time...
Flat metric for tensor calculation.
Scalar & set_expa()
Sets the expansion function on the surface at time t (considering to protect this function) ...
Tensor field of valence 0 (or component of a tensorial field).
void operator=(const Excision_surf &)
Assignment to another Excision_surf.
Spheroid sph
The associated Spheroid object.
void coef_i() const
Computes the physical value of *this.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Tensor field of valence 1.
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Scalar & set_dt_expa()
Sets the time derivative of the expansion function on the surface at time t (considering to protect t...
Scalar * p_get_BC_lapse_3
Source of Dirichlet condtion on , based on einstein equations.
Scalar * p_get_BC_Npsi_5
Source of Neumann boundary condition on .
void set_dzpuis(int)
Modifies the dzpuis flag.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
void annule_hard()
Sets the Scalar to zero in a hard way.
Scalar lapse
The lapse defined on the 3 slice.
int get_etat() const
Returns the logical state.
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Excision_surf(const Scalar &h_in, const Metric &gij, const Sym_tensor &Kij2, const Scalar &ppsi, const Scalar &nn, const Vector &beta, double timestep, int int_nos)
Constructor of an excision surface embedded in a 3-slice (Time_slice ) of 3+1 formalism.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Scalar * p_get_BC_conf_fact_4
Source of Birichlet boundary condition on ,.
Scalar * p_get_BC_conf_fact_3
Source of Neumann boundary condition on ,.
virtual void del_deriv() const
Deletes all the derived quantities.
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
void get_evol_params_from_ID(double alpha, double beta, double gamma, Scalar &Ee, Vector &Jj, Sym_tensor &Ss)
Computes the parameters for the hyperbolic evolution in set_expa_hyperb(), so that the expansion has ...
Scalar * p_get_BC_Npsi_3
Source of Dirichlet boundary condition on .
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Scalar * p_derive_t_expa
Computation of an updated expansion scalar.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
virtual void sauve(FILE *) const
Save in a file.
const Scalar & get_BC_Npsi_2(double value, double c_npsi_lap, double c_npsi_fin) const
Source for the Dirichlet BC on (N*Psi1), based on a parabolic driver.
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Vector * p_get_BC_shift_3
Source of Dirichlet BC for the shift vector , partly derived from kinematical relation.
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Active physical coordinates and mapping derivatives.
const Scalar & get_BC_lapse_3(Scalar &dttheta, Scalar &Ee, Vector &Jj, Sym_tensor &Sij, bool sph_sym=true) const
Source for Dirichlet BC on the lapse, based on einstein equations.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Scalar * p_get_BC_lapse_1
Source of Dirichlet boundary condition of .
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
double no_of_steps
The internal number of timesteps for one iteration.
Surface where boundary conditions for quantities in the bulk will be calculated It relies on geometri...
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
virtual ~Excision_surf()
Destructor.
Scalar dt_expa
The time derivative of the expansion, derived from Einstein equations and arbitrary evolution...
Scalar conf_fact
The value of the conformal factor on the 3-slice.
const Vector & get_ll() const
Returns the vector .
const Metric & get_qab() const
Returns the metric .
Coefficients storage for the multi-domain spectral method.
const Scalar & get_BC_Npsi_3(double n_0, double beta) const
Source for the Dirichlet BC on (N*Psi1), with Kerr_Schild-like form for the lapse boundary...
const Scalar & dsdr() const
Returns of *this .
const Vector & get_BC_shift_2(double c_bb_lap, double c_bb_fin, double c_V_lap, double epsilon) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; no assumptions are made except a...
Scalar * p_get_BC_lapse_2
Source of Dirichlet boundary condition of .
const Scalar & get_BC_lapse_2(double lapse_fin, double c_lapse_lap, double c_lapse_fi) const
Source for Dirichlet BC on the lapse, based on a parabolic driver towards arbitrary constant value...
Valeur & set_spectral_va()
Returns va (read/write version)
Scalar & set(int)
Read/write access to a component.
const Scalar & get_BC_conf_fact_4() const
Source for the Dirchlet BC on the conformal factor, based on the consistency condition derived from t...
Vector * p_get_BC_shift_1
Source of Dirichlet BC for the shift vector .
Metric gamij
The 3-d metric on the slice.
Vector * p_get_BC_shift_2
Source of Dirichlet BC for the shift vector .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Map & get_mp() const
Returns the mapping.
Scalar * p_get_BC_Npsi_1
Source of Neumann boundary condition on .
const Vector & get_BC_shift_4(Scalar &dttheta, Scalar &Ee, Vector &Jj, Sym_tensor &Sij, double c_V_lap, double epsilon, bool sph_sym=true) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; Radial part is dealt with using ...
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Vector shift
The Shift 3-vector on the slice.
Class intended to describe valence-2 symmetric tensors.
const Scalar & get_BC_conf_fact_1(bool isMOTS=false) const
Source for a Neumann BC on the conformal factor. If boolean isMOTS is false, it is based on expansion...
const Valeur & get_spectral_va() const
Returns va (read only version)
const Scalar & get_hsurf() const
Returns the field h_surf.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Scalar * p_get_BC_conf_fact_1
Source of Neumann boundary condition on ,.
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.