LORENE
blackhole_static.C
1 /*
2  * Method of class Black_hole to set metric quantities to a spherical,
3  * static, analytic solution
4  *
5  * (see file blackhole.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005-2007 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: blackhole_static.C,v 1.4 2016/12/05 16:17:48 j_novak Exp $
33  * $Log: blackhole_static.C,v $
34  * Revision 1.4 2016/12/05 16:17:48 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.3 2014/10/13 08:52:46 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.2 2008/05/15 19:31:17 k_taniguchi
41  * Change of some parameters.
42  *
43  * Revision 1.1 2007/06/22 01:20:50 k_taniguchi
44  * *** empty log message ***
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_static.C,v 1.4 2016/12/05 16:17:48 j_novak Exp $
48  *
49  */
50 
51 // C++ headers
52 //#include <>
53 
54 // C headers
55 //#include <math.h>
56 
57 // Lorene headers
58 #include "blackhole.h"
59 #include "unites.h"
60 #include "utilitaires.h"
61 
62 namespace Lorene {
63 void Black_hole::static_bh(bool neumann, bool first) {
64 
65  // Fundamental constants and units
66  // -------------------------------
67  using namespace Unites ;
68 
69  double mass = ggrav * mass_bh ;
70 
71  Scalar rr(mp) ;
72  rr = mp.r ;
73  rr.std_spectral_base() ;
74 
75  //-------------------------------------//
76  // Metric quantities //
77  //-------------------------------------//
78 
79  Scalar st(mp) ;
80  st = mp.sint ;
81  st.std_spectral_base() ;
82  Scalar ct(mp) ;
83  ct = mp.cost ;
84  ct.std_spectral_base() ;
85  Scalar sp(mp) ;
86  sp = mp.sinp ;
87  sp.std_spectral_base() ;
88  Scalar cp(mp) ;
89  cp = mp.cosp ;
90  cp.std_spectral_base() ;
91 
92  Vector ll(mp, CON, mp.get_bvect_cart()) ;
93  ll.set_etat_qcq() ;
94  ll.set(1) = st * cp ;
95  ll.set(2) = st * sp ;
96  ll.set(3) = ct ;
97  ll.std_spectral_base() ;
98 
99  if (kerrschild) {
100 
101  // Lapconf function
102  // ----------------
103  lapconf = 1./sqrt(1.+2.*mass/rr) ;
106  lapconf.raccord(1) ;
107 
108  // Conformal factor
109  // ----------------
110  confo = 1. ;
112 
113  // Shift vector
114  // ------------
115  for (int i=1; i<=3; i++)
116  shift.set(i) = 2. * mass * lapconf % lapconf % ll(i) / rr ;
117 
118  shift.annule_domain(0) ;
120 
121  }
122  else { // Isotropic coordinates with the maximal slicing
123 
124  // Sets C/M^2 for each case of the lapse boundary condition
125  // --------------------------------------------------------
126  double cc ;
127 
128  if (neumann) { // Neumann boundary condition
129  if (first) { // First condition
130  // d(\alpha \psi)/dr = 0
131  // ---------------------
132  cc = 2. * (sqrt(13.) - 1.) / 3. ;
133  }
134  else { // Second condition
135  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
136  // -----------------------------------------
137  cc = 4. / 3. ;
138  }
139  }
140  else { // Dirichlet boundary condition
141  if (first) { // First condition
142  // (\alpha \psi) = 1/2
143  // -------------------
144  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
145  abort() ;
146  }
147  else { // Second condition
148  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
149  // ----------------------------------
150  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
151  abort() ;
152  }
153  }
154 
155  Scalar r_are(mp) ;
156  r_are = r_coord(neumann, first) ;
157  r_are.std_spectral_base() ;
158 
159  // Lapconf function
160  // ----------------
161  lapconf = sqrt(1. - 2.*mass/r_are/rr
162  + cc*cc*pow(mass/r_are/rr,4.)) * sqrt(r_are) ;
163 
166  lapconf.raccord(1) ;
167 
168  // Conformal factor
169  // ----------------
170  confo = sqrt(r_are) ;
172  confo.annule_domain(0) ;
173  confo.raccord(1) ;
174 
175  // Lapse function
176  // --------------
177  lapse = lapconf / confo ;
179  lapse.annule_domain(0) ;
180  lapse.raccord(1) ;
181 
182  // Shift vector
183  // ------------
184  for (int i=1; i<=3; i++) {
185  shift.set(i) = mass * mass * cc * ll(i) / rr / rr
186  / pow(r_are,3.) ;
187  }
188 
190 
191  for (int i=1; i<=3; i++) {
192  shift.set(i).annule_domain(0) ;
193  shift.set(i).raccord(1) ;
194  }
195 
196  }
197 
198 }
199 }
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
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
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.
void static_bh(bool neumann, bool first)
Sets the metric quantities to a spherical, static blach-hole analytic solution.
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
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Coord sint
Definition: map.h:739
Coord sinp
Definition: map.h:741
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
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
Scalar lapse
Lapse function generated by the black hole.
Definition: blackhole.h:106
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Coord cosp
Definition: map.h:742
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition: blackhole.h:97
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord cost
Definition: map.h:740