29 char star_bhns_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $" ;
58 #include "star_bhns.h"
62 #include "utilitaires.h"
68 const double& sepa,
bool kerrschild,
69 int mer,
int mermax_ns,
int mermax_potvit,
70 int mermax_poisson,
int filter_r,
71 int filter_r_s,
int filter_p_s,
72 double relax_poisson,
double relax_potvit,
73 double thres_adapt,
double resize_ns,
74 const Tbl& fact_resize,
Tbl& diff) {
92 int i_b = mg->
get_nr(l_b) - 1 ;
102 double& diff_ent = diff.
set(0) ;
103 double& diff_vel_pot = diff.
set(1) ;
104 double& diff_lapconf = diff.
set(2) ;
105 double& diff_confo = diff.
set(3) ;
106 double& diff_shift_x = diff.
set(4) ;
107 double& diff_shift_y = diff.
set(5) ;
108 double& diff_shift_z = diff.
set(6) ;
109 double& diff_dHdr = diff.
set(7) ;
110 double& diff_dHdr_min = diff.
set(8) ;
111 double& diff_phi_min = diff.
set(9) ;
112 double& diff_radius = diff.
set(10) ;
125 int nz_search =
nzet + 1 ;
128 double precis_secant = 1.e-14 ;
130 double reg_map = 1. ;
134 par_adapt.
add_int(nitermax, 0) ;
139 par_adapt.
add_int(nz_search, 2) ;
141 par_adapt.
add_int(adapt_flag, 3) ;
157 par_adapt.
add_tbl(ent_limit, 0) ;
166 double precis_poisson = 1.e-14 ;
173 par_lapconf.
add_int(mermax_poisson, 0) ;
184 par_confo.
add_int(mermax_poisson, 0) ;
195 par_shift2.
add_int(mermax_poisson, 0) ;
205 for (
int i=0; i<3; i++) {
229 for (
int mer_ns=0; mer_ns<mermax_ns; mer_ns++) {
231 cout <<
"-----------------------------------------------" << endl ;
232 cout <<
"step: " << mer_ns << endl ;
233 cout <<
"diff_ent = " << diff_ent << endl ;
242 mermax_potvit, precis_poisson,
267 double hhh_c =
exp(ent_c) ;
268 double hhh_b =
exp(ent_b) ;
274 double alpha_r2 = 0. ;
276 for (
int k=0; k<mg->
get_np(l_b); k++) {
277 for (
int j=0; j<mg->
get_nt(l_b); j++) {
279 double lapconf_auto_b =
281 double lapconf_comp_b =
289 double aaa = (gam0_c*gam_b*hhh_b*confo_c)
290 / (gam0_b*gam_c*hhh_c*confo_b) ;
293 double alpha_r2_jk = (aaa * lapconf_comp_b - lapconf_comp_c
295 / (lapconf_auto_c - aaa * lapconf_auto_b ) ;
297 if (alpha_r2_jk > alpha_r2) {
298 alpha_r2 = alpha_r2_jk ;
305 alpha_r =
sqrt(alpha_r2) ;
307 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
334 double log_lapconf_c =
log(lapconf_tot_tmp.val_grid_point(0,0,0,0)) ;
339 ent = (ent_c + log_lapconf_c - log_confo_c + loggam_c + pot_centri_c)
352 cout <<
"dH/dx|_center = " << dentdx << endl ;
353 cout <<
"dH/dy|_center = " << dentdy << endl ;
355 double dec_fact_x = 1. ;
356 double dec_fact_y = 1. ;
359 func_in = 1. - dec_fact_x * (dentdx/ent_c) *
mp.
x
360 - dec_fact_y * (dentdy/ent_c) *
mp.
y ;
372 ent =
ent * (func_in + func_ex) ;
378 cout <<
"dH/dx|_new = " << dentdx_new << endl ;
379 cout <<
"dH/dy|_new = " << dentdy_new << endl ;
387 double rap_dent = fabs( dent_eq / dent_pole ) ;
388 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
389 diff_dHdr = rap_dent ;
391 if ( rap_dent < thres_adapt ) {
393 cout <<
"******* FROZEN MAPPING *********" << endl ;
401 for (
int l=0; l<
nzet; l++) {
404 ent_limit.
set(
nzet-1) = ent_b ;
415 double rr_in_1 =
mp.
val_r(1, -1., M_PI/2., 0.) ;
418 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
419 mp.
resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
422 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
423 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
430 double rr_out_1 =
mp.
val_r(1, 1., M_PI/2., 0.) ;
431 mp.
resize(1, rr_in_1/rr_out_1 * resize_ns) ;
436 double rr_out_nm3_new =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
438 for (
int i=1; i<nz-4; i++) {
440 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
442 double rr_mid = rr_out_i
443 + (rr_out_nm3_new - rr_out_i) /
double(nz - 3 - i) ;
445 double rr_2timesi = 2. * rr_out_i ;
447 if (rr_2timesi < rr_mid) {
449 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
450 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
455 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
456 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
478 double ent_s_max = -1 ;
482 for (
int k=0; k<mg->
get_np(l_b); k++) {
483 for (
int j=0; j<mg->
get_nt(l_b); j++) {
485 if (xx > ent_s_max) {
492 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1"
493 <<
" and nzet : " << endl ;
494 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max
495 <<
" and j = " << j_s_max << endl ;
517 double rad_chi_min =
radius_p(azimu_min) ;
518 double chi_min =
chi_rp(rad_chi_min, azimu_min) ;
520 cout <<
"| dH/dr_eq / dH/dr_pole | (minimum) = " << chi_min << endl ;
521 cout <<
" phi = " << azimu_min / M_PI <<
" [M_PI]" << endl ;
522 cout <<
" radius = " << rad_chi_min / km <<
" [km]" << endl ;
524 diff_dHdr_min = chi_min ;
525 diff_phi_min = azimu_min ;
526 diff_radius = rad_chi_min ;
548 source_lapconf = sou_lap1 + sou_lap2 ;
554 if (source_lapconf.
get_etat() != ETATZERO) {
555 source_lapconf.
filtre(filter_r) ;
567 source_lapconf.
poisson(par_lapconf, lapconf_m1) ;
575 "Relative error in the resolution of the equation for lapconf_auto : "
577 for (
int l=0; l<nz; l++) {
578 cout << tdiff_lapconf(l) <<
" " ;
581 diff_lapconf =
max(
abs(tdiff_lapconf)) ;
607 source_confo = sou_con1 + sou_con2 ;
613 if (source_confo.
get_etat() != ETATZERO) {
614 source_confo.
filtre(filter_r) ;
626 source_confo.
poisson(par_confo, confo_m1) ;
634 "Relative error in the resolution of the equation for confo_auto : "
636 for (
int l=0; l<nz; l++) {
637 cout << tdiff_confo(l) <<
" " ;
640 diff_confo =
max(
abs(tdiff_confo)) ;
658 for (
int i=1; i<=3; i++) {
659 sou_shif1.
set(i) = 4.*qpig * lapconf_tot_tmp
667 for (
int i=1; i<=3; i++) {
668 (sou_shif1.
set(i)).inc_dzpuis(4) ;
673 for (
int i=1; i<=3; i++) {
674 sou_shif2.
set(i) = 2. *
688 source_shift = sou_shif1 + sou_shif2 ;
698 if (filter_r_s != 0) {
699 for (
int i=1; i<=3; i++) {
700 if (source_shift(i).get_etat() != ETATZERO) {
707 if (filter_p_s != 0) {
708 for (
int i=1; i<=3; i++) {
709 if (source_shift(i).get_etat() != ETATZERO) {
710 (source_shift.
set(i)).filtre_phi(filter_p_s, nz-1) ;
716 for (
int i=1; i<=3; i++) {
717 if(source_shift(i).dz_nonzero()) {
718 assert( source_shift(i).get_dzpuis() == 4 ) ;
721 (source_shift.
set(i)).set_dzpuis(4) ;
725 double lambda = double(1) / double(3) ;
729 for (
int i=0; i<3; i++) {
730 source_p.
set(i) =
Cmp(source_shift(i+1)) ;
735 for (
int i=0; i<3 ;i++) {
736 vect_auxi.
set(i) = 0. ;
746 source_p.
poisson_vect(lambda, par_shift2, resu_p, vect_auxi,
749 for (
int i=1; i<=3; i++)
754 for (
int i=0; i<3; i++) {
766 Tbl tdiff_shift_x =
diffrel(lap_shift(1), source_shift(1)) ;
767 Tbl tdiff_shift_y =
diffrel(lap_shift(2), source_shift(2)) ;
768 Tbl tdiff_shift_z =
diffrel(lap_shift(3), source_shift(3)) ;
771 "Relative error in the resolution of the equation for shift_auto : "
773 cout <<
"x component : " ;
774 for (
int l=0; l<nz; l++) {
775 cout << tdiff_shift_x(l) <<
" " ;
778 cout <<
"y component : " ;
779 for (
int l=0; l<nz; l++) {
780 cout << tdiff_shift_y(l) <<
" " ;
783 cout <<
"z component : " ;
784 for (
int l=0; l<nz; l++) {
785 cout << tdiff_shift_z(l) <<
" " ;
789 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
790 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
791 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
799 diff_ent = diff_ent_tbl(0) ;
800 for (
int l=0; l<
nzet; l++) {
801 diff_ent += diff_ent_tbl(l) ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void annule_hard()
Sets the Cmp to zero in a hard way.
Radial mapping of rather general form.
virtual void homothetie(double lambda)
Sets a new radial scale.
Coord y
y coordinate centered on the grid
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 resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Coord x
x coordinate centered on the grid
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
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_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).
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
int get_dzpuis() const
Returns dzpuis.
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
const Scalar & dsdy() const
Returns of *this , where .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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 inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & dsdx() const
Returns of *this , where .
const Scalar & dsdr() const
Returns of *this .
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Valeur & set_spectral_va()
Returns va (read/write version)
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star ).
Sym_tensor taij_auto
Part of the extrinsic curvature tensor generated by shift_auto , lapse_auto , and confo_auto .
Scalar lapconf_auto
Lapconf function generated by the star.
double velo_pot_bhns(const double &mass_bh, const double &sepa, bool kerrschild, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
void equilibrium_bhns(double ent_c, const double &mass_bh, const double &sepa, bool kerrschild, int mer, int mermax_ns, int mermax_potvit, int mermax_poisson, int filter_r, int filter_r_s, int filter_p_s, double relax_poisson, double relax_potvit, double thres_adapt, double resize_ns, const Tbl &fact_resize, Tbl &diff)
Computes an equilibrium configuration.
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Scalar confo_tot
Total conformal factor.
Vector d_lapconf_auto
Derivative of the lapconf function generated by the star .
Vector d_confo_auto
Derivative of the conformal factor generated by the star .
double radius_p(double phi)
Radius of the star to the direction of and .
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
void hydro_euler_bhns(bool kerrschild, const double &mass_bh, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Scalar gam
Lorentz factor between the fluid and the co-orbiting observer.
Scalar ssjm1_lapconf
Effective source at the previous step for the resolution of the Poisson equation for lapconf_auto .
bool irrotational
true for an irrotational star, false for a corotating one
Scalar pot_centri
Centrifugal potential.
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Scalar lapconf_comp
Lapconf function generated by the companion black hole.
Scalar taij_quad_auto
Part of the scalar generated by .
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Vector shift_auto
Shift vector generated by the star.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Scalar confo_auto
Conformal factor generated by the star.
Scalar ssjm1_confo
Effective source at the previous step for the resolution of the Poisson equation for confo_auto .
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Scalar ener_euler
Total energy density in the Eulerian frame.
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Scalar press
Fluid pressure.
double ray_eq_pi() const
Coordinate radius at , [r_unit].
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].
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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.
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
void smooth(int nzet, Valeur &uuva) const
Changes the function *this as a smooth one when there exists a discontinuity between the nucleus and ...
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Scalar & set(int)
Read/write access to a component.
Cmp sqrt(const Cmp &)
Square root.
Cmp exp(const Cmp &)
Exponential.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
Cmp abs(const Cmp &)
Absolute value.
Cmp log(const Cmp &)
Neperian logarithm.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
const Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Standard units of space, time and mass.