LORENE
star_QI.C
1 /*
2  * Methods of the class Star_QI
3  *
4  * (see file compobj.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2012 Claire Some, Eric Gourgoulhon
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_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $" ;
29 
30 /*
31  * $Id: star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
32  * $Log: star_QI.C,v $
33  * Revision 1.5 2014/10/13 08:52:49 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.4 2012/12/03 15:27:30 c_some
37  * Small changes
38  *
39  * Revision 1.3 2012/11/22 16:04:51 c_some
40  * Minor modifications
41  *
42  * Revision 1.2 2012/11/21 14:55:27 c_some
43  * Added methods fait_shift and fait_nphi
44  *
45  * Revision 1.1 2012/11/20 16:29:50 c_some
46  * New class Star_QI
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
50  *
51  */
52 
53 
54 // C headers
55 #include <cassert>
56 
57 // Lorene headers
58 #include "compobj.h"
59 #include "unites.h"
60 
61  //--------------//
62  // Constructors //
63  //--------------//
64 
65 // Standard constructor
66 // --------------------
67 namespace Lorene {
69  Compobj_QI(mpi) ,
70  logn(mpi),
71  tnphi(mpi),
72  nuf(mpi),
73  nuq(mpi),
74  dzeta(mpi),
75  tggg(mpi),
76  w_shift(mpi, CON, mp.get_bvect_cart()),
77  khi_shift(mpi),
78  ssjm1_nuf(mpi),
79  ssjm1_nuq(mpi),
80  ssjm1_dzeta(mpi),
81  ssjm1_tggg(mpi),
82  ssjm1_khi(mpi),
83  ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
84 {
85  // Pointers of derived quantities initialized to zero :
86  set_der_0x0() ;
87 
88  // Initialization to a flat metric :
89  logn = 0 ;
90  tnphi = 0 ;
91  nuf = 0 ;
92  nuq = 0 ;
93  dzeta = 0 ;
94  tggg = 0 ;
95 
97  khi_shift = 0 ;
98 
99  beta.set_etat_zero() ;
100 
101  ssjm1_nuf = 0 ;
102  ssjm1_nuq = 0 ;
103  ssjm1_dzeta = 0 ;
104  ssjm1_tggg = 0 ;
105  ssjm1_khi = 0 ;
106 
108 
109 }
110 
111 // Copy constructor
112 // --------------------
114  Compobj_QI(st),
115  logn(st.logn),
116  tnphi(st.tnphi),
117  nuf(st.nuf),
118  nuq(st.nuq),
119  dzeta(st.dzeta),
120  tggg(st.tggg),
121  w_shift(st.w_shift),
122  khi_shift(st.khi_shift),
123  ssjm1_nuf(st.ssjm1_nuf),
124  ssjm1_nuq(st.ssjm1_nuq),
126  ssjm1_tggg(st.ssjm1_tggg),
127  ssjm1_khi(st.ssjm1_khi),
129 {
130  // Pointers of derived quantities initialized to zero :
131  set_der_0x0() ;
132 }
133 
134 
135 // Constructor from a file
136 // -----------------------
137 Star_QI::Star_QI(Map& mpi, FILE* fich) :
138  Compobj_QI(mpi) ,
139  logn(mpi),
140  tnphi(mpi),
141  nuf(mpi),
142  nuq(mpi),
143  dzeta(mpi),
144  tggg(mpi),
145  w_shift(mpi, CON, mp.get_bvect_cart()),
146  khi_shift(mpi),
147  ssjm1_nuf(mpi),
148  ssjm1_nuq(mpi),
149  ssjm1_dzeta(mpi),
150  ssjm1_tggg(mpi),
151  ssjm1_khi(mpi),
152  ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
153 {
154  // Pointers of derived quantities initialized to zero :
155  set_der_0x0() ;
156 
157  // Read of the saved fields:
158  // ------------------------
159 
160  Scalar nuf_file(mp, *(mp.get_mg()), fich) ;
161  nuf = nuf_file ;
162 
163  Scalar nuq_file(mp, *(mp.get_mg()), fich) ;
164  nuq = nuq_file ;
165 
166  logn = nuf + nuq ; //## to be checked !
167 
168  Scalar dzeta_file(mp, *(mp.get_mg()), fich) ;
169  dzeta = dzeta_file ;
170 
171  Scalar tggg_file(mp, *(mp.get_mg()), fich) ;
172  tggg = tggg_file ;
173 
174  Vector w_shift_file(mp, mp.get_bvect_cart(), fich) ;
175  w_shift = w_shift_file ;
176 
177  Scalar khi_shift_file(mp, *(mp.get_mg()), fich) ;
178  khi_shift = khi_shift_file ;
179 
180  fait_shift() ; // constructs shift from w_shift and khi_shift
181  fait_nphi() ; // constructs N^phi from (N^x,N^y,N^z)
182 
183  Scalar ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ;
184  ssjm1_nuf = ssjm1_nuf_file ;
185 
186  Scalar ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ;
187  ssjm1_nuq = ssjm1_nuq_file ;
188 
189  Scalar ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ;
190  ssjm1_dzeta = ssjm1_dzeta_file ;
191 
192  Scalar ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ;
193  ssjm1_tggg = ssjm1_tggg_file ;
194 
195  Scalar ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
196  ssjm1_khi = ssjm1_khi_file ;
197 
198  Vector ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
199  ssjm1_wshift = ssjm1_wshift_file ;
200 
201 
202  // Initialization of N, A^2, B^2, gamma_ij, tkij and ak_car
203  update_metric() ;
204 
205 }
206 
207  //------------//
208  // Destructor //
209  //------------//
210 
212 
213  del_deriv() ;
214 
215 }
216 
217 
218  //----------------------------------//
219  // Management of derived quantities //
220  //----------------------------------//
221 
222 void Star_QI::del_deriv() const {
223 
225 
226  if (p_grv2 != 0x0) delete p_grv2 ;
227  if (p_grv3 != 0x0) delete p_grv3 ;
228  if (p_mom_quad != 0x0) delete p_mom_quad ;
229  if (p_mass_g != 0x0) delete p_mass_g ;
230 
232 }
233 
234 
235 void Star_QI::set_der_0x0() const {
236 
237  p_grv2 = 0x0 ;
238  p_grv3 = 0x0 ;
239  p_mom_quad = 0x0 ;
240  p_mass_g = 0x0 ;
241 
242 }
243 
244  //--------------//
245  // Assignment //
246  //--------------//
247 
248 // Assignment to another Star_QI
249 // --------------------------------
250 void Star_QI::operator=(const Star_QI& st) {
251 
252  // Assignment of quantities common to all the derived classes of Compobj_QI
253  Compobj_QI::operator=(st) ;
254 
255  logn = st.logn ;
256  tnphi = st.tnphi ;
257  nuf = st.nuf ;
258  nuq = st.nuq ;
259  dzeta = st.dzeta ;
260  tggg = st.tggg ;
261  w_shift = st.w_shift ;
262  khi_shift = st.khi_shift ;
263  ssjm1_nuf = st.ssjm1_nuf ;
264  ssjm1_nuq = st.ssjm1_nuq ;
265  ssjm1_dzeta = st.ssjm1_dzeta ;
266  ssjm1_tggg = st.ssjm1_tggg ;
267  ssjm1_khi = st.ssjm1_khi ;
268  ssjm1_wshift = st.ssjm1_wshift ;
269 
270  del_deriv() ; // Deletes all derived quantities
271 }
272 
273  //--------------//
274  // Outputs //
275  //--------------//
276 
277 // Save in a file
278 // --------------
279 void Star_QI::sauve(FILE* fich) const {
280 
281  nuf.sauve(fich) ;
282  nuq.sauve(fich) ;
283  dzeta.sauve(fich) ;
284  tggg.sauve(fich) ;
285  w_shift.sauve(fich) ;
286  khi_shift.sauve(fich) ;
287 
288  ssjm1_nuf.sauve(fich) ;
289  ssjm1_nuq.sauve(fich) ;
290  ssjm1_dzeta.sauve(fich) ;
291  ssjm1_tggg.sauve(fich) ;
292  ssjm1_khi.sauve(fich) ;
293  ssjm1_wshift.sauve(fich) ;
294 
295 }
296 
297 // Printing
298 // --------
299 
300 ostream& Star_QI::operator>>(ostream& ost) const {
301 
302  using namespace Unites ;
303 
304  Compobj_QI::operator>>(ost) ;
305 
306  ost << endl << "Axisymmetric stationary compact star in quasi-isotropic coordinates (class Star_QI) " << endl ;
307 
308  ost << "Central values of various fields : " << endl ;
309  ost << "-------------------------------- " << endl ;
310  ost << " ln(N) : " << logn.val_grid_point(0,0,0,0) << endl ;
311  ost << " nuf : " << nuf.val_grid_point(0,0,0,0) << endl ;
312  ost << " nuq : " << nuq.val_grid_point(0,0,0,0) << endl ;
313  ost << " zeta = ln(AN): " << dzeta.val_grid_point(0,0,0,0) << endl << endl ;
314 
315  ost << "Error on the virial identity GRV2 : " << grv2() << endl ;
316  ost << "Error on the virial identity GRV3 : " << grv3(&ost) << endl ;
317 
318  double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
319  / double(1.e38) ) ;
320  ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
321  << endl ;
322  ost << "c^4 Q / (G^2 M^3) : "
323  << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
324  << endl ;
325 
326  ost << "Total angular momentum J : "
327  << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
328  << endl ;
329  ost << "c J / (G M^2) : "
330  << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;
331 
332  return ost ;
333 
334 }
335 
336  //-------------------------//
337  // Computational methods //
338  //-------------------------//
339 
341 
342  Vector d_khi = khi_shift.derive_con( mp.flat_met_cart() ) ;
343 
344  d_khi.dec_dzpuis(2) ; // divide by r^2 in the external compactified domain
345 
346  // x_k dW^k/dx_i
347  Scalar xx(mp) ;
348  Scalar yy(mp) ;
349  Scalar zz(mp) ;
350  Scalar sintcosp(mp) ;
351  Scalar sintsinp(mp) ;
352  Scalar cost(mp) ;
353  xx = mp.x ;
354  yy = mp.y ;
355  zz = mp.z ;
356  sintcosp = mp.sint * mp.cosp ;
357  sintsinp = mp.sint * mp.sinp ;
358  cost = mp.cost ;
359 
360  int nz = mp.get_mg()->get_nzone() ;
361  Vector xk(mp, COV, mp.get_bvect_cart()) ;
362  xk.set(1) = xx ;
363  xk.set(1).set_domain(nz-1) = sintcosp.domain(nz-1) ;
364  xk.set(1).set_dzpuis(-1) ;
365  xk.set(2) = yy ;
366  xk.set(2).set_domain(nz-1) = sintsinp.domain(nz-1) ;
367  xk.set(2).set_dzpuis(-1) ;
368  xk.set(3) = zz ;
369  xk.set(3).set_domain(nz-1) = cost.domain(nz-1) ;
370  xk.set(3).set_dzpuis(-1) ;
371  xk.std_spectral_base() ;
372 
374 
375  Vector x_d_w = contract(xk, 0, d_w, 0) ;
376  x_d_w.dec_dzpuis() ;
377 
378  double lambda = double(1) / double(3) ;
379 
380  beta = - (lambda+2)/2./(lambda+1) * w_shift
381  + (lambda/2./(lambda+1)) * (d_khi + x_d_w) ;
382 
383 }
384 
385 
387 
388  if ( (beta(1).get_etat() == ETATZERO) && (beta(2).get_etat() == ETATZERO) ) {
389  tnphi = 0 ;
390  nphi = 0 ;
391  return ;
392  }
393 
394  // Computation of tnphi
395  // --------------------
396 
397  mp.comp_p_from_cartesian( -beta(1), -beta(2), tnphi ) ;
398 
399  // Computation of nphi
400  // -------------------
401 
402  nphi = tnphi ;
403  nphi.div_rsint() ;
404 
405 }
406 
407 
408 // Updates the 3-metric and the shift
409 
411 
412  // Lapse function N
413  // ----------------
414 
415  nn = exp( logn ) ;
416 
417  nn.std_spectral_base() ; // set the bases for spectral expansions
418 
419 
420  // Metric factor A^2
421  // -----------------
422 
423  a_car = exp( 2*( dzeta - logn ) ) ;
424 
425  a_car.std_spectral_base() ; // set the bases for spectral expansions
426 
427  // Metric factor B
428  // ---------------
429 
430  Scalar tmp = tggg ;
431  tmp.div_rsint() ; //... Division of tG by r sin(theta)
432 
433  bbb = (1 + tmp) / nn ;
434 
435  bbb.std_spectral_base() ; // set the bases for spectral expansions
436 
437  b_car = bbb * bbb ;
438 
439 
440  Compobj_QI::update_metric() ; // updates gamma_{ij}
441 
442 
443  // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
444  // -------------------------------------------------
445 
446  //## extrinsic_curvature() ; // should be done by Compobj_QI::update_metric()
447 
448 
449  // The derived quantities are no longer up to date :
450  // -----------------------------------------------
451 
452  del_deriv() ;
453 
454 }
455 
456 
457 
458 }
Base class for axisymmetric stationary compact objects in Quasi-Isotropic coordinates (under developm...
Definition: compobj.h:280
Scalar logn
Logarithm of the lapse N .
Definition: compobj.h:495
Scalar ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta ...
Definition: compobj.h:556
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: compobj.h:539
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Scalar tnphi
Component of the shift vector.
Definition: compobj.h:500
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition: scalar.h:625
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
virtual void comp_p_from_cartesian(const Scalar &v_x, const Scalar &v_y, Scalar &v_p) const =0
Computes the Spherical component (with respect to bvect_spher ) of a vector given by its cartesian c...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
virtual double mass_g() const
Gravitational mass.
Base class for coordinate mappings.
Definition: map.h:670
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: compobj.h:510
virtual double mom_quad() const
Quadrupole moment.
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
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: compobj.h:551
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: compobj.h:569
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
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Scalar a_car
Square of the metric factor A.
Definition: compobj.h:287
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:331
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
double * p_grv2
Error on the virial identity GRV2.
Definition: compobj.h:585
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: compobj.h:505
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
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
Base class for axisymmetric stationary compact stars in Quasi-Isotropic coordinates (under developmen...
Definition: compobj.h:487
Scalar nphi
Metric coefficient .
Definition: compobj.h:296
Scalar b_car
Square of the metric factor B.
Definition: compobj.h:293
Coord sint
Definition: map.h:721
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:686
Scalar tggg
Metric potential .
Definition: compobj.h:516
Scalar bbb
Metric factor B.
Definition: compobj.h:290
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata&#39;s prescription [Prog.
Definition: star_QI.C:340
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg ...
Definition: compobj.h:561
Vector beta
Shift vector .
Definition: compobj.h:138
Star_QI(Map &mp_i)
Standard constructor.
Definition: star_QI.C:68
void operator=(const Compobj_QI &)
Assignment to another Compobj_QI.
Definition: compobj_QI.C:195
Coord sinp
Definition: map.h:723
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Definition: compobj.h:578
Scalar dzeta
Metric potential .
Definition: compobj.h:513
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Tensor handling.
Definition: tensor.h:288
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: star_QI.C:235
double * p_grv3
Error on the virial identity GRV3.
Definition: compobj.h:586
double * p_mass_g
Gravitational mass (ADM mass as a volume integral)
Definition: compobj.h:588
double * p_mom_quad
Quadrupole moment.
Definition: compobj.h:587
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Coord y
y coordinate centered on the grid
Definition: map.h:727
virtual void update_metric()
Updates the 3-metric from A and B and the shift vector from .
Definition: compobj_QI.C:302
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: compobj_QI.C:264
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: compobj.h:545
Coord cosp
Definition: map.h:724
Scalar nn
Lapse function N .
Definition: compobj.h:135
Coord x
x coordinate centered on the grid
Definition: map.h:726
virtual double grv2() const
Error on the virial identity GRV2.
void div_rsint()
Division by everywhere; dzpuis is not changed.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
virtual double angu_mom() const
Angular momentum.
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_QI.C:222
Coord z
z coordinate centered on the grid
Definition: map.h:728
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj_QI.C:163
void operator=(const Star_QI &)
Assignment to another Star_QI.
Definition: star_QI.C:250
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift ...
Definition: star_QI.C:386
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: star_QI.C:300
void update_metric()
Computes metric coefficients from known potentials.
Definition: star_QI.C:410
Vector w_shift
Vector used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: compobj.h:529
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:132
virtual void sauve(FILE *) const
Save in a file.
Definition: star_QI.C:279
virtual ~Star_QI()
Destructor.
Definition: star_QI.C:211
Coord cost
Definition: map.h:722