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 char 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 
63 namespace Lorene {
64 void 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 
111  Scalar source_lapconf(mp) ;
112  Scalar source_confo(mp) ;
113  Scalar lapconf_auto_m1(mp) ;
114  Scalar confo_auto_m1(mp) ;
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  // ----------------------
140  source_lapconf = qpig * lapconf_auto * pow(confo_auto,4.)
141  * (0.5*ener + 3.*press) ;
142 
143  source_lapconf.inc_dzpuis(4-source_lapconf.get_dzpuis()) ;
144  source_lapconf.std_spectral_base() ;
145 
146  Cmp sou_lapconf(source_lapconf) ;
147  Cmp lapconf_cmp(lapconf_auto_m1) ;
148  lapconf_cmp.set_etat_qcq() ;
149 
150  mpaff.poisson(sou_lapconf, par_nul, lapconf_cmp) ;
151 
152  // Re-construction of a scalar
153  lapconf_auto_m1 = lapconf_cmp ;
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) ;
166  double confo_b = confo_auto.val_grid_point(l_b,k_b,j_b,i_b) ;
167 
168  double alpha_r2 = (exp_ent_b*confo_c - exp_ent_c*confo_b)
169  / ( exp_ent_c*confo_b*lap_auto_c
170  - exp_ent_b*confo_c*lap_auto_b ) ;
171 
172  alpha_r = sqrt(alpha_r2) ;
173 
174  // New radial scale
175  mpaff.homothetie( alpha_r ) ;
176 
177  //----------------
178  // First integral
179  //----------------
180 
181  // Lapconf function
182  lapconf_auto_m1 = alpha_r2 * lapconf_auto_m1 ;
183  lapconf_auto = lapconf_auto_m1 + 1. ;
184 
185  // Enthalpy in all space
186  double lapconfo_c = lapconf_auto.val_grid_point(0,0,0,0) ;
187  confo_c = confo_auto.val_grid_point(0,0,0,0) ;
188  ent = ent_c + log(lapconfo_c) - log(confo_c)
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 
206  Cmp sou_confo(source_confo) ;
207  Cmp cnf_auto_cmp(confo_auto_m1) ;
208  cnf_auto_cmp.set_etat_qcq() ;
209 
210  mpaff.poisson(sou_confo, par_nul, cnf_auto_cmp) ;
211 
212  // Re-construction of a scalr
213  confo_auto_m1 = cnf_auto_cmp ;
214 
215  confo_auto = confo_auto_m1 + 1. ;
216 
217  // Relative difference with enthalphy at the previous step
218  // -------------------------------------------------------
219 
220  diff_ent = norme( diffrel(ent, ent_jm1) ) / nzet ;
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) ;
299  tmp2 = confo_auto * lapconf_auto.dsdr() / lapconf_auto - tmp1 ;
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 }
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:296
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Map & mp
Mapping associated with the star.
Definition: star.h:180
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_af.C:537
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:61
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
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:504
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:637
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:481
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:557
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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:348
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
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:452
Multi-domain grid.
Definition: grilles.h:273
Affine radial mapping.
Definition: map.h:2027
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:457
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
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