29 char boson_star_equil_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Compobj/boson_star_equil.C,v 1.6 2014/10/13 08:52:49 j_novak Exp $" ;
61 #include "boson_star.h"
65 #include "graphique.h"
66 #include "utilitaires.h"
71 int nzadapt,
const Tbl& phi_limit,
const Itbl& icontrol,
81 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
82 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
99 int j_b = mg->
get_nt(l_b) - 1 ;
108 int mer_max = icontrol(0) ;
113 int mermax_poisson = icontrol(5) ;
114 int mer_triax = icontrol(6) ;
118 double precis = control(0) ;
120 double relax = control(2) ;
121 double relax_prev = double(1) - relax ;
122 double relax_poisson = control(3) ;
124 double ampli_triax = control(5) ;
125 double precis_adapt = control(6) ;
132 double& diff_phi = diff.
set(0) ;
133 double& diff_nuf = diff.
set(1) ;
134 double& diff_nuq = diff.
set(2) ;
137 double& diff_shift_x = diff.
set(5) ;
138 double& diff_shift_y = diff.
set(6) ;
139 double& vit_triax = diff.
set(7) ;
150 int nz_search = nzet + 1 ;
153 double reg_map = 1. ;
155 par_adapt.
add_int(nitermax, 0) ;
157 par_adapt.
add_int(nzadapt, 1) ;
160 par_adapt.
add_int(nz_search, 2) ;
162 par_adapt.
add_int(adapt_flag, 3) ;
181 par_adapt.
add_tbl(phi_limit, 0) ;
187 double precis_poisson = 1.e-16 ;
202 for (
int i=1; i<=3; i++) {
208 for (
int i=1; i<=3; i++) {
214 for (
int i=1; i<=3; i++) {
222 Param par_poisson_nuf ;
223 par_poisson_nuf.
add_int(mermax_poisson, 0) ;
224 par_poisson_nuf.
add_double(relax_poisson, 0) ;
225 par_poisson_nuf.
add_double(precis_poisson, 1) ;
229 Param par_poisson_nuq ;
230 par_poisson_nuq.
add_int(mermax_poisson, 0) ;
231 par_poisson_nuq.
add_double(relax_poisson, 0) ;
232 par_poisson_nuq.
add_double(precis_poisson, 1) ;
236 Param par_poisson_tggg ;
237 par_poisson_tggg.
add_int(mermax_poisson, 0) ;
238 par_poisson_tggg.
add_double(relax_poisson, 0) ;
239 par_poisson_tggg.
add_double(precis_poisson, 1) ;
245 Param par_poisson_dzeta ;
252 Param par_poisson_vect ;
254 par_poisson_vect.
add_int(mermax_poisson, 0) ;
255 par_poisson_vect.
add_double(relax_poisson, 0) ;
256 par_poisson_vect.
add_double(precis_poisson, 1) ;
288 ofstream fichconv(
"convergence.d") ;
289 fichconv <<
"# diff_phi GRV2 max_triax vit_triax" << endl ;
292 ofstream fichevol(
"evolution.d") ;
294 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq rphi_c"
298 double err_grv2 = 1 ;
299 double max_triax_prev = 0 ;
305 for(
int mer=0 ; (diff_phi > precis) && (mer<mer_max) ; mer++ ) {
307 cout <<
"-----------------------------------------------" << endl ;
308 cout <<
"step: " << mer << endl ;
309 cout <<
"diff_phi = " << display_bold << diff_phi << display_normal
311 cout <<
"err_grv2 = " << err_grv2 << endl ;
331 source_nuq =
ak_car - d_logn(1)*(d_logn(1)+d_bet(1))
332 - d_logn(2)*(d_logn(2)+d_bet(2))
333 - d_logn(3)*(d_logn(3)+d_bet(3)) ;
336 source_nuq.std_spectral_base() ;
344 - d_logn(1)*d_logn(1) - d_logn(2)*d_logn(2) - d_logn(3)*d_logn(3) ;
362 source_shift = (-4*qpig) *
nn *
a_car * mom_euler_cart ;
379 cout <<
"Test of the Poisson equation for nuf :" << endl ;
381 diff_nuf = err(0, 0) ;
387 if (mer == mer_triax) {
389 if ( mg->
get_np(0) == 1 ) {
391 "Boson_star::equilibrium: np must be stricly greater than 1"
392 << endl <<
" to set a triaxial perturbation !" << endl ;
399 perturb = 1 + ampli_triax * sint*sint *
cos(2*phi) ;
413 double max_triax = 0 ;
415 if ( mg->
get_np(0) > 1 ) {
417 for (
int l=0; l<nz; l++) {
418 for (
int j=0; j<mg->
get_nt(l); j++) {
419 for (
int i=0; i<mg->
get_nr(l); i++) {
422 double xcos2p = (*(va_nuf.
c_cf))(l, 2, j, i) ;
425 double xsin2p = (*(va_nuf.
c_cf))(l, 3, j, i) ;
427 double xx =
sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
429 max_triax = ( xx > max_triax ) ? xx : max_triax ;
436 cout <<
"Triaxial part of nuf : " << max_triax << endl ;
442 source_nuq.poisson(par_poisson_nuq,
nuq) ;
444 cout <<
"Test of the Poisson equation for nuq :" << endl ;
445 err = source_nuq.test_poisson(
nuq, cout,
true) ;
446 diff_nuq = err(0, 0) ;
453 for (
int i=1; i<=3; i++) {
454 if(source_shift(i).get_etat() != ETATZERO) {
455 if(source_shift(i).dz_nonzero()) {
456 assert( source_shift(i).get_dzpuis() == 4 ) ;
459 (source_shift.
set(i)).set_dzpuis(4) ;
464 double lambda_shift = double(1) / double(3) ;
466 if ( mg->
get_np(0) == 1 ) {
472 for (
int i=1; i<=3; i++) {
473 csource_shift.
set(i-1) = source_shift(i) ;
477 csource_shift.
poisson_vect(lambda_shift, par_poisson_vect,
478 cshift, cw_shift, ckhi_shift) ;
480 for (
int i=1; i<=3; i++) {
487 cout <<
"Test of the Poisson equation for shift_x :" << endl ;
488 err = source_shift(1).test_poisson(-
beta(1), cout,
true) ;
489 diff_shift_x = err(0, 0) ;
491 cout <<
"Test of the Poisson equation for shift_y :" << endl ;
492 err = source_shift(2).test_poisson(-
beta(2), cout,
true) ;
493 diff_shift_y = err(0, 0) ;
532 Cmp csource_tggg(source_tggg) ;
543 Cmp csource_dzf(source_dzf) ;
544 Cmp csource_dzq(source_dzq) ;
546 mp.
poisson2d(csource_dzf, csource_dzq, par_poisson_dzeta,
550 err_grv2 = lbda_grv2 - 1;
551 cout <<
"GRV2: " << err_grv2 << endl ;
561 logn = relax *
logn + relax_prev * logn_prev ;
563 dzeta = relax *
dzeta + relax_prev * dzeta_prev ;
578 cout << *
this << endl ;
586 diff_phi = diff_phi_tbl(0) ;
587 for (
int l=1; l<nzet; l++) {
588 diff_phi += diff_phi_tbl(l) ;
592 fichconv <<
" " <<
log10( fabs(diff_phi) + 1.e-16 ) ;
593 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
594 fichconv <<
" " <<
log10( fabs(max_triax) + 1.e-16 ) ;
597 if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
598 vit_triax = (max_triax - max_triax_prev) / max_triax_prev ;
601 fichconv <<
" " << vit_triax ;
610 max_triax_prev = max_triax ;
627 for (
int i=1; i<=3; i++) {
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Radial mapping of rather general form.
double & set(int i)
Read/write of a particular element (index i) (1D case)
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Scalar dzeta
Metric potential .
Values and coefficients of a (real-value) function.
Scalar bbb
Metric factor B.
Map & mp
Mapping describing the coordinate system (r,theta,phi)
virtual void equilibrium(double rphi_c, double iphi_c, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, Tbl &diff, Param *=0x0)
Solves the equation satisfied by the scalar field.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Cmp log(const Cmp &)
Neperian logarithm.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar rphi
Real part of the scalar field Phi.
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.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Basic integer array class.
Vector beta
Shift vector .
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Tensor field of valence 0 (or component of a tensorial field).
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void coef() const
Computes the coeffcients of *this.
const Valeur & get_spectral_va() const
Returns va (read only version)
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Standard units of space, time and mass.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Sym_tensor kk
Extrinsic curvature tensor
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Active physical coordinates and mapping derivatives.
Scalar logn
Logarithm of the lapse N .
void update_metric()
Computes metric coefficients from known potentials.
Scalar & set(int)
Read/write access to a component.
Scalar ener_euler
Total energy density E in the Eulerian frame.
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.
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
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_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Cmp log10(const Cmp &)
Basis 10 logarithm.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord phi
coordinate centered on the grid
Cmp sqrt(const Cmp &)
Square root.
void set_dzpuis(int)
Modifies the dzpuis flag.
Tensor field of valence 1.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
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.
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Cmp cos(const Cmp &)
Cosine.
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Scalar b_car
Square of the metric factor B.
Scalar tggg
Metric potential .
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
void update_ener_mom()
Computes the 3+1 components of the energy-momentum tensor (E, P_i and S_{ij}) from the values of the ...
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Scalar a_car
Square of the metric factor A.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Scalar nn
Lapse function N .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.