32char et_bin_equil_regu_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_equil_regu.C,v 1.8 2014/10/13 08:52:55 j_novak Exp $" ;
137#include "utilitaires.h"
145 double relax_poisson,
int mermax_potvit,
146 double relax_potvit,
double thres_adapt,
147 const Tbl& fact_resize,
Tbl& diff) {
166 int i_b = mg->
get_nr(l_b) - 1 ;
176 double& diff_ent = diff.
set(0) ;
177 double& diff_vel_pot = diff.
set(1) ;
178 double& diff_logn = diff.
set(2) ;
179 double& diff_beta = diff.
set(3) ;
180 double& diff_shift_x = diff.
set(4) ;
181 double& diff_shift_y = diff.
set(5) ;
182 double& diff_shift_z = diff.
set(6) ;
194 int nz_search =
nzet ;
197 double precis_secant = 1.e-14 ;
199 double reg_map = 1. ;
203 par_adapt.
add_int(nitermax, 0) ;
208 par_adapt.
add_int(nz_search, 2) ;
210 par_adapt.
add_int(adapt_flag, 3) ;
226 par_adapt.
add_tbl(ent_limit, 0) ;
232 double precis_poisson = 1.e-16 ;
236 par_poisson1.
add_int(mermax_poisson, 0) ;
248 par_poisson2.
add_int(mermax_poisson, 0) ;
258 Param par_poisson_vect ;
260 par_poisson_vect.
add_int(mermax_poisson, 0) ;
262 par_poisson_vect.
add_double(relax_poisson, 0) ;
263 par_poisson_vect.
add_double(precis_poisson, 1) ;
285 Cmp source_regu(
mp) ;
292 for(
int mer=0 ; mer<mermax ; mer++ ) {
294 cout <<
"-----------------------------------------------" << endl ;
295 cout <<
"step: " << mer << endl ;
296 cout <<
"diff_ent = " << diff_ent << endl ;
318 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
319 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
325 double alpha_r2 = 0 ;
326 for (
int k=0; k<mg->
get_np(l_b); k++) {
327 for (
int j=0; j<mg->
get_nt(l_b); j++) {
329 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
330 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
332 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
333 ( logn_auto_b - logn_auto_c ) ;
338 if (alpha_r2_jk > alpha_r2) {
339 alpha_r2 = alpha_r2_jk ;
347 alpha_r =
sqrt(alpha_r2) ;
349 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
372 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
383 double dent_eq =
ent().dsdr().val_point(
ray_eq(),M_PI/2.,0.) ;
384 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
385 double rap_dent = fabs( dent_eq / dent_pole ) ;
386 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
388 if ( rap_dent < thres_adapt ) {
390 cout <<
"******* FROZEN MAPPING *********" << endl ;
399 for (
int l=0; l<
nzet; l++) {
400 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
402 ent_limit.
set(
nzet-1) = ent_b ;
415 double rr_in_1 =
mp.
val_r(
nzet, -1., M_PI/2., 0.) ;
418 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
420 mp.
resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
423 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
425 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
430 double rr_out_nm3_new =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
432 for (
int i=
nzet-1; i<nz-4; i++) {
434 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
436 double rr_mid = rr_out_i
437 + (rr_out_nm3_new - rr_out_i) /
double(nz - 3 - i) ;
439 double rr_2timesi = 2. * rr_out_i ;
441 if (rr_2timesi < rr_mid) {
443 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
445 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
450 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
452 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
473 double ent_s_max = -1 ;
476 for (
int k=0; k<mg->
get_np(l_b); k++) {
477 for (
int j=0; j<mg->
get_nt(l_b); j++) {
478 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
479 if (xx > ent_s_max) {
486 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1"
487 <<
" and nzet : " << endl ;
488 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max <<
489 " and j = " << j_s_max << endl ;
526 source = qpig *
nbar ;
544 source().poisson_regular(
k_div,
nzet, unsgam1, par_poisson1,
548 source_regu, source_div) ;
555 "Relative error in the resolution of the equation for logn_auto : "
557 for (
int l=0; l<nz; l++) {
558 cout << tdiff_logn(l) <<
" " ;
561 diff_logn =
max(
abs(tdiff_logn)) ;
593 "Relative error in the resolution of the equation for beta_auto : "
595 for (
int l=0; l<nz; l++) {
596 cout << tdiff_beta(l) <<
" " ;
599 diff_beta =
max(
abs(tdiff_beta)) ;
622 for (
int i=0; i<3; i++) {
624 if (source_shift(i).get_etat() != ETATZERO)
635 for (
int i=0; i<3; i++) {
636 if(source_shift(i).dz_nonzero()) {
637 assert( source_shift(i).get_dzpuis() == 4 ) ;
640 (source_shift.
set(i)).set_dzpuis(4) ;
647 double lambda_shift = double(1) / double(3) ;
650 lambda_shift, par_poisson_vect,
668 for (
int i=0; i<3; i++) {
670 + lambda_shift * graddivn(i) ;
673 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
674 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
675 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
678 "Relative error in the resolution of the equation "
681 cout <<
"x component : " ;
682 for (
int l=0; l<nz; l++) {
683 cout << tdiff_shift_x(l) <<
" " ;
686 cout <<
"y component : " ;
687 for (
int l=0; l<nz; l++) {
688 cout << tdiff_shift_y(l) <<
" " ;
691 cout <<
"z component : " ;
692 for (
int l=0; l<nz; l++) {
693 cout << tdiff_shift_z(l) <<
" " ;
697 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
698 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
699 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
718 diff_ent = diff_ent_tbl(0) ;
719 for (
int l=1; l<
nzet; l++) {
720 diff_ent += diff_ent_tbl(l) ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy with extra parameters (virtual function im...
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy and extra parameters (virtual function imp...
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Tenseur akcar_comp
Part of the scalar generated by shift_auto and shift_comp , i.e.
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
bool irrotational
true for an irrotational star, false for a corotating one
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Tenseur pot_centri
Centrifugal potential.
void equil_regular(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration by regularizing the diverging source.
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Tenseur d_beta_comp
Gradient of beta_comp (Cartesian components with respect to ref_triad )
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by t...
int nzet
Number of domains of *mp occupied by the star.
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
int k_div
Index of regularity of the gravitational potential logn_auto .
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Tenseur nbar
Baryon density in the fluid frame.
const Eos & eos
Equation of state of the stellar matter.
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
Tenseur press
Fluid pressure.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Tenseur a_car
Total conformal factor .
double ray_pole() const
Coordinate radius at [r_unit].
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 resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
virtual void reevaluate_symy(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.
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
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.
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 dec2_dzpuis()
dzpuis -= 2 ;
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void inc2_dzpuis()
dzpuis += 2 ;
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void poisson_vect_regu(int k_div, int nzet, double unsgam1, double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Cmp sqrt(const Cmp &)
Square root.
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 abs(const Cmp &)
Absolute value.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Standard units of space, time and mass.