30char gravastar_equil_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Gravastar/gravastar_equil.C,v 1.3 2014/10/13 08:52:59 j_novak Exp $" ;
43#include "utilitaires.h"
48 int nzadapt,
const Tbl& ent_limit,
49 const Itbl& icontrol,
const Tbl& control,
59 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
60 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
75 int i_b = mg->
get_nr(l_b) - 1 ;
76 int j_b = mg->
get_nt(l_b) - 1 ;
80 double ent_b = ent_limit(
nzet-1) ;
85 int mer_max = icontrol(0) ;
86 int mer_rot = icontrol(1) ;
87 int mer_change_omega = icontrol(2) ;
88 int mer_fix_omega = icontrol(3) ;
89 int mermax_poisson = icontrol(4) ;
90 int delta_mer_kep = icontrol(5) ;
93 if (mer_change_omega < mer_rot) {
94 cout <<
"Gravastar::equilibrium: mer_change_omega < mer_rot !" << endl ;
95 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
96 cout <<
" mer_rot = " << mer_rot << endl ;
99 if (mer_fix_omega < mer_change_omega) {
100 cout <<
"Gravastar::equilibrium: mer_fix_omega < mer_change_omega !"
102 cout <<
" mer_fix_omega = " << mer_fix_omega << endl ;
103 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
117 double precis = control(0) ;
118 double omega_ini = control(1) ;
119 double relax = control(2) ;
120 double relax_prev = double(1) - relax ;
121 double relax_poisson = control(3) ;
122 double thres_adapt = control(4) ;
123 double precis_adapt = control(5) ;
130 double& diff_ent = diff.
set(0) ;
131 double& diff_nuf = diff.
set(1) ;
132 double& diff_nuq = diff.
set(2) ;
135 double& diff_shift_x = diff.
set(5) ;
136 double& diff_shift_y = diff.
set(6) ;
148 int nz_search =
nzet + 1 ;
151 double reg_map = 1. ;
153 par_adapt.
add_int(nitermax, 0) ;
155 par_adapt.
add_int(nzadapt, 1) ;
158 par_adapt.
add_int(nz_search, 2) ;
160 par_adapt.
add_int(adapt_flag, 3) ;
179 par_adapt.
add_tbl(ent_limit, 0) ;
185 double precis_poisson = 1.e-16 ;
193 for (
int i=1; i<=3; i++) {
199 for (
int i=1; i<=3; i++) {
205 for (
int i=1; i<=3; i++) {
213 Param par_poisson_nuf ;
214 par_poisson_nuf.
add_int(mermax_poisson, 0) ;
215 par_poisson_nuf.
add_double(relax_poisson, 0) ;
216 par_poisson_nuf.
add_double(precis_poisson, 1) ;
220 Param par_poisson_nuq ;
221 par_poisson_nuq.
add_int(mermax_poisson, 0) ;
222 par_poisson_nuq.
add_double(relax_poisson, 0) ;
223 par_poisson_nuq.
add_double(precis_poisson, 1) ;
227 Param par_poisson_tggg ;
228 par_poisson_tggg.
add_int(mermax_poisson, 0) ;
229 par_poisson_tggg.
add_double(relax_poisson, 0) ;
230 par_poisson_tggg.
add_double(precis_poisson, 1) ;
236 Param par_poisson_dzeta ;
243 Param par_poisson_vect ;
245 par_poisson_vect.
add_int(mermax_poisson, 0) ;
246 par_poisson_vect.
add_double(relax_poisson, 0) ;
247 par_poisson_vect.
add_double(precis_poisson, 1) ;
259 double accrois_omega = (omega0 - omega_ini) /
260 double(mer_fix_omega - mer_change_omega) ;
299 ofstream fichconv(
"convergence.d") ;
300 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" << endl ;
302 ofstream fichfreq(
"frequency.d") ;
303 fichfreq <<
"# f [Hz]" << endl ;
305 ofstream fichevol(
"evolution.d") ;
307 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_cr"
311 double err_grv2 = 1 ;
318 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
320 cout <<
"-----------------------------------------------" << endl ;
321 cout <<
"step: " << mer << endl ;
322 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
324 cout <<
"err_grv2 = " << err_grv2 << endl ;
329 if (mer >= mer_rot) {
331 if (mer < mer_change_omega) {
335 if (mer <= mer_fix_omega) {
336 omega = omega_ini + accrois_omega *
337 (mer - mer_change_omega) ;
360 source_nuq =
ak_car - d_logn(1)*(d_logn(1)+d_bet(1))
361 - d_logn(2)*(d_logn(2)+d_bet(2))
362 - d_logn(3)*(d_logn(3)+d_bet(3)) ;
366 source_nuf = qpig *
nbar ;
393 - d_logn(1)*d_logn(1) - d_logn(2)*d_logn(2) - d_logn(3)*d_logn(3) ;
427 cout <<
"$$$???$$$ Test of the Poisson equation for nuf :" << endl ;
429 diff_nuf = err(0, 0) ;
494 cout <<
"Test of the Poisson equation for nuq :" << endl ;
496 diff_nuq = err(0, 0) ;
503 for (
int i=1; i<=3; i++) {
504 if(source_shift(i).get_etat() != ETATZERO) {
505 if(source_shift(i).dz_nonzero()) {
506 assert( source_shift(i).get_dzpuis() == 4 ) ;
509 (source_shift.
set(i)).set_dzpuis(4) ;
514 double lambda_shift = double(1) / double(3) ;
516 if ( mg->
get_np(0) == 1 ) {
522 for (
int i=1; i<=3; i++) {
523 csource_shift.
set(i-1) = source_shift(i) ;
527 csource_shift.
poisson_vect(lambda_shift, par_poisson_vect,
528 cshift, cw_shift, ckhi_shift) ;
530 for (
int i=1; i<=3; i++) {
537 cout <<
"Test of the Poisson equation for shift_x :" << endl ;
538 err = source_shift(1).test_poisson(-
beta(1), cout,
true) ;
539 diff_shift_x = err(0, 0) ;
541 cout <<
"Test of the Poisson equation for shift_y :" << endl ;
542 err = source_shift(2).test_poisson(-
beta(2), cout,
true) ;
543 diff_shift_y = err(0, 0) ;
557 if (mer > mer_fix_omega + delta_mer_kep) {
559 omega *= fact_omega ;
562 bool omega_trop_grand = false ;
569 bool superlum = true ;
592 for (
int l=0; l<
nzet; l++) {
593 for (
int i=0; i<mg->
get_nr(l); i++) {
598 cout <<
"U > c for l, i : " << l <<
" " << i
599 <<
" U = " << u1 << endl ;
604 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
605 omega /= fact_omega ;
606 cout <<
"New rotation frequency : "
607 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
608 omega_trop_grand = true ;
629 mlngamma = - 0.5 *
uuu*
uuu ;
654 int j_cr = mg->
get_nt(l_cr) - 1 ;
658 double mlngamma_cr = mlngamma.
val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
673 double alpha_r2 = ( ent_cr - ent_b + mlngamma_cr - mlngamma_b
674 + nuq_cr - nuq_b) / ( nuf_b - nuf_cr ) ;
675 alpha_r =
sqrt(alpha_r2) ;
676 cout <<
"ent_cr,ent_b,mlngamma_cr,mlngamma_b,nuq_cr,nuq_b " << ent_cr <<
" " << ent_b <<
" " << mlngamma_cr <<
" " << mlngamma_b <<
" " << nuq_cr <<
" " << nuq_b << endl;
677 cout <<
"nuf_b, nuf_cr= " << nuf_b <<
" " << nuf_cr << endl;
678 cout <<
"num= " << ent_cr - ent_b + mlngamma_cr - mlngamma_b
679 + nuq_cr - nuq_b << endl;
680 cout <<
"deno= " << nuf_b - nuf_cr << endl;
681 cout <<
"alpha_r = " << alpha_r << endl ;
682 if (alpha_r != alpha_r){
683 cout <<
"alpha_r nan!" << endl;
695 cout <<
"ent_cr, ent_crbis, nu_cr= " << ent_cr <<
" " <<
ent.
val_grid_point(l_cr,k_cr,j_cr,i_cr) <<
" " << nu_cr << endl;
698 ent = (ent_cr + nu_cr + mlngamma_cr) -
logn - mlngamma ;
717 for (
int l=0; l<
nzet; l++) {
718 int imax = mg->
get_nr(l) - 1 ;
719 if (l == l_b) imax-- ;
720 for (
int i=0; i<imax; i++) {
723 cout <<
"ent < 0 for l, i : " << l <<
" " << i
730 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" << endl ;
731 omega /= fact_omega ;
732 cout <<
"New rotation frequency : "
733 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
734 omega_trop_grand = true ;
739 if ( omega_trop_grand ) {
741 fact_omega =
sqrt( fact_omega ) ;
742 cout <<
"**** New fact_omega : " << fact_omega << endl ;
754 double rap_dent = fabs( dent_eq / dent_pole ) ;
755 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
757 if ( rap_dent < thres_adapt ) {
759 cout <<
"******* FROZEN MAPPING *********" << endl ;
787 des_profile(
ent,0.,10.,0.,0.,
"enthalpy apres mapping (used for eos)");
818 Cmp csource_tggg(source_tggg) ;
829 Cmp csource_dzf(source_dzf) ;
830 Cmp csource_dzq(source_dzq) ;
832 mp.
poisson2d(csource_dzf, csource_dzq, par_poisson_dzeta,
836 err_grv2 = lbda_grv2 - 1;
837 cout <<
"GRV2: " << err_grv2 << endl ;
852 logn = relax *
logn + relax_prev * logn_prev ;
854 dzeta = relax *
dzeta + relax_prev * dzeta_prev ;
866 fichfreq <<
" " <<
omega / (2*M_PI) * f_unit ;
867 fichevol <<
" " << rap_dent ;
911 diff_ent = diff_ent_tbl(0) ;
912 for (
int l=1; l<
nzet; l++) {
913 diff_ent += diff_ent_tbl(l) ;
917 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
918 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
954 for (
int i=1; i<=3; i++) {
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
void equilibrium(double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, Tbl &diff, Param *=0x0)
Computes an equilibrium configuration.
void equation_of_state()
Allows to computes the proper baryon and energy density, as well as pressure from the enthalpy,...
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.
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
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.
Tensor field of valence 0 (or component of a tensorial field).
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
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.
Tbl test_poisson(const Scalar &uu, ostream &ostr, bool detail=false) const
Checks if a Poisson equation with *this as a source has been correctly solved.
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
Scalar tggg
Metric potential .
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Scalar bbb
Metric factor B.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double omega
Rotation angular velocity ([f_unit] )
Scalar uuu
Norm of u_euler.
Scalar nphi
Metric coefficient .
void update_metric()
Computes metric coefficients from known potentials.
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Scalar dzeta
Metric potential .
Scalar a_car
Square of the metric factor A.
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
virtual double grv2() const
Error on the virial identity GRV2.
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Scalar logn
Logarithm of the lapse N .
Scalar nn
Lapse function N .
Scalar nbar
Baryon density in the fluid frame.
Scalar ener_euler
Total energy density in the Eulerian frame.
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
double ray_eq() const
Coordinate radius at , [r_unit].
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Scalar press
Fluid pressure.
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.
double ray_pole() const
Coordinate radius at [r_unit].
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 poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Tensor field of valence 1.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
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 log(const Cmp &)
Neperian logarithm.
void des_profile(const float *uutab, int nx, float xmin, float xmax, const char *nomx, const char *nomy, const char *title, const char *device=0x0, int nbound=0, float *xbound=0x0)
Basic routine for drawing a single profile with uniform x sampling.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Standard units of space, time and mass.