30char et_rot_diff_equil_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff_equil.C,v 1.9 2014/10/13 08:52:57 j_novak Exp $" ;
80#include "et_rot_diff.h"
84#include "utilitaires.h"
89 int nzadapt,
const Tbl& ent_limit,
91 const Tbl& control,
double mbar_wanted,
92 double aexp_mass,
Tbl& diff,
Param*) {
100 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
101 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
116 int i_b = mg->
get_nr(l_b) - 1 ;
117 int j_b = mg->
get_nt(l_b) - 1 ;
121 double ent_b = ent_limit(
nzet-1) ;
126 int mer_max = icontrol(0) ;
127 int mer_rot = icontrol(1) ;
128 int mer_change_omega = icontrol(2) ;
129 int mer_fix_omega = icontrol(3) ;
130 int mer_mass = icontrol(4) ;
131 int mermax_poisson = icontrol(5) ;
132 int mer_triax = icontrol(6) ;
133 int delta_mer_kep = icontrol(7) ;
136 if (mer_change_omega < mer_rot) {
137 cout <<
"Et_rot_diff::equilibrium: mer_change_omega < mer_rot !" << endl ;
138 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
139 cout <<
" mer_rot = " << mer_rot << endl ;
142 if (mer_fix_omega < mer_change_omega) {
143 cout <<
"Et_rot_diff::equilibrium: mer_fix_omega < mer_change_omega !"
145 cout <<
" mer_fix_omega = " << mer_fix_omega << endl ;
146 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
152 bool change_ent = true ;
155 mer_mass =
abs(mer_mass) ;
158 double precis = control(0) ;
159 double omega_ini = control(1) ;
160 double relax = control(2) ;
161 double relax_prev = double(1) - relax ;
162 double relax_poisson = control(3) ;
163 double thres_adapt = control(4) ;
164 double ampli_triax = control(5) ;
165 double precis_adapt = control(6) ;
172 double& diff_ent = diff.
set(0) ;
173 double& diff_nuf = diff.
set(1) ;
174 double& diff_nuq = diff.
set(2) ;
177 double& diff_shift_x = diff.
set(5) ;
178 double& diff_shift_y = diff.
set(6) ;
179 double& vit_triax = diff.
set(7) ;
190 int nz_search =
nzet + 1 ;
193 double reg_map = 1. ;
195 par_adapt.
add_int(nitermax, 0) ;
197 par_adapt.
add_int(nzadapt, 1) ;
200 par_adapt.
add_int(nz_search, 2) ;
202 par_adapt.
add_int(adapt_flag, 3) ;
221 par_adapt.
add_tbl(ent_limit, 0) ;
227 double precis_poisson = 1.e-16 ;
229 Param par_poisson_nuf ;
230 par_poisson_nuf.
add_int(mermax_poisson, 0) ;
231 par_poisson_nuf.
add_double(relax_poisson, 0) ;
232 par_poisson_nuf.
add_double(precis_poisson, 1) ;
236 Param par_poisson_nuq ;
237 par_poisson_nuq.
add_int(mermax_poisson, 0) ;
238 par_poisson_nuq.
add_double(relax_poisson, 0) ;
239 par_poisson_nuq.
add_double(precis_poisson, 1) ;
243 Param par_poisson_tggg ;
244 par_poisson_tggg.
add_int(mermax_poisson, 0) ;
245 par_poisson_tggg.
add_double(relax_poisson, 0) ;
246 par_poisson_tggg.
add_double(precis_poisson, 1) ;
252 Param par_poisson_dzeta ;
259 Param par_poisson_vect ;
261 par_poisson_vect.
add_int(mermax_poisson, 0) ;
262 par_poisson_vect.
add_double(relax_poisson, 0) ;
263 par_poisson_vect.
add_double(precis_poisson, 1) ;
275 double accrois_omega = (omega_c0 - omega_ini) /
276 double(mer_fix_omega - mer_change_omega) ;
326 ofstream fichconv(
"convergence.d") ;
327 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" << endl ;
329 ofstream fichfreq(
"frequency.d") ;
330 fichfreq <<
"# f [Hz]" << endl ;
332 ofstream fichevol(
"evolution.d") ;
334 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
338 double err_grv2 = 1 ;
339 double max_triax_prev = 0 ;
345 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
347 cout <<
"-----------------------------------------------" << endl ;
348 cout <<
"step: " << mer << endl ;
349 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
351 cout <<
"err_grv2 = " << err_grv2 << endl ;
356 if (mer >= mer_rot) {
358 if (mer < mer_change_omega) {
359 omega_c = omega_ini ;
362 if (mer <= mer_fix_omega) {
363 omega_c = omega_ini + accrois_omega *
364 (mer - mer_change_omega) ;
386 source_nuf = qpig *
nbar ;
408 (source_tggg.
set()).mult_rsint() ;
429 for (
int i=0; i<3; i++) {
430 source_shift.
set(i) += squad(i) ;
440 source_nuf().poisson(par_poisson_nuf,
nuf.
set()) ;
442 cout <<
"Test of the Poisson equation for nuf :" << endl ;
443 Tbl err = source_nuf().test_poisson(
nuf(), cout,
true) ;
444 diff_nuf = err(0, 0) ;
450 if (mer == mer_triax) {
452 if ( mg->
get_np(0) == 1 ) {
454 "Et_rot_diff::equilibrium: np must be stricly greater than 1"
455 << endl <<
" to set a triaxial perturbation !" << endl ;
462 perturb = 1 + ampli_triax * sint*sint *
cos(2*phi) ;
476 double max_triax = 0 ;
478 if ( mg->
get_np(0) > 1 ) {
480 for (
int l=0; l<nz; l++) {
481 for (
int j=0; j<mg->
get_nt(l); j++) {
482 for (
int i=0; i<mg->
get_nr(l); i++) {
485 double xcos2p = (*(va_nuf.
c_cf))(l, 2, j, i) ;
488 double xsin2p = (*(va_nuf.
c_cf))(l, 3, j, i) ;
490 double xx =
sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
492 max_triax = ( xx > max_triax ) ? xx : max_triax ;
499 cout <<
"Triaxial part of nuf : " << max_triax << endl ;
507 source_nuq().poisson(par_poisson_nuq,
nuq.
set()) ;
509 cout <<
"Test of the Poisson equation for nuq :" << endl ;
510 err = source_nuq().test_poisson(
nuq(), cout,
true) ;
511 diff_nuq = err(0, 0) ;
518 if (source_shift.
get_etat() != ETATZERO) {
520 for (
int i=0; i<3; i++) {
521 if(source_shift(i).dz_nonzero()) {
522 assert( source_shift(i).get_dzpuis() == 4 ) ;
525 (source_shift.
set(i)).set_dzpuis(4) ;
533 double lambda_shift = double(1) / double(3) ;
535 if ( mg->
get_np(0) == 1 ) {
539 source_shift.
poisson_vect(lambda_shift, par_poisson_vect,
542 cout <<
"Test of the Poisson equation for shift_x :" << endl ;
543 err = source_shift(0).test_poisson(
shift(0), cout,
true) ;
544 diff_shift_x = err(0, 0) ;
546 cout <<
"Test of the Poisson equation for shift_y :" << endl ;
547 err = source_shift(1).test_poisson(
shift(1), cout,
true) ;
548 diff_shift_y = err(0, 0) ;
562 if (mer > mer_fix_omega + delta_mer_kep) {
564 omega_c *= fact_omega ;
568 bool omega_trop_grand = false ;
575 bool superlum = true ;
589 double omeg_min = 0 ;
590 double omeg_max = omega_c ;
591 double precis1 = 1.e-14 ;
592 int nitermax1 = 100 ;
609 ((
uuu.
set()).va).set_base( (tmp.
va).base ) ;
617 for (
int l=0; l<
nzet; l++) {
618 for (
int i=0; i<mg->
get_nr(l); i++) {
620 double u1 =
uuu()(l, 0, j_b, i) ;
623 cout <<
"U > c for l, i : " << l <<
" " << i
624 <<
" U = " << u1 << endl ;
629 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
630 omega_c /= fact_omega ;
631 cout <<
"New central rotation frequency : "
632 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
633 omega_trop_grand = true ;
654 mlngamma = - 0.5 *
uuu*
uuu ;
658 double nuf_b =
nuf()(l_b, k_b, j_b, i_b) ;
659 double nuq_b =
nuq()(l_b, k_b, j_b, i_b) ;
660 double mlngamma_b = mlngamma()(l_b, k_b, j_b, i_b) ;
661 double primf_b =
prim_field()(l_b, k_b, j_b, i_b) ;
665 double nuf_c =
nuf()(0,0,0,0) ;
666 double nuq_c =
nuq()(0,0,0,0) ;
667 double mlngamma_c = 0 ;
672 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
673 + nuq_c - nuq_b + primf_c - primf_b)
674 / ( nuf_b - nuf_c ) ;
675 alpha_r =
sqrt(alpha_r2) ;
676 cout <<
"alpha_r = " << alpha_r << endl ;
682 double nu_c =
logn()(0,0,0,0) ;
694 for (
int l=0; l<
nzet; l++) {
695 int imax = mg->
get_nr(l) - 1 ;
696 if (l == l_b) imax-- ;
697 for (
int i=0; i<imax; i++) {
698 if (
ent()(l, 0, j_b, i) < 0. ) {
700 cout <<
"ent < 0 for l, i : " << l <<
" " << i
701 <<
" ent = " <<
ent()(l, 0, j_b, i) << endl ;
707 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" << endl ;
708 omega_c /= fact_omega ;
709 cout <<
"New central rotation frequency : "
710 << omega_c/(2.*M_PI) * f_unit <<
" Hz" << endl ;
711 omega_trop_grand = true ;
716 if ( omega_trop_grand ) {
718 fact_omega =
sqrt( fact_omega ) ;
719 cout <<
"**** New fact_omega : " << fact_omega << endl ;
730 double dent_eq =
ent().dsdr()(l_b, k_b, j_b, i_b) ;
731 double dent_pole =
ent().dsdr()(l_b, k_b, 0, i_b) ;
732 double rap_dent = fabs( dent_eq / dent_pole ) ;
733 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
735 if ( rap_dent < thres_adapt ) {
737 cout <<
"******* FROZEN MAPPING *********" << endl ;
789 mp.
poisson2d(source_dzf(), source_dzq(), par_poisson_dzeta,
792 err_grv2 = lbda_grv2 - 1;
793 cout <<
"GRV2: " << err_grv2 << endl ;
808 logn = relax *
logn + relax_prev * logn_prev ;
810 dzeta = relax *
dzeta + relax_prev * dzeta_prev ;
822 fichfreq <<
" " << omega_c / (2*M_PI) * f_unit ;
823 fichevol <<
" " << rap_dent ;
825 fichevol <<
" " << ent_c ;
831 if (mer > mer_mass) {
834 if (mbar_wanted > 0.) {
835 xx =
mass_b() / mbar_wanted - 1. ;
836 cout <<
"Discrep. baryon mass <-> wanted bar. mass : " << xx
840 xx =
mass_g() / fabs(mbar_wanted) - 1. ;
841 cout <<
"Discrep. grav. mass <-> wanted grav. mass : " << xx
844 double xprog = ( mer > 2*mer_mass) ? 1. :
845 double(mer-mer_mass)/double(mer_mass) ;
847 double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
848 double fact =
pow(ax, aexp_mass) ;
849 cout <<
" xprog, xx, ax, fact : " << xprog <<
" " <<
850 xx <<
" " << ax <<
" " << fact << endl ;
856 if (mer%4 == 0) omega_c *= fact ;
866 diff_ent = diff_ent_tbl(0) ;
867 for (
int l=1; l<
nzet; l++) {
868 diff_ent += diff_ent_tbl(l) ;
872 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
873 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
874 fichconv <<
" " <<
log10( fabs(max_triax) + 1.e-16 ) ;
877 if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
878 vit_triax = (max_triax - max_triax_prev) / max_triax_prev ;
881 fichconv <<
" " << vit_triax ;
890 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.
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
Tenseur prim_field
Field .
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur omega_field
Field .
Tbl par_frot
Parameters of the function .
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 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.
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...
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.