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
28char 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// --------------------
81namespace 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) ,
116 stress_euler(co.stress_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// -----------------------
126Compobj::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
155void Compobj::del_deriv() const {
156
157 if (p_adm_mass != 0x0) delete p_adm_mass ;
158
160}
161
162
164
165 p_adm_mass = 0x0 ;
166
167}
168
169 //--------------//
170 // Assignment //
171 //--------------//
172
173// Assignment to another Compobj
174// ----------------------------
176
177 assert( &(co.mp) == &mp ) ; // Same mapping
178
179 nn = co.nn ;
180 beta = co.beta ;
181 gamma = co.gamma ;
183 mom_euler = co.mom_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// --------------
196void 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
210void 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
233ostream& operator<<(ostream& ost, const Compobj& co) {
234 co >> ost ;
235 return ost ;
236}
237
238
239ostream& 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
310double 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}
Base class for stationary compact objects (under development).
Definition compobj.h:126
virtual double adm_mass() const
ADM mass (computed as a surface integral at spatial infinity)
Definition compobj.C:310
Sym_tensor kk
Extrinsic curvature tensor
Definition compobj.h:153
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Definition compobj.h:147
virtual void del_deriv() const
Deletes all the derived quantities.
Definition compobj.C:155
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Definition compobj.h:150
Scalar ener_euler
Total energy density E in the Eulerian frame.
Definition compobj.h:144
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition compobj.C:163
void operator=(const Compobj &)
Assignment to another Compobj.
Definition compobj.C:175
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition compobj.C:239
Metric gamma
3-metric
Definition compobj.h:141
Compobj(Map &map_i)
Standard constructor.
Definition compobj.C:82
virtual void sauve(FILE *) const
Save in a file.
Definition compobj.C:196
Scalar nn
Lapse function N .
Definition compobj.h:135
Vector beta
Shift vector .
Definition compobj.h:138
void gyoto_data(const char *file_name) const
Save in a file for GYOTO.
Definition compobj.C:210
double * p_adm_mass
ADM mass.
Definition compobj.h:158
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition compobj.C:290
virtual ~Compobj()
Destructor.
Definition compobj.C:144
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition compobj.h:132
Base class for coordinate mappings.
Definition map.h:670
virtual void sauve(FILE *) const
Save in a file.
Definition map.C:224
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 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 void sauve(FILE *) const
Save in a file.
Definition metric.C:414
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
Definition mg3d.C:371
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
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
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
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition vector.C:807
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
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
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
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
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 set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene prototypes.
Definition app_hor.h:64