LORENE
bin_ns_bh_coal.C
1 /*
2  * Copyright (c) 2004 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 bin_ns_bh_coal_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $" ;
24 
25 /*
26  * $Id: bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $
27  * $Log: bin_ns_bh_coal.C,v $
28  * Revision 1.12 2014/10/13 08:52:43 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.11 2014/10/06 15:13:01 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.10 2007/04/24 20:13:53 f_limousin
35  * Implementation of Dirichlet and Neumann BC for the lapse
36  *
37  * Revision 1.9 2006/09/05 13:39:43 p_grandclement
38  * update of the bin_ns_bh project
39  *
40  * Revision 1.8 2006/06/23 07:09:24 p_grandclement
41  * Addition of spinning black hole
42  *
43  * Revision 1.7 2006/06/01 12:47:52 p_grandclement
44  * update of the Bin_ns_bh project
45  *
46  * Revision 1.6 2006/04/27 09:12:33 p_grandclement
47  * First try at irrotational black holes
48  *
49  * Revision 1.5 2006/04/25 07:21:57 p_grandclement
50  * Various changes for the NS_BH project
51  *
52  * Revision 1.4 2006/03/30 07:33:45 p_grandclement
53  * *** empty log message ***
54  *
55  * Revision 1.3 2005/12/01 12:59:10 p_grandclement
56  * Files for bin_ns_bh project
57  *
58  * Revision 1.2 2005/10/18 13:12:32 p_grandclement
59  * update of the mixted binary codes
60  *
61  * Revision 1.1 2005/08/29 15:10:15 p_grandclement
62  * Addition of things needed :
63  * 1) For BBH with different masses
64  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
65  * WORKING YET !!!
66  *
67  *
68  * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $
69  *
70  */
71 
72 //standard
73 #include <cstdlib>
74 
75 // Lorene
76 #include "tenseur.h"
77 #include "bin_ns_bh.h"
78 #include "unites.h"
79 #include "graphique.h"
80 
81 namespace Lorene {
82 void Bin_ns_bh::coal (double precis, double relax, int itemax_equil,
83  int itemax_mp_et, double ent_c_init, double seuil_masses,
84  double dist, double m1, double m2, double spin_cible,
85  double scale_ome_local, const int sortie, int bound_nn,
86  double lim_nn) {
87 
88  using namespace Unites ;
89 
90  int nz_bh = hole.mp.get_mg()->get_nzone() ;
91  int nz_ns = star.mp.get_mg()->get_nzone() ;
92 
93  // Distance initiale
94  double distance = fabs(hole.mp.get_ori_x()-star.mp.get_ori_x()) ;
95  double mass_ns = star.mass_g() * ggrav;
96  double mass_bh = hole.masse_adm_seul() ;
97  double axe_pos = star.mp.get_ori_x() ;
98  double scale_linear = (mass_ns+mass_bh)/2.*distance*omega ;
99 
100  char name_iteration[40] ;
101  char name_correction[40] ;
102  char name_viriel[40] ;
103  char name_ome [40] ;
104  char name_linear[40] ;
105  char name_axe[40] ;
106  char name_error_m1[40] ;
107  char name_error_m2[40] ;
108  char name_hor[40] ;
109  char name_entc[40] ;
110  char name_dist[40] ;
111  char name_spin[40] ;
112  char name_ome_local[40] ;
113 
114  sprintf(name_iteration, "ite.dat") ;
115  sprintf(name_correction, "cor.dat") ;
116  sprintf(name_viriel, "vir.dat") ;
117  sprintf(name_ome, "ome.dat") ;
118  sprintf(name_linear, "linear.dat") ;
119  sprintf(name_axe, "axe.dat") ;
120  sprintf(name_error_m1, "error_m_bh.dat") ;
121  sprintf(name_error_m2, "error_m_ns.dat") ;
122  sprintf(name_hor, "hor.dat") ;
123  sprintf(name_entc, "entc.dat") ;
124  sprintf(name_dist, "distance.dat") ;
125  sprintf(name_spin, "spin.dat") ;
126  sprintf(name_ome_local, "ome_local.dat") ;
127 
128  ofstream fiche_iteration(name_iteration) ;
129  fiche_iteration.precision(8) ;
130 
131  ofstream fiche_correction(name_correction) ;
132  fiche_correction.precision(8) ;
133 
134  ofstream fiche_viriel(name_viriel) ;
135  fiche_viriel.precision(8) ;
136 
137  ofstream fiche_ome(name_ome) ;
138  fiche_ome.precision(8) ;
139 
140  ofstream fiche_linear(name_linear) ;
141  fiche_linear.precision(8) ;
142 
143  ofstream fiche_axe(name_axe) ;
144  fiche_axe.precision(8) ;
145 
146  ofstream fiche_error_m1 (name_error_m1) ;
147  fiche_error_m1.precision(8) ;
148 
149  ofstream fiche_error_m2 (name_error_m2) ;
150  fiche_error_m2.precision(8) ;
151 
152  ofstream fiche_hor (name_hor) ;
153  fiche_hor.precision(8) ;
154 
155  ofstream fiche_entc (name_entc) ;
156  fiche_entc.precision(8) ;
157 
158  ofstream fiche_dist (name_dist) ;
159  fiche_dist.precision(8) ;
160 
161  ofstream fiche_spin (name_spin) ;
162  fiche_spin.precision(8) ;
163 
164  ofstream fiche_ome_local (name_ome_local) ;
165  fiche_ome_local.precision(8) ;
166 
167  bool loop = true ;
168  bool search_masses = false ;
169  double ent_c = ent_c_init ;
170 
171  Cmp shift_bh_old (hole.mp) ;
172  Cmp shift_ns_old (star.mp) ;
173 
174  double erreur = 1 ;
175 
176  int conte = 0 ;
177 
178  double old_mass_ns ;
179  mass_ns = star.mass_b() ;
180  double spin_old ;
181  double spin = 1;
182  bool adapt = true ;
183 
184  while (loop) {
185 
186  spin_old = spin ;
187  spin = hole.local_momentum() ;
188  if (sortie !=0) {
189  fiche_ome_local << conte << " " << hole.omega_local << endl ;
190  fiche_spin << conte << " " << spin/m1/m1 << endl ;
191  }
192 
193  double conv_spin = fabs(spin-spin_old) ;
194  double error_spin = spin - spin_cible ;
195  double rel_diff_spin = (spin_cible==0) ? fabs(error_spin) :
196  fabs(error_spin)/spin_cible ;
197  if ((conv_spin*2<rel_diff_spin) && (search_masses) &&
198  hole.get_rot_state() != COROT) {
199  double func = scale_ome_local*log((2+error_spin)/(2+2*error_spin));
201  }
202 
203  old_mass_ns = mass_ns ;
204 
205  if (hole.get_shift_auto().get_etat() != ETATZERO)
206  shift_bh_old = hole.get_shift_auto()(0) ;
207  else
208  shift_bh_old = 0 ;
209 
210  if (star.get_shift_auto().get_etat() != ETATZERO)
211  shift_ns_old = star.get_shift_auto()(0) ;
212  else
213  shift_ns_old = 0 ;
214 
216  star.fait_d_psi() ;
217  star.hydro_euler() ;
218 
219  Tbl diff (7) ;
220  diff.set_etat_qcq() ;
221  int ite ;
222 
223 
224  star.equilibrium_nsbh (adapt, ent_c, ite, itemax_equil, itemax_mp_et,
225  relax, itemax_mp_et, relax, diff) ;
226 
227  hole.update_metric(star) ;
228 
229  hole.equilibrium (star, precis, relax, bound_nn, lim_nn) ;
230  cout << "Apres equilibrium" << endl ;
231 
233  cout << "Apres star::update_metric" << endl ;
234 
236  cout << "Apres star::update_metric_der_comp" << endl ;
237  fait_tkij(bound_nn, lim_nn) ;
238  cout << "Apres Bin_ns_bh::fait_tkij" << endl ;
239 
240  erreur = 0 ;
241  Tbl diff_bh (diffrelmax (shift_bh_old, hole.get_shift_auto()(0))) ;
242  for (int i=1 ; i<nz_bh ; i++)
243  if (diff_bh(i) > erreur)
244  erreur = diff_bh(i) ;
245  Tbl diff_ns (diffrelmax (shift_ns_old, star.get_shift_auto()(0))) ;
246  for (int i=0 ; i<nz_ns ; i++)
247  if (diff_ns(i) > erreur)
248  erreur = diff_ns(i) ;
249 
250  if (erreur<seuil_masses)
251  search_masses = true ;
252 
253  mass_ns = star.mass_b() ;
254 
255  cout << "Avant viriel" << endl ;
256  double error_viriel = viriel() ;
257  cout << "Apres viriel" << endl ;
258  double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
259  cout << "Apres linear" << endl ;
260  double error_m1 = 1.-sqrt(hole.area()/16./M_PI)/m1 ;
261  cout << "Apres Mbh" << endl ;
262  double error_m2 = 1 - mass_ns/m2 ;
263  cout << "Apres Mns" << endl ;
264 
265  if (sortie != 0) {
266  fiche_iteration << conte << " " << erreur << endl ;
267  fiche_correction << conte << " " << hole.regul << endl ;
268  fiche_viriel << conte << " " << error_viriel << endl ;
269  fiche_linear << conte << " " << error_linear << endl ;
270  fiche_error_m1 << conte << " " << error_m1 << endl ;
271  fiche_error_m2 << conte << " " << error_m2 << endl ;
272  fiche_hor << conte << " " << hole.mp.val_r(0, 1, 0,0) << endl ;
273  fiche_entc << conte << " " << ent_c << endl ;
274  fiche_dist << conte << " " << distance << endl ;
275  fiche_ome << conte << " " << omega << endl ;
276  fiche_axe << conte << " " << axe_pos << endl ;
277  }
278 
279  // The axis position
280  double scaling_axe = pow((2+error_linear)/(2+2*error_linear), 0.1) ;
281  axe_pos *= scaling_axe ;
282  star.set_mp().set_ori (axe_pos, 0, 0) ;
283  hole.set_mp().set_ori (-distance + axe_pos, 0, 0) ;
284 
285  // Value of omega
286  double new_ome = star.compute_angul() ;
287  if (new_ome !=0) {
288  set_omega(new_ome) ;
289  if (hole.get_rot_state() == COROT)
290  hole.set_omega_local(new_ome) ;
291  }
292 
293  // The right distance
294  double error_dist = (distance-dist)/dist ;
295  double scale_d = pow((2+error_dist)/(2+2*error_dist), 0.2) ;
296  distance *= scale_d ;
297 
298 
299  // Converge to the right masses :
300  if (search_masses) {
301 
302  double scaling_r = pow((2-error_m1)/(2-2*error_m1), 0.25) ;
303  hole.mp.homothetie_interne(scaling_r) ;
304  hole.set_rayon(hole.get_rayon()*scaling_r) ;
305  }
306 
307  // The case of the NS :
308  double convergence = fabs(mass_ns - old_mass_ns)/mass_ns ;
309  double rel_diff = fabs(error_m2) ;
310  if ((search_masses) && (convergence*2 < rel_diff)) {
311  double scaling_ent = pow((2-error_m2)/(2-2*error_m2), 1) ;
312  ent_c *= scaling_ent ;
313 
314  }
315 
316 
317 
318  cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
319  //if (erreur < 1e-4)
320  //adapt = false ;
321 
322  if (erreur < precis)
323  loop = false ;
324  conte ++ ;
325  }
326 
327 
328  fiche_iteration.close() ;
329  fiche_correction.close() ;
330  fiche_viriel.close() ;
331  fiche_ome.close() ;
332  fiche_linear.close() ;
333  fiche_axe.close() ;
334  fiche_error_m1.close() ;
335  fiche_error_m2.close() ;
336  fiche_hor.close() ;
337  fiche_entc.close() ;
338  fiche_dist.close() ;
339  fiche_spin.close() ;
340  fiche_ome_local.close() ;
341 }
342 }
Lorene::Map_af::homothetie_interne
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
Lorene::Bhole::mp
Map_af & mp
Affine mapping.
Definition: bhole.h:273
Lorene::Bhole::get_rayon
double get_rayon() const
Returns the radius of the horizon.
Definition: bhole.h:347
Lorene::Map::get_ori_x
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Lorene::Bhole::omega_local
double omega_local
local angular velocity
Definition: bhole.h:276
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Lorene::Etoile_bin::fait_d_psi
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition: etoile_bin.C:764
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Bin_ns_bh::x_axe
double x_axe
Absolute X coordinate of the rotation axis.
Definition: bin_ns_bh.h:140
Lorene::Bhole::set_mp
Map_af & set_mp()
Read/write of the mapping.
Definition: bhole.h:342
Lorene::Map::set_ori
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Definition: map.C:253
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:539
Lorene::Bhole::get_rot_state
double get_rot_state() const
Returns the state of rotation.
Definition: bhole.h:373
Lorene::Bin_ns_bh::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_ns_bh.h:136
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Map_af::val_r
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:96
Lorene::Bhole::area
double area() const
Computes the area of the throat.
Definition: bhole.C:545
Lorene::Etoile_bin::get_shift_auto
const Tenseur & get_shift_auto() const
Returns the part of the shift vector generated principaly by the star.
Definition: etoile.h:1152
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
Unites
Standard units of space, time and mass.
Lorene::Et_bin_nsbh::equilibrium_nsbh
virtual void equilibrium_nsbh(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration in a NS-BH binary system.
Definition: et_bin_nsbh_equilibrium.C:494
Lorene::Bhole::get_shift_auto
const Tenseur & get_shift_auto() const
Returns the part of generated by the hole.
Definition: bhole.h:440
Lorene::Etoile_bin::hydro_euler
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition: et_bin_hydro.C:106
Lorene::Bin_ns_bh::hole
Bhole hole
The black hole.
Definition: bin_ns_bh.h:131
Lorene::Etoile_bin::mass_g
virtual double mass_g() const
Gravitational mass.
Definition: et_bin_global.C:137
Lorene::Bhole::set_rayon
void set_rayon(double ray)
Sets the radius of the horizon to ray .
Definition: bhole.h:352
Lorene::Bin_ns_bh::fait_tkij
void fait_tkij(int bound_nn=-1, double lim_nn=0)
Computation of the extrinsic curvature tensor for both {\tt star} and {\tt bhole}.
Definition: bin_ns_bh_kij.C:257
Lorene::Et_bin_nsbh::kinematics
virtual void kinematics(double omega, double x_axe)
Computes the quantities bsn and pot_centri .
Definition: et_bin_nsbh_kinema.C:51
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Etoile_bin::mass_b
virtual double mass_b() const
Baryon mass.
Definition: et_bin_global.C:100
Lorene::Bhole::set_omega_local
void set_omega_local(double ome)
Sets the local angular velocity to ome .
Definition: bhole.h:369
Lorene::Et_bin_nsbh::update_metric
void update_metric(const Bhole &comp)
Computes metric coefficients from known potentials, when the companion is a black hole.
Definition: et_bin_nsbh_upmetr.C:64
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Bin_ns_bh::star
Et_bin_nsbh star
The neutron star.
Definition: bin_ns_bh.h:128
Lorene::Bhole::regul
double regul
Intensity of the correction on the shift vector.
Definition: bhole.h:284
Lorene::Etoile::set_mp
Map & set_mp()
Read/write of the mapping.
Definition: etoile.h:601
Lorene::Bhole::masse_adm_seul
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
Definition: bhole_pseudo_kerr.C:344
Lorene::Bin_ns_bh::set_omega
void set_omega(double)
Sets the orbital angular velocity [{\tt f_unit}].
Definition: bin_ns_bh.C:218
Lorene::Et_bin_nsbh::update_metric_der_comp
void update_metric_der_comp(const Bhole &comp)
Computes the derivative of metric functions related to the companion black hole.
Definition: et_bin_nsbh_upmetr_der.C:71