28char strot_dirac_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_equilibrium.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $" ;
89#include "star_rot_dirac.h"
91#include "utilitaires.h"
96 double fact_omega,
int ,
const Tbl& ent_limit,
97 const Itbl& icontrol,
const Tbl& control,
98 double mbar_wanted,
double aexp_mass,
Tbl& diff){
107 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
108 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
120 assert( ( type_t == SYM) || (type_t == NONSYM) ) ;
122 int i_b = mg->
get_nr(l_b) - 1 ;
123 int j_b = (type_t == SYM ? mg->
get_nt(l_b) - 1 : mg->
get_nt(l_b)/2 ) ;
127 double ent_b = ent_limit(
nzet-1) ;
132 int mer_max = icontrol(0) ;
133 int mer_rot = icontrol(1) ;
134 int mer_change_omega = icontrol(2) ;
135 int mer_fix_omega = icontrol(3) ;
136 int mer_mass = icontrol(4) ;
137 int delta_mer_kep = icontrol(5) ;
138 int mer_hij = icontrol(6) ;
141 if (mer_change_omega < mer_rot) {
142 cout <<
"Star_rot_Dirac::equilibrium: mer_change_omega < mer_rot !"
144 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
145 cout <<
" mer_rot = " << mer_rot << endl ;
148 if (mer_fix_omega < mer_change_omega) {
149 cout <<
"Star_rot_Dirac::equilibrium: mer_fix_omega < mer_change_omega !"
151 cout <<
" mer_fix_omega = " << mer_fix_omega << endl ;
152 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
158 bool change_ent = true ;
161 mer_mass =
abs(mer_mass) ;
165 double precis = control(0) ;
166 double omega_ini = control(1) ;
167 double relax = control(2) ;
168 double relax_prev = double(1) - relax ;
174 double& diff_ent = diff.
set(0) ;
184 double accrois_omega = (omega0 - omega_ini) /
185 double(mer_fix_omega - mer_change_omega) ;
205 ofstream fichconv(
"convergence.d") ;
206 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" << endl ;
208 ofstream fichfreq(
"frequency.d") ;
209 fichfreq <<
"# f [Hz]" << endl ;
211 ofstream fichevol(
"evolution.d") ;
213 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
217 double err_grv2 = 1 ;
225 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
227 cout <<
"-----------------------------------------------" << endl ;
228 cout <<
"step: " << mer << endl ;
229 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
231 cout <<
"err_grv2 = " << err_grv2 << endl ;
238 if (mer >= mer_rot) {
240 if (mer < mer_change_omega) {
244 if (mer <= mer_fix_omega) {
245 omega = omega_ini + accrois_omega *
246 (mer - mer_change_omega) ;
282 if (mer > mer_fix_omega + delta_mer_kep) {
284 omega *= fact_omega ;
287 bool omega_trop_grand = false ;
294 bool superlum = true ;
323 for (
int l=0; l<
nzet; l++) {
324 for (
int i=0; i<mg->
get_nr(l); i++) {
329 cout <<
"U > c for l, i : " << l <<
" " << i
330 <<
" U = " <<
sqrt( u2 ) << endl ;
335 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
336 omega /= fact_omega ;
337 cout <<
"New rotation frequency : "
338 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
339 omega_trop_grand = true ;
367 double mlngamma_c = 0 ;
371 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
372 + ln_q_c - ln_q_b) / ( ln_f_b - ln_f_c ) ;
374 alpha_r =
sqrt(alpha_r2) ;
376 cout <<
"alpha_r = " << alpha_r << endl ;
385 logn = alpha_r2 * ln_f_new + ln_q_new ;
392 ent = (ent_c + logn_c + mlngamma_c) -
logn - mlngamma ;
401 for (
int l=0; l<
nzet; l++) {
402 int imax = mg->
get_nr(l) - 1 ;
403 if (l == l_b) imax-- ;
404 for (
int i=0; i<imax; i++) {
407 cout <<
"ent < 0 for l, i : " << l <<
" " << i
414 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" << endl ;
415 omega /= fact_omega ;
416 cout <<
"New rotation frequency : "
417 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
418 omega_trop_grand = true ;
423 if ( omega_trop_grand ) {
425 fact_omega =
sqrt( fact_omega ) ;
426 cout <<
"**** New fact_omega : " << fact_omega << endl ;
478 logn = relax *
logn + relax_prev * logn_prev ;
480 qqq = relax *
qqq + relax_prev * qqq_prev ;
494 fichfreq <<
" " <<
omega / (2*M_PI) * f_unit ;
495 fichevol <<
" " << ent_c ;
502 if (mer > mer_mass) {
505 if (mbar_wanted > 0.) {
506 xx =
mass_b() / mbar_wanted - 1. ;
507 cout <<
"Discrep. baryon mass <-> wanted bar. mass : " << xx
511 xx =
mass_g() / fabs(mbar_wanted) - 1. ;
512 cout <<
"Discrep. grav. mass <-> wanted grav. mass : " << xx
515 double xprog = ( mer > 2*mer_mass) ? 1. :
516 double(mer-mer_mass)/double(mer_mass) ;
518 double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
519 double fact =
pow(ax, aexp_mass) ;
520 cout <<
" xprog, xx, ax, fact : " << xprog <<
" " <<
521 xx <<
" " << ax <<
" " << fact << endl ;
527 if (mer%4 == 0)
omega *= fact ;
538 diff_ent = diff_ent_tbl(0) ;
539 for (
int l=1; l<
nzet; l++) {
540 diff_ent += diff_ent_tbl(l) ;
544 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
545 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
Basic integer array class.
virtual void homothetie(double lambda)=0
Sets a new radial scale.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
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.
Tensor field of valence 0 (or component of a tensorial field).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
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.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
virtual double mass_g() const
Gravitational mass.
Sym_tensor_trans hh
is defined by .
double omega
Rotation angular velocity ([f_unit] )
const Metric_flat & flat
flat metric (spherical components)
void solve_logn_f(Scalar &ln_f_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.
void update_metric()
Computes metric quantities from known potentials.
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
virtual double mass_b() const
Baryonic mass.
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.
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
virtual double grv2() const
Error on the virial identity GRV2.
Scalar logn
Logarithm of the lapse N .
Scalar nn
Lapse function N .
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
int nzet
Number of domains of *mp occupied by the star.
Transverse symmetric tensors of rank 2.
void annule_hard()
Sets the Tbl to zero in a hard way.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 1.
Scalar & set(int)
Read/write access to a component.
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 abs(const Cmp &)
Absolute value.
Cmp log(const Cmp &)
Neperian logarithm.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Standard units of space, time and mass.