LORENE
bin_hor.C
1 /*
2  * Methods of class Bin_hor
3  *
4  * (see file bin_hor.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004-2005 Francois Limousin
10  * Jose Luis Jaramillo
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 as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: bin_hor.C,v 1.13 2016/12/05 16:17:46 j_novak Exp $
35  * $Log: bin_hor.C,v $
36  * Revision 1.13 2016/12/05 16:17:46 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.12 2014/10/13 08:52:42 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.11 2014/10/06 15:13:00 j_novak
43  * Modified #include directives to use c++ syntax.
44  *
45  * Revision 1.10 2007/04/13 15:28:55 f_limousin
46  * Lots of improvements, generalisation to an arbitrary state of
47  * rotation, implementation of the spatial metric given by Samaya.
48  *
49  * Revision 1.9 2006/08/01 14:37:19 f_limousin
50  * New version
51  *
52  * Revision 1.8 2006/06/29 08:51:00 f_limousin
53  * *** empty log message ***
54  *
55  * Revision 1.7 2006/06/28 13:36:09 f_limousin
56  * Convergence to a given irreductible mass
57  *
58  * Revision 1.6 2006/05/24 16:56:37 f_limousin
59  * Many small modifs.
60  *
61  * Revision 1.5 2005/06/13 15:47:29 jl_jaramillo
62  * Add some quatities in write_global()
63  *
64  * Revision 1.4 2005/06/09 16:12:04 f_limousin
65  * Implementation of amg_mom_adm().
66  *
67  * Revision 1.3 2005/04/29 14:02:44 f_limousin
68  * Important changes : manage the dependances between quantities (for
69  * instance psi and psi4). New function write_global(ost).
70  *
71  * Revision 1.2 2005/03/04 09:38:41 f_limousin
72  * Implement the constructor from a file, operator>>, operator<<
73  * and function sauve.
74  *
75  * Revision 1.1 2004/12/29 16:11:02 f_limousin
76  * First version
77  *
78  *
79  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/bin_hor.C,v 1.13 2016/12/05 16:17:46 j_novak Exp $
80  *
81  */
82 
83 //standard
84 #include <cstdlib>
85 #include <cmath>
86 
87 // Lorene
88 #include "nbr_spx.h"
89 #include "tenseur.h"
90 #include "tensor.h"
91 #include "isol_hor.h"
92 #include "proto.h"
93 #include "utilitaires.h"
94 //#include "graphique.h"
95 
96 // Standard constructor
97 // --------------------
98 
99 namespace Lorene {
101  hole1(mp1), hole2(mp2), omega(0){
102 
103  holes[0] = &hole1 ;
104  holes[1] = &hole2 ;
105 }
106 
107 // Copy constructor
108 // ----------------
109 
110 Bin_hor::Bin_hor (const Bin_hor& source) :
111  hole1(source.hole1), hole2(source.hole2), omega(source.omega) {
112 
113  holes[0] = &hole1 ;
114  holes[1] = &hole2 ;
115  }
116 
117 // Constructor from a file
118 // -----------------------
119 
120 Bin_hor::Bin_hor(Map_af& mp1, Map_af& mp2, FILE* fich)
121  : hole1(mp1, fich),
122  hole2(mp2, fich),
123  omega(0) {
124 
125  fread_be(&omega, sizeof(double), 1, fich) ;
126  holes[0] = &hole1 ;
127  holes[1] = &hole2 ;
128 
129 }
130 
131  //--------------//
132  // Destructor //
133  //--------------//
134 
136 }
137 
138  //-----------------------//
139  // Mutators / assignment //
140  //-----------------------//
141 
142 void Bin_hor::operator= (const Bin_hor& source) {
143  hole1 = source.hole1 ;
144  hole2 = source.hole2 ;
145 
146  omega = source.omega ;
147 }
148 
149 
150  //--------------------------//
151  // Save in a file //
152  //--------------------------//
153 
154 void Bin_hor::sauve(FILE* fich) const {
155 
156  hole1.sauve(fich) ;
157  hole2.sauve(fich) ;
158  fwrite_be (&omega, sizeof(double), 1, fich) ;
159 
160 }
161 
162 
163 //Initialisation : Sum of two static BH
165  set_omega (0) ;
166  hole1.init_bhole() ;
167  hole2.init_bhole() ;
168 
171 
174 
175  decouple() ;
177 
178 }
179 
180 
181 void Bin_hor::write_global(ostream& ost, double lim_nn, int bound_nn,
182  int bound_psi, int bound_beta, double alpha) const {
183 
184  double distance = hole1.get_mp().get_ori_x() - hole2.get_mp().get_ori_x() ;
185  double mass_adm = adm_mass() ;
186  double mass_komar = komar_mass() ;
187  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
188  sqrt(hole2.area_hor()/16/M_PI) ;
189  double J_adm = ang_mom_adm() ;
190  double J_hor = ang_mom_hor() ; //hole1.ang_mom_hor() + hole2.ang_mom_hor() ;
191  double j1 = hole1.ang_mom_hor() ;
192  double j2 = hole2.ang_mom_hor() ;
193  double mass_ih1 = hole1.mass_hor() ;
194  double mass_ih2 = hole2.mass_hor() ;
195  double mass_ih = mass_ih1 + mass_ih2 ;
196  double omega1 = hole1.omega_hor() ;
197  double omega2 = hole2.omega_hor() ;
198 
199  // Verification of Smarr :
200  // -----------------------
201 
202  // Les alignemenents pour le signe des CL.
203  double orientation1 = hole1.mp.get_rot_phi() ;
204  assert ((orientation1 == 0) || (orientation1 == M_PI)) ;
205  int aligne1 = (orientation1 == 0) ? 1 : -1 ;
206 
207  Vector angular1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
208  Scalar yya1 (hole1.mp) ;
209  yya1 = hole1.mp.ya ;
210  Scalar xxa1 (hole1.mp) ;
211  xxa1 = hole1.mp.xa ;
212 
213  angular1.set(1) = aligne1 * omega * yya1 ;
214  angular1.set(2) = - aligne1 * omega * xxa1 ;
215  angular1.set(3).annule_hard() ;
216 
217  angular1.set(1).set_spectral_va()
218  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[0])) ;
219  angular1.set(2).set_spectral_va()
220  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[1])) ;
221  angular1.set(3).set_spectral_va()
222  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[2])) ;
223 
224  angular1.change_triad(hole1.mp.get_bvect_spher()) ;
225 
226 
227  double orientation2 = hole2.mp.get_rot_phi() ;
228  assert ((orientation2 == 0) || (orientation2 == M_PI)) ;
229  int aligne2 = (orientation2 == 0) ? 1 : -1 ;
230 
231  Vector angular2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
232  Scalar yya2 (hole2.mp) ;
233  yya2 = hole2.mp.ya ;
234  Scalar xxa2 (hole2.mp) ;
235  xxa2 = hole2.mp.xa ;
236 
237  angular2.set(1) = aligne2 * omega * yya2 ;
238  angular2.set(2) = - aligne2 * omega * xxa2 ;
239  angular2.set(3).annule_hard() ;
240 
241  angular2.set(1).set_spectral_va()
242  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[0])) ;
243  angular2.set(2).set_spectral_va()
244  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[1])) ;
245  angular2.set(3).set_spectral_va()
246  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[2])) ;
247 
248  angular2.change_triad(hole2.mp.get_bvect_spher()) ;
249 
250 
251  Scalar btilde1 (hole1.b_tilde() -
252  contract(angular1, 0, hole1.tgam.radial_vect().up_down(hole1.tgam), 0)) ;
253  Scalar btilde2 (hole2.b_tilde() -
254  contract(angular2, 0, hole2.tgam.radial_vect().up_down(hole2.tgam), 0)) ;
255 
256 
257 
258 
259  Vector integrand_un (hole1.mp, COV, hole1.mp.get_bvect_spher()) ;
260  integrand_un = hole1.nn.derive_cov(hole1.ff)*pow(hole1.psi, 2)
261  - btilde1*contract(hole1.get_k_dd(), 1,
262  hole1.tgam.radial_vect(), 0)*pow(hole1.psi, 2) ;
263  integrand_un.std_spectral_base() ;
264 
265  Vector integrand_deux (hole2.mp, COV, hole2.mp.get_bvect_spher()) ;
266  integrand_deux = hole2.nn.derive_cov(hole2.ff)*pow(hole2.psi, 2)
267  - btilde2*contract(hole2.get_k_dd(), 1,
268  hole2.tgam.radial_vect(), 0)*pow(hole2.psi, 2) ;
269  integrand_deux.std_spectral_base() ;
270 
271  double horizon = hole1.mp.integrale_surface(integrand_un(1),
272  hole1.get_radius())+
273  hole2.mp.integrale_surface(integrand_deux(1), hole2.get_radius()) ;
274 
275  horizon /= 4*M_PI ;
276 
277  double J_smarr = (mass_komar - horizon) / 2. / omega ;
278 
279  ost.precision(8) ;
280  ost << "# Grid : " << hole1.mp.get_mg()->get_nr(1) << "x"
281  << hole1.mp.get_mg()->get_nt(1) << "x"
282  << hole1.mp.get_mg()->get_np(1) << " R_out(l) : " ;
283 
284  for (int ll=0; ll<hole1.mp.get_mg()->get_nzone(); ll++) {
285  ost << " " << hole1.mp.val_r(ll, 1., M_PI/2, 0) ;
286  }
287  ost << endl ;
288  ost << "# bound N, lim N : " << bound_nn << " " << lim_nn
289  << " - bound Psi : " << bound_psi << " - bound shift : " << bound_beta
290  << " alpha = " << alpha << endl ;
291 
292  ost << "# distance omega Mass_ADM Mass_K M_area J_ADM J_hor" << endl ;
293  ost << distance << " " ;
294  ost << omega << " " ;
295  ost << mass_adm << " " ;
296  ost << mass_komar << " " ;
297  ost << mass_area << " " ;
298  ost << J_adm << " " ;
299  ost << J_hor << endl ;
300  ost << "# mass_ih1 mass_ih2 mass_ih j1 J2 omega1 omega2" << endl ;
301  ost << mass_ih1 << " " ;
302  ost << mass_ih2 << " " ;
303  ost << mass_ih << " " ;
304  ost << j1 << " " ;
305  ost << j2 << " " ;
306  ost << omega1 << " " ;
307  ost << omega2 << endl ;
308  ost << "# ADM_mass/M_area J/M_area2 omega*M_area" << endl ;
309  ost << mass_adm / mass_area << " " ;
310  ost << J_adm /mass_area / mass_area << " " ;
311  ost << omega * mass_area << endl ;
312  ost << "# Diff J_hor/J_ADM Diff J_ADM/J_Smarr Diff J_hor/J_smarr"
313  << endl ;
314  ost << fabs(J_adm - J_hor) / J_adm << " " << fabs(J_adm - J_smarr) / J_adm
315  << " " << fabs(J_hor - J_smarr) / J_hor << endl ;
316 
317 
318 }
319 
320 }
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition: binhor_kij.C:476
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 mass_hor() const
Mass computed at the horizon.
Definition: single_param.C:138
double omega_hor() const
Orbital velocity.
Definition: single_param.C:163
double omega
Angular velocity.
Definition: isol_hor.h:1348
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
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
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
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
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:365
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
Tensor field of valence 1.
Definition: vector.h:188
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
const Map_af & get_mp() const
Returns the mapping (readonly).
Definition: isol_hor.h:1038
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition: isol_hor.h:1409
Single_hor * holes[2]
Array on the black holes.
Definition: isol_hor.h:1346
void operator=(const Bin_hor &)
Affectation operator.
Definition: bin_hor.C:142
Scalar psi
Conformal factor .
Definition: isol_hor.h:930
void init_bin_hor()
Initialisation of the system.
Definition: bin_hor.C:164
double get_radius() const
Returns the radius of the horizon.
Definition: isol_hor.h:1046
void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition: bin_hor.C:154
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void init_bhole()
Sets the values of the fields to :
Definition: single_hor.C:487
void write_global(ostream &, double lim_nn, int bound_nn, int bound_psi, int bound_beta, double alpha) const
Write global quantities in a formatted file.
Definition: bin_hor.C:181
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:73
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition: single_hor.C:247
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
double komar_mass() const
Calculates the Komar mass of the system using : .
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:72
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
Affine radial mapping.
Definition: map.h:2048
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
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Definition: single_hor.C:393
Bin_hor(Map_af &mp1, Map_af &mp2)
Standard constructor.
Definition: bin_hor.C:100
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Definition: single_hor.C:427
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
double ang_mom_adm() const
Calculates the angular momentum of the black hole.
double ang_mom_hor() const
Angular momentum (modulo)
Definition: single_param.C:110
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
double adm_mass() const
Calculates the ADM mass of the system.
Definition: binhor_global.C:86
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: single_param.C:71
Scalar nn
Lapse function .
Definition: isol_hor.h:921
virtual ~Bin_hor()
Destructor.
Definition: bin_hor.C:135