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