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