31char et_bin_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_equilibrium.C,v 1.14 2014/10/13 08:52:55 j_novak Exp $" ;
179#include "graphique.h"
180#include "utilitaires.h"
185 int mermax,
int mermax_poisson,
186 double relax_poisson,
int mermax_potvit,
187 double relax_potvit,
double thres_adapt,
188 const Tbl& fact_resize,
Tbl& diff,
const Tbl* pent_limit ) {
207 int i_b = mg->
get_nr(l_b) - 1 ;
217 double& diff_ent = diff.
set(0) ;
218 double& diff_vel_pot = diff.
set(1) ;
219 double& diff_logn = diff.
set(2) ;
220 double& diff_beta = diff.
set(3) ;
221 double& diff_shift_x = diff.
set(4) ;
222 double& diff_shift_y = diff.
set(5) ;
223 double& diff_shift_z = diff.
set(6) ;
235 int nz_search =
nzet ;
238 double precis_secant = 1.e-14 ;
240 double reg_map = 1. ;
242 par_adapt.
add_int(nitermax, 0) ;
247 par_adapt.
add_int(nz_search, 2) ;
249 par_adapt.
add_int(adapt_flag, 3) ;
270 if (pent_limit != 0x0) ent_limit = *pent_limit ;
272 par_adapt.
add_tbl(ent_limit, 0) ;
278 double precis_poisson = 1.e-16 ;
282 par_poisson1.
add_int(mermax_poisson, 0) ;
293 par_poisson2.
add_int(mermax_poisson, 0) ;
303 Param par_poisson_vect ;
305 par_poisson_vect.
add_int(mermax_poisson, 0) ;
306 par_poisson_vect.
add_double(relax_poisson, 0) ;
307 par_poisson_vect.
add_double(precis_poisson, 1) ;
334 for(
int mer=0 ; mer<mermax ; mer++ ) {
336 cout <<
"-----------------------------------------------" << endl ;
337 cout <<
"step: " << mer << endl ;
338 cout <<
"diff_ent = " << diff_ent << endl ;
347 precis_poisson, relax_potvit) ;
359 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
360 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
366 double alpha_r2 = 0 ;
367 for (
int k=0; k<mg->
get_np(l_b); k++) {
368 for (
int j=0; j<mg->
get_nt(l_b); j++) {
370 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
371 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
375 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
376 ( logn_auto_b - logn_auto_c ) ;
381 if (alpha_r2_jk > alpha_r2) {
382 alpha_r2 = alpha_r2_jk ;
390 alpha_r =
sqrt(alpha_r2) ;
392 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
416 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
427 double dent_eq =
ent().dsdr().val_point(
ray_eq(),M_PI/2.,0.) ;
428 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
429 double rap_dent = fabs( dent_eq / dent_pole ) ;
430 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
432 if ( rap_dent < thres_adapt ) {
434 cout <<
"******* FROZEN MAPPING *********" << endl ;
442 if (pent_limit == 0x0) {
444 for (
int l=0; l<
nzet; l++) {
445 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
447 ent_limit.
set(
nzet-1) = ent_b ;
498 double rr_in_1 =
mp.
val_r(
nzet, -1., M_PI/2., 0.) ;
501 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
503 mp.
resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
506 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
508 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
513 double rr_out_nm3_new =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
515 for (
int i=
nzet-1; i<nz-4; i++) {
517 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
519 double rr_mid = rr_out_i
520 + (rr_out_nm3_new - rr_out_i) /
double(nz - 3 - i) ;
522 double rr_2timesi = 2. * rr_out_i ;
524 if (rr_2timesi < rr_mid) {
526 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
528 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
533 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
535 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
554 double ent_s_max = -1 ;
557 for (
int k=0; k<mg->
get_np(l_b); k++) {
558 for (
int j=0; j<mg->
get_nt(l_b); j++) {
559 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
560 if (xx > ent_s_max) {
567 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1"
568 <<
" and nzet : " << endl ;
569 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max <<
570 " and j = " << j_s_max << endl ;
602 source = qpig *
nbar ;
622 "Relative error in the resolution of the equation for logn_auto : "
624 for (
int l=0; l<nz; l++) {
625 cout << tdiff_logn(l) <<
" " ;
628 diff_logn =
max(
abs(tdiff_logn)) ;
660 "Relative error in the resolution of the equation for beta_auto : "
662 for (
int l=0; l<nz; l++) {
663 cout << tdiff_beta(l) <<
" " ;
666 diff_beta =
max(
abs(tdiff_beta)) ;
690 for (
int i=0; i<3; i++) {
692 if (source_shift(i).get_etat() != ETATZERO)
702 for (
int i=0; i<3; i++) {
703 if(source_shift(i).dz_nonzero()) {
704 assert( source_shift(i).get_dzpuis() == 4 ) ;
707 (source_shift.
set(i)).set_dzpuis(4) ;
714 double lambda_shift = double(1) / double(3) ;
716 source_shift.
poisson_vect(lambda_shift, par_poisson_vect,
734 for (
int i=0; i<3; i++) {
736 + lambda_shift * graddivn(i) ;
739 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
740 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
741 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
744 "Relative error in the resolution of the equation for shift_auto : "
746 cout <<
"x component : " ;
747 for (
int l=0; l<nz; l++) {
748 cout << tdiff_shift_x(l) <<
" " ;
751 cout <<
"y component : " ;
752 for (
int l=0; l<nz; l++) {
753 cout << tdiff_shift_y(l) <<
" " ;
756 cout <<
"z component : " ;
757 for (
int l=0; l<nz; l++) {
758 cout << tdiff_shift_z(l) <<
" " ;
762 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
763 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
764 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
781 for (
int ltrans = 0; ltrans <
nzet-1; ltrans++) {
782 cout << endl <<
"Values at boundary between domains no. " << ltrans <<
" and " << ltrans+1 <<
" for theta = pi/2 and phi = 0 :" << endl ;
784 double rt1 =
mp.
val_r(ltrans, 1., M_PI/2, 0.) ;
785 double rt2 =
mp.
val_r(ltrans+1, -1., M_PI/2, 0.) ;
786 cout <<
" Coord. r [km] (left, right, rel. diff) : "
787 << rt1 / km <<
" " << rt2 / km <<
" " << (rt2 - rt1)/rt1 << endl ;
789 int ntm1 = mg->
get_nt(ltrans) - 1;
790 int nrm1 = mg->
get_nr(ltrans) - 1 ;
791 double ent1 =
ent()(ltrans, 0, ntm1, nrm1) ;
792 double ent2 =
ent()(ltrans+1, 0, ntm1, 0) ;
793 cout <<
" Enthalpy (left, right, rel. diff) : "
794 << ent1 <<
" " << ent2 <<
" " << (ent2-ent1)/ent1 << endl ;
796 double press1 =
press()(ltrans, 0, ntm1, nrm1) ;
797 double press2 =
press()(ltrans+1, 0, ntm1, 0) ;
798 cout <<
" Pressure (left, right, rel. diff) : "
799 << press1 <<
" " << press2 <<
" " << (press2-press1)/press1 << endl ;
801 double nb1 =
nbar()(ltrans, 0, ntm1, nrm1) ;
802 double nb2 =
nbar()(ltrans+1, 0, ntm1, 0) ;
803 cout <<
" Baryon density (left, right, rel. diff) : "
804 << nb1 <<
" " << nb2 <<
" " << (nb2-nb1)/nb1 << endl ;
821 diff_ent = diff_ent_tbl(0) ;
822 for (
int l=1; l<
nzet; l++) {
823 diff_ent += diff_ent_tbl(l) ;
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Tenseur akcar_comp
Part of the scalar generated by shift_auto and shift_comp , i.e.
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
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.
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
bool irrotational
true for an irrotational star, false for a corotating one
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
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...
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Tenseur pot_centri
Centrifugal potential.
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Tenseur d_beta_comp
Gradient of beta_comp (Cartesian components with respect to ref_triad )
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
int nzet
Number of domains of *mp occupied by the star.
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Tenseur nbar
Baryon density in the fluid frame.
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
Tenseur press
Fluid pressure.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Tenseur a_car
Total conformal factor .
double ray_pole() const
Coordinate radius at [r_unit].
Radial mapping of rather general form.
virtual void homothetie(double lambda)
Sets a new radial scale.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
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.
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void dec2_dzpuis()
dzpuis -= 2 ;
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void inc2_dzpuis()
dzpuis += 2 ;
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Cmp sqrt(const Cmp &)
Square root.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp abs(const Cmp &)
Absolute value.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Standard units of space, time and mass.