LORENE
compobj.C
1 /*
2  * Methods of the class Compobj
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: compobj.C,v 1.10 2016/12/05 16:17:49 j_novak Exp $
32  * $Log: compobj.C,v $
33  * Revision 1.10 2016/12/05 16:17:49 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.9 2014/10/13 08:52:49 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.8 2014/05/16 11:55:18 o_straub
40  * fixed: GYOTO output from compobj & compobj_QI
41  *
42  * Revision 1.7 2014/04/11 17:22:07 o_straub
43  * Risco and Rms output for GYOTO
44  *
45  * Revision 1.6 2013/07/25 19:44:11 o_straub
46  * calculation of the marginally bound radius
47  *
48  * Revision 1.5 2013/04/03 12:10:13 e_gourgoulhon
49  * Added member kk to Compobj; suppressed tkij
50  *
51  * Revision 1.4 2012/12/03 15:27:30 c_some
52  * Small changes
53  *
54  * Revision 1.3 2012/11/22 16:04:51 c_some
55  * Minor modifications
56  *
57  * Revision 1.2 2012/11/20 16:24:09 c_some
58  * Added computation of ADM mass (method mass_q())
59  *
60  * Revision 1.1 2012/11/15 16:20:51 c_some
61  * New class Compobj
62  *
63  *
64  * $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj.C,v 1.10 2016/12/05 16:17:49 j_novak Exp $
65  *
66  */
67 
68 
69 // C headers
70 #include <cassert>
71 #include <cmath>
72 
73 // Lorene headers
74 #include "compobj.h"
75 #include "nbr_spx.h"
76 #include "utilitaires.h"
77 
78  //--------------//
79  // Constructors //
80  //--------------//
81 
82 // Standard constructor
83 // --------------------
84 namespace Lorene {
86  mp(map_i) ,
87  nn(map_i) ,
88  beta(map_i, CON, map_i.get_bvect_spher()) ,
89  gamma(map_i.flat_met_spher()) ,
90  ener_euler(map_i) ,
91  mom_euler(map_i, CON, map_i.get_bvect_spher()) ,
92  stress_euler(map_i, COV, map_i.get_bvect_spher()) ,
93  kk(map_i, COV, map_i.get_bvect_spher())
94 {
95  // Pointers of derived quantities initialized to zero :
96  set_der_0x0() ;
97 
98  // Some initialisations:
99  nn = 1 ;
100  nn.std_spectral_base() ;
101 
102  beta.set_etat_zero() ;
103  ener_euler = 0 ;
106  kk.set_etat_zero() ;
107 
108 }
109 
110 // Copy constructor
111 // --------------------
113  mp(co.mp) ,
114  nn(co.nn) ,
115  beta(co.beta) ,
116  gamma(co.gamma) ,
117  ener_euler(co.ener_euler) ,
118  mom_euler(co.mom_euler) ,
119  stress_euler(co.stress_euler) ,
120  kk(co.kk)
121 {
122  // Pointers of derived quantities initialized to zero :
123  set_der_0x0() ;
124 }
125 
126 
127 // Constructor from a file
128 // -----------------------
129 Compobj::Compobj(Map& map_i, FILE* fich) :
130  mp(map_i) ,
131  nn(map_i, *(map_i.get_mg()), fich) ,
132  beta(map_i, map_i.get_bvect_spher(), fich) ,
133  gamma(map_i, fich) ,
134  ener_euler(map_i, *(map_i.get_mg()), fich) ,
135  mom_euler(map_i, map_i.get_bvect_spher(), fich) ,
136  stress_euler(map_i, map_i.get_bvect_spher(), fich) ,
137  kk(map_i, COV, map_i.get_bvect_spher())
138 {
139  // Pointers of derived quantities initialized to zero :
140  set_der_0x0() ;
141 }
142 
143  //------------//
144  // Destructor //
145  //------------//
146 
148 
149  del_deriv() ;
150 
151 }
152 
153 
154  //----------------------------------//
155  // Management of derived quantities //
156  //----------------------------------//
157 
158 void Compobj::del_deriv() const {
159 
160  if (p_adm_mass != 0x0) delete p_adm_mass ;
161 
163 }
164 
165 
166 void Compobj::set_der_0x0() const {
167 
168  p_adm_mass = 0x0 ;
169 
170 }
171 
172  //--------------//
173  // Assignment //
174  //--------------//
175 
176 // Assignment to another Compobj
177 // ----------------------------
178 void Compobj::operator=(const Compobj& co) {
179 
180  assert( &(co.mp) == &mp ) ; // Same mapping
181 
182  nn = co.nn ;
183  beta = co.beta ;
184  gamma = co.gamma ;
185  ener_euler = co.ener_euler ;
186  mom_euler = co.mom_euler ;
187  stress_euler = co.stress_euler ;
188  kk = co.kk ;
189 
190  del_deriv() ; // Deletes all derived quantities
191 }
192 
193  //--------------//
194  // Outputs //
195  //--------------//
196 
197 // Save in a file
198 // --------------
199 void Compobj::sauve(FILE* fich) const {
200 
201  nn.sauve(fich) ;
202  beta.sauve(fich) ;
203  gamma.sauve(fich) ;
204  ener_euler.sauve(fich) ;
205  mom_euler.sauve(fich) ;
206  stress_euler.sauve(fich) ;
207 
208 }
209 
210 // Save in a file for GYOTO input
211 // ------------------------------
212 
213 void Compobj::gyoto_data(const char* file_name) const {
214 
215  FILE* file_out = fopen(file_name, "w") ;
216  double total_time = 0. ; // for compatibility
217 
218  fwrite_be(&total_time, sizeof(double), 1, file_out) ;
219  mp.get_mg()->sauve(file_out) ;
220  mp.sauve(file_out) ;
221  nn.sauve(file_out) ;
222  beta.sauve(file_out) ;
223  gamma.cov().sauve(file_out) ;
224  gamma.con().sauve(file_out) ;
225  kk.sauve(file_out) ;
226 
227  fclose(file_out) ;
228 
229 
230  cout << "WRITING TO GYOTO FILE - OK: " << endl ;
231 }
232 
233 // Printing
234 // --------
235 
236 ostream& operator<<(ostream& ost, const Compobj& co) {
237  co >> ost ;
238  return ost ;
239 }
240 
241 
242 ostream& Compobj::operator>>(ostream& ost) const {
243 
244  ost << endl << "Compact object (class Compobj) " << endl ;
245  ost << "Mapping : " << mp << endl ;
246  ost << "Central values of various fields : " << endl ;
247  ost << "-------------------------------- " << endl ;
248  ost << " lapse function : N_c = " << nn.val_grid_point(0,0,0,0) << endl ;
249  ost << " metric components gamma_{ij} : " << endl
250  << " ( " << gamma.cov()(1,1).val_grid_point(0,0,0,0) << " "
251  << gamma.cov()(1,2).val_grid_point(0,0,0,0) << " "
252  << gamma.cov()(1,3).val_grid_point(0,0,0,0) << " )" << endl
253  << " ( " << gamma.cov()(2,1).val_grid_point(0,0,0,0) << " "
254  << gamma.cov()(2,2).val_grid_point(0,0,0,0) << " "
255  << gamma.cov()(2,3).val_grid_point(0,0,0,0) << " )" << endl
256  << " ( " << gamma.cov()(3,1).val_grid_point(0,0,0,0) << " "
257  << gamma.cov()(3,2).val_grid_point(0,0,0,0) << " "
258  << gamma.cov()(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
259  ost << " components of the extrinsic curvature K_{ij} : " << endl
260  << " ( " << kk(1,1).val_grid_point(0,0,0,0) << " "
261  << kk(1,2).val_grid_point(0,0,0,0) << " "
262  << kk(1,3).val_grid_point(0,0,0,0) << " )" << endl
263  << " ( " << kk(2,1).val_grid_point(0,0,0,0) << " "
264  << kk(2,2).val_grid_point(0,0,0,0) << " "
265  << kk(2,3).val_grid_point(0,0,0,0) << " )" << endl
266  << " ( " << kk(3,1).val_grid_point(0,0,0,0) << " "
267  << kk(3,2).val_grid_point(0,0,0,0) << " "
268  << kk(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
269  ost << " energy density / Eulerian observer : E_c = " << ener_euler.val_grid_point(0,0,0,0) << endl ;
270  ost << " components of the stress tensor S_{ij} / Eulerian observer : " << endl
271  << " ( " << stress_euler(1,1).val_grid_point(0,0,0,0) << " "
272  << stress_euler(1,2).val_grid_point(0,0,0,0) << " "
273  << stress_euler(1,3).val_grid_point(0,0,0,0) << " )" << endl
274  << " ( " << stress_euler(2,1).val_grid_point(0,0,0,0) << " "
275  << stress_euler(2,2).val_grid_point(0,0,0,0) << " "
276  << stress_euler(2,3).val_grid_point(0,0,0,0) << " )" << endl
277  << " ( " << stress_euler(3,1).val_grid_point(0,0,0,0) << " "
278  << stress_euler(3,2).val_grid_point(0,0,0,0) << " "
279  << stress_euler(3,3).val_grid_point(0,0,0,0) << " )" << endl ;
280 
281 //## ost << endl << "ADM mass : " << adm_mass() << endl ;
282 
283  return ost ;
284 
285 }
286 
287 
288  //-------------------------//
289  // Computational methods //
290  //-------------------------//
291 
292 // Extrinsic curvature
294 
295  cout << "WARNING: Compobj::extrinsic_curvature() NOT TESTED !" << endl ;
296 
297  // Gradient of the shift D_j beta_i
298  Vector cobeta = beta.down(0, gamma) ;
299 
300  Tensor dn = cobeta.derive_cov(gamma) ;
301 
302  kk.set_etat_qcq() ;
303  for (int i=1; i<=3; i++) {
304  for (int j=i; j<=3; j++) {
305  kk.set(i, j) = (dn(i, j) + dn(j, i))/(2*nn) ;
306  }
307  }
308 
309 }
310 
311 
312 // Gravitational mass
313 double Compobj::adm_mass() const {
314 
315  if (p_adm_mass == 0x0) { // a new computation is required
316 
317  const Sym_tensor& gam_dd = gamma.cov() ; // components \gamma_{ij} of the 3-metric
318  Metric_flat ff(mp, *(gam_dd.get_triad())) ;
319 
320  Vector ww = gam_dd.derive_con(ff).trace(1,2).up(0,ff)
321  - gam_dd.trace(ff).derive_con(ff) ;
322 
323  p_adm_mass = new double( ww.flux(__infinity, ff) / (16.* M_PI) ) ;
324  }
325 
326  return *p_adm_mass ;
327 
328 }
329 
330 
331 }
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:375
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Definition: compobj.h:150
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition: vector.C:813
virtual void sauve(FILE *) const
Save in a file.
Definition: map.C:227
virtual double adm_mass() const
ADM mass (computed as a surface integral at spatial infinity)
Definition: compobj.C:313
Lorene prototypes.
Definition: app_hor.h:67
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Flat metric for tensor calculation.
Definition: metric.h:261
virtual ~Compobj()
Destructor.
Definition: compobj.C:147
double * p_adm_mass
ADM mass.
Definition: compobj.h:161
Base class for coordinate mappings.
Definition: map.h:688
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
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
Sym_tensor kk
Extrinsic curvature tensor .
Definition: compobj.h:156
Tensor field of valence 1.
Definition: vector.h:188
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition: compobj.C:293
void gyoto_data(const char *file_name) const
Save in a file for GYOTO.
Definition: compobj.C:213
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 sauve(FILE *) const
Save in a file.
Definition: metric.C:417
Base class for stationary compact objects (under development).
Definition: compobj.h:129
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:692
void operator=(const Compobj &)
Assignment to another Compobj.
Definition: compobj.C:178
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
Definition: mg3d.C:500
Vector beta
Shift vector .
Definition: compobj.h:141
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj.C:158
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:73
Tensor handling.
Definition: tensor.h:294
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: compobj.C:242
Compobj(Map &map_i)
Standard constructor.
Definition: compobj.C:85
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
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
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Scalar ener_euler
Total energy density E in the Eulerian frame.
Definition: compobj.h:147
Scalar nn
Lapse function N .
Definition: compobj.h:138
Metric gamma
3-metric
Definition: compobj.h:144
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
virtual void sauve(FILE *) const
Save in a file.
Definition: compobj.C:199
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: compobj.C:166
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:135
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Definition: compobj.h:153
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.