117 #include "et_rot_bifluid.h" 119 #include "utilitaires.h" 129 Etoile_rot(mpi, nzet_i, relat, *eos_i.trans2Eos()),
138 j_euler(mpi, 1, CON, mp.get_bvect_cart()),
139 j_euler1 (mpi, 1, CON, mp.get_bvect_cart()),
140 j_euler2(mpi, 1, CON, mp.get_bvect_cart()),
182 alpha_eos(et.alpha_eos),
183 sphph_euler(et.sphph_euler),
185 j_euler1(et.j_euler1),
186 j_euler2(et.j_euler2),
187 enerps_euler(et.enerps_euler),
189 gam_euler2(et.gam_euler2),
190 delta_car(et.delta_car)
211 j_euler(mpi, 1, CON, mp.get_bvect_cart()),
212 j_euler1(mpi, 1, CON, mp.get_bvect_cart()),
213 j_euler2(mpi, 1, CON, mp.get_bvect_cart()),
357 assert( &(et.
eos) == &
eos ) ;
407 ost <<
"Bifluid rotating star" << endl ;
408 ost <<
"-------------" << endl ;
409 ost << setprecision(16);
410 double freq =
omega / (2.*M_PI) ;
411 ost <<
"Omega1 : " <<
omega * f_unit
412 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
413 ost <<
"Rotation period 1: " << 1000. / (freq * f_unit) <<
" ms" 416 double freq2 =
omega2 / (2.*M_PI) ;
417 ost <<
"Omega2 : " <<
omega2 * f_unit
418 <<
" rad/s f : " << freq2 * f_unit <<
" Hz" << endl ;
419 ost <<
"Rotation period 2: " << 1000. / (freq2 * f_unit) <<
" ms" 423 ost <<
"Total angular momentum J : " 424 <<
angu_mom()/( qpig / (4* M_PI) * msol*msol) <<
" G M_sol^2 / c" 426 ost <<
"c J / (G M^2) : " 430 ost <<
"Total moment of inertia I = J/Omega2 : " 431 << mom_iner <<
" Lorene units" 433 double mom_iner_38si = mom_iner * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
434 ost <<
"Total moment of inertia I = J/Omega2 : " 435 << mom_iner_38si <<
" 10^38 kg m^2" 439 ost <<
"Angular momentum J of fluid 1 : " 440 <<
angu_mom_1()/( qpig / (4* M_PI) * msol*msol) <<
" G M_sol^2 / c" 442 ost <<
"Angular momentum J of fluid 2 : " 443 <<
angu_mom_2()/( qpig / (4* M_PI) * msol*msol) <<
" G M_sol^2 / c" 452 double mom_iner_1 = 0.;
453 double mom_iner_2 = 0.;
455 if (
omega != __infinity) {
457 ost <<
"Partial moments of inertia (defined in corotation only) :" << endl;
460 ost <<
"Moment of inertia of fluid 1 : " << mom_iner_1 <<
" Lorene units " << endl ;
461 double mom_iner_1_38si = mom_iner_1 * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
462 ost <<
"Moment of inertia of fluid 1 : " << mom_iner_1_38si <<
" 10^38 kg m^2" 466 ost <<
"Moment of inertia of fluid 2 : " << mom_iner_2 <<
" Lorene units " << endl ;
467 double mom_iner_2_38si = mom_iner_2 * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
468 ost <<
"Moment of inertia of fluid 2 : " << mom_iner_2_38si <<
" 10^38 kg m^2" 473 ost <<
"***** Fluid coupling quantities *****" << endl;
475 double mominert_1_tilde_si =
coupling_mominert_1() * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
476 double mominert_2_tilde_si = coupling_mominert_2() * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
478 ost <<
"tilde{I}_n : " <<
coupling_mominert_1() <<
" SI: " << mominert_1_tilde_si <<
" 10^38 kg m^2" 480 ost <<
"tilde{I}_p : " << coupling_mominert_2() <<
" SI: " << mominert_2_tilde_si <<
" 10^38 kg m^2" 483 ost <<
"tilde{varepsilon}_p : " << coupling_entr() / coupling_mominert_2() << endl;
485 ost <<
"tilde{omega}_p : " << coupling_LT_2() / coupling_mominert_2() << endl;
487 <<
" Jp = " << coupling_mominert_2() *
omega2 - coupling_LT_2() + coupling_entr() * (
omega -
omega2)
492 double nphi_c =
nphi()(0, 0, 0, 0) ;
493 if ( (
omega==0) && (nphi_c==0) ) {
494 ost <<
"Central N^phi : " << nphi_c << endl ;
497 ost <<
"Central N^phi/Omega : " << nphi_c /
omega << endl ;
500 ost <<
"Central N^phi/Omega2 : " << nphi_c /
omega2 << endl ;
502 ost <<
"Error on the virial identity GRV2 : " << endl ;
503 ost <<
"GRV2 = " <<
grv2() << endl ;
504 ost <<
"Error on the virial identity GRV3 : " << endl ;
505 double xgrv3 =
grv3(&ost) ;
506 ost <<
"GRV3 = " << xgrv3 << endl ;
508 ost <<
"Circumferential equatorial radius R_circ : " 509 <<
r_circ()/km <<
" km" << endl ;
510 ost <<
"Mean radius R_mean : " 512 ost <<
"Coordinate equatorial radius r_eq : " <<
ray_eq()/km <<
" km" 514 ost <<
"Flattening r_pole/r_eq : " <<
aplat() << endl ;
515 ost <<
"Circumferential equatorial radius R_circ2 : " 516 <<
r_circ2()/km <<
" km" << endl ;
517 ost <<
"Mean radius R_mean2 : " 519 ost <<
"Coordinate equatorial radius r_eq2 : " <<
ray_eq2()/km <<
" km" 521 ost <<
"Flattening r_pole2/r_eq2 : " <<
aplat2() << endl ;
523 int lsurf =
nzet - 1;
526 ost <<
"Equatorial value of the velocity U: " 527 <<
uuu()(
nzet-1, 0, nt-1, nr-1) <<
" c" << endl ;
528 ost <<
"Equatorial value of the velocity U2: " 529 <<
uuu2()(
nzet-1, 0, nt-1, nr-1) <<
" c" << endl ;
530 ost <<
"Redshift at the equator (forward) : " <<
z_eqf() << endl ;
531 ost <<
"Redshift at the equator (backward): " <<
z_eqb() << endl ;
532 ost <<
"Redshift at the pole : " <<
z_pole() << endl ;
535 ost <<
"Central value of log(N) : " 536 <<
logn()(0, 0, 0, 0) << endl ;
538 ost <<
"Central value of dzeta=log(AN) : " 539 <<
dzeta()(0, 0, 0, 0) << endl ;
541 ost <<
"Central value of B^2 : " <<
b_car()(0,0,0,0) << endl ;
545 "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : " 547 for (
int l=0; l<diff_a_b.get_taille(); l++) {
548 ost << diff_a_b(l) <<
" " ;
551 ost <<
"Quadrupole moment : " <<
mom_quad() << endl ;
552 double mom_quad_38si =
mom_quad() * rho_unit * (
pow(r_unit,
double(5.)) / double(1.e38) ) ;
553 ost <<
"Quadrupole moment Q : " << mom_quad_38si <<
" 10^38 kg m^2" 555 ost <<
"Old quadrupole moment : " <<
mom_quad_old() << endl ;
557 ost <<
"q = c^4 Q / (G^2 M^3) : " 559 ost <<
"j = c J / (G M^2) : " 563 ost <<
"Baryon mass 1 : " <<
mass_b1() / msol <<
" Msol" << endl ;
564 ost <<
"Baryon mass 2 : " <<
mass_b2() / msol <<
" Msol" << endl ;
577 double freq =
omega / (2.*M_PI) ;
578 ost <<
"Omega : " <<
omega * f_unit
579 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
580 ost <<
"Rotation period : " << 1000. / (freq * f_unit) <<
" ms" 582 ost << endl <<
"Central enthalpy : " <<
ent()(0,0,0,0) <<
" c^2" << endl ;
583 ost <<
"Central proper baryon density : " <<
nbar()(0,0,0,0)
584 <<
" x 0.1 fm^-3" << endl ;
585 double freq2 =
omega2 / (2.*M_PI) ;
586 ost <<
"Omega2 : " <<
omega2 * f_unit
587 <<
" rad/s f : " << freq2 * f_unit <<
" Hz" << endl ;
588 ost <<
"Rotation period 2: " << 1000. / (freq2 * f_unit) <<
" ms" 590 ost << endl <<
"Central enthalpy 2: " <<
ent2()(0,0,0,0) <<
" c^2" << endl ;
591 ost <<
"Central proper baryon density 2: " <<
nbar2()(0,0,0,0)
592 <<
" x 0.1 fm^-3" << endl ;
593 ost <<
"Central proper energy density : " <<
ener()(0,0,0,0)
594 <<
" rho_nuc c^2" << endl ;
595 ost <<
"Central pressure : " <<
press()(0,0,0,0)
596 <<
" rho_nuc c^2" << endl ;
598 ost <<
"Central value of log(N) : " 599 <<
logn()(0, 0, 0, 0) << endl ;
600 ost <<
"Central lapse N : " <<
nnn()(0,0,0,0) << endl ;
601 ost <<
"Central value of dzeta=log(AN) : " 602 <<
dzeta()(0, 0, 0, 0) << endl ;
603 ost <<
"Central value of A^2 : " <<
a_car()(0,0,0,0) << endl ;
604 ost <<
"Central value of B^2 : " <<
b_car()(0,0,0,0) << endl ;
606 double nphi_c =
nphi()(0, 0, 0, 0) ;
607 if ( (
omega==0) && (nphi_c==0) ) {
608 ost <<
"Central N^phi : " << nphi_c << endl ;
611 ost <<
"Central N^phi/Omega : " << nphi_c /
omega << endl ;
615 int lsurf =
nzet - 1;
618 ost <<
"Equatorial value of the velocity U: " 619 <<
uuu()(
nzet-1, 0, nt-1, nr-1) <<
" c" << endl ;
621 ost <<
"Equatorial value of the velocity U2: " 622 <<
uuu2()(
nzet-1, 0, nt-1, nr-1) <<
" c" << endl ;
625 <<
"Coordinate equatorial radius r_eq = " 626 <<
ray_eq()/km <<
" km" << endl ;
627 ost <<
"Flattening r_pole/r_eq : " <<
aplat() << endl ;
629 <<
"Coordinate equatorial radius r_eq2 = " 630 <<
ray_eq2()/km <<
" km" << endl ;
631 ost <<
"Flattening r_pole2/r_eq2 : " <<
aplat2() << endl ;
650 cout <<
"Et_rot_bifluid::equation_of_state: not ready yet for nzet > 2 !" << endl ;
653 double epsilon = 1.e-12 ;
660 for (
int l=0; l<nz; l++) {
662 for (
int k=0; k<mg->
get_np(l); k++) {
663 for (
int j=0; j<mg->
get_nt(l); j++) {
664 for (
int i=0; i<mg->
get_nr(l); i++) {
676 fact_ent.
set(0) = 1 + epsilon * xi(0) * xi(0) ;
677 fact_ent.
set(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
679 for (
int l=
nzet; l<nz; l++) {
680 fact_ent.
set(l) = 1 ;
683 ent_eos = fact_ent * ent_eos ;
684 ent2_eos = fact_ent * ent2_eos ;
686 ent2_eos.std_base_scal() ;
756 (
uuu.
set()).std_base_scal() ;
782 int j_b = mg->
get_nt(l_b) - 1 ;
784 bool superlum = false ;
786 for (
int l=0; l<
nzet; l++) {
787 for (
int i=0; i<mg->
get_nr(l); i++) {
789 double u1 =
uuu()(l, 0, j_b, i) ;
790 double u2 =
uuu2()(l, 0, j_b, i) ;
791 if ((u1 >= 1.) || (u2>=1.)) {
793 cout <<
"U > c for l, i : " << l <<
" " << i
794 <<
" U1 = " << u1 << endl ;
795 cout <<
" U2 = " << u2 << endl ;
800 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
842 E_euler = Ann + 2. * Anp + App -
press ;
virtual double z_eqf() const
Forward redshift factor at equator.
Cmp get_Kpp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 2) .
double get_m2() const
Return the individual particule mass .
double * p_mass_b2
Baryon mass of fluid 2.
virtual double r_circ() const
Circumferential equatorial radius.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
void annule(int l)
Sets the Tenseur to zero in a given domain.
void calcule_interpol(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, Cmp &K_nn, Cmp &K_np, Cmp &K_pp, Cmp &alpha_eos, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
virtual double mass_g() const
Gravitational mass.
Tenseur j_euler2
To compute .
Tenseur j_euler1
To compute .
Tenseur K_pp
Coefficient .
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Cmp sqrt(const Cmp &)
Square root.
double mass_b1() const
Baryon mass of fluid 1.
void set_enthalpies(const Cmp &, const Cmp &)
Sets both enthalpy profiles.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
double * p_r_circ2
Circumferential radius of fluid no.2.
Tbl * p_xi_surf2
Description of the surface of fluid 2: 2-D Tbl containing the values of the radial coordinate on the...
Standard units of space, time and mass.
Tenseur uuu2
Norm of the (fluid no.2) 3-velocity with respect to the eulerian observer.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tenseur nnn
Total lapse function.
virtual void equation_of_state()
Computes the proper baryon and energy densities, as well as pressure and the coefficients Knn...
Tenseur nphi
Metric coefficient .
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
double * p_angu_mom_1
Angular momentum of fluid 1.
double ray_eq() const
Coordinate radius at , [r_unit].
Base class for coordinate mappings.
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Tenseur b_car
Square of the metric factor B.
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
virtual void sauve(FILE *) const
Save in a file.
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
double * p_ray_eq2
Coordinate radius at , .
Tenseur press
Fluid pressure.
double * x
Array of values of at the nr collocation points.
Class for two-fluid rotating relativistic stars.
virtual double mean_radius() const
Mean radius.
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
virtual void del_deriv() const
Deletes all the derived quantities.
2-fluids equation of state base class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tenseur j_euler
Total angular momentum (flat-space!) 3-vector , which is related to of the "3+1" decomposition...
Cmp get_Knn(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 1) .
double * p_mass_b1
Baryon mass of fluid 1.
Tenseur gam_euler2
Lorentz factor between the fluid 2 and Eulerian observers.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void sauve(FILE *) const
Save in a file.
const Eos_bifluid & eos
Equation of state for two-fluids model.
virtual double z_pole() const
Redshift factor at North pole.
Tenseur nbar
Baryon density in the fluid frame.
double * p_coupling_mominert_1
Quantities used to describe the different couplings between the fluids.
double omega2
Rotation angular velocity for fluid 2 ([f_unit] )
virtual void sauve(FILE *) const
Save in a file.
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
double * p_coupling_mominert_2
virtual double angu_mom_2() const
Angular momentum of fluid 2.
Map & mp
Mapping associated with the star.
double * p_area2
Surface area of fluid no.2.
int get_nzone() const
Returns the number of domains.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
virtual double z_eqb() const
Backward redshift factor at equator.
int get_etat() const
Returns the logical state.
Tenseur enerps_euler
the combination : useful because in the Newtonian limit .
Tenseur bbb
Metric factor B.
double * p_ray_eq2_pi
Coordinate radius at , .
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Tenseur K_np
Coefficient .
Cmp pow(const Cmp &, int)
Power .
void operator=(const Et_rot_bifluid &)
Assignment to another Et_rot_bifluid.
virtual double mom_quad_old() const
Part of the quadrupole moment.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual double grv2() const
Error on the virial identity GRV2.
Tenseur uuu
Norm of u_euler.
Tbl & set(int l)
Read/write of the Tbl in a given domain.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
double omega
Rotation angular velocity ([f_unit] )
virtual double mean_radius2() const
Mean radius for fluid 2.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int nzet
Number of domains of *mp occupied by the star.
virtual void del_deriv() const
Deletes all the derived quantities.
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
virtual double angu_mom_1() const
Angular momentum of fluid 1.
Tenseur a_car
Total conformal factor .
double * p_angu_mom_2
Angular momentum of fluid 2.
Tenseur ent2
Log-enthalpy for the second fluid.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Tenseur sphph_euler
The component of the stress tensor .
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Tbl & set(int l)
Read/write of the value in a given domain.
virtual double mom_quad() const
Quadrupole moment.
Tenseur ener
Total energy density in the fluid frame.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
virtual void calcule_tout(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
double * p_ray_eq2_pis2
Coordinate radius at , .
double * p_ray_pole2
Coordinate radius at .
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
virtual double r_circ2() const
Circumferential radius for fluid 2.
Tenseur nbar2
Baryon density in the fluid frame, for fluid 2.
Tenseur K_nn
Coefficient .
Class for a two-fluid (tabulated) equation of state.
Tenseur & logn
Metric potential = logn_auto.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
virtual double coupling_mominert_1() const
Quantities used to study the different fluid couplings: , and .
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
virtual double angu_mom() const
Angular momentum.
double mass_b2() const
Baryon mass of fluid 2.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual ~Et_rot_bifluid()
Destructor.
Cmp get_Knp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy with respect to .
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Tenseur alpha_eos
Coefficient relative to entrainment effects.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Tenseur delta_car
The "relative velocity" (squared) of the two fluids.
Tenseur ener_euler
Total energy density in the Eulerian frame.
double ray_eq2() const
Coordinate radius for fluid 2 at , [r_unit].
Tenseur & dzeta
Metric potential = beta_auto.
virtual double aplat2() const
Flatening r_pole/r_eq for fluid 2.
double * p_aplat2
Flatening r_pole/r_eq of fluid no.2.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double get_m1() const
Return the individual particule mass .
Tensor handling *** DEPRECATED : use class Tensor instead ***.
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Et_rot_bifluid(Map &mp_i, int nzet_i, bool relat, const Eos_bifluid &eos_i)
Standard constructor.
Itbl * p_l_surf2
Description of the surface of fluid 2: 2-D Itbl containing the values of the domain index l on the su...
virtual double aplat() const
Flatening r_pole/r_eq.