30char et_rot_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_equilibrium.C,v 1.8 2014/10/13 08:52:57 j_novak Exp $" ;
142#include "graphique.h"
143#include "utilitaires.h"
148 int nzadapt,
const Tbl& ent_limit,
const Itbl& icontrol,
149 const Tbl& control,
double mbar_wanted,
150 double aexp_mass,
Tbl& diff,
Param*) {
159 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
160 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
175 int i_b = mg->
get_nr(l_b) - 1 ;
176 int j_b = mg->
get_nt(l_b) - 1 ;
180 double ent_b = ent_limit(
nzet-1) ;
185 int mer_max = icontrol(0) ;
186 int mer_rot = icontrol(1) ;
187 int mer_change_omega = icontrol(2) ;
188 int mer_fix_omega = icontrol(3) ;
189 int mer_mass = icontrol(4) ;
190 int mermax_poisson = icontrol(5) ;
191 int mer_triax = icontrol(6) ;
192 int delta_mer_kep = icontrol(7) ;
195 if (mer_change_omega < mer_rot) {
196 cout <<
"Etoile_rot::equilibrium: mer_change_omega < mer_rot !" << endl ;
197 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
198 cout <<
" mer_rot = " << mer_rot << endl ;
201 if (mer_fix_omega < mer_change_omega) {
202 cout <<
"Etoile_rot::equilibrium: mer_fix_omega < mer_change_omega !"
204 cout <<
" mer_fix_omega = " << mer_fix_omega << endl ;
205 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
211 bool change_ent = true ;
214 mer_mass =
abs(mer_mass) ;
217 double precis = control(0) ;
218 double omega_ini = control(1) ;
219 double relax = control(2) ;
220 double relax_prev = double(1) - relax ;
221 double relax_poisson = control(3) ;
222 double thres_adapt = control(4) ;
223 double ampli_triax = control(5) ;
224 double precis_adapt = control(6) ;
231 double& diff_ent = diff.
set(0) ;
232 double& diff_nuf = diff.
set(1) ;
233 double& diff_nuq = diff.
set(2) ;
236 double& diff_shift_x = diff.
set(5) ;
237 double& diff_shift_y = diff.
set(6) ;
238 double& vit_triax = diff.
set(7) ;
249 int nz_search =
nzet + 1 ;
252 double reg_map = 1. ;
254 par_adapt.
add_int(nitermax, 0) ;
256 par_adapt.
add_int(nzadapt, 1) ;
259 par_adapt.
add_int(nz_search, 2) ;
261 par_adapt.
add_int(adapt_flag, 3) ;
280 par_adapt.
add_tbl(ent_limit, 0) ;
286 double precis_poisson = 1.e-16 ;
288 Param par_poisson_nuf ;
289 par_poisson_nuf.
add_int(mermax_poisson, 0) ;
290 par_poisson_nuf.
add_double(relax_poisson, 0) ;
291 par_poisson_nuf.
add_double(precis_poisson, 1) ;
295 Param par_poisson_nuq ;
296 par_poisson_nuq.
add_int(mermax_poisson, 0) ;
297 par_poisson_nuq.
add_double(relax_poisson, 0) ;
298 par_poisson_nuq.
add_double(precis_poisson, 1) ;
302 Param par_poisson_tggg ;
303 par_poisson_tggg.
add_int(mermax_poisson, 0) ;
304 par_poisson_tggg.
add_double(relax_poisson, 0) ;
305 par_poisson_tggg.
add_double(precis_poisson, 1) ;
311 Param par_poisson_dzeta ;
318 Param par_poisson_vect ;
320 par_poisson_vect.
add_int(mermax_poisson, 0) ;
321 par_poisson_vect.
add_double(relax_poisson, 0) ;
322 par_poisson_vect.
add_double(precis_poisson, 1) ;
334 double accrois_omega = (omega0 - omega_ini) /
335 double(mer_fix_omega - mer_change_omega) ;
385 ofstream fichconv(
"convergence.d") ;
386 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" << endl ;
388 ofstream fichfreq(
"frequency.d") ;
389 fichfreq <<
"# f [Hz]" << endl ;
391 ofstream fichevol(
"evolution.d") ;
393 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
397 double err_grv2 = 1 ;
398 double max_triax_prev = 0 ;
404 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
406 cout <<
"-----------------------------------------------" << endl ;
407 cout <<
"step: " << mer << endl ;
408 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
410 cout <<
"err_grv2 = " << err_grv2 << endl ;
415 if (mer >= mer_rot) {
417 if (mer < mer_change_omega) {
421 if (mer <= mer_fix_omega) {
422 omega = omega_ini + accrois_omega *
423 (mer - mer_change_omega) ;
445 source_nuf = qpig *
nbar ;
467 (source_tggg.
set()).mult_rsint() ;
488 for (
int i=0; i<3; i++) {
489 source_shift.
set(i) += squad(i) ;
499 source_nuf().poisson(par_poisson_nuf,
nuf.
set()) ;
501 cout <<
"Test of the Poisson equation for nuf :" << endl ;
502 Tbl err = source_nuf().test_poisson(
nuf(), cout,
true) ;
503 diff_nuf = err(0, 0) ;
509 if (mer == mer_triax) {
511 if ( mg->
get_np(0) == 1 ) {
513 "Etoile_rot::equilibrium: np must be stricly greater than 1"
514 << endl <<
" to set a triaxial perturbation !" << endl ;
521 perturb = 1 + ampli_triax * sint*sint *
cos(2*phi) ;
535 double max_triax = 0 ;
537 if ( mg->
get_np(0) > 1 ) {
539 for (
int l=0; l<nz; l++) {
540 for (
int j=0; j<mg->
get_nt(l); j++) {
541 for (
int i=0; i<mg->
get_nr(l); i++) {
544 double xcos2p = (*(va_nuf.
c_cf))(l, 2, j, i) ;
547 double xsin2p = (*(va_nuf.
c_cf))(l, 3, j, i) ;
549 double xx =
sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
551 max_triax = ( xx > max_triax ) ? xx : max_triax ;
558 cout <<
"Triaxial part of nuf : " << max_triax << endl ;
566 source_nuq().poisson(par_poisson_nuq,
nuq.
set()) ;
568 cout <<
"Test of the Poisson equation for nuq :" << endl ;
569 err = source_nuq().test_poisson(
nuq(), cout,
true) ;
570 diff_nuq = err(0, 0) ;
577 if (source_shift.
get_etat() != ETATZERO) {
579 for (
int i=0; i<3; i++) {
580 if(source_shift(i).dz_nonzero()) {
581 assert( source_shift(i).get_dzpuis() == 4 ) ;
584 (source_shift.
set(i)).set_dzpuis(4) ;
592 double lambda_shift = double(1) / double(3) ;
594 if ( mg->
get_np(0) == 1 ) {
598 source_shift.
poisson_vect(lambda_shift, par_poisson_vect,
601 cout <<
"Test of the Poisson equation for shift_x :" << endl ;
602 err = source_shift(0).test_poisson(
shift(0), cout,
true) ;
603 diff_shift_x = err(0, 0) ;
605 cout <<
"Test of the Poisson equation for shift_y :" << endl ;
606 err = source_shift(1).test_poisson(
shift(1), cout,
true) ;
607 diff_shift_y = err(0, 0) ;
626 if (mer > mer_fix_omega + delta_mer_kep) {
628 omega *= fact_omega ;
631 bool omega_trop_grand = false ;
638 bool superlum = true ;
654 ((
uuu.
set()).va).set_base( (tmp.
va).base ) ;
662 for (
int l=0; l<
nzet; l++) {
663 for (
int i=0; i<mg->
get_nr(l); i++) {
665 double u1 =
uuu()(l, 0, j_b, i) ;
668 cout <<
"U > c for l, i : " << l <<
" " << i
669 <<
" U = " << u1 << endl ;
674 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
675 omega /= fact_omega ;
676 cout <<
"New rotation frequency : "
677 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
678 omega_trop_grand = true ;
699 mlngamma = - 0.5 *
uuu*
uuu ;
703 double nuf_b =
nuf()(l_b, k_b, j_b, i_b) ;
704 double nuq_b =
nuq()(l_b, k_b, j_b, i_b) ;
705 double mlngamma_b = mlngamma()(l_b, k_b, j_b, i_b) ;
708 double nuf_c =
nuf()(0,0,0,0) ;
709 double nuq_c =
nuq()(0,0,0,0) ;
710 double mlngamma_c = 0 ;
714 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
715 + nuq_c - nuq_b) / ( nuf_b - nuf_c ) ;
716 alpha_r =
sqrt(alpha_r2) ;
717 cout <<
"alpha_r = " << alpha_r << endl ;
723 double nu_c =
logn()(0,0,0,0) ;
728 ent = (ent_c + nu_c + mlngamma_c) -
logn - mlngamma ;
735 for (
int l=0; l<
nzet; l++) {
736 int imax = mg->
get_nr(l) - 1 ;
737 if (l == l_b) imax-- ;
738 for (
int i=0; i<imax; i++) {
739 if (
ent()(l, 0, j_b, i) < 0. ) {
741 cout <<
"ent < 0 for l, i : " << l <<
" " << i
742 <<
" ent = " <<
ent()(l, 0, j_b, i) << endl ;
748 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" << endl ;
749 omega /= fact_omega ;
750 cout <<
"New rotation frequency : "
751 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
752 omega_trop_grand = true ;
757 if ( omega_trop_grand ) {
759 fact_omega =
sqrt( fact_omega ) ;
760 cout <<
"**** New fact_omega : " << fact_omega << endl ;
775 double dent_eq =
ent().dsdr()(l_b, k_b, j_b, i_b) ;
776 double dent_pole =
ent().dsdr()(l_b, k_b, 0, i_b) ;
777 double rap_dent = fabs( dent_eq / dent_pole ) ;
778 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
780 if ( rap_dent < thres_adapt ) {
782 cout <<
"******* FROZEN MAPPING *********" << endl ;
844 mp.
poisson2d(source_dzf(), source_dzq(), par_poisson_dzeta,
847 err_grv2 = lbda_grv2 - 1;
848 cout <<
"GRV2: " << err_grv2 << endl ;
863 logn = relax *
logn + relax_prev * logn_prev ;
865 dzeta = relax *
dzeta + relax_prev * dzeta_prev ;
877 fichfreq <<
" " <<
omega / (2*M_PI) * f_unit ;
878 fichevol <<
" " << rap_dent ;
880 fichevol <<
" " << ent_c ;
886 if (mer > mer_mass) {
889 if (mbar_wanted > 0.) {
890 xx =
mass_b() / mbar_wanted - 1. ;
891 cout <<
"Discrep. baryon mass <-> wanted bar. mass : " << xx
895 xx =
mass_g() / fabs(mbar_wanted) - 1. ;
896 cout <<
"Discrep. grav. mass <-> wanted grav. mass : " << xx
899 double xprog = ( mer > 2*mer_mass) ? 1. :
900 double(mer-mer_mass)/double(mer_mass) ;
902 double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
903 double fact =
pow(ax, aexp_mass) ;
904 cout <<
" xprog, xx, ax, fact : " << xprog <<
" " <<
905 xx <<
" " << ax <<
" " << fact << endl ;
911 if (mer%4 == 0)
omega *= fact ;
921 diff_ent = diff_ent_tbl(0) ;
922 for (
int l=1; l<
nzet; l++) {
923 diff_ent += diff_ent_tbl(l) ;
927 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
928 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
929 fichconv <<
" " <<
log10( fabs(max_triax) + 1.e-16 ) ;
932 if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
933 vit_triax = (max_triax - max_triax_prev) / max_triax_prev ;
936 fichconv <<
" " << vit_triax ;
945 max_triax_prev = max_triax ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void mult_rsint()
Multiplication by .
Valeur va
The numerical value of the Cmp
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void annule(int l)
Sets the Cmp to zero in a given domain.
Active physical coordinates and mapping derivatives.
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Tenseur uuu
Norm of u_euler.
double omega
Rotation angular velocity ([f_unit] )
Tenseur & logn
Metric potential = logn_auto.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur nuq
Part of the Metric potential = logn generated by the quadratic terms.
virtual double mass_g() const
Gravitational mass.
Tenseur khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Tenseur tggg
Metric potential .
Tenseur nuf
Part of the Metric potential = logn generated by the matter terms.
Cmp ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Tenseur nphi
Metric coefficient .
virtual double mass_b() const
Baryon mass.
Tenseur bbb
Metric factor B.
void update_metric()
Computes metric coefficients from known potentials.
Tenseur & dzeta
Metric potential = beta_auto.
Cmp ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
virtual double grv2() const
Error on the virial identity GRV2.
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
virtual void equilibrium(double ent_c, double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, double mbar_wanted, double aexp_mass, Tbl &diff, Param *=0x0)
Computes an equilibrium configuration.
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Cmp ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Tenseur w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
int nzet
Number of domains of *mp occupied by the star.
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
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.
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
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 shift
Total shift vector.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Tenseur a_car
Total conformal factor .
double ray_pole() const
Coordinate radius at [r_unit].
Basic integer array class.
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 reevaluate(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.
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const =0
Computes the solution of a 2-D Poisson equation.
Coord phi
coordinate centered on the grid
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
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_double_mod(double &x, int position=0)
Adds the address of a new modifiable double 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 set_std_base()
Set the standard spectal basis of decomposition for each component.
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
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 : .
int get_etat() const
Returns the logical state.
Values and coefficients of a (real-value) function.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
Cmp sqrt(const Cmp &)
Square root.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
Cmp abs(const Cmp &)
Absolute value.
Cmp log(const Cmp &)
Neperian logarithm.
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...
Standard units of space, time and mass.