140 #include "utilitaires.h" 148 double relax_poisson,
int mermax_potvit,
149 double relax_potvit,
double thres_adapt,
150 const Tbl& fact_resize,
Tbl& diff) {
169 int i_b = mg->
get_nr(l_b) - 1 ;
179 double& diff_ent = diff.
set(0) ;
180 double& diff_vel_pot = diff.
set(1) ;
181 double& diff_logn = diff.
set(2) ;
182 double& diff_beta = diff.
set(3) ;
183 double& diff_shift_x = diff.
set(4) ;
184 double& diff_shift_y = diff.
set(5) ;
185 double& diff_shift_z = diff.
set(6) ;
197 int nz_search =
nzet ;
200 double precis_secant = 1.e-14 ;
202 double reg_map = 1. ;
206 par_adapt.
add_int(nitermax, 0) ;
211 par_adapt.
add_int(nz_search, 2) ;
213 par_adapt.
add_int(adapt_flag, 3) ;
229 par_adapt.
add_tbl(ent_limit, 0) ;
235 double precis_poisson = 1.e-16 ;
239 par_poisson1.
add_int(mermax_poisson, 0) ;
251 par_poisson2.
add_int(mermax_poisson, 0) ;
261 Param par_poisson_vect ;
263 par_poisson_vect.
add_int(mermax_poisson, 0) ;
265 par_poisson_vect.
add_double(relax_poisson, 0) ;
266 par_poisson_vect.
add_double(precis_poisson, 1) ;
288 Cmp source_regu(
mp) ;
295 for(
int mer=0 ; mer<mermax ; mer++ ) {
297 cout <<
"-----------------------------------------------" << endl ;
298 cout <<
"step: " << mer << endl ;
299 cout <<
"diff_ent = " << diff_ent << endl ;
321 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
322 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
328 double alpha_r2 = 0 ;
329 for (
int k=0; k<mg->
get_np(l_b); k++) {
330 for (
int j=0; j<mg->
get_nt(l_b); j++) {
332 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
333 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
335 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
336 ( logn_auto_b - logn_auto_c ) ;
341 if (alpha_r2_jk > alpha_r2) {
342 alpha_r2 = alpha_r2_jk ;
350 alpha_r =
sqrt(alpha_r2) ;
352 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" " 375 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
386 double dent_eq =
ent().dsdr().val_point(
ray_eq(),M_PI/2.,0.) ;
387 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
388 double rap_dent = fabs( dent_eq / dent_pole ) ;
389 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
391 if ( rap_dent < thres_adapt ) {
393 cout <<
"******* FROZEN MAPPING *********" << endl ;
402 for (
int l=0; l<
nzet; l++) {
403 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
405 ent_limit.
set(
nzet-1) = ent_b ;
418 double rr_in_1 =
mp.
val_r(
nzet, -1., M_PI/2., 0.) ;
421 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
423 mp.
resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
426 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
428 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
433 double rr_out_nm3_new =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
435 for (
int i=
nzet-1; i<nz-4; i++) {
437 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
439 double rr_mid = rr_out_i
440 + (rr_out_nm3_new - rr_out_i) /
double(nz - 3 - i) ;
442 double rr_2timesi = 2. * rr_out_i ;
444 if (rr_2timesi < rr_mid) {
446 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
448 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
453 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
455 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
476 double ent_s_max = -1 ;
479 for (
int k=0; k<mg->
get_np(l_b); k++) {
480 for (
int j=0; j<mg->
get_nt(l_b); j++) {
481 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
482 if (xx > ent_s_max) {
489 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1" 490 <<
" and nzet : " << endl ;
491 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max <<
492 " and j = " << j_s_max << endl ;
529 source = qpig *
nbar ;
547 source().poisson_regular(
k_div,
nzet, unsgam1, par_poisson1,
551 source_regu, source_div) ;
558 "Relative error in the resolution of the equation for logn_auto : " 560 for (
int l=0; l<nz; l++) {
561 cout << tdiff_logn(l) <<
" " ;
564 diff_logn =
max(
abs(tdiff_logn)) ;
596 "Relative error in the resolution of the equation for beta_auto : " 598 for (
int l=0; l<nz; l++) {
599 cout << tdiff_beta(l) <<
" " ;
602 diff_beta =
max(
abs(tdiff_beta)) ;
625 for (
int i=0; i<3; i++) {
627 if (source_shift(i).get_etat() != ETATZERO)
638 for (
int i=0; i<3; i++) {
639 if(source_shift(i).dz_nonzero()) {
640 assert( source_shift(i).get_dzpuis() == 4 ) ;
643 (source_shift.
set(i)).set_dzpuis(4) ;
650 double lambda_shift = double(1) / double(3) ;
653 lambda_shift, par_poisson_vect,
671 for (
int i=0; i<3; i++) {
673 + lambda_shift * graddivn(i) ;
676 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
677 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
678 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
681 "Relative error in the resolution of the equation " 684 cout <<
"x component : " ;
685 for (
int l=0; l<nz; l++) {
686 cout << tdiff_shift_x(l) <<
" " ;
689 cout <<
"y component : " ;
690 for (
int l=0; l<nz; l++) {
691 cout << tdiff_shift_y(l) <<
" " ;
694 cout <<
"z component : " ;
695 for (
int l=0; l<nz; l++) {
696 cout << tdiff_shift_z(l) <<
" " ;
700 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
701 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
702 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
721 diff_ent = diff_ent_tbl(0) ;
722 for (
int l=1; l<
nzet; l++) {
723 diff_ent += diff_ent_tbl(l) ;
int k_div
Index of regularity of the gravitational potential logn_auto .
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void dec2_dzpuis()
dzpuis -= 2 ;
Radial mapping of rather general form.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy and extra parameters (virtual function imp...
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.
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star...
Standard units of space, time and mass.
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)
Tenseur nnn
Total lapse function.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
double ray_eq() const
Coordinate radius at , [r_unit].
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 d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
void inc2_dzpuis()
dzpuis += 2 ;
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by t...
bool irrotational
true for an irrotational star, false for a corotating one
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog...
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Tenseur d_beta_comp
Gradient of beta_comp (Cartesian components with respect to ref_triad )
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Tenseur nbar
Baryon density in the fluid frame.
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Tenseur akcar_comp
Part of the scalar generated by shift_auto and shift_comp , i.e.
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
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.
const Eos & eos
Equation of state of the stellar matter.
virtual void homothetie(double lambda)
Sets a new radial scale.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog...
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int nzet
Number of domains of *mp occupied by the star.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tenseur a_car
Total conformal factor .
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
double ray_pole() const
Coordinate radius at [r_unit].
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
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.
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
void equil_regular(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 by regularizing the diverging source.
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.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Tenseur ener_euler
Total energy density in the Eulerian frame.
void poisson_vect_regu(int k_div, int nzet, double unsgam1, double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy with extra parameters (virtual function im...
Tensor handling *** DEPRECATED : use class Tensor instead ***.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
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...