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 
25 
26 /*
27  * $Id: binhor_coal.C,v 1.16 2016/12/05 16:17:46 j_novak Exp $
28  * $Log: binhor_coal.C,v $
29  * Revision 1.16 2016/12/05 16:17:46 j_novak
30  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31  *
32  * Revision 1.15 2014/10/13 08:52:42 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.14 2014/10/06 15:13:01 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.13 2007/04/13 15:28:55 f_limousin
39  * Lots of improvements, generalisation to an arbitrary state of
40  * rotation, implementation of the spatial metric given by Samaya.
41  *
42  * Revision 1.12 2006/08/01 14:37:19 f_limousin
43  * New version
44  *
45  * Revision 1.11 2006/06/28 13:36:09 f_limousin
46  * Convergence to a given irreductible mass
47  *
48  * Revision 1.10 2006/05/24 16:56:37 f_limousin
49  * Many small modifs.
50  *
51  * Revision 1.9 2005/09/13 18:33:15 f_limousin
52  * New function vv_bound_cart_bin(double) for computing binaries with
53  * berlin condition for the shift vector.
54  * Suppress all the symy and asymy in the importations.
55  *
56  * Revision 1.8 2005/07/11 08:21:57 f_limousin
57  * Implementation of a new boundary condition for the lapse in the binary
58  * case : boundary_nn_Dir_lapl().
59  *
60  * Revision 1.7 2005/03/10 17:21:52 f_limousin
61  * Add the Berlin boundary condition for the shift.
62  * Some changes to avoid warnings.
63  *
64  * Revision 1.6 2005/03/10 17:09:05 f_limousin
65  * Display the logarithm of viriel and convergence.
66  *
67  * Revision 1.5 2005/03/10 16:57:00 f_limousin
68  * Improve the convergence of the code coal_bh.
69  *
70  * Revision 1.4 2005/02/24 17:24:26 f_limousin
71  * The boundary conditions for psi, N and beta are now parameters in
72  * par_init.d and par_coal.d.
73  *
74  * Revision 1.3 2005/02/07 10:43:36 f_limousin
75  * Add the printing of the regularisation of the shift in the case N=0
76  * on the horizon.
77  *
78  * Revision 1.2 2004/12/31 15:40:21 f_limousin
79  * Improve the initialisation of several quantities in set_statiques().
80  *
81  * Revision 1.1 2004/12/29 16:11:19 f_limousin
82  * First version
83  *
84  *
85  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.16 2016/12/05 16:17:46 j_novak Exp $
86  *
87  */
88 
89 //standard
90 #include <cstdlib>
91 
92 // Lorene
93 #include "tensor.h"
94 #include "isol_hor.h"
95 #include "graphique.h"
96 
97 
98 namespace Lorene {
99 void Bin_hor::set_statiques (double precis, double relax, int bound_nn,
100  double lim_nn, int bound_psi) {
101 
102  int nz = hole1.mp.get_mg()->get_nzone() ;
103 
104  set_omega(0) ;
105  hole1.init_met_trK() ;
106  hole2.init_met_trK() ;
107  init_bin_hor() ;
109 
110  int indic = 1 ;
111  int conte = 0 ;
112 
113  cout << "Static black holes : " << endl ;
114  while (indic == 1) {
115  Scalar lapse_un_old (hole1.n_auto) ;
116 
117  solve_psi (precis, relax, bound_psi) ;
118  solve_lapse (precis, relax, bound_nn, lim_nn) ;
119 
120  // des_profile(hole1.nn(), 0, 20, M_PI/2, M_PI) ;
121 
122  double erreur = 0 ;
123  Tbl diff (diffrelmax (lapse_un_old, hole1.n_auto)) ;
124  for (int i=1 ; i<nz ; i++)
125  if (diff(i) > erreur)
126  erreur = diff(i) ;
127 
128  cout << "Step : " << conte << " Difference : " << erreur << endl ;
129 
130  if (erreur < precis)
131  indic = -1 ;
132  conte ++ ;
133  }
134 }
135 
136 double Bin_hor::coal (double angu_vel, double relax, int nb_ome,
137  int nb_it, int bound_nn, double lim_nn,
138  int bound_psi, int bound_beta, double omega_eff,
139  double alpha,
140  ostream& fich_iteration, ostream& fich_correction,
141  ostream& fich_viriel, ostream& fich_kss,
142  int step, int search_mass, double mass_irr,
143  const int sortie) {
144 
145  int nz = hole1.mp.get_mg()->get_nzone() ;
146 
147  double precis = 1e-7 ;
148 
149  // LOOP INCREASING OMEGA :
150  cout << "OMEGA INCREASED STEP BY STEP." << endl ;
151  double homme = get_omega() ;
152  double inc_homme = (angu_vel - homme)/nb_ome ;
153  for (int pas = 0 ; pas <nb_ome ; pas ++) {
154 
155  bool verif = false ;
156  if (omega_eff == alpha*homme ) verif = true ;
157 
158  homme += inc_homme ;
159  set_omega (homme) ;
160  if (verif)
161  omega_eff = alpha*homme ;
162  Scalar beta_un_old (hole1.beta_auto(1)) ;
163 
164  solve_shift (precis, relax, bound_beta, omega_eff) ;
166 
167  solve_psi (precis, relax, bound_psi) ;
168  solve_lapse (precis, relax, bound_nn, lim_nn) ;
169 
170  // Convergence to the given irreductible mass
171  if (search_mass == 1 && step >= 30) {
172  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
173  sqrt(hole2.area_hor()/16/M_PI) ;
174  double error_m = (mass_irr - mass_area) / mass_irr ;
175  double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
176  hole1.mp.homothetie_interne(scaling_r) ;
177  hole1.radius = hole1.radius *scaling_r ;
178  hole2.mp.homothetie_interne(scaling_r) ;
179  hole2.radius = hole2.radius *scaling_r ;
180 
181  // Update of the different metrics (another possibility would
182  // be to set all derived quantities to 0x0, especially
183  // the connection p_connect
188 
189  }
190 
191  cout << "Angular momentum computed at the horizon : " << ang_mom_hor()
192  << endl ;
193 
194  double erreur = 0 ;
195  Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
196  for (int i=1 ; i<nz ; i++)
197  if (diff(i) > erreur)
198  erreur = diff(i) ;
199 
200  // Saving ok K_{ij}s^is^j
201  // -----------------------
202 
203  Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
205  hole1.get_gam().radial_vect(), 0, 1)) ;
206  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
207  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
208  int nnp = hole1.mp.get_mg()->get_np(1) ;
209  int nnt = hole2.mp.get_mg()->get_nt(1) ;
210  for (int k=0 ; k<nnp ; k++)
211  for (int j=0 ; j<nnt ; j++){
212  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
213  max_kss = kkss.val_grid_point(1, k, j, 0) ;
214  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
215  min_kss = kkss.val_grid_point(1, k, j, 0) ;
216  }
217 
218  if (sortie != 0) {
219  fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
220  fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
221  // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
222  fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
223  fich_kss << step << " " << max_kss << " " << min_kss << endl ;
224  }
225 
226  cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
227  step ++ ;
228  }
229 
230  // LOOP WITH FIXED OMEGA :
231 
232  if (nb_it !=0)
233  cout << "OMEGA FIXED" << endl ;
234  double erreur ;
235 
236  for (int pas = 0 ; pas <nb_it ; pas ++) {
237 
238  Scalar beta_un_old (hole1.beta_auto(1)) ;
239 
240  solve_shift (precis, relax, bound_beta, omega_eff) ;
242 
243  solve_psi (precis, relax, bound_psi) ;
244  solve_lapse (precis, relax, bound_nn, lim_nn) ;
245 
246  // Convergence to the given irreductible mass
247  if (search_mass == 1 && step >= 30) {
248  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
249  sqrt(hole2.area_hor()/16/M_PI) ;
250  double error_m = (mass_irr - mass_area) / mass_irr ;
251  double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
252 
253  hole1.mp.homothetie_interne(scaling_r) ;
254  hole1.radius = hole1.radius *scaling_r ;
255  hole2.mp.homothetie_interne(scaling_r) ;
256  hole2.radius = hole2.radius *scaling_r ;
257  }
258 
259  erreur = 0 ;
260  Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
261  for (int i=1 ; i<nz ; i++)
262  if (diff(i) > erreur)
263  erreur = diff(i) ;
264 
265  // Saving ok K_{ij}s^is^j
266  // -----------------------
267 
268  Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
270  hole1.get_gam().radial_vect(), 0, 1)) ;
271  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
272  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
273  int nnp = hole1.mp.get_mg()->get_np(1) ;
274  int nnt = hole2.mp.get_mg()->get_nt(1) ;
275  for (int k=0 ; k<nnp ; k++)
276  for (int j=0 ; j<nnt ; j++){
277  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
278  max_kss = kkss.val_grid_point(1, k, j, 0) ;
279  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
280  min_kss = kkss.val_grid_point(1, k, j, 0) ;
281  }
282 
283 
284  if (sortie != 0) {
285  fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
286  fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
287  // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
288  fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
289  fich_kss << step << " " << max_kss << " " << min_kss << endl ;
290  }
291 
292  cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
293  step ++ ;
294  }
295 
296  if (nb_it != 0){
297  fich_iteration << "#----------------------------" << endl ;
298  fich_correction << "#-----------------------------" << endl ;
299  fich_viriel << "#------------------------------" << endl ;
300  }
301 
302  return viriel() ;
303 }
304 }
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:342
double omega_hor() const
Orbital velocity.
Definition: single_param.C:163
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:91
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
double area_hor() const
Area of the horizon.
Definition: single_param.C:91
const Sym_tensor & get_k_dd() const
k_dd
Definition: single_hor.C:351
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
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:365
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:99
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:74
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
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:136
void init_bin_hor()
Initialisation of the system.
Definition: bin_hor.C:164
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:465
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
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:325
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 homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
Definition: map_af.C:744
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:539
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:324
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542