LORENE
bhole_pseudo_kerr.C
1 /*
2  * Copyright (c) 2000-2001 Philippe Grandclement
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 
24 
25 /*
26  * $Id: bhole_pseudo_kerr.C,v 1.8 2016/12/05 16:17:45 j_novak Exp $
27  * $Log: bhole_pseudo_kerr.C,v $
28  * Revision 1.8 2016/12/05 16:17:45 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.7 2014/10/13 08:52:40 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.6 2014/10/06 15:12:58 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.5 2008/08/19 06:41:59 j_novak
38  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
39  * cast-type operations, and constant strings that must be defined as const char*
40  *
41  * Revision 1.4 2004/03/25 10:28:57 j_novak
42  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
43  *
44  * Revision 1.3 2003/10/03 15:58:43 j_novak
45  * Cleaning of some headers
46  *
47  * Revision 1.2 2002/10/16 14:36:32 j_novak
48  * Reorganization of #include instructions of standard C++, in order to
49  * use experimental version 3 of gcc.
50  *
51  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
52  * LORENE
53  *
54  * Revision 2.9 2001/05/16 14:51:36 phil
55  * correction calcul de kij
56  *
57  * Revision 2.8 2001/05/07 12:24:30 phil
58  * correction calcul de J
59  *
60  * Revision 2.7 2001/05/07 09:28:33 phil
61  * *** empty log message ***
62  *
63  * Revision 2.6 2001/02/12 15:36:58 phil
64  * ajout calcul de J a l infini
65  *
66  * Revision 2.5 2000/12/14 14:11:54 phil
67  * correction cl sur psi
68  *
69  * Revision 2.4 2000/12/14 12:41:54 phil
70  * on met les bases dans les sources
71  *
72  * Revision 2.3 2000/12/14 10:57:47 phil
73  * corections diverses et sans importances
74  *
75  * Revision 2.2 2000/12/14 10:45:00 phil
76  * ATTENTION : PASSAGE DE PHI A PSI
77  *
78  * Revision 2.1 2000/11/03 12:56:33 phil
79  * ajout de const
80  *
81  * Revision 2.0 2000/10/20 09:19:09 phil
82  * *** empty log message ***
83  *
84  *
85  * $Header: /cvsroot/Lorene/C++/Source/Bhole/bhole_pseudo_kerr.C,v 1.8 2016/12/05 16:17:45 j_novak Exp $
86  *
87  */
88 
89 //standard
90 #include <cstdlib>
91 #include <cmath>
92 
93 // Lorene
94 #include "tenseur.h"
95 #include "bhole.h"
96 #include "proto.h"
97 
98 //Resolution pour le lapse pour 1 seul trou
99 namespace Lorene {
100 void Bhole::solve_lapse_seul (double relax) {
101 
102  assert ((relax>0) && (relax<=1)) ;
103 
104  cout << "Resolution LAPSE" << endl ;
105 
106  // Pour la relaxation ...
107  Cmp lapse_old (n_auto()) ;
109  Tenseur kk (mp) ;
110  kk = 0 ;
111  Tenseur work(mp) ;
112  work.set_etat_qcq() ;
113  for (int i=0 ; i<3 ; i++) {
114  work.set() = auxi(i, i) ;
115  kk = kk + work ;
116  }
117 
118  // La source
119  Cmp source
121  +pow(psi_auto(), 4.)*n_auto()*kk()) ;
122  source.std_base_scal() ;
123 
124  // On resout pour N-1 :
125  Valeur limite (mp.get_mg()->get_angu()) ;
126  limite = -1 ;
127  limite.std_base_scal() ;
128 
129  Cmp soluce (source.poisson_dirichlet(limite, 0)) ;
130  soluce = soluce + 1 ; // Permet de trouver N
131  soluce.raccord(1) ;
132 
133  n_auto.set() = relax*soluce + (1-relax)*lapse_old ;
134 }
135 
136 
137 // Resolution sur Psi :
138 void Bhole::solve_psi_seul (double relax) {
139 
140  assert ((relax>0) && (relax<=1)) ;
141 
142  cout << "Resolution PSI" << endl ;
143 
144  Cmp psi_old (psi_auto()) ;
146  Tenseur kk (mp) ;
147  kk = 0 ;
148  Tenseur work(mp) ;
149  work.set_etat_qcq() ;
150  for (int i=0 ; i<3 ; i++) {
151  work.set() = auxi(i, i) ;
152  kk = kk + work ;
153  }
154 
155  // La source :
156  Cmp source (-pow(psi_auto(), 5.)*kk()/8.) ;
157  source.std_base_scal() ;
158 
159  // La condition limite de type neumann :
160  int np = mp.get_mg()->get_np(1) ;
161  int nt = mp.get_mg()->get_nt(1) ;
162  Valeur limite (mp.get_mg()->get_angu()) ;
163  limite = 1 ;
164  for (int k=0 ; k<np ; k++)
165  for (int j=0 ; j<nt ; j++)
166  limite.set(0, k, j, 0) = -0.5/rayon*psi_auto()(1, k, j, 0) ;
167  limite.std_base_scal() ;
168 
169  Cmp soluce (source.poisson_neumann(limite, 0)) ;
170  soluce = soluce + 1 ;
171  soluce.raccord(1) ;
172 
173  psi_auto.set() = relax*soluce + (1-relax)*psi_old ;
174 
175 }
176 
177 
178 // Le shift. Processus iteratif pour cause de CL.
179 void Bhole::solve_shift_seul (double precision, double relax) {
180 
181  assert (precision > 0) ;
182  assert ((relax>0) && (relax<=1)) ;
183 
184  cout << "resolution SHIFT" << endl ;
185 
186  Tenseur shift_old (shift_auto) ;
187 
190  source.set_std_base() ;
191 
192  // On verifie si les 3 composantes ne sont pas nulles :
193  if (source.get_etat() == ETATQCQ) {
194  int indic = 0 ;
195  for (int i=0 ; i<3 ; i++)
196  if (source(i).get_etat() == ETATQCQ)
197  indic = 1 ;
198  if (indic ==0)
199  source.set_etat_zero() ;
200  }
201 
202  // On filtre les hautes frequences pour raison de stabilite :
203  if (source.get_etat() == ETATQCQ)
204  for (int i=0 ; i<3 ; i++)
205  source.set(i).filtre(4) ;
206 
207 
208  // On determine les conditions limites en fonction du omega et du boost :
209  int np = mp.get_mg()->get_np(1) ;
210  int nt = mp.get_mg()->get_nt(1) ;
211 
212  Mtbl x_mtbl (mp.get_mg()) ;
213  x_mtbl.set_etat_qcq() ;
214  Mtbl y_mtbl (mp.get_mg()) ;
215  y_mtbl.set_etat_qcq() ;
216  x_mtbl = mp.x ;
217  y_mtbl = mp.y ;
218 
219  // Les bases pour les conditions limites :
220  Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
221 
222  Valeur lim_x (mp.get_mg()->get_angu()) ;
223  lim_x = 1 ;
224  for (int k=0 ; k<np ; k++)
225  for (int j=0 ; j<nt ; j++)
226  lim_x.set(0, k, j, 0) = omega*y_mtbl(1, k, j, 0)-boost[0] ;
227  lim_x.base = *bases[0] ;
228 
229  Valeur lim_y (mp.get_mg()->get_angu()) ;
230  lim_y = 1 ;
231  for (int k=0 ; k<np ; k++)
232  for (int j=0 ; j<nt ; j++)
233  lim_y.set(0, k, j, 0) = - omega*x_mtbl(1, k, j, 0)-boost[1] ;
234  lim_y.base = *bases[1] ;
235 
236  Valeur lim_z (mp.get_mg()->get_angu()) ;
237  lim_z = 1 ;
238  for (int k=0 ; k<np ; k++)
239  for (int j=0 ; j<nt ; j++)
240  lim_z.set(0, k, j, 0) = -boost[2] ;
241  lim_z.base = *bases[2] ;
242 
243  // On n'en a plus besoin
244  for (int i=0 ; i<3 ; i++)
245  delete bases[i] ;
246  delete [] bases ;
247 
248  // On resout :
249  poisson_vect_frontiere(1./3., source, shift_auto, lim_x, lim_y,
250  lim_z, 0, precision, 20) ;
251 
252  shift_auto = relax*shift_auto + (1-relax)*shift_old ;
253 }
254 
255 
256 // La regularisation si un seul trou noir :
258 
259  // Le vecteur B (non tournant et non boostant)
260  Tenseur tbi (shift_auto) ;
261  for (int i=0 ; i<3 ; i++) {
262  tbi.set(i).va.coef_i() ;
263  tbi.set(i).va.set_etat_c_qcq() ;
264  }
265 
266  for (int i=0 ; i<3 ; i++)
267  shift_auto(i).va.coef_i() ;
268 
269  tbi.set(0) = *shift_auto(0).va.c - omega*shift_auto.get_mp()->y + boost[0];
270  tbi.set(1) = *shift_auto(1).va.c + omega*shift_auto.get_mp()->x + boost[1];
271  tbi.set(2) = *shift_auto(2).va.c + boost[2];
272  tbi.set_std_base() ;
273 
274  // On evite soucis a l'infini (on a besoin que de la valeur sur horizon
275  tbi.set(0).annule(mp.get_mg()->get_nzone()-1) ;
276  tbi.set(1).annule(mp.get_mg()->get_nzone()-1) ;
277 
278  Tenseur derive_r (mp, 1, CON, mp.get_bvect_cart()) ;
279  derive_r.set_etat_qcq() ;
280  for (int i=0 ; i<3 ; i++)
281  derive_r.set(i) = tbi(i).dsdr() ;
282 
283  Valeur val_hor (mp.get_mg()) ;
284  Valeur fonction_radiale (mp.get_mg()) ;
285  Cmp enleve (mp) ;
286 
287  double erreur = 0 ;
288  int nz = mp.get_mg()->get_nzone() ;
289  int np = mp.get_mg()->get_np(1) ;
290  int nt = mp.get_mg()->get_nt(1) ;
291  int nr = mp.get_mg()->get_nr(1) ;
292 
293  double r_0 = mp.val_r(1, -1, 0, 0) ;
294  double r_1 = mp.val_r(1, 1, 0, 0) ;
295 
296  for (int comp=0 ; comp<3 ; comp++) {
297  val_hor.annule_hard() ;
298  for (int k=0 ; k<np ; k++)
299  for (int j=0 ; j<nt ; j++)
300  for (int i=0 ; i<nr ; i++)
301  val_hor.set(1, k, j, i) = derive_r(comp)(1, k, j, 0) ;
302 
303  fonction_radiale = pow(r_1-mp.r, 3.)* (mp.r-r_0)/pow(r_1-r_0, 3.) ;
304  fonction_radiale.annule(0) ;
305  fonction_radiale.annule(2, nz-1) ;
306 
307  enleve = fonction_radiale*val_hor ;
308  enleve.va.base = shift_auto(comp).va.base ;
309 
310  // Ca devrai annuler la derivee de B sur H et donc rendre K regulier
311  Cmp copie (shift_auto(comp)) ;
312  shift_auto.set(comp) = shift_auto(comp)-enleve ;
313 
314  assert (shift_auto(comp).check_dzpuis(0)) ;
315 
316  // On regarde l'intensite de la correction si non nul !
317  Tbl norm (norme(shift_auto(comp))) ;
318  if (norm(1) > 1e-5) {
319  Tbl diff (diffrelmax (copie, shift_auto(comp))) ;
320  if (erreur<diff(1))
321  erreur = diff(1) ;
322  }
323  }
324  regul = erreur ;
325 }
326 
327 // On calcul Kij sachant que la regulatisation sur le shift doit
328 // etre faite avant.
330 
331  fait_taij_auto() ;
332 
333  Cmp lapse_non_sing (division_xpun(n_auto(), 0)) ;
334  Cmp auxi (mp) ;
335 
337  for (int i=0 ; i<3 ; i++)
338  for (int j=i ; j<3 ; j++) {
339  auxi = taij_auto(i, j) ;
340  auxi = division_xpun (auxi, 0) ;
341  tkij_auto.set(i, j) = auxi/lapse_non_sing/2. ;
342  tkij_auto.set(i, j).raccord(1) ;
343  }
344 }
345 
346 // Calcul masse ADM via valeur asymptotiques
347 double Bhole::masse_adm_seul () const {
348 
349  Cmp integrant (psi_auto().dsdr()) ;
350  double masse = mp.integrale_surface_infini (integrant) ;
351  masse /= -2*M_PI ;
352  return masse ;
353 }
354 
355 double Bhole::masse_komar_seul() const {
356 
357  Cmp integrant (n_auto().dsdr()) ;
358  double masse = mp.integrale_surface_infini (integrant) ;
359  masse /= 4*M_PI ;
360  return masse ;
361 }
362 
363 // Calcul du moment angulaire via integrale a l infini
364 // Non valable si le boost != 0 ;
365 
366 double Bhole::moment_seul_inf() const {
367 
368  // On verifie si le boost est bien nul :
369  double indic = 0 ;
370  for (int i=0 ; i<3 ; i++)
371  if (boost[i] != 0)
372  indic = 1 ;
373  if (indic == 1) {
374  cout << "Calcul du moment non valable pour un boost != 0" << endl ;
375  abort() ;
376  }
377 
378  if (omega == 0)
379  return 0 ;
380  else {
381  Tenseur vecteur_un (mp, 1, CON, mp.get_bvect_cart()) ;
382  vecteur_un.set_etat_qcq() ;
383  for (int i=0 ; i<3 ; i++)
384  vecteur_un.set(i) = tkij_auto(0, i) ;
385  vecteur_un.change_triad (mp.get_bvect_spher()) ;
386  Cmp integrant_un (vecteur_un(0)) ;
387  integrant_un.mult_r_zec() ;
388 
389  Tenseur vecteur_deux (mp, 1, CON, mp.get_bvect_cart()) ;
390  vecteur_deux.set_etat_qcq() ;
391  for (int i=0 ; i<3 ; i++)
392  vecteur_deux.set(i) = tkij_auto(1, i) ;
393  vecteur_deux.change_triad (mp.get_bvect_spher()) ;
394  Cmp integrant_deux (vecteur_deux(0)) ;
395  integrant_deux.mult_r_zec() ;
396 
397  // On multiplie par y et x :
398  integrant_un.va = integrant_un.va.mult_st() ;
399  integrant_un.va = integrant_un.va.mult_sp() ;
400 
401  integrant_deux.va = integrant_deux.va.mult_st() ;
402  integrant_deux.va = integrant_deux.va.mult_cp() ;
403 
404  double moment = mp.integrale_surface_infini (-integrant_un+integrant_deux) ;
405  moment /= 8*M_PI ;
406  return moment ;
407  }
408 }
409 
410 // Calcul du moment angulaire via integrale sur l horizon
411 // Non valable si le boost != 0 ;
412 
413 double Bhole::moment_seul_hor() const {
414 
415  // On verifie si le boost est bien nul :
416  double indic = 0 ;
417  for (int i=0 ; i<3 ; i++)
418  if (boost[i] != 0)
419  indic = 1 ;
420  if (indic == 1) {
421  cout << "Calcul du moment non valable pour un boost != 0" << endl ;
422  abort() ;
423  }
424 
425  if (omega == 0)
426  return 0 ;
427  else {
428  // Integrale sur l horizon :
429  Cmp xa (mp) ;
430  xa = mp.xa ;
431  xa.std_base_scal() ;
432 
433  Cmp ya (mp) ;
434  ya = mp.ya ;
435  ya.std_base_scal() ;
436 
437  Tenseur vecteur (mp, 1, CON, mp.get_bvect_cart()) ;
438  vecteur.set_etat_qcq() ;
439  for (int i=0 ; i<3 ; i++)
440  vecteur.set(i) = (-ya*tkij_auto(0, i)+xa * tkij_auto(1, i)) ;
441  vecteur.set_std_base() ;
442  vecteur.annule(mp.get_mg()->get_nzone()-1) ;
443  vecteur.change_triad (mp.get_bvect_spher()) ;
444 
445  Cmp integrant (pow(psi_auto(), 6)*vecteur(0)) ;
446  integrant.std_base_scal() ;
447  double moment = mp.integrale_surface (integrant, rayon)/8/M_PI ;
448  return moment ;
449  }
450 }
451 
452 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
double moment_seul_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and H de...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Cmp ** c
The components.
Definition: tenseur.h:325
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
const Valeur & mult_sp() const
Returns applied to *this.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
Multi-domain array.
Definition: mtbl.h:118
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
void coef_i() const
Computes the physical value of *this.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Definition: cmp_r_manip.C:106
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Definition: map_af_radius.C:99
void solve_psi_seul(double relax)
Solves the equation for ~: with the condition that on the horizon, where f is the value of on the ...
void solve_lapse_seul(double relax)
Solves the equation for N ~: with the condition that N =0 on the horizon.
Tenseur psi_auto
Part of generated by the hole.
Definition: bhole.h:290
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
void fait_taij_auto()
Calculates the part of generated by shift_auto .
Definition: bhole.C:382
Tenseur_sym taij_auto
Part of generated by the hole.
Definition: bhole.h:299
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:827
Tenseur shift_auto
Part of generated by the hole.
Definition: bhole.h:297
void regularise_seul()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon...
void solve_shift_seul(double precis, double relax)
Solves the equation for ~: with on the horizon, being the boost and .
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:702
void fait_tkij()
Calculates the total .
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: cmp_manip.C:77
Map_af & mp
Affine mapping.
Definition: bhole.h:273
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Valeur & mult_st() const
Returns applied to *this.
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
double rayon
Radius of the horizon in LORENE&#39;s units.
Definition: bhole.h:274
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur_sym tkij_auto
Auto .
Definition: bhole.h:307
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
double * boost
Lapse on the horizon.
Definition: bhole.h:283
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Coord ya
Absolute y coordinate.
Definition: map.h:749
Bases of the spectral expansions.
Definition: base_val.h:325
const Valeur & mult_cp() const
Returns applied to *this.
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
Coord y
y coordinate centered on the grid
Definition: map.h:745
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: cmp_raccord.C:173
double masse_komar_seul() const
Calculates the Komar mass of the black hole using : .
Coord x
x coordinate centered on the grid
Definition: map.h:744
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
double omega
Angular velocity in LORENE&#39;s units.
Definition: bhole.h:275
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:704
double moment_seul_inf() const
Calculates the angular momentum of the black hole using the formula at infinity : where ...
Tenseur n_auto
Part of N generated by the hole.
Definition: bhole.h:286
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
double regul
Intensity of the correction on the shift vector.
Definition: bhole.h:284
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558
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