LORENE
blackhole_rk_phi.C
1 /*
2  * Methods of class Black_hole to compute a forth-order Runge-Kutta
3  * integration to the phi direction for the solution of the Killing vectors
4  *
5  * (see file blackhole.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 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_rk_phi.C,v 1.5 2016/12/05 16:17:48 j_novak Exp $
33  * $Log: blackhole_rk_phi.C,v $
34  * Revision 1.5 2016/12/05 16:17:48 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.4 2014/10/13 08:52:46 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.3 2014/10/06 15:13:02 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.2 2008/07/02 20:43:54 k_taniguchi
44  * Typos removed.
45  *
46  * Revision 1.1 2008/05/15 19:33:32 k_taniguchi
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_rk_phi.C,v 1.5 2016/12/05 16:17:48 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 "unites.h"
63 #include "utilitaires.h"
64 
65  //--------------------------------------------------//
66  // Forth-order Runge-Kutta on the equator //
67  //--------------------------------------------------//
68 
69 namespace Lorene {
70 Tbl Black_hole::runge_kutta_phi_bh(const Tbl& xi_i, const double& phi_i,
71  const int& nrk_phi) const {
72 
73  using namespace Unites ;
74 
75  const Mg3d* mg = mp.get_mg() ;
76  int np = mg->get_np(1) ;
77 
78  Tbl xi_f(3) ; // xi_f(0)=xi_hat{theta}, xi_f(1)=xi_hat{phi}, xi_f(2)=L
79  xi_f.set_etat_qcq() ;
80 
81  if (kerrschild) {
82 
83  cout << "Not yet prepared!!!" << endl ;
84  abort() ;
85 
86  }
87  else { // Isotropic coordinates
88 
89  // Initial data at phi=0 on the equator
90  double xi_t0 = xi_i(0) ; // xi_hat{theta}
91  double xi_p0 = xi_i(1) ; // xi_hat{phi}
92  double xi_l0 = xi_i(2) ; // L
93  double phi0 = phi_i ;
94 
95  double dp = 2. * M_PI / double(np) / double(nrk_phi) ;
96 
97  double rah = rad_ah() ;
98 
99  Scalar dlnconfo(mp) ;
100  dlnconfo = confo.dsdt() / confo ;
101  dlnconfo.std_spectral_base() ;
102 
103  Scalar laplnconfo(mp) ;
104  laplnconfo = confo.lapang() / confo ;
105  laplnconfo.std_spectral_base() ;
106 
107  Scalar confo2(mp) ;
108  confo2 = confo * confo ;
109  confo2.std_spectral_base() ;
110 
111  double xi_t1, xi_t2, xi_t3, xi_t4, xi_tf ;
112  double xi_p1, xi_p2, xi_p3, xi_p4, xi_pf ;
113  double xi_l1, xi_l2, xi_l3, xi_l4, xi_lf ;
114  double f1, f2, f3, f4 ;
115  double g1, g2, g3, g4 ;
116  double h1, h2, h3, h4 ;
117 
118  // Forth-order Runge-Kutta
119  // (nrk_phi times steps between two collocation points)
120  // ----------------------------------------------------
121 
122  for (int i=0; i<nrk_phi; i++) {
123 
124  // First
125  f1 = - xi_l0 * rah * confo2.val_point(rah, M_PI/2., phi0)
126  + 2. * xi_p0 * dlnconfo.val_point(rah, M_PI/2., phi0) ;
127  g1 = -2. * xi_t0 * dlnconfo.val_point(rah, M_PI/2., phi0) ;
128  h1 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0)) * xi_t0
129  / rah / confo2.val_point(rah, M_PI/2., phi0) ;
130 
131  xi_t1 = dp * f1 ;
132  xi_p1 = dp * g1 ;
133  xi_l1 = dp * h1 ;
134 
135  // Second
136  f2 = - (xi_l0+0.5*xi_l1) * rah
137  * confo2.val_point(rah, M_PI/2., phi0+0.5*dp)
138  + 2. * (xi_p0+0.5*xi_p1)
139  * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
140  g2 = -2. * (xi_t0+0.5*xi_t1)
141  * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
142  h2 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+0.5*dp))
143  * (xi_t0+0.5*xi_t1) / rah
144  / confo2.val_point(rah, M_PI/2., phi0+0.5*dp) ;
145 
146  xi_t2 = dp * f2 ;
147  xi_p2 = dp * g2 ;
148  xi_l2 = dp * h2 ;
149 
150  // Third
151  f3 = - (xi_l0+0.5*xi_l2) * rah
152  * confo2.val_point(rah, M_PI/2., phi0+0.5*dp)
153  + 2. * (xi_p0+0.5*xi_p2)
154  * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
155  g3 = -2. * (xi_t0+0.5*xi_t2)
156  * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
157  h3 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+0.5*dp))
158  * (xi_t0+0.5*xi_t2) / rah
159  / confo2.val_point(rah, M_PI/2., phi0+0.5*dp) ;
160 
161  xi_t3 = dp * f3 ;
162  xi_p3 = dp * g3 ;
163  xi_l3 = dp * h3 ;
164 
165  // Forth
166  f4 = - (xi_l0+xi_l3) * rah * confo2.val_point(rah, M_PI/2., phi0+dp)
167  + 2. * (xi_p0+xi_p3) * dlnconfo.val_point(rah, M_PI/2., phi0+dp) ;
168  g4 = -2. * (xi_t0+xi_t3) * dlnconfo.val_point(rah, M_PI/2., phi0+dp) ;
169  h4 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+dp))
170  * (xi_t0+xi_t3) / rah / confo2.val_point(rah, M_PI/2., phi0+dp) ;
171 
172  xi_t4 = dp * f4 ;
173  xi_p4 = dp * g4 ;
174  xi_l4 = dp * h4 ;
175 
176  // Final results
177  // -------------
178  xi_tf = xi_t0 + (xi_t1 + 2.*xi_t2 + 2.*xi_t3 + xi_t4) / 6. ;
179  xi_pf = xi_p0 + (xi_p1 + 2.*xi_p2 + 2.*xi_p3 + xi_p4) / 6. ;
180  xi_lf = xi_l0 + (xi_l1 + 2.*xi_l2 + 2.*xi_l3 + xi_l4) / 6. ;
181 
182  // Final results are put into the initial data
183  // in order for the next step
184  // -------------------------------------------
185  xi_t0 = xi_tf ;
186  xi_p0 = xi_pf ;
187  xi_l0 = xi_lf ;
188 
189  } // End of the loop
190 
191  xi_f.set(0) = xi_tf ;
192  xi_f.set(1) = xi_pf ;
193  xi_f.set(2) = xi_lf ;
194 
195  }
196 
197  return xi_f ;
198 
199 }
200 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
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
Standard units of space, time and mass.
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
virtual double rad_ah() const
Radius of the apparent horizon.
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition: scalar.C:896
Multi-domain grid.
Definition: grilles.h:279
Basic array class.
Definition: tbl.h:164
Tbl runge_kutta_phi_bh(const Tbl &xi_i, const double &phi_i, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the phi direction for the solution of the Killing ve...