182 #include "graphique.h" 183 #include "utilitaires.h" 188 int mermax,
int mermax_poisson,
189 double relax_poisson,
int mermax_potvit,
190 double relax_potvit,
double thres_adapt,
191 const Tbl& fact_resize,
Tbl& diff,
const Tbl* pent_limit ) {
210 int i_b = mg->
get_nr(l_b) - 1 ;
220 double& diff_ent = diff.
set(0) ;
221 double& diff_vel_pot = diff.
set(1) ;
222 double& diff_logn = diff.
set(2) ;
223 double& diff_beta = diff.
set(3) ;
224 double& diff_shift_x = diff.
set(4) ;
225 double& diff_shift_y = diff.
set(5) ;
226 double& diff_shift_z = diff.
set(6) ;
238 int nz_search =
nzet ;
241 double precis_secant = 1.e-14 ;
243 double reg_map = 1. ;
245 par_adapt.
add_int(nitermax, 0) ;
250 par_adapt.
add_int(nz_search, 2) ;
252 par_adapt.
add_int(adapt_flag, 3) ;
273 if (pent_limit != 0x0) ent_limit = *pent_limit ;
275 par_adapt.
add_tbl(ent_limit, 0) ;
281 double precis_poisson = 1.e-16 ;
285 par_poisson1.
add_int(mermax_poisson, 0) ;
296 par_poisson2.
add_int(mermax_poisson, 0) ;
306 Param par_poisson_vect ;
308 par_poisson_vect.
add_int(mermax_poisson, 0) ;
309 par_poisson_vect.
add_double(relax_poisson, 0) ;
310 par_poisson_vect.
add_double(precis_poisson, 1) ;
337 for(
int mer=0 ; mer<mermax ; mer++ ) {
339 cout <<
"-----------------------------------------------" << endl ;
340 cout <<
"step: " << mer << endl ;
341 cout <<
"diff_ent = " << diff_ent << endl ;
350 precis_poisson, relax_potvit) ;
362 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
363 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
369 double alpha_r2 = 0 ;
370 for (
int k=0; k<mg->
get_np(l_b); k++) {
371 for (
int j=0; j<mg->
get_nt(l_b); j++) {
373 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
374 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
378 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
379 ( logn_auto_b - logn_auto_c ) ;
384 if (alpha_r2_jk > alpha_r2) {
385 alpha_r2 = alpha_r2_jk ;
393 alpha_r =
sqrt(alpha_r2) ;
395 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" " 419 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
430 double dent_eq =
ent().dsdr().val_point(
ray_eq(),M_PI/2.,0.) ;
431 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
432 double rap_dent = fabs( dent_eq / dent_pole ) ;
433 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
435 if ( rap_dent < thres_adapt ) {
437 cout <<
"******* FROZEN MAPPING *********" << endl ;
445 if (pent_limit == 0x0) {
447 for (
int l=0; l<
nzet; l++) {
448 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
450 ent_limit.
set(
nzet-1) = ent_b ;
501 double rr_in_1 =
mp.
val_r(
nzet, -1., M_PI/2., 0.) ;
504 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
506 mp.
resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
509 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
511 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
516 double rr_out_nm3_new =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
518 for (
int i=
nzet-1; i<nz-4; i++) {
520 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
522 double rr_mid = rr_out_i
523 + (rr_out_nm3_new - rr_out_i) /
double(nz - 3 - i) ;
525 double rr_2timesi = 2. * rr_out_i ;
527 if (rr_2timesi < rr_mid) {
529 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
531 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
536 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
538 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
557 double ent_s_max = -1 ;
560 for (
int k=0; k<mg->
get_np(l_b); k++) {
561 for (
int j=0; j<mg->
get_nt(l_b); j++) {
562 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
563 if (xx > ent_s_max) {
570 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1" 571 <<
" and nzet : " << endl ;
572 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max <<
573 " and j = " << j_s_max << endl ;
605 source = qpig *
nbar ;
625 "Relative error in the resolution of the equation for logn_auto : " 627 for (
int l=0; l<nz; l++) {
628 cout << tdiff_logn(l) <<
" " ;
631 diff_logn =
max(
abs(tdiff_logn)) ;
663 "Relative error in the resolution of the equation for beta_auto : " 665 for (
int l=0; l<nz; l++) {
666 cout << tdiff_beta(l) <<
" " ;
669 diff_beta =
max(
abs(tdiff_beta)) ;
693 for (
int i=0; i<3; i++) {
695 if (source_shift(i).get_etat() != ETATZERO)
705 for (
int i=0; i<3; i++) {
706 if(source_shift(i).dz_nonzero()) {
707 assert( source_shift(i).get_dzpuis() == 4 ) ;
710 (source_shift.
set(i)).set_dzpuis(4) ;
717 double lambda_shift = double(1) / double(3) ;
719 source_shift.
poisson_vect(lambda_shift, par_poisson_vect,
737 for (
int i=0; i<3; i++) {
739 + lambda_shift * graddivn(i) ;
742 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
743 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
744 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
747 "Relative error in the resolution of the equation for shift_auto : " 749 cout <<
"x component : " ;
750 for (
int l=0; l<nz; l++) {
751 cout << tdiff_shift_x(l) <<
" " ;
754 cout <<
"y component : " ;
755 for (
int l=0; l<nz; l++) {
756 cout << tdiff_shift_y(l) <<
" " ;
759 cout <<
"z component : " ;
760 for (
int l=0; l<nz; l++) {
761 cout << tdiff_shift_z(l) <<
" " ;
765 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
766 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
767 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
784 for (
int ltrans = 0; ltrans <
nzet-1; ltrans++) {
785 cout << endl <<
"Values at boundary between domains no. " << ltrans <<
" and " << ltrans+1 <<
" for theta = pi/2 and phi = 0 :" << endl ;
787 double rt1 =
mp.
val_r(ltrans, 1., M_PI/2, 0.) ;
788 double rt2 =
mp.
val_r(ltrans+1, -1., M_PI/2, 0.) ;
789 cout <<
" Coord. r [km] (left, right, rel. diff) : " 790 << rt1 / km <<
" " << rt2 / km <<
" " << (rt2 - rt1)/rt1 << endl ;
792 int ntm1 = mg->
get_nt(ltrans) - 1;
793 int nrm1 = mg->
get_nr(ltrans) - 1 ;
794 double ent1 =
ent()(ltrans, 0, ntm1, nrm1) ;
795 double ent2 =
ent()(ltrans+1, 0, ntm1, 0) ;
796 cout <<
" Enthalpy (left, right, rel. diff) : " 797 << ent1 <<
" " << ent2 <<
" " << (ent2-ent1)/ent1 << endl ;
799 double press1 =
press()(ltrans, 0, ntm1, nrm1) ;
800 double press2 =
press()(ltrans+1, 0, ntm1, 0) ;
801 cout <<
" Pressure (left, right, rel. diff) : " 802 << press1 <<
" " << press2 <<
" " << (press2-press1)/press1 << endl ;
804 double nb1 =
nbar()(ltrans, 0, ntm1, nrm1) ;
805 double nb2 =
nbar()(ltrans+1, 0, ntm1, 0) ;
806 cout <<
" Baryon density (left, right, rel. diff) : " 807 << nb1 <<
" " << nb2 <<
" " << (nb2-nb1)/nb1 << endl ;
824 diff_ent = diff_ent_tbl(0) ;
825 for (
int l=1; l<
nzet; l++) {
826 diff_ent += diff_ent_tbl(l) ;
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...
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.
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
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).
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.
virtual void homothetie(double lambda)
Sets a new radial scale.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
void equilibrium(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, const Tbl *ent_limit=0x0)
Computes an equilibrium configuration.
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.
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.
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.
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
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...