LORENE
strot_dirac_global.C
1/*
2 * Methods for computing global quantities within the class Star_rot_Dirac
3 *
4 * (see file star.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005 Lap-Ming Lin & 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
28char strot_dirac_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $" ;
29
30/*
31 * $Id: strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
32 * $Log: strot_dirac_global.C,v $
33 * Revision 1.13 2014/10/13 08:53:40 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.12 2014/10/06 15:13:18 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.11 2010/11/17 11:21:52 j_novak
40 * Corrected minor problems in the case nz=1 and nt=1.
41 *
42 * Revision 1.10 2010/10/22 08:08:40 j_novak
43 * Removal of the method Star_rot_dirac::lambda_grv2() and call to the C++ version of integrale2d.
44 *
45 * Revision 1.9 2009/10/26 10:54:33 j_novak
46 * Added the case of a NONSYM base in theta.
47 *
48 * Revision 1.8 2008/05/30 08:27:38 j_novak
49 * New global quantities rp_circ and ellipt (circumferential polar coordinate and
50 * ellipticity).
51 *
52 * Revision 1.7 2008/02/18 13:51:20 j_novak
53 * Change of a dzpuis
54 *
55 * Revision 1.6 2005/03/25 14:11:49 j_novak
56 * In the 1D case, GRV2 returns -1 (because of a problem in integral2d).
57 *
58 * Revision 1.5 2005/02/17 17:30:42 f_limousin
59 * Change the name of some quantities to be consistent with other classes
60 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
61 *
62 * Revision 1.4 2005/02/09 13:36:01 lm_lin
63 *
64 * Add the calculations of GRV2, T/W, R_circ, and flattening.
65 *
66 * Revision 1.3 2005/02/02 10:11:24 j_novak
67 * Better calculation of the GRV3 identity.
68 *
69 *
70 *
71 * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
72 *
73 */
74
75
76// C headers
77#include <cmath>
78#include <cassert>
79
80// Lorene headers
81#include "star_rot_dirac.h"
82#include "unites.h"
83#include "utilitaires.h"
84#include "proto.h"
85
86 //-------------------------------------------------------//
87 // Baryonic mass //
88 // //
89 // Note: In Lorene units, mean particle mass is unity //
90 //-------------------------------------------------------//
91
92
93namespace Lorene {
94double Star_rot_Dirac::mass_b() const {
95
96 if (p_mass_b == 0x0) { // a new computation is required
97
99
100 dens.std_spectral_base() ;
101
102 p_mass_b = new double( dens.integrale() ) ;
103
104 }
105
106 return *p_mass_b ;
107
108}
109
110 //---------------------------------------------------------//
111 // Gravitational mass //
112 // //
113 // Note: This is the Komar mass for stationary and //
114 // asymptotically flat spacetime (see, eg, Wald) //
115 //---------------------------------------------------------//
116
118
119 if (p_mass_g == 0x0) { // a new computation is required
120
122 0, beta, 0) ;
123
125 ( nn * (ener_euler + s_euler) - j_source ) ;
126
127 dens.std_spectral_base() ;
128
129 p_mass_g = new double( dens.integrale() ) ;
130
131 }
132
133 return *p_mass_g ;
134
135}
136
137 //--------------------------------------//
138 // Angular momentum //
139 // //
140 // Komar-type integral (see, eg, Wald) //
141 //--------------------------------------//
142
144
145 if (p_angu_mom == 0x0) { // a new computation is required
146
147 // phi_kill = axial killing vector
148
150
151 phi_kill.set(1).set_etat_zero() ;
152 phi_kill.set(2).set_etat_zero() ;
153 phi_kill.set(3) = 1. ;
154 phi_kill.set(3).std_spectral_base() ;
155 phi_kill.set(3).mult_rsint() ;
156
158 0, phi_kill, 0) ;
159
161
162 dens.std_spectral_base() ;
163
164 p_angu_mom = new double( dens.integrale() ) ;
165
166
167 }
168
169 return *p_angu_mom ;
170
171}
172
173 //---------------------//
174 // T/W //
175 //---------------------//
176
177double Star_rot_Dirac::tsw() const {
178
179 if (p_tsw == 0x0) { // a new computation is required
180
181 double tcin = 0.5 * omega * angu_mom() ;
182
184
185 dens.std_spectral_base() ;
186
187 double mass_p = dens.integrale() ;
188
189 p_tsw = new double( tcin / ( mass_p + tcin - mass_g() ) ) ;
190
191 }
192
193 return *p_tsw ;
194
195}
196
197
198 //--------------------------------------------------------------//
199 // GRV2 //
200 // cf. Eq. (28) of Bonazzola & Gourgoulhon CQG, 11, 1775 (1994) //
201 // //
202 //--------------------------------------------------------------//
203
204double Star_rot_Dirac::grv2() const {
205
206 using namespace Unites ;
207
208 if (p_grv2 == 0x0) { // a new computation is required
209
210 // determinant of the 2-metric k_{ab}
211
212 Scalar k_det = gamma.cov()(1,1)*gamma.cov()(2,2) -
213 gamma.cov()(1,2)*gamma.cov()(1,2) ;
214
215
216 //**
217 // sou_m = 8\pi T_{\mu\nu} m^{\mu}m^{\nu}
218 // => sou_m = 8\pi [ (E+P) U^2 + P ], where v^2 = v_i v^i
219 //
220 //**
221
222 Scalar sou_m = 2 * qpig * ( (ener_euler + press)*v2 + press ) ;
223
224 sou_m = sqrt( k_det )*sou_m ;
225
226 sou_m.std_spectral_base() ;
227
228
229 // This is the term 3k_a k^a.
230
231 Scalar sou_q = 3 *( taa(1,3) * aa(1,3) + taa(2,3)*aa(2,3) ) ;
232
233
234 // This is the term \nu_{|| a}\nu^{|| a}.
235 //
236
237 Scalar sou_tmp = gamma.con()(1,1) * logn.dsdr() * logn.dsdr() ;
238
239 Scalar term_2 = 2 * gamma.con()(1,2) * logn.dsdr() * logn.dsdt() ;
240
241 term_2.div_r_dzpuis(4) ;
242
243 Scalar term_3 = gamma.con()(2,2) * logn.dsdt() * logn.dsdt() ;
244
245 term_3.div_r_dzpuis(2) ;
246 term_3.div_r_dzpuis(4) ;
247
248 sou_tmp += term_2 + term_3 ;
249
250
251 // Source of the gravitational part
252
253 sou_q -= sou_tmp ;
254
255 sou_q = sqrt( k_det )*sou_q ;
256
257 sou_q.std_spectral_base() ;
258
259 p_grv2 = new double( double(1) + integrale2d(sou_m)/integrale2d(sou_q) ) ;
260
261 }
262
263 return *p_grv2 ;
264
265}
266
267
268 //-------------------------------------------------------------//
269 // GRV3 //
270 // cf. Eq. (29) of Gourgoulhon & Bonazzola CQG, 11, 443 (1994) //
271 //-------------------------------------------------------------//
272
273double Star_rot_Dirac::grv3() const {
274
275 using namespace Unites ;
276
277 if (p_grv3 == 0x0) { // a new computation is required
278
279 // Gravitational term
280 // -------------------
281
283 logn.derive_con(gamma), 0) ;
284
285
286 Tensor t_tmp = contract(gamma.connect().get_delta(), 2,
287 gamma.connect().get_delta(), 0) ;
288
289 Scalar tmp_1 = 0.25* contract( gamma.con(), 0, 1,
290 contract(t_tmp, 0, 3), 0, 1 ) ;
291
292 Scalar tmp_2 = 0.25* contract( gamma.con(), 0, 1,
293 contract( contract( gamma.connect().get_delta(), 0, 1),
294 0, gamma.connect().get_delta(), 0), 0, 1) ;
295
296 sou_q = sou_q + tmp_1 - tmp_2 ;
297
298 sou_q = sqrt( gamma.determinant() ) * sou_q ;
299
300 sou_q.std_spectral_base() ;
301
302 double int_grav = sou_q.integrale() ;
303
304
305 // Matter term
306 // --------------
307
308 Scalar sou_m = qpig*s_euler ;
309
310 sou_m = sqrt( gamma.determinant() ) * sou_m ;
311
312 sou_m.std_spectral_base() ;
313
314 double int_mat = sou_m.integrale() ;
315
316 p_grv3 = new double( (int_grav + int_mat) / int_mat ) ;
317
318
319
320 }
321
322 return *p_grv3 ;
323
324}
325
326
327 //--------------------//
328 // R_circ //
329 //--------------------//
330
332
333 if (p_r_circ ==0x0) { // a new computation is required
334
335 // Index of the point at phi=0, theta=pi/2 at the surface of the star:
336 const Mg3d* mg = mp.get_mg() ;
337 int l_b = nzet - 1 ;
338 int i_b = mg->get_nr(l_b) - 1 ;
339 int j_b = (mg->get_type_t() == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b) / 2) ;
340 int k_b = 0 ;
341
342 double gamma_phi = gamma.cov()(3,3).val_grid_point(l_b, k_b, j_b, i_b) ;
343
344 p_r_circ = new double( sqrt( gamma_phi ) * ray_eq() ) ;
345
346 }
347
348 return *p_r_circ ;
349
350}
351
352
353 //--------------------//
354 // R_circ //
355 //--------------------//
356
358
359 if (p_rp_circ ==0x0) { // a new computation is required
360 const Mg3d& mg = *mp.get_mg() ;
361 int nz = mg.get_nzone() ;
362 assert(nz>2) ;
363 int np = mg.get_np(0) ;
364 if (np != 1) {
365 cout << "The polar circumferential radius is only well defined\n"
366 << "with np = 1!" << endl ;
367 abort() ;
368 }
369 int nt = mg.get_nt(0) ;
370 Sym_tensor gam = gamma.cov() ;
371 const Coord& tet = mp.tet ;
372 Scalar rrr(mp) ;
373 rrr.annule_hard() ;
374 rrr.annule(nzet,nz-1) ;
375 double phi = 0 ;
376 for (int j=0; j<nt; j++) {
377 double theta = (+tet)(0, 0, j, 0) ;
378 double erre =
379 mp.val_r(l_surf()(0,j), xi_surf()(0, j), theta, phi) ;
380 for (int lz=0; lz<nzet; lz++) {
381 int nrz = mg.get_nr(lz) ;
382 for (int i=0; i<nrz; i++) {
383 rrr.set_grid_point(lz, 0, j, i) = erre ;
384 }
385 }
386 }
387 rrr.std_spectral_base() ;
388 Scalar drrr = rrr.dsdt() ;
389
390 Scalar fff(mp) ;
391 fff.annule_hard() ;
392 fff.annule(nzet,nz-1) ;
393 for (int j=0; j<nt; j++) {
394 double theta = (+tet)(0, 0, j, 0) ;
395 int ls = l_surf()(0, j) ;
396 double xs = xi_surf()(0, j) ;
397 double grr = gam(1,1).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
398 double grt = gam(1,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
399 double gtt = gam(2,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
400 double rr = mp.val_r(ls, xs, theta, phi) ;
401 double dr = drrr.get_spectral_va().val_point_jk(ls, xs, j, 0) ;
402 for (int i=0; i< mg.get_nr(nzet-1); i++) {
403 fff.set_grid_point(nzet-1, 0, j, i)
404 = sqrt(grr*dr*dr + 2*grt*rr*dr + gtt*rr*rr) ;
405 }
406 }
407 fff.std_spectral_base() ;
408 fff.set_spectral_va().coef() ;
409 p_rp_circ = new double((*fff.get_spectral_va().c_cf)(nzet-1, 0, 0, 0)) ;
410 }
411 return *p_rp_circ ;
412}
413
414 //--------------------------//
415 // Flattening //
416 //--------------------------//
417
418double Star_rot_Dirac::aplat() const {
419
420 return ray_pole() / ray_eq() ;
421
422}
423
424 //--------------------------//
425 // Ellipticity //
426 //--------------------------//
427
429
430 return sqrt(1. - (rp_circ()*rp_circ())/(r_circ()*r_circ())) ;
431
432}
433}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Coord tet
coordinate centered on the grid
Definition map.h:719
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
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
virtual const Connection & connect() const
Returns the connection.
Definition metric.C:301
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
Multi-domain grid.
Definition grilles.h:273
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition grilles.h:485
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Scalar & dsdt() const
Returns of *this .
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
const Scalar & dsdr() const
Returns of *this .
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
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:684
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
virtual double mass_g() const
Gravitational mass.
virtual double ellipt() const
Ellipticity e.
virtual double grv3() const
Error on the virial identity GRV3.
virtual double angu_mom() const
Angular momentum.
double * p_grv3
Error on the virial identity GRV3.
virtual double r_circ() const
Circumferential equatorial radius.
double omega
Rotation angular velocity ([f_unit] )
double * p_tsw
Ratio T/W.
double * p_r_circ
Circumferential equatorial radius.
virtual double tsw() const
Ratio T/W.
Vector j_euler
Momentum density 3-vector with respect to the Eulerian observer.
virtual double rp_circ() const
Circumferential polar radius.
double * p_grv2
Error on the virial identity GRV2.
virtual double mass_b() const
Baryonic mass.
double * p_angu_mom
Angular momentum.
double * p_rp_circ
Circumferential polar radius.
virtual double aplat() const
Flattening r_pole/r_eq.
virtual double grv2() const
Error on the virial identity GRV2.
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Scalar nn
Lapse function N .
Definition star.h:225
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:63
double * p_mass_b
Baryon mass.
Definition star.h:268
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
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:89
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
double ray_eq() const
Coordinate radius at , [r_unit].
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
double * p_mass_g
Gravitational mass.
Definition star.h:269
Metric gamma
3-metric
Definition star.h:235
Scalar press
Fluid pressure.
Definition star.h:194
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Vector beta
Shift vector.
Definition star.h:228
double ray_pole() const
Coordinate radius at [r_unit].
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor handling.
Definition tensor.h:288
void coef() const
Computes the coeffcients of *this.
Tensor field of valence 1.
Definition vector.h:188
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.