app_hor_finder.C
/*
* Function ah_finder
*
* (see file app_hor.h for documentation)
*/
/*
* Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
*
*
* 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 app_hor_finder_C[] = "$Header: /cvsroot/Lorene/C++/Source/App_hor/app_hor_finder.C,v 1.11 2014/10/13 08:52:38 j_novak Exp $" ;
/*
* $Id: app_hor_finder.C,v 1.11 2014/10/13 08:52:38 j_novak Exp $
* $Log: app_hor_finder.C,v $
* Revision 1.11 2014/10/13 08:52:38 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.10 2014/10/06 15:12:56 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.9 2012/01/02 13:52:57 j_novak
* New parameter 'verbose' to get less output if needed.
*
* Revision 1.8 2008/01/08 13:56:54 j_novak
* Minor modif.
*
* Revision 1.7 2007/10/23 12:26:08 j_novak
* Added a test for the case where there is no AH, h(theta,phi) is then going out of the grid
*
* Revision 1.6 2005/12/09 09:35:06 lm_lin
*
* Add more information to screen output if no convergence.
*
* Revision 1.5 2005/12/07 14:16:36 lm_lin
*
* Add option to turn off screen output if no horizon is found
* (for performance reason in hydrodynamics simulation).
*
* Revision 1.4 2005/12/07 11:11:09 lm_lin
*
* Add option to turn off screen output during iterations.
*
* Revision 1.3 2005/11/17 15:53:28 lm_lin
*
* A tiny fix.
*
* Revision 1.2 2005/11/17 14:20:43 lm_lin
*
* Check the expansion function evaluated on the apparent horizon after the
* iteration of the 2-surface converges.
*
* Revision 1.1 2005/10/13 08:51:15 j_novak
* New stuff for apparent horizon finder. For the moment, there is only an
* external function. A class should come soon...
*
*
*
* $Header: /cvsroot/Lorene/C++/Source/App_hor/app_hor_finder.C,v 1.11 2014/10/13 08:52:38 j_novak Exp $
*
*/
// C headers
#include <cmath>
#include <cassert>
// Lorene headers
#include "app_hor.h"
#include "graphique.h"
namespace Lorene {
bool ah_finder(const Metric& gamma, const Sym_tensor& k_dd, Valeur& h, Scalar& ex_fcn,
double a_axis, double b_axis, double c_axis, bool verbose, bool print,
double precis, double precis_exp, int step_max, int step_relax,
double relax)
{
bool ah_flag = false ;
//Get the mapping, grid, base vector...etc
const Map& map = gamma.get_mp() ;
const Mg3d* mg = map.get_mg() ;
const Mg3d* g_angu = mg->get_angu() ;
int nz = mg->get_nzone() ;
const Base_vect_spher& bspher = map.get_bvect_spher() ;
const Coord& rr = map.r ;
const Coord& theta = map.tet ;
const Coord& phi = map.phi ;
const Coord& costh = map.cost ;
const Coord& cosph = map.cosp ;
const Coord& sinth = map.sint ;
const Coord& sinph = map.sinp ;
double r_min = min(+rr)(0) ;
double r_max = max(+rr)(nz-1) ;
// Set up a triaxial ellipsoidal surface as the initial guess for h
//------------------------------------------------------------------
double aa = a_axis ;
double bb = b_axis ;
double cc = c_axis ;
Scalar ct(map) ;
ct = costh ;
ct.std_spectral_base() ;
Scalar st(map) ;
st = sinth ;
st.std_spectral_base() ;
Scalar cp(map) ;
cp = cosph ;
cp.std_spectral_base() ;
Scalar sp(map) ;
sp = sinph ;
sp.std_spectral_base() ;
Scalar h_tmp(map) ;
h_tmp = st*st*( cp*cp/aa/aa + sp*sp/bb/bb ) + ct*ct/cc/cc ;
h_tmp = 1./sqrt( h_tmp ) ;
h_tmp.std_spectral_base() ;
Valeur h_guess(g_angu) ;
h_guess.annule_hard() ;
for (int l=0; l<nz; l++) {
int jmax = mg->get_nt(l) ;
int kmax = mg->get_np(l) ;
for (int k=0; k<kmax; k++) {
for (int j=0; j<jmax; j++) {
h_guess.set(l,k,j,0) = h_tmp.val_grid_point(l,k,j,0) ;
}
}
}
h_guess.std_base_scal() ;
h = h_guess ; // initialize h to h_guess
h.std_base_scal() ;
//End setting initial guess for h
//------------------------------------
const Metric_flat& fmets = map.flat_met_spher() ;
// Define the conformal factor
double one_third = double(1) / double(3) ;
Scalar psi4 = gamma.determinant() / fmets.determinant() ;
psi4 = pow( psi4, one_third ) ;
psi4.std_spectral_base() ;
// The expansion function at the n-1 iteration step
Scalar ex_fcn_old(map) ;
ex_fcn_old.set_etat_zero() ;
// The expansion function evaulated on the 2 surface h
Valeur ex_AH(g_angu) ;
ex_AH.annule_hard() ;
// Normal unit vector of the level set surface F = r - h(\theta,phi)
Vector s_unit(map, CON, bspher) ;
double relax_prev = double(1) - relax ;
double diff_exfcn = 1. ;
Tbl diff_h(nz) ;
diff_h = 1. ;
bool no_AH_in_grid = false ;
//--------------------------------------------------------
// Start of iteration
//--------------------------------------------------------
for (int step=0 ;
(max(diff_h) > precis) && (step < step_max) && (!no_AH_in_grid);
step++) {
// ***To be fixed: the function "set_grid_point" does not delete the derived
// quantities of F.
// Temporary fix: Define the level set F inside the iteration loop...
//----------------------------------------------------------------
Scalar F(map) ; // level set function: F = r - h(theta,phi)
F.allocate_all() ;
for (int l=0; l<nz; l++) {
int imax = mg->get_nr(l) ;
int jmax = mg->get_nt(l) ;
int kmax = mg->get_np(l) ;
for (int k=0; k<kmax; k++) {
for (int j=0; j<jmax; j++) {
for (int i=0; i<imax; i++) {
// (+rr) converts rr to Mtbl
F.set_grid_point(l,k,j,i) = (+rr)(l,k,j,i) - h(l,k,j,0) ;
}
}
}
}
F.std_spectral_base() ;
// Construct the unit normal vector s^i of the surface F
Scalar dF_norm(map) ;
dF_norm = contract( contract(gamma.con(), 0, F.derive_cov(gamma), 0),
0, F.derive_cov(gamma), 0) ;
dF_norm = sqrt( dF_norm ) ;
dF_norm.std_spectral_base() ;
s_unit = F.derive_con(gamma) / dF_norm ;
// The expansion function
ex_fcn = s_unit.divergence(gamma) - k_dd.trace(gamma) +
contract( s_unit, 0, contract(s_unit, 0, k_dd, 1), 0) ;
// Construct the source term for the angular Laplace equation
//---------------------------------------------------------
Sym_tensor sou_1(map, CON, bspher) ;
sou_1 = gamma.con() - fmets.con()/psi4 - s_unit*s_unit ;
Scalar source_tmp(map) ;
source_tmp = contract( sou_1, 0, 1, F.derive_cov(fmets).derive_cov(fmets), 0, 1 ) ;
source_tmp = source_tmp / dF_norm ;
Sym_tensor d_gam(map, COV, bspher) ;
d_gam = contract( gamma.connect().get_delta(), 0, F.derive_cov(fmets), 0) ;
source_tmp -= contract( gamma.con() - s_unit*s_unit, 0, 1,
d_gam, 0, 1) / dF_norm ;
source_tmp = psi4*dF_norm*( source_tmp - k_dd.trace(gamma) +
contract(s_unit, 0, contract(s_unit, 0, k_dd, 1), 0) ) ;
source_tmp.std_spectral_base() ;
Valeur sou_angu(g_angu) ; // source defined on the angular grid
// S(theta, phi) = S(h(theta,phi),theta,phi)
sou_angu.annule_hard() ;
double h_min = min(h)(0) ;
double h_max = max(h)(0) ;
if ( (r_min < h_min) && (h_max < r_max) ) {
for (int l=0; l<nz; l++) {
int jmax = mg->get_nt(l) ;
int kmax = mg->get_np(l) ;
for (int k=0; k<kmax; k++) {
for (int j=0; j<jmax; j++) {
sou_angu.set(l,k,j,0) = source_tmp.val_point(h(l,k,j,0)
,(+theta)(l,k,j,0) ,(+phi)(l,k,j,0)) ;
}
}
}
sou_angu = h*h*sou_angu ; // Final source term: psi4*dF_norm*h^2*(source_tmp)
}
else {
no_AH_in_grid = true ;
break ;
}
sou_angu.std_base_scal() ;
// Done with the source term
//-----------------------------------------
// Start solving the equation L^2h - 2h = source
//-----------------------------------------------
sou_angu.ylm() ;
Valeur h_new = sou_angu ;
h_new.c_cf->poisson_angu(-2.) ;
h_new.ylm_i() ;
if (h_new.c != 0x0)
delete h_new.c ;
h_new.c = 0x0 ;
h_new.coef_i() ;
// Convergence condition:
diff_h = max(abs(h - h_new)) ;
// Relaxations
if (step >= step_relax) {
h_new = relax * h_new + relax_prev * h ;
}
// Recycling for the next step
h = h_new ;
if (print)
{
cout << "-------------------------------------" << endl ;
cout << "App_hor iteration step: " << step << endl ;
cout << " " << endl ;
cout << "Difference in h : " << diff_h << endl ;
// Check: calculate the difference between ex_fcn and ex_fcn_old
Tbl diff_exfcn_tbl = diffrel( ex_fcn, ex_fcn_old ) ;
diff_exfcn = diff_exfcn_tbl(0) ;
for (int l=1; l<nz; l++) {
diff_exfcn += diff_exfcn_tbl(l) ;
}
diff_exfcn /= nz ;
cout << "diff_exfcn : " << diff_exfcn << endl ;
ex_fcn_old = ex_fcn ; // recycling
// End check
}
if ( (step == step_max-1) && (max(diff_h) > precis) ) {
//Check: Evaluate the expansion function on the 2-surface
for (int l=0; l<nz; l++) {
int jmax = mg->get_nt(l) ;
int kmax = mg->get_np(l) ;
for (int k=0; k<kmax; k++) {
for (int j=0; j<jmax; j++) {
ex_AH.set(l,k,j,0) = ex_fcn.val_point(h(l,k,j,0),(+theta)(l,k,j,0)
,(+phi)(l,k,j,0)) ;
}
}
}
if (verbose) {
cout << " " << endl ;
cout << "###############################################" << endl ;
cout << "AH finder: maximum number of iteration reached!" << endl ;
cout << " No convergence in the 2-surface h! " << endl ;
cout << " max( difference in h ) > prescribed tolerance " << endl ;
cout << " " << endl ;
cout << " prescribed tolerance = " << precis << endl ;
cout << " max( difference in h ) = " << max(diff_h) << endl ;
cout << " max( expansion function on h ) = " << max(abs(ex_AH(0))) << endl ;
cout << "###############################################" << endl ;
cout << " " << endl ;
}
}
} // End of iteration
//Done with the AH finder
//Check: Evaluate the expansion function on the 2-surface
if (no_AH_in_grid) {
if (print) {
cout << " " << endl ;
cout << "###############################################" << endl ;
cout << " AH finder: no horizon found inside the grid!" << endl ;
cout << " Grid: rmin= " << r_min << ", rmax= " << r_max << endl ;
cout << "###############################################" << endl ;
cout << " " << endl ;
}
}
else {
for (int l=0; l<nz; l++) {
int jmax = mg->get_nt(l) ;
int kmax = mg->get_np(l) ;
for (int k=0; k<kmax; k++) {
for (int j=0; j<jmax; j++) {
ex_AH.set(l,k,j,0) = ex_fcn.val_point(h(l,k,j,0),(+theta)(l,k,j,0)
,(+phi)(l,k,j,0)) ;
}
}
}
if ( (max(diff_h) < precis) && (max(abs(ex_AH(0))) < precis_exp) ) {
ah_flag = true ;
if (verbose) {
cout << " " << endl ;
cout << "################################################" << endl ;
cout << " AH finder: Apparent horizon found!!! " << endl ;
cout << " Max error of the expansion function on h: " << endl ;
cout << " max( expansion function on AH ) = " << max(abs(ex_AH(0))) << endl ;
cout << "################################################" << endl ;
cout << " " << endl ;
}
}
if ( (max(diff_h) < precis) && (max(abs(ex_AH(0))) > precis_exp) ) {
if (print) {
cout << " " << endl ;
cout << "#############################################" << endl ;
cout << " AH finder: convergence in the 2 surface h. " << endl ;
cout << " But max error of the expansion function evaulated on h > precis_exp" << endl ;
cout << " max( expansion function on AH ) = " << max(abs(ex_AH(0))) << endl ;
cout << " Probably not an apparent horizon! " << endl ;
cout << "#############################################" << endl ;
cout << " " << endl ;
}
}
}
return ah_flag ;
} // End ah_finder
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Connection::get_delta
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition connection.h:271
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Metric_flat::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric_flat.C:153
Lorene::Metric_flat::determinant
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric_flat.C:214
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Lorene::Metric::get_mp
const Map & get_mp() const
Returns the mapping.
Definition metric.h:202
Lorene::Metric::connect
virtual const Connection & connect() const
Returns the connection.
Definition metric.C:301
Lorene::Metric::determinant
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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::Mtbl_cf::poisson_angu
void poisson_angu(double lambda=0)
Resolution of the generalized angular Poisson equation.
Definition mtbl_cf_pde.C:83
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition scalar_deriv.C:390
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::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::allocate_all
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
Lorene::Scalar::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:890
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::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c
Mtbl * c
Values of the function at the points of the multi-grid.
Definition valeur.h:299
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Lorene::Valeur::coef_i
void coef_i() const
Computes the physical value of *this.
Definition valeur_coef_i.C:140
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
Lorene::Valeur::std_base_scal
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
Lorene::Valeur::annule_hard
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
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::ah_finder
bool ah_finder(const Metric &gamma, const Sym_tensor &k_dd_in, Valeur &h, Scalar &ex_fcn, double a_axis, double b_axis, double c_axis, bool verbose=true, bool print=false, double precis=1.e-8, double precis_exp=1.e-6, int it_max=200, int it_relax=200, double relax_fac=1.)
Class for apparent horizon (under heavy development).
Definition app_hor_finder.C:91
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition tensor_calculus.C:94
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::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720