LORENE
bhole_equations_bin.C
1 /*
2  * Copyright (c) 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 char bhole_equations_bin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $" ;
24 
25 /*
26  * $Id: bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
27  * $Log: bhole_equations_bin.C,v $
28  * Revision 1.5 2014/10/13 08:52:40 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.4 2014/10/06 15:12:58 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.3 2006/04/27 09:12:32 p_grandclement
35  * First try at irrotational black holes
36  *
37  * Revision 1.2 2002/10/16 14:36:33 j_novak
38  * Reorganization of #include instructions of standard C++, in order to
39  * use experimental version 3 of gcc.
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.6 2001/05/07 09:11:56 phil
45  * *** empty log message ***
46  *
47  * Revision 2.5 2001/04/30 09:30:50 phil
48  * filtre pour poisson vectoriel
49  *
50  * Revision 2.4 2001/04/26 12:04:34 phil
51  * *** empty log message ***
52  *
53  * Revision 2.3 2001/04/25 15:08:48 phil
54  * vire fait_tkij
55  *
56  * Revision 2.2 2001/04/02 12:16:11 phil
57  * *** empty log message ***
58  *
59  * Revision 2.1 2001/03/22 10:41:36 phil
60  * modification prototypage
61  *
62  * Revision 2.0 2001/03/01 08:18:04 phil
63  * *** empty log message ***
64  *
65  *
66  * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
67  *
68  */
69 
70 //standard
71 #include <cstdlib>
72 #include <cmath>
73 
74 // Lorene
75 #include "nbr_spx.h"
76 #include "tenseur.h"
77 #include "bhole.h"
78 #include "proto.h"
79 #include "utilitaires.h"
80 #include "graphique.h"
81 
82 // Resolution pour le lapse
83 namespace Lorene {
84 void Bhole_binaire::solve_lapse (double precision, double relax) {
85 
86  assert ((relax >0) && (relax<=1)) ;
87 
88  cout << "-----------------------------------------------" << endl ;
89  cout << "Resolution LAPSE" << endl ;
90 
91  Tenseur lapse_un_old (hole1.n_auto) ;
92  Tenseur lapse_deux_old (hole2.n_auto) ;
93 
95  Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
96 
98  Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
99 
100  // Les sources
101 
102  Cmp source_un
104  +hole1.n_tot()*pow(hole1.psi_tot(), 4.)*kk_un) ;
105  source_un.std_base_scal() ;
106 
107  Cmp source_deux
109  +hole2.n_tot()*pow(hole2.psi_tot(), 4.)*kk_deux) ;
110  source_deux.std_base_scal() ;
111 
112  //On resout
113  hole1.n_auto.set() = hole1.n_auto() - 1./2. ;
114  hole2.n_auto.set() = hole2.n_auto() - 1./2. ;
115 
116  dirichlet_binaire (source_un, source_deux, -1., -1., hole1.n_auto.set(),
117  hole2.n_auto.set(), 0, precision) ;
118 
119  hole1.n_auto.set() = hole1.n_auto() + 1./2. ;
120  hole2.n_auto.set() = hole2.n_auto() + 1./2. ;
121 
122  hole1.n_auto.set().raccord(1) ;
123  hole2.n_auto.set().raccord(1) ;
124 
125  // La relaxation :
126  hole1.n_auto.set() = relax*hole1.n_auto() + (1-relax)*lapse_un_old() ;
127  hole2.n_auto.set() = relax*hole2.n_auto() + (1-relax)*lapse_deux_old() ;
128 
131 }
132 
133 //Resolution sur Psi
134 void Bhole_binaire::solve_psi (double precision, double relax) {
135 
136  assert ((relax>0) && (relax<=1)) ;
137 
138  cout << "-----------------------------------------------" << endl ;
139  cout << "Resolution PSI" << endl ;
140 
141  Tenseur psi_un_old (hole1.psi_auto) ;
142  Tenseur psi_deux_old (hole2.psi_auto) ;
143 
145  Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
146 
148  Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
149 
150  // Les sources
151  Cmp source_un (-pow(hole1.psi_tot(), 5.)*kk_un/8.) ;
152  source_un.std_base_scal() ;
153 
154  Cmp source_deux (-pow(hole2.psi_tot(), 5.)*kk_deux/8.) ;
155  source_deux.std_base_scal() ;
156 
157  // Les valeurs limites :
158  int np_un = hole1.mp.get_mg()->get_np(1) ;
159  int nt_un = hole1.mp.get_mg()->get_nt(1) ;
160  Valeur lim_un (hole1.mp.get_mg()->get_angu()) ;
161  lim_un = 1 ;
162  for (int k=0 ; k<np_un ; k++)
163  for (int j=0 ; j<nt_un ; j++)
164  lim_un.set(0, k, j, 0) = -0.5/hole1.rayon*hole1.psi_tot()(1, k, j, 0) ;
165  lim_un.std_base_scal() ;
166 
167  int np_deux = hole2.mp.get_mg()->get_np(1) ;
168  int nt_deux = hole2.mp.get_mg()->get_nt(1) ;
169  Valeur lim_deux (hole2.mp.get_mg()->get_angu()) ;
170  lim_deux = 1 ;
171  for (int k=0 ; k<np_deux ; k++)
172  for (int j=0 ; j<nt_deux ; j++)
173  lim_deux.set(0, k, j, 0) = -0.5/hole2.rayon*hole2.psi_tot()(1, k, j, 0) ;
174  lim_deux.std_base_scal() ;
175 
176  //On resout
177  neumann_binaire (source_un, source_deux, lim_un, lim_deux,
178  hole1.psi_auto.set(), hole2.psi_auto.set(), 0, precision) ;
179 
180  hole1.psi_auto.set() = hole1.psi_auto() + 1./2. ;
181  hole2.psi_auto.set() = hole2.psi_auto() + 1./2. ;
182 
183  hole1.psi_auto.set().raccord(1) ;
184  hole2.psi_auto.set().raccord(1) ;
185 
186  // La relaxation :
187  hole1.psi_auto.set() = relax*hole1.psi_auto() + (1-relax)*psi_un_old() ;
188  hole2.psi_auto.set() = relax*hole2.psi_auto() + (1-relax)*psi_deux_old() ;
189 
192 }
193 
194 
195 // Resolution sur shift a omega bloque.
196 void Bhole_binaire::solve_shift (double precision, double relax) {
197 
198  cout << "------------------------------------------------" << endl ;
199  cout << "Resolution shift : Omega = " << omega << endl ;
200 
201  // On determine les sources
202  Tenseur source_un (
205  if (source_un.get_etat() == ETATZERO) {
206  source_un.set_etat_qcq() ;
207  for (int i=0 ; i<3 ; i++)
208  source_un.set(i).set_etat_zero() ;
209  }
210  source_un.set_std_base() ;
211 
212  Tenseur source_deux (
215  if (source_deux.get_etat() == ETATZERO) {
216  source_deux.set_etat_qcq() ;
217  for (int i=0 ; i<3 ; i++)
218  source_deux.set(i).set_etat_zero() ;
219  }
220  source_deux.set_std_base() ;
221 
222  // On filtre les hautes frequences.
223  for (int i=0 ; i<3 ; i++) {
224  if (source_un(i).get_etat() != ETATZERO)
225  source_un.set(i).filtre(4) ;
226  if (source_deux(i).get_etat() != ETATZERO)
227  source_deux.set(i).filtre(4) ;
228  }
229 
230  // Les alignemenents pour le signe des CL.
231  double orientation_un = hole1.mp.get_rot_phi() ;
232  assert ((orientation_un==0) || (orientation_un == M_PI)) ;
233 
234  double orientation_deux = hole2.mp.get_rot_phi() ;
235  assert ((orientation_deux==0) || (orientation_deux == M_PI)) ;
236 
237  int aligne_un = (orientation_un == 0) ? 1 : -1 ;
238  int aligne_deux = (orientation_deux == 0) ? 1 : -1 ;
239 
240  // On regarde si toutes les composantes sont nulles :
241  int ind1 = 0 ;
242  int ind2 = 0 ;
243  for (int i=0 ; i<3 ; i++) {
244  if (source_un(i).get_etat() == ETATQCQ)
245  ind1 = 1 ;
246  if (source_deux(i).get_etat() == ETATQCQ)
247  ind2 = 1 ;
248  }
249 
250  if (ind1==0)
251  source_un.set_etat_zero() ;
252  if (ind2==0)
253  source_deux.set_etat_zero() ;
254 
255  // On determine les Cl en fonction de omega :
256  int np_un = hole1.mp.get_mg()->get_np (1) ;
257  int nt_un = hole1.mp.get_mg()->get_nt (1) ;
258 
259  Mtbl xa_mtbl_un (source_un.get_mp()->get_mg()) ;
260  xa_mtbl_un.set_etat_qcq() ;
261  Mtbl ya_mtbl_un (source_un.get_mp()->get_mg()) ;
262  ya_mtbl_un.set_etat_qcq() ;
263 
264  xa_mtbl_un = source_un.get_mp()->xa ;
265  ya_mtbl_un = source_un.get_mp()->ya ;
266 
267  Mtbl x_mtbl_un (source_un.get_mp()->get_mg()) ;
268  x_mtbl_un.set_etat_qcq() ;
269  Mtbl y_mtbl_un (source_un.get_mp()->get_mg()) ;
270  y_mtbl_un.set_etat_qcq() ;
271 
272  x_mtbl_un = source_un.get_mp()->x ;
273  y_mtbl_un = source_un.get_mp()->y ;
274 
275  // Les bases
276  Base_val** bases_un = hole1.mp.get_mg()->std_base_vect_cart() ;
277  Base_val** bases_deux = hole2.mp.get_mg()->std_base_vect_cart() ;
278 
279  Valeur lim_x_un (*hole1.mp.get_mg()->get_angu()) ;
280  lim_x_un = 1 ; // Juste pour affecter dans espace des configs ;
281  lim_x_un.set_etat_c_qcq() ;
282  for (int k=0 ; k<np_un ; k++)
283  for (int j=0 ; j<nt_un ; j++)
284  lim_x_un.set(0, k, j, 0) = aligne_un*omega*ya_mtbl_un(0, 0, 0, 0) + aligne_un*hole1.omega_local*y_mtbl_un(1,k,j,0) ;
285  lim_x_un.base = *bases_un[0] ;
286 
287  Valeur lim_y_un (*hole1.mp.get_mg()->get_angu()) ;
288  lim_y_un = 1 ; // Juste pour affecter dans espace des configs ;
289  lim_y_un.set_etat_c_qcq() ;
290  for (int k=0 ; k<np_un ; k++)
291  for (int j=0 ; j<nt_un ; j++)
292  lim_y_un.set(0, k, j, 0) = -aligne_un*omega*xa_mtbl_un(0, 0, 0, 0) - aligne_un*hole1.omega_local*x_mtbl_un(1,k,j,0) ;;
293  lim_y_un.base = *bases_un[1] ;
294 
295  Valeur lim_z_un (*hole1.mp.get_mg()->get_angu()) ;
296  lim_z_un = 1 ;
297  for (int k=0 ; k<np_un ; k++)
298  for (int j=0 ; j<nt_un ; j++)
299  lim_z_un.set(0, k, j, 0) = 0 ;
300  lim_z_un.base = *bases_un[2] ;
301 
302  // On determine les Cl en fonction de omega :
303  int np_deux = hole2.mp.get_mg()->get_np (1) ;
304  int nt_deux = hole2.mp.get_mg()->get_nt (1) ;
305 
306  Mtbl xa_mtbl_deux (source_deux.get_mp()->get_mg()) ;
307  xa_mtbl_deux.set_etat_qcq() ;
308  Mtbl ya_mtbl_deux (source_deux.get_mp()->get_mg()) ;
309  ya_mtbl_deux.set_etat_qcq() ;
310 
311  xa_mtbl_deux = source_deux.get_mp()->xa ;
312  ya_mtbl_deux = source_deux.get_mp()->ya ;
313 
314  Mtbl x_mtbl_deux (source_deux.get_mp()->get_mg()) ;
315  x_mtbl_deux.set_etat_qcq() ;
316  Mtbl y_mtbl_deux (source_deux.get_mp()->get_mg()) ;
317  y_mtbl_deux.set_etat_qcq() ;
318 
319  x_mtbl_deux = source_deux.get_mp()->x ;
320  y_mtbl_deux = source_deux.get_mp()->y ;
321 
322  Valeur lim_x_deux (*hole2.mp.get_mg()->get_angu()) ;
323  lim_x_deux = 1 ; // Juste pour affecter dans espace des configs ;
324  lim_x_deux.set_etat_c_qcq() ;
325  for (int k=0 ; k<np_deux ; k++)
326  for (int j=0 ; j<nt_deux ; j++)
327  lim_x_deux.set(0, k, j, 0) = aligne_deux*omega*ya_mtbl_deux(0, 0, 0, 0) + aligne_deux*hole2.omega_local*y_mtbl_deux(1,k,j,0) ;
328  lim_x_deux.base = *bases_deux[0] ;
329 
330  Valeur lim_y_deux (*hole2.mp.get_mg()->get_angu()) ;
331  lim_y_deux = 1 ; // Juste pour affecter dans espace des configs ;
332  lim_y_deux.set_etat_c_qcq() ;
333  for (int k=0 ; k<np_deux ; k++)
334  for (int j=0 ; j<nt_deux ; j++)
335  lim_y_deux.set(0, k, j, 0) = -aligne_deux*omega*xa_mtbl_deux(0, 0, 0, 0) - aligne_deux*hole2.omega_local*x_mtbl_deux(1,k,j,0) ;
336  lim_y_deux.base = *bases_deux[1] ;
337 
338  Valeur lim_z_deux (*hole2.mp.get_mg()->get_angu()) ;
339  lim_z_deux = 1 ;
340  for (int k=0 ; k<np_deux ; k++)
341  for (int j=0 ; j<nt_deux ; j++)
342  lim_z_deux.set(0, k, j, 0) = 0 ;
343  lim_z_deux.base = *bases_deux[2] ;
344 
345  for (int i=0 ; i<3 ; i++) {
346  delete bases_un[i] ;
347  delete bases_deux[i] ;
348  }
349  delete [] bases_un ;
350  delete [] bases_deux ;
351 
352  // On resout le truc :
353  Tenseur shift_un_old (hole1.shift_auto) ;
354  Tenseur shift_deux_old (hole2.shift_auto) ;
355 
356  poisson_vect_binaire (1./3., source_un, source_deux,
357  lim_x_un, lim_y_un, lim_z_un,
358  lim_x_deux, lim_y_deux, lim_z_deux,
359  hole1.shift_auto, hole2.shift_auto, 0, precision) ;
360 
361  for (int i=0 ; i<3 ; i++) {
362  hole1.shift_auto.set(i).raccord(1) ;
363  hole2.shift_auto.set(i).raccord(1) ;
364  }
365 
366  // On regularise les shift.
367  hole1.shift_auto = relax*hole1.shift_auto +
368  (1-relax)*shift_un_old ;
369  hole2.shift_auto = relax*hole2.shift_auto +
370  (1-relax)*shift_deux_old ;
371 
372  double diff_un = regle (hole1.shift_auto, hole2.shift_auto, omega, hole1.omega_local) ;
373  double diff_deux = regle (hole2.shift_auto, hole1.shift_auto, omega, hole2.omega_local) ;
374  hole1.regul = diff_un ;
375  hole2.regul = diff_deux ;
376 
377  cout << "Difference relatives : " << diff_un << " " << diff_deux << endl ;
378 }
379 }
Lorene::Bhole::tkij_tot
Tenseur_sym tkij_tot
Total .
Definition: bhole.h:308
Lorene::Bhole::mp
Map_af & mp
Affine mapping.
Definition: bhole.h:273
Lorene::Tenseur::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
Lorene::Base_val
Bases of the spectral expansions.
Definition: base_val.h:322
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Bhole::shift_auto
Tenseur shift_auto
Part of generated by the hole.
Definition: bhole.h:297
Lorene::Valeur::set_etat_c_qcq
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition: valeur.C:701
Lorene::Bhole::omega_local
double omega_local
local angular velocity
Definition: bhole.h:276
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Lorene::Bhole::rayon
double rayon
Radius of the horizon in LORENE's units.
Definition: bhole.h:274
Lorene::Map::get_rot_phi
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Lorene::Bhole::fait_psi_comp
void fait_psi_comp(const Bhole &comp)
Imports the part of due to the companion hole comp .
Definition: bhole.C:280
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
Lorene::Bhole_binaire::hole2
Bhole hole2
Black hole two.
Definition: bhole.h:763
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
Lorene::Mtbl
Multi-domain array.
Definition: mtbl.h:118
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Bhole_binaire::solve_lapse
void solve_lapse(double precis, double relax)
Solves the equation for the lapse~:
Definition: bhole_equations_bin.C:84
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Lorene::Bhole::n_tot
Tenseur n_tot
Total N .
Definition: bhole.h:288
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Bhole_binaire::hole1
Bhole hole1
Black hole one.
Definition: bhole.h:762
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
Lorene::Bhole_binaire::solve_psi
void solve_psi(double precis, double relax)
Solves the equation for the conformal factor~:
Definition: bhole_equations_bin.C:134
Lorene::Bhole::psi_tot
Tenseur psi_tot
Total .
Definition: bhole.h:292
Lorene::Bhole_binaire::solve_shift
void solve_shift(double precis, double relax)
Solves the equation for the shift, using the Oohara-Nakarmure scheme~:
Definition: bhole_equations_bin.C:196
Lorene::Bhole::n_auto
Tenseur n_auto
Part of N generated by the hole.
Definition: bhole.h:286
Lorene::Valeur::std_base_scal
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:824
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Cmp::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: cmp_manip.C:74
Lorene::flat_scalar_prod
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...
Definition: tenseur_operateur.C:653
Lorene::Bhole::fait_n_comp
void fait_n_comp(const Bhole &comp)
Imports the part of N due to the companion hole comp .
Definition: bhole.C:254
Lorene::Mg3d::std_base_vect_cart
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Definition: mg3d_std_base.C:177
Lorene::Bhole::tkij_auto
Tenseur_sym tkij_auto
Auto .
Definition: bhole.h:307
Lorene::Tenseur::get_mp
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:699
Lorene::Map::x
Coord x
x coordinate centered on the grid
Definition: map.h:726
Lorene::Map::xa
Coord xa
Absolute x coordinate.
Definition: map.h:730
Lorene::Bhole::regul
double regul
Intensity of the correction on the shift vector.
Definition: bhole.h:284
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Lorene::Bhole_binaire::omega
double omega
Position of the axis of rotation.
Definition: bhole.h:769
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363
Lorene::Map::y
Coord y
y coordinate centered on the grid
Definition: map.h:727
Lorene::Bhole::psi_auto
Tenseur psi_auto
Part of generated by the hole.
Definition: bhole.h:290
Lorene::Cmp::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: cmp_raccord.C:170
Lorene::Bhole::taij_tot
Tenseur_sym taij_tot
Total , which must be zero on the horizon of the regularisation on the shift has been done.
Definition: bhole.h:305
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene::Bhole::grad_n_tot
Tenseur grad_n_tot
Total gradient of N .
Definition: bhole.h:294
Lorene::Map::ya
Coord ya
Absolute y coordinate.
Definition: map.h:731