LORENE
et_equil_spher_regu.C
1/*
2 * Method of class Etoile to compute a static spherical configuration
3 * by regularizing source.
4 *
5 * (see file etoile.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2000-2001 Eric Gourgoulhon
11 * Copyright (c) 2000-2001 Keisuke Taniguchi
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 et_equil_spher_regu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_equil_spher_regu.C,v 1.4 2014/10/13 08:52:56 j_novak Exp $" ;
33
34/*
35 * $Id: et_equil_spher_regu.C,v 1.4 2014/10/13 08:52:56 j_novak Exp $
36 * $Log: et_equil_spher_regu.C,v $
37 * Revision 1.4 2014/10/13 08:52:56 j_novak
38 * Lorene classes and functions now belong to the namespace Lorene.
39 *
40 * Revision 1.3 2014/10/06 15:13:08 j_novak
41 * Modified #include directives to use c++ syntax.
42 *
43 * Revision 1.2 2004/03/25 10:29:04 j_novak
44 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
45 *
46 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
47 * LORENE
48 *
49 * Revision 2.14 2001/03/06 16:34:04 keisuke
50 * Change the regularization degree k_div=1 to the arbitrary one.
51 *
52 * Revision 2.13 2001/02/07 09:46:11 eric
53 * unsgam1 est desormais donne par Eos::der_nbar_ent (cas newtonien)
54 * ou Eos::der_ener_ent (cas relativiste).
55 *
56 * Revision 2.12 2000/09/26 15:42:35 keisuke
57 * Correction erreur: la triade de duu_div doit etre celle de mp et non celle
58 * de l'objet temporaire mpaff.
59 *
60 * Revision 2.11 2000/09/26 06:56:04 keisuke
61 * Suppress "int" from the declaration of k_div.
62 *
63 * Revision 2.10 2000/09/22 15:53:36 keisuke
64 * d_logn_auto_div est desormais un membre de la classe Etoile.
65 *
66 * Revision 2.9 2000/09/08 12:23:17 keisuke
67 * Correct a typological error in the equation.
68 *
69 * Revision 2.8 2000/09/07 15:37:23 keisuke
70 * Add a new argument in poisson_regular
71 * and suppress logn_auto_total.
72 *
73 * Revision 2.7 2000/09/06 12:47:35 keisuke
74 * Suppress #include "graphique.h".
75 *
76 * Revision 2.6 2000/09/06 12:39:59 keisuke
77 * Save the map in the every iterative step after operating "homothetie".
78 *
79 * Revision 2.5 2000/09/04 16:15:05 keisuke
80 * Change the argument of Map_af::poisson_regular.
81 *
82 * Revision 2.4 2000/08/31 16:05:26 keisuke
83 * Modify the arguments of Map_af::poisson_regular.
84 *
85 * Revision 2.3 2000/08/29 11:38:25 eric
86 * Ajout des membres k_div et logn_auto_div a la classe Etoile.
87 *
88 * Revision 2.2 2000/08/28 16:10:39 keisuke
89 * Add "nzet" in the argumant of poisson_regular.
90 *
91 * Revision 2.1 2000/08/25 14:58:15 keisuke
92 * Modif (Virial theorem).
93 *
94 * Revision 2.0 2000/08/25 09:01:33 keisuke
95 * *** empty log message ***
96 *
97 *
98 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_equil_spher_regu.C,v 1.4 2014/10/13 08:52:56 j_novak Exp $
99 *
100 */
101
102// Headers C
103#include <cmath>
104
105// Headers Lorene
106#include "etoile.h"
107#include "eos.h"
108#include "param.h"
109#include "unites.h"
110
111//********************************************************************
112
113namespace Lorene {
114
115void Etoile::equil_spher_regular(double ent_c, double precis){
116
117 // Fundamental constants and units
118 // -------------------------------
119 using namespace Unites ;
120
121 // Initializations
122 // ---------------
123
124 // k_div = 1 ; // Regularization index
125
126 cout << "Input the regularization degree (k_div) : " ;
127 cin >> k_div ; // Regularization index
128
129 const Mg3d* mg = mp.get_mg() ;
130 int nz = mg->get_nzone() ; // total number of domains
131
132 // Index of the point at phi=0, theta=pi/2 at the surface of the star:
133 int l_b = nzet - 1 ;
134 int i_b = mg->get_nr(l_b) - 1 ;
135 int j_b = mg->get_nt(l_b) - 1 ;
136 int k_b = 0 ;
137
138 // Value of the enthalpy defining the surface of the star
139 double ent_b = 0 ;
140
141 // Initialization of the enthalpy field to the constant value ent_c :
142
143 ent = ent_c ;
144 ent.annule(nzet, nz-1) ;
145
146 // Corresponding profiles of baryon density, energy density and pressure
147
149
150 // Initial metric
151 a_car = 1 ; // this value will remain unchanged in the Newtonian case
152 beta_auto = 0 ; // this value will remain unchanged in the Newtonian case
153
154 // Auxiliary quantities
155 // --------------------
156
157 // Affine mapping for solving the Poisson equations
158 Map_af mpaff(mp) ;
159
160 Param par_nul ; // Param (null) for Map_af::poisson.
161
162 Tenseur ent_jm1(mp) ; // Enthalpy at previous step
163 ent_jm1 = 0 ;
164
165 Tenseur source(mp) ;
166 Tenseur logn(mp) ;
167 Tenseur logn_quad(mp) ;
168 logn = 0 ;
169 logn_quad = 0 ;
170
171 Tenseur dlogn_auto_regu(mp, 1, COV, mp.get_bvect_spher()) ;
172
173 Cmp source_regu(mp) ;
174 Cmp source_div(mp) ;
175
176 Cmp dlogn(mp) ;
177 dlogn = 0 ;
178 Cmp dbeta(mp) ;
179
180 Tenseur dlogn_auto(mp, 1, COV, mp.get_bvect_spher()) ;
181 Tenseur dlogn_quad(mp) ;
182 dlogn_quad = 0 ;
183
184 double diff_ent = 1 ;
185 int mermax = 200 ; // Max number of iterations
186
187 double alpha_r = 1 ;
188
189 //=========================================================================
190 // Start of iteration
191 //=========================================================================
192
193 for(int mer=0 ; (diff_ent > precis) && (mer<mermax) ; mer++ ) {
194
195 cout << "-----------------------------------------------" << endl ;
196 cout << "step: " << mer << endl ;
197 cout << "alpha_r: " << alpha_r << endl ;
198 cout << "diff_ent = " << diff_ent << endl ;
199
200 //-----------------------------------------------------
201 // Resolution of Poisson equation for ln(N)
202 //-----------------------------------------------------
203
204 // Matter part of ln(N)
205 // --------------------
206
207 double unsgam1 ; // effective power of H in the source
208 // close to the surface
209
210 if (relativistic) {
211
212 source = a_car * (ener + 3*press) ;
213
214 // 1/(gam-1) = dln(e)/dln(H) close to the surface :
215 unsgam1 = eos.der_ener_ent_p(ent_b + 1e-10*(ent_c-ent_b)) ;
216 }
217 else {
218
219 source = nbar ;
220
221 // 1/(gam-1) = dln(n)/dln(H) close to the surface :
222 unsgam1 = eos.der_nbar_ent_p(ent_b + 1e-10*(ent_c-ent_b)) ;
223 }
224
225 (source.set()).set_dzpuis(4) ;
226
227 source.set_std_base() ; // Sets the standard spectral bases.
228
232
233 source_regu.std_base_scal() ;
234 source_div.std_base_scal() ;
235
236 mpaff.poisson_regular(source(), k_div, nzet, unsgam1, par_nul,
239 d_logn_auto_div, source_regu, source_div) ;
240
241 dlogn_auto_regu = logn_auto_regu.gradient_spher() ;
242
243 dlogn_auto = dlogn_auto_regu + d_logn_auto_div ;
244
245 // NB: at this stage logn_auto is in machine units, not in c^2
246
247 // Quadratic part of ln(N)
248 // -----------------------
249
250 if (relativistic) {
251
252 mpaff.dsdr(beta_auto(), dbeta) ;
253
254 source = - dlogn * dbeta ;
255
256 logn_quad.set_etat_qcq() ;
257
258 mpaff.poisson(source(), par_nul, logn_quad.set()) ;
259
260 dlogn_quad.set_etat_qcq() ;
261
262 mpaff.dsdr(logn_quad(), dlogn_quad.set()) ;
263
264 }
265
266 //-----------------------------------------------------
267 // Computation of the new radial scale
268 //-----------------------------------------------------
269
270 // alpha_r (r = alpha_r r') is determined so that the enthalpy
271 // takes the requested value ent_b at the stellar surface
272
273 double nu_mat0_b = logn_auto()(l_b, k_b, j_b, i_b) ;
274 double nu_mat0_c = logn_auto()(0, 0, 0, 0) ;
275
276 double nu_quad0_b = logn_quad()(l_b, k_b, j_b, i_b) ;
277 double nu_quad0_c = logn_quad()(0, 0, 0, 0) ;
278
279 double alpha_r2 = ( ent_c - ent_b - nu_quad0_b + nu_quad0_c )
280 / ( qpig*(nu_mat0_b - nu_mat0_c) ) ;
281
282 alpha_r = sqrt(alpha_r2) ;
283
284 // New radial scale
285 // -----------------
286 mpaff.homothetie( alpha_r ) ;
287
288 // The mapping is transfered to that of the star:
289 // ----------------------------------------------
290 mp = mpaff ;
291
292 d_logn_auto_div.set_triad( mp.get_bvect_spher() ) ; // Absolutely necessary !!!
293
294 //--------------------
295 // First integral
296 //--------------------
297
298 // Gravitation potential in units c^2 :
299 logn_auto = alpha_r2*qpig * logn_auto ;
300 logn = logn_auto + logn_quad ;
301
302 // Enthalpy in all space
303 double logn_c = logn()(0, 0, 0, 0) ;
304 ent = ent_c - logn() + logn_c ;
305
306 //---------------------
307 // Equation of state
308 //---------------------
309
311
312 // derivative of gravitation potential in units c^2 :
313 dlogn_auto = alpha_r*qpig * dlogn_auto ;
314 dlogn = dlogn_auto(0) + dlogn_quad() ;
315
316 if (relativistic) {
317
318 //----------------------------
319 // Equation for beta = ln(AN)
320 //----------------------------
321
322 mpaff.dsdr(beta_auto(), dbeta) ;
323
324 source = 3 * qpig * a_car * press ;
325
326 source = source()
327 - 0.5 * ( dlogn * dlogn + dbeta * dbeta ) ;
328
329 source.set_std_base() ; // Sets the standard spectral bases.
330
332
333 mpaff.poisson(source(), par_nul, beta_auto.set()) ;
334
335 // Metric coefficient A^2 update
336
337 a_car = exp(2*(beta_auto - logn)) ;
338
339 }
340
341 // Relative difference with enthalpy at the previous step
342 // ------------------------------------------------------
343
344 diff_ent = norme( diffrel(ent(), ent_jm1()) ) / nzet ;
345
346 // Next step
347 // ---------
348
349 ent_jm1 = ent ;
350
351
352 } // End of iteration loop
353
354 //=========================================================================
355 // End of iteration
356 //=========================================================================
357
358
359 // Sets value to all the Tenseur's of the star
360 // -------------------------------------------
361
362 // ... hydro
363 ent.annule(nzet, nz-1) ; // enthalpy set to zero at the exterior of
364 // the star
365 ener_euler = ener ;
366 s_euler = 3 * press ;
367 gam_euler = 1 ;
368 u_euler = 0 ;
369
370 // ... metric
371 nnn = exp( unsurc2 * logn ) ;
372 shift = 0 ;
373
374 // Info printing
375 // -------------
376
377 cout << endl
378 << "Characteristics of the star obtained by Etoile::equil_spher_regular : "
379 << endl
380 << "-----------------------------------------------------------------"
381 << endl ;
382
383 double ray = mp.val_r(l_b, 1., M_PI/2., 0) ;
384 cout << "Coordinate radius : " << ray / km << " km" << endl ;
385
386 double rcirc = ray * sqrt( a_car()(l_b, k_b, j_b, i_b) ) ;
387
388 double compact = qpig/(4.*M_PI) * mass_g() / rcirc ;
389
390 cout << "Circumferential radius R : " << rcirc/km << " km" << endl ;
391 cout << "Baryon mass : " << mass_b()/msol << " Mo" << endl ;
392 cout << "Gravitational mass M : " << mass_g()/msol << " Mo" << endl ;
393 cout << "Compacity parameter GM/(c^2 R) : " << compact << endl ;
394
395
396 //-----------------
397 // Virial theorem
398 //-----------------
399
400 //... Pressure term
401
402 source = qpig * a_car * sqrt(a_car) * s_euler ;
403 source.set_std_base() ;
404 double vir_mat = source().integrale() ;
405
406 //... Gravitational term
407
408 Cmp tmp = beta_auto().dsdr() - dlogn ;
409
410 source = - ( dlogn * dlogn - 0.5 * tmp * tmp ) * sqrt(a_car()) ;
411
412 source.set_std_base() ;
413 double vir_grav = source().integrale() ;
414
415 //... Relative error on the virial identity GRV3
416
417 double grv3 = ( vir_mat + vir_grav ) / vir_mat ;
418
419 cout << "Virial theorem GRV3 : " << endl ;
420 cout << " 3P term : " << vir_mat << endl ;
421 cout << " grav. term : " << vir_grav << endl ;
422 cout << " relative error : " << grv3 << endl ;
423
424
425}
426
427}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy with extra parameters (virtual function im...
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const =0
Computes the logarithmic derivative from the log-enthalpy and extra parameters (virtual function imp...
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by t...
Definition etoile.h:497
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:491
int k_div
Index of regularity of the gravitational potential logn_auto .
Definition etoile.h:449
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
const Eos & eos
Equation of state of the stellar matter.
Definition etoile.h:451
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
virtual double mass_b() const
Baryon mass.
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
void equil_spher_regular(double ent_c, double precis=1.e-14)
Computes a spherical static configuration.
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
Definition etoile.h:501
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
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
virtual void poisson_regular(const Cmp &source, int k_div, int nzet, double unsgam1, Param &par, Cmp &uu, Cmp &uu_regu, Cmp &uu_div, Tenseur &duu_div, Cmp &source_regu, Cmp &source_div) const
Computes the solution of a scalar Poisson equation.
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.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
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 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
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition tenseur.C:1548
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:674
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.