et_bin_bhns_extr_eq_ylm.C
/*
* Method of class Et_bin_bhns_extr to compute an equilibrium configuration
* of a BH-NS binary system with an extreme mass ratio
*
* (see file et_bin_bhns_extr.h for documentation).
*
*/
/*
* Copyright (c) 2004-2005 Keisuke Taniguchi
* 2004 Joshua A. Faber
*
* 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 version 2
* as published by the Free Software Foundation.
*
* 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 et_bin_bhns_extr_eq_ylm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_eq_ylm.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $" ;
/*
* $Id: et_bin_bhns_extr_eq_ylm.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $
* $Log: et_bin_bhns_extr_eq_ylm.C,v $
* Revision 1.7 2014/10/13 08:52:54 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2014/10/06 15:13:07 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.5 2005/02/28 23:12:28 k_taniguchi
* Modification to include the case of the conformally flat background metric
*
* Revision 1.4 2005/01/31 20:28:14 k_taniguchi
* Addition of a number of coefficients which are deleted in the filtre_phi
* in the argument.
*
* Revision 1.3 2005/01/25 17:42:02 k_taniguchi
* Suppression of the filter for the source term of the shift vector.
*
* Revision 1.2 2005/01/03 18:03:39 k_taniguchi
* Addition of the method to fix the position of the neutron star
* in the coordinate system.
*
* Revision 1.1 2004/12/29 16:31:24 k_taniguchi
* *** empty log message ***
*
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_eq_ylm.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $
*
*/
// C headers
#include <cmath>
// Lorene headers
#include "et_bin_bhns_extr.h"
#include "etoile.h"
#include "map.h"
#include "coord.h"
#include "param.h"
#include "eos.h"
#include "graphique.h"
#include "utilitaires.h"
#include "unites.h"
namespace Lorene {
void Et_bin_bhns_extr::equil_bhns_extr_ylm_ks(double ent_c,
const double& mass,
const double& sepa,
double* nu_int,
double* beta_int,
double* shift_int,
int mermax,
int mermax_poisson,
double relax_poisson,
double relax_ylm,
int mermax_potvit,
double relax_potvit,
int np_filter,
double thres_adapt,
Tbl& diff) {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
assert( kerrschild == true ) ;
// Initializations
// ---------------
const Mg3d* mg = mp.get_mg() ;
int nz = mg->get_nzone() ; // total number of domains
// The following is required to initialize mp_prev as a Map_et:
Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
// Domain and radial indices of points at the surface of the star:
int l_b = nzet - 1 ;
int i_b = mg->get_nr(l_b) - 1 ;
int k_b ;
int j_b ;
// Value of the enthalpy defining the surface of the star
double ent_b = 0 ;
// Error indicators
// ----------------
double& diff_ent = diff.set(0) ;
double& diff_vel_pot = diff.set(1) ;
double& diff_logn = diff.set(2) ;
double& diff_beta = diff.set(3) ;
double& diff_shift_x = diff.set(4) ;
double& diff_shift_y = diff.set(5) ;
double& diff_shift_z = diff.set(6) ;
// Parameters for the function Map_et::adapt
// -----------------------------------------
Param par_adapt ;
int nitermax = 100 ;
int niter ;
int adapt_flag = 1 ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
int nz_search = nzet ; // Number of domains for searching
// the enthalpy isosurfaces
double precis_secant = 1.e-14 ;
double alpha_r ;
double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
Tbl ent_limit(nz) ;
par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
// locate zeros by the secant method
par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
// to the isosurfaces of ent is to be
// performed
par_adapt.add_int(nz_search, 2) ; // number of domains to search for
// the enthalpy isosurface
par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
// used in the secant method
par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
// the determination of zeros by
// the secant method
par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
// 0 = contracting mapping
par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
// distances will be multiplied
par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
// to define the isosurfaces
// Parameters for the function Map_et::poisson for logn_auto
// ---------------------------------------------------------
double precis_poisson = 1.e-16 ;
Param par_poisson1 ;
par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson1.add_double(precis_poisson, 1) ; // required precision
par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
// used
par_poisson1.add_cmp_mod( ssjm1_logn ) ;
// Parameters for the function Map_et::poisson for beta_auto
// ---------------------------------------------------------
Param par_poisson2 ;
par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson2.add_double(precis_poisson, 1) ; // required precision
par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
// used
par_poisson2.add_cmp_mod( ssjm1_beta ) ;
// Parameters for the function Tenseur::poisson_vect
// -------------------------------------------------
Param par_poisson_vect ;
par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
// iterations
par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
par_poisson_vect.add_int_mod(niter, 0) ;
// External potential
// See Eq (99) from Gourgoulhon et al. (2001)
// -----------------------------------------
Tenseur pot_ext = logn_comp + pot_centri + loggam ;
Tenseur ent_jm1 = ent ; // Enthalpy at previous step
Tenseur source(mp) ; // source term in the equation for logn_auto
// and beta_auto
Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
// equation for shift_auto
// maximum l to evaluate multipole moments y_lm
int l_max = 4 ;
if (l_max > 6) {
cout << "We don't have those ylm programmed yet!!!!" << endl ;
abort() ;
}
int nylm = (l_max+1) * (l_max+1) ;
Cmp** ylmvec = new Cmp* [nylm] ;
//==========================================================//
// Start of iteration //
//==========================================================//
for(int mer=0 ; mer<mermax ; mer++ ) {
cout << "-----------------------------------------------" << endl ;
cout << "step: " << mer << endl ;
cout << "diff_ent = " << diff_ent << endl ;
//------------------------------------------------------
// Resolution of the elliptic equation for the velocity
// scalar potential
//------------------------------------------------------
if (irrotational) {
diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
precis_poisson, relax_potvit) ;
}
//-------------------------------------
// Computation of the new radial scale
//-------------------------------------
// alpha_r (r = alpha_r r') is determined so that the enthalpy
// takes the requested value ent_b at the stellar surface
// Values at the center of the star:
double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
// Search for the reference point (theta_*, phi_*)
// [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
// at the surface of the star
// ------------------------------------------------------------------
double alpha_r2 = 0 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
// See Eq (100) from Gourgoulhon et al. (2001)
double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
/ ( logn_auto_b - logn_auto_c ) ;
if (alpha_r2_jk > alpha_r2) {
alpha_r2 = alpha_r2_jk ;
k_b = k ;
j_b = j ;
}
}
}
alpha_r = sqrt(alpha_r2) ;
cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
<< alpha_r << endl ;
// New value of logn_auto
// ----------------------
logn_auto = alpha_r2 * logn_auto ;
logn_auto_regu = alpha_r2 * logn_auto_regu ;
logn_auto_c = logn_auto()(0, 0, 0, 0) ;
//------------------------------------------------------------
// Change the values of the inner points of the second domain
// by those of the outer points of the first domain
//------------------------------------------------------------
(logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
//--------------------------------------------
// First integral --> enthalpy in all space
// See Eq (98) from Gourgoulhon et al. (2001)
//--------------------------------------------
ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
//----------------------------------------------------------
// Change the enthalpy field to be set its maximum position
// at the coordinate center
//----------------------------------------------------------
double dentdx = ent().dsdx()(0, 0, 0, 0) ;
double dentdy = ent().dsdy()(0, 0, 0, 0) ;
cout << "dH/dx|_center = " << dentdx << endl ;
cout << "dH/dy|_center = " << dentdy << endl ;
double dec_fact = 1. ;
Tenseur func_in(mp) ;
func_in.set_etat_qcq() ;
func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x
- dec_fact * (dentdy/ent_c) * mp.y ;
func_in.set().annule(nzet, nz-1) ;
func_in.set_std_base() ;
Tenseur func_ex(mp) ;
func_ex.set_etat_qcq() ;
func_ex.set() = 1. ;
func_ex.set().annule(0, nzet-1) ;
func_ex.set_std_base() ;
// New enthalpy field
// ------------------
ent.set() = ent() * (func_in() + func_ex()) ;
(ent().va).smooth(nzet, (ent.set()).va) ;
double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
double dentdy_new = ent().dsdy()(0, 0, 0, 0) ;
cout << "dH/dx|_new = " << dentdx_new << endl ;
cout << "dH/dy|_new = " << dentdy_new << endl ;
//-----------------------------------------------------
// Adaptation of the mapping to the new enthalpy field
//----------------------------------------------------
// Shall the adaptation be performed (cusp) ?
// ------------------------------------------
double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
double rap_dent = fabs( dent_eq / dent_pole ) ;
cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
if ( rap_dent < thres_adapt ) {
adapt_flag = 0 ; // No adaptation of the mapping
cout << "******* FROZEN MAPPING *********" << endl ;
}
else{
adapt_flag = 1 ; // The adaptation of the mapping is to be
// performed
}
ent_limit.set_etat_qcq() ;
for (int l=0; l<nzet; l++) { // loop on domains inside the star
ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
}
ent_limit.set(nzet-1) = ent_b ;
Map_et mp_prev = mp_et ;
mp.adapt(ent(), par_adapt) ;
//----------------------------------------------------
// Computation of the enthalpy at the new grid points
//----------------------------------------------------
mp_prev.homothetie(alpha_r) ;
mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
double fact_resize = 1. / alpha_r ;
for (int l=nzet; l<nz-1; l++) {
mp_et.resize(l, fact_resize) ;
}
mp_et.resize_extr(fact_resize) ;
double ent_s_max = -1 ;
int k_s_max = -1 ;
int j_s_max = -1 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double xx = fabs( ent()(l_b, k, j, i_b) ) ;
if (xx > ent_s_max) {
ent_s_max = xx ;
k_s_max = k ;
j_s_max = j ;
}
}
}
cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
<< " and nzet : " << endl ;
cout << " " << ent_s_max << " reached for k = " << k_s_max
<< " and j = " << j_s_max << endl ;
//******************************************
// Call this each time the mapping changes,
// but no more often than that
//******************************************
int oldindex = -1 ;
for (int l=0; l<=l_max; l++) {
for (int m=0; m<= l; m++) {
int index = l*l + 2*m ;
if(m > 0) index -= 1 ;
if(index != oldindex+1) {
cout << "error!! " << l << " " << m << " " << index
<< " " << oldindex << endl ;
abort() ;
}
// set up Cmp's in ylmvec with proper parity
if ((l+m) % 2 == 0) {
ylmvec[index] = new Cmp (logn_auto()) ;
} else {
ylmvec[index] = new Cmp (shift_auto(2)) ;
}
oldindex = index ;
if (m > 0) {
index += 1 ;
if((l+m) % 2 == 0) {
ylmvec[index] = new Cmp (logn_auto()) ;
} else {
ylmvec[index] = new Cmp (shift_auto(2)) ;
}
oldindex = index ;
}
}
}
if(oldindex+1 != nylm) {
cout << "ERROR!!! " << oldindex << " "<< nylm << endl ;
abort() ;
}
// This may need to be in some class definition,
// it references the relevant map_et
get_ylm(nylm, ylmvec);
//-------------------
// Equation of state
//-------------------
equation_of_state() ; // computes new values for nbar (n), ener (e)
// and press (p) from the new ent (H)
//----------------------------------------------------------
// Matter source terms in the gravitational field equations
//---------------------------------------------------------
hydro_euler_extr(mass, sepa) ; // computes new values for
// ener_euler (E), s_euler (S)
// and u_euler (U^i)
//----------------------------------------------------
// Auxiliary terms for the source of Poisson equation
//----------------------------------------------------
const Coord& xx = mp.x ;
const Coord& yy = mp.y ;
const Coord& zz = mp.z ;
Tenseur r_bh(mp) ;
r_bh.set_etat_qcq() ;
r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
r_bh.set_std_base() ;
Tenseur xx_cov(mp, 1, COV, ref_triad) ;
xx_cov.set_etat_qcq() ;
xx_cov.set(0) = xx + sepa ;
xx_cov.set(1) = yy ;
xx_cov.set(2) = zz ;
xx_cov.set_std_base() ;
Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
xsr_cov = xx_cov / r_bh ;
xsr_cov.set_std_base() ;
Tenseur msr(mp) ;
msr = ggrav * mass / r_bh ;
msr.set_std_base() ;
Tenseur lapse_bh(mp) ;
lapse_bh = 1. / sqrt( 1.+2.*msr ) ;
lapse_bh.set_std_base() ;
Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
lapse_bh2 = 1. / (1.+2.*msr) ;
lapse_bh2.set_std_base() ;
Tenseur ldnu(mp) ;
ldnu = flat_scalar_prod_desal(xsr_cov, d_logn_auto) ;
ldnu.set_std_base() ;
Tenseur ldbeta(mp) ;
ldbeta = flat_scalar_prod_desal(xsr_cov, d_beta_auto) ;
ldbeta.set_std_base() ;
Tenseur lltkij(mp) ;
lltkij.set_etat_qcq() ;
lltkij.set() = 0 ;
lltkij.set_std_base() ;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
lltkij.set() += xsr_cov(i) * xsr_cov(j) * tkij_auto(i, j) ;
Tenseur lshift(mp) ;
lshift = flat_scalar_prod_desal(xsr_cov, shift_auto) ;
lshift.set_std_base() ;
Tenseur d_ldnu(mp, 1, COV, ref_triad) ;
d_ldnu = ldnu.gradient() ; // (d/dx, d/dy, d/dz)
d_ldnu.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
Tenseur ldldnu(mp) ;
ldldnu = flat_scalar_prod_desal(xsr_cov, d_ldnu) ;
ldldnu.set_std_base() ;
Tenseur d_ldbeta(mp, 1, COV, ref_triad) ;
d_ldbeta = ldbeta.gradient() ; // (d/dx, d/dy, d/dz)
d_ldbeta.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
Tenseur ldldbeta(mp) ;
ldldbeta = flat_scalar_prod_desal(xsr_cov, d_ldbeta) ;
ldldbeta.set_std_base() ;
//------------------------------------------
// Poisson equation for logn_auto (nu_auto)
//------------------------------------------
// Source
// ------
if (relativistic) {
source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
- flat_scalar_prod_desal(d_logn_auto, d_beta_auto)
+ 2.*lapse_bh2 % msr % (ldnu % ldbeta + ldldnu)
+ lapse_bh2 % lapse_bh2 % msr % (2.*(ldnu + 4.*msr % ldnu)
- ldbeta) / r_bh
- (4.*a_car % lapse_bh2 % lapse_bh2 % msr / 3. / nnn / r_bh)
* (2.+3.*msr) * (3.+4.*msr) * lltkij
+ (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
/ nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
+ (4.*pow(lapse_bh2, 3.) % msr % msr / 3. / r_bh / r_bh)
% (2.*(a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
+ (a_car - 1.) % pow(1.+3.*msr, 2.)
- 3.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
}
else {
cout << "The computation of BH-NS binary systems"
<< " should be relativistic !!!" << endl ;
abort() ;
}
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
double* intvec_logn = new double[nylm] ;
// get_integrals needs access to map_et
get_integrals(nylm, intvec_logn, ylmvec, source()) ;
// relax if you have a previously calculated set of moments
if(nu_int[0] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_logn[i] *= relax_ylm ;
intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
}
}
source().poisson_ylm(par_poisson1, logn_auto.set(), nylm,
intvec_logn) ;
for (int i=0; i<nylm; i++) {
nu_int[i] = intvec_logn[i] ;
}
delete [] intvec_logn ;
// Construct logn_auto_regu for et_bin_upmetr_extr.C
// -------------------------------------------------
logn_auto_regu = logn_auto ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
cout <<
"Relative error in the resolution of the equation for logn_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_logn(l) << " " ;
}
cout << endl ;
diff_logn = max(abs(tdiff_logn)) ;
if (relativistic) {
//--------------------------------
// Poisson equation for beta_auto
//--------------------------------
// Source
// ------
source = qpig * a_car % s_euler + 0.75 * akcar_auto
- 0.5 * flat_scalar_prod_desal(d_logn_auto, d_logn_auto)
- 0.5 * flat_scalar_prod_desal(d_beta_auto, d_beta_auto)
+ lapse_bh2 % msr % (ldnu%ldnu + ldbeta%ldbeta + 2.*ldldbeta)
+ lapse_bh2 % lapse_bh2 % msr % (2.*(1.+4.*msr) * ldbeta
- ldnu) / r_bh
- (a_car % lapse_bh2 % lapse_bh2 % msr / nnn / r_bh)
* (2.+3.*msr) * (3.+4.*msr) * lltkij
+ (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
/ nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
+ (2.*pow(lapse_bh2, 3.) % msr % msr / r_bh / r_bh)
% ((a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
+ (a_car - 1.) * pow(1.+3.*msr, 2.)
- 2.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
double* intvec_beta = new double[nylm] ;
// get_integrals needs access to map_et
get_integrals(nylm, intvec_beta, ylmvec, source()) ;
// relax if you have a previously calculated set of moments
if(beta_int[0] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_beta[i] *= relax_ylm ;
intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
}
}
source().poisson_ylm(par_poisson2, beta_auto.set(),
nylm, intvec_beta) ;
for (int i=0; i<nylm; i++) {
beta_int[i] = intvec_beta[i] ;
}
delete [] intvec_beta ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
cout << "Relative error in the resolution of the equation for "
<< "beta_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta(l) << " " ;
}
cout << endl ;
diff_beta = max(abs(tdiff_beta)) ;
//----------------------------------------
// Vector Poisson equation for shift_auto
//----------------------------------------
// Some auxiliary terms for the source
// -----------------------------------
Tenseur xx_con(mp, 1, CON, ref_triad) ;
xx_con.set_etat_qcq() ;
xx_con.set(0) = xx + sepa ;
xx_con.set(1) = yy ;
xx_con.set(2) = zz ;
xx_con.set_std_base() ;
Tenseur xsr_con(mp, 1, CON, ref_triad) ;
xsr_con = xx_con / r_bh ;
xsr_con.set_std_base() ;
// Components of shift_auto with respect to the Cartesian triad
// (d/dx, d/dy, d/dz) of the mapping :
Tenseur shift_auto_local = shift_auto ;
shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
// Gradient (partial derivatives with respect to the Cartesian
// coordinates of the mapping)
// dn(i, j) = D_i N^j
Tenseur dn = shift_auto_local.gradient() ;
// Return to the absolute reference frame
dn.change_triad(ref_triad) ;
// Trace of D_i N^j = divergence of N^j :
Tenseur divn = contract(dn, 0, 1) ;
// l^j D_j N^i
Tenseur ldn_con = contract(xsr_con, 0, dn, 0) ;
// D_j (l^k D_k N^i): dldn(j, i)
Tenseur ldn_local = ldn_con ;
ldn_local.change_triad( mp.get_bvect_cart() ) ;
Tenseur dldn = ldn_local.gradient() ;
dldn.change_triad(ref_triad) ;
// l^j D_j (l^k D_k N^i)
Tenseur ldldn = contract(xsr_con, 0, dldn, 0) ;
// l_k D_j N^k
Tenseur ldn_cov = contract(xsr_cov, 0, dn, 1) ;
// l^j l_k D_j N^k
Tenseur lldn_cov = contract(xsr_con, 0, ldn_cov, 0) ;
// eta^{ij} l_k D_j N^k
Tenseur eldn(mp, 1, CON, ref_triad) ;
eldn.set_etat_qcq() ;
eldn.set(0) = ldn_cov(0) ;
eldn.set(1) = ldn_cov(1) ;
eldn.set(2) = ldn_cov(2) ;
eldn.set_std_base() ;
// l^i D_j N^j
Tenseur ldivn = xsr_con % divn ;
// D_j (l^i D_k N^k): dldivn(j, i)
Tenseur ldivn_local = ldivn ;
ldivn_local.change_triad( mp.get_bvect_cart() ) ;
Tenseur dldivn = ldivn_local.gradient() ;
dldivn.change_triad(ref_triad) ;
// l^j D_j (l^i D_k N^k)
Tenseur ldldivn = contract(xsr_con, 0, dldivn, 0) ;
// l_j N^j
Tenseur ln = contract(xsr_cov, 0, shift_auto, 0) ;
Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
Tenseur lvtmp = contract(xsr_con, 0, vtmp, 0) ;
// eta^{ij} vtmp_j
Tenseur evtmp(mp, 1, CON, ref_triad) ;
evtmp.set_etat_qcq() ;
evtmp.set(0) = vtmp(0) ;
evtmp.set(1) = vtmp(1) ;
evtmp.set(2) = vtmp(2) ;
evtmp.set_std_base() ;
// lapse_ns
Tenseur lapse_ns(mp) ;
lapse_ns = exp(logn_auto) ;
lapse_ns.set_std_base() ;
// Source
// ------
source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
% u_euler
+ nnn % flat_scalar_prod_desal(tkij_auto, vtmp)
- 2.*nnn % lapse_bh2 * msr / r_bh
% flat_scalar_prod_desal(tkij_auto, xsr_cov)
+ 2.*lapse_bh2 * msr * (3.*ldldn + ldldivn) / 3.
- lapse_bh2 * msr / r_bh
* (4.*ldivn - lapse_bh2 % (3.*ldn_con + 8.*msr * ldn_con)
- (eldn + 2.*lapse_bh2*(9.+11.*msr)*lldn_cov%xsr_con) / 3.)
- 2.*lapse_bh2 % lapse_bh2 * msr / r_bh / r_bh
* ( (4.+11.*msr) * shift_auto
- lapse_bh2 * (12.+51.*msr+46.*msr*msr) * ln % xsr_con )
/ 3.
+ 8.*pow(lapse_bh2, 4.) * msr / r_bh / r_bh
% (lapse_ns - 1.) * (2.+10.*msr+9.*msr*msr) * xsr_con / 3.
+ 2.*pow(lapse_bh2, 3.) * msr / r_bh * (2.+3.*msr)
* ( (1.+2.*msr) * evtmp - (3.+2.*msr) * lvtmp * xsr_con) / 3. ;
source_shift.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
// Filter for the source of shift vector :
for (int i=0; i<3; i++) {
for (int l=0; l<nz; l++) {
if (source_shift(i).get_etat() != ETATZERO)
source_shift.set(i).filtre_phi(np_filter, l) ;
}
}
// For Tenseur::poisson_vect, the triad must be the mapping
// triad, not the reference one:
source_shift.change_triad( mp.get_bvect_cart() ) ;
/*
for (int i=0; i<3; i++) {
if(source_shift(i).dz_nonzero()) {
assert( source_shift(i).get_dzpuis() == 4 ) ;
}
else {
(source_shift.set(i)).set_dzpuis(4) ;
}
}
source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
*/
double lambda_shift = double(1) / double(3) ;
double* intvec_shift = new double[4*nylm] ;
for (int i=0; i<3; i++) {
double* intvec = new double[nylm];
get_integrals(nylm, intvec, ylmvec, source_shift(i));
for (int j=0; j<nylm; j++) {
intvec_shift[i*nylm+j] = intvec[j] ;
}
if(shift_int[i*nylm] != -1.0) {
for (int j=0; j<nylm; j++) {
intvec_shift[i*nylm+j] *= relax_ylm ;
intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
* shift_int[i*nylm+j] ;
}
}
delete [] intvec ;
}
Tenseur source_scal (-skxk(source_shift)) ;
Cmp soscal (source_scal()) ;
double* intvec2 = new double[nylm] ;
get_integrals(nylm, intvec2, ylmvec, soscal) ;
for (int j=0; j<nylm; j++) {
intvec_shift[3*nylm+j] = intvec2[j] ;
}
if(shift_int[3*nylm] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_shift[3*nylm+i] *= relax_ylm ;
intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
*shift_int[3*nylm+i] ;
}
}
delete [] intvec2 ;
source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
shift_auto, w_shift,
khi_shift, nylm, intvec_shift) ;
for (int i=0; i<4*nylm; i++) {
shift_int[i] = intvec_shift[i] ;
}
delete[] intvec_shift ;
// Check: has the equation for shift_auto been correctly solved ?
// --------------------------------------------------------------
// Divergence of shift_auto :
Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
// divna.dec2_dzpuis() ; // dzpuis 2 -> 0
// Grad(div) :
Tenseur graddivn = divna.gradient() ;
// graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
// Full operator :
Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
lap_shift.set_etat_qcq() ;
for (int i=0; i<3; i++) {
lap_shift.set(i) = shift_auto(i).laplacien()
+ lambda_shift * graddivn(i) ;
}
Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
cout << "Relative error in the resolution of the equation for "
<< "shift_auto : " << endl ;
cout << "x component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_x(l) << " " ;
}
cout << endl ;
cout << "y component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_y(l) << " " ;
}
cout << endl ;
cout << "z component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_z(l) << " " ;
}
cout << endl ;
diff_shift_x = max(abs(tdiff_shift_x)) ;
diff_shift_y = max(abs(tdiff_shift_y)) ;
diff_shift_z = max(abs(tdiff_shift_z)) ;
// Final result
// ------------
// The output of Tenseur::poisson_vect_falloff is on the mapping
// triad, it should therefore be transformed to components on the
// reference triad :
shift_auto.change_triad( ref_triad ) ;
} // End of relativistic equations
//------------------------------
// Relative change in enthalpy
//------------------------------
Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
diff_ent = diff_ent_tbl(0) ;
for (int l=1; l<nzet; l++) {
diff_ent += diff_ent_tbl(l) ;
}
diff_ent /= nzet ;
ent_jm1 = ent ;
for (int i=0; i<nylm; i++) {
delete ylmvec[i];
}
} // End of main loop
//========================================================//
// End of iteration //
//========================================================//
delete[] ylmvec ;
}
void Et_bin_bhns_extr::equil_bhns_extr_ylm_cf(double ent_c,
const double& mass,
const double& sepa,
double* nu_int,
double* beta_int,
double* shift_int,
int mermax,
int mermax_poisson,
double relax_poisson,
double relax_ylm,
int mermax_potvit,
double relax_potvit,
int np_filter,
double thres_adapt,
Tbl& diff) {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
assert( kerrschild == false ) ;
// Initializations
// ---------------
const Mg3d* mg = mp.get_mg() ;
int nz = mg->get_nzone() ; // total number of domains
// The following is required to initialize mp_prev as a Map_et:
Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
// Domain and radial indices of points at the surface of the star:
int l_b = nzet - 1 ;
int i_b = mg->get_nr(l_b) - 1 ;
int k_b ;
int j_b ;
// Value of the enthalpy defining the surface of the star
double ent_b = 0 ;
// Error indicators
// ----------------
double& diff_ent = diff.set(0) ;
double& diff_vel_pot = diff.set(1) ;
double& diff_logn = diff.set(2) ;
double& diff_beta = diff.set(3) ;
double& diff_shift_x = diff.set(4) ;
double& diff_shift_y = diff.set(5) ;
double& diff_shift_z = diff.set(6) ;
// Parameters for the function Map_et::adapt
// -----------------------------------------
Param par_adapt ;
int nitermax = 100 ;
int niter ;
int adapt_flag = 1 ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
int nz_search = nzet ; // Number of domains for searching
// the enthalpy isosurfaces
double precis_secant = 1.e-14 ;
double alpha_r ;
double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
Tbl ent_limit(nz) ;
par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
// locate zeros by the secant method
par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
// to the isosurfaces of ent is to be
// performed
par_adapt.add_int(nz_search, 2) ; // number of domains to search for
// the enthalpy isosurface
par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
// 0 = performs only the rescaling by
// the factor alpha_r
par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
// (theta_*, phi_*)
par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
// used in the secant method
par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
// the determination of zeros by
// the secant method
par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
// 0 = contracting mapping
par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
// distances will be multiplied
par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
// to define the isosurfaces
// Parameters for the function Map_et::poisson for logn_auto
// ---------------------------------------------------------
double precis_poisson = 1.e-16 ;
Param par_poisson1 ;
par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson1.add_double(precis_poisson, 1) ; // required precision
par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
// used
par_poisson1.add_cmp_mod( ssjm1_logn ) ;
// Parameters for the function Map_et::poisson for beta_auto
// ---------------------------------------------------------
Param par_poisson2 ;
par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson2.add_double(precis_poisson, 1) ; // required precision
par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
// used
par_poisson2.add_cmp_mod( ssjm1_beta ) ;
// Parameters for the function Tenseur::poisson_vect
// -------------------------------------------------
Param par_poisson_vect ;
par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
// iterations
par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
par_poisson_vect.add_int_mod(niter, 0) ;
// External potential
// See Eq (99) from Gourgoulhon et al. (2001)
// -----------------------------------------
Tenseur pot_ext = logn_comp + pot_centri + loggam ;
Tenseur ent_jm1 = ent ; // Enthalpy at previous step
Tenseur source(mp) ; // source term in the equation for logn_auto
// and beta_auto
Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
// equation for shift_auto
// maximum l to evaluate multipole moments y_lm
int l_max = 4 ;
if (l_max > 6) {
cout << "We don't have those ylm programmed yet!!!!" << endl ;
abort() ;
}
int nylm = (l_max+1) * (l_max+1) ;
Cmp** ylmvec = new Cmp* [nylm] ;
//==========================================================//
// Start of iteration //
//==========================================================//
for(int mer=0 ; mer<mermax ; mer++ ) {
cout << "-----------------------------------------------" << endl ;
cout << "step: " << mer << endl ;
cout << "diff_ent = " << diff_ent << endl ;
//------------------------------------------------------
// Resolution of the elliptic equation for the velocity
// scalar potential
//------------------------------------------------------
if (irrotational) {
diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
precis_poisson, relax_potvit) ;
}
//-------------------------------------
// Computation of the new radial scale
//-------------------------------------
// alpha_r (r = alpha_r r') is determined so that the enthalpy
// takes the requested value ent_b at the stellar surface
// Values at the center of the star:
double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
// Search for the reference point (theta_*, phi_*)
// [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
// at the surface of the star
// ------------------------------------------------------------------
double alpha_r2 = 0 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
// See Eq (100) from Gourgoulhon et al. (2001)
double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
/ ( logn_auto_b - logn_auto_c ) ;
if (alpha_r2_jk > alpha_r2) {
alpha_r2 = alpha_r2_jk ;
k_b = k ;
j_b = j ;
}
}
}
alpha_r = sqrt(alpha_r2) ;
cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
<< alpha_r << endl ;
// New value of logn_auto
// ----------------------
logn_auto = alpha_r2 * logn_auto ;
logn_auto_regu = alpha_r2 * logn_auto_regu ;
logn_auto_c = logn_auto()(0, 0, 0, 0) ;
//------------------------------------------------------------
// Change the values of the inner points of the second domain
// by those of the outer points of the first domain
//------------------------------------------------------------
(logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
//--------------------------------------------
// First integral --> enthalpy in all space
// See Eq (98) from Gourgoulhon et al. (2001)
//--------------------------------------------
ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
//---------------------------------------------------------
// Change the enthalpy field to accelerate the convergence
//---------------------------------------------------------
double dentdx = ent().dsdx()(0, 0, 0, 0) ;
double dentdy = ent().dsdy()(0, 0, 0, 0) ;
cout << "dH/dx|_center = " << dentdx << endl ;
cout << "dH/dy|_center = " << dentdy << endl ;
double dec_fact = 1. ;
Tenseur func_in(mp) ;
func_in.set_etat_qcq() ;
func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x ;
func_in.set().annule(nzet, nz-1) ;
func_in.set_std_base() ;
Tenseur func_ex(mp) ;
func_ex.set_etat_qcq() ;
func_ex.set() = 1. ;
func_ex.set().annule(0, nzet-1) ;
func_ex.set_std_base() ;
// New enthalpy field
// ------------------
ent.set() = ent() * (func_in() + func_ex()) ;
(ent().va).smooth(nzet, (ent.set()).va) ;
double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
cout << "dH/dx|_new = " << dentdx_new << endl ;
//-----------------------------------------------------
// Adaptation of the mapping to the new enthalpy field
//----------------------------------------------------
// Shall the adaptation be performed (cusp) ?
// ------------------------------------------
double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
double rap_dent = fabs( dent_eq / dent_pole ) ;
cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
if ( rap_dent < thres_adapt ) {
adapt_flag = 0 ; // No adaptation of the mapping
cout << "******* FROZEN MAPPING *********" << endl ;
}
else{
adapt_flag = 1 ; // The adaptation of the mapping is to be
// performed
}
ent_limit.set_etat_qcq() ;
for (int l=0; l<nzet; l++) { // loop on domains inside the star
ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
}
ent_limit.set(nzet-1) = ent_b ;
Map_et mp_prev = mp_et ;
mp.adapt(ent(), par_adapt) ;
//----------------------------------------------------
// Computation of the enthalpy at the new grid points
//----------------------------------------------------
mp_prev.homothetie(alpha_r) ;
mp.reevaluate_symy(&mp_prev, nzet+1, ent.set()) ;
double fact_resize = 1. / alpha_r ;
for (int l=nzet; l<nz-1; l++) {
mp_et.resize(l, fact_resize) ;
}
mp_et.resize_extr(fact_resize) ;
double ent_s_max = -1 ;
int k_s_max = -1 ;
int j_s_max = -1 ;
for (int k=0; k<mg->get_np(l_b); k++) {
for (int j=0; j<mg->get_nt(l_b); j++) {
double xx = fabs( ent()(l_b, k, j, i_b) ) ;
if (xx > ent_s_max) {
ent_s_max = xx ;
k_s_max = k ;
j_s_max = j ;
}
}
}
cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
<< " and nzet : " << endl ;
cout << " " << ent_s_max << " reached for k = " << k_s_max
<< " and j = " << j_s_max << endl ;
//******************************************
// Call this each time the mapping changes,
// but no more often than that
//******************************************
int oldindex = -1 ;
for (int l=0; l<=l_max; l++) {
for (int m=0; m<= l; m++) {
int index = l*l + 2*m ;
if(m > 0) index -= 1 ;
if(index != oldindex+1) {
cout << "error!! " << l << " " << m << " " << index
<< " " << oldindex << endl ;
abort() ;
}
// set up Cmp's in ylmvec with proper parity
if ((l+m) % 2 == 0) {
ylmvec[index] = new Cmp (logn_auto()) ;
} else {
ylmvec[index] = new Cmp (shift_auto(2)) ;
}
oldindex = index ;
if (m > 0) {
index += 1 ;
if((l+m) % 2 == 0) {
ylmvec[index] = new Cmp (logn_auto()) ;
} else {
ylmvec[index] = new Cmp (shift_auto(2)) ;
}
oldindex = index ;
}
}
}
if(oldindex+1 != nylm) {
cout << "ERROR!!! " << oldindex << " "<< nylm << endl ;
abort() ;
}
// This may need to be in some class definition,
// it references the relevant map_et
get_ylm(nylm, ylmvec);
//-------------------
// Equation of state
//-------------------
equation_of_state() ; // computes new values for nbar (n), ener (e)
// and press (p) from the new ent (H)
//----------------------------------------------------------
// Matter source terms in the gravitational field equations
//---------------------------------------------------------
hydro_euler_extr(mass, sepa) ; // computes new values for
// ener_euler (E), s_euler (S)
// and u_euler (U^i)
//----------------------------------------------------
// Auxiliary terms for the source of Poisson equation
//----------------------------------------------------
const Coord& xx = mp.x ;
const Coord& yy = mp.y ;
const Coord& zz = mp.z ;
Tenseur r_bh(mp) ;
r_bh.set_etat_qcq() ;
r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
r_bh.set_std_base() ;
Tenseur xx_cov(mp, 1, COV, ref_triad) ;
xx_cov.set_etat_qcq() ;
xx_cov.set(0) = xx + sepa ;
xx_cov.set(1) = yy ;
xx_cov.set(2) = zz ;
xx_cov.set_std_base() ;
Tenseur msr(mp) ;
msr = ggrav * mass / r_bh ;
msr.set_std_base() ;
Tenseur tmp(mp) ;
tmp = 1. / ( 1. - 0.25*msr*msr ) ;
tmp.set_std_base() ;
Tenseur xdnu(mp) ;
xdnu = flat_scalar_prod_desal(xx_cov, d_logn_auto) ;
xdnu.set_std_base() ;
Tenseur xdbeta(mp) ;
xdbeta = flat_scalar_prod_desal(xx_cov, d_beta_auto) ;
xdbeta.set_std_base() ;
//------------------------------------------
// Poisson equation for logn_auto (nu_auto)
//------------------------------------------
// Source
// ------
if (relativistic) {
source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
- flat_scalar_prod_desal(d_logn_auto, d_beta_auto)
- 0.5 * tmp % msr % msr % xdnu / r_bh / r_bh
- tmp % msr % xdbeta / r_bh / r_bh ;
}
else {
cout << "The computation of BH-NS binary systems"
<< " should be relativistic !!!" << endl ;
abort() ;
}
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
double* intvec_logn = new double[nylm] ;
// get_integrals needs access to map_et
get_integrals(nylm, intvec_logn, ylmvec, source()) ;
// relax if you have a previously calculated set of moments
if(nu_int[0] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_logn[i] *= relax_ylm ;
intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
}
}
source().poisson_ylm(par_poisson1, logn_auto.set(), nylm,
intvec_logn) ;
for (int i=0; i<nylm; i++) {
nu_int[i] = intvec_logn[i] ;
}
delete [] intvec_logn ;
// Construct logn_auto_regu for et_bin_upmetr_extr.C
// -------------------------------------------------
logn_auto_regu = logn_auto ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
cout <<
"Relative error in the resolution of the equation for logn_auto : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_logn(l) << " " ;
}
cout << endl ;
diff_logn = max(abs(tdiff_logn)) ;
if (relativistic) {
//--------------------------------
// Poisson equation for beta_auto
//--------------------------------
// Source
// ------
source = qpig * a_car % s_euler + 0.75 * akcar_auto
- 0.5 * flat_scalar_prod_desal(d_logn_auto, d_logn_auto)
- 0.5 * flat_scalar_prod_desal(d_beta_auto, d_beta_auto)
- tmp % msr % xdnu / r_bh / r_bh
- 0.5 * tmp % msr %msr % xdbeta / r_bh / r_bh ;
source.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
double* intvec_beta = new double[nylm] ;
// get_integrals needs access to map_et
get_integrals(nylm, intvec_beta, ylmvec, source()) ;
// relax if you have a previously calculated set of moments
if(beta_int[0] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_beta[i] *= relax_ylm ;
intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
}
}
source().poisson_ylm(par_poisson2, beta_auto.set(),
nylm, intvec_beta) ;
for (int i=0; i<nylm; i++) {
beta_int[i] = intvec_beta[i] ;
}
delete [] intvec_beta ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
cout << "Relative error in the resolution of the equation for "
<< "beta_auto : " << endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta(l) << " " ;
}
cout << endl ;
diff_beta = max(abs(tdiff_beta)) ;
//----------------------------------------
// Vector Poisson equation for shift_auto
//----------------------------------------
// Some auxiliary terms for the source
// -----------------------------------
Tenseur bhtmp(mp, 1, COV, ref_triad) ;
bhtmp.set_etat_qcq() ;
bhtmp = tmp % msr % (3.*msr-8.) % xx_cov / r_bh / r_bh ;
bhtmp.set_std_base() ;
Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
// Source
// ------
source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
% u_euler
+ nnn % flat_scalar_prod_desal(tkij_auto, vtmp+bhtmp) ;
source_shift.set_std_base() ;
// Resolution of the Poisson equation
// ----------------------------------
// Filter for the source of shift vector :
for (int i=0; i<3; i++) {
for (int l=0; l<nz; l++) {
if (source_shift(i).get_etat() != ETATZERO)
source_shift.set(i).filtre_phi(np_filter, l) ;
}
}
// For Tenseur::poisson_vect, the triad must be the mapping
// triad, not the reference one:
source_shift.change_triad( mp.get_bvect_cart() ) ;
/*
for (int i=0; i<3; i++) {
if(source_shift(i).dz_nonzero()) {
assert( source_shift(i).get_dzpuis() == 4 ) ;
}
else {
(source_shift.set(i)).set_dzpuis(4) ;
}
}
source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
*/
double lambda_shift = double(1) / double(3) ;
double* intvec_shift = new double[4*nylm] ;
for (int i=0; i<3; i++) {
double* intvec = new double[nylm];
get_integrals(nylm, intvec, ylmvec, source_shift(i));
for (int j=0; j<nylm; j++) {
intvec_shift[i*nylm+j] = intvec[j] ;
}
if(shift_int[i*nylm] != -1.0) {
for (int j=0; j<nylm; j++) {
intvec_shift[i*nylm+j] *= relax_ylm ;
intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
* shift_int[i*nylm+j] ;
}
}
delete [] intvec ;
}
Tenseur source_scal (-skxk(source_shift)) ;
Cmp soscal (source_scal()) ;
double* intvec2 = new double[nylm] ;
get_integrals(nylm, intvec2, ylmvec, soscal) ;
for (int j=0; j<nylm; j++) {
intvec_shift[3*nylm+j] = intvec2[j] ;
}
if(shift_int[3*nylm] != -1.0) {
for (int i=0; i<nylm; i++) {
intvec_shift[3*nylm+i] *= relax_ylm ;
intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
*shift_int[3*nylm+i] ;
}
}
delete [] intvec2 ;
source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
shift_auto, w_shift,
khi_shift, nylm, intvec_shift) ;
for (int i=0; i<4*nylm; i++) {
shift_int[i] = intvec_shift[i] ;
}
delete[] intvec_shift ;
// Check: has the equation for shift_auto been correctly solved ?
// --------------------------------------------------------------
// Divergence of shift_auto :
Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
// divna.dec2_dzpuis() ; // dzpuis 2 -> 0
// Grad(div) :
Tenseur graddivn = divna.gradient() ;
// graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
// Full operator :
Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
lap_shift.set_etat_qcq() ;
for (int i=0; i<3; i++) {
lap_shift.set(i) = shift_auto(i).laplacien()
+ lambda_shift * graddivn(i) ;
}
Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
cout << "Relative error in the resolution of the equation for "
<< "shift_auto : " << endl ;
cout << "x component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_x(l) << " " ;
}
cout << endl ;
cout << "y component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_y(l) << " " ;
}
cout << endl ;
cout << "z component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_shift_z(l) << " " ;
}
cout << endl ;
diff_shift_x = max(abs(tdiff_shift_x)) ;
diff_shift_y = max(abs(tdiff_shift_y)) ;
diff_shift_z = max(abs(tdiff_shift_z)) ;
// Final result
// ------------
// The output of Tenseur::poisson_vect_falloff is on the mapping
// triad, it should therefore be transformed to components on the
// reference triad :
shift_auto.change_triad( ref_triad ) ;
} // End of relativistic equations
//------------------------------
// Relative change in enthalpy
//------------------------------
Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
diff_ent = diff_ent_tbl(0) ;
for (int l=1; l<nzet; l++) {
diff_ent += diff_ent_tbl(l) ;
}
diff_ent /= nzet ;
ent_jm1 = ent ;
for (int i=0; i<nylm; i++) {
delete ylmvec[i];
}
} // End of main loop
//========================================================//
// End of iteration //
//========================================================//
delete[] ylmvec ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::filtre_phi
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
Definition cmp_manip.C:101
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Et_bin_bhns_extr::equil_bhns_extr_ylm_ks
void equil_bhns_extr_ylm_ks(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the Kerr-Schild background metric u...
Definition et_bin_bhns_extr_eq_ylm.C:79
Lorene::Et_bin_bhns_extr::kerrschild
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
Definition et_bin_bhns_extr.h:78
Lorene::Et_bin_bhns_extr::get_ylm
void get_ylm(int nylm, Cmp **ylmvec) const
Constructs spherical harmonics.
Definition et_bin_bhns_extr_ylm.C:63
Lorene::Et_bin_bhns_extr::get_integrals
void get_integrals(int nylm, double *intvec, Cmp **ylmvec, Cmp source) const
Computes multipole moments.
Definition et_bin_bhns_extr_ylm.C:255
Lorene::Et_bin_bhns_extr::velocity_pot_extr
double velocity_pot_extr(const double &mass, const double &sepa, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Definition et_bin_bhns_extr_velpot.C:67
Lorene::Et_bin_bhns_extr::hydro_euler_extr
void hydro_euler_extr(const double &mass, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition et_bin_bhns_extr_hydro.C:56
Lorene::Et_bin_bhns_extr::equil_bhns_extr_ylm_cf
void equil_bhns_extr_ylm_cf(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the conformally flat background met...
Definition et_bin_bhns_extr_eq_ylm.C:969
Lorene::Etoile_bin::d_beta_auto
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad ).
Definition etoile.h:879
Lorene::Etoile_bin::d_logn_auto
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad ).
Definition etoile.h:859
Lorene::Etoile_bin::ssjm1_khi
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:973
Lorene::Etoile_bin::ref_triad
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:828
Lorene::Etoile_bin::irrotational
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
Lorene::Etoile_bin::w_shift
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:908
Lorene::Etoile_bin::logn_comp
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Definition etoile.h:854
Lorene::Etoile_bin::ssjm1_wshift
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition etoile.h:983
Lorene::Etoile_bin::ssjm1_logn
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Definition etoile.h:959
Lorene::Etoile_bin::pot_centri
Tenseur pot_centri
Centrifugal potential.
Definition etoile.h:953
Lorene::Etoile_bin::loggam
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:849
Lorene::Etoile_bin::shift_auto
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition etoile.h:889
Lorene::Etoile_bin::tkij_auto
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition etoile.h:925
Lorene::Etoile_bin::akcar_auto
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition etoile.h:938
Lorene::Etoile_bin::ssjm1_beta
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Definition etoile.h:965
Lorene::Etoile_bin::khi_shift
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:918
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::logn_auto_regu
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:491
Lorene::Etoile::ray_eq_pi
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:256
Lorene::Etoile::nnn
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Lorene::Etoile::logn_auto
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:484
Lorene::Etoile::equation_of_state
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Lorene::Etoile::u_euler
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Etoile::press
Tenseur press
Fluid pressure.
Definition etoile.h:461
Lorene::Etoile::relativistic
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Lorene::Etoile::ener_euler
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Lorene::Etoile::s_euler
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Lorene::Etoile::ent
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case).
Definition etoile.h:457
Lorene::Etoile::beta_auto
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:506
Lorene::Etoile::a_car
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
Lorene::Etoile::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition etoile_global.C:415
Lorene::Map_et
Radial mapping of rather general form.
Definition map.h:2752
Lorene::Map_et::resize_extr
void resize_extr(double lambda)
Rescales the outer boundary of the outermost domain in the case of non-compactified external domain.
Definition map_et_resize_extr.C:55
Lorene::Map_et::homothetie
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
Lorene::Map_et::resize
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Definition map_et.C:928
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
Lorene::Param::add_cmp_mod
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition param.C:1004
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
Lorene::Param::add_tenseur_mod
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition param.C:1142
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Param::add_tbl
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition param.C:522
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates).
Definition tenseur.C:1542
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
Lorene::Tenseur::change_triad
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:668
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Lorene::skxk
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Definition tenseur_operateur.C:580
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::flat_scalar_prod_desal
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Definition tenseur_operateur.C:735
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::laplacien
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const =0
Computes the Laplacian of a scalar field.
Unites
Standard units of space, time and mass.