LORENE
excised_slice.C
1/*
2 * Methods of class Excised_slice
3 *
4 * (see file excised_slice.h for documentation)
5 * Disclaimer: the class Isol_hole() is redundant with this class under a set of parameters;
6 * therefore it has to go at some point.
7 */
8
9
10/*
11 * Copyright (c) 2010 Nicolas Vasset
12
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or
18 * (at your option) any later version.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31
32// Headers C
33#include "math.h"
34
35// Headers Lorene
36#include "excised_slice.h"
37#include "spheroid.h"
38#include "excision_surf.h"
39#include "excision_hor.h"
40#include "utilitaires.h"
41#include "param.h"
42#include "unites.h"
43#include "proto.h"
44
45namespace Lorene {
46
47 // Fundamental constants and units
48 // -------------------------------
49 using namespace Unites ;
50
51 //--------------//
52 // Constructors //
53 //--------------//
54
55
56// Standard constructor
57// --------------------
58Excised_slice::Excised_slice (const Map& mpi, int hor_type, int scheme_type)
59 : mp(mpi),
60 type_hor(hor_type),
61 field_set(scheme_type),
62 lapse(mpi),
63 conf_fact(mpi),
64 shift(mpi, CON, mpi.get_bvect_spher()),
65 hij(mpi, CON, mpi.get_bvect_spher()),
66 hatA(mpi, CON, mpi.get_bvect_spher()),
67 Xx(mpi, CON, mpi.get_bvect_spher()){
68
69 // Pointers of derived quantities initialized to zero :
70 set_der_0x0() ;
71
72 // Initializing primary quantities.
73
74 lapse = 1. ;
75 conf_fact = 1. ;
80}
81
82// Copy constructor
83// ----------------
85 : mp(ih.mp),
86 type_hor(ih.type_hor),
87 field_set(ih.field_set),
88 lapse(ih.lapse),
89 conf_fact(ih.conf_fact),
90 shift(ih.shift),
91 hij(ih.hij),
92 hatA(ih.hatA),
93 Xx(ih.Xx){
94
95 set_der_0x0() ;
96
97}
98
99// Constructor from a file
100// -----------------------
101Excised_slice::Excised_slice(const Map& mpi, int hor_type, int scheme_type, FILE* fich)
102 : mp(mpi),
103 type_hor(hor_type),
104 field_set(scheme_type),
105 lapse(mpi,*(mpi.get_mg()), fich),
106 conf_fact(mpi,*(mpi.get_mg()), fich),
107 shift(mpi, mpi.get_bvect_spher(), fich),
108 hij(mpi, mpi.get_bvect_spher(), fich),
109 hatA(mpi, mpi.get_bvect_spher(), fich),
110 Xx(mpi, mpi.get_bvect_spher(), fich){
111
112
113
114 // Pointers of derived quantities initialized to zero
115 // --------------------------------------------------
116 set_der_0x0() ;
117
118}
119
120 //------------//
121 // Destructor //
122 //------------//
123
129
130
131 //----------------------------------//
132 // Management of derived quantities //
133 //----------------------------------//
134
136 if (p_hor != 0x0) delete p_hor ;
137 if (p_adm_mass != 0x0) delete p_adm_mass ;
138 if (p_komar_angmom != 0x0) delete p_komar_angmom ;
139 if (p_virial_residue != 0x0) delete p_virial_residue ;
141}
142
143
145 p_hor = 0x0;
146 p_adm_mass = 0x0 ;
147 p_komar_angmom = 0x0 ;
148 p_virial_residue = 0x0 ;
149}
150
151
152
153 //--------------//
154 // Assignment //
155 //--------------//
156
157// Assignment to another Excised_slice
158// ----------------------------
160
161 assert( &(ih.mp) == &mp ) ; // Same mapping
162 type_hor = ih.type_hor ;
163 field_set = ih.field_set ;
164 lapse = ih.lapse ;
165 conf_fact = ih.conf_fact ;
166 shift = ih.shift ;
167 hij = ih.hij ;
168 hatA = ih.hatA ;
169 Xx = ih.Xx;
170
171 del_deriv() ; // Deletes all derived quantities
172
173}
174
175 //--------------//
176 // Outputs //
177 //--------------//
178
179// Save in a file
180// --------------
181void Excised_slice::sauve(FILE* fich) const {
182
183 lapse.sauve(fich) ;
184 conf_fact.sauve(fich) ;
185 shift.sauve(fich);
186 hij.sauve(fich);
187 hatA.sauve(fich);
188 Xx.sauve(fich);}
189
190
191// Prints out maximal errors in Einstein equations for the obtained metric fields
192
194
195const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
196 const Metric_flat& mets = (*map).flat_met_spher() ;
197
198 Sym_tensor gamtcon = mets.con() + hij;
199 Metric gamt(gamtcon);
200 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
201 gamcon.std_spectral_base();
202 Metric gam(gamcon);
203 Sym_tensor k_uu = hatA/(pow(conf_fact,10)) ;
204 k_uu.std_spectral_base();
205 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
206 Sym_tensor k_dd = k_uu.up_down(gam);
207
208 Scalar TrK3 = k_uu.trace(gam);
209
210 // Hamiltonian constraint
211 //-----------------------
212 Scalar ham_constr = gam.ricci_scal() ;
213 ham_constr.dec_dzpuis(3) ; // To check
214 ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
215 maxabs(ham_constr, "Hamiltonian constraint: ") ;
216
217 // Momentum constraint
218 //-------------------
219 Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
220 mom_constr.dec_dzpuis(2) ; // To check
221 maxabs(mom_constr, "Momentum constraint: ") ;
222
223 // Evolution equations
224 //--------------------
225 Sym_tensor evol_eq = lapse*gam.ricci()
226 - lapse.derive_cov(gam).derive_cov(gam);
227 evol_eq.dec_dzpuis() ;
228 evol_eq += k_dd.derive_lie(shift) ;
229 evol_eq.dec_dzpuis(2) ; // To check
230 evol_eq += lapse*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
231 maxabs(evol_eq, "Evolution equations: ") ;
232
233 return;
234}
235
236
237
238
239
240 //----------------------------//
241 // Accessors/ Derived data //
242 //----------------------------//
243
244// Computation of the Spheroid corresponding to the black hole excision surface.
245
247 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
248 const Mg3d* mgrid = (*map).get_mg();
249
250 // Construct angular grid for h(theta,phi)
251 const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
252
253 const Coord& rr = (*map).r;
254 Scalar rrr (*map) ;
255 rrr = rr ;
256 rrr.std_spectral_base();
257 assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9); // For now the code handles only horizons at r=1, corresponding to the first shell inner boundary. This test assures this is the case with our mapping.
258
259 // Angular mapping defined for one domain (argument of spheroid Class)
260 //--------------------------------------------------------------------
261
262 double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ;
263 const Map_af map_2(*g_angu, r_limits2);
264
265 //Full 3-metric and extrrinsic curvature
266 const Metric_flat& mets = (*map).flat_met_spher() ;
267
268 Sym_tensor gamtcon = mets.con() + hij;
269 Metric gamt(gamtcon);
270 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
271 Metric gam(gamcon);
272
273
274 Sym_tensor kuu = hatA/pow(conf_fact,10) ;
275 kuu.std_spectral_base();
276 Sym_tensor kdd = kuu.up_down(gam);
277
278 //---------------------------------------------------------
279 // Construction of the spheroid associated with those data
280 //--------------------------------------------------------
281 double hor_posd = rrr.val_grid_point(1,0,0,0);
282 Scalar hor_pos(map_2); hor_pos = hor_posd; hor_pos.std_spectral_base();
283 Spheroid hor_loc(hor_pos, gam, kdd);
284 return hor_loc;
285}
286
287// Computation of the ADM mass of the BH spacetime
289 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
290 const Metric_flat& mets = (*map).flat_met_spher() ;
291
292 Sym_tensor gamtcon = mets.con() + hij;
293 Metric gamt(gamtcon);
294 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
295 Metric gam(gamcon);
296
297 Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
298 detgam.std_spectral_base();
299 Vector corr_adm = - (0.125*contract(gamt.cov().derive_con(mets),1,2));
300 Scalar admintegr = conf_fact.dsdr() + corr_adm(1);
301
302 double M_ADM = - (1/(2.*3.1415927*ggrav))*(*map).integrale_surface_infini(admintegr*detgam);
303 return M_ADM;
304}
305
306
307// Computation of the Komar angular momentum w.r.t. assumed rotational symmetry
309const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
310 const Metric_flat& mets = (*map).flat_met_spher() ;
311
312 Sym_tensor gamtcon = mets.con() + hij;
313 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
314 gamcon.std_spectral_base();
315 Metric gam(gamcon);
316 Sym_tensor k_uu = hatA/(pow(conf_fact,10));
317 k_uu.std_spectral_base();
318 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
319 Sym_tensor k_dd = k_uu.up_down(gam);
320
321 Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
322 detgam.std_spectral_base();
323 Scalar contraction = k_dd(1,3); contraction.mult_r_dzpuis(2); contraction.mult_sint();
324 double angu_komar = - (1/(8.*3.1415927*ggrav))*(*map).integrale_surface_infini(detgam*contraction);
325
326 return angu_komar;
327}
328
329
330// Computation of the Virial residual, as rescaled difference at infinity
331// between the ADM mass and the Komar integral associated to the mass
332
334 const Mg3d* mgrid = mp.get_mg();
335 const int nz = (*mgrid).get_nzone(); // Number of domains
336 Valeur** devel_psi (conf_fact.asymptot(1)) ;
337 Valeur** devel_n (lapse.asymptot(1)) ;
338
339
340 double erreur = (2*(*devel_psi[1])(nz-1, 0, 0, 0)
341 + (*devel_n[1])(nz-1, 0, 0, 0))/fabs ((*devel_n[1])(nz-1, 0, 0, 0)) ;
342
343 return erreur;
344}
345
346
347
348
349
350}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Class to compute single black hole spacetime excised slices.
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: rescaled same way as Cordero et al.
double * p_adm_mass
Computation of the ADM mass of the BH spacetime.
double komar_angmom()
Computation of the Komar angular momentum w.r.t.
int type_hor
Chosen horizon type.
void Einstein_errors()
Prints out errors in Einstein equations for the data obtained.
double * p_komar_angmom
Computation of the Komar angular momentum w.r.t.
virtual void del_deriv() const
Deletes all the derived quantities.
Spheroid * p_hor
Computation of the spheroid associated with the black hole horizon.
double virial_residue()
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Scalar lapse
Lapse function.
double * p_virial_residue
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
void operator=(const Excised_slice &)
Assignment to another Excised_slice.
virtual ~Excised_slice()
Destructor.
const Map & mp
Mapping associated with the slice.
int field_set
Chose field set type.
Vector shift
Shift vector.
virtual void sauve(FILE *) const
Save in a file.
Excised_slice(const Map &mp_i, int hor_type, int scheme_type)
Standard constructor.
double adm_mass()
Computation of the ADM mass of the BH spacetime.
Vector Xx
Longitudinal part of the extrinsic curvature.
Spheroid hor()
Spheroid at the excised boundary associated with the black hole MOTS on the slice.
Affine radial mapping.
Definition map.h:2027
Base class for coordinate mappings.
Definition map.h:670
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.
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition metric.C:338
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition metric.C:350
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition mg3d.C:494
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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
const Scalar & dsdr() const
Returns of *this .
void mult_sint()
Multiplication by .
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Definition spheroid.h:84
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:349
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tensor field of valence 1.
Definition vector.h:188
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
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
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
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...
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
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.