LORENE
hole_bhns_global.C
1 /*
2  * Methods of class Hole_bhns to compute global quantities
3  *
4  * (see file hole_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005,2007 Keisuke Taniguchi
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 
29 
30 /*
31  * $Id: hole_bhns_global.C,v 1.6 2016/12/05 16:17:55 j_novak Exp $
32  * $Log: hole_bhns_global.C,v $
33  * Revision 1.6 2016/12/05 16:17:55 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.5 2014/10/13 08:53:00 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.4 2014/10/06 15:13:10 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.3 2008/07/02 21:10:15 k_taniguchi
43  * A bug removed.
44  *
45  * Revision 1.2 2008/05/15 19:07:26 k_taniguchi
46  * Introduction of the quasilocal spin angular momentum.
47  *
48  * Revision 1.1 2007/06/22 01:25:15 k_taniguchi
49  * *** empty log message ***
50  *
51  *
52  * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_global.C,v 1.6 2016/12/05 16:17:55 j_novak Exp $
53  *
54  */
55 
56 // C++ headers
57 //#include <>
58 
59 // C headers
60 #include <cmath>
61 
62 // Lorene headers
63 #include "hole_bhns.h"
64 #include "unites.h"
65 #include "utilitaires.h"
66 
67  //-----------------------------------------//
68  // Irreducible mass of BH //
69  //-----------------------------------------//
70 
71 namespace Lorene {
72 double Hole_bhns::mass_irr_bhns() const {
73 
74  // Fundamental constants and units
75  // -------------------------------
76  using namespace Unites ;
77 
78  if (p_mass_irr_bhns == 0x0) { // a new computation is required
79 
80  Scalar psi4(mp) ;
81  psi4 = pow(confo_tot, 4.) ;
82  psi4.std_spectral_base() ;
83  psi4.annule_domain(0) ;
84  psi4.raccord(1) ;
85 
86  double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
87 
88  Map_af& mp_aff= dynamic_cast<Map_af&>(mp) ;
89 
90  double a_ah = mp_aff.integrale_surface(psi4, radius_ah) ;
91  double mirr = sqrt(a_ah/16./M_PI) / ggrav ;
92 
93  p_mass_irr_bhns = new double( mirr ) ;
94 
95  }
96 
97  return *p_mass_irr_bhns ;
98 
99 }
100 
101  //----------------------------------------------------------//
102  // Quasilocal spin angular momentum of BH //
103  //----------------------------------------------------------//
104 
105 double Hole_bhns::spin_am_bhns(const Tbl& xi_i, const double& phi_i,
106  const double& theta_i, const int& nrk_phi,
107  const int& nrk_theta) const {
108 
109  // Fundamental constants and units
110  // -------------------------------
111  using namespace Unites ;
112 
113  if (p_spin_am_bhns == 0x0) { // a new computation is required
114 
115  double mass = ggrav * mass_bh ;
116 
117  Scalar rr(mp) ;
118  rr = mp.r ;
119  rr.std_spectral_base() ;
120 
121  Scalar st(mp) ;
122  st = mp.sint ;
123  st.std_spectral_base() ;
124  Scalar ct(mp) ;
125  ct = mp.cost ;
126  ct.std_spectral_base() ;
127  Scalar sp(mp) ;
128  sp = mp.sinp ;
129  sp.std_spectral_base() ;
130  Scalar cp(mp) ;
131  cp = mp.cosp ;
132  cp.std_spectral_base() ;
133 
134  Vector ll(mp, CON, mp.get_bvect_cart()) ;
135  ll.set_etat_qcq() ;
136  ll.set(1) = st % cp ;
137  ll.set(2) = st % sp ;
138  ll.set(3) = ct ;
139  ll.std_spectral_base() ;
140 
141  double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
142 
143  if (kerrschild) {
144 
145  cout << "Not yet prepared!!!" << endl ;
146  abort() ;
147 
148  }
149  else { // Isotropic coordinates
150 
151  // Sets C/M^2 for each case of the lapse boundary condition
152  // --------------------------------------------------------
153  double cc ;
154 
155  if (bc_lapconf_nd) { // Neumann boundary condition
156  if (bc_lapconf_fs) { // First condition
157  // d(\alpha \psi)/dr = 0
158  // ---------------------
159  cc = 2. * (sqrt(13.) - 1.) / 3. ;
160  }
161  else { // Second condition
162  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
163  // -----------------------------------------
164  cc = 4. / 3. ;
165  }
166  }
167  else { // Dirichlet boundary condition
168  if (bc_lapconf_fs) { // First condition
169  // (\alpha \psi) = 1/2
170  // -------------------
171  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
172  abort() ;
173  }
174  else { // Second condition
175  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
176  // ----------------------------------
177  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
178  abort() ;
179  // cc = 2. * sqrt(2.) ;
180  }
181  }
182 
183  Scalar r_are(mp) ;
185  r_are.std_spectral_base() ;
186 
187  // Killing vector of the spherical components
188  Vector killing_spher(mp, COV, mp.get_bvect_spher()) ;
189  killing_spher.set_etat_qcq() ;
190  killing_spher = killing_vect(xi_i, phi_i, theta_i,
191  nrk_phi, nrk_theta) ;
192  killing_spher.std_spectral_base() ;
193 
194  killing_spher.set(2) = confo_tot * confo_tot * radius_ah
195  * killing_spher(2) ;
196  killing_spher.set(3) = confo_tot * confo_tot * radius_ah
197  * killing_spher(3) ;
198  // killing_spher(3) is already divided by sin(theta)
199  killing_spher.std_spectral_base() ;
200 
201  // Killing vector of the Cartesian components
202  Vector killing(mp, COV, mp.get_bvect_cart()) ;
203  killing.set_etat_qcq() ;
204  killing.set(1) = (killing_spher(2) * ct * cp - killing_spher(3) * sp)
205  / radius_ah ;
206  killing.set(2) = (killing_spher(2) * ct * sp + killing_spher(3) * cp)
207  / radius_ah ;
208  killing.set(3) = - killing_spher(2) * st / radius_ah ;
209  killing.std_spectral_base() ;
210 
211  // Surface integral <- dzpuis should be 0
212  // --------------------------------------
213  // Source terms in the surface integral
214  Scalar source_1(mp) ;
215  source_1 = (ll(1) * (taij_tot_rs(1,1) * killing(1)
216  + taij_tot_rs(1,2) * killing(2)
217  + taij_tot_rs(1,3) * killing(3))
218  + ll(2) * (taij_tot_rs(2,1) * killing(1)
219  + taij_tot_rs(2,2) * killing(2)
220  + taij_tot_rs(2,3) * killing(3))
221  + ll(3) * (taij_tot_rs(3,1) * killing(1)
222  + taij_tot_rs(3,2) * killing(2)
223  + taij_tot_rs(3,3) * killing(3)))
224  / pow(confo_tot, 4.) ;
225  source_1.std_spectral_base() ;
226  source_1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
227 
228  Scalar source_2(mp) ;
229  source_2 = -2. * pow(confo_tot, 3.) * mass * mass * cc
230  * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
231  * (ll(1)*killing(1) + ll(2)*killing(2) + ll(3)*killing(3))
232  / lapconf_tot / pow(r_are*rr, 3.) ;
233  source_2.std_spectral_base() ;
234 
235  Scalar source_surf(mp) ;
236  source_surf = source_1 + source_2 ;
237  source_surf.std_spectral_base() ;
238  source_surf.annule_domain(0) ;
239  source_surf.raccord(1) ;
240 
241  Map_af& mp_aff= dynamic_cast<Map_af&>(mp) ;
242 
243  double spin = mp_aff.integrale_surface(source_surf, radius_ah) ;
244  double spin_angmom = 0.5 * spin / qpig ;
245 
246  p_spin_am_bhns = new double( spin_angmom ) ;
247 
248  }
249 
250  }
251 
252  return *p_spin_am_bhns ;
253 
254 }
255 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Standard units of space, time and mass.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
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
Tensor field of valence 1.
Definition: vector.h:188
double * p_spin_am_bhns
Irreducible mass of BH.
Definition: hole_bhns.h:248
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Coord sint
Definition: map.h:739
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Vector killing_vect(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Compute the Killing vector of a black hole normalized so that its affine length is 2 M_PI...
Coord sinp
Definition: map.h:741
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Scalar confo_tot
Total conformal factor.
Definition: hole_bhns.h:169
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Definition: hole_bhns.h:190
Coord cosp
Definition: map.h:742
virtual double mass_irr_bhns() const
Irreducible mass of the black hole.
bool bc_lapconf_fs
true for the first type BC for the lapconf function, false for the second type BH ...
Definition: hole_bhns.h:78
bool bc_lapconf_nd
true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH ...
Definition: hole_bhns.h:73
Basic array class.
Definition: tbl.h:164
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
double spin_am_bhns(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Spin angular momentum of the black hole.
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...
Scalar lapconf_tot
Total lapconf function.
Definition: hole_bhns.h:101
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord cost
Definition: map.h:740