LORENE
isol_hole.C
1/*
2 * Methods of class Isol_hole
3 *
4 * (see file isol_hole.h for documentation)
5 */
6
7/*
8 * Copyright (c) 2009 Nicolas Vasset
9
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
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// Headers Lorene
29#include "isol_hole.h"
30#include "spheroid.h"
31#include "excision_surf.h"
32#include "excision_hor.h"
33#include "utilitaires.h"
34#include "param.h"
35#include "unites.h"
36#include "proto.h"
37
38namespace Lorene {
39
40 // Fundamental constants and units
41 // -------------------------------
42 using namespace Unites ;
43
44 //--------------//
45 // Constructors //
46 //--------------//
47
48
49// Standard constructor
50// --------------------
51Isol_hole::Isol_hole (const Map& mpi, double Omega_i, bool NorKappa_i, Scalar NoK_i, bool isCF_i)
52 : mp(mpi),
53 Omega(Omega_i),
54 NorKappa(NorKappa_i),
55 boundNoK(NoK_i),
56 isCF (isCF_i),
57 lapse(mpi),
58 conf_fact(mpi),
59 shift(mpi, CON, mpi.get_bvect_spher()),
60 hij(mpi, CON, mpi.get_bvect_spher()),
61 hatA(mpi, CON, mpi.get_bvect_spher()){
62
63
64
65
66
67 assert (boundNoK.get_mp() == mpi); // Check if this is not too strong a condition
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. ;
79}
80
81// Copy constructor
82// ----------------
84 : mp(ih.mp),
85 Omega(ih.Omega),
86 NorKappa(ih.NorKappa),
87 boundNoK(ih.boundNoK),
88 isCF (ih.isCF),
89 lapse(ih.lapse),
90 conf_fact(ih.conf_fact),
91 shift(ih.shift),
92 hij(ih.hij),
93 hatA(ih.hatA){
94
95 set_der_0x0() ;
96
97}
98
99// Constructor from a file
100// -----------------------
101Isol_hole::Isol_hole(const Map& mpi, double Omega_i, bool NorKappa_i, Scalar NoK_i, bool isCF_i, FILE* fich)
102 : mp(mpi),
103 Omega(Omega_i),
104 NorKappa(NorKappa_i),
105 boundNoK(NoK_i),
106 isCF(isCF_i),
107 lapse(mpi,*(mpi.get_mg()), fich),
108 conf_fact(mpi,*(mpi.get_mg()), fich),
109 shift(mpi, mpi.get_bvect_spher(), fich),
110 hij(mpi, mpi.get_bvect_spher(), fich),
111 hatA(mpi, mpi.get_bvect_spher(), fich){
112
113
114
115 // Pointers of derived quantities initialized to zero
116 // --------------------------------------------------
117 set_der_0x0() ;
118
119}
120
121 //------------//
122 // Destructor //
123 //------------//
124
126
127 del_deriv() ;
128
129}
130
131
132 //----------------------------------//
133 // Management of derived quantities //
134 //----------------------------------//
135
137 if (p_hor != 0x0) delete p_hor ;
138 if (p_adm_mass != 0x0) delete p_adm_mass ;
139 if (p_komar_angmom != 0x0) delete p_komar_angmom ;
140 if (p_virial_residue != 0x0) delete p_virial_residue ;
142}
143
144
146 p_hor = 0x0;
147 p_adm_mass = 0x0 ;
148 p_komar_angmom = 0x0 ;
149 p_virial_residue = 0x0 ;
150}
151
152
153
154 //--------------//
155 // Assignment //
156 //--------------//
157
158// Assignment to another Isol_hole
159// ----------------------------
161
162 assert( &(ih.mp) == &mp ) ; // Same mapping
163 Omega = ih.Omega;
164 NorKappa = ih.NorKappa ;
165 boundNoK = ih.boundNoK ;
166 isCF = ih.isCF ;
167 lapse = ih.lapse ;
168 conf_fact = ih.conf_fact ;
169 shift = ih.shift ;
170 hij = ih.hij ;
171 hatA = ih.hatA ;
172
173 del_deriv() ; // Deletes all derived quantities
174
175}
176
177 //--------------//
178 // Outputs //
179 //--------------//
180
181// Save in a file
182// --------------
183void Isol_hole::sauve(FILE* fich) const {
184
185 lapse.sauve(fich) ;
186 conf_fact.sauve(fich) ;
187 shift.sauve(fich);
188 hij.sauve(fich);
189 hatA.sauve(fich);}
190
191
192// Prints out maximal errors in Einstein equations for the obtained metric fields
193
195
196const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
197 const Metric_flat& mets = (*map).flat_met_spher() ;
198
199 Sym_tensor gamtcon = mets.con() + hij;
200 Metric gamt(gamtcon);
201 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
202 gamcon.std_spectral_base();
203 Metric gam(gamcon);
204 Sym_tensor k_uu = hatA/(pow(conf_fact,10)) ;
205 k_uu.std_spectral_base();
206 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
207 Sym_tensor k_dd = k_uu.up_down(gam);
208
209 Scalar TrK3 = k_uu.trace(gam);
210
211 // Hamiltonian constraint
212 //-----------------------
213 Scalar ham_constr = gam.ricci_scal() ;
214 ham_constr.dec_dzpuis(3) ; // To check
215 ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
216 maxabs(ham_constr, "Hamiltonian constraint: ") ;
217
218 // Momentum constraint
219 //-------------------
220 Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
221 mom_constr.dec_dzpuis(2) ; // To check
222 maxabs(mom_constr, "Momentum constraint: ") ;
223
224 // Evolution equations
225 //--------------------
226 Sym_tensor evol_eq = lapse*gam.ricci()
227 - lapse.derive_cov(gam).derive_cov(gam);
228 evol_eq.dec_dzpuis() ;
229 evol_eq += k_dd.derive_lie(shift) ;
230 evol_eq.dec_dzpuis(2) ; // To check
231 evol_eq += lapse*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
232 maxabs(evol_eq, "Evolution equations: ") ;
233
234 return;
235}
236
237
238
239
240
241 //----------------------------//
242 // Accessors/ Derived data //
243 //----------------------------//
244
245// Computation of the Spheroid corresponding to the black hole horizon
246
248 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
249 const Mg3d* mgrid = (*map).get_mg();
250
251 // Construct angular grid for h(theta,phi)
252 const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
253
254 const Coord& rr = (*map).r;
255 Scalar rrr (*map) ;
256 rrr = rr ;
257 rrr.std_spectral_base();
258 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.
259
260 // Angular mapping defined for one domain (argument of spheroid Class)
261 //--------------------------------------------------------------------
262
263 double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ;
264 const Map_af map_2(*g_angu, r_limits2);
265
266 //Full 3-metric and extrrinsic curvature
267 const Metric_flat& mets = (*map).flat_met_spher() ;
268
269 Sym_tensor gamtcon = mets.con() + hij;
270 Metric gamt(gamtcon);
271 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
272 Metric gam(gamcon);
273
274
275 Sym_tensor kuu = hatA/pow(conf_fact,10) ;
276 kuu.std_spectral_base();
277 Sym_tensor kdd = kuu.up_down(gam);
278
279 //---------------------------------------------------------
280 // Construction of the spheroid associated with those data
281 //--------------------------------------------------------
282 double hor_posd = rrr.val_grid_point(1,0,0,0);
283 Scalar hor_pos(map_2); hor_pos = hor_posd; hor_pos.std_spectral_base();
284 Spheroid hor_loc(hor_pos, gam, kdd);
285 return hor_loc;
286}
287
288// Computation of the ADM mass of the BH spacetime
290 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
291 const Metric_flat& mets = (*map).flat_met_spher() ;
292
293 Sym_tensor gamtcon = mets.con() + hij;
294 Metric gamt(gamtcon);
295 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
296 Metric gam(gamcon);
297
298 Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
299 detgam.std_spectral_base();
300 Vector corr_adm = - (0.125*contract(gamt.cov().derive_con(mets),1,2));
301 Scalar admintegr = conf_fact.dsdr() + corr_adm(1);
302
303 double M_ADM = - (1/(2.*3.1415927*ggrav))*(*map).integrale_surface_infini(admintegr*detgam);
304 return M_ADM;
305}
306
307
308// Computation of the Komar angular momentum w.r.t. assumed rotational symmetry
310const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
311 const Metric_flat& mets = (*map).flat_met_spher() ;
312
313 Sym_tensor gamtcon = mets.con() + hij;
314 Sym_tensor gamcon = gamtcon/(conf_fact*conf_fact*conf_fact*conf_fact);
315 gamcon.std_spectral_base();
316 Metric gam(gamcon);
317 Sym_tensor k_uu = hatA/(pow(conf_fact,10));
318 k_uu.std_spectral_base();
319 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); //WTF?
320 Sym_tensor k_dd = k_uu.up_down(gam);
321
322 Scalar detgam = sqrt((gam.cov())(2,2)*(gam.cov())(3,3) - (gam.cov())(2,3)*(gam.cov())(3,2));
323 detgam.std_spectral_base();
324 Scalar contraction = k_dd(1,3); contraction.mult_r_dzpuis(2); contraction.mult_sint();
325 double angu_komar = - (1/(8.*3.1415927*ggrav))*(*map).integrale_surface_infini(detgam*contraction);
326
327 return angu_komar;
328}
329
330
331// Computation of the Virial residual, as rescaled difference at infinity
332// between the ADM mass and the Komar integral associated to the mass
333
335 const Mg3d* mgrid = mp.get_mg();
336 const int nz = (*mgrid).get_nzone(); // Number of domains
337 Valeur** devel_psi (conf_fact.asymptot(1)) ;
338 Valeur** devel_n (lapse.asymptot(1)) ;
339
340
341 double erreur = (2*(*devel_psi[1])(nz-1, 0, 0, 0)
342 + (*devel_n[1])(nz-1, 0, 0, 0))/fabs ((*devel_n[1])(nz-1, 0, 0, 0)) ;
343
344 return erreur;
345}
346
347
348
349
350
351}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Class to compute quasistationary single black hole spacetimes in vacuum.
Definition isol_hole.h:50
Spheroid hor()
Spheroid at the excised boundary associated with the black hole MOTS on the slice.
Definition isol_hole.C:247
double * p_adm_mass
Computation of the ADM mass of the BH spacetime.
Definition isol_hole.h:103
double * p_komar_angmom
Computation of the Komar angular momentum w.r.t.
Definition isol_hole.h:108
virtual ~Isol_hole()
Destructor.
Definition isol_hole.C:125
void Einstein_errors()
Prints out errors in Einstein equations for the data obtained.
Definition isol_hole.C:194
Scalar boundNoK
Indicates if the boundary value for the lapse or the surface gravity is used.
Definition isol_hole.h:65
virtual void del_deriv() const
Deletes all the derived quantities.
Definition isol_hole.C:136
Scalar lapse
Lapse function.
Definition isol_hole.h:76
bool NorKappa
Rotation rate of the horizon in the azimuthal direction.
Definition isol_hole.h:61
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
Definition isol_hole.h:87
double * p_virial_residue
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Definition isol_hole.h:113
void operator=(const Isol_hole &)
Assignment to another Isol_hole.
Definition isol_hole.C:160
Vector shift
Shift vector.
Definition isol_hole.h:82
virtual void sauve(FILE *) const
Save in a file.
Definition isol_hole.C:183
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: see Cordero et al.
Definition isol_hole.h:92
const Map & mp
Mapping associated with the star.
Definition isol_hole.h:55
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition isol_hole.C:145
bool isCF
Stores the boundary value of the lapse or surface gravity.
Definition isol_hole.h:69
double adm_mass()
Computation of the ADM mass of the BH spacetime.
Definition isol_hole.C:289
Spheroid * p_hor
Computation of the spheroid associated with the black hole horizon.
Definition isol_hole.h:101
Isol_hole(const Map &mp_i, double Omega_i, bool NorKappa_i, Scalar NoK_i, bool isCF_i=false)
Standard constructor.
Definition isol_hole.C:51
double komar_angmom()
Computation of the Komar angular momentum w.r.t.
Definition isol_hole.C:309
double virial_residue()
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Definition isol_hole.C:334
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
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
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.