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
31char bin_hor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $" ;
32
33/*
34 * $Id: bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $
35 * $Log: bin_hor.C,v $
36 * Revision 1.12 2014/10/13 08:52:42 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.11 2014/10/06 15:13:00 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.10 2007/04/13 15:28:55 f_limousin
43 * Lots of improvements, generalisation to an arbitrary state of
44 * rotation, implementation of the spatial metric given by Samaya.
45 *
46 * Revision 1.9 2006/08/01 14:37:19 f_limousin
47 * New version
48 *
49 * Revision 1.8 2006/06/29 08:51:00 f_limousin
50 * *** empty log message ***
51 *
52 * Revision 1.7 2006/06/28 13:36:09 f_limousin
53 * Convergence to a given irreductible mass
54 *
55 * Revision 1.6 2006/05/24 16:56:37 f_limousin
56 * Many small modifs.
57 *
58 * Revision 1.5 2005/06/13 15:47:29 jl_jaramillo
59 * Add some quatities in write_global()
60 *
61 * Revision 1.4 2005/06/09 16:12:04 f_limousin
62 * Implementation of amg_mom_adm().
63 *
64 * Revision 1.3 2005/04/29 14:02:44 f_limousin
65 * Important changes : manage the dependances between quantities (for
66 * instance psi and psi4). New function write_global(ost).
67 *
68 * Revision 1.2 2005/03/04 09:38:41 f_limousin
69 * Implement the constructor from a file, operator>>, operator<<
70 * and function sauve.
71 *
72 * Revision 1.1 2004/12/29 16:11:02 f_limousin
73 * First version
74 *
75 *
76 * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $
77 *
78 */
79
80//standard
81#include <cstdlib>
82#include <cmath>
83
84// Lorene
85#include "nbr_spx.h"
86#include "tenseur.h"
87#include "tensor.h"
88#include "isol_hor.h"
89#include "proto.h"
90#include "utilitaires.h"
91//#include "graphique.h"
92
93// Standard constructor
94// --------------------
95
96namespace Lorene {
98 hole1(mp1), hole2(mp2), omega(0){
99
100 holes[0] = &hole1 ;
101 holes[1] = &hole2 ;
102}
103
104// Copy constructor
105// ----------------
106
107Bin_hor::Bin_hor (const Bin_hor& source) :
108 hole1(source.hole1), hole2(source.hole2), omega(source.omega) {
109
110 holes[0] = &hole1 ;
111 holes[1] = &hole2 ;
112 }
113
114// Constructor from a file
115// -----------------------
116
117Bin_hor::Bin_hor(Map_af& mp1, Map_af& mp2, FILE* fich)
118 : hole1(mp1, fich),
119 hole2(mp2, fich),
120 omega(0) {
121
122 fread_be(&omega, sizeof(double), 1, fich) ;
123 holes[0] = &hole1 ;
124 holes[1] = &hole2 ;
125
126}
127
128 //--------------//
129 // Destructor //
130 //--------------//
131
134
135 //-----------------------//
136 // Mutators / assignment //
137 //-----------------------//
138
139void Bin_hor::operator= (const Bin_hor& source) {
140 hole1 = source.hole1 ;
141 hole2 = source.hole2 ;
142
143 omega = source.omega ;
144}
145
146
147 //--------------------------//
148 // Save in a file //
149 //--------------------------//
150
151void Bin_hor::sauve(FILE* fich) const {
152
153 hole1.sauve(fich) ;
154 hole2.sauve(fich) ;
155 fwrite_be (&omega, sizeof(double), 1, fich) ;
156
157}
158
159
160//Initialisation : Sum of two static BH
176
177
178void Bin_hor::write_global(ostream& ost, double lim_nn, int bound_nn,
179 int bound_psi, int bound_beta, double alpha) const {
180
181 double distance = hole1.get_mp().get_ori_x() - hole2.get_mp().get_ori_x() ;
182 double mass_adm = adm_mass() ;
183 double mass_komar = komar_mass() ;
184 double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
185 sqrt(hole2.area_hor()/16/M_PI) ;
186 double J_adm = ang_mom_adm() ;
187 double J_hor = ang_mom_hor() ; //hole1.ang_mom_hor() + hole2.ang_mom_hor() ;
188 double j1 = hole1.ang_mom_hor() ;
189 double j2 = hole2.ang_mom_hor() ;
190 double mass_ih1 = hole1.mass_hor() ;
191 double mass_ih2 = hole2.mass_hor() ;
192 double mass_ih = mass_ih1 + mass_ih2 ;
193 double omega1 = hole1.omega_hor() ;
194 double omega2 = hole2.omega_hor() ;
195
196 // Verification of Smarr :
197 // -----------------------
198
199 // Les alignemenents pour le signe des CL.
200 double orientation1 = hole1.mp.get_rot_phi() ;
201 assert ((orientation1 == 0) || (orientation1 == M_PI)) ;
202 int aligne1 = (orientation1 == 0) ? 1 : -1 ;
203
204 Vector angular1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
205 Scalar yya1 (hole1.mp) ;
206 yya1 = hole1.mp.ya ;
207 Scalar xxa1 (hole1.mp) ;
208 xxa1 = hole1.mp.xa ;
209
210 angular1.set(1) = aligne1 * omega * yya1 ;
211 angular1.set(2) = - aligne1 * omega * xxa1 ;
212 angular1.set(3).annule_hard() ;
213
214 angular1.set(1).set_spectral_va()
216 angular1.set(2).set_spectral_va()
218 angular1.set(3).set_spectral_va()
220
222
223
224 double orientation2 = hole2.mp.get_rot_phi() ;
225 assert ((orientation2 == 0) || (orientation2 == M_PI)) ;
226 int aligne2 = (orientation2 == 0) ? 1 : -1 ;
227
228 Vector angular2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
229 Scalar yya2 (hole2.mp) ;
230 yya2 = hole2.mp.ya ;
231 Scalar xxa2 (hole2.mp) ;
232 xxa2 = hole2.mp.xa ;
233
234 angular2.set(1) = aligne2 * omega * yya2 ;
235 angular2.set(2) = - aligne2 * omega * xxa2 ;
236 angular2.set(3).annule_hard() ;
237
238 angular2.set(1).set_spectral_va()
240 angular2.set(2).set_spectral_va()
242 angular2.set(3).set_spectral_va()
244
246
247
248 Scalar btilde1 (hole1.b_tilde() -
249 contract(angular1, 0, hole1.tgam.radial_vect().up_down(hole1.tgam), 0)) ;
250 Scalar btilde2 (hole2.b_tilde() -
251 contract(angular2, 0, hole2.tgam.radial_vect().up_down(hole2.tgam), 0)) ;
252
253
254
255
256 Vector integrand_un (hole1.mp, COV, hole1.mp.get_bvect_spher()) ;
257 integrand_un = hole1.nn.derive_cov(hole1.ff)*pow(hole1.psi, 2)
258 - btilde1*contract(hole1.get_k_dd(), 1,
259 hole1.tgam.radial_vect(), 0)*pow(hole1.psi, 2) ;
260 integrand_un.std_spectral_base() ;
261
262 Vector integrand_deux (hole2.mp, COV, hole2.mp.get_bvect_spher()) ;
263 integrand_deux = hole2.nn.derive_cov(hole2.ff)*pow(hole2.psi, 2)
264 - btilde2*contract(hole2.get_k_dd(), 1,
265 hole2.tgam.radial_vect(), 0)*pow(hole2.psi, 2) ;
266 integrand_deux.std_spectral_base() ;
267
268 double horizon = hole1.mp.integrale_surface(integrand_un(1),
269 hole1.get_radius())+
270 hole2.mp.integrale_surface(integrand_deux(1), hole2.get_radius()) ;
271
272 horizon /= 4*M_PI ;
273
274 double J_smarr = (mass_komar - horizon) / 2. / omega ;
275
276 ost.precision(8) ;
277 ost << "# Grid : " << hole1.mp.get_mg()->get_nr(1) << "x"
278 << hole1.mp.get_mg()->get_nt(1) << "x"
279 << hole1.mp.get_mg()->get_np(1) << " R_out(l) : " ;
280
281 for (int ll=0; ll<hole1.mp.get_mg()->get_nzone(); ll++) {
282 ost << " " << hole1.mp.val_r(ll, 1., M_PI/2, 0) ;
283 }
284 ost << endl ;
285 ost << "# bound N, lim N : " << bound_nn << " " << lim_nn
286 << " - bound Psi : " << bound_psi << " - bound shift : " << bound_beta
287 << " alpha = " << alpha << endl ;
288
289 ost << "# distance omega Mass_ADM Mass_K M_area J_ADM J_hor" << endl ;
290 ost << distance << " " ;
291 ost << omega << " " ;
292 ost << mass_adm << " " ;
293 ost << mass_komar << " " ;
294 ost << mass_area << " " ;
295 ost << J_adm << " " ;
296 ost << J_hor << endl ;
297 ost << "# mass_ih1 mass_ih2 mass_ih j1 J2 omega1 omega2" << endl ;
298 ost << mass_ih1 << " " ;
299 ost << mass_ih2 << " " ;
300 ost << mass_ih << " " ;
301 ost << j1 << " " ;
302 ost << j2 << " " ;
303 ost << omega1 << " " ;
304 ost << omega2 << endl ;
305 ost << "# ADM_mass/M_area J/M_area2 omega*M_area" << endl ;
306 ost << mass_adm / mass_area << " " ;
307 ost << J_adm /mass_area / mass_area << " " ;
308 ost << omega * mass_area << endl ;
309 ost << "# Diff J_hor/J_ADM Diff J_ADM/J_Smarr Diff J_hor/J_smarr"
310 << endl ;
311 ost << fabs(J_adm - J_hor) / J_adm << " " << fabs(J_adm - J_smarr) / J_adm
312 << " " << fabs(J_hor - J_smarr) / J_hor << endl ;
313
314
315}
316
317}
double omega
Angular velocity.
Definition isol_hor.h:1348
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Single_hor * holes[2]
Array on the black holes.
Definition isol_hor.h:1346
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
void init_bin_hor()
Initialisation of the system.
Definition bin_hor.C:161
double adm_mass() const
Calculates the ADM mass of the system.
double ang_mom_adm() const
Calculates the angular momentum of the black hole.
void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition bin_hor.C:151
double komar_mass() const
Calculates the Komar mass of the system using : .
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition binhor_kij.C:473
void operator=(const Bin_hor &)
Affectation operator.
Definition bin_hor.C:139
Bin_hor(Map_af &mp1, Map_af &mp2)
Standard constructor.
Definition bin_hor.C:97
virtual ~Bin_hor()
Destructor.
Definition bin_hor.C:132
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition isol_hor.h:1409
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:178
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition binhor_kij.C:88
Affine radial mapping.
Definition map.h:2027
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
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.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord ya
Absolute y coordinate.
Definition map.h:731
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
double get_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition metric.C:362
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
double mass_hor() const
Mass computed at the horizon
const Sym_tensor & get_k_dd() const
k_dd
Definition single_hor.C:348
Metric_flat ff
3 metric flat
Definition isol_hor.h:980
Scalar psi
Conformal factor .
Definition isol_hor.h:930
Map_af & mp
Affine mapping.
Definition isol_hor.h:900
double area_hor() const
Area of the horizon.
void init_bhole()
Sets the values of the fields to :
Definition single_hor.C:484
const Map_af & get_mp() const
Returns the mapping (readonly).
Definition isol_hor.h:1038
double ang_mom_hor() const
Angular momentum (modulo)
double omega_hor() const
Orbital velocity
virtual void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition single_hor.C:244
double get_radius() const
Returns the radius of the horizon.
Definition isol_hor.h:1046
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Definition single_hor.C:390
Scalar nn
Lapse function .
Definition isol_hor.h:921
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Definition single_hor.C:424
Metric tgam
3 metric tilde
Definition isol_hor.h:977
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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:69
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:70
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64