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