init_data.C
/*
* Method of class Hor_isol to compute valid initial data for standard boundary
* conditions
*
* (see file isol_hor.h for documentation).
*
*/
/*
* Copyright (c) 2004 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 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 init_data_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $" ;
/*
* $Id: init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $
* $Log: init_data.C,v $
* Revision 1.30 2014/10/13 08:53:01 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.29 2014/10/06 15:13:10 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.28 2008/08/19 06:42:00 j_novak
* Minor modifications to avoid warnings with gcc 4.3. Most of them concern
* cast-type operations, and constant strings that must be defined as const char*
*
* Revision 1.27 2006/02/22 17:02:04 f_limousin
* Removal of warnings
*
* Revision 1.26 2006/02/22 16:29:55 jl_jaramillo
* corrections on the relaxation and boundary conditions
*
* Revision 1.24 2006/01/18 09:04:27 f_limousin
* Minor modifs (warnings and errors at the compilation with gcc-3.4)
*
* Revision 1.23 2006/01/16 17:13:40 jl_jaramillo
* function for solving the spherical case
*
* Revision 1.22 2005/11/02 16:09:44 jl_jaramillo
* changes in boundary_nn_Dir_lapl
*
* Revision 1.21 2005/10/24 16:44:40 jl_jaramillo
* Cook boundary condition ans minot bound of kss
*
* Revision 1.20 2005/10/21 16:20:55 jl_jaramillo
* Version for the paper JaramL05
*
* Revision 1.19 2005/07/08 13:15:23 f_limousin
* Improvements of boundary_vv_cart(), boundary_nn_lapl().
* Add a fonction to compute the departure of axisymmetry.
*
* Revision 1.18 2005/06/09 08:05:32 f_limousin
* Implement a new function sol_elliptic_boundary() and
* Vector::poisson_boundary(...) which solve the vectorial poisson
* equation (method 6) with an inner boundary condition.
*
* Revision 1.17 2005/05/12 14:48:07 f_limousin
* New boundary condition for the lapse : boundary_nn_lapl().
*
* Revision 1.16 2005/04/08 12:16:52 f_limousin
* Function set_psi(). And dependance in phi.
*
* Revision 1.15 2005/04/03 19:48:22 f_limousin
* Implementation of set_psi(psi_in). And minor changes to avoid warnings.
*
* Revision 1.14 2005/04/02 15:49:21 f_limousin
* New choice (Lichnerowicz) for aaquad. New member data nz.
*
* Revision 1.13 2005/03/31 09:45:31 f_limousin
* New functions compute_ww(...) and aa_kerr_ww().
*
* Revision 1.12 2005/03/24 16:50:28 f_limousin
* Add parameters solve_shift and solve_psi in par_isol.d and in function
* init_dat(...). Implement Isolhor::kerr_perturb().
*
* Revision 1.11 2005/03/22 13:25:36 f_limousin
* Small changes. The angular velocity and A^{ij} are computed
* with a differnet sign.
*
* Revision 1.10 2005/03/09 10:18:08 f_limousin
* Save K_{ij}s^is^j in a file. Add solve_lapse in a file
*
* Revision 1.9 2005/03/06 16:56:13 f_limousin
* The computation of A^{ij} is no more necessary here thanks to the new
* function Isol_hor::aa().
*
* Revision 1.8 2005/03/04 17:04:57 jl_jaramillo
* Addition of boost to the shift after solving the shift equation
*
* Revision 1.7 2005/03/03 10:03:55 f_limousin
* The boundary conditions for the lapse, psi and shift are now
* parameters (in file par_hor.d).
*
* Revision 1.6 2004/12/22 18:15:30 f_limousin
* Many different changes.
*
* Revision 1.5 2004/11/08 14:51:21 f_limousin
* A regularisation for the computation of A^{ij } is done in the
* case lapse equal to zero on the horizon.
*
* Revision 1.1 2004/10/29 12:54:53 jl_jaramillo
* First version
*
* Revision 1.4 2004/10/01 16:47:51 f_limousin
* Case \alpha=0 included
*
* Revision 1.3 2004/09/28 16:10:05 f_limousin
* Many improvements. Now the resolution for the shift is working !
*
* Revision 1.1 2004/09/09 16:41:50 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Isol_hor/init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $
*
*/
// C++ headers
#include "headcpp.h"
// C headers
#include <cstdlib>
#include <cassert>
// Lorene headers
#include "isol_hor.h"
#include "metric.h"
#include "unites.h"
#include "graphique.h"
#include "cmp.h"
#include "tenseur.h"
#include "utilitaires.h"
#include "param.h"
namespace Lorene {
void Isol_hor::init_data(int bound_nn, double lim_nn, int bound_psi,
int bound_beta, int solve_lapse, int solve_psi,
int solve_shift, double precis,
double relax_nn, double relax_psi, double relax_beta, int niter) {
using namespace Unites ;
// Initialisations
// ---------------
double ttime = the_time[jtime] ;
ofstream conv("resconv.d") ;
ofstream kss("kss.d") ;
conv << " # diff_nn diff_psi diff_beta " << endl ;
// Iteration
// ---------
// double relax_nn_fin = relax_nn ;
// double relax_psi_fin = relax_psi ;
// double relax_beta_fin = relax_beta ;
for (int mer=0; mer<niter; mer++) {
//=========================================================
// Boundary conditions and resolution of elliptic equations
//=========================================================
// Resolution of the Poisson equation for the lapse
// ------------------------------------------------
double relax_init = 0.05 ;
double relax_speed = 0.005 ;
double corr = 1 - (1 - relax_init) * exp (- relax_speed * mer) ;
// relax_nn = relax_nn_fin - ( relax_nn_fin - relax_init ) * exp (- relax_speed * mer) ;
// relax_psi = relax_psi_fin - ( relax_psi_fin - relax_init ) * exp (- relax_speed * mer) ;
// relax_beta = relax_beta_fin - ( relax_beta_fin - relax_init ) * exp (- relax_speed * mer) ;
cout << "nn = " << mer << " corr = " << corr << endl ;
cout << " relax_nn = " << relax_nn << endl ;
cout << " relax_psi = " << relax_psi << endl ;
cout << " relax_beta = " << relax_beta << endl ;
Scalar sou_nn (source_nn()) ;
Scalar nn_jp1 (mp) ;
if (solve_lapse == 1) {
Valeur nn_bound (mp.get_mg()-> get_angu()) ;
switch (bound_nn) {
case 0 : {
nn_bound = boundary_nn_Dir(lim_nn) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 1 : {
nn_bound = boundary_nn_Neu_eff(lim_nn) ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
case 2 : {
nn_bound = boundary_nn_Dir_eff(lim_nn) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 3 : {
nn_bound = boundary_nn_Neu_kk(mer) ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
case 4 : {
nn_bound = boundary_nn_Dir_kk() ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 5 : {
nn_bound = boundary_nn_Dir_lapl(mer) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 6 : {
nn_bound = boundary_nn_Neu_Cook() ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for the lapse!"
<< endl
<< " bound_nn = " << bound_nn << endl ;
abort() ;
break ;
}
} // End of switch
// Test:
maxabs(nn_jp1.laplacian() - sou_nn,
"Absolute error in the resolution of the equation for N") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
if (mer==0)
n_evol.update(nn_jp1, jtime, ttime) ;
else
nn_jp1 = relax_nn * nn_jp1 + (1 - relax_nn) * nn() ;
}
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar sou_psi (source_psi()) ;
Scalar psi_jp1 (mp) ;
if (solve_psi == 1) {
Valeur psi_bound (mp.get_mg()-> get_angu()) ;
switch (bound_psi) {
case 0 : {
psi_bound = boundary_psi_app_hor() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 1 : {
psi_bound = boundary_psi_Neu_spat() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 2 : {
psi_bound = boundary_psi_Dir_spat() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
case 3 : {
psi_bound = boundary_psi_Neu_evol() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 4 : {
psi_bound = boundary_psi_Dir_evol() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
case 5 : {
psi_bound = boundary_psi_Dir() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for psi!"
<< endl
<< " bound_psi = " << bound_psi << endl ;
abort() ;
break ;
}
} // End of switch
// Test:
maxabs(psi_jp1.laplacian() - sou_psi,
"Absolute error in the resolution of the equation for Psi") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi() ;
}
// Resolution of the vector Poisson equation for the shift
//---------------------------------------------------------
// Source
Vector beta_jp1(beta()) ;
if (solve_shift == 1) {
Vector source_vector ( source_beta() ) ;
double lambda = 0; //1./3.;
Vector source_reg = - (1./3. - lambda) * beta().divergence(ff)
.derive_con(ff) ;
source_reg.inc_dzpuis() ;
source_vector = source_vector + source_reg ;
// CARTESIAN CASE
// #################################
// Boundary values
Valeur boundary_x (mp.get_mg()-> get_angu()) ;
Valeur boundary_y (mp.get_mg()-> get_angu()) ;
Valeur boundary_z (mp.get_mg()-> get_angu()) ;
switch (bound_beta) {
case 0 : {
boundary_x = boundary_beta_x(omega) ;
boundary_y = boundary_beta_y(omega) ;
boundary_z = boundary_beta_z() ;
break ;
}
case 1 : {
boundary_x = boundary_vv_x(omega) ;
boundary_y = boundary_vv_y(omega) ;
boundary_z = boundary_vv_z(omega) ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for psi!"
<< endl
<< " bound_psi = " << bound_psi << endl ;
abort() ;
break ;
}
} // End of switch
if (boost_x != 0.)
boundary_x -= beta_boost_x() ;
if (boost_z != 0.)
boundary_z -= beta_boost_z() ;
// Resolution
//-----------
double precision = 1e-8 ;
poisson_vect_boundary(lambda, source_vector, beta_jp1, boundary_x,
boundary_y, boundary_z, 0, precision, 20) ;
/*
// SPHERICAL CASE
// #################################
// Boundary values
Valeur boundary_r (mp.get_mg()-> get_angu()) ;
Valeur boundary_t (mp.get_mg()-> get_angu()) ;
Valeur boundary_p (mp.get_mg()-> get_angu()) ;
switch (bound_beta) {
case 0 : {
boundary_r = boundary_beta_r() ;
boundary_t = boundary_beta_theta() ;
boundary_p = boundary_beta_phi(omega) ;
break ;
}
case 1 : {
boundary_r = boundary_vv_x(omega) ;
boundary_t = boundary_vv_y(omega) ;
boundary_p = boundary_vv_z(omega) ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for psi!"
<< endl
<< " bound_psi = " << bound_psi << endl ;
abort() ;
break ;
}
} // End of switch
// Resolution
//-----------
beta_jp1 = source_vector.poisson_dirichlet(lambda, boundary_r,
boundary_t, boundary_p, 0) ;
des_meridian(beta_jp1(1), 1.0000001, 10., "beta_r", 0) ;
des_meridian(beta_jp1(2), 1.0000001, 10., "beta_t", 1) ;
des_meridian(beta_jp1(3), 1.0000001, 10., "beta_p", 2) ;
arrete() ;
// #########################
// End of spherical case
// #########################
*/
// Test
source_vector.dec_dzpuis() ;
maxabs(beta_jp1.derive_con(ff).divergence(ff)
+ lambda * beta_jp1.divergence(ff)
.derive_con(ff) - source_vector,
"Absolute error in the resolution of the equation for beta") ;
cout << endl ;
// Boost
// -----
Vector boost_vect(mp, CON, mp.get_bvect_cart()) ;
if (boost_x != 0.) {
boost_vect.set(1) = boost_x ;
boost_vect.set(2) = 0. ;
boost_vect.set(3) = 0. ;
boost_vect.std_spectral_base() ;
boost_vect.change_triad(mp.get_bvect_spher()) ;
beta_jp1 = beta_jp1 + boost_vect ;
}
if (boost_z != 0.) {
boost_vect.set(1) = boost_z ;
boost_vect.set(2) = 0. ;
boost_vect.set(3) = 0. ;
boost_vect.std_spectral_base() ;
boost_vect.change_triad(mp.get_bvect_spher()) ;
beta_jp1 = beta_jp1 + boost_vect ;
}
// Relaxation (relax=1 -> new ; relax=0 -> old )
beta_jp1 = relax_beta * beta_jp1 + (1 - relax_beta) * beta() ;
}
//===========================================
// Convergence control
//===========================================
double diff_nn, diff_psi, diff_beta ;
diff_nn = 1.e-16 ;
diff_psi = 1.e-16 ;
diff_beta = 1.e-16 ;
if (solve_lapse == 1)
diff_nn = max( diffrel(nn(), nn_jp1) ) ;
if (solve_psi == 1)
diff_psi = max( diffrel(psi(), psi_jp1) ) ;
if (solve_shift == 1)
diff_beta = max( maxabs(beta_jp1 - beta()) ) ;
cout << "step = " << mer << " : diff_psi = " << diff_psi
<< " diff_nn = " << diff_nn
<< " diff_beta = " << diff_beta << endl ;
cout << "----------------------------------------------" << endl ;
if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
break ;
if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
<< " " << log10(diff_beta) << endl ; } ;
//=============================================
// Updates for next step
//=============================================
if (solve_psi == 1)
set_psi(psi_jp1) ;
if (solve_lapse == 1)
n_evol.update(nn_jp1, jtime, ttime) ;
if (solve_shift == 1)
beta_evol.update(beta_jp1, jtime, ttime) ;
if (solve_shift == 1)
update_aa() ;
// Saving ok K_{ij}s^is^j
// -----------------------
Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
gam().radial_vect(), 0, 1)) ;
double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
Scalar aaLss (pow(psi(), 6) * kkss) ;
double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
int nnp = mp.get_mg()->get_np(1) ;
int nnt = mp.get_mg()->get_nt(1) ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
if (kkss.val_grid_point(1, k, j, 0) > max_kss)
max_kss = kkss.val_grid_point(1, k, j, 0) ;
if (kkss.val_grid_point(1, k, j, 0) < min_kss)
min_kss = kkss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
}
kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
<< " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
}
conv.close() ;
kss.close() ;
}
/*
void Isol_hor::init_data_loop(int bound_nn, double lim_nn, int bound_psi,
int bound_beta, int solve_lapse, int solve_psi,
int solve_shift, double precis, double precis_loop,
double relax_nn, double relax_psi, double relax_beta, double relax_loop, int niter) {
using namespace Unites ;
// Initialisations
// ---------------
double ttime = the_time[jtime] ;
ofstream conv("resconv.d") ;
ofstream kss("kss.d") ;
conv << " # diff_nn diff_psi diff_beta " << endl ;
// Iteration
// ---------
for (int mer=0; mer<niter; mer++) {
//=========================================================
// Boundary conditions and resolution of elliptic equations
//=========================================================
// Resolution of the Poisson equation for the lapse
// ------------------------------------------------
Scalar sou_nn (source_nn()) ;
Scalar nn_jp1 (mp) ;
if (solve_lapse == 1) {
Valeur nn_bound (mp.get_mg()-> get_angu()) ;
switch (bound_nn) {
case 0 : {
nn_bound = boundary_nn_Dir(lim_nn) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 1 : {
nn_bound = boundary_nn_Neu_eff(lim_nn) ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
case 2 : {
nn_bound = boundary_nn_Dir_eff(lim_nn) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 3 : {
nn_bound = boundary_nn_Neu_kk() ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
case 4 : {
nn_bound = boundary_nn_Dir_kk() ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 5 : {
nn_bound = boundary_nn_Dir_lapl(mer) ;
nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
break ;
}
case 6 : {
nn_bound = boundary_nn_Neu_Cook() ;
nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for the lapse!"
<< endl
<< " bound_nn = " << bound_nn << endl ;
abort() ;
break ;
}
} // End of switch
// Test:
maxabs(nn_jp1.laplacian() - sou_nn,
"Absolute error in the resolution of the equation for N") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
if (mer==0)
n_evol.update(nn_jp1, jtime, ttime) ;
else
nn_jp1 = relax_nn * nn_jp1 + (1 - relax_nn) * nn() ;
}
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar sou_psi (source_psi()) ;
Scalar psi_jp1 (mp) ;
if (solve_psi == 1) {
Valeur psi_bound (mp.get_mg()-> get_angu()) ;
switch (bound_psi) {
case 0 : {
psi_bound = boundary_psi_app_hor() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 1 : {
psi_bound = boundary_psi_Neu_spat() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 2 : {
psi_bound = boundary_psi_Dir_spat() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
case 3 : {
psi_bound = boundary_psi_Neu_evol() ;
psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
break ;
}
case 4 : {
psi_bound = boundary_psi_Dir_evol() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
case 5 : {
psi_bound = boundary_psi_Dir() ;
psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for psi!"
<< endl
<< " bound_psi = " << bound_psi << endl ;
abort() ;
break ;
}
} // End of switch
// Test:
maxabs(psi_jp1.laplacian() - sou_psi,
"Absolute error in the resolution of the equation for Psi") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi() ;
}
// Resolution of the vector Poisson equation for the shift
//---------------------------------------------------------
// Source
Vector beta_j(beta()) ;
if (solve_shift == 1) {
double lambda = 1./3.;
Vector beta_jp1 (beta()) ;
double thresh_loop = 1;
int n_loop = 0 ;
while( thresh_loop > precis_loop ){
Vector source_vector ( source_beta() ) ;
Vector source_reg = - (1./3. - lambda) * beta().divergence(ff)
.derive_con(ff) ;
source_reg.inc_dzpuis() ;
source_vector = source_vector + source_reg ;
// Boundary values
// ===============
Valeur boundary_x (mp.get_mg()-> get_angu()) ;
Valeur boundary_y (mp.get_mg()-> get_angu()) ;
Valeur boundary_z (mp.get_mg()-> get_angu()) ;
switch (bound_beta) {
case 0 : {
boundary_x = boundary_beta_x(omega) ;
boundary_y = boundary_beta_y(omega) ;
boundary_z = boundary_beta_z() ;
break ;
}
case 1 : {
boundary_x = boundary_vv_x(omega) ;
boundary_y = boundary_vv_y(omega) ;
boundary_z = boundary_vv_z(omega) ;
break ;
}
default : {
cout <<"Unexpected type of boundary conditions for beta!"
<< endl
<< " bound_beta = " << bound_beta << endl ;
abort() ;
break ;
}
} // End of switch
if (boost_x != 0.)
boundary_x -= beta_boost_x() ;
if (boost_z != 0.)
boundary_z -= beta_boost_z() ;
// Resolution
//-----------
double precision = 1e-8 ;
poisson_vect_boundary(lambda, source_vector, beta_jp1, boundary_x,
boundary_y, boundary_z, 0, precision, 20) ;
// Test
source_vector.dec_dzpuis() ;
maxabs(beta_jp1.derive_con(ff).divergence(ff)
+ lambda * beta_jp1.divergence(ff)
.derive_con(ff) - source_vector,
"Absolute error in the resolution of the equation for beta") ;
cout << endl ;
// Relaxation_loop (relax=1 -> new ; relax=0 -> old )
beta_jp1 = relax_loop * beta_jp1 + (1 - relax_loop) * beta() ;
// Convergence loop
//=================
double diff_beta_loop ;
diff_beta_loop = 1.e-16 ;
if (solve_shift == 1)
diff_beta_loop = max( maxabs(beta_jp1 - beta()) ) ;
cout << "step_loop = " << n_loop <<
<< " diff_beta_loop = " << diff_beta_loop << endl ;
cout << "----------------------------------------------" << endl ;
thresh_loop = diff_beta_loop ;
//Update loop
//===========
beta_evol.update(beta_jp1, jtime, ttime) ;
update_aa() ;
n_loop += 1 ;
// End internal loop
}
// Test for resolution of beta at this setp mer is already done in the internal loop
// Relaxation beta (relax=1 -> new ; relax=0 -> old )
beta_jp1 = relax_beta * beta_jp1 + (1 - relax_loop) * beta_j ;
}
//===========================================
// Convergence control
//===========================================
double diff_nn, diff_psi, diff_beta ;
diff_nn = 1.e-16 ;
diff_psi = 1.e-16 ;
diff_beta = 1.e-16 ;
if (solve_lapse == 1)
diff_nn = max( diffrel(nn(), nn_jp1) ) ;
if (solve_psi == 1)
diff_psi = max( diffrel(psi(), psi_jp1) ) ;
if (solve_shift == 1)
diff_beta = max( maxabs(beta_jp1 - beta_j) ) ;
cout << "step = " << mer << " : diff_psi = " << diff_psi
<< " diff_nn = " << diff_nn
<< " diff_beta = " << diff_beta << endl ;
cout << "----------------------------------------------" << endl ;
if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
break ;
if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
<< " " << log10(diff_beta) << endl ; } ;
//=============================================
// Updates for next step
//=============================================
if (solve_psi == 1)
set_psi(psi_jp1) ;
if (solve_lapse == 1)
n_evol.update(nn_jp1, jtime, ttime) ;
if (solve_shift == 1)
beta_evol.update(beta_jp1, jtime, ttime) ;
if (solve_shift == 1)
update_aa() ;
// Saving ok K_{ij}s^is^j
// -----------------------
Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
gam().radial_vect(), 0, 1)) ;
double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
Scalar aaLss (pow(psi(), 6) * kkss) ;
double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
int nnp = mp.get_mg()->get_np(1) ;
int nnt = mp.get_mg()->get_nt(1) ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
if (kkss.val_grid_point(1, k, j, 0) > max_kss)
max_kss = kkss.val_grid_point(1, k, j, 0) ;
if (kkss.val_grid_point(1, k, j, 0) < min_kss)
min_kss = kkss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
}
kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
<< " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
}
conv.close() ;
kss.close() ;
}
*/
void Isol_hor::init_data_spher(int bound_nn, double lim_nn, int bound_psi,
int bound_beta, int solve_lapse, int solve_psi,
int solve_shift, double precis,
double relax, int niter) {
using namespace Unites ;
// Initialisations
// ---------------
double ttime = the_time[jtime] ;
ofstream conv("resconv.d") ;
ofstream kss("kss.d") ;
conv << " # diff_nn diff_psi diff_beta " << endl ;
// Iteration
// ---------
for (int mer=0; mer<niter; mer++) {
// des_meridian(psi(), 1, 10., "psi", 0) ;
// des_meridian(b_tilde(), 1, 10., "b_tilde", 1) ;
// des_meridian(nn(), 1, 10., "nn", 2) ;
// arrete() ;
//========
// Sources
//========
// Useful functions
// ----------------
Vector tem_vect (beta() ) ;
Scalar dif_b = b_tilde() - tem_vect.set(1) ;
// cout << "dif_b = " << dif_b << endl ;
// arrete() ;
Scalar dbdr ( b_tilde().dsdr() ) ;
Scalar bsr (b_tilde()) ;
bsr.div_r() ;
bsr.inc_dzpuis(2) ;
Scalar bsr2 ( bsr) ;
bsr2.div_r() ;
bsr2.inc_dzpuis(2) ;
Scalar psisr (psi()) ;
psisr.div_r() ;
psisr.inc_dzpuis(2) ;
// Source Psi
// ----------
Scalar source_psi_spher(mp) ;
source_psi_spher = -1./12. * psi4()*psi()/(nn() * nn()) * (dbdr - bsr) * (dbdr - bsr) ;
// Source N
//---------
Scalar source_nn_spher(mp) ;
source_nn_spher = 2./3. * psi4() /nn() * (dbdr - bsr) * (dbdr - bsr)
- 2 * ln_psi().dsdr() * nn().dsdr() ;
// Source b_tilde
//---------------
Scalar source_btilde_spher(mp) ;
Scalar tmp ( -1./3. * (dbdr + 2 * bsr).dsdr() ) ;
tmp.std_spectral_base() ;
tmp.inc_dzpuis() ;
source_btilde_spher = tmp + 2 * bsr2
+ 4./3. * (dbdr - bsr) * ( nn().dsdr()/nn() - 6 * psi().dsdr()/psi() ) ;
Scalar source_btilde_trun(mp) ;
source_btilde_trun = tmp +
4./3. * (dbdr - bsr) * ( nn().dsdr()/nn() - 6 * psi().dsdr()/psi() ) ;
// Scalar diff_dbeta ( (dbdr + 2 * bsr).dsdr() - beta().divergence(ff).derive_con(ff)(1) ) ;
// Parallel calculation
//---------------------
Scalar sourcepsi (source_psi()) ;
Scalar sourcenn (source_nn()) ;
Vector sourcebeta (source_beta()) ;
Vector source_reg = 1./3. * beta().divergence(ff).derive_con(ff) ;
source_reg.inc_dzpuis() ;
sourcebeta -= source_reg ;
Scalar source_btilde (sourcebeta(1) ) ;
// Scalar diff_div = source_reg(1) + tmp ; ;
Scalar mag_sou_psi ( source_psi_spher ) ;
mag_sou_psi.dec_dzpuis(4) ;
Scalar mag_sou_nn ( source_nn_spher ) ;
mag_sou_nn.dec_dzpuis(4) ;
Scalar mag_sou_btilde ( source_btilde_trun ) ;
mag_sou_btilde.dec_dzpuis(4) ;
Scalar diff_sou_psi ( source_psi_spher - sourcepsi) ;
diff_sou_psi.dec_dzpuis(4) ;
Scalar diff_sou_nn ( source_nn_spher - sourcenn) ;
diff_sou_nn.dec_dzpuis(4) ;
Scalar diff_sou_btilde ( source_btilde_trun - source_btilde) ;
diff_sou_btilde.dec_dzpuis(4) ;
/*
cout << "dzpuis mag_btilde =" << mag_sou_btilde.get_dzpuis()<<endl ;
des_meridian(diff_sou_psi, 1, 10., "diff_psi", 0) ;
des_meridian(diff_sou_nn, 1, 10., "diff_nn", 1) ;
des_meridian(diff_sou_btilde, 1, 10., "diff_btilde", 2) ;
des_meridian(mag_sou_psi, 1, 10., "mag_psi", 3) ;
des_meridian(mag_sou_nn, 1, 10., "mag_nn", 4) ;
des_meridian(mag_sou_btilde, 1, 10., "mag_btilde", 5) ;
// des_meridian(diff_dbeta, 1, 10., "diff_dbeta", 6) ;
arrete() ;
*/
//====================
// Boundary conditions
//====================
// To avoid warnings;
bound_nn = 1 ; lim_nn = 1. ; bound_psi = 1 ; bound_beta = 1 ;
double kappa_0 = 0.2 - 1. ;
Scalar kappa (mp) ;
kappa = kappa_0 ;
kappa.std_spectral_base() ;
kappa.inc_dzpuis(2) ;
int nnp = mp.get_mg()->get_np(1) ;
int nnt = mp.get_mg()->get_nt(1) ;
Valeur psi_bound (mp.get_mg()-> get_angu()) ;
Valeur nn_bound (mp.get_mg()-> get_angu()) ;
Valeur btilde_bound (mp.get_mg()-> get_angu()) ;
psi_bound = 1. ; // Juste pour affecter dans espace des configs ;
nn_bound = 1. ; // Juste pour affecter dans espace des configs ;
btilde_bound = 1. ; // Juste pour affecter dans espace des configs ;
Scalar tmp_psi = -1./4. * (2 * psisr +
2./3. * psi4()/(psi() * nn()) * (dbdr - bsr) ) ;
Scalar tmp_nn = kappa ; //+ 2./3. * psi() * psi() * (dbdr - bsr) ;
Scalar tmp_btilde = nn() / (psi() * psi()) ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
psi_bound.set(0, k, j, 0) = tmp_psi.val_grid_point(1, k, j, 0) ; // BC Psi
nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ; // BC N
btilde_bound.set(0, k, j, 0) = tmp_btilde.val_grid_point(1, k, j, 0) ; // BC b_tilde
}
psi_bound.std_base_scal() ;
nn_bound.std_base_scal() ;
btilde_bound.std_base_scal() ;
//=================================
// Resolution of elliptic equations
//=================================
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar psi_jp1 (mp) ;
if (solve_psi == 1) {
psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
// Test:
maxabs(psi_jp1.laplacian() - source_psi_spher,
"Absolute error in the resolution of the equation for Psi") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
psi_jp1 = relax * psi_jp1 + (1 - relax) * psi() ;
}
// Resolution of the Poisson equation for the lapse
// ------------------------------------------------
Scalar nn_jp1 (mp) ;
if (solve_lapse == 1) {
nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
// Test:
maxabs(nn_jp1.laplacian() - source_nn_spher,
"Absolute error in the resolution of the equation for N") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
if (mer==0)
n_evol.update(nn_jp1, jtime, ttime) ;
else
nn_jp1 = relax * nn_jp1 + (1 - relax) * nn() ;
}
// Resolution of the Poisson equation for b_tilde
// ----------------------------------------------
Scalar btilde_jp1 (mp) ;
if (solve_shift == 1) {
btilde_jp1 = source_btilde_spher.poisson_dirichlet(btilde_bound, 0) ;
// Test:
maxabs(btilde_jp1.laplacian() - source_btilde_spher,
"Absolute error in the resolution of the equation for btilde") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
btilde_jp1 = relax * btilde_jp1 + (1 - relax) * b_tilde() ;
}
//===========================================
// Convergence control
//===========================================
double diff_nn, diff_psi, diff_btilde ;
diff_nn = 1.e-16 ;
diff_psi = 1.e-16 ;
diff_btilde = 1.e-16 ;
if (solve_lapse == 1)
diff_nn = max( diffrel(nn(), nn_jp1) ) ;
if (solve_psi == 1)
diff_psi = max( diffrel(psi(), psi_jp1) ) ;
if (solve_shift == 1)
diff_btilde = max( diffrel(btilde_jp1, b_tilde()) ) ;
cout << "step = " << mer << " : diff_psi = " << diff_psi
<< " diff_nn = " << diff_nn
<< " diff_btilde = " << diff_btilde << endl ;
cout << "----------------------------------------------" << endl ;
if ((diff_psi<precis) && (diff_nn<precis) && (diff_btilde<precis))
break ;
if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
<< " " << log10(diff_btilde) << endl ; } ;
//=============================================
// Updates for next step
//=============================================
if (solve_psi == 1)
set_psi(psi_jp1) ;
if (solve_lapse == 1)
n_evol.update(nn_jp1, jtime, ttime) ;
if (solve_shift == 1)
{
Vector beta_jp1 (btilde_jp1 * tgam().radial_vect()) ;
cout << tgam().radial_vect() << endl ;
beta_evol.update(beta_jp1, jtime, ttime) ;
}
if (solve_shift == 1 || solve_lapse == 1)
{
update_aa() ;
}
// Saving ok K_{ij}s^is^j
// -----------------------
Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
gam().radial_vect(), 0, 1)) ;
double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
Scalar aaLss (pow(psi(), 6) * kkss) ;
double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
if (kkss.val_grid_point(1, k, j, 0) > max_kss)
max_kss = kkss.val_grid_point(1, k, j, 0) ;
if (kkss.val_grid_point(1, k, j, 0) < min_kss)
min_kss = kkss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
}
kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
<< " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
}
conv.close() ;
kss.close() ;
}
void Isol_hor::init_data_alt(int, double, int,
int, int solve_lapse, int solve_psi,
int solve_shift, double precis,
double relax, int niter) {
using namespace Unites ;
// Initialisations
// ---------------
double ttime = the_time[jtime] ;
ofstream conv("resconv.d") ;
ofstream kss("kss.d") ;
conv << " # diff_nn diff_psi diff_beta " << endl ;
Scalar psi_j (psi()) ;
Scalar nn_j (nn()) ;
Scalar btilde_j (b_tilde()) ;
Scalar diffb ( btilde_j - b_tilde() ) ;
Scalar theta_j (beta().divergence(ff)) ;
theta_j.dec_dzpuis(2) ;
Scalar chi_j (b_tilde()) ;
chi_j.mult_r() ;
// Iteration
// ---------
for (int mer=0; mer<niter; mer++) {
des_meridian(psi_j, 1, 10., "psi", 0) ;
des_meridian(nn_j, 1, 10., "nn", 1) ;
des_meridian(theta_j, 1, 10., "Theta", 2) ;
des_meridian(chi_j, 1, 10., "chi", 3) ;
arrete() ;
//========
// Sources
//========
// Useful functions
// ----------------
Scalar psisr (psi_j) ;
psisr.div_r() ;
psisr.inc_dzpuis(2) ;
Scalar dchidr ( chi_j.dsdr() ) ;
Scalar chisr (chi_j) ;
chisr.div_r() ;
chisr.inc_dzpuis(2) ;
Scalar rdthetadr (theta_j.dsdr() ) ;
rdthetadr.mult_r() ;
rdthetadr.inc_dzpuis(2) ;
Scalar theta_dz4 (theta_j) ;
theta_dz4.inc_dzpuis(4) ;
Scalar dbmb (dchidr - 2 * chisr) ;
dbmb.div_r() ;
// Source Psi
// ----------
Scalar source_psi_spher(mp) ;
source_psi_spher = -1./12. * psi_j*psi_j*psi_j*psi_j*psi_j/(nn_j * nn_j)
* dbmb *dbmb ;
// Source N
//---------
Scalar source_nn_spher(mp) ;
source_nn_spher = 2./3. * psi_j*psi_j*psi_j*psi_j/nn_j * dbmb *dbmb
- 2 * psi_j.dsdr()/psi_j * nn_j.dsdr() ;
//====================
// Boundary conditions
//====================
double kappa_0 = 0.2 - 1. ;
Scalar kappa (mp) ;
kappa = kappa_0 ;
kappa.std_spectral_base() ;
kappa.inc_dzpuis(2) ;
int nnp = mp.get_mg()->get_np(1) ;
int nnt = mp.get_mg()->get_nt(1) ;
Valeur psi_bound (mp.get_mg()-> get_angu()) ;
Valeur nn_bound (mp.get_mg()-> get_angu()) ;
psi_bound = 1. ; // Juste pour affecter dans espace des configs ;
nn_bound = 1. ; // Juste pour affecter dans espace des configs ;
//psi
Scalar tmp_psi = -1./4. * (2 * psisr +
2./3. * psi_j*psi_j*psi_j/ nn_j * dbmb ) ;
// tmp_psi = 2./3. * psi_j*psi_j*psi_j/ nn_j * (dchidr - 2 * chisr) ;
// tmp_psi.div_r() ;
// tmp_psi = -1./4. * (2 * psisr + tmp_psi) ;
//nn
Scalar tmp_nn = kappa ; //+ 2./3. * psi_j*psi_j * dbmb ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
psi_bound.set(0, k, j, 0) = tmp_psi.val_grid_point(1, k, j, 0) ; // BC Psi
nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ; // BC N
}
psi_bound.std_base_scal() ;
nn_bound.std_base_scal() ;
//=================================
// Resolution of elliptic equations
//=================================
// Resolution of the Poisson equation for Psi
// ------------------------------------------
Scalar psi_jp1 (mp) ;
if (solve_psi == 1) {
psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
// Test:
maxabs(psi_jp1.laplacian() - source_psi_spher,
"Absolute error in the resolution of the equation for Psi") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
psi_jp1 = relax * psi_jp1 + (1 - relax) * psi_j ;
}
// Resolution of the Poisson equation for the lapse
// ------------------------------------------------
Scalar nn_jp1 (mp) ;
if (solve_lapse == 1) {
nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
// Test:
maxabs(nn_jp1.laplacian() - source_nn_spher,
"Absolute error in the resolution of the equation for N") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
if (mer==0)
n_evol.update(nn_jp1, jtime, ttime) ;
else
nn_jp1 = relax * nn_jp1 + (1 - relax) * nn_j ;
}
// Resolution for chi and Theta
// ----------------------------
Scalar theta_jp1 (mp) ;
Scalar chi_jp1 (mp) ;
if (solve_shift == 1) {
// Initialisations loop on theta/chi
Scalar theta_i(theta_j) ;
Scalar chi_i(chi_j) ;
// Iteration in theta/chi
for (int i=0 ; i<niter ; i++) {
des_meridian(theta_i, 1, 10., "Theta", 2) ;
des_meridian(chi_i, 1, 10., "chi", 3) ;
arrete() ;
//Sources
// Source_theta
//-------------
Scalar source_theta_spher(mp) ;
source_theta_spher = (dbmb * (nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j)).dsdr() ;
source_theta_spher.dec_dzpuis() ;
// Source chi
//-----------
Scalar source_chi_spher(mp) ;
source_chi_spher = 4./3. * (dchidr - 2 * chisr) * ( nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j )
- 1./3. * rdthetadr + 2 * theta_dz4 ;
//Boundaries
Valeur theta_bound (mp.get_mg()-> get_angu()) ;
Valeur chi_bound (mp.get_mg()-> get_angu()) ;
theta_bound = 1. ; // Juste pour affecter dans espace des configs ;
chi_bound = 1. ; // Juste pour affecter dans espace des configs ;
//theta
Scalar tmp_theta = dchidr ;
tmp_theta.dec_dzpuis(2) ;
tmp_theta += nn_j/(psi_j*psi_j) ;
tmp_theta.div_r() ;
//chi
Scalar tmp_chi = nn_j/(psi_j*psi_j) ;
tmp_chi.mult_r() ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
theta_bound.set(0, k, j, 0) = tmp_theta.val_grid_point(1, k, j, 0) ; // BC Theta
chi_bound.set(0, k, j, 0) = tmp_chi.val_grid_point(1, k, j, 0) ; // BC chi
}
theta_bound.std_base_scal() ;
chi_bound.std_base_scal() ;
//Resolution equations
Scalar theta_ip1(mp) ;
Scalar chi_ip1(mp) ;
theta_ip1 = source_theta_spher.poisson_dirichlet(theta_bound, 0) ;
chi_ip1 = source_chi_spher.poisson_dirichlet(chi_bound, 0) ;
// Test:
maxabs(theta_ip1.laplacian() - source_theta_spher,
"Absolute error in the resolution of the equation for Theta") ;
maxabs(chi_ip1.laplacian() - source_chi_spher,
"Absolute error in the resolution of the equation for chi") ;
// Relaxation (relax=1 -> new ; relax=0 -> old )
theta_ip1 = relax * theta_ip1 + (1 - relax) * theta_i ;
chi_ip1 = relax * chi_ip1 + (1 - relax) * chi_i ;
// Convergence control of loop in theta/chi
double diff_theta_int, diff_chi_int, int_precis ;
diff_theta_int = 1.e-16 ;
diff_chi_int = 1.e-16 ;
int_precis = 1.e-3 ;
diff_theta_int = max( diffrel(theta_ip1, theta_i) ) ;
diff_chi_int = max( diffrel(chi_ip1, chi_i) ) ;
cout << "internal step = " << i
<< " diff_theta_int = " << diff_theta_int
<< " diff_chi_int = " << diff_chi_int << endl ;
cout << "----------------------------------------------" << endl ;
if ((diff_theta_int<int_precis) && (diff_chi_int<int_precis))
{
theta_jp1 = theta_ip1 ;
chi_jp1 = chi_ip1 ;
break ;
}
// Updates of internal loop in theta/chi
theta_i = theta_ip1 ;
chi_i = chi_ip1 ;
}
}
//===========================================
// Convergence control
//===========================================
double diff_nn, diff_psi, diff_theta, diff_chi ;
diff_nn = 1.e-16 ;
diff_psi = 1.e-16 ;
diff_theta = 1.e-16 ;
diff_chi = 1.e-16 ;
if (solve_lapse == 1)
diff_nn = max( diffrel(nn_j, nn_jp1) ) ;
if (solve_psi == 1)
diff_psi = max( diffrel(psi_j, psi_jp1) ) ;
if (solve_shift == 1)
diff_theta = max( diffrel(theta_jp1, theta_j) ) ;
if (solve_shift == 1)
diff_chi = max( diffrel(chi_jp1, chi_j) ) ;
cout << "step = " << mer << " : diff_psi = " << diff_psi
<< " diff_nn = " << diff_nn
<< " diff_theta = " << diff_theta
<< " diff_chi = " << diff_chi << endl ;
cout << "----------------------------------------------" << endl ;
if ((diff_psi<precis) && (diff_nn<precis) && (diff_theta<precis) && (diff_chi<precis))
break ;
if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
<< " " << log10(diff_theta)
<< " " << log10(diff_chi) << endl ; } ;
//=============================================
// Updates for next step
//=============================================
if (solve_psi == 1)
set_psi(psi_jp1) ;
psi_j = psi_jp1 ;
if (solve_lapse == 1)
n_evol.update(nn_jp1, jtime, ttime) ;
nn_j = nn_jp1 ;
if (solve_shift == 1)
{
theta_j = theta_jp1 ;
chi_j = chi_jp1 ;
chi_jp1.mult_r() ;
Vector beta_jp1 (chi_jp1 * tgam().radial_vect()) ;
beta_evol.update(beta_jp1, jtime, ttime) ;
}
if (solve_shift == 1 || solve_lapse == 1)
{
update_aa() ;
}
// Saving ok K_{ij}s^is^j
// -----------------------
Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
gam().radial_vect(), 0, 1)) ;
double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
Scalar aaLss (pow(psi(), 6) * kkss) ;
double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
for (int k=0 ; k<nnp ; k++)
for (int j=0 ; j<nnt ; j++){
if (kkss.val_grid_point(1, k, j, 0) > max_kss)
max_kss = kkss.val_grid_point(1, k, j, 0) ;
if (kkss.val_grid_point(1, k, j, 0) < min_kss)
min_kss = kkss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
}
kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
<< " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
}
conv.close() ;
kss.close() ;
}
}
Lorene::Isol_hor::tgam
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Definition isol_hor.h:514
Lorene::Isol_hor::boundary_vv_y
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
Definition bound_hor.C:1521
Lorene::Isol_hor::boundary_psi_Dir_evol
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution).
Definition bound_hor.C:174
Lorene::Isol_hor::boundary_psi_Dir_spat
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial).
Definition bound_hor.C:231
Lorene::Isol_hor::source_beta
const Vector source_beta() const
Source for .
Definition sources_hor.C:190
Lorene::Isol_hor::omega
double omega
Angular velocity in LORENE's units.
Definition isol_hor.h:269
Lorene::Isol_hor::boundary_nn_Neu_eff
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif) .
Definition bound_hor.C:622
Lorene::Isol_hor::boundary_vv_x
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
Definition bound_hor.C:1493
Lorene::Isol_hor::boundary_beta_y
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
Definition bound_hor.C:1071
Lorene::Isol_hor::boundary_psi_app_hor
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial).
Definition bound_hor.C:297
Lorene::Isol_hor::boundary_beta_x
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Definition bound_hor.C:1021
Lorene::Isol_hor::met_gamt
Metric met_gamt
3 metric tilde
Definition isol_hor.h:326
Lorene::Isol_hor::boundary_nn_Neu_Cook
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
Definition bound_hor.C:489
Lorene::Isol_hor::b_tilde
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition phys_param.C:136
Lorene::Isol_hor::source_psi
const Scalar source_psi() const
Source for .
Definition sources_hor.C:108
Lorene::Isol_hor::beta_boost_x
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
Definition bound_hor.C:1144
Lorene::Isol_hor::boundary_vv_z
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
Definition bound_hor.C:1547
Lorene::Isol_hor::set_psi
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition isol_hor.C:515
Lorene::Isol_hor::boundary_nn_Dir_lapl
const Valeur boundary_nn_Dir_lapl(int mer=1) const
Dirichlet boundary condition for N fixing the divergence of the connection form .
Definition bound_hor.C:693
Lorene::Isol_hor::boost_z
double boost_z
Boost velocity in z-direction.
Definition isol_hor.h:275
Lorene::Isol_hor::beta_boost_z
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
Definition bound_hor.C:1155
Lorene::Isol_hor::update_aa
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition isol_hor.C:842
Lorene::Isol_hor::boundary_beta_z
const Valeur boundary_beta_z() const
Component z of boundary value of .
Definition bound_hor.C:1121
Lorene::Isol_hor::boundary_nn_Dir_kk
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
Definition bound_hor.C:371
Lorene::Isol_hor::mp
Map_af & mp
Affine mapping.
Definition isol_hor.h:260
Lorene::Isol_hor::boundary_psi_Neu_spat
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial).
Definition bound_hor.C:267
Lorene::Isol_hor::boundary_nn_Dir_eff
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif) .
Definition bound_hor.C:586
Lorene::Isol_hor::boundary_nn_Dir
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
Definition bound_hor.C:647
Lorene::Isol_hor::boundary_nn_Neu_kk
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
Definition bound_hor.C:420
Lorene::Isol_hor::boost_x
double boost_x
Boost velocity in x-direction.
Definition isol_hor.h:272
Lorene::Isol_hor::boundary_psi_Dir
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial).
Definition bound_hor.C:324
Lorene::Isol_hor::boundary_psi_Neu_evol
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution).
Definition bound_hor.C:203
Lorene::Isol_hor::source_nn
const Scalar source_nn() const
Source for N.
Definition sources_hor.C:145
Lorene::Metric::radial_vect
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition metric.C:362
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
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::Time_slice_conf::k_dd
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ).
Definition time_slice_conf.C:630
Lorene::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime ).
Definition time_slice_conf.C:591
Lorene::Time_slice_conf::ln_psi
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Definition time_slice_conf.C:719
Lorene::Time_slice_conf::psi
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Definition time_slice_conf.C:693
Lorene::Time_slice_conf::psi4
const Scalar & psi4() const
Factor at the current time step (jtime ).
Definition time_slice_conf.C:707
Lorene::Time_slice_conf::ff
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition time_slice.h:507
Lorene::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition time_slice_access.C:87
Lorene::Time_slice::gam
const Metric & gam() const
Induced metric at the current time step (jtime ).
Definition time_slice_access.C:95
Lorene::Time_slice::n_evol
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition time_slice.h:216
Lorene::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
Lorene::Time_slice::beta_evol
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:219
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
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
des_meridian
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
Lorene::arrete
void arrete(int a=0)
Setting a stop point in a code.
Definition arrete.C:61
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::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::maxabs
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Definition tensor_calculus_ext.C:606
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::dsdr
virtual void dsdr(const Cmp &ci, Cmp &resu) const =0
Computes of a Cmp .