LORENE
binhor_coal.C
1 /*
2  * Copyright (c) 2004-2005 Francois Limousin
3  * Jose-Luis Jaramillo
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 char binhor_coal_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $" ;
25 
26 /*
27  * $Id: binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $
28  * $Log: binhor_coal.C,v $
29  * Revision 1.15 2014/10/13 08:52:42 j_novak
30  * Lorene classes and functions now belong to the namespace Lorene.
31  *
32  * Revision 1.14 2014/10/06 15:13:01 j_novak
33  * Modified #include directives to use c++ syntax.
34  *
35  * Revision 1.13 2007/04/13 15:28:55 f_limousin
36  * Lots of improvements, generalisation to an arbitrary state of
37  * rotation, implementation of the spatial metric given by Samaya.
38  *
39  * Revision 1.12 2006/08/01 14:37:19 f_limousin
40  * New version
41  *
42  * Revision 1.11 2006/06/28 13:36:09 f_limousin
43  * Convergence to a given irreductible mass
44  *
45  * Revision 1.10 2006/05/24 16:56:37 f_limousin
46  * Many small modifs.
47  *
48  * Revision 1.9 2005/09/13 18:33:15 f_limousin
49  * New function vv_bound_cart_bin(double) for computing binaries with
50  * berlin condition for the shift vector.
51  * Suppress all the symy and asymy in the importations.
52  *
53  * Revision 1.8 2005/07/11 08:21:57 f_limousin
54  * Implementation of a new boundary condition for the lapse in the binary
55  * case : boundary_nn_Dir_lapl().
56  *
57  * Revision 1.7 2005/03/10 17:21:52 f_limousin
58  * Add the Berlin boundary condition for the shift.
59  * Some changes to avoid warnings.
60  *
61  * Revision 1.6 2005/03/10 17:09:05 f_limousin
62  * Display the logarithm of viriel and convergence.
63  *
64  * Revision 1.5 2005/03/10 16:57:00 f_limousin
65  * Improve the convergence of the code coal_bh.
66  *
67  * Revision 1.4 2005/02/24 17:24:26 f_limousin
68  * The boundary conditions for psi, N and beta are now parameters in
69  * par_init.d and par_coal.d.
70  *
71  * Revision 1.3 2005/02/07 10:43:36 f_limousin
72  * Add the printing of the regularisation of the shift in the case N=0
73  * on the horizon.
74  *
75  * Revision 1.2 2004/12/31 15:40:21 f_limousin
76  * Improve the initialisation of several quantities in set_statiques().
77  *
78  * Revision 1.1 2004/12/29 16:11:19 f_limousin
79  * First version
80  *
81  *
82  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $
83  *
84  */
85 
86 //standard
87 #include <cstdlib>
88 
89 // Lorene
90 #include "tensor.h"
91 #include "isol_hor.h"
92 #include "graphique.h"
93 
94 
95 namespace Lorene {
96 void Bin_hor::set_statiques (double precis, double relax, int bound_nn,
97  double lim_nn, int bound_psi) {
98 
99  int nz = hole1.mp.get_mg()->get_nzone() ;
100 
101  set_omega(0) ;
102  hole1.init_met_trK() ;
103  hole2.init_met_trK() ;
104  init_bin_hor() ;
106 
107  int indic = 1 ;
108  int conte = 0 ;
109 
110  cout << "Static black holes : " << endl ;
111  while (indic == 1) {
112  Scalar lapse_un_old (hole1.n_auto) ;
113 
114  solve_psi (precis, relax, bound_psi) ;
115  solve_lapse (precis, relax, bound_nn, lim_nn) ;
116 
117  // des_profile(hole1.nn(), 0, 20, M_PI/2, M_PI) ;
118 
119  double erreur = 0 ;
120  Tbl diff (diffrelmax (lapse_un_old, hole1.n_auto)) ;
121  for (int i=1 ; i<nz ; i++)
122  if (diff(i) > erreur)
123  erreur = diff(i) ;
124 
125  cout << "Step : " << conte << " Difference : " << erreur << endl ;
126 
127  if (erreur < precis)
128  indic = -1 ;
129  conte ++ ;
130  }
131 }
132 
133 double Bin_hor::coal (double angu_vel, double relax, int nb_ome,
134  int nb_it, int bound_nn, double lim_nn,
135  int bound_psi, int bound_beta, double omega_eff,
136  double alpha,
137  ostream& fich_iteration, ostream& fich_correction,
138  ostream& fich_viriel, ostream& fich_kss,
139  int step, int search_mass, double mass_irr,
140  const int sortie) {
141 
142  int nz = hole1.mp.get_mg()->get_nzone() ;
143 
144  double precis = 1e-7 ;
145 
146  // LOOP INCREASING OMEGA :
147  cout << "OMEGA INCREASED STEP BY STEP." << endl ;
148  double homme = get_omega() ;
149  double inc_homme = (angu_vel - homme)/nb_ome ;
150  for (int pas = 0 ; pas <nb_ome ; pas ++) {
151 
152  bool verif = false ;
153  if (omega_eff == alpha*homme ) verif = true ;
154 
155  homme += inc_homme ;
156  set_omega (homme) ;
157  if (verif)
158  omega_eff = alpha*homme ;
159  Scalar beta_un_old (hole1.beta_auto(1)) ;
160 
161  solve_shift (precis, relax, bound_beta, omega_eff) ;
163 
164  solve_psi (precis, relax, bound_psi) ;
165  solve_lapse (precis, relax, bound_nn, lim_nn) ;
166 
167  // Convergence to the given irreductible mass
168  if (search_mass == 1 && step >= 30) {
169  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
170  sqrt(hole2.area_hor()/16/M_PI) ;
171  double error_m = (mass_irr - mass_area) / mass_irr ;
172  double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
173  hole1.mp.homothetie_interne(scaling_r) ;
174  hole1.radius = hole1.radius *scaling_r ;
175  hole2.mp.homothetie_interne(scaling_r) ;
176  hole2.radius = hole2.radius *scaling_r ;
177 
178  // Update of the different metrics (another possibility would
179  // be to set all derived quantities to 0x0, especially
180  // the connection p_connect
185 
186  }
187 
188  cout << "Angular momentum computed at the horizon : " << ang_mom_hor()
189  << endl ;
190 
191  double erreur = 0 ;
192  Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
193  for (int i=1 ; i<nz ; i++)
194  if (diff(i) > erreur)
195  erreur = diff(i) ;
196 
197  // Saving ok K_{ij}s^is^j
198  // -----------------------
199 
200  Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
202  hole1.get_gam().radial_vect(), 0, 1)) ;
203  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
204  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
205  int nnp = hole1.mp.get_mg()->get_np(1) ;
206  int nnt = hole2.mp.get_mg()->get_nt(1) ;
207  for (int k=0 ; k<nnp ; k++)
208  for (int j=0 ; j<nnt ; j++){
209  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
210  max_kss = kkss.val_grid_point(1, k, j, 0) ;
211  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
212  min_kss = kkss.val_grid_point(1, k, j, 0) ;
213  }
214 
215  if (sortie != 0) {
216  fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
217  fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
218  // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
219  fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
220  fich_kss << step << " " << max_kss << " " << min_kss << endl ;
221  }
222 
223  cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
224  step ++ ;
225  }
226 
227  // LOOP WITH FIXED OMEGA :
228 
229  if (nb_it !=0)
230  cout << "OMEGA FIXED" << endl ;
231  double erreur ;
232 
233  for (int pas = 0 ; pas <nb_it ; pas ++) {
234 
235  Scalar beta_un_old (hole1.beta_auto(1)) ;
236 
237  solve_shift (precis, relax, bound_beta, omega_eff) ;
239 
240  solve_psi (precis, relax, bound_psi) ;
241  solve_lapse (precis, relax, bound_nn, lim_nn) ;
242 
243  // Convergence to the given irreductible mass
244  if (search_mass == 1 && step >= 30) {
245  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
246  sqrt(hole2.area_hor()/16/M_PI) ;
247  double error_m = (mass_irr - mass_area) / mass_irr ;
248  double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
249 
250  hole1.mp.homothetie_interne(scaling_r) ;
251  hole1.radius = hole1.radius *scaling_r ;
252  hole2.mp.homothetie_interne(scaling_r) ;
253  hole2.radius = hole2.radius *scaling_r ;
254  }
255 
256  erreur = 0 ;
257  Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
258  for (int i=1 ; i<nz ; i++)
259  if (diff(i) > erreur)
260  erreur = diff(i) ;
261 
262  // Saving ok K_{ij}s^is^j
263  // -----------------------
264 
265  Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
267  hole1.get_gam().radial_vect(), 0, 1)) ;
268  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
269  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
270  int nnp = hole1.mp.get_mg()->get_np(1) ;
271  int nnt = hole2.mp.get_mg()->get_nt(1) ;
272  for (int k=0 ; k<nnp ; k++)
273  for (int j=0 ; j<nnt ; j++){
274  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
275  max_kss = kkss.val_grid_point(1, k, j, 0) ;
276  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
277  min_kss = kkss.val_grid_point(1, k, j, 0) ;
278  }
279 
280 
281  if (sortie != 0) {
282  fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
283  fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
284  // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
285  fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
286  fich_kss << step << " " << max_kss << " " << min_kss << endl ;
287  }
288 
289  cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
290  step ++ ;
291  }
292 
293  if (nb_it != 0){
294  fich_iteration << "#----------------------------" << endl ;
295  fich_correction << "#-----------------------------" << endl ;
296  fich_viriel << "#------------------------------" << endl ;
297  }
298 
299  return viriel() ;
300 }
301 }
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
const Metric & get_gam() const
metric
Definition: single_hor.C:339
double omega_hor() const
Orbital velocity.
Definition: single_param.C:160
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition: binhor_kij.C:88
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
double area_hor() const
Area of the horizon.
Definition: single_param.C:88
const Sym_tensor & get_k_dd() const
k_dd
Definition: single_hor.C:348
Lorene prototypes.
Definition: app_hor.h:64
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
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:362
void set_statiques(double precis, double relax, int bound_nn, double lim_nn, int bound_psi)
Initialize the systeme to Misner Lindquist solution, that is solving for N and in the case ...
Definition: binhor_coal.C:96
double viriel() const
Computes the viriel error, that is the difference between the ADM and the Komar masses, calculated by the asymptotic behaviours of respectively and N .
Definition: binhor_viriel.C:71
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
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition: isol_hor.h:1409
double coal(double ang_vel, double relax, int nb_om, int nb_it, int bound_nn, double lim_nn, int bound_psi, int bound_beta, double omega_eff, double alpha, ostream &fich_iteration, ostream &fich_correction, ostream &fich_viriel, ostream &fich_kss, int step, int search_mass, double mass_irr, const int sortie=0)
Solves the equation for a particular angular velocity, the systeme being initialized to Misner-Lindqu...
Definition: binhor_coal.C:133
void init_bin_hor()
Initialisation of the system.
Definition: bin_hor.C:161
double regul
Intensity of the correction on the shift vector.
Definition: isol_hor.h:912
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Scalar n_auto
Lapse function .
Definition: isol_hor.h:915
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ...
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
Definition: map_af.C:614
void init_met_trK()
Sets the 3-metric tilde to the flat metric and gamt_point, trK and trK_point to zero.
Definition: single_hor.C:536
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:906
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
double get_omega() const
Returns the angular velocity.
Definition: isol_hor.h:1422
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:321
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539