LORENE
etoile_equil_spher.C
1 /*
2  * Method of class Etoile to compute a static spherical configuration.
3  *
4  * (see file etoile.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2000-2001 Eric Gourgoulhon
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: etoile_equil_spher.C,v 1.7 2016/12/05 16:17:54 j_novak Exp $
34  * $Log: etoile_equil_spher.C,v $
35  * Revision 1.7 2016/12/05 16:17:54 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.6 2014/10/13 08:52:59 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.5 2008/11/14 13:48:06 e_gourgoulhon
42  * Added parameter pent_limit to force the enthalpy values at the
43  * boundaries between the domains, in case of more than one domain inside
44  * the star.
45  *
46  * Revision 1.4 2004/05/07 12:13:15 k_taniguchi
47  * Change the position of the initialization of alpha_r.
48  *
49  * Revision 1.3 2004/03/25 10:29:06 j_novak
50  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
51  *
52  * Revision 1.2 2003/04/23 15:09:38 j_novak
53  * Standard basis is set to a_car and nnn before exiting.
54  *
55  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
56  * LORENE
57  *
58  * Revision 2.4 2000/02/02 09:23:51 eric
59  * Ajout du theoreme du viriel.
60  * Affichage quantites globales a la fin.
61  *
62  * Revision 2.3 2000/01/28 17:18:36 eric
63  * Modifs mineures.
64  *
65  * Revision 2.2 2000/01/27 16:47:16 eric
66  * Premiere version qui tourne !
67  *
68  * Revision 2.1 2000/01/26 13:18:19 eric
69  * *** empty log message ***
70  *
71  * Revision 2.0 2000/01/24 17:13:56 eric
72  * *** empty log message ***
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_equil_spher.C,v 1.7 2016/12/05 16:17:54 j_novak Exp $
76  *
77  */
78 
79 // Headers C
80 #include "math.h"
81 
82 // Headers Lorene
83 #include "etoile.h"
84 #include "param.h"
85 #include "unites.h"
86 #include "nbr_spx.h"
87 #include "graphique.h"
88 
89 namespace Lorene {
90 void Etoile::equilibrium_spher(double ent_c, double precis, const Tbl* pent_limit){
91 
92  // Fundamental constants and units
93  // -------------------------------
94  using namespace Unites ;
95 
96  // Initializations
97  // ---------------
98 
99  const Mg3d* mg = mp.get_mg() ;
100  int nz = mg->get_nzone() ; // total number of domains
101 
102  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
103  int l_b = nzet - 1 ;
104  int i_b = mg->get_nr(l_b) - 1 ;
105  int j_b = mg->get_nt(l_b) - 1 ;
106  int k_b = 0 ;
107 
108  // Value of the enthalpy defining the surface of the star
109  double ent_b = 0 ;
110 
111  // Initialization of the enthalpy field to the constant value ent_c :
112 
113  ent = ent_c ;
114  ent.annule(nzet, nz-1) ;
115 
116  // Corresponding profiles of baryon density, energy density and pressure
117 
118  equation_of_state() ;
119 
120  // Initial metric
121  a_car = 1 ; // this value will remain unchanged in the Newtonian case
122  beta_auto = 0 ; // this value will remain unchanged in the Newtonian case
123 
124 
125  // Auxiliary quantities
126  // --------------------
127 
128  // Affine mapping for solving the Poisson equations
129  Map_af mpaff(mp);
130 
131  Param par_nul ; // Param (null) for Map_af::poisson.
132 
133  Tenseur ent_jm1(mp) ; // Enthalpy at previous step
134  ent_jm1 = 0 ;
135 
136  Tenseur source(mp) ;
137  Tenseur logn(mp) ;
138  Tenseur logn_quad(mp) ;
139  logn = 0 ;
140  logn_quad = 0 ;
141 
142  Cmp dlogn(mp) ;
143  Cmp dbeta(mp) ;
144 
145  double diff_ent = 1 ;
146  int mermax = 200 ; // Max number of iterations
147 
148  double alpha_r = 1 ;
149 
150  //=========================================================================
151  // Start of iteration
152  //=========================================================================
153 
154  for(int mer=0 ; (diff_ent > precis) && (mer<mermax) ; mer++ ) {
155 
156  cout << "-----------------------------------------------" << endl ;
157  cout << "step: " << mer << endl ;
158  cout << "alpha_r: " << alpha_r << endl ;
159  cout << "diff_ent = " << diff_ent << endl ;
160 
161  //-----------------------------------------------------
162  // Resolution of Poisson equation for ln(N)
163  //-----------------------------------------------------
164 
165  // Matter part of ln(N)
166  // --------------------
167  if (relativistic) {
168  source = a_car * (ener + 3*press) ;
169  }
170  else {
171  source = nbar ;
172  }
173 
174  (source.set()).set_dzpuis(4) ;
175 
176  source.set_std_base() ; // Sets the standard spectral bases.
177 
179 
180  mpaff.poisson(source(), par_nul, logn_auto.set()) ;
181 
182  // NB: at this stage logn_auto is in machine units, not in c^2
183 
184  // Quadratic part of ln(N)
185  // -----------------------
186 
187  if (relativistic) {
188 
189  mpaff.dsdr(logn(), dlogn) ;
190  mpaff.dsdr(beta_auto(), dbeta) ;
191 
192  source = - dlogn * dbeta ;
193 
194  logn_quad.set_etat_qcq() ;
195 
196  mpaff.poisson(source(), par_nul, logn_quad.set()) ;
197 
198  }
199 
200  //-----------------------------------------------------
201  // Computation of the new radial scale
202  //-----------------------------------------------------
203 
204  // alpha_r (r = alpha_r r') is determined so that the enthalpy
205  // takes the requested value ent_b at the stellar surface
206 
207  double nu_mat0_b = logn_auto()(l_b, k_b, j_b, i_b) ;
208  double nu_mat0_c = logn_auto()(0, 0, 0, 0) ;
209 
210  double nu_quad0_b = logn_quad()(l_b, k_b, j_b, i_b) ;
211  double nu_quad0_c = logn_quad()(0, 0, 0, 0) ;
212 
213  double alpha_r2 = ( ent_c - ent_b - nu_quad0_b + nu_quad0_c )
214  / ( qpig*(nu_mat0_b - nu_mat0_c) ) ;
215 
216  alpha_r = sqrt(alpha_r2) ;
217 
218 
219  // One domain inside the star:
220  // ---------------------------
221  if(nzet==1) {
222 
223  mpaff.homothetie( alpha_r ) ;
224 
225  }
226 
227  //--------------------
228  // First integral
229  //--------------------
230 
231  // Gravitation potential in units c^2 :
232  logn_auto = alpha_r2*qpig * logn_auto ;
233  logn = logn_auto + logn_quad ;
234 
235  // Enthalpy in all space
236  double logn_c = logn()(0, 0, 0, 0) ;
237  ent = ent_c - logn() + logn_c ;
238 
239  // Two or more domains inside the star:
240  // ------------------------------------
241  if(nzet>1) {
242 
243  // Parameters for the function Map_et::adapt
244  // -----------------------------------------
245 
246  Param par_adapt ;
247  int nitermax = 100 ;
248  int niter ;
249  int adapt_flag = 1 ; // 1 = performs the full computation,
250  // 0 = performs only the rescaling by
251  // the factor alpha_r
252  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
253  // isosurfaces
254 
255  int nzadapt = nzet ;
256 
257  //cout << "no. of domains where the ent adjustment will be done: " << nzet << endl ;
258  //cout << "ent limits: " << ent_limit << endl ;
259 
260  double precis_adapt = 1.e-14 ;
261 
262  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
263 
264  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
265  // locate zeros by the secant method
266  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
267  // to the isosurfaces of ent is to be
268  // performed
269  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
270  // the enthalpy isosurface
271  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
272  // 0 = performs only the rescaling by
273  // the factor alpha_r
274  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
275  // (theta_*, phi_*)
276  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
277  // (theta_*, phi_*)
278 
279  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
280  // the secant method
281 
282  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
283  // the determination of zeros by
284  // the secant method
285  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
286 
287  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
288  // distances will be multiplied
289 
290  // Enthalpy values for the adaptation
291  Tbl ent_limit(nzet) ;
292  if (pent_limit != 0x0) ent_limit = *pent_limit ;
293 
294  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
295  // to define the isosurfaces.
296 
297  double* bornes = new double[nz+1] ;
298  bornes[0] = 0. ;
299 
300  for(int l=0; l<nz; l++) {
301 
302  bornes[l+1] = mpaff.get_alpha()[l] + mpaff.get_beta()[l] ;
303 
304  }
305  bornes[nz] = __infinity ;
306 
307  Map_et mp0(*mg, bornes) ;
308 
309  mp0 = mpaff;
310  mp0.adapt(ent(), par_adapt) ;
311 
312  //Map_af mpaff_prev (mpaff) ;
313 
314  double alphal, betal ;
315 
316  for(int l=0; l<nz; l++) {
317 
318  alphal = mp0.get_alpha()[l] ;
319  betal = mp0.get_beta()[l] ;
320 
321  mpaff.set_alpha(alphal, l) ;
322  mpaff.set_beta(betal, l) ;
323 
324  }
325 
326 
327  //mbtest
328  int num_r1 = mg->get_nr(0) - 1;
329 
330  cout << "Pressure difference:" << get_press()()(0,0,0,num_r1) - get_press()()(1,0,0,0) << endl ;
331  cout << "Difference in enthalpies at the domain boundary:" << endl ;
332  cout << get_ent()()(0,0,0,num_r1) << endl ;
333  cout << get_ent()()(1,0,0,0) << endl ;
334 
335  cout << "Enthalpy difference: " << get_ent()()(0,0,0,num_r1) - get_ent()()(1,0,0,0) << endl ;
336 
337  // Computation of the enthalpy at the new grid points
338  //----------------------------------------------------
339 
340  //mpaff.reevaluate(&mpaff_prev, nzet+1, ent.set()) ;
341 
342  }
343 
344  //---------------------
345  // Equation of state
346  //---------------------
347 
348  equation_of_state() ;
349 
350  if (relativistic) {
351 
352  //----------------------------
353  // Equation for beta = ln(AN)
354  //----------------------------
355 
356  mpaff.dsdr(logn(), dlogn) ;
357  mpaff.dsdr(beta_auto(), dbeta) ;
358 
359  source = 3 * qpig * a_car * press ;
360 
361  source = source()
362  - 0.5 * ( dlogn * dlogn + dbeta * dbeta ) ;
363 
364  source.set_std_base() ; // Sets the standard spectral bases.
365 
367 
368  mpaff.poisson(source(), par_nul, beta_auto.set()) ;
369 
370 
371  // Metric coefficient A^2 update
372 
373  a_car = exp(2*(beta_auto - logn)) ;
374 
375  }
376 
377  // Relative difference with enthalpy at the previous step
378  // ------------------------------------------------------
379 
380  diff_ent = norme( diffrel(ent(), ent_jm1()) ) / nzet ;
381 
382  // Next step
383  // ---------
384 
385  ent_jm1 = ent ;
386 
387 
388  } // End of iteration loop
389 
390  //=========================================================================
391  // End of iteration
392  //=========================================================================
393 
394 
395  // The mapping is transfered to that of the star:
396  // ----------------------------------------------
397  mp = mpaff ;
398 
399 
400  // Sets value to all the Tenseur's of the star
401  // -------------------------------------------
402 
403  // ... hydro
404  ent.annule(nzet, nz-1) ; // enthalpy set to zero at the exterior of
405  // the star
406  ener_euler = ener ;
407  s_euler = 3 * press ;
408  gam_euler = 1 ;
409  u_euler = 0 ;
410 
411  // ... metric
412  nnn = exp( unsurc2 * logn ) ;
413  nnn.set_std_base() ;
414  shift = 0 ;
415  a_car.set_std_base() ;
416 
417  // Info printing
418  // -------------
419 
420  cout << endl
421  << "Characteristics of the star obtained by Etoile::equilibrium_spher : "
422  << endl
423  << "-----------------------------------------------------------------"
424  << endl ;
425 
426  double ray = mp.val_r(l_b, 1., M_PI/2., 0) ;
427  cout << "Coordinate radius : " << ray / km << " km" << endl ;
428 
429  double rcirc = ray * sqrt( a_car()(l_b, k_b, j_b, i_b) ) ;
430 
431  double compact = qpig/(4.*M_PI) * mass_g() / rcirc ;
432 
433  cout << "Circumferential radius R : " << rcirc/km << " km" << endl ;
434  cout << "Baryon mass : " << mass_b()/msol << " Mo" << endl ;
435  cout << "Gravitational mass M : " << mass_g()/msol << " Mo" << endl ;
436  cout << "Compacity parameter GM/(c^2 R) : " << compact << endl ;
437 
438 
439  //-----------------
440  // Virial theorem
441  //-----------------
442 
443  //... Pressure term
444 
445  source = qpig * a_car * sqrt(a_car) * s_euler ;
446  source.set_std_base() ;
447  double vir_mat = source().integrale() ;
448 
449  //... Gravitational term
450 
451  Cmp tmp = beta_auto() - logn() ;
452 
453  source = - ( logn().dsdr() * logn().dsdr()
454  - 0.5 * tmp.dsdr() * tmp.dsdr() )
455  * sqrt(a_car()) ;
456 
457  source.set_std_base() ;
458  double vir_grav = source().integrale() ;
459 
460  //... Relative error on the virial identity GRV3
461 
462  double grv3 = ( vir_mat + vir_grav ) / vir_mat ;
463 
464  cout << "Virial theorem GRV3 : " << endl ;
465  cout << " 3P term : " << vir_mat << endl ;
466  cout << " grav. term : " << vir_grav << endl ;
467  cout << " relative error : " << grv3 << endl ;
468 
469  if (nzet > 1) {
470  cout.precision(10) ;
471 
472  for (int ltrans = 0; ltrans < nzet-1; ltrans++) {
473  cout << endl << "Values at boundary between domains no. " << ltrans << " and " << ltrans+1 << " for theta = pi/2 and phi = 0 :" << endl ;
474 
475  double rt1 = mp.val_r(ltrans, 1., M_PI/2, 0.) ;
476  double rt2 = mp.val_r(ltrans+1, -1., M_PI/2, 0.) ;
477  cout << " Coord. r [km] (left, right, rel. diff) : "
478  << rt1 / km << " " << rt2 / km << " " << (rt2 - rt1)/rt1 << endl ;
479 
480  int ntm1 = mg->get_nt(ltrans) - 1;
481  int nrm1 = mg->get_nr(ltrans) - 1 ;
482  double ent1 = ent()(ltrans, 0, ntm1, nrm1) ;
483  double ent2 = ent()(ltrans+1, 0, ntm1, 0) ;
484  cout << " Enthalpy (left, right, rel. diff) : "
485  << ent1 << " " << ent2 << " " << (ent2-ent1)/ent1 << endl ;
486 
487  double press1 = press()(ltrans, 0, ntm1, nrm1) ;
488  double press2 = press()(ltrans+1, 0, ntm1, 0) ;
489  cout << " Pressure (left, right, rel. diff) : "
490  << press1 << " " << press2 << " " << (press2-press1)/press1 << endl ;
491 
492  double nb1 = nbar()(ltrans, 0, ntm1, nrm1) ;
493  double nb2 = nbar()(ltrans+1, 0, ntm1, 0) ;
494  cout << " Baryon density (left, right, rel. diff) : "
495  << nb1 << " " << nb2 << " " << (nb2-nb1)/nb1 << endl ;
496  }
497  }
498 
499 
500 /* double r_max = 1.2 * ray_eq() ;
501  des_profile(nbar(), 0., r_max, M_PI/2, 0., "n", "Baryon density") ;
502  des_profile(ener(), 0., r_max, M_PI/2, 0., "e", "Energy density") ;
503  des_profile(press(), 0., r_max, M_PI/2, 0., "p", "Pressure") ;
504  des_profile(ent(), 0., r_max, M_PI/2, 0., "H", "Enthalpy") ;
505 */
506 
507 }
508 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:87
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:916
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
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
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 homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:667
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
virtual double mass_g() const
Gravitational mass.
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
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
virtual void equilibrium_spher(double ent_c, double precis=1.e-14, const Tbl *ent_limit=0x0)
Computes a spherical static configuration.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:471
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: etoile.h:445
Tenseur press
Fluid pressure.
Definition: etoile.h:464
Tenseur shift
Total shift vector.
Definition: etoile.h:515
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
const Tenseur & get_press() const
Returns the fluid pressure.
Definition: etoile.h:685
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:771
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.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:477
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.
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:462
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:474
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
virtual double mass_b() const
Baryon mass.
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:569
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:760
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:435
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Tenseur a_car
Total conformal factor .
Definition: etoile.h:518
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:440
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
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:463
Affine radial mapping.
Definition: map.h:2048
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition: etoile.h:487
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:460
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
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
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition: etoile.h:509
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:468
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:282
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
const Tenseur & get_ent() const
Returns the enthalpy field.
Definition: etoile.h:676