118 #include "isol_hor.h" 120 #include "utilitaires.h" 130 assert ((relax >0) && (relax<=1)) ;
132 cout <<
"-----------------------------------------------" << endl ;
133 cout <<
"Resolution LAPSE" << endl ;
181 source_un += tmp_un ;
217 source_deux += tmp_deux ;
219 cout <<
"source lapse" << endl <<
norme(source_un) << endl ;
234 lim_un =
hole1.boundary_nn_Dir(lim_nn) ;
235 lim_deux =
hole2.boundary_nn_Dir(lim_nn) ;
237 n_un_temp = n_un_temp - 1./2. ;
238 n_deux_temp = n_deux_temp - 1./2. ;
240 dirichlet_binaire (source_un, source_deux, lim_un, lim_deux,
241 n_un_temp, n_deux_temp, 0, precision) ;
247 lim_un =
hole1.boundary_nn_Neu(lim_nn) ;
248 lim_deux =
hole2.boundary_nn_Neu(lim_nn) ;
250 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
251 n_un_temp, n_deux_temp, 0, precision) ;
256 cout <<
"Unexpected type of boundary conditions for the lapse!" 258 <<
" bound_nn = " << bound_nn << endl ;
265 n_un_temp = n_un_temp + 1./2. ;
266 n_deux_temp = n_deux_temp + 1./2. ;
275 cout <<
"lapse auto" << endl <<
norme (n_un_temp) << endl ;
276 Tbl tdiff_nn =
diffrel(n_un_temp.laplacian(), source_un) ;
279 "Relative error in the resolution of the equation for the lapse : " 281 for (
int l=0; l<nz; l++) {
282 cout << tdiff_nn(l) <<
" " ;
289 n_un_temp = relax*n_un_temp + (1-relax)*lapse_un_old ;
290 n_deux_temp = relax*n_deux_temp + (1-relax)*lapse_deux_old ;
306 assert ((relax>0) && (relax<=1)) ;
308 cout <<
"-----------------------------------------------" << endl ;
309 cout <<
"Resolution PSI" << endl ;
379 cout <<
"source psi" << endl <<
norme(source_un) << endl ;
397 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
398 psi_un_temp, psi_deux_temp, 0, precision) ;
403 cout <<
"Unexpected type of boundary conditions for psi!" 405 <<
" bound_psi = " << bound_psi << endl ;
412 psi_un_temp = psi_un_temp + 1./2. ;
413 psi_deux_temp = psi_deux_temp + 1./2. ;
422 cout <<
"psi auto" << endl <<
norme (psi_un_temp) << endl ;
423 Tbl tdiff_psi =
diffrel(psi_un_temp.laplacian(), source_un) ;
426 "Relative error in the resolution of the equation for psi : " 428 for (
int l=0; l<nz; l++) {
429 cout << tdiff_psi(l) <<
" " ;
436 psi_un_temp = relax*psi_un_temp + (1-relax)*psi_un_old ;
437 psi_deux_temp = relax*psi_deux_temp + (1-relax)*psi_deux_old ;
459 cout <<
"------------------------------------------------" << endl ;
460 cout <<
"Resolution shift : Omega = " <<
omega << endl ;
494 source_un += 2.*
hole1.
nn * ( tmp_vect_un
507 source_un -= vtmp_un ;
537 source_deux += 2.*
hole2.
nn * ( tmp_vect_deux
550 source_deux -= vtmp_deux ;
558 Vector source_1 (source_un) ;
559 Vector source_2 (source_deux) ;
563 cout <<
"source shift_x" << endl <<
norme(source_1(1)) << endl ;
564 cout <<
"source shift_y" << endl <<
norme(source_1(2)) << endl ;
565 cout <<
"source shift_z" << endl <<
norme(source_1(3)) << endl ;
568 for (
int i=1 ; i<=3 ; i++) {
569 source_un.
set(i).filtre(4) ;
570 source_deux.set(i).filtre(4) ;
584 switch (bound_beta) {
588 lim_x_un =
hole1.boundary_beta_x(
omega, omega_eff) ;
589 lim_y_un =
hole1.boundary_beta_y(
omega, omega_eff) ;
590 lim_z_un =
hole1.boundary_beta_z() ;
592 lim_x_deux =
hole2.boundary_beta_x(
omega, omega_eff) ;
593 lim_y_deux =
hole2.boundary_beta_y(
omega, omega_eff) ;
594 lim_z_deux =
hole2.boundary_beta_z() ;
599 cout <<
"Unexpected type of boundary conditions for beta!" 601 <<
" bound_beta = " << bound_beta << endl ;
617 poisson_vect_binaire (1./3., source_un, source_deux,
618 lim_x_un, lim_y_un, lim_z_un,
619 lim_x_deux, lim_y_deux, lim_z_deux,
620 beta1, beta2, 0, precision) ;
626 for (
int i=1 ; i<=3 ; i++) {
631 cout <<
"shift_auto x" << endl <<
norme(beta1(1)) << endl ;
632 cout <<
"shift_auto y" << endl <<
norme(beta1(2)) << endl ;
633 cout <<
"shift_auto z" << endl <<
norme(beta1(3)) << endl ;
646 Tbl tdiff_beta_r =
diffrel(lap_beta(1), source_un(1)) ;
647 Tbl tdiff_beta_t =
diffrel(lap_beta(2), source_un(2)) ;
648 Tbl tdiff_beta_p =
diffrel(lap_beta(3), source_un(3)) ;
651 "Relative error in the resolution of the equation for beta : " 653 cout <<
"r component : " ;
654 for (
int l=0; l<nz; l++) {
655 cout << tdiff_beta_r(l) <<
" " ;
658 cout <<
"t component : " ;
659 for (
int l=0; l<nz; l++) {
660 cout << tdiff_beta_t(l) <<
" " ;
663 cout <<
"p component : " ;
664 for (
int l=0; l<nz; l++) {
665 cout << tdiff_beta_p(l) <<
" " ;
688 omdsdp1.set(2) =
omega * xxa1 ;
692 omdsdp1.set(1) =
omega * yya1 ;
693 omdsdp1.set(2) = -
omega * xxa1 ;
697 omdsdp1.set(1).set_spectral_va()
699 omdsdp1.set(2).set_spectral_va()
701 omdsdp1.set(3).set_spectral_va()
704 omdsdp1.annule_domain(nz-1) ;
716 omdsdp2.set(2) =
omega * xxa2 ;
720 omdsdp2.set(1) =
omega * yya2 ;
721 omdsdp2.set(2) = -
omega * xxa2 ;
725 omdsdp2.set(1).set_spectral_va()
727 omdsdp2.set(2).set_spectral_va()
729 omdsdp2.set(3).set_spectral_va()
732 omdsdp2.annule_domain(nz-1) ;
736 beta1_new = relax*(beta1+
hole1.
decouple*omdsdp1) + (1-relax)*beta_un_old ;
737 beta2_new = relax*(beta2+
hole2.
decouple*omdsdp2) + (1-relax)*beta_deux_old ;
753 for (
int k=0; k<nnp; k++)
754 for (
int j=0; j<nnt; j++){
Coord xa
Absolute x coordinate.
void beta_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
double omega
Angular velocity.
Vector dn
Covariant derivative of the lapse with respect to the flat metric .
Metric tgam
3 metric tilde
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Scalar decouple
Function used to construct from the total .
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Scalar n_comp
Lapse function .
double regularisation(const Vector &shift_auto, const Vector &shift_comp, double ang_vel)
Corrects shift_auto in such a way that the total is equal to zero in the horizon, which should ensure the regularity of .
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Tensor field of valence 0 (or component of a tensorial field).
Single_hor hole1
Black hole one.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Scalar & get_psi4() const
Conformal factor .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Values and coefficients of a (real-value) function.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Sym_tensor aa_auto
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Tensor field of valence 1.
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
void annule_hard()
Sets the Scalar to zero in a hard way.
Sym_tensor aa
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for.
Scalar psi
Conformal factor .
Scalar trK
Trace of the extrinsic curvature.
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...
double regul
Intensity of the correction on the shift vector.
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
Sym_tensor hh
Deviation metric.
virtual const Connection & connect() const
Returns the connection.
Vector beta_auto
Shift function .
int get_nzone() const
Returns the number of domains.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Scalar n_auto
Lapse function .
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 .
Single_hor hole2
Black hole two.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Vector dpsi
Covariant derivative of the conformal factor .
Coord ya
Absolute y coordinate.
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ...
Map_af & mp
Affine mapping.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Scalar psi_auto
Conformal factor .
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
Metric_flat ff
3 metric flat
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Class intended to describe valence-2 symmetric tensors.
Scalar nn
Lapse function .