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 
45 namespace Lorene {
46 
47  // Fundamental constants and units
48  // -------------------------------
49  using namespace Unites ;
50 
51  //--------------//
52  // Constructors //
53  //--------------//
54 
55 
56 // Standard constructor
57 // --------------------
58 Excised_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. ;
79  Xx.set_etat_zero();
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 // -----------------------
101 Excised_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 
125 
126  del_deriv() ;
127 
128 }
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 // --------------
181 void 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 
195 const 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
309 const 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 }
Vector shift
Shift vector.
Definition: excised_slice.h:64
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: rescaled same way as Cordero et al.
Definition: excised_slice.h:74
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.
Metric for tensor calculation.
Definition: metric.h:90
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:375
void Einstein_errors()
Prints out errors in Einstein equations for the data obtained.
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
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:777
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:682
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
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...
void mult_sint()
Multiplication by .
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
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition: metric.C:341
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:915
virtual void sauve(FILE *) const
Save in a file.
Class to compute single black hole spacetime excised slices.
Definition: excised_slice.h:41
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
int field_set
Chose field set type.
Definition: excised_slice.h:55
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 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:817
double * p_adm_mass
Computation of the ADM mass of the BH spacetime.
Definition: excised_slice.h:91
void operator=(const Excised_slice &)
Assignment to another Excised_slice.
virtual void del_deriv() const
Deletes all the derived quantities.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:692
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Definition: spheroid.h:84
double * p_virial_residue
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
double komar_angmom()
Computation of the Komar angular momentum w.r.t.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:353
Vector Xx
Longitudinal part of the extrinsic curvature.
Definition: excised_slice.h:79
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al...
Definition: excised_slice.h:69
Spheroid * p_hor
Computation of the spheroid associated with the black hole horizon.
Definition: excised_slice.h:89
const Map & mp
Mapping associated with the slice.
Definition: excised_slice.h:46
Spheroid hor()
Spheroid at the excised boundary associated with the black hole MOTS on the slice.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int type_hor
Chosen horizon type.
Definition: excised_slice.h:52
Scalar lapse
Lapse function.
Definition: excised_slice.h:58
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
double virial_residue()
Computation of the Virial residual, as difference at infinity between the ADM mass and the Komar inte...
Excised_slice(const Map &mp_i, int hor_type, int scheme_type)
Standard constructor.
Definition: excised_slice.C:58
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
double adm_mass()
Computation of the ADM mass of the BH spacetime.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Multi-domain grid.
Definition: grilles.h:279
Affine radial mapping.
Definition: map.h:2042
virtual ~Excised_slice()
Destructor.
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:625
double * p_komar_angmom
Computation of the Komar angular momentum w.r.t.
Definition: excised_slice.h:96
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...
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:935
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:363
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r ...
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.