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