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
30char 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
86namespace Lorene {
87void 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
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
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 ;
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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
virtual void equilibrium_spher(double ent_c, double precis=1.e-14, const Tbl *ent_limit=0x0)
Computes a spherical static configuration.
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:484
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:471
Map & mp
Mapping associated with the star.
Definition etoile.h:429
const Tenseur & get_ent() const
Returns the enthalpy field.
Definition etoile.h:673
virtual double mass_b() const
Baryon mass.
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Tenseur press
Fluid pressure.
Definition etoile.h:461
virtual double mass_g() const
Gravitational mass.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Tenseur shift
Total shift vector.
Definition etoile.h:512
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:506
const Tenseur & get_press() const
Returns the fluid pressure.
Definition etoile.h:682
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition etoile.h:442
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.
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.
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:900
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
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
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.