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 
28 
29 
30 /*
31  * $Id: star_bhns_spher.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
32  * $Log: star_bhns_spher.C,v $
33  * Revision 1.5 2016/12/05 16:18:16 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.4 2014/10/13 08:53:41 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.3 2014/10/06 15:13:16 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.2 2008/05/15 19:16:54 k_taniguchi
43  * Change of some parameters.
44  *
45  * Revision 1.1 2007/06/22 01:32:19 k_taniguchi
46  * *** empty log message ***
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_spher.C,v 1.5 2016/12/05 16:18:16 j_novak Exp $
50  *
51  */
52 
53 // C++ headers
54 //#include <>
55 
56 // C headers
57 #include <cmath>
58 
59 // Lorene headers
60 #include "star_bhns.h"
61 #include "param.h"
62 #include "cmp.h"
63 #include "tenseur.h"
64 #include "unites.h"
65 
66 namespace Lorene {
67 void Star_bhns::equil_spher_bhns(double ent_c, double precis) {
68 
69  // Fundamental constants and units
70  // -------------------------------
71  using namespace Unites ;
72 
73  // Initializations
74  // ---------------
75 
76  const Mg3d* mg = mp.get_mg() ;
77  int nz = mg->get_nzone() ;
78 
79  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
80  int l_b = nzet - 1 ;
81  int i_b = mg->get_nr(l_b) - 1 ;
82  int j_b = mg->get_nt(l_b) - 1 ;
83  int k_b = 0 ;
84 
85  // Value of the enthalpy defining the surface of the star
86  double ent_b = 0. ;
87 
88  // Initialization of the enthalpy field to the constant value ent_c :
89  ent = ent_c ;
90  ent.annule(nzet, nz-1) ;
91 
92  // Corresponding profiles of baryon density, energy density and pressure
94 
95  // Initial metric
96  lapconf_auto = 1. ;
98  confo_auto = 1. ;
100 
101  // Auxiliary quantities
102  // --------------------
103 
104  // Affine mapping for solving the Poisson equations
105  Map_af mpaff(mp) ;
106 
107  Param par_nul ; // Param (null) for Map_af::poisson
108 
109  Scalar ent_jm1(mp) ; // Enthalpy at previous step
110  ent_jm1 = ent_c ;
111  ent_jm1.std_spectral_base() ;
112  ent_jm1.annule(nzet, nz-1) ;
113 
114  Scalar source_lapconf(mp) ;
115  Scalar source_confo(mp) ;
116  Scalar lapconf_auto_m1(mp) ;
117  Scalar confo_auto_m1(mp) ;
118  lapconf_auto_m1 = 0. ;
119  confo_auto_m1 = 0. ;
120 
121  double diff_ent = 1. ;
122  int mermax = 200 ; // Maximum number of iterations
123 
124  double alpha_r = 1. ;
125 
126  //==========================================================//
127  // Start of iteration //
128  //==========================================================//
129 
130  for (int mer=0 ; (diff_ent > precis) && (mer<mermax) ; mer++ ) {
131 
132  cout << "-----------------------------------------------" << endl ;
133  cout << "step: " << mer << endl ;
134  cout << "alpha_r: " << alpha_r << endl ;
135  cout << "diff_ent = " << diff_ent << endl ;
136 
137  //----------------------------------------------------
138  // Resolution of Poisson equation for lapconf function
139  //----------------------------------------------------
140 
141  // Matter part of lapconf
142  // ----------------------
143  source_lapconf = qpig * lapconf_auto * pow(confo_auto,4.)
144  * (0.5*ener + 3.*press) ;
145 
146  source_lapconf.inc_dzpuis(4-source_lapconf.get_dzpuis()) ;
147  source_lapconf.std_spectral_base() ;
148 
149  Cmp sou_lapconf(source_lapconf) ;
150  Cmp lapconf_cmp(lapconf_auto_m1) ;
151  lapconf_cmp.set_etat_qcq() ;
152 
153  mpaff.poisson(sou_lapconf, par_nul, lapconf_cmp) ;
154 
155  // Re-construction of a scalar
156  lapconf_auto_m1 = lapconf_cmp ;
157 
158  //-------------------------------------
159  // Computation of the new radial scale
160  //-------------------------------------
161 
162  double exp_ent_c = exp(ent_c) ;
163  double exp_ent_b = exp(ent_b) ;
164 
165  double lap_auto_c = lapconf_auto_m1.val_grid_point(0,0,0,0) ;
166  double lap_auto_b = lapconf_auto_m1.val_grid_point(l_b,k_b,j_b,i_b) ;
167 
168  double confo_c = confo_auto.val_grid_point(0,0,0,0) ;
169  double confo_b = confo_auto.val_grid_point(l_b,k_b,j_b,i_b) ;
170 
171  double alpha_r2 = (exp_ent_b*confo_c - exp_ent_c*confo_b)
172  / ( exp_ent_c*confo_b*lap_auto_c
173  - exp_ent_b*confo_c*lap_auto_b ) ;
174 
175  alpha_r = sqrt(alpha_r2) ;
176 
177  // New radial scale
178  mpaff.homothetie( alpha_r ) ;
179 
180  //----------------
181  // First integral
182  //----------------
183 
184  // Lapconf function
185  lapconf_auto_m1 = alpha_r2 * lapconf_auto_m1 ;
186  lapconf_auto = lapconf_auto_m1 + 1. ;
187 
188  // Enthalpy in all space
189  double lapconfo_c = lapconf_auto.val_grid_point(0,0,0,0) ;
190  confo_c = confo_auto.val_grid_point(0,0,0,0) ;
191  ent = ent_c + log(lapconfo_c) - log(confo_c)
194 
195  //-------------------
196  // Equation of state
197  //-------------------
198 
200 
201  //-----------------------------------------------------
202  // Resolution of Poisson equation for conformal factor
203  //-----------------------------------------------------
204 
205  source_confo = - 0.5 * qpig * pow(confo_auto,5.) * ener ;
206  source_confo.inc_dzpuis(4-source_confo.get_dzpuis()) ;
207  source_confo.std_spectral_base() ;
208 
209  Cmp sou_confo(source_confo) ;
210  Cmp cnf_auto_cmp(confo_auto_m1) ;
211  cnf_auto_cmp.set_etat_qcq() ;
212 
213  mpaff.poisson(sou_confo, par_nul, cnf_auto_cmp) ;
214 
215  // Re-construction of a scalr
216  confo_auto_m1 = cnf_auto_cmp ;
217 
218  confo_auto = confo_auto_m1 + 1. ;
219 
220  // Relative difference with enthalphy at the previous step
221  // -------------------------------------------------------
222 
223  diff_ent = norme( diffrel(ent, ent_jm1) ) / nzet ;
224 
225  // Next step
226  // ---------
227 
228  ent_jm1 = ent ;
229 
230  } // End of iteration loop
231 
232  //========================================================//
233  // End of iteration //
234  //========================================================//
235 
236  // The mapping is transfered to that of the star
237  // ---------------------------------------------
238  mp = mpaff ;
239 
240  // Sets values
241  // -----------
242 
243  // ... hydro
244  ent.annule(nzet, nz-1) ;
245 
246  ener_euler = ener ;
247  s_euler = 3. * press ;
248  gam_euler = 1. ;
249  for(int i=1; i<=3; i++)
250  u_euler.set(i) = 0 ;
251 
252  // ... metric
257  psi4 = pow(confo_auto, 4.) ;
258  for (int i=1; i<=3; i++)
259  shift_auto.set(i) = 0. ;
260 
261  // Info printing
262  // -------------
263 
264  cout << endl
265  << "Characteristics of the star obtained by Star_bhns::equil_spher_bhns : "
266  << endl
267  << "------------------------------------------------------------------- "
268  << endl ;
269 
270  cout.precision(16) ;
271  double ray = mp.val_r(l_b, 1., M_PI/2., 0) ;
272  cout << "Coordinate radius : "
273  << ray / km << " [km]" << endl ;
274 
275  double rcirc = ray * sqrt(psi4.val_grid_point(l_b, k_b, j_b, i_b) ) ;
276  double compact = qpig/(4.*M_PI) * mass_g_bhns() / rcirc ;
277 
278  cout << "Circumferential radius R : "
279  << rcirc/km << " [km]" << endl ;
280  cout << "Baryon mass : "
281  << mass_b_bhns(0,0.,1.)/msol << " [Mo]" << endl ;
282  cout << "Gravitational mass M : "
283  << mass_g_bhns()/msol << " [Mo]" << endl ;
284  cout << "Compaction parameter GM/(c^2 R) : " << compact << endl ;
285 
286  //----------------
287  // Virial theorem
288  //----------------
289 
290  //... Pressure term
291  Scalar source(mp) ;
292  source = qpig * pow(confo_auto,6.) * s_euler ;
293  source.std_spectral_base() ;
294  double vir_mat = source.integrale() ;
295 
296  //... Gravitational term
297  Scalar tmp1(mp) ;
298  tmp1 = confo_auto.dsdr() ;
299  tmp1.std_spectral_base() ;
300 
301  Scalar tmp2(mp) ;
302  tmp2 = confo_auto * lapconf_auto.dsdr() / lapconf_auto - tmp1 ;
303  tmp2.std_spectral_base() ;
304 
305  source = 2. * tmp1 * tmp1 - tmp2 * tmp2 ;
306 
307  /*
308  Scalar tmp1(mp) ;
309  tmp1 = log(lapse_auto) ;
310  tmp1.std_spectral_base() ;
311 
312  Scalar tmp2(mp) ;
313  tmp2 = log(confo_auto) ;
314  tmp2.std_spectral_base() ;
315 
316  source = confo_auto * confo_auto
317  * ( 2. * tmp2.dsdr() * tmp2.dsdr() - tmp1.dsdr() * tmp1.dsdr() ) ;
318  */
319  source.std_spectral_base() ;
320  double vir_grav = source.integrale() ;
321 
322  //... Relative error on the virial identity GRV3
323  double grv3 = ( vir_mat + vir_grav ) / vir_mat ;
324 
325  cout << "Virial theorem GRV3 : " << endl ;
326  cout << " 3P term : " << vir_mat << endl ;
327  cout << " grav. term : " << vir_grav << endl ;
328  cout << " relative error : " << grv3 << endl ;
329 
330 }
331 }
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Map & mp
Mapping associated with the star.
Definition: star.h:180
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:667
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:64
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Scalar ent
Log-enthalpy.
Definition: star.h:190
Vector shift_auto
Shift vector generated by the star.
Definition: star_bhns.h:138
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Scalar lapse_tot
Total lapse function.
Definition: star_bhns.h:125
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:643
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar confo_tot
Total conformal factor.
Definition: star_bhns.h:163
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
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.
Scalar ener
Total energy density in the fluid frame.
Definition: star.h:193
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
virtual void poisson(const Cmp &source, Param &par, Cmp &uu) const
Computes the solution of a scalar Poisson equation.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Scalar press
Fluid pressure.
Definition: star.h:194
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void equil_spher_bhns(double ent_c, double precis)
Computes a spherical configuration.
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
Scalar lapconf_auto
Lapconf function generated by the star.
Definition: star_bhns.h:113
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
Affine radial mapping.
Definition: map.h:2048
Scalar confo_auto
Conformal factor generated by the star.
Definition: star_bhns.h:157
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:465
Scalar lapse_auto
Lapse function generated by the "star".
Definition: star_bhns.h:122
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198