LORENE
blackhole_bc.C
1 /*
2  * Methods of class Black_hole to compute the inner boundary condition
3  * at the excised surface
4  *
5  * (see file blackhole.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005-2007 Keisuk 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 char blackhole_bc_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $" ;
30 
31 /*
32  * $Id: blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $
33  * $Log: blackhole_bc.C,v $
34  * Revision 1.5 2014/10/13 08:52:45 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.4 2014/10/06 15:13:02 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.3 2008/07/03 14:53:47 k_taniguchi
41  * Modification of a signature in bc_shift_x and bc_shift_y.
42  *
43  * Revision 1.2 2008/05/15 19:25:43 k_taniguchi
44  * Change of some parameters.
45  *
46  * Revision 1.1 2007/06/22 01:18:23 k_taniguchi
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $
51  *
52  */
53 
54 // C++ headers
55 //#include <>
56 
57 // C headers
58 #include <cmath>
59 
60 // Lorene headers
61 #include "blackhole.h"
62 #include "valeur.h"
63 #include "grilles.h"
64 #include "unites.h"
65 
66  //----------------------------------//
67  // Inner boundary condition //
68  //----------------------------------//
69 
70 namespace Lorene {
71 const Valeur Black_hole::bc_lapconf(bool neumann, bool first) const {
72 
73  // Fundamental constants and units
74  // -------------------------------
75  using namespace Unites ;
76 
77  const Mg3d* mg = mp.get_mg() ;
78  const Mg3d* mg_angu = mg->get_angu() ;
79  Valeur bc(mg_angu) ;
80 
81  if (kerrschild) {
82 
83  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
84  abort() ;
85 
86  /*
87  if (neumann) {
88 
89  if (first) {
90 
91  Scalar rr(mp) ;
92  rr = mp.r ;
93  rr.std_spectral_base() ;
94 
95  int nt = mg->get_nt(0) ;
96  int np = mg->get_np(0) ;
97 
98  Scalar tmp(mp) ;
99 
100  tmp = - pow(lapse_bh,3.) * ggrav * mass_bh / rr / rr ;
101  // dlapse/dr = 0
102 
103  bc = 1. ;
104  for (int j=0; j<nt; j++) {
105  for (int k=0; k<np; k++) {
106  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
107  }
108  }
109  }
110  else {
111 
112  bc = 0. ; // dlapse/dr = 0.25*lapse/rr
113 
114  }
115  }
116  else {
117  if (first) { // The poisson solver in LORENE assumes the
118  // asymptotic behavior of the function -> 0
119  bc = 0.5 - 1./sqrt(2.) ; // <- bc of the real function = 0.5
120 
121  }
122  else {
123 
124  bc = 0. ; // <- bc of the real function = 1./sqrt(2.)
125 
126  }
127  }
128  */
129  }
130  else { // Isotropic coordinates with the maximal slicing
131 
132  if (neumann) {
133 
134  if (first) {
135 
136  bc = 0. ; // d(lapconf)/dr = 0
137 
138  }
139  else {
140 
141  Scalar rr(mp) ;
142  rr = mp.r ;
143  rr.std_spectral_base() ;
144 
145  int nt = mg->get_nt(0) ;
146  int np = mg->get_np(0) ;
147 
148  Scalar tmp(mp) ;
149 
150  tmp = 0.5 * lapconf / rr ;
151  // d(lapconf)/dr = lapconf/2/rah
152 
153  bc = 1. ;
154  for (int j=0; j<nt; j++) {
155  for (int k=0; k<np; k++) {
156  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
157  }
158  }
159 
160 
161  }
162 
163  }
164  else {
165 
166  if (first) {
167 
168  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
169  abort() ;
170  // bc = - 0.5 ; // lapconf = 0.5
171 
172  }
173  else {
174 
175  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
176  abort() ;
177  /*
178  Scalar r_are(mp) ;
179  r_are = r_coord(neumann, first) ;
180  r_are.std_spectral_base() ;
181  r_are.annule_domain(0) ;
182  r_are.raccord(1) ;
183 
184  int nt = mg->get_nt(0) ;
185  int np = mg->get_np(0) ;
186 
187  Scalar tmp(mp) ;
188 
189  tmp = sqrt(0.5*r_are) - 1. ; // lapse = 1./sqrt(2.)
190 
191  bc = 1. ;
192  for (int j=0; j<nt; j++) {
193  for (int k=0; k<np; k++) {
194  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
195  }
196  }
197  */
198 
199  }
200 
201  }
202 
203  }
204 
205  bc.std_base_scal() ;
206  return bc ;
207 
208 }
209 
210 
211 const Valeur Black_hole::bc_shift_x(double omega_r) const {
212 
213  // Fundamental constants and units
214  // -------------------------------
215  using namespace Unites ;
216 
217  const Mg3d* mg = mp.get_mg() ;
218  const Mg3d* mg_angu = mg->get_angu() ;
219  Valeur bc(mg_angu) ;
220 
221  Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
222 
223  Scalar rr(mp) ;
224  rr = mp.r ;
225  rr.std_spectral_base() ;
226  Scalar st(mp) ;
227  st = mp.sint ;
228  st.std_spectral_base() ;
229  Scalar cp(mp) ;
230  cp = mp.cosp ;
231  cp.std_spectral_base() ;
232  Scalar yy(mp) ;
233  yy = mp.y ;
234  yy.std_spectral_base() ;
235 
236  Scalar tmp(mp) ;
237 
238  if (kerrschild) {
239 
240  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
241  abort() ;
242  /*
243  tmp = lapse_bh * (lapse / confo / confo) * st * cp
244  - omega_r * yy - shift_bh(1) ;
245  */
246 
247  // tmp = lap_bh * lap_bh * st * cp - omega_r * yy ;
248  }
249  else { // Isotropic coordinates with the maximal slicing
250 
251  // Note: the signature of omega_r is opposite to that in the
252  // binary case because of the direction of the spin
253  tmp = lapconf / pow(confo, 3.) * st * cp + omega_r * yy ;
254  // tmp = lapconf / pow(confo, 3.) * st * cp - omega_r * yy ;
255 
256  }
257 
258  int nt = mg->get_nt(0) ;
259  int np = mg->get_np(0) ;
260 
261  bc = 1. ;
262  for (int j=0; j<nt; j++) {
263  for (int k=0; k<np; k++) {
264  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
265  }
266  }
267 
268  bc.base = *bases[0] ;
269  // bc.std_base_scal() ;
270 
271  for (int i=0; i<3; i++)
272  delete bases[i] ;
273 
274  delete [] bases ;
275 
276  return bc ;
277 
278 }
279 
280 
281 const Valeur Black_hole::bc_shift_y(double omega_r) const {
282 
283  // Fundamental constants and units
284  // -------------------------------
285  using namespace Unites ;
286 
287  const Mg3d* mg = mp.get_mg() ;
288  const Mg3d* mg_angu = mg->get_angu() ;
289  Valeur bc(mg_angu) ;
290 
291  Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
292 
293  Scalar rr(mp) ;
294  rr = mp.r ;
295  rr.std_spectral_base() ;
296  Scalar st(mp) ;
297  st = mp.sint ;
298  st.std_spectral_base() ;
299  Scalar sp(mp) ;
300  sp = mp.sinp ;
301  sp.std_spectral_base() ;
302  Scalar xx(mp) ;
303  xx = mp.x ;
304  xx.std_spectral_base() ;
305 
306  Scalar tmp(mp) ;
307 
308  if (kerrschild) {
309 
310  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
311  abort() ;
312  /*
313  tmp = lapse_bh * (lapse / confo / confo) * st * sp
314  + omega_r * xx - shift_bh(2) ;
315  */
316  // tmp = lap_bh * lap_bh * st * sp + omega_r * xx ;
317  }
318  else {
319 
320  // Note: the signature of omega_r is opposite to that in the
321  // binary case because of the direction of the spin
322  tmp = lapconf / pow(confo, 3.) * st * sp - omega_r * xx ;
323  // tmp = lapconf / pow(confo, 3.) * st * sp + omega_r * xx ;
324 
325  }
326 
327  int nt = mg->get_nt(0) ;
328  int np = mg->get_np(0) ;
329 
330  bc = 1. ;
331  for (int j=0; j<nt; j++) {
332  for (int k=0; k<np; k++) {
333  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
334  }
335  }
336 
337  bc.base = *bases[1] ;
338  // bc.std_base_scal() ;
339 
340  for (int i=0; i<3; i++)
341  delete bases[i] ;
342 
343  delete [] bases ;
344 
345  return bc ;
346 
347 }
348 
349 
351 
352  // Fundamental constants and units
353  // -------------------------------
354  using namespace Unites ;
355 
356  const Mg3d* mg = mp.get_mg() ;
357  const Mg3d* mg_angu = mg->get_angu() ;
358  Valeur bc(mg_angu) ;
359 
360  Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
361 
362  Scalar rr(mp) ;
363  rr = mp.r ;
364  rr.std_spectral_base() ;
365  Scalar ct(mp) ;
366  ct = mp.cost ;
367  ct.std_spectral_base() ;
368 
369  Scalar tmp(mp) ;
370 
371  if (kerrschild) {
372 
373  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
374  abort() ;
375  /*
376  tmp = lapse_bh * (lapse / confo / confo) * ct
377  - shift_bh(3) ;
378  */
379  // tmp = lap_bh * lap_bh * ct ;
380  }
381  else {
382 
383  tmp = lapconf / pow(confo, 3.) * ct ;
384 
385  }
386 
387  int nt = mg->get_nt(0) ;
388  int np = mg->get_np(0) ;
389 
390  bc = 1. ;
391  for (int j=0; j<nt; j++) {
392  for (int k=0; k<np; k++) {
393  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
394  }
395  }
396 
397  bc.base = *bases[2] ;
398  // bc.std_base_scal() ;
399 
400  for (int i=0; i<3; i++)
401  delete bases[i] ;
402 
403  delete [] bases ;
404 
405  return bc ;
406 
407 }
408 
409 
411 
412  // Fundamental constants and units
413  // -------------------------------
414  using namespace Unites ;
415 
416  const Mg3d* mg = mp.get_mg() ;
417  const Mg3d* mg_angu = mg->get_angu() ;
418  Valeur bc(mg_angu) ;
419 
420  double mass = ggrav * mass_bh ;
421 
422  Scalar rr(mp) ;
423  rr = mp.r ;
424  rr.std_spectral_base() ;
425  Scalar st(mp) ;
426  st = mp.sint ;
427  st.std_spectral_base() ;
428  Scalar ct(mp) ;
429  ct = mp.cost ;
430  ct.std_spectral_base() ;
431  Scalar sp(mp) ;
432  sp = mp.sinp ;
433  sp.std_spectral_base() ;
434  Scalar cp(mp) ;
435  cp = mp.cosp ;
436  cp.std_spectral_base() ;
437 
438  int nt = mg->get_nt(0) ;
439  int np = mg->get_np(0) ;
440 
441  Scalar tmp(mp) ;
442 
443  if (kerrschild) { // Assumes that r_BH = 1.
444 
445  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
446  abort() ;
447  /*
448  Scalar divshift(mp) ;
449  divshift = shift_rs(1).deriv(1) + shift_rs(2).deriv(2)
450  + shift_rs(3).deriv(3) ;
451  divshift.std_spectral_base() ;
452 
453  Scalar llshift(mp) ;
454  llshift = st*cp*shift_rs(1) + st*sp*shift_rs(2) + ct*shift_rs(3) ;
455  llshift.std_spectral_base() ;
456 
457  Scalar lldllsh = llshift.dsdr() ;
458  lldllsh.std_spectral_base() ;
459 
460  Scalar tmp1 = divshift ;
461  Scalar tmp2 = -3.*lldllsh ;
462 
463  Scalar tmp5 = 0.5*confo*(lapse_bh*confo*confo/lapse - 1.)/rr ;
464  tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
465  tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
466 
467  Scalar tmp3 = 2. * lapse_bh * lapse_bh * mass * llshift / rr / rr ;
468  Scalar tmp4 = 4. * pow(lapse_bh,3.) * mass * (1.+3.*mass/rr)
469  * lapse_rs / rr / rr ;
470 
471  tmp3.set_dzpuis(tmp5.get_dzpuis()) ;
472  tmp4.set_dzpuis(tmp5.get_dzpuis()) ;
473 
474  tmp = tmp5 + pow(confo,3.)*(tmp1+tmp2+tmp3+tmp4)/12./lapse/lapse_bh ;
475  */
476  // tmp = -0.5 * (1. - 2. * mass / rr) / rr ;
477 
478  }
479  else { // Isotropic coordinates with the maximal slicing
480 
481  Scalar divshift(mp) ;
482  divshift = shift(1).deriv(1) + shift(2).deriv(2)
483  + shift(3).deriv(3) ;
484  divshift.std_spectral_base() ;
485 
486  Scalar llshift(mp) ;
487  llshift = st*cp*shift(1) + st*sp*shift(2) + ct*shift(3) ;
488  llshift.std_spectral_base() ;
489 
490  Scalar lldllsh = llshift.dsdr() ;
491  lldllsh.std_spectral_base() ;
492 
493  Scalar tmp1 = divshift ;
494  Scalar tmp2 = -3.*lldllsh ;
495 
496  Scalar tmp5 = - 0.5 * confo / rr ;
497 
498  tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
499  tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
500 
501  tmp = tmp5 + pow(confo, 4.) * (tmp1 + tmp2) / 12. / lapconf ;
502 
503  }
504 
505  bc = 1. ;
506  for (int j=0; j<nt; j++) {
507  for (int k=0; k<np; k++) {
508  bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
509  }
510  }
511 
512  bc.std_base_scal() ;
513  return bc ;
514 
515 }
516 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
const Valeur bc_lapconf(bool neumann, bool first) const
Boundary condition on the apparent horizon of the black hole for the lapse function: 2-D Valeur...
Definition: blackhole_bc.C:71
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene prototypes.
Definition: app_hor.h:64
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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:784
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:824
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
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
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Coord sint
Definition: map.h:721
const Valeur bc_shift_x(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: blackhole_bc.C:211
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
const Valeur bc_shift_y(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: blackhole_bc.C:281
Coord sinp
Definition: map.h:723
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: blackhole_bc.C:350
Multi-domain grid.
Definition: grilles.h:273
Bases of the spectral expansions.
Definition: base_val.h:322
Coord y
y coordinate centered on the grid
Definition: map.h:727
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Coord cosp
Definition: map.h:724
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Coord x
x coordinate centered on the grid
Definition: map.h:726
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition: blackhole.h:97
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur...
Definition: blackhole_bc.C:410
Coord r
r coordinate centered on the grid
Definition: map.h:718
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
Coord cost
Definition: map.h:722