LORENE
star_equil_spher.C
1 /*
2  * Method of class Star to compute a static spherical configuration.
3  *
4  * (see file star.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Francois Limousin
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: star_equil_spher.C,v 1.17 2017/10/20 13:54:54 j_novak Exp $
34  * $Log: star_equil_spher.C,v $
35  * Revision 1.17 2017/10/20 13:54:54 j_novak
36  * Adaptation of the method to be used within nrotstar.
37  *
38  * Revision 1.16 2016/12/05 16:18:15 j_novak
39  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40  *
41  * Revision 1.15 2014/10/13 08:53:39 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.14 2010/10/18 20:16:10 m_bejger
45  * Commented-out the reevaluation of the mapping for the case of many domains inside the star
46  *
47  * Revision 1.13 2010/10/18 18:56:33 m_bejger
48  * Changed to allow initial data with more than one domain in the star
49  *
50  * Revision 1.12 2005/09/14 12:30:52 f_limousin
51  * Saving of fields lnq and logn in class Star.
52  *
53  * Revision 1.11 2005/09/13 19:38:31 f_limousin
54  * Reintroduction of the resolution of the equations in cartesian coordinates.
55  *
56  * Revision 1.10 2005/02/17 17:32:35 f_limousin
57  * Change the name of some quantities to be consistent with other classes
58  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
59  *
60  * Revision 1.9 2004/06/22 12:45:31 f_limousin
61  * Improve the convergence
62  *
63  * Revision 1.8 2004/06/07 16:27:47 f_limousin
64  * Add the computation of virial error.
65  *
66  * Revision 1.6 2004/04/08 16:33:42 f_limousin
67  * The new variable is ln(Q) instead of Q=psi^2*N. It improves the
68  * convergence of the code.
69  *
70  * Revision 1.5 2004/03/25 10:29:26 j_novak
71  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
72  *
73  * Revision 1.4 2004/02/27 09:57:55 f_limousin
74  * We now update the metric gamma at the end of this routine for
75  * the calculus of mass_b and mass_g.
76  *
77  * Revision 1.3 2004/02/21 17:05:13 e_gourgoulhon
78  * Method Scalar::point renamed Scalar::val_grid_point.
79  * Method Scalar::set_point renamed Scalar::set_grid_point.
80  *
81  * Revision 1.2 2004/01/20 15:20:35 f_limousin
82  * First version
83  *
84  *
85  * $Header: /cvsroot/Lorene/C++/Source/Star/star_equil_spher.C,v 1.17 2017/10/20 13:54:54 j_novak Exp $
86  *
87  */
88 
89 // Headers C
90 #include "math.h"
91 
92 // Headers Lorene
93 #include "tenseur.h"
94 #include "star_rot.h"
95 #include "param.h"
96 #include "graphique.h"
97 #include "nbr_spx.h"
98 #include "unites.h"
99 
100 namespace Lorene {
101 void Star::equilibrium_spher(double ent_c, double precis, const Tbl* pent_limit){
102 
103  // Fundamental constants and units
104  // -------------------------------
105  using namespace Unites ;
106 
107  // Initializations
108  // ---------------
109 
110  const Mg3d* mg = mp.get_mg() ;
111  int nz = mg->get_nzone() ; // total number of domains
112 
113  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
114  int l_b = nzet - 1 ;
115  int i_b = mg->get_nr(l_b) - 1 ;
116  int j_b = mg->get_nt(l_b) - 1 ;
117  int k_b = 0 ;
118 
119  // Value of the enthalpy defining the surface of the star
120  double ent_b = 0 ;
121 
122  // Initialization of the enthalpy field to the constant value ent_c :
123 
124  ent = ent_c ;
125  ent.annule(nzet, nz-1) ;
126 
127 
128  // Corresponding profiles of baryon density, energy density and pressure
129 
130  equation_of_state() ;
131 
132  Scalar a_car(mp) ;
133  a_car = 1 ;
134  lnq = 0 ;
135  lnq.std_spectral_base() ;
136 
137  // Auxiliary quantities
138  // --------------------
139 
140  // Affine mapping for solving the Poisson equations
141  Map_af mpaff(mp);
142 
143  Param par_nul ; // Param (null) for Map_af::poisson.
144 
145  Scalar ent_jm1(mp) ; // Enthalpy at previous step
146  ent_jm1 = 0 ;
147 
148  Scalar source(mp) ;
149  Scalar logn_quad(mp) ;
150  Scalar logn_mat(mp) ;
151  logn_mat = 0 ;
152  logn = 0 ;
153  logn_quad = 0 ;
154 
155  Scalar dlogn(mp) ;
156  Scalar dlnq(mp) ;
157 
158  double diff_ent = 1 ;
159  int mermax = 200 ; // Max number of iterations
160 
161  //=========================================================================
162  // Start of iteration
163  //=========================================================================
164 
165  for(int mer=0 ; (diff_ent > precis) && (mer<mermax) ; mer++ ) {
166 
167  double alpha_r = 1 ;
168 
169  cout << "-----------------------------------------------" << endl ;
170  cout << "step: " << mer << endl ;
171  cout << "alpha_r: " << alpha_r << endl ;
172  cout << "diff_ent = " << diff_ent << endl ;
173 
174  //-----------------------------------------------------
175  // Resolution of Poisson equation for ln(N)
176  //-----------------------------------------------------
177 
178  // Matter part of ln(N)
179  // --------------------
180 
181  source = a_car * (ener + 3.*press) ;
182 
183  source.set_dzpuis(4) ;
184 
185  Cmp source_logn_mat (source) ;
186  Cmp logn_mat_cmp (logn_mat) ;
187  logn_mat_cmp.set_etat_qcq() ;
188 
189  mpaff.poisson(source_logn_mat, par_nul, logn_mat_cmp) ;
190 
191  logn_mat = logn_mat_cmp ;
192 
193  // NB: at this stage logn is in machine units, not in c^2
194 
195  // Quadratic part of ln(N)
196  // -----------------------
197 
198  mpaff.dsdr(logn, dlogn) ;
199  mpaff.dsdr(lnq, dlnq) ;
200 
201  source = - dlogn * dlnq ;
202 
203  Cmp source_logn_quad (source) ;
204  Cmp logn_quad_cmp (logn_quad) ;
205  logn_quad_cmp.set_etat_qcq() ;
206 
207  mpaff.poisson(source_logn_quad, par_nul, logn_quad_cmp) ;
208 
209  logn_quad = logn_quad_cmp ;
210 
211 
212  //-----------------------------------------------------
213  // Computation of the new radial scale
214  //-----------------------------------------------------
215 
216  // alpha_r (r = alpha_r r') is determined so that the enthalpy
217  // takes the requested value ent_b at the stellar surface
218 
219  double nu_mat0_b = logn_mat.val_grid_point(l_b, k_b, j_b, i_b) ;
220  double nu_mat0_c = logn_mat.val_grid_point(0, 0, 0, 0) ;
221 
222  double nu_quad0_b = logn_quad.val_grid_point(l_b, k_b, j_b, i_b) ;
223  double nu_quad0_c = logn_quad.val_grid_point(0, 0, 0, 0) ;
224 
225  double alpha_r2 = ( ent_c - ent_b - nu_quad0_b + nu_quad0_c )
226  / ( qpig*(nu_mat0_b - nu_mat0_c) ) ;
227 
228  alpha_r = sqrt(alpha_r2) ;
229 
230  // One domain inside the star:
231  // ---------------------------
232  if(nzet==1) mpaff.homothetie( alpha_r ) ;
233 
234  //--------------------
235  // First integral
236  //--------------------
237 
238  // Gravitation potential in units c^2 :
239  logn_mat = alpha_r2*qpig * logn_mat ;
240  logn = logn_mat + logn_quad ;
241 
242  // Enthalpy in all space
243  double logn_c = logn.val_grid_point(0, 0, 0, 0) ;
244  ent = ent_c - logn + logn_c ;
245 
246  // Two or more domains inside the star:
247  // ------------------------------------
248 
249  if(nzet>1) {
250 
251  // Parameters for the function Map_et::adapt
252  // -----------------------------------------
253 
254  Param par_adapt ;
255  int nitermax = 100 ;
256  int niter ;
257  int adapt_flag = 1 ; // 1 = performs the full computation,
258  // 0 = performs only the rescaling by
259  // the factor alpha_r
260  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
261  // isosurfaces
262 
263  int nzadapt = nzet ;
264 
265  //cout << "no. of domains where the ent adjustment will be done: " << nzet << endl ;
266  //cout << "ent limits: " << ent_limit << endl ;
267 
268  double precis_adapt = 1.e-14 ;
269 
270  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
271 
272  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
273  // locate zeros by the secant method
274  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
275  // to the isosurfaces of ent is to be
276  // performed
277  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
278  // the enthalpy isosurface
279  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
280  // 0 = performs only the rescaling by
281  // the factor alpha_r
282  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
283  // (theta_*, phi_*)
284  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
285  // (theta_*, phi_*)
286 
287  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
288  // the secant method
289 
290  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
291  // the determination of zeros by
292  // the secant method
293  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
294 
295  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
296  // distances will be multiplied
297 
298  // Enthalpy values for the adaptation
299  Tbl ent_limit(nzet) ;
300  if (pent_limit != 0x0) ent_limit = *pent_limit ;
301 
302  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
303  // to define the isosurfaces.
304 
305  double* bornes = new double[nz+1] ;
306  bornes[0] = 0. ;
307 
308  for(int l=0; l<nz; l++) {
309 
310  bornes[l+1] = mpaff.get_alpha()[l] + mpaff.get_beta()[l] ;
311 
312  }
313  bornes[nz] = __infinity ;
314 
315  Map_et mp0(*mg, bornes) ;
316 
317  mp0 = mpaff;
318  mp0.adapt(ent, par_adapt) ;
319 
320  //Map_af mpaff_prev (mpaff) ;
321 
322  double alphal, betal ;
323 
324  for(int l=0; l<nz; l++) {
325 
326  alphal = mp0.get_alpha()[l] ;
327  betal = mp0.get_beta()[l] ;
328 
329  mpaff.set_alpha(alphal, l) ;
330  mpaff.set_beta(betal, l) ;
331 
332  }
333 
334 
335  // testing printout
336  //int num_r1 = mg->get_nr(0) - 1;
337  //cout << "Pressure difference:" << press.val_grid_point(0,0,0,num_r1)
338  //- press.val_grid_point(1,0,0,0) << endl ;
339  //cout << "Difference in enthalpies at the domain boundary:" << endl ;
340  //cout << ent.val_grid_point(0,0,0,num_r1) << endl ;
341  //cout << ent.val_grid_point(1,0,0,0) << endl ;
342 
343  //cout << "Enthalpy difference: " << ent.val_grid_point(0,0,0,num_r1)
344  //- ent.val_grid_point(1,0,0,0) << endl ;
345 
346  // Computation of the enthalpy at the new grid points
347  //----------------------------------------------------
348  //mpaff.reevaluate(&mpaff_prev, nzet+1, ent) ;
349 
350  }
351 
352 
353  //---------------------
354  // Equation of state
355  //---------------------
356 
357  equation_of_state() ;
358 
359  //---------------------
360  // Equation for lnq_auto
361  //---------------------
362 
363  mpaff.dsdr(logn, dlogn) ;
364  mpaff.dsdr(lnq, dlnq) ;
365 
366  source = 3 * qpig * a_car * press ;
367 
368  source = source - 0.5 * ( dlnq * dlnq + dlogn * dlogn ) ;
369 
370  source.std_spectral_base() ;
371  Cmp source_lnq (source) ;
372  Cmp lnq_cmp (logn_quad) ;
373  lnq_cmp.set_etat_qcq() ;
374 
375  mpaff.poisson(source_lnq, par_nul, lnq_cmp) ;
376 
377  lnq = lnq_cmp ;
378 
379  // Metric coefficient psi4 update
380 
381  nn = exp( logn ) ;
383 
384  Scalar qq = exp( lnq ) ;
385  qq.std_spectral_base() ;
386 
387  a_car = qq * qq / ( nn * nn ) ;
388  a_car.std_spectral_base() ;
389 
390  // Relative difference with enthalpy at the previous step
391  // ------------------------------------------------------
392 
393  diff_ent = norme( diffrel(ent, ent_jm1) ) / nzet ;
394 
395  // Next step
396  // ---------
397 
398  ent_jm1 = ent ;
399 
400 
401  } // End of iteration loop
402 
403  //=========================================================================
404  // End of iteration
405  //=========================================================================
406 
407  // The mapping is transfered to that of the star:
408  // ----------------------------------------------
409  mp = mpaff ;
410 
411  // Sets value to all the Tenseur's of the star
412  // -------------------------------------------
413 
414  nn = exp( logn ) ;
416 
417  Scalar qq = exp( lnq ) ;
418  qq.std_spectral_base() ;
419 
420  a_car = qq * qq / ( nn * nn ) ;
421 
422  Sym_tensor gamma_cov(mp, COV, mp.get_bvect_cart()) ;
423  gamma_cov.set_etat_zero() ;
424 
425  for (int i=1; i<=3; i++){
426  gamma_cov.set(i,i) = a_car ;
427  }
428  gamma = gamma_cov ;
429 
430  // ... hydro
431  ent.annule(nzet, nz-1) ; // enthalpy set to zero at the exterior of
432  // the star
433  ener_euler = ener ;
434  s_euler = 3 * press ;
435  gam_euler = 1 ;
436  for(int i=1; i<=3; i++) u_euler.set(i) = 0 ;
437 
438  Star_rot* p_star_rot = dynamic_cast<Star_rot*>(this) ;
439  if ( p_star_rot != 0x0) {
440  p_star_rot->a_car = a_car ;
441  p_star_rot->bbb = qq / nn ;
442  p_star_rot->b_car = p_star_rot->bbb * p_star_rot->bbb ;
443  p_star_rot->dzeta = lnq ;
444  }
445 
446  // Info printing
447  // -------------
448 
449  cout << endl
450  << "Characteristics of the star obtained by Etoile::equilibrium_spher : "
451  << endl
452  << "-----------------------------------------------------------------"
453  << endl ;
454 
455  double ray = mp.val_r(l_b, 1., M_PI/2., 0) ;
456  cout << "Coordinate radius : " << ray / km << " km" << endl ;
457 
458  double rcirc = ray * sqrt(a_car.val_grid_point(l_b, k_b, j_b, i_b) ) ;
459 
460  double compact = qpig/(4.*M_PI) * mass_g() / rcirc ;
461 
462  cout << "Circumferential radius R : " << rcirc/km << " km" << endl ;
463  cout << "Baryon mass : " << mass_b()/msol << " Mo" << endl ;
464  cout << "Gravitational mass M : " << mass_g()/msol << " Mo" << endl ;
465  cout << "Compacity parameter GM/(c^2 R) : " << compact << endl ;
466 
467  //-----------------
468  // Virial theorem
469  //-----------------
470 
471  //... Pressure term
472 
473  source = qpig * a_car * sqrt(a_car) * s_euler ;
474  source.std_spectral_base() ;
475  double vir_mat = source.integrale() ;
476 
477  //... Gravitational term
478 
479  Scalar tmp = lnq - logn ;
480 
481  source = - ( logn.dsdr() * logn.dsdr()
482  - 0.5 * tmp.dsdr() * tmp.dsdr() )
483  * sqrt(a_car) ;
484 
485  source.std_spectral_base() ;
486  double vir_grav = source.integrale() ;
487 
488  //... Relative error on the virial identity GRV3
489 
490  double grv3 = ( vir_mat + vir_grav ) / vir_mat ;
491 
492  cout << "Virial theorem GRV3 : " << endl ;
493  cout << " 3P term : " << vir_mat << endl ;
494  cout << " grav. term : " << vir_grav << endl ;
495  cout << " relative error : " << grv3 << endl ;
496 
497  // To be compatible with class Etoile :
498  //logn = logn - logn_quad ;
499 
500 
501 }
502 }
Scalar dzeta
Metric potential .
Definition: star_rot.h:146
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:607
Scalar a_car
Square of the metric factor A.
Definition: star_rot.h:116
Radial mapping of rather general form.
Definition: map.h:2783
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
const double * get_alpha() const
Returns a pointer on the array alpha (values of in each domain)
Definition: map_et.C:1049
Map & mp
Mapping associated with the star.
Definition: star.h:180
const double * get_beta() const
Returns a pointer on the array beta (values of in each domain)
Definition: map_et.C:1053
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Metric gamma
3-metric
Definition: star.h:235
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:667
Scalar bbb
Metric factor B.
Definition: star_rot.h:119
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:64
virtual void equilibrium_spher(double ent_c, double precis=1.e-14, const Tbl *pent_limit=0x0)
Computes a spherical static configuration.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Scalar ent
Log-enthalpy.
Definition: star.h:190
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Class for isolated rotating stars.
Definition: star_rot.h:98
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:771
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
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.
Scalar ener
Total energy density in the fluid frame.
Definition: star.h:193
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
virtual void poisson(const Cmp &source, Param &par, Cmp &uu) const
Computes the solution of a scalar Poisson equation.
Scalar press
Fluid pressure.
Definition: star.h:194
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:611
Parameter storage.
Definition: param.h:125
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:525
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:760
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Scalar b_car
Square of the metric factor B.
Definition: star_rot.h:122
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
virtual void adapt(const Cmp &ent, const Param &par, int nbr_filtre=0)
Adaptation of the mapping to a given scalar field.
Definition: map_et_adapt.C:111
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar nn
Lapse function N .
Definition: star.h:225
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
virtual double mass_b() const =0
Baryon mass.
Definition: star_global.C:310
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:465
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:282
virtual double mass_g() const =0
Gravitational mass.
Definition: star_global.C:330
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388