39#include "excision_surf.h"
40#include "excision_hor.h"
41#include "param_elliptic.h"
58 cout <<
"================================================================" << endl;
59 cout <<
"STARTING THE MAIN ITERATION FOR COMPUTING METRIC FIELDS" << endl;
60 cout <<
" iteration parameters are the following: " << endl;
61 cout <<
" convergence precision required:" << precis << endl;
62 cout <<
" max number of global steps :" << mer_max << endl;
63 cout <<
" relaxation parameter :" << relax << endl;
64 cout <<
"================================================================" << endl;
70 const Mg3d* mgrid = (*map).get_mg();
76 int nt = (*mgrid).get_nt(1);
77 const int np = (*mgrid).get_np(1);
78 const Coord& rr = (*map).r;
87 const Map_af map_2(*g_angu, r_limits2);
100 Scalar logpsi(*map) ; logpsi =
log(conf_fact) ;
102 Scalar psi4 (*map) ; psi4 = conf_fact*conf_fact*conf_fact*conf_fact ;
109 Vector shift_new (*map, CON, (*map).get_bvect_spher());
110 for(
int i=1; i<=3; i++){
126 Sym_tensor aa(*map, CON, (*map).get_bvect_spher());
127 for (
int iii= 1; iii<=3; iii++){
128 for(
int j=1; j<=3; j++){
133 Scalar aa_quad_scal(*map) ; aa_quad_scal = 0. ;
135 Sym_tensor aa_hat(*map, CON, (*map).get_bvect_spher());
136 for (
int iii= 1; iii<=3; iii++){
137 for(
int j=1; j<=3; j++){
138 aa_hat.
set(iii,j)= 0;
152 for (
int i=1; i<=3; i++) {
153 for (
int j=1; j<=3; j++) {
154 for (
int k=1; k<=3; k++) {
157 for (
int l=1; l<=3; l++) {
160 delta.
set(k,i,j) += tmp ;
177 if (isvoid ==
false){
179 cout <<
"FAIL: case of non-void spacetime not treated yet" << endl;
188 double diff_ent = 1 ;
197 for(
int mer=0 ;(diff_ent > precis) && (mer<mer_max) ; mer++) {
233 ssalt = ssalt/ssnormalt;
234 ssaltcon = ssaltcon/ssnormalt;
235 Vector ssconalt = ssaltcon*conf_fact*conf_fact;
241 bound3bis += -conf_fact*ssconalt.
divergence(gamt);
243 bound3bis = 0.25*bound3bis;
244 bound3bis += -
contract(conf_fact.derive_cov(gamt),0,ssconalt,0) + conf_fact.dsdr();
258 Scalar source_conf_fact(*map) ; source_conf_fact=3. ;
261 Scalar d2logpsi =
contract(conf_fact.derive_cov(mets).derive_cov(mets), 0, 1,
hij, 0,1);
264 source_conf_fact = -(0.125* aa_quad_scal )/(psi4*conf_fact*conf_fact*conf_fact) + conf_fact* 0.125* Rstar - d2logpsi;
268 if (source_conf_fact.
get_etat() == ETATZERO) {
289 Scalar baba2 = (conf_fact_new-1).laplacian();
295 psinewbis = psinewbis.
dsdr();
296 Scalar psinewfin2 (map_2) ;
301 for (
int k=0; k<np; k++)
302 for (
int j=0; j<nt; j++) {
315 conf_fact = conf_fact_new* (1-relax) + conf_fact* relax ;
316 psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
317 logpsi =
log(conf_fact) ;
354 Scalar source_npsi = npsi*(aa_quad_scal*(7./8.)/(psi4*psi4) + Rstar/8.) - d2lognpsi;
356 if (source_npsi.
get_etat() == ETATZERO) {
369 npsi_new = npsi_new +1;
381 Scalar npsibound2 (map_2) ;
385 for (
int k=0; k<np; k++)
386 for (
int j=0; j<nt; j++) {
400 npsi = npsi_new*(1-relax) + npsi* relax;
401 lapse = npsi/conf_fact;
417 bound = (
boundNoK)/(conf_fact*conf_fact) ;
426 Vector ephi(*map, CON, (*map).get_bvect_spher());
433 limit = bound*ssconalt + hor_rot*ephi;
451 Vector sourcevect2(*map,CON, (*map).get_bvect_spher());
462 double lam = (1./3.);
469 sourcevect2.
poisson_boundary2(lam, shift_new, Vrb, etab, mmub, 1., 0., 1. ,0. ,1. ,0.) ;
488 for (
int k=0; k<np; k++)
489 for (
int j=0; j<nt; j++) {
498 Scalar brfin = shift_new(1);
506 for (
int k=0; k<np; k++)
507 for (
int j=0; j<nt; j++) {
519 for (
int ii=1; ii <=3; ii++){
524 diff_ent =
max(
maxabs(npsi_new - npsi ));
540 if (diff_ent <=5.e-3) {
551 double precis2 = 1.e5*precis ;
553 double diff_ent2 = 1 ;
559 for(
int mer2=0 ;(diff_ent2 > precis2) && (mer2<mer_max2) ; mer2++) {
576 hij_new = boundfree_tensBC (sourcehij,
shift , conf_fact,
lapse,
hij, precis2);
578 cout <<
"maximum of convergence marker for the subiteration" << endl;
581 hij = relax2*hij_new + (1 - relax2)*
hij;
583 cout <<
"mer2, diffent2" << endl;
585 cout << mer2 << endl;
586 cout << diff_ent2 << endl;
606 Sym_tensor youps = test - sourcehij/((
lapse/(conf_fact*conf_fact))*(
lapse/(conf_fact*conf_fact)));
631 gamtuu = mets.
con() +
hij;
633 gam = gamt.
cov()*psi4;
637 for (
int i=1; i<=3; i++) {
638 for (
int j=1; j<=3; j++) {
649 aa_hat = aa*psi4*
sqrt(psi4);
650 aa_hat.std_spectral_base();
662 for (
int i=1; i<=3; i++) {
663 for (
int j=1; j<=3; j++) {
664 for (
int k=1; k<=3; k++) {
667 for (
int l=1; l<=3; l++) {
670 delta.
set(k,i,j) += tmp ;
689 cout <<
"diffent" << endl;
690 cout<< diff_ent << endl;
691 cout <<
"mer" << mer << endl;
715 Sym_tensor gamt2(*map, COV, (*map).get_bvect_spher());
716 for (
int i=1; i<=3; i++)
717 for (
int j=1; j<=3; j++)
718 { gamt2.
set(i,j)=gam.
cov()(i,j);
732 TrK3 = k_uu.
trace(gam);
739 ham_constr += TrK3*TrK3 -
contract(k_uu, 0, 1, k_dd, 0, 1) ;
764 evol_eq += lapse2*(TrK3*k_dd - 2*
contract(k_dd, 1, k_dd.
up(0, gam), 0) ) ;
778 cout <<
"================================================================" << endl;
779 cout <<
" THE ITERATION HAS NOW CONVERGED" << endl;
780 cout <<
"================================================================" << endl;
Active physical coordinates and mapping derivatives.
Scalar boundNoK
Indicates if the boundary value for the lapse or the surface gravity is used.
Scalar lapse
Lapse function.
bool NorKappa
Rotation rate of the horizon in the azimuthal direction.
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
void secmembre_kerr(Sym_tensor &source_hh)
Computes the rhs of hyperbolic equation for conformal metric assuming statioarity; WARNING; up to now...
Vector shift
Shift vector.
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: see Cordero et al.
const Map & mp
Mapping associated with the star.
bool isCF
Stores the boundary value of the lapse or surface gravity.
void compute_stat_metric(double precis, double relax, int mer_max, int mer_max2, bool isvoid=true)
Computes a quasi-stationary 3-slice from the chosen parameters.
Flat metric for tensor calculation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Metric for tensor calculation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
virtual const Scalar & determinant() const
Returns the determinant.
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
int get_nzone() const
Returns the number of domains.
Coefficients storage for the multi-domain spectral method.
This class contains the parameters needed to call the general elliptic solver.
Tensor field of valence 0 (or component of a tensorial field).
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
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 inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
void annule_hard()
Sets the Scalar to zero in a hard way.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
Scalar sol_elliptic_boundary(Param_elliptic ¶ms, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Class intended to describe valence-2 symmetric tensors.
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
void ylm()
Computes the coefficients of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
Scalar & set(int)
Read/write access to a component.
Cmp sqrt(const Cmp &)
Square root.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp log(const Cmp &)
Neperian logarithm.
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Standard units of space, time and mass.