bin_ns_bh_coal.C
/*
* Copyright (c) 2004 Philippe Grandclement
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $
* $Log: bin_ns_bh_coal.C,v $
* Revision 1.12 2014/10/13 08:52:43 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.11 2014/10/06 15:13:01 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.10 2007/04/24 20:13:53 f_limousin
* Implementation of Dirichlet and Neumann BC for the lapse
*
* Revision 1.9 2006/09/05 13:39:43 p_grandclement
* update of the bin_ns_bh project
*
* Revision 1.8 2006/06/23 07:09:24 p_grandclement
* Addition of spinning black hole
*
* Revision 1.7 2006/06/01 12:47:52 p_grandclement
* update of the Bin_ns_bh project
*
* Revision 1.6 2006/04/27 09:12:33 p_grandclement
* First try at irrotational black holes
*
* Revision 1.5 2006/04/25 07:21:57 p_grandclement
* Various changes for the NS_BH project
*
* Revision 1.4 2006/03/30 07:33:45 p_grandclement
* *** empty log message ***
*
* Revision 1.3 2005/12/01 12:59:10 p_grandclement
* Files for bin_ns_bh project
*
* Revision 1.2 2005/10/18 13:12:32 p_grandclement
* update of the mixted binary codes
*
* Revision 1.1 2005/08/29 15:10:15 p_grandclement
* Addition of things needed :
* 1) For BBH with different masses
* 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
* WORKING YET !!!
*
*
* $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 $
*
*/
//standard
#include <cstdlib>
// Lorene
#include "tenseur.h"
#include "bin_ns_bh.h"
#include "unites.h"
#include "graphique.h"
namespace Lorene {
void Bin_ns_bh::coal (double precis, double relax, int itemax_equil,
int itemax_mp_et, double ent_c_init, double seuil_masses,
double dist, double m1, double m2, double spin_cible,
double scale_ome_local, const int sortie, int bound_nn,
double lim_nn) {
using namespace Unites ;
int nz_bh = hole.mp.get_mg()->get_nzone() ;
int nz_ns = star.mp.get_mg()->get_nzone() ;
// Distance initiale
double distance = fabs(hole.mp.get_ori_x()-star.mp.get_ori_x()) ;
double mass_ns = star.mass_g() * ggrav;
double mass_bh = hole.masse_adm_seul() ;
double axe_pos = star.mp.get_ori_x() ;
double scale_linear = (mass_ns+mass_bh)/2.*distance*omega ;
char name_iteration[40] ;
char name_correction[40] ;
char name_viriel[40] ;
char name_ome [40] ;
char name_linear[40] ;
char name_axe[40] ;
char name_error_m1[40] ;
char name_error_m2[40] ;
char name_hor[40] ;
char name_entc[40] ;
char name_dist[40] ;
char name_spin[40] ;
char name_ome_local[40] ;
sprintf(name_iteration, "ite.dat") ;
sprintf(name_correction, "cor.dat") ;
sprintf(name_viriel, "vir.dat") ;
sprintf(name_ome, "ome.dat") ;
sprintf(name_linear, "linear.dat") ;
sprintf(name_axe, "axe.dat") ;
sprintf(name_error_m1, "error_m_bh.dat") ;
sprintf(name_error_m2, "error_m_ns.dat") ;
sprintf(name_hor, "hor.dat") ;
sprintf(name_entc, "entc.dat") ;
sprintf(name_dist, "distance.dat") ;
sprintf(name_spin, "spin.dat") ;
sprintf(name_ome_local, "ome_local.dat") ;
ofstream fiche_iteration(name_iteration) ;
fiche_iteration.precision(8) ;
ofstream fiche_correction(name_correction) ;
fiche_correction.precision(8) ;
ofstream fiche_viriel(name_viriel) ;
fiche_viriel.precision(8) ;
ofstream fiche_ome(name_ome) ;
fiche_ome.precision(8) ;
ofstream fiche_linear(name_linear) ;
fiche_linear.precision(8) ;
ofstream fiche_axe(name_axe) ;
fiche_axe.precision(8) ;
ofstream fiche_error_m1 (name_error_m1) ;
fiche_error_m1.precision(8) ;
ofstream fiche_error_m2 (name_error_m2) ;
fiche_error_m2.precision(8) ;
ofstream fiche_hor (name_hor) ;
fiche_hor.precision(8) ;
ofstream fiche_entc (name_entc) ;
fiche_entc.precision(8) ;
ofstream fiche_dist (name_dist) ;
fiche_dist.precision(8) ;
ofstream fiche_spin (name_spin) ;
fiche_spin.precision(8) ;
ofstream fiche_ome_local (name_ome_local) ;
fiche_ome_local.precision(8) ;
bool loop = true ;
bool search_masses = false ;
double ent_c = ent_c_init ;
Cmp shift_bh_old (hole.mp) ;
Cmp shift_ns_old (star.mp) ;
double erreur = 1 ;
int conte = 0 ;
double old_mass_ns ;
mass_ns = star.mass_b() ;
double spin_old ;
double spin = 1;
bool adapt = true ;
while (loop) {
spin_old = spin ;
spin = hole.local_momentum() ;
if (sortie !=0) {
fiche_ome_local << conte << " " << hole.omega_local << endl ;
fiche_spin << conte << " " << spin/m1/m1 << endl ;
}
double conv_spin = fabs(spin-spin_old) ;
double error_spin = spin - spin_cible ;
double rel_diff_spin = (spin_cible==0) ? fabs(error_spin) :
fabs(error_spin)/spin_cible ;
if ((conv_spin*2<rel_diff_spin) && (search_masses) &&
hole.get_rot_state() != COROT) {
double func = scale_ome_local*log((2+error_spin)/(2+2*error_spin));
hole.set_omega_local(hole.omega_local+func) ;
}
old_mass_ns = mass_ns ;
if (hole.get_shift_auto().get_etat() != ETATZERO)
shift_bh_old = hole.get_shift_auto()(0) ;
else
shift_bh_old = 0 ;
if (star.get_shift_auto().get_etat() != ETATZERO)
shift_ns_old = star.get_shift_auto()(0) ;
else
shift_ns_old = 0 ;
star.kinematics(omega, x_axe) ;
star.fait_d_psi() ;
star.hydro_euler() ;
Tbl diff (7) ;
diff.set_etat_qcq() ;
int ite ;
star.equilibrium_nsbh (adapt, ent_c, ite, itemax_equil, itemax_mp_et,
relax, itemax_mp_et, relax, diff) ;
hole.update_metric(star) ;
hole.equilibrium (star, precis, relax, bound_nn, lim_nn) ;
cout << "Apres equilibrium" << endl ;
star.update_metric(hole) ;
cout << "Apres star::update_metric" << endl ;
star.update_metric_der_comp(hole) ;
cout << "Apres star::update_metric_der_comp" << endl ;
fait_tkij(bound_nn, lim_nn) ;
cout << "Apres Bin_ns_bh::fait_tkij" << endl ;
erreur = 0 ;
Tbl diff_bh (diffrelmax (shift_bh_old, hole.get_shift_auto()(0))) ;
for (int i=1 ; i<nz_bh ; i++)
if (diff_bh(i) > erreur)
erreur = diff_bh(i) ;
Tbl diff_ns (diffrelmax (shift_ns_old, star.get_shift_auto()(0))) ;
for (int i=0 ; i<nz_ns ; i++)
if (diff_ns(i) > erreur)
erreur = diff_ns(i) ;
if (erreur<seuil_masses)
search_masses = true ;
mass_ns = star.mass_b() ;
cout << "Avant viriel" << endl ;
double error_viriel = viriel() ;
cout << "Apres viriel" << endl ;
double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
cout << "Apres linear" << endl ;
double error_m1 = 1.-sqrt(hole.area()/16./M_PI)/m1 ;
cout << "Apres Mbh" << endl ;
double error_m2 = 1 - mass_ns/m2 ;
cout << "Apres Mns" << endl ;
if (sortie != 0) {
fiche_iteration << conte << " " << erreur << endl ;
fiche_correction << conte << " " << hole.regul << endl ;
fiche_viriel << conte << " " << error_viriel << endl ;
fiche_linear << conte << " " << error_linear << endl ;
fiche_error_m1 << conte << " " << error_m1 << endl ;
fiche_error_m2 << conte << " " << error_m2 << endl ;
fiche_hor << conte << " " << hole.mp.val_r(0, 1, 0,0) << endl ;
fiche_entc << conte << " " << ent_c << endl ;
fiche_dist << conte << " " << distance << endl ;
fiche_ome << conte << " " << omega << endl ;
fiche_axe << conte << " " << axe_pos << endl ;
}
// The axis position
double scaling_axe = pow((2+error_linear)/(2+2*error_linear), 0.1) ;
axe_pos *= scaling_axe ;
star.set_mp().set_ori (axe_pos, 0, 0) ;
hole.set_mp().set_ori (-distance + axe_pos, 0, 0) ;
// Value of omega
double new_ome = star.compute_angul() ;
if (new_ome !=0) {
set_omega(new_ome) ;
if (hole.get_rot_state() == COROT)
hole.set_omega_local(new_ome) ;
}
// The right distance
double error_dist = (distance-dist)/dist ;
double scale_d = pow((2+error_dist)/(2+2*error_dist), 0.2) ;
distance *= scale_d ;
// Converge to the right masses :
if (search_masses) {
double scaling_r = pow((2-error_m1)/(2-2*error_m1), 0.25) ;
hole.mp.homothetie_interne(scaling_r) ;
hole.set_rayon(hole.get_rayon()*scaling_r) ;
}
// The case of the NS :
double convergence = fabs(mass_ns - old_mass_ns)/mass_ns ;
double rel_diff = fabs(error_m2) ;
if ((search_masses) && (convergence*2 < rel_diff)) {
double scaling_ent = pow((2-error_m2)/(2-2*error_m2), 1) ;
ent_c *= scaling_ent ;
}
cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
//if (erreur < 1e-4)
//adapt = false ;
if (erreur < precis)
loop = false ;
conte ++ ;
}
fiche_iteration.close() ;
fiche_correction.close() ;
fiche_viriel.close() ;
fiche_ome.close() ;
fiche_linear.close() ;
fiche_axe.close() ;
fiche_error_m1.close() ;
fiche_error_m2.close() ;
fiche_hor.close() ;
fiche_entc.close() ;
fiche_dist.close() ;
fiche_spin.close() ;
fiche_ome_local.close() ;
}
}
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::Bin_ns_bh::set_omega
void set_omega(double)
Sets the orbital angular velocity [{\tt f_unit}].
Definition bin_ns_bh.C:218
Lorene::Bin_ns_bh::x_axe
double x_axe
Absolute X coordinate of the rotation axis.
Definition bin_ns_bh.h:140
Lorene::Bin_ns_bh::star
Et_bin_nsbh star
The neutron star.
Definition bin_ns_bh.h:128
Lorene::Bin_ns_bh::hole
Bhole hole
The black hole.
Definition bin_ns_bh.h:131
Lorene::Bin_ns_bh::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition bin_ns_bh.h:136
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::diffrelmax
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::adapt
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.