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