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 
32 char 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 
113 namespace Lorene {
114 
115 void 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,
238  logn_auto_div.set(),
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 }
int k_div
Index of regularity of the gravitational potential logn_auto .
Definition: etoile.h:449
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:900
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1548
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...
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:537
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
virtual double mass_g() const
Gravitational mass.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
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
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:468
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.
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: etoile.h:442
Tenseur press
Fluid pressure.
Definition: etoile.h:461
Tenseur shift
Total shift vector.
Definition: etoile.h:512
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by the...
Definition: etoile.h:497
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.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
virtual void poisson(const Cmp &source, Param &par, Cmp &uu) const
Computes the solution of a scalar Poisson equation.
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
Parameter storage.
Definition: param.h:125
void equil_spher_regular(double ent_c, double precis=1.e-14)
Computes a spherical static configuration.
virtual double mass_b() const
Baryon mass.
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:566
const Eos & eos
Equation of state of the stellar matter.
Definition: etoile.h:451
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:437
Multi-domain grid.
Definition: grilles.h:273
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:460
Affine radial mapping.
Definition: map.h:2027
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition: etoile.h:484
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:457
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition: etoile.h:506
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:465
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
Definition: etoile.h:501
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 void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:279
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298