85 #include "utilitaires.h" 86 #include "et_bin_nsbh.h" 87 #include "graphique.h" 94 assert ((relax>0) && (relax<=1)) ;
96 cout <<
"Resolution LAPSE" << endl ;
105 for (
int i=0 ; i<3 ; i++) {
106 work.
set() = auxi(i, i) ;
115 +psiq*
n_tot()*kk()) ;
123 limite = -0.5 + lim_nn ;
126 for (
int k=0 ; k<np ; k++)
127 for (
int j=0 ; j<nt ; j++)
128 limite.set(0,k,j,0) -=
n_comp() (1, k, j, 0) ;
129 limite.std_base_scal() ;
134 assert(bound_nn == 1);
140 for (
int k=0 ; k<np ; k++)
141 for (
int j=0 ; j<nt ; j++)
142 limite.set(0,k,j,0) -=
n_tot()(1, k, j, 0)/
psi_tot()(1,k,j,0)*
144 limite.std_base_scal() ;
149 soluce = soluce + 0.5 ;
151 n_auto.
set() = relax*soluce + (1-relax)*lapse_old ;
158 assert ((relax>0) && (relax<=1)) ;
160 cout <<
"Resolution PSI" << endl ;
168 for (
int i=0 ; i<3 ; i++) {
169 work.
set() = auxi(i, i) ;
176 Cmp source (-psic*kk()/8.) ;
187 double* vec_s =
new double[3] ;
193 for (
int k=0 ; k<np ; k++)
194 for (
int j=0 ; j<nt ; j++) {
196 double tet = tet_mtbl(1,k,j,0) ;
197 double phi = phi_mtbl(1,k,j,0) ;
198 vec_s[0] =
cos(phi)*
sin(tet) ;
199 vec_s[1] =
sin(phi)*
sin(tet) ;
200 vec_s[2] =
cos(tet) ;
203 for (
int m=0 ; m<3 ; m++)
204 for (
int n=0 ; n<3 ; n++)
205 part_ss += vec_s[m]*vec_s[n]*
tkij_tot(m,n)(1,k,j,0) ;
210 psi_comp().dsdr()(1, k, j, 0) - part_ss ;
214 limite.std_base_scal() ;
217 soluce = soluce + 1./2. ;
219 psi_auto.
set() = relax*soluce + (1-relax)*psi_old ;
226 double precision,
double relax,
227 int bound_nn,
double lim_nn) {
229 assert (precision > 0) ;
230 assert ((relax>0) && (relax<=1)) ;
232 cout <<
"resolution SHIFT" << endl ;
243 for (
int i=0 ; i<3 ; i++)
244 if (source(i).
get_etat() == ETATQCQ)
247 for (
int i=0 ; i<3 ; i++)
254 for (
int i=0 ; i<3 ; i++)
269 double air, theta, phi, xabs, yabs, zabs ;
292 for (
int k=0 ; k<np ; k++)
293 for (
int j=0 ; j<nt ; j++) {
295 double tet = tet_mtbl(1,k,j,0) ;
296 double phy = phi_mtbl(1,k,j,0) ;
298 xabs = Xabs (1, k, j, 0) ;
299 yabs = Yabs (1, k, j, 0) ;
300 zabs = Zabs (1, k, j, 0) ;
304 lim_x.set(0, k, j, 0) =
omega*Yabs(0, 0, 0, 0) +
309 lim_x.base = *bases[0] ;
312 lim_y.set(0, k, j, 0) = -
omega*Xabs(0, 0, 0, 0) -
318 lim_z.set(0, k, j, 0) = -
323 lim_x.base = *bases[0] ;
324 lim_y.base = *bases[1] ;
325 lim_z.base = *bases[2] ;
328 for (
int i=0 ; i<3 ; i++)
333 poisson_vect_frontiere(1./3., source,
shift_auto, lim_x, lim_y,
334 lim_z, 0, precision, 20) ;
338 for (
int i=0; i<3; i++)
344 if (bound_nn == 0 && lim_nn == 0)
353 double precision,
double relax,
354 int bound_nn,
double lim_nn) {
369 void Bhole::update_metric (
const Et_bin_nsbh& comp) {
Coord xa
Absolute x coordinate.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Map & get_mp() const
Returns the mapping.
Tenseur grad_n_tot
Total gradient of N .
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void annule_hard()
Sets the Valeur to zero in a hard way.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
double omega_local
local angular velocity
Tenseur psi_auto
Part of generated by the hole.
Values and coefficients of a (real-value) function.
Tenseur n_comp
Part of N generated by the companion hole.
Tenseur_sym taij_tot
Total , which must be zero on the horizon of the regularisation on the shift has been done...
void fait_n_comp(const Bhole &comp)
Imports the part of N due to the companion hole comp .
Cmp cos(const Cmp &)
Cosine.
Coord tet
coordinate centered on the grid
Tenseur shift_auto
Part of generated by the hole.
Coord phi
coordinate centered on the grid
void solve_psi_with_ns(double relax)
Solves the equation for ~: with the condition that on the horizon, where f is the value of on th...
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tenseur psi_comp
Part of generated by the companion hole.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Map_af & mp
Affine mapping.
Cmp poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Cmp::poisson() .
int get_etat() const
Returns the logical state.
double rayon
Radius of the horizon in LORENE's units.
Cmp pow(const Cmp &, int)
Power .
Tenseur_sym tkij_auto
Auto .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Tenseur_sym tkij_tot
Total .
Coord ya
Absolute y coordinate.
Bases of the spectral expansions.
Tbl & set(int l)
Read/write of the value in a given domain.
Class for a star in a NS-BH binary system.
Coord y
y coordinate centered on the grid
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
void solve_lapse_with_ns(double relax, int bound_nn, double lim_nn)
Solves the equation for N ~: with the condition that N =0 on the horizon.
Coord za
Absolute z coordinate.
void solve_shift_with_ns(const Et_bin_nsbh &ns, double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for ~: with on the horizon, being the boost and .
Coord x
x coordinate centered on the grid
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
double omega
Angular velocity in LORENE's units.
Cmp sin(const Cmp &)
Sine.
Tenseur n_auto
Part of N generated by the hole.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void fait_psi_comp(const Bhole &comp)
Imports the part of due to the companion hole comp .
const Tenseur & get_shift_auto() const
Returns the part of the shift vector generated principaly by the star.
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X...
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Tensor handling *** DEPRECATED : use class Tensor instead ***.
double regul
Intensity of the correction on the shift vector.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)