92 #include "et_bin_nsbh.h" 95 #include "graphique.h" 96 #include "utilitaires.h" 101 int mermax_poisson,
double relax_poisson,
102 int mermax_potvit,
double relax_potvit,
116 Map_et& mp_et =
dynamic_cast<Map_et&
>(
mp) ;
120 double& diff_ent = diff.set(0) ;
121 double& diff_vel_pot = diff.set(1) ;
122 double& diff_lapse = diff.set(2) ;
123 double& diff_confpsi = diff.set(3) ;
124 double& diff_shift_x = diff.set(4) ;
125 double& diff_shift_y = diff.set(5) ;
126 double& diff_shift_z = diff.set(6) ;
131 double precis_poisson = 1.e-16 ;
134 par_poisson1.add_int(mermax_poisson, 0) ;
135 par_poisson1.add_double(relax_poisson, 0) ;
136 par_poisson1.add_double(precis_poisson, 1) ;
137 par_poisson1.add_int_mod(niter, 0) ;
144 par_poisson2.add_int(mermax_poisson, 0) ;
145 par_poisson2.add_double(relax_poisson, 0) ;
146 par_poisson2.add_double(precis_poisson, 1) ;
147 par_poisson2.add_int_mod(niter, 0) ;
153 Param par_poisson_vect ;
154 par_poisson_vect.add_int(mermax_poisson, 0) ;
156 par_poisson_vect.add_double(relax_poisson, 0) ;
157 par_poisson_vect.add_double(precis_poisson, 1) ;
158 par_poisson_vect.add_cmp_mod(
ssjm1_khi ) ;
160 par_poisson_vect.add_int_mod(niter, 0) ;
166 int adapt_flag = (adapt) ? 1 : 0 ;
167 int nz_search =
nzet + 1 ;
168 double precis_secant = 1.e-14 ;
170 double reg_map = 1. ;
173 Tbl ent_limit(
nzet) ;
175 par_adapt.add_int(nitermax, 0) ;
176 par_adapt.add_int(
nzet, 1) ;
177 par_adapt.add_int(nz_search, 2) ;
178 par_adapt.add_int(adapt_flag, 3) ;
179 par_adapt.add_int(j_b, 4) ;
180 par_adapt.add_int(k_b, 5) ;
181 par_adapt.add_int_mod(niter_adapt, 0) ;
182 par_adapt.add_double(precis_secant, 0) ;
183 par_adapt.add_double(reg_map, 1) ;
184 par_adapt.add_double(alpha_r, 2) ;
185 par_adapt.add_tbl(ent_limit, 0) ;
192 Tenseur ent_jm1 =
ent ;
202 for(
int mer=0 ; mer<mermax ; mer++ ) {
204 cout <<
"-----------------------------------------------" << endl ;
205 cout <<
"step: " << mer << endl ;
206 cout <<
"diff_ent = " << diff_ent << endl ;
221 int nt = mg->get_nt(
nzet-1) ;
222 int np = mg->get_np(
nzet-1) ;
223 int nr = mg->get_nr(
nzet-1) ;
226 double hc =
exp(ent_c) ;
227 double gamma_c =
exp(
loggam())(0,0,0,0) ;
229 double n_auto_c =
n_auto()(0,0,0,0) ;
230 double n_comp_c =
n_comp()(0,0,0,0) ;
232 double alpha_square = 0 ;
233 double constante = 0;
234 for (
int k=0; k<np; k++) {
235 for (
int j=0; j<nt; j++) {
244 double alpha_square_courant = (gamma_0_c*gamma_b*n_comp_b - hc*gamma_c*gamma_0_b*n_comp_c) /
245 (hc*gamma_c*gamma_0_b*n_auto_c-gamma_0_c*gamma_b*n_auto_b) ;
246 double constante_courant = gamma_b*(n_comp_b+alpha_square_courant*n_auto_b)/gamma_0_b ;
248 if (alpha_square_courant > alpha_square) {
249 alpha_square = alpha_square_courant ;
252 constante = constante_courant ;
257 alpha_r =
sqrt(alpha_square) ;
258 cout <<
"Adaptation : " << k_b <<
" " << j_b <<
" " << alpha_r << endl ;
262 potentiel.set_std_base() ;
263 for (
int l=
nzet+1 ; l<nz ; l++)
264 potentiel.set().va.set(l) = 1 ;
266 Map_et mp_prev = mp_et ;
272 for (
int l=0; l<
nzet; l++) {
273 ent_limit.set(l) =
ent()(l, k_b, j_b, nr-1) ;
278 mp_prev.homothetie(alpha_r) ;
280 for (
int l=
nzet ; l<nz-1 ; l++)
302 Tbl diff_ent_tbl =
diffrel(
ent(), ent_jm1() ) ;
303 diff_ent = diff_ent_tbl(0) ;
304 for (
int l=1; l<
nzet; l++) {
305 diff_ent += diff_ent_tbl(l) ;
328 tmp2.set_etat_qcq() ;
329 for (
int i=0 ; i<3 ; i++) {
330 tmp2.set() = tmp(i, i) ;
335 +
nnn * confpsi_q * kk
341 "WARNING : Et_bin_nsbh is for the relativistic calculation" 346 source.set_std_base() ;
350 Cmp n_auto_old (
n_auto()) ;
351 source().poisson(par_poisson1,
n_auto.
set()) ;
359 "Relative difference on n_auto : " 361 for (
int l=0; l<nz; l++) {
362 cout << tdiff_lapse(l) <<
" " ;
365 diff_lapse =
max(
abs(tdiff_lapse)) ;
382 tmp2.set_etat_qcq() ;
383 for (
int i=0 ; i<3 ; i++) {
384 tmp2.set() = tmp(i, i) ;
389 - 0.125 * confpsi_c * kk ;
405 "Relative difference on confpsi_auto : " 407 for (
int l=0; l<nz; l++) {
408 cout << tdiff_confpsi(l) <<
" " ;
411 diff_confpsi =
max(
abs(tdiff_confpsi)) ;
429 for (
int i=0 ; i<3 ; i++)
430 if (source_shift(i).get_etat() != ETATZERO)
431 source_shift.set(i).va.coef_i() ;
433 for (
int i=0; i<3; i++)
434 if ((source_shift(i).get_etat() != ETATZERO) && (source_shift(i).va.c->t[nz-1]->get_etat() != ETATZERO))
435 source_shift.set(i).filtre(4) ;
436 for (
int i=0; i<3; i++) {
437 if(source_shift(i).dz_nonzero()) {
438 assert( source_shift(i).get_dzpuis() == 4 ) ;
441 (source_shift.set(i)).set_dzpuis(4) ;
447 double lambda_shift = double(1) / double(3) ;
451 source_shift.poisson_vect_oohara(lambda_shift, par_poisson_vect,
465 "Relative difference on shift_auto : " 467 cout <<
"x component : " ;
468 for (
int l=0; l<nz; l++) {
469 cout << tdiff_shift_x(l) <<
" " ;
472 cout <<
"y component : " ;
473 for (
int l=0; l<nz; l++) {
474 cout << tdiff_shift_y(l) <<
" " ;
477 cout <<
"z component : " ;
478 for (
int l=0; l<nz; l++) {
479 cout << tdiff_shift_z(l) <<
" " ;
483 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
484 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
485 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
498 int,
double,
double,
const Tbl&,
Tbl&) {
500 cout <<
"Not implemented !" << endl ;
Cmp log(const Cmp &)
Neperian logarithm.
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Cmp exp(const Cmp &)
Exponential.
Cmp sqrt(const Cmp &)
Square root.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Tenseur pot_centri
Centrifugal potential.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tenseur nnn
Total lapse function.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
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...
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Tenseur press
Fluid pressure.
Tenseur n_auto
Part of the lapse { N} generated principaly by the star.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tenseur confpsi
Total conformal factor $$.
Tenseur d_confpsi_comp
Gradient of { confpsi_comp} (Cartesian components with respect to { ref_triad})
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
bool irrotational
true for an irrotational star, false for a corotating one
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Map & mp
Mapping associated with the star.
int get_nzone() const
Returns the number of domains.
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
int get_etat() const
Returns the logical state.
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog...
Cmp pow(const Cmp &, int)
Power .
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Tenseur d_confpsi_auto
Gradient of { confpsi_auto} (Cartesian components with respect to { ref_triad})
Tenseur d_n_auto
Gradient of { n_auto} (Cartesian components with respect to { ref_triad})
int nzet
Number of domains of *mp occupied by the star.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Tenseur confpsi_auto
Part of the conformal factor $$ generated principaly by the star.
virtual void reevaluate_symy(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Cmp ssjm1_confpsi
Effective source at the previous step for the resolution of the Poisson equation for { confpsi_auto} ...
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Cmp abs(const Cmp &)
Absolute value.
Cmp ssjm1_lapse
Effective source at the previous step for the resolution of the Poisson equation for { n_auto} by mea...
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Tenseur_sym tkij_tot
Total extrinsic curvature tensor $K^{ij}$ generated by { shift_auto} and { shift_comp}.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor $K^{ij}$ generated by { shift_auto}.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur n_comp
Part of the lapse { N} generated principaly by the companion star.
virtual void equilibrium_nsbh(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration in a NS-BH binary system.
Valeur va
The numerical value of the Cmp.
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...