29char 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"
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 rphi
Real part of the scalar field 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.
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 ...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Scalar b_car
Square of the metric factor B.
Scalar bbb
Metric factor B.
Scalar a_car
Square of the metric factor A.
Sym_tensor kk
Extrinsic curvature tensor
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Scalar ener_euler
Total energy density E in the Eulerian frame.
Scalar nn
Lapse function N .
Vector beta
Shift vector .
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Active physical coordinates and mapping derivatives.
Basic integer array class.
Radial mapping of rather general form.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
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.
Coord phi
coordinate centered on the grid
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.
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 Valeur & get_spectral_va() const
Returns va (read only version)
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...
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Scalar logn
Logarithm of the lapse N .
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
void update_metric()
Computes metric coefficients from known potentials.
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Scalar tggg
Metric potential .
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Scalar dzeta
Metric potential .
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 : .
Values and coefficients of a (real-value) function.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
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 cos(const Cmp &)
Cosine.
Cmp log(const Cmp &)
Neperian logarithm.
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Standard units of space, time and mass.