41 #include "star_rot_dirac_diff.h" 43 #include "utilitaires.h" 48 double fact_omega,
int ,
const Tbl& ent_limit,
49 const Itbl& icontrol,
const Tbl& control,
50 double,
double,
Tbl& diff){
59 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
60 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
73 int i_b = mg->
get_nr(l_b) - 1 ;
74 int j_b = mg->
get_nt(l_b) - 1 ;
78 double ent_b = ent_limit(
nzet-1) ;
83 int mer_max = icontrol(0) ;
84 int mer_rot = icontrol(1) ;
85 int mer_change_omega = icontrol(2) ;
86 int mer_fix_omega = icontrol(3) ;
88 int delta_mer_kep = icontrol(5) ;
91 if (mer_change_omega < mer_rot) {
92 cout <<
"Star_rot_Dirac_diff::equilibrium: mer_change_omega < mer_rot !" 94 cout <<
" mer_change_omega = " << mer_change_omega <<
'\n' ;
95 cout <<
" mer_rot = " << mer_rot <<
'\n' ;
98 if (mer_fix_omega < mer_change_omega) {
99 cout <<
"Star_rot_Dirac_diff::equilibrium: mer_fix_omega < mer_change_omega !" 101 cout <<
" mer_fix_omega = " << mer_fix_omega <<
'\n' ;
102 cout <<
" mer_change_omega = " << mer_change_omega <<
'\n' ;
106 double precis = control(0) ;
107 double omega_ini = control(1) ;
108 double relax = control(2) ;
109 double relax_prev = double(1) - relax ;
115 double& diff_ent = diff.
set(0) ;
125 double accrois_omega = (omega_c0 - omega_ini) /
126 double(mer_fix_omega - mer_change_omega) ;
146 ofstream fichconv(
"convergence.d") ;
147 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" <<
'\n' ;
149 ofstream fichfreq(
"frequency.d") ;
150 fichfreq <<
"# f [Hz]" <<
'\n' ;
152 ofstream fichevol(
"evolution.d") ;
154 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c" 158 double err_grv2 = 1 ;
166 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
168 cout <<
"-----------------------------------------------" <<
'\n' ;
169 cout <<
"step: " << mer <<
'\n' ;
170 cout <<
"ent_c = " << display_bold << ent_c << display_normal
172 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
174 cout <<
"err_grv2 = " << err_grv2 <<
'\n' ;
181 if (mer >= mer_rot) {
183 if (mer < mer_change_omega) {
184 omega_c = omega_ini ;
187 if (mer <= mer_fix_omega) {
188 omega_c = omega_ini + accrois_omega *
189 (mer - mer_change_omega) ;
225 if (mer > mer_fix_omega + delta_mer_kep) {
227 omega_c *= fact_omega ;
230 bool omega_trop_grand = false ;
237 bool superlum = true ;
251 double omeg_min = 0 ;
252 double omeg_max = omega_c ;
253 double precis1 = 1.e-14 ;
254 int nitermax1 = 200 ;
284 for (
int l=0; l<
nzet; l++) {
285 for (
int i=0; i<mg->
get_nr(l); i++) {
290 cout <<
"U > c for l, i : " << l <<
" " << i
291 <<
" U = " <<
sqrt( u2 ) <<
'\n' ;
296 cout <<
"**** VELOCITY OF LIGHT REACHED ****" <<
'\n' ;
297 omega /= fact_omega ;
298 cout <<
"New rotation frequency : " 299 <<
omega/(2.*M_PI) * f_unit <<
" Hz" <<
'\n' ;
300 omega_trop_grand = true ;
329 double mlngamma_c = 0 ;
334 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
335 + ln_q_c - ln_q_b + primf_c - primf_b)
336 / ( ln_f_b - ln_f_c ) ;
338 alpha_r =
sqrt(alpha_r2) ;
340 cout <<
"alpha_r = " << alpha_r <<
'\n' ;
349 logn = alpha_r2 * ln_f_new + ln_q_new ;
365 for (
int l=0; l<
nzet; l++) {
366 int imax = mg->
get_nr(l) - 1 ;
367 if (l == l_b) imax-- ;
368 for (
int i=0; i<imax; i++) {
371 cout <<
"ent < 0 for l, i : " << l <<
" " << i
378 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" <<
'\n' ;
379 omega /= fact_omega ;
380 cout <<
"New central rotation frequency : " 381 << omega_c/(2.*M_PI) * f_unit <<
" Hz" <<
'\n' ;
382 omega_trop_grand = true ;
387 if ( omega_trop_grand ) {
389 fact_omega =
sqrt( fact_omega ) ;
390 cout <<
"**** New fact_omega : " << fact_omega <<
'\n' ;
439 logn = relax *
logn + relax_prev * logn_prev ;
441 qqq = relax *
qqq + relax_prev * qqq_prev ;
455 fichfreq <<
" " << omega_c / (2*M_PI) * f_unit ;
456 fichevol <<
" " << ent_c ;
464 diff_ent = diff_ent_tbl(0) ;
465 for (
int l=1; l<
nzet; l++) {
466 diff_ent += diff_ent_tbl(l) ;
470 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
471 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
void solve_logn_f(Scalar &ln_f_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Cmp log(const Cmp &)
Neperian logarithm.
double omega
Rotation angular velocity ([f_unit] )
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Map & mp
Mapping associated with the star.
Cmp sqrt(const Cmp &)
Square root.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
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)
Computes an equilibrium configuration.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
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)
Tensor field of valence 0 (or component of a tensorial field).
Scalar omega_field
Field .
Sym_tensor_trans hh
is defined by .
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Basic integer array class.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Tensor field of valence 1.
void update_metric()
Computes metric quantities from known potentials.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl par_frot
Parameters of the function .
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int nzet
Number of domains of *mp occupied by the star.
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual double grv2() const
Error on the virial identity GRV2.
int get_nzone() const
Returns the number of domains.
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Transverse symmetric tensors of rank 2.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar nn
Lapse function N .
Cmp log10(const Cmp &)
Basis 10 logarithm.
double ray_eq() const
Coordinate radius at , [r_unit].
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
void solve_qqq(Scalar &q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
const Metric_flat & flat
flat metric (spherical components)