binaire_constr.C
/*
* Methods of class Binaire for estimating the error in the Hamiltionian
* and momentum constraints
*
* (see file binaire.h for documentation).
*/
/*
* Copyright (c) 2000-2001 Eric Gourgoulhon
*
* 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 binaire_constr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $" ;
/*
* $Id: binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
* $Log: binaire_constr.C,v $
* Revision 1.3 2014/10/13 08:52:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.2 2004/03/25 10:28:59 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
* LORENE
*
* Revision 2.1 2000/03/13 17:05:34 eric
* *** empty log message ***
*
* Revision 2.0 2000/03/13 14:26:08 eric
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "binaire.h"
#include "unites.h"
//----------------------------------------------//
// Hamiltonian constraint //
//----------------------------------------------//
namespace Lorene {
double Binaire::ham_constr() const {
using namespace Unites ;
if (p_ham_constr == 0x0) { // A new computation is required
Tenseur lap_alpha1( star1.get_mp() ) ;
Tenseur lap_alpha2( star2.get_mp() ) ;
Tenseur source1( star1.get_mp() ) ;
Tenseur source2( star2.get_mp() ) ;
Tenseur* p_lap_alpha[2] ;
Tenseur* p_source[2] ;
p_lap_alpha[0] = &lap_alpha1 ;
p_lap_alpha[1] = &lap_alpha2 ;
p_source[0] = &source1 ;
p_source[1] = &source2 ;
// Computation of the l.h.s. and r.h.s. of the Hamiltonian
// constraint in each star.
// -------------------------------------------------------
double som = 0 ;
for (int i=0; i<2; i++) {
// Laplacian of alpha = ln(A)
// --------------------------
Tenseur alpha_auto = et[i]->get_beta_auto()
- et[i]->get_logn_auto() ;
*(p_lap_alpha[i]) = alpha_auto().laplacien() ;
// Right-hand-side of the Hamiltonian constraint
// ---------------------------------------------
const Tenseur& a_car = et[i]->get_a_car() ;
const Tenseur& ener_euler = et[i]->get_ener_euler() ;
Tenseur d_alpha_auto = et[i]->get_d_beta_auto()
- et[i]->get_d_logn_auto() ;
Tenseur d_alpha_comp = et[i]->get_d_beta_comp()
- et[i]->get_d_logn_comp() ;
const Tenseur& akcar_auto = et[i]->get_akcar_auto() ;
const Tenseur& akcar_comp = et[i]->get_akcar_comp() ;
*(p_source[i]) = - qpig * a_car * ener_euler
- 0.25 * ( akcar_auto + akcar_comp )
- 0.5 *
( flat_scalar_prod(d_alpha_auto, d_alpha_auto)
+ flat_scalar_prod(d_alpha_auto, d_alpha_comp)
) ;
// Relative difference
// -------------------
Tbl diff = diffrel( (*(p_lap_alpha[i]))(), (*(p_source[i]))() ) ;
cout <<
"Binaire::ham_constr : relative difference Lap(alpha) <-> source : "
<< endl << diff << endl ;
som += max( abs(diff) ) ;
}
// Total error
// -----------
p_ham_constr = new double ;
*p_ham_constr = 0.5 * som ;
}
return *p_ham_constr ;
}
//----------------------------------------------//
// Momentum constraint //
//----------------------------------------------//
const Tbl& Binaire::mom_constr() const {
using namespace Unites ;
if (p_mom_constr == 0x0) { // A new computation is required
Tenseur divk1( star1.get_mp(), 1, CON, ref_triad ) ;
Tenseur divk2( star2.get_mp(), 1, CON, ref_triad ) ;
Tenseur source1( star1.get_mp(), 1, CON, ref_triad ) ;
Tenseur source2( star2.get_mp(), 1, CON, ref_triad ) ;
Tenseur* p_divk[2] ;
Tenseur* p_source[2] ;
p_divk[0] = &divk1 ;
p_divk[1] = &divk2 ;
p_source[0] = &source1 ;
p_source[1] = &source2 ;
// Computation of the l.h.s. and r.h.s. of the momentum
// constraint in each star.
// -------------------------------------------------------
double somx = 0 ;
double somy = 0 ;
double somz = 0 ;
for (int i=0; i<2; i++) {
// (flat space) divergence of K^{ij}
// ---------------------------------
const Tenseur& a_car = et[i]->get_a_car() ;
Tenseur kij_auto = et[i]->get_tkij_auto() / a_car ;
kij_auto.dec2_dzpuis() ; // dzpuis : 2 --> 0
// so that in the external domain, kij_auto
// contains now exactly K^{ij}
// The gradient of K^{ij} is computed on the local triad:
kij_auto.change_triad( (et[i]->get_mp()).get_bvect_cart() ) ;
*(p_divk[i]) = contract( kij_auto.gradient(), 0, 1) ;
// Back to the Reference triad :
p_divk[i]->change_triad( ref_triad ) ;
kij_auto.change_triad( ref_triad ) ;
// Right-hand-side of the momentum constraint
// ------------------------------------------
const Tenseur& u_euler = et[i]->get_u_euler() ;
const Tenseur& ener_euler = et[i]->get_ener_euler() ;
const Tenseur& press = et[i]->get_press() ;
Tenseur d_alpha = et[i]->get_d_beta_auto()
- et[i]->get_d_logn_auto()
+ et[i]->get_d_beta_comp()
- et[i]->get_d_logn_comp() ;
*(p_source[i]) = 2 * qpig * (ener_euler + press) * u_euler
- 5 * contract(kij_auto, 1, d_alpha, 0) ;
// Relative differences
// --------------------
Tbl diffx = diffrel( (*(p_divk[i]))(0), (*(p_source[i]))(0)) ;
Tbl diffy = diffrel( (*(p_divk[i]))(1), (*(p_source[i]))(1)) ;
Tbl diffz = diffrel( (*(p_divk[i]))(2), (*(p_source[i]))(2)) ;
cout << "Binaire::mom_constr : norme div(K) : " << endl ;
cout << "X component : " << norme( (*(p_divk[i]))(0) ) << endl ;
cout << "Y component : " << norme( (*(p_divk[i]))(1) ) << endl ;
cout << "Z component : " << norme( (*(p_divk[i]))(2) ) << endl ;
cout << "Binaire::mom_constr : norme source : " << endl ;
cout << "X component : " << norme( (*(p_source[i]))(0) ) << endl ;
cout << "Y component : " << norme( (*(p_source[i]))(1) ) << endl ;
cout << "Z component : " << norme( (*(p_source[i]))(2) ) << endl ;
cout <<
"Binaire::mom_constr : rel. diff. div(K) <-> source : "
<< endl ;
cout << "X component : " << diffx << endl ;
cout << "Y component : " << diffy << endl ;
cout << "Z component : " << diffz << endl ;
somx += max( abs(diffx) ) ;
somy += max( abs(diffy) ) ;
somz += max( abs(diffz) ) ;
}
// Total error
// -----------
p_mom_constr = new Tbl(3) ;
p_mom_constr->set_etat_qcq() ;
p_mom_constr->set(0) = 0.5 * somx ;
p_mom_constr->set(1) = 0.5 * somy ;
p_mom_constr->set(2) = 0.5 * somz ;
}
return *p_mom_constr ;
}
}
Lorene::Binaire::p_ham_constr
double * p_ham_constr
Relative error on the Hamiltonian constraint.
Definition binaire.h:160
Lorene::Binaire::et
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
Definition binaire.h:124
Lorene::Binaire::p_mom_constr
Tbl * p_mom_constr
Relative error on the momentum constraint.
Definition binaire.h:163
Lorene::Binaire::star2
Etoile_bin star2
Second star of the system.
Definition binaire.h:118
Lorene::Binaire::ref_triad
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition binaire.h:112
Lorene::Binaire::mom_constr
const Tbl & mom_constr() const
Estimates the relative error on the momentum constraint equation by comparing ${\overline\nabla}...
Definition binaire_constr.C:158
Lorene::Binaire::star1
Etoile_bin star1
First star of the system.
Definition binaire.h:115
Lorene::Binaire::ham_constr
double ham_constr() const
Estimates the relative error on the Hamiltonian constraint equation by comparing $\underline\Delta\ln...
Definition binaire_constr.C:69
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::dec2_dzpuis
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates).
Definition tenseur.C:1542
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::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Lorene::flat_scalar_prod
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Definition tenseur_operateur.C:653
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_bvect_cart
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
Unites
Standard units of space, time and mass.