LORENE
star.C
1/*
2 * Methods of class Star
3 *
4 * (see file star.h for documentation)
5 */
6
7/*
8 * Copyright (c) 2004 Francois Limousin
9 *
10 * Copyright (c) 2000-2001 Eric Gourgoulhon (for preceding class Etoile)
11 * Copyright (c) 2000-2001 Keisuke Taniguchi (for preceding class Etoile)
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or
18 * (at your option) any later version.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31
32char star_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star.C,v 1.21 2014/10/13 08:53:37 j_novak Exp $" ;
33
34/*
35 * $Id: star.C,v 1.21 2014/10/13 08:53:37 j_novak Exp $
36 * $Log: star.C,v $
37 * Revision 1.21 2014/10/13 08:53:37 j_novak
38 * Lorene classes and functions now belong to the namespace Lorene.
39 *
40 * Revision 1.20 2013/04/20 20:56:15 m_bejger
41 * Fix for three domains in star in Star::equation_of_state from Etoile/etoile.C
42 *
43 * Revision 1.19 2010/02/02 12:45:16 e_gourgoulhon
44 * Improved the display (operator>>)
45 *
46 * Revision 1.18 2010/01/26 16:49:03 e_gourgoulhon
47 * Commented the test on the relativistic character of the EOS: the
48 * relativity parameter is not defined (yet !) in the base class Star.
49 *
50 * Revision 1.17 2007/11/06 16:22:03 j_novak
51 * The data member stress_euler is now a Sym_tensor instead of a Tensor.
52 *
53 * Revision 1.16 2007/06/21 19:53:47 k_taniguchi
54 * Addition of p_ray_eq_3pis2
55 *
56 * Revision 1.15 2005/09/14 12:30:52 f_limousin
57 * Saving of fields lnq and logn in class Star.
58 *
59 * Revision 1.14 2005/09/13 19:38:31 f_limousin
60 * Reintroduction of the resolution of the equations in cartesian coordinates.
61 *
62 * Revision 1.13 2005/02/17 17:29:04 f_limousin
63 * Change the name of some quantities to be consistent with other classes
64 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
65 *
66 * Revision 1.12 2005/01/05 17:43:03 f_limousin
67 * u_euler is now constructed in the spherical triad to be consistent
68 * with all the others vectors ans tensors.
69 *
70 * Revision 1.11 2004/11/11 16:29:49 j_novak
71 * set_der_0x0 is no longer virtual (to be coherent with Tensor/Scalar classes).
72 *
73 * Revision 1.10 2004/06/22 12:48:08 f_limousin
74 * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
75 *
76 * Revision 1.9 2004/06/07 16:21:35 f_limousin
77 * Add outputs
78 *
79 * Revision 1.8 2004/04/08 16:32:10 f_limousin
80 * The new variable is ln(Q) instead of Q=psi^2*N. It improves the
81 * convergence of the code.
82 *
83 * Revision 1.7 2004/03/25 10:29:26 j_novak
84 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
85 *
86 * Revision 1.6 2004/03/08 11:48:00 f_limousin
87 * Error in del_deriv() and set_der_0x0() : p_mass_b and p_mass_g were
88 * missing. And so they were never recomputed.
89 *
90 * Revision 1.5 2004/02/27 09:36:46 f_limousin
91 * u_euler is now constructed on a cartesian basis instead
92 * of a spherical one.
93 *
94 * Revision 1.4 2004/02/21 17:05:13 e_gourgoulhon
95 * Method Scalar::point renamed Scalar::val_grid_point.
96 * Method Scalar::set_point renamed Scalar::set_grid_point.
97 *
98 * Revision 1.3 2004/01/20 15:16:58 f_limousin
99 * First version
100 *
101 *
102 * $Header: /cvsroot/Lorene/C++/Source/Star/star.C,v 1.21 2014/10/13 08:53:37 j_novak Exp $
103 *
104 */
105
106// Headers C
107#include "math.h"
108
109// Headers Lorene
110#include "star.h"
111#include "eos.h"
112#include "utilitaires.h"
113#include "param.h"
114#include "unites.h"
115
116
117 //--------------//
118 // Constructors //
119 //--------------//
120
121// Standard constructor
122// --------------------
123namespace Lorene {
125 : mp(mpi),
126 nzet(nzet_i),
127 eos(eos_i),
128 ent(mpi),
129 nbar(mpi),
130 ener(mpi),
131 press(mpi),
132 ener_euler(mpi),
133 s_euler(mpi),
134 gam_euler(mpi),
135 u_euler(mpi, CON, mp.get_bvect_spher()),
136 stress_euler(mpi, CON, mp.get_bvect_spher()),
137 logn(mpi),
138 nn(mpi),
139 beta(mpi, CON, mp.get_bvect_spher()),
140 lnq(mpi),
141 gamma(mp.flat_met_spher()){
142
143
144 // Check of the EOS
145// const Eos_poly_newt* p_eos_poly_newt =
146// dynamic_cast<const Eos_poly_newt*>( &eos ) ;
147//
148// const Eos_incomp_newt* p_eos_incomp_newt =
149// dynamic_cast<const Eos_incomp_newt*>( &eos ) ;
150//
151//
152//
153// if (p_eos_poly_newt != 0x0) {
154// cout <<
155// "Star::Star : the EOS Eos_poly_newt must not be employed"
156// << " for a relativistic star ! " << endl ;
157// cout << "(Use Eos_poly instead)" << endl ;
158// abort() ;
159// }
160// if (p_eos_incomp_newt != 0x0) {
161// cout <<
162// "Star::Star : the EOS Eos_incomp_newt must not be employed"
163// << " for a relativistic star ! " << endl ;
164// cout << "(Use Eos_incomp instead)" << endl ;
165// abort() ;
166// }
167//
168 // Pointers of derived quantities initialized to zero :
169 set_der_0x0() ;
170
171 // All the matter quantities are initialized to zero :
172 nbar = 0 ;
173 ener = 0 ;
174 press = 0 ;
175 ent = 0 ;
176 ener_euler = 0 ;
177 s_euler = 0 ;
178 gam_euler = 1. ;
182
183 // The metric is initialized to the flat one :
184 Metric flat(mp.flat_met_spher()) ;
185 flat.cov() ;
186 gamma = flat ;
187
188 logn = 0 ;
189 nn = 1. ;
192 lnq = 0 ;
193
194}
195
196// Copy constructor
197// ----------------
198Star::Star(const Star& et)
199 : mp(et.mp),
200 nzet(et.nzet),
201 eos(et.eos),
202 ent(et.ent),
203 nbar(et.nbar),
204 ener(et.ener),
205 press(et.press),
206 ener_euler(et.ener_euler),
207 s_euler(et.s_euler),
208 gam_euler(et.gam_euler),
209 u_euler(et.u_euler),
210 stress_euler(et.stress_euler),
211 logn(et.logn),
212 nn(et.nn),
213 beta(et.beta),
214 lnq(et.lnq),
215 gamma(et.gamma){
216
217 set_der_0x0() ;
218
219}
220
221// Constructor from a file
222// -----------------------
224 : mp(mpi),
225 eos(eos_i),
226 ent(mpi),
227 nbar(mpi),
228 ener(mpi),
229 press(mpi),
230 ener_euler(mpi),
231 s_euler(mpi),
232 gam_euler(mpi),
233 u_euler(mpi, CON, mp.get_bvect_spher()),
234 stress_euler(mpi, CON, mp.get_bvect_spher()),
235 logn(mpi, *(mpi.get_mg()), fich),
236 nn(mpi),
237 beta(mpi, CON, mp.get_bvect_spher()),
238 lnq(mpi, *(mpi.get_mg()), fich),
239 gamma(mpi.flat_met_spher()){
240
241 // Star parameters
242 // -----------------
243
244 // nzet is read in the file:
245 int xx ;
246 fread_be(&xx, sizeof(int), 1, fich) ;
247 nzet = xx ;
248
249 // Equation of state
250 // -----------------
251
252 // Read of the saved EOS
254
255 // Comparison with the assigned EOS:
256 if (eos != *p_eos_file) {
257 cout <<
258 "Star::Star(const Map&, const Eos&, FILE*) : the EOS given in "
259 << endl <<
260 " argument and that read in the file are different !" << endl ;
261 abort() ;
262 }
263
264 // p_eos_file is no longer required (it was used only for checking the
265 // EOS compatibility)
266 delete p_eos_file ;
267
268 // Read of the saved fields:
269 // ------------------------
270 Scalar ent_file(mp, *(mp.get_mg()), fich) ;
271 ent = ent_file ;
274 nn = 1 ;
276
277 // Pointers of derived quantities initialized to zero
278 // --------------------------------------------------
279 set_der_0x0() ;
280
281}
282
283 //------------//
284 // Destructor //
285 //------------//
286
288
289 del_deriv() ;
290
291}
292
293
294 //----------------------------------//
295 // Management of derived quantities //
296 //----------------------------------//
297
298void Star::del_deriv() const {
299
300 if (p_mass_b != 0x0) delete p_mass_b ;
301 if (p_mass_g != 0x0) delete p_mass_g ;
302 if (p_ray_eq != 0x0) delete p_ray_eq ;
303 if (p_ray_eq_pis2 != 0x0) delete p_ray_eq_pis2 ;
304 if (p_ray_eq_pi != 0x0) delete p_ray_eq_pi ;
305 if (p_ray_eq_3pis2 != 0x0) delete p_ray_eq_3pis2 ;
306 if (p_ray_pole != 0x0) delete p_ray_pole ;
307 if (p_l_surf != 0x0) delete p_l_surf ;
308 if (p_xi_surf != 0x0) delete p_xi_surf ;
309
311}
312
313
314
315
316void Star::set_der_0x0() const {
317
318 p_mass_b = 0x0 ;
319 p_mass_g = 0x0 ;
320 p_ray_eq = 0x0 ;
321 p_ray_eq_pis2 = 0x0 ;
322 p_ray_eq_pi = 0x0 ;
323 p_ray_eq_3pis2 = 0x0 ;
324 p_ray_pole = 0x0 ;
325 p_l_surf = 0x0 ;
326 p_xi_surf = 0x0 ;
327
328}
329
341
342
343
344
345 //--------------//
346 // Assignment //
347 //--------------//
348
349// Assignment to another Star
350// ----------------------------
351void Star::operator=(const Star& et) {
352
353 assert( &(et.mp) == &mp ) ; // Same mapping
354 assert( &(et.eos) == &eos ) ; // Same EOS
355
356 nzet = et.nzet ;
357 ent = et.ent ;
358 nbar = et.nbar ;
359 ener = et.ener ;
360 press = et.press ;
362 s_euler = et.s_euler ;
363 gam_euler = et.gam_euler ;
364 u_euler = et.u_euler ;
366 logn = et.logn ;
367 nn = et.nn ;
368 beta = et.beta ;
369 lnq = et.lnq ;
370 gamma = et.gamma ;
371
372 del_deriv() ; // Deletes all derived quantities
373
374}
375
376// Assignment of the enthalpy field
377// --------------------------------
378
380
381 ent = ent_i ;
382
383 // Update of (nbar, ener, press) :
385
386 // The derived quantities are obsolete:
387 del_deriv() ;
388
389}
390
391 //--------------//
392 // Outputs //
393 //--------------//
394
395// Save in a file
396// --------------
397void Star::sauve(FILE* fich) const {
398
399 logn.sauve(fich) ;
400 lnq.sauve(fich) ;
401
402 int xx = nzet ;
403 fwrite_be(&xx, sizeof(int), 1, fich) ;
404
405 eos.sauve(fich) ;
406 ent.sauve(fich) ;
407}
408
409// Printing
410// --------
411
413 et >> ost ;
414 return ost ;
415}
416
418
419 using namespace Unites ;
420
421 ost << endl ;
422
423 ost << "Number of domains occupied by the star : " << nzet << endl ;
424
425 ost << "Equation of state : " << endl ;
426 ost << eos << endl ;
427
428 ost << endl << "Central enthalpy : " << ent.val_grid_point(0,0,0,0) << " c^2" << endl ;
429 ost << "Central proper baryon density : " << nbar.val_grid_point(0,0,0,0)
430 << " x 0.1 fm^-3" << endl ;
431 ost << "Central proper energy density : " << ener.val_grid_point(0,0,0,0)
432 << " rho_nuc c^2" << endl ;
433 ost << "Central pressure : " << press.val_grid_point(0,0,0,0)
434 << " rho_nuc c^2" << endl ;
435
436 ost << endl ;
437 ost << "Central lapse N : " << nn.val_grid_point(0,0,0,0) << endl ;
438// ost << "Central value of lnq : " << lnq.val_grid_point(0,0,0,0) << endl ;
439
440 ost << endl
441 << "Coordinate equatorial radius (phi=0) a1 = "
442 << ray_eq()/km << " km" << endl ;
443 ost << "Coordinate equatorial radius (phi=pi/2) a2 = "
444 << ray_eq_pis2()/km << " km" << endl ;
445 ost << "Coordinate equatorial radius (phi=pi): "
446 << ray_eq_pi()/km << " km" << endl ;
447 ost << "Coordinate polar radius a3 = "
448 << ray_pole()/km << " km" << endl ;
449 ost << "Axis ratio a2/a1 = " << ray_eq_pis2() / ray_eq()
450 << " a3/a1 = " << ray_pole() / ray_eq() << endl ;
451 ost << endl << "Baryon mass : " << mass_b() / msol << " M_sol" << endl ;
452 ost << "Gravitational mass : " << mass_g() / msol << " M_sol" << endl ;
453
454
455 return ost ;
456}
457
458 //-----------------------------------------//
459 // Computation of hydro quantities //
460 //-----------------------------------------//
461
463
464 Scalar ent_eos = ent ;
465
466
467 // Slight rescale of the enthalpy field in case of 2 domains inside the
468 // star
469
470
471 double epsilon = 1.e-12 ;
472
473 const Mg3d* mg = mp.get_mg() ;
474 int nz = mg->get_nzone() ;
475
476 Mtbl xi(mg) ;
477 xi.set_etat_qcq() ;
478 for (int l=0; l<nz; l++) {
479 xi.t[l]->set_etat_qcq() ;
480 for (int k=0; k<mg->get_np(l); k++) {
481 for (int j=0; j<mg->get_nt(l); j++) {
482 for (int i=0; i<mg->get_nr(l); i++) {
483 xi.set(l,k,j,i) =
484 mg->get_grille3d(l)->x[i] ;
485 }
486 }
487 }
488
489 }
490
492 fact_ent.allocate_all() ;
493
494 fact_ent.set_domain(0) = 1 + epsilon * xi(0) * xi(0) ;
495 fact_ent.set_domain(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
496
497 for (int l=nzet; l<nz; l++) {
498 fact_ent.set_domain(l) = 1 ;
499 }
500
501 if (nzet > 1) {
502
503 if(nzet == 3) {
504 fact_ent.set_domain(1) = 1 - 0.5 * epsilon * (xi(1) - 0.5) * (xi(1) - 0.5) ;
505 fact_ent.set_domain(2) = 1 - 0.25 * epsilon * (xi(2) - 1) * (xi(2) - 1) ;
506 }
507
508 if (nzet > 3) {
509
510 cout << "Star::equation_of_state: not ready yet for nzet > 3 !"
511 << endl ;
512 }
513
515 ent_eos.std_spectral_base() ;
516 }
517
518
519
520
521
522 // Call to the EOS (the EOS is called domain by domain in order to
523 // allow for the use of MEos)
524
525 Scalar tempo(mp) ;
526
528 nbar = 0 ;
529 for (int l=0; l<nzet; l++) {
530
531 Param par ; // Paramater for multi-domain equation of state
532 par.add_int(l) ;
533
534 tempo = eos.nbar_ent(ent_eos, 1, l, &par) ;
535
536 nbar = nbar + tempo ;
537
538 }
539
541 ener = 0 ;
542 for (int l=0; l<nzet; l++) {
543
544 Param par ; // Paramater for multi-domain equation of state
545 par.add_int(l) ;
546
547 tempo = eos.ener_ent(ent_eos, 1, l, &par) ;
548
549 ener = ener + tempo ;
550
551 }
552
554 press = 0 ;
555 for (int l=0; l<nzet; l++) {
556
557 Param par ; // Paramater for multi-domain equation of state
558 par.add_int(l) ;
559
560 tempo = eos.press_ent(ent_eos, 1, l, &par) ;
561
562 press = press + tempo ;
563
564 }
565
566
567 // Set the bases for spectral expansion
571
572 // The derived quantities are obsolete
573 del_deriv() ;
574
575}
576
578
579 cout <<
580 "Star::hydro_euler : hydro_euler must be called via a derived class"
581 << endl << " of Star !" << endl ;
582
583 abort() ;
584
585}
586}
Equation of state base class.
Definition eos.h:190
static Eos * eos_from_file(FILE *)
Construction of an EOS from a binary file.
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Multi-domain grid.
Definition grilles.h:273
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
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
Multi-domain array.
Definition mtbl.h:118
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition mtbl.h:221
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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
virtual void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition scalar.C:344
Base class for stars.
Definition star.h:175
virtual ~Star()
Destructor.
Definition star.C:287
Star(Map &mp_i, int nzet_i, const Eos &eos_i)
Standard constructor.
Definition star.C:124
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 del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition star.C:330
double * p_mass_b
Baryon mass.
Definition star.h:268
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition star.C:577
double * p_ray_eq
Coordinate radius at , .
Definition star.h:242
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
double * p_ray_eq_pi
Coordinate radius at , .
Definition star.h:248
Itbl * p_l_surf
Description of the stellar surface: 2-D Itbl containing the values of the domain index l on the surfa...
Definition star.h:260
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition star.C:462
double * p_ray_eq_pis2
Coordinate radius at , .
Definition star.h:245
double * p_ray_eq_3pis2
Coordinate radius at , .
Definition star.h:251
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star.C:298
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition star.C:417
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition star.C:316
double ray_eq() const
Coordinate radius at , [r_unit].
void set_enthalpy(const Scalar &)
Assignment of the enthalpy field.
Definition star.C:379
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
double * p_mass_g
Gravitational mass.
Definition star.h:269
Metric gamma
3-metric
Definition star.h:235
Tbl * p_xi_surf
Description of the stellar surface: 2-D Tbl containing the values of the radial coordinate on the su...
Definition star.h:266
virtual void sauve(FILE *) const
Save in a file.
Definition star.C:397
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition star.h:212
Scalar press
Fluid pressure.
Definition star.h:194
double ray_eq_pi() const
Coordinate radius at , [r_unit].
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
void operator=(const Star &)
Assignment to another Star.
Definition star.C:351
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
double * p_ray_pole
Coordinate radius at .
Definition star.h:254
Vector beta
Shift vector.
Definition star.h:228
double ray_pole() const
Coordinate radius at [r_unit].
virtual double mass_b() const =0
Baryon mass.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
virtual void set_etat_nondef()
Sets the logical state of all components to ETATNONDEF (undefined state).
Definition tensor.C:489
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.