binhor_equations.C
/*
* Copyright- (c) 2004-2005 Francois Limousin
* Jose-Luis Jaramillo
*
* 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 binhor_equations_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $" ;
/*
* $Id: binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $
* $Log: binhor_equations.C,v $
* Revision 1.21 2014/10/13 08:52:42 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.20 2014/10/06 15:13:01 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.19 2008/02/06 18:20:02 f_limousin
* Fixed an error in the triad
*
* Revision 1.18 2007/08/22 16:12:33 f_limousin
* Changed the name of the function computing \tilde{\gamma}_{ij}
*
* Revision 1.17 2007/04/13 15:28:55 f_limousin
* Lots of improvements, generalisation to an arbitrary state of
* rotation, implementation of the spatial metric given by Samaya.
*
* Revision 1.16 2006/08/01 14:37:19 f_limousin
* New version
*
* Revision 1.15 2006/06/29 08:51:00 f_limousin
* *** empty log message ***
*
* Revision 1.14 2006/05/24 16:56:37 f_limousin
* Many small modifs.
*
* Revision 1.13 2005/11/15 14:04:00 f_limousin
* Minor change to control the resolution of the equation for psi.
*
* Revision 1.12 2005/10/23 16:39:43 f_limousin
* Simplification of the equation in the case of a conformally
* flat metric and maximal slicing
*
* Revision 1.11 2005/09/13 18:33:15 f_limousin
* New function vv_bound_cart_bin(double) for computing binaries with
* berlin condition for the shift vector.
* Suppress all the symy and asymy in the importations.
*
* Revision 1.10 2005/07/11 08:21:57 f_limousin
* Implementation of a new boundary condition for the lapse in the binary
* case : boundary_nn_Dir_lapl().
*
* Revision 1.9 2005/06/09 16:12:04 f_limousin
* Implementation of amg_mom_adm().
*
* Revision 1.8 2005/04/29 14:02:44 f_limousin
* Important changes : manage the dependances between quantities (for
* instance psi and psi4). New function write_global(ost).
*
* Revision 1.7 2005/03/10 17:21:52 f_limousin
* Add the Berlin boundary condition for the shift.
* Some changes to avoid warnings.
*
* Revision 1.6 2005/03/03 10:29:02 f_limousin
* Delete omega in the parameters of the function boundary_beta_z().
*
* Revision 1.5 2005/02/24 17:25:23 f_limousin
* The boundary conditions for psi, N and beta are now parameters in
* par_init.d and par_coal.d.
*
* Revision 1.4 2005/02/11 18:21:38 f_limousin
* Dirichlet_binaire and neumann_binaire are taking Scalars as arguments
* instead of Cmps.
*
* Revision 1.3 2005/02/07 10:46:28 f_limousin
* Many changes !! The sources are written differently to minimize the
* numerical error, the boundary conditions are changed...
*
* Revision 1.2 2004/12/31 15:41:26 f_limousin
* Correction of an error
*
* Revision 1.1 2004/12/29 16:11:34 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $
*
*/
//standard
#include <cstdlib>
#include <cmath>
// Lorene
#include "nbr_spx.h"
#include "tensor.h"
#include "tenseur.h"
#include "isol_hor.h"
#include "proto.h"
#include "utilitaires.h"
//#include "graphique.h"
// Resolution for the lapse
// ------------------------
namespace Lorene {
void Bin_hor::solve_lapse (double precision, double relax, int bound_nn,
double lim_nn) {
assert ((relax >0) && (relax<=1)) ;
cout << "-----------------------------------------------" << endl ;
cout << "Resolution LAPSE" << endl ;
Scalar lapse_un_old (hole1.n_auto) ;
Scalar lapse_deux_old (hole2.n_auto) ;
Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
// Source 1
// --------
Scalar source_un (hole1.mp) ;
// Conformally flat
/*
source_un = hole1.get_psi4()*hole1.nn*aa_quad_un
-2.*contract(hole1.dn, 0, hole1.psi_auto
.derive_con(hole1.ff), 0)/hole1.psi ;
*/
source_un = hole1.get_psi4()*( hole1.nn*( aa_quad_un + 0.3333333333333333 *
hole1.trK*hole1.trK*hole1.decouple)
- hole1.trK_point*hole1.decouple )
-2.*contract(contract(hole1.hh, 0, hole1.n_auto
.derive_cov(hole1.ff), 0), 0, hole1.dpsi, 0)/hole1.psi
-2.*contract(hole1.dn, 0, hole1.psi_auto
.derive_con(hole1.ff), 0)/hole1.psi ;
- contract(hdirac1, 0, hole1.n_auto.derive_cov(hole1.ff), 0) ;
Scalar tmp_un (hole1.mp) ;
tmp_un = hole1.get_psi4()* contract(hole1.beta_auto, 0, hole1.trK.
derive_cov(hole1.ff), 0)
- contract( hole1.hh, 0, 1, hole1.n_auto.derive_cov(hole1.ff)
.derive_cov(hole1.ff), 0, 1 ) ;
tmp_un.inc_dzpuis() ; // dzpuis: 3 -> 4
source_un += tmp_un ;
// Source 2
// ---------
Scalar source_deux (hole2.mp) ;
// Conformally flat
/*
source_deux = hole2.get_psi4()*hole2.nn*aa_quad_deux
-2.*contract(hole2.dn, 0, hole2.psi_auto
.derive_con(hole2.ff), 0)/hole2.psi ;
*/
source_deux = hole2.get_psi4()*( hole2.nn*( aa_quad_deux + 0.3333333333333333
* hole2.trK*hole2.trK*hole2.decouple)
- hole2.trK_point*hole2.decouple )
-2.*contract(contract(hole2.hh, 0, hole2.n_auto
.derive_cov(hole2.ff), 0), 0, hole2.dpsi, 0)/hole2.psi
-2.*contract(hole2.dn, 0, hole2.psi_auto
.derive_con(hole2.ff), 0)/hole2.psi ;
- contract(hdirac2, 0, hole2.n_auto.derive_cov(hole2.ff), 0) ;
Scalar tmp_deux (hole2.mp) ;
tmp_deux = hole2.get_psi4()* contract(hole2.beta_auto, 0, hole2.trK.
derive_cov(hole2.ff), 0)
- contract( hole2.hh, 0, 1, hole2.n_auto.derive_cov(hole2.ff)
.derive_cov(hole2.ff), 0, 1 ) ;
tmp_deux.inc_dzpuis() ; // dzpuis: 3 -> 4
source_deux += tmp_deux ;
cout << "source lapse" << endl << norme(source_un) << endl ;
// Boundary conditions and resolution
// -----------------------------------
Valeur lim_un (hole1.mp.get_mg()-> get_angu()) ;
Valeur lim_deux (hole1.mp.get_mg()-> get_angu()) ;
Scalar n_un_temp (hole1.n_auto) ;
Scalar n_deux_temp (hole2.n_auto) ;
switch (bound_nn) {
case 0 : {
lim_un = hole1.boundary_nn_Dir(lim_nn) ;
lim_deux = hole2.boundary_nn_Dir(lim_nn) ;
n_un_temp = n_un_temp - 1./2. ;
n_deux_temp = n_deux_temp - 1./2. ;
dirichlet_binaire (source_un, source_deux, lim_un, lim_deux,
n_un_temp, n_deux_temp, 0, precision) ;
break ;
}
case 1 : {
lim_un = hole1.boundary_nn_Neu(lim_nn) ;
lim_deux = hole2.boundary_nn_Neu(lim_nn) ;
neumann_binaire (source_un, source_deux, lim_un, lim_deux,
n_un_temp, n_deux_temp, 0, precision) ;
break ;
}
default : {
cout << "Unexpected type of boundary conditions for the lapse!"
<< endl
<< " bound_nn = " << bound_nn << endl ;
abort() ;
break ;
}
} // End of switch
n_un_temp = n_un_temp + 1./2. ;
n_deux_temp = n_deux_temp + 1./2. ;
n_un_temp.raccord(3) ;
n_deux_temp.raccord(3) ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
int nz = hole1.mp.get_mg()->get_nzone() ;
cout << "lapse auto" << endl << norme (n_un_temp) << endl ;
Tbl tdiff_nn = diffrel(n_un_temp.laplacian(), source_un) ;
cout <<
"Relative error in the resolution of the equation for the lapse : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_nn(l) << " " ;
}
cout << endl ;
// Relaxation :
// -------------
n_un_temp = relax*n_un_temp + (1-relax)*lapse_un_old ;
n_deux_temp = relax*n_deux_temp + (1-relax)*lapse_deux_old ;
hole1.n_auto = n_un_temp ;
hole2.n_auto = n_deux_temp ;
hole1.n_comp_import(hole2) ;
hole2.n_comp_import(hole1) ;
}
// Resolution for Psi
// -------------------
void Bin_hor::solve_psi (double precision, double relax, int bound_psi) {
assert ((relax>0) && (relax<=1)) ;
cout << "-----------------------------------------------" << endl ;
cout << "Resolution PSI" << endl ;
Scalar psi_un_old (hole1.psi_auto) ;
Scalar psi_deux_old (hole2.psi_auto) ;
Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
// Source 1
// ---------
Scalar source_un (hole1.mp) ;
/*
// Conformally flat source
source_un.annule_hard() ;
source_un.set_dzpuis(4) ;
source_un += - hole1.psi*hole1.get_psi4()* 0.125* aa_quad_un ;
source_un.std_spectral_base() ;
*/
Scalar tmp_un (hole1.mp) ;
tmp_un = 0.125* hole1.psi_auto * (hole1.tgam).ricci_scal()
- contract(hole1.hh, 0, 1, hole1.psi_auto.derive_cov(hole1.ff)
.derive_cov(hole1.ff), 0, 1 ) ;
tmp_un.inc_dzpuis() ; // dzpuis : 3 -> 4
tmp_un -= contract(hdirac1, 0, hole1.psi_auto
.derive_cov(hole1.ff), 0) ;
source_un = tmp_un - hole1.psi*hole1.get_psi4()* ( 0.125* aa_quad_un
- 8.33333333333333e-2* hole1.trK*hole1.trK*hole1.decouple ) ;
source_un.std_spectral_base() ;
// Source 2
// ---------
Scalar source_deux (hole2.mp) ;
/*
// Conformally flat source
source_deux.annule_hard() ;
source_deux.set_dzpuis(4) ;
source_deux += - hole2.psi*hole2.get_psi4()* 0.125* aa_quad_deux ;
source_deux.std_spectral_base() ;
*/
Scalar tmp_deux (hole2.mp) ;
tmp_deux = 0.125* hole2.psi_auto * (hole2.tgam).ricci_scal()
- contract(hole2.hh, 0, 1, hole2.psi_auto.derive_cov(hole2.ff)
.derive_cov(hole2.ff), 0, 1 ) ;
tmp_deux.inc_dzpuis() ; // dzpuis : 3 -> 4
tmp_deux -= contract(hdirac2, 0, hole2.psi_auto
.derive_cov(hole2.ff), 0) ;
source_deux = tmp_deux - hole2.psi*hole2.get_psi4()* ( 0.125* aa_quad_deux
- 8.33333333333333e-2* hole2.trK*hole2.trK*hole2.decouple ) ;
source_deux.std_spectral_base() ;
cout << "source psi" << endl << norme(source_un) << endl ;
// Boundary conditions and resolution :
// ------------------------------------
Valeur lim_un (hole1.mp.get_mg()-> get_angu()) ;
Valeur lim_deux (hole1.mp.get_mg()-> get_angu()) ;
Scalar psi_un_temp (hole1.psi_auto) ;
Scalar psi_deux_temp (hole2.psi_auto) ;
switch (bound_psi) {
case 0 : {
lim_un = hole1.boundary_psi_app_hor() ;
lim_deux = hole2.boundary_psi_app_hor() ;
neumann_binaire (source_un, source_deux, lim_un, lim_deux,
psi_un_temp, psi_deux_temp, 0, precision) ;
break ;
}
default : {
cout << "Unexpected type of boundary conditions for psi!"
<< endl
<< " bound_psi = " << bound_psi << endl ;
abort() ;
break ;
}
} // End of switch
psi_un_temp = psi_un_temp + 1./2. ;
psi_deux_temp = psi_deux_temp + 1./2. ;
psi_un_temp.raccord(3) ;
psi_deux_temp.raccord(3) ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
int nz = hole1.mp.get_mg()->get_nzone() ;
cout << "psi auto" << endl << norme (psi_un_temp) << endl ;
Tbl tdiff_psi = diffrel(psi_un_temp.laplacian(), source_un) ;
cout <<
"Relative error in the resolution of the equation for psi : "
<< endl ;
for (int l=0; l<nz; l++) {
cout << tdiff_psi(l) << " " ;
}
cout << endl ;
// Relaxation :
// -------------
psi_un_temp = relax*psi_un_temp + (1-relax)*psi_un_old ;
psi_deux_temp = relax*psi_deux_temp + (1-relax)*psi_deux_old ;
hole1.psi_auto = psi_un_temp ;
hole2.psi_auto = psi_deux_temp ;
hole1.psi_comp_import(hole2) ;
hole2.psi_comp_import(hole1) ;
hole1.set_der_0x0() ;
hole2.set_der_0x0() ;
//set_hh_Samaya() ;
}
// Resolution for shift with omega fixed.
// --------------------------------------
void Bin_hor::solve_shift (double precision, double relax, int bound_beta,
double omega_eff) {
cout << "------------------------------------------------" << endl ;
cout << "Resolution shift : Omega = " << omega << endl ;
Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
// Source 1
// ---------
Vector source_un (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
/*
// Conformally flat source
source_un = 2.* contract(hole1.aa_auto, 1, hole1.dn, 0)
- 12.*hole1.nn*contract(hole1.aa_auto, 1, hole1.dpsi, 0)
/hole1.psi;
*/
Vector tmp_vect_un (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
source_un = 2.* contract(hole1.aa_auto, 1, hole1.dn, 0)
- 12.*hole1.nn*contract(hole1.aa_auto, 1, hole1.dpsi, 0)
/hole1.psi;
tmp_vect_un = 2./3.* hole1.trK.derive_con(hole1.tgam)
* hole1.decouple ;
tmp_vect_un.inc_dzpuis() ;
source_un += 2.* hole1.nn * ( tmp_vect_un
- contract(hole1.tgam.connect().get_delta(), 1, 2,
hole1.aa_auto, 0, 1) ) ;
Vector vtmp_un = contract(hole1.hh, 0, 1,
hole1.beta_auto.derive_cov(hole1.ff)
.derive_cov(hole1.ff), 1, 2)
+ 1./3.*contract(hole1.hh, 1, hole1.beta_auto
.divergence(hole1.ff).derive_cov(hole1.ff), 0)
- hdirac1.derive_lie(hole1.beta_auto)
+ hole1.gamt_point.divergence(hole1.ff)*hole1.decouple ;
vtmp_un.inc_dzpuis() ; // dzpuis: 3 -> 4
source_un -= vtmp_un ;
source_un += 2./3.* hole1.beta_auto.divergence(hole1.ff)
* hdirac1 ;
source_un.std_spectral_base() ;
// Source 2
// ---------
Vector source_deux (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
/*
// Conformally flat source
source_deux = 2.* contract(hole2.aa_auto, 1, hole2.dn, 0)
- 12.*hole2.nn*contract(hole2.aa_auto, 1, hole2.dpsi, 0)
/hole2.psi;
*/
Vector tmp_vect_deux (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
source_deux = 2.* contract(hole2.aa_auto, 1, hole2.dn, 0)
- 12.*hole2.nn*contract(hole2.aa_auto, 1, hole2.dpsi, 0)
/hole2.psi;
tmp_vect_deux = 2./3.* hole2.trK.derive_con(hole2.tgam)
* hole2.decouple ;
tmp_vect_deux.inc_dzpuis() ;
source_deux += 2.* hole2.nn * ( tmp_vect_deux
- contract(hole2.tgam.connect().get_delta(), 1, 2,
hole2.aa*hole2.decouple, 0, 1) ) ;
Vector vtmp_deux = contract(hole2.hh, 0, 1,
hole2.beta_auto.derive_cov(hole2.ff)
.derive_cov(hole2.ff), 1, 2)
+ 1./3.*contract(hole2.hh, 1, hole2.beta_auto
.divergence(hole2.ff).derive_cov(hole2.ff), 0)
- hdirac2.derive_lie(hole2.beta_auto)
+ hole2.gamt_point.divergence(hole2.ff)*hole2.decouple ;
vtmp_deux.inc_dzpuis() ; // dzpuis: 3 -> 4
source_deux -= vtmp_deux ;
source_deux += 2./3.* hole2.beta_auto.divergence(hole2.ff)
* hdirac2 ;
source_deux.std_spectral_base() ;
Vector source_1 (source_un) ;
Vector source_2 (source_deux) ;
source_1.change_triad(hole1.mp.get_bvect_cart()) ;
source_2.change_triad(hole2.mp.get_bvect_cart()) ;
cout << "source shift_x" << endl << norme(source_1(1)) << endl ;
cout << "source shift_y" << endl << norme(source_1(2)) << endl ;
cout << "source shift_z" << endl << norme(source_1(3)) << endl ;
// Filter for high frequencies.
for (int i=1 ; i<=3 ; i++) {
source_un.set(i).filtre(4) ;
source_deux.set(i).filtre(4) ;
}
// Boundary conditions
// --------------------
Valeur lim_x_un (hole1.mp.get_mg()-> get_angu()) ;
Valeur lim_y_un (hole1.mp.get_mg()-> get_angu()) ;
Valeur lim_z_un (hole1.mp.get_mg()-> get_angu()) ;
Valeur lim_x_deux (hole2.mp.get_mg()-> get_angu()) ;
Valeur lim_y_deux (hole2.mp.get_mg()-> get_angu()) ;
Valeur lim_z_deux (hole2.mp.get_mg()-> get_angu()) ;
switch (bound_beta) {
case 0 : {
lim_x_un = hole1.boundary_beta_x(omega, omega_eff) ;
lim_y_un = hole1.boundary_beta_y(omega, omega_eff) ;
lim_z_un = hole1.boundary_beta_z() ;
lim_x_deux = hole2.boundary_beta_x(omega, omega_eff) ;
lim_y_deux = hole2.boundary_beta_y(omega, omega_eff) ;
lim_z_deux = hole2.boundary_beta_z() ;
break ;
}
default : {
cout << "Unexpected type of boundary conditions for beta!"
<< endl
<< " bound_beta = " << bound_beta << endl ;
abort() ;
break ;
}
} // End of switch
// We solve :
// -----------
Vector beta_un_old (hole1.beta_auto) ;
Vector beta_deux_old (hole2.beta_auto) ;
Vector beta1 (hole1.beta_auto) ;
Vector beta2 (hole2.beta_auto) ;
poisson_vect_binaire (1./3., source_un, source_deux,
lim_x_un, lim_y_un, lim_z_un,
lim_x_deux, lim_y_deux, lim_z_deux,
beta1, beta2, 0, precision) ;
beta1.change_triad(hole1.mp.get_bvect_cart()) ;
beta2.change_triad(hole2.mp.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++) {
beta1.set(i).raccord(3) ;
beta2.set(i).raccord(3) ;
}
cout << "shift_auto x" << endl << norme(beta1(1)) << endl ;
cout << "shift_auto y" << endl << norme(beta1(2)) << endl ;
cout << "shift_auto z" << endl << norme(beta1(3)) << endl ;
beta1.change_triad(hole1.mp.get_bvect_spher()) ;
beta2.change_triad(hole2.mp.get_bvect_spher()) ;
// Check: has the Poisson equation been correctly solved ?
// -----------------------------------------------------
int nz = hole1.mp.get_mg()->get_nzone() ;
Vector lap_beta = (beta1.derive_con(hole1.ff)).divergence(hole1.ff)
+ 1./3.* beta1.divergence(hole1.ff).derive_con(hole1.ff) ;
source_un.dec_dzpuis() ;
Tbl tdiff_beta_r = diffrel(lap_beta(1), source_un(1)) ;
Tbl tdiff_beta_t = diffrel(lap_beta(2), source_un(2)) ;
Tbl tdiff_beta_p = diffrel(lap_beta(3), source_un(3)) ;
cout <<
"Relative error in the resolution of the equation for beta : "
<< endl ;
cout << "r component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_r(l) << " " ;
}
cout << endl ;
cout << "t component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_t(l) << " " ;
}
cout << endl ;
cout << "p component : " ;
for (int l=0; l<nz; l++) {
cout << tdiff_beta_p(l) << " " ;
}
cout << endl ;
// Relaxation
// -----------
Vector beta1_new (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
Vector beta2_new (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
// Construction of Omega d/dphi
// ----------------------------
Vector omdsdp1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
Scalar yya1 (hole1.mp) ;
yya1 = hole1.mp.ya ;
Scalar xxa1 (hole1.mp) ;
xxa1 = hole1.mp.xa ;
if (fabs(hole1.mp.get_rot_phi()) < 1e-10){
omdsdp1.set(1) = - omega * yya1 ;
omdsdp1.set(2) = omega * xxa1 ;
omdsdp1.set(3).annule_hard() ;
}
else{
omdsdp1.set(1) = omega * yya1 ;
omdsdp1.set(2) = - omega * xxa1 ;
omdsdp1.set(3).annule_hard() ;
}
omdsdp1.set(1).set_spectral_va()
.set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[0])) ;
omdsdp1.set(2).set_spectral_va()
.set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[1])) ;
omdsdp1.set(3).set_spectral_va()
.set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[2])) ;
omdsdp1.annule_domain(nz-1) ;
omdsdp1.change_triad(hole1.mp.get_bvect_spher()) ;
Vector omdsdp2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
Scalar yya2 (hole2.mp) ;
yya2 = hole2.mp.ya ;
Scalar xxa2 (hole2.mp) ;
xxa2 = hole2.mp.xa ;
if (fabs(hole2.mp.get_rot_phi()) < 1e-10){
omdsdp2.set(1) = - omega * yya2 ;
omdsdp2.set(2) = omega * xxa2 ;
omdsdp2.set(3).annule_hard() ;
}
else{
omdsdp2.set(1) = omega * yya2 ;
omdsdp2.set(2) = - omega * xxa2 ;
omdsdp2.set(3).annule_hard() ;
}
omdsdp2.set(1).set_spectral_va()
.set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[0])) ;
omdsdp2.set(2).set_spectral_va()
.set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[1])) ;
omdsdp2.set(3).set_spectral_va()
.set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[2])) ;
omdsdp2.annule_domain(nz-1) ;
omdsdp2.change_triad(hole2.mp.get_bvect_spher()) ;
// New shift
beta1_new = relax*(beta1+hole1.decouple*omdsdp1) + (1-relax)*beta_un_old ;
beta2_new = relax*(beta2+hole2.decouple*omdsdp2) + (1-relax)*beta_deux_old ;
hole1.beta_auto = beta1_new ;
hole2.beta_auto = beta2_new ;
hole1.beta_comp_import(hole2) ;
hole2.beta_comp_import(hole1) ;
// Regularisation of the shifts if necessary
// -----------------------------------------
int nnt = hole1.mp.get_mg()->get_nt(1) ;
int nnp = hole1.mp.get_mg()->get_np(1) ;
int check ;
check = 0 ;
for (int k=0; k<nnp; k++)
for (int j=0; j<nnt; j++){
if (fabs((hole1.n_auto+hole1.n_comp).val_grid_point(1, k, j , 0)) < 1e-8){
check = 1 ;
break ;
}
}
if (check == 1){
double diff_un = hole1.regularisation (hole1.beta_auto,
hole2.beta_auto, omega) ;
double diff_deux = hole2.regularisation (hole2.beta_auto,
hole1.beta_auto, omega) ;
hole1.regul = diff_un ;
hole2.regul = diff_deux ;
}
else {
hole1.regul = 0. ;
hole2.regul = 0. ;
}
}
}
Lorene::Bin_hor::solve_psi
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
Definition binhor_equations.C:301
Lorene::Bin_hor::omega
double omega
Angular velocity.
Definition isol_hor.h:1348
Lorene::Bin_hor::hole1
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Lorene::Bin_hor::hole2
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
Lorene::Bin_hor::solve_shift
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
Definition binhor_equations.C:453
Lorene::Bin_hor::solve_lapse
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ,...
Definition binhor_equations.C:124
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition scalar_raccord.C:62
Lorene::Scalar::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition scalar_manip.C:151
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition scalar_deriv.C:436
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
Lorene::Scalar::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:470
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition scalar_deriv.C:402
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:810
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:75
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
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::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
Lorene::Tensor::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
Lorene::Tensor::derive_lie
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition tensor_calculus.C:481
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
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