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
30char 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
94namespace Lorene {
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
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) ;
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
183 logn_mat_cmp.set_etat_qcq() ;
184
186
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
201 logn_quad_cmp.set_etat_qcq() ;
202
204
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
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 ;
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
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
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() ;
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
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
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
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
Radial mapping of rather general form.
Definition map.h:2752
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.
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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 .
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
virtual double mass_g() const =0
Gravitational mass.
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.
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
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.