LORENE
strot_cfc_global.C
1 /*
2  * Methods for computing global quantities within the class Star_rot_CFC
3  *
4  * (see file star_rot_cfc.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2024 Jerome Novak
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 // Lorene headers
31 #include "star_rot_cfc.h"
32 #include "unites.h"
33 #include "utilitaires.h"
34 #include "proto.h"
35 
36  //-------------------------------------------------------//
37  // Baryon mass //
38  // //
39  // Note: In Lorene units, mean particle mass is unity //
40  //-------------------------------------------------------//
41 
42 
43 namespace Lorene {
44 double Star_rot_CFC::mass_b() const {
45 
46  if (p_mass_b == 0x0) { // a new computation is required
47 
48  Scalar dens = psi4*psi2 * gam_euler * nbar ;
49 
50  dens.std_spectral_base() ;
51 
52  p_mass_b = new double( dens.integrale() ) ;
53 
54  }
55 
56  return *p_mass_b ;
57 
58 }
59 
60  //---------------------------------------------------------//
61  // Gravitational mass //
62  // //
63  // Note: This is the Komar mass for stationary and //
64  // asymptotically flat spacetime (see, eg, Wald) //
65  //---------------------------------------------------------//
66 
67 double Star_rot_CFC::mass_g() const {
68 
69  if (p_mass_g == 0x0) { // a new computation is required
70 
71  Scalar j_source = 2.* contract(contract(gamma.cov(), 0, j_euler, 0),
72  0, beta, 0) ;
73 
74  Scalar dens = psi4*psi2 * ( nn * (ener_euler + s_euler) - j_source ) ;
75 
76  dens.std_spectral_base() ;
77 
78  p_mass_g = new double( dens.integrale() ) ;
79 
80  }
81 
82  return *p_mass_g ;
83 
84 }
85 
86  //--------------------------------------//
87  // Angular momentum //
88  // //
89  // Komar-type integral (see, eg, Wald) //
90  //--------------------------------------//
91 
92 double Star_rot_CFC::angu_mom() const {
93 
94  if (p_angu_mom == 0x0) { // a new computation is required
95 
96  // phi_kill = axial killing vector
97 
98  Vector phi_kill(mp, CON, mp.get_bvect_spher()) ;
99 
100  phi_kill.set(1).set_etat_zero() ;
101  phi_kill.set(2).set_etat_zero() ;
102  phi_kill.set(3) = 1. ;
103  phi_kill.set(3).std_spectral_base() ;
104  phi_kill.set(3).mult_rsint() ;
105 
106  Scalar j_source = contract(contract(gamma.cov(), 0, j_euler, 0),
107  0, phi_kill, 0) ;
108 
109  Scalar dens = psi4*psi2 * j_source ;
110 
111  dens.std_spectral_base() ;
112 
113  p_angu_mom = new double( dens.integrale() ) ;
114 
115 
116  }
117 
118  return *p_angu_mom ;
119 
120 }
121 
122  //---------------------//
123  // T/W //
124  //---------------------//
125 
126 double Star_rot_CFC::tsw() const {
127 
128  if (p_tsw == 0x0) { // a new computation is required
129 
130  double tcin = 0.5 * omega * angu_mom() ;
131 
132  Scalar dens = psi4*psi2 * gam_euler * ener ;
133 
134  dens.std_spectral_base() ;
135 
136  double mass_p = dens.integrale() ;
137 
138  p_tsw = new double( tcin / ( mass_p + tcin - mass_g() ) ) ;
139 
140  }
141 
142  return *p_tsw ;
143 
144 }
145 
146 
147  //--------------------------------------------------------------//
148  // GRV2 //
149  // cf. Eq. (28) of Bonazzola & Gourgoulhon CQG, 11, 1775 (1994) //
150  // //
151  //--------------------------------------------------------------//
152 
153 double Star_rot_CFC::grv2() const {
154 
155  using namespace Unites ;
156 
157  if (p_grv2 == 0x0) { // a new computation is required
158 
159  // determinant of the 2-metric k_{ab}
160 
161  Scalar k_det = gamma.cov()(1,1)*gamma.cov()(2,2) -
162  gamma.cov()(1,2)*gamma.cov()(1,2) ;
163 
164 
165  //**
166  // sou_m = 8\pi T_{\mu\nu} m^{\mu}m^{\nu}
167  // => sou_m = 8\pi [ (E+P) U^2 + P ], where v^2 = v_i v^i
168  //
169  //**
170 
171  Scalar sou_m = 2 * qpig * ( (ener_euler + press)*v2 + press ) ;
172 
173  sou_m = sqrt( k_det )*sou_m ;
174 
175  sou_m.std_spectral_base() ;
176 
177 
178  // This is the term 3k_a k^a.
179 
180  Scalar sou_q = 3 *( hatA(1,3) * hatA(1,3) + hatA(2,3)*hatA(2,3) )
181  / (psi4*psi4*psi4) ;
182 
183 
184  // This is the term \nu_{|| a}\nu^{|| a}.
185  //
186 
187  Scalar sou_tmp = gamma.con()(1,1) * logn.dsdr() * logn.dsdr() ;
188 
189  Scalar term_2 = 2 * gamma.con()(1,2) * logn.dsdr() * logn.dsdt() ;
190 
191  term_2.div_r_dzpuis(4) ;
192 
193  Scalar term_3 = gamma.con()(2,2) * logn.dsdt() * logn.dsdt() ;
194 
195  term_3.div_r_dzpuis(2) ;
196  term_3.div_r_dzpuis(4) ;
197 
198  sou_tmp += term_2 + term_3 ;
199 
200 
201  // Source of the gravitational part
202 
203  sou_q -= sou_tmp ;
204 
205  sou_q = sqrt( k_det )*sou_q ;
206 
207  sou_q.std_spectral_base() ;
208 
209  p_grv2 = new double( double(1) + integrale2d(sou_m)/integrale2d(sou_q) ) ;
210 
211  }
212 
213  return *p_grv2 ;
214 
215 }
216 
217 
218  //-------------------------------------------------------------//
219  // GRV3 //
220  // cf. Eq. (29) of Gourgoulhon & Bonazzola CQG, 11, 443 (1994) //
221  //-------------------------------------------------------------//
222 
223 double Star_rot_CFC::grv3() const {
224 
225  using namespace Unites ;
226 
227  if (p_grv3 == 0x0) { // a new computation is required
228 
229  // Gravitational term
230  // -------------------
231 
232  Scalar psi6 = psi4*psi2 ;
233  Scalar sou_q = 0.75*hatA_quad/psi6
234  - psi6*contract(logn.derive_cov(gamma), 0,logn.derive_con(gamma), 0) ;
235 
236 
237  Tensor t_tmp = contract(gamma.connect().get_delta(), 2,
238  gamma.connect().get_delta(), 0) ;
239 
240  Scalar tmp_1 = 0.25* contract( gamma.con(), 0, 1,
241  contract(t_tmp, 0, 3), 0, 1 ) ;
242 
243  Scalar tmp_2 = 0.25* contract( gamma.con(), 0, 1,
244  contract( contract( gamma.connect().get_delta(), 0, 1),
245  0, gamma.connect().get_delta(), 0), 0, 1) ;
246 
247  sou_q = sou_q + psi6*(tmp_1 - tmp_2) ;
248 
249  // sou_q = psi4*psi2 * sou_q ;
250 
251  sou_q.std_spectral_base() ;
252 
253  double int_grav = sou_q.integrale() ;
254 
255 
256  // Matter term
257  // --------------
258 
259  Scalar sou_m = qpig*s_euler ;
260 
261  sou_m = psi6 * sou_m ;
262 
263  sou_m.std_spectral_base() ;
264 
265  double int_mat = sou_m.integrale() ;
266 
267  p_grv3 = new double( (int_grav + int_mat) / int_mat ) ;
268 
269 
270 
271  }
272 
273  return *p_grv3 ;
274 
275 }
276 
277 
278  //--------------------//
279  // R_circ //
280  //--------------------//
281 
282 double Star_rot_CFC::r_circ() const {
283 
284  if (p_r_circ ==0x0) { // a new computation is required
285 
286  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
287  const Mg3d* mg = mp.get_mg() ;
288  int l_b = nzet - 1 ;
289  int i_b = mg->get_nr(l_b) - 1 ;
290  int j_b = (mg->get_type_t() == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b) / 2) ;
291  int k_b = 0 ;
292 
293  double gamma_phi = gamma.cov()(3,3).val_grid_point(l_b, k_b, j_b, i_b) ;
294 
295  p_r_circ = new double( sqrt( gamma_phi ) * ray_eq() ) ;
296 
297  }
298 
299  return *p_r_circ ;
300 
301 }
302 
303 
304  //--------------------//
305  // R_circ //
306  //--------------------//
307 
308 double Star_rot_CFC::rp_circ() const {
309 
310  if (p_rp_circ ==0x0) { // a new computation is required
311  const Mg3d& mg = *mp.get_mg() ;
312  int nz = mg.get_nzone() ;
313  assert(nz>2) ;
314  int np = mg.get_np(0) ;
315  if (np != 1) {
316  cout << "The polar circumferential radius is only well defined\n"
317  << "with np = 1!" << endl ;
318  abort() ;
319  }
320  int nt = mg.get_nt(0) ;
321  Sym_tensor gam = gamma.cov() ;
322  const Coord& tet = mp.tet ;
323  Scalar rrr(mp) ;
324  rrr.annule_hard() ;
325  rrr.annule(nzet,nz-1) ;
326  double phi = 0 ;
327  for (int j=0; j<nt; j++) {
328  double theta = (+tet)(0, 0, j, 0) ;
329  double erre =
330  mp.val_r(l_surf()(0,j), xi_surf()(0, j), theta, phi) ;
331  for (int lz=0; lz<nzet; lz++) {
332  int nrz = mg.get_nr(lz) ;
333  for (int i=0; i<nrz; i++) {
334  rrr.set_grid_point(lz, 0, j, i) = erre ;
335  }
336  }
337  }
338  rrr.std_spectral_base() ;
339  Scalar drrr = rrr.dsdt() ;
340 
341  Scalar fff(mp) ;
342  fff.annule_hard() ;
343  fff.annule(nzet,nz-1) ;
344  for (int j=0; j<nt; j++) {
345  double theta = (+tet)(0, 0, j, 0) ;
346  int ls = l_surf()(0, j) ;
347  double xs = xi_surf()(0, j) ;
348  double grr = gam(1,1).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
349  double grt = gam(1,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
350  double gtt = gam(2,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
351  double rr = mp.val_r(ls, xs, theta, phi) ;
352  double dr = drrr.get_spectral_va().val_point_jk(ls, xs, j, 0) ;
353  for (int i=0; i< mg.get_nr(nzet-1); i++) {
354  fff.set_grid_point(nzet-1, 0, j, i)
355  = sqrt(grr*dr*dr + 2*grt*rr*dr + gtt*rr*rr) ;
356  }
357  }
358  fff.std_spectral_base() ;
359  fff.set_spectral_va().coef() ;
360  p_rp_circ = new double((*fff.get_spectral_va().c_cf)(nzet-1, 0, 0, 0)) ;
361  }
362  return *p_rp_circ ;
363 }
364 
365  //--------------------------//
366  // Flattening //
367  //--------------------------//
368 
369 double Star_rot_CFC::aplat() const {
370 
371  return ray_pole() / ray_eq() ;
372 
373 }
374 
375  //--------------------------//
376  // Ellipticity //
377  //--------------------------//
378 
379 double Star_rot_CFC::ellipt() const {
380 
381  return sqrt(1. - (rp_circ()*rp_circ())/(r_circ()*r_circ())) ;
382 
383 }
384 }
double * p_grv2
Error on the virial identity GRV2.
Definition: star_rot_cfc.h:98
double * p_mass_b
Baryon mass.
Definition: star.h:268
virtual double tsw() const
Ratio T/W.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Metric gamma
3-metric
Definition: star.h:235
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
virtual double grv3() const
Error on the virial identity GRV3.
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
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
virtual double angu_mom() const
Angular momentum.
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:64
virtual double r_circ() const
Circumferential equatorial radius.
double * p_rp_circ
Circumferential polar radius.
Definition: star_rot_cfc.h:102
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
Scalar psi4
Conformal factor .
Definition: star_rot_cfc.h:81
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 nbar
Baryon density in the fluid frame.
Definition: star.h:192
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Vector beta
Shift vector.
Definition: star.h:228
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition: connection.h:271
Coord tet
coordinate centered on the grid
Definition: map.h:737
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinate ...
Definition: star_global.C:92
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
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
virtual double aplat() const
Flattening r_pole/r_eq.
double omega
Rotation angular velocity ([f_unit] )
Definition: star_rot_cfc.h:60
Scalar press
Fluid pressure.
Definition: star.h:194
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:304
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l on...
Definition: star_global.C:66
virtual double rp_circ() const
Circumferential polar radius.
double * p_angu_mom
Angular momentum.
Definition: star_rot_cfc.h:97
double * p_r_circ
Circumferential equatorial radius.
Definition: star_rot_cfc.h:101
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual double mass_b() const
Baryonic mass.
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Definition: valeur.C:903
Tensor handling.
Definition: tensor.h:294
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:281
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Multi-domain grid.
Definition: grilles.h:279
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
Scalar nn
Lapse function N .
Definition: star.h:225
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
double * p_tsw
Ratio T/W.
Definition: star_rot_cfc.h:100
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
double * p_grv3
Error on the virial identity GRV3.
Definition: star_rot_cfc.h:99
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
virtual double grv2() const
Error on the virial identity GRV2.
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Vector j_euler
Momentum density 3-vector with respect to the Eulerian observer.
Definition: star_rot_cfc.h:75
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
double * p_mass_g
Gravitational mass.
Definition: star.h:269
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
virtual double mass_g() const
Gravitational mass.
virtual double ellipt() const
Ellipticity e.