LORENE
star_bhns_spher.C
1/*
2 * Method of class Star_bhns to compute a spherical star configuration
3 *
4 * (see file star_bhns.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005,2007 Keisuke Taniguchi
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 version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char star_bhns_spher_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_spher.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $" ;
29
30/*
31 * $Id: star_bhns_spher.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
32 * $Log: star_bhns_spher.C,v $
33 * Revision 1.4 2014/10/13 08:53:41 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2014/10/06 15:13:16 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.2 2008/05/15 19:16:54 k_taniguchi
40 * Change of some parameters.
41 *
42 * Revision 1.1 2007/06/22 01:32:19 k_taniguchi
43 * *** empty log message ***
44 *
45 *
46 * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_spher.C,v 1.4 2014/10/13 08:53:41 j_novak Exp $
47 *
48 */
49
50// C++ headers
51//#include <>
52
53// C headers
54#include <cmath>
55
56// Lorene headers
57#include "star_bhns.h"
58#include "param.h"
59#include "cmp.h"
60#include "tenseur.h"
61#include "unites.h"
62
63namespace Lorene {
64void Star_bhns::equil_spher_bhns(double ent_c, double precis) {
65
66 // Fundamental constants and units
67 // -------------------------------
68 using namespace Unites ;
69
70 // Initializations
71 // ---------------
72
73 const Mg3d* mg = mp.get_mg() ;
74 int nz = mg->get_nzone() ;
75
76 // Index of the point at phi=0, theta=pi/2 at the surface of the star:
77 int l_b = nzet - 1 ;
78 int i_b = mg->get_nr(l_b) - 1 ;
79 int j_b = mg->get_nt(l_b) - 1 ;
80 int k_b = 0 ;
81
82 // Value of the enthalpy defining the surface of the star
83 double ent_b = 0. ;
84
85 // Initialization of the enthalpy field to the constant value ent_c :
86 ent = ent_c ;
87 ent.annule(nzet, nz-1) ;
88
89 // Corresponding profiles of baryon density, energy density and pressure
91
92 // Initial metric
93 lapconf_auto = 1. ;
95 confo_auto = 1. ;
97
98 // Auxiliary quantities
99 // --------------------
100
101 // Affine mapping for solving the Poisson equations
102 Map_af mpaff(mp) ;
103
104 Param par_nul ; // Param (null) for Map_af::poisson
105
106 Scalar ent_jm1(mp) ; // Enthalpy at previous step
107 ent_jm1 = ent_c ;
108 ent_jm1.std_spectral_base() ;
109 ent_jm1.annule(nzet, nz-1) ;
110
115 lapconf_auto_m1 = 0. ;
116 confo_auto_m1 = 0. ;
117
118 double diff_ent = 1. ;
119 int mermax = 200 ; // Maximum number of iterations
120
121 double alpha_r = 1. ;
122
123 //==========================================================//
124 // Start of iteration //
125 //==========================================================//
126
127 for (int mer=0 ; (diff_ent > precis) && (mer<mermax) ; mer++ ) {
128
129 cout << "-----------------------------------------------" << endl ;
130 cout << "step: " << mer << endl ;
131 cout << "alpha_r: " << alpha_r << endl ;
132 cout << "diff_ent = " << diff_ent << endl ;
133
134 //----------------------------------------------------
135 // Resolution of Poisson equation for lapconf function
136 //----------------------------------------------------
137
138 // Matter part of lapconf
139 // ----------------------
141 * (0.5*ener + 3.*press) ;
142
143 source_lapconf.inc_dzpuis(4-source_lapconf.get_dzpuis()) ;
144 source_lapconf.std_spectral_base() ;
145
148 lapconf_cmp.set_etat_qcq() ;
149
151
152 // Re-construction of a scalar
154
155 //-------------------------------------
156 // Computation of the new radial scale
157 //-------------------------------------
158
159 double exp_ent_c = exp(ent_c) ;
160 double exp_ent_b = exp(ent_b) ;
161
162 double lap_auto_c = lapconf_auto_m1.val_grid_point(0,0,0,0) ;
163 double lap_auto_b = lapconf_auto_m1.val_grid_point(l_b,k_b,j_b,i_b) ;
164
165 double confo_c = confo_auto.val_grid_point(0,0,0,0) ;
167
171
173
174 // New radial scale
175 mpaff.homothetie( alpha_r ) ;
176
177 //----------------
178 // First integral
179 //----------------
180
181 // Lapconf function
184
185 // Enthalpy in all space
186 double lapconfo_c = lapconf_auto.val_grid_point(0,0,0,0) ;
191
192 //-------------------
193 // Equation of state
194 //-------------------
195
197
198 //-----------------------------------------------------
199 // Resolution of Poisson equation for conformal factor
200 //-----------------------------------------------------
201
202 source_confo = - 0.5 * qpig * pow(confo_auto,5.) * ener ;
203 source_confo.inc_dzpuis(4-source_confo.get_dzpuis()) ;
204 source_confo.std_spectral_base() ;
205
208 cnf_auto_cmp.set_etat_qcq() ;
209
211
212 // Re-construction of a scalr
214
216
217 // Relative difference with enthalphy at the previous step
218 // -------------------------------------------------------
219
221
222 // Next step
223 // ---------
224
225 ent_jm1 = ent ;
226
227 } // End of iteration loop
228
229 //========================================================//
230 // End of iteration //
231 //========================================================//
232
233 // The mapping is transfered to that of the star
234 // ---------------------------------------------
235 mp = mpaff ;
236
237 // Sets values
238 // -----------
239
240 // ... hydro
241 ent.annule(nzet, nz-1) ;
242
243 ener_euler = ener ;
244 s_euler = 3. * press ;
245 gam_euler = 1. ;
246 for(int i=1; i<=3; i++)
247 u_euler.set(i) = 0 ;
248
249 // ... metric
254 psi4 = pow(confo_auto, 4.) ;
255 for (int i=1; i<=3; i++)
256 shift_auto.set(i) = 0. ;
257
258 // Info printing
259 // -------------
260
261 cout << endl
262 << "Characteristics of the star obtained by Star_bhns::equil_spher_bhns : "
263 << endl
264 << "------------------------------------------------------------------- "
265 << endl ;
266
267 cout.precision(16) ;
268 double ray = mp.val_r(l_b, 1., M_PI/2., 0) ;
269 cout << "Coordinate radius : "
270 << ray / km << " [km]" << endl ;
271
272 double rcirc = ray * sqrt(psi4.val_grid_point(l_b, k_b, j_b, i_b) ) ;
273 double compact = qpig/(4.*M_PI) * mass_g_bhns() / rcirc ;
274
275 cout << "Circumferential radius R : "
276 << rcirc/km << " [km]" << endl ;
277 cout << "Baryon mass : "
278 << mass_b_bhns(0,0.,1.)/msol << " [Mo]" << endl ;
279 cout << "Gravitational mass M : "
280 << mass_g_bhns()/msol << " [Mo]" << endl ;
281 cout << "Compaction parameter GM/(c^2 R) : " << compact << endl ;
282
283 //----------------
284 // Virial theorem
285 //----------------
286
287 //... Pressure term
288 Scalar source(mp) ;
289 source = qpig * pow(confo_auto,6.) * s_euler ;
290 source.std_spectral_base() ;
291 double vir_mat = source.integrale() ;
292
293 //... Gravitational term
294 Scalar tmp1(mp) ;
295 tmp1 = confo_auto.dsdr() ;
296 tmp1.std_spectral_base() ;
297
298 Scalar tmp2(mp) ;
300 tmp2.std_spectral_base() ;
301
302 source = 2. * tmp1 * tmp1 - tmp2 * tmp2 ;
303
304 /*
305 Scalar tmp1(mp) ;
306 tmp1 = log(lapse_auto) ;
307 tmp1.std_spectral_base() ;
308
309 Scalar tmp2(mp) ;
310 tmp2 = log(confo_auto) ;
311 tmp2.std_spectral_base() ;
312
313 source = confo_auto * confo_auto
314 * ( 2. * tmp2.dsdr() * tmp2.dsdr() - tmp1.dsdr() * tmp1.dsdr() ) ;
315 */
316 source.std_spectral_base() ;
317 double vir_grav = source.integrale() ;
318
319 //... Relative error on the virial identity GRV3
320 double grv3 = ( vir_mat + vir_grav ) / vir_mat ;
321
322 cout << "Virial theorem GRV3 : " << endl ;
323 cout << " 3P term : " << vir_mat << endl ;
324 cout << " grav. term : " << vir_grav << endl ;
325 cout << " relative error : " << grv3 << endl ;
326
327}
328}
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
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
void equil_spher_bhns(double ent_c, double precis)
Computes a spherical configuration.
Scalar lapconf_auto
Lapconf function generated by the star.
Definition star_bhns.h:113
Scalar lapse_tot
Total lapse function.
Definition star_bhns.h:125
Scalar confo_tot
Total conformal factor.
Definition star_bhns.h:163
Scalar psi4
Fourth power of the total conformal factor.
Definition star_bhns.h:176
Scalar lapconf_tot
Total lapconf function.
Definition star_bhns.h:119
Scalar lapse_auto
Lapse function generated by the "star".
Definition star_bhns.h:122
Vector shift_auto
Shift vector generated by the star.
Definition star_bhns.h:138
Scalar confo_auto
Conformal factor generated by the star.
Definition star_bhns.h:157
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
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
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
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
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.