sym_tensor_trans_aux.C
/*
* Manipulation of auxiliary potentials for Sym_tensor_trans objects.
*
* (see file sym_tensor.h for documentation).
*
*/
/*
* Copyright (c) 2005-2006 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 sym_tensor_trans_aux_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $" ;
/*
* $Id: sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $
* $Log: sym_tensor_trans_aux.C,v $
* Revision 1.21 2014/10/13 08:53:43 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.20 2014/10/06 15:13:19 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.19 2010/10/11 10:23:03 j_novak
* Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
*
* Revision 1.18 2008/12/05 08:46:19 j_novak
* New method Sym_tensor_trans_aux::set_tt_part_det_one.
*
* Revision 1.17 2007/04/27 13:52:55 j_novak
* In method Sym_tensor_trans::set_AtBtt_det_one(), the member p_tt (the TT-part)
* is updated.
*
* Revision 1.16 2007/03/20 12:20:56 j_novak
* In Sym_tensor_trans::set_AtBtt_det_one(), the trace is stored in the resulting
* object.
*
* Revision 1.15 2006/10/24 13:03:19 j_novak
* New methods for the solution of the tensor wave equation. Perhaps, first
* operational version...
*
* Revision 1.14 2006/09/15 08:48:01 j_novak
* Suppression of some messages in the optimized version.
*
* Revision 1.13 2006/08/31 12:13:22 j_novak
* Added an argument of type Param to Sym_tensor_trans::sol_Dirac_A().
*
* Revision 1.12 2006/06/26 09:28:13 j_novak
* Added a forgotten initialisation in set_AtB_trace_zero().
*
* Revision 1.11 2006/06/20 12:07:15 j_novak
* Improved execution speed for sol_Dirac_tildeB...
*
* Revision 1.10 2006/06/14 10:04:21 j_novak
* New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
*
* Revision 1.9 2006/01/17 15:50:53 j_novak
* Slight re-arrangment of the det=1 formula.
*
* Revision 1.8 2005/11/28 14:45:17 j_novak
* Improved solution of the Poisson tensor equation in the case of a transverse
* tensor.
*
* Revision 1.7 2005/09/16 13:34:44 j_novak
* Back to dzpuis 1 for the source for mu. eta is computed the same way as hrr.
*
* Revision 1.6 2005/09/08 16:00:23 j_novak
* Change of dzpuis for source for mu (temporary?).
*
* Revision 1.5 2005/09/07 16:47:43 j_novak
* Removed method Sym_tensor_trans::T_from_det_one
* Modified Sym_tensor::set_auxiliary, so that it takes eta/r and mu/r as
* arguments.
* Modified Sym_tensor_trans::set_hrr_mu.
* Added new protected method Sym_tensor_trans::solve_hrr
*
* Revision 1.4 2005/04/08 08:22:04 j_novak
* New methods set_hrr_mu_det_one() and set_WX_det_one(). Not tested yet...
*
* Revision 1.3 2005/04/07 07:56:22 j_novak
* Better handling of dzpuis (first try).
*
* Revision 1.2 2005/04/06 15:49:46 j_novak
* Error corrected.
*
* Revision 1.1 2005/04/06 15:43:59 j_novak
* New method Sym_tensor_trans::T_from_det_one(...).
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $
*
*/
// C headers
#include <cassert>
// Lorene headers
#include "metric.h"
#include "param.h"
namespace Lorene {
void Sym_tensor_trans::set_hrr_mu_det_one(const Scalar& hrr, const Scalar& mu_in,
double precis, int it_max ) {
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
assert(hrr.check_dzpuis(0)) ;
assert(mu_in.check_dzpuis(0)) ;
assert(&mu_in != p_mu) ;
assert( (precis > 0.) && (it_max > 0) ) ;
Sym_tensor_tt tens_tt(*mp, *triad, *met_div) ;
tens_tt.set_rr_mu(hrr, mu_in) ;
tens_tt.inc_dzpuis(2) ;
trace_from_det_one(tens_tt, precis, it_max) ;
dec_dzpuis(2) ;
return ;
}
void Sym_tensor_trans::set_AtBtt_det_one(const Scalar& a_in, const Scalar& tbtt_in,
const Scalar* h_prev, Param* par_bc,
Param* par_mat, double precis,
int it_max ) {
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
assert(a_in.check_dzpuis(2)) ;
assert(tbtt_in.check_dzpuis(2)) ;
assert(&a_in != p_aaa) ;
assert(&tbtt_in != p_tilde_b) ;
assert( (precis > 0.) && (it_max > 0) ) ;
//Computation of mu and X from A
//-------------------------------
Scalar mu_over_r(*mp) ;
Scalar x_new(*mp) ;
sol_Dirac_A(a_in, mu_over_r, x_new, par_bc) ;
// Preparation for the iteration
//------------------------------
Scalar h_old(*mp) ;
if (h_prev != 0x0)
h_old = *h_prev ;
else
h_old.set_etat_zero() ;
double lambda = 1. ;
Scalar zero(*mp) ;
zero.set_etat_zero() ;
Scalar hrr_tt(*mp) ;
Scalar eta_sr_tt(*mp) ;
Scalar w_tt(*mp) ;
sol_Dirac_tilde_B(tbtt_in, zero, hrr_tt, eta_sr_tt, w_tt, par_bc, par_mat) ;
Sym_tensor_tt hijtt(*mp, *triad, *met_div) ;
hijtt.set_auxiliary(hrr_tt, eta_sr_tt, mu_over_r, w_tt, x_new, zero) ;
Scalar hrr_new(*mp) ;
Scalar eta_over_r(*mp) ;
Scalar w_new(*mp) ;
for (int it=0; it<=it_max; it++) {
Scalar tilde_B = get_tilde_B_from_TT_trace(zero, h_old) ;
sol_Dirac_tilde_B(tilde_B, h_old, hrr_new, eta_over_r, w_new, 0x0, par_mat) ;
set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
const Sym_tensor_trans& hij = *this ;
Scalar h_new = (1 + hij(1,1))*( hij(2,3)*hij(2,3) - hij(2,2)*hij(3,3) )
+ hij(1,2)*hij(1,2)*(1 + hij(3,3))
+ hij(1,3)*hij(1,3)*(1 + hij(2,2))
- hij(1,1)*(hij(2,2) + hij(3,3)) - 2*hij(1,2)*hij(1,3)*hij(2,3) ;
h_new.set_spectral_base(hrr_tt.get_spectral_base()) ;
double diff = max(max(abs(h_new - h_old))) ;
#ifndef NDEBUG
cout << "Sym_tensor_trans::set_AtB_det_one : "
<< "iteration : " << it << " convergence on h: "
<< diff << endl ;
#endif
if (diff < precis) {
set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
if (p_aaa != 0x0) delete p_aaa ;
p_aaa = new Scalar(a_in) ;
if (p_tilde_b != 0x0) delete p_tilde_b ;
p_tilde_b = new Scalar(tilde_B + tbtt_in) ;
if (p_trace != 0x0) delete p_trace ;
p_trace = new Scalar(h_old) ;
if (p_tt != 0x0) delete p_tt ;
p_tt = new Sym_tensor_tt(hijtt) ;
p_tt->inc_dzpuis(2) ;
break ;
}
else {
h_old = lambda*h_new +(1-lambda)*h_old ;
}
if (it == it_max) {
cout << "Sym_tensor_trans:::set_AtBtt_det_one : convergence not reached \n" ;
cout << " for the required accuracy (" << precis << ") ! "
<< endl ;
abort() ;
}
}
return ;
}
void Sym_tensor_trans::set_tt_part_det_one(const Sym_tensor_tt& hijtt, const
Scalar* h_prev, Param* par_mat,
double precis, int it_max ) {
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
assert( (precis > 0.) && (it_max > 0) ) ;
Scalar mu_over_r = hijtt.mu() ;
mu_over_r.div_r() ;
const Scalar& x_new = hijtt.xxx() ;
// Preparation for the iteration
//------------------------------
Scalar h_old(*mp) ;
if (h_prev != 0x0)
h_old = *h_prev ;
else
h_old.set_etat_zero() ;
double lambda = 1. ;
Scalar zero(*mp) ;
zero.set_etat_zero() ;
const Scalar& hrr_tt = hijtt( 1, 1 ) ;
Scalar eta_sr_tt = hijtt.eta() ;
eta_sr_tt.div_r() ;
const Scalar w_tt = hijtt.www() ;
Scalar hrr_new(*mp) ;
Scalar eta_over_r(*mp) ;
Scalar w_new(*mp) ;
for (int it=0; it<=it_max; it++) {
Scalar tilde_B = get_tilde_B_from_TT_trace(zero, h_old) ;
sol_Dirac_tilde_B(tilde_B, h_old, hrr_new, eta_over_r, w_new, 0x0, par_mat) ;
set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
const Sym_tensor_trans& hij = *this ;
Scalar h_new = (1 + hij(1,1))*( hij(2,3)*hij(2,3) - hij(2,2)*hij(3,3) )
+ hij(1,2)*hij(1,2)*(1 + hij(3,3))
+ hij(1,3)*hij(1,3)*(1 + hij(2,2))
- hij(1,1)*(hij(2,2) + hij(3,3)) - 2*hij(1,2)*hij(1,3)*hij(2,3) ;
h_new.set_spectral_base(hrr_tt.get_spectral_base()) ;
double diff = max(max(abs(h_new - h_old))) ;
#ifndef NDEBUG
cout << "Sym_tensor_trans::set_tt_part_det_one : "
<< "iteration : " << it << " convergence on h: "
<< diff << endl ;
#endif
if (diff < precis) {
set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
if (p_trace != 0x0) delete p_trace ;
p_trace = new Scalar(h_old) ;
if (p_tt != 0x0) delete p_tt ;
p_tt = new Sym_tensor_tt(hijtt) ;
p_tt->inc_dzpuis(2) ;
break ;
}
else {
h_old = lambda*h_new +(1-lambda)*h_old ;
}
if (it == it_max) {
cout << "Sym_tensor_trans:::set_AtBtt_det_one : convergence not reached \n" ;
cout << " for the required accuracy (" << precis << ") ! "
<< endl ;
abort() ;
}
}
return ;
}
void Sym_tensor_trans::set_AtB_trace(const Scalar& a_in, const Scalar& tb_in,
const Scalar& hh, Param* par_bc, Param* par_mat) {
// All this has a meaning only for spherical components:
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
assert(a_in.check_dzpuis(2)) ;
assert(tb_in.check_dzpuis(2)) ;
assert(hh.check_dzpuis(0)) ;
assert(&a_in != p_aaa) ;
assert(&tb_in != p_tilde_b) ;
//Computation of mu and X from A
//-------------------------------
Scalar mu_over_r(*mp) ;
Scalar x_new(*mp) ;
sol_Dirac_A(a_in, mu_over_r, x_new, par_bc) ;
// Computation of the other potentials
//------------------------------------
Scalar hrr_new(*mp) ;
Scalar eta_over_r(*mp) ;
Scalar w_new(*mp) ;
sol_Dirac_tilde_B(tb_in, hh, hrr_new, eta_over_r, w_new, par_bc, par_mat) ;
set_auxiliary(hrr_new, eta_over_r, mu_over_r, w_new, x_new, hh - hrr_new) ;
if (p_aaa != 0x0) delete p_aaa ;
p_aaa = new Scalar(a_in) ;
if (p_tilde_b != 0x0) delete p_tilde_b ;
p_tilde_b = new Scalar(tb_in) ;
return ;
}
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
Lorene::Scalar::div_r
void div_r()
Division by r everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:133
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition scalar.h:1294
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition scalar.C:797
Lorene::Sym_tensor_trans::p_tt
Sym_tensor_tt * p_tt
Traceless part with respect to the metric *met_div.
Definition sym_tensor.h:620
Lorene::Sym_tensor_trans::set_AtBtt_det_one
void set_AtBtt_det_one(const Scalar &a_in, const Scalar &tbtt_in, const Scalar *h_prev=0x0, Param *par_bc=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assigns the derived member A and computes from its TT-part (see Sym_tensor::compute_tilde_B_tt() ).
Definition sym_tensor_trans_aux.C:137
Lorene::Sym_tensor_trans::sol_Dirac_tilde_B
void sol_Dirac_tilde_B(const Scalar &tilde_b, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &www, Param *par_bc=0x0, Param *par_mat=0x0) const
Solves a system of three coupled first-order PDEs obtained from divergence-free conditions (Dirac gau...
Definition sym_tensor_trans_dirac.C:583
Lorene::Sym_tensor_trans::Sym_tensor_trans
Sym_tensor_trans(const Map &map, const Base_vect &triad_i, const Metric &met)
Standard constructor.
Definition sym_tensor_trans.C:114
Lorene::Sym_tensor_trans::trace_from_det_one
void trace_from_det_one(const Sym_tensor_tt &htt, double precis=1.e-14, int it_max=100)
Assigns the derived member p_tt and computes the trace so that *this + the flat metric has a determin...
Definition sym_tensor_trans.C:315
Lorene::Sym_tensor_trans::met_div
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition sym_tensor.h:614
Lorene::Sym_tensor_trans::p_trace
Scalar * p_trace
Trace with respect to the metric *met_div.
Definition sym_tensor.h:617
Lorene::Sym_tensor_trans::set_AtB_trace
void set_AtB_trace(const Scalar &a_in, const Scalar &tb_in, const Scalar &trace, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A , and the trace.
Definition sym_tensor_trans_aux.C:302
Lorene::Sym_tensor_trans::set_tt_part_det_one
void set_tt_part_det_one(const Sym_tensor_tt &hijtt, const Scalar *h_prev=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assignes the TT-part of the tensor.
Definition sym_tensor_trans_aux.C:226
Lorene::Sym_tensor_trans::sol_Dirac_A
void sol_Dirac_A(const Scalar &aaa, Scalar &tilde_mu, Scalar &xxx, const Param *par_bc=0x0) const
Solves a system of two coupled first-order PDEs obtained from the divergence-free condition (Dirac ga...
Definition sym_tensor_trans_dirac.C:82
Lorene::Sym_tensor_trans::set_hrr_mu_det_one
void set_hrr_mu_det_one(const Scalar &hrr, const Scalar &mu_in, double precis=1.e-14, int it_max=100)
Assigns the rr component and the derived member .
Definition sym_tensor_trans_aux.C:117
Lorene::Sym_tensor_tt
Transverse and traceless symmetric tensors of rank 2.
Definition sym_tensor.h:938
Lorene::Sym_tensor_tt::set_rr_mu
void set_rr_mu(const Scalar &hrr, const Scalar &mu_i)
Sets the component , as well as the angular potential (see member p_mu ).
Definition sym_tensor_tt_etamu.C:238
Lorene::Sym_tensor_tt::eta
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
Definition sym_tensor_tt_etamu.C:140
Lorene::Sym_tensor::xxx
const Scalar & xxx() const
Gives the field X (see member p_xxx ).
Definition sym_tensor_aux.C:240
Lorene::Sym_tensor::p_aaa
Scalar * p_aaa
Field A defined from X and insensitive to the longitudinal part of the Sym_tensor (only for ).
Definition sym_tensor.h:322
Lorene::Sym_tensor::p_tilde_b
Scalar * p_tilde_b
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition sym_tensor.h:334
Lorene::Sym_tensor::p_mu
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:274
Lorene::Sym_tensor::www
const Scalar & www() const
Gives the field W (see member p_www ).
Definition sym_tensor_aux.C:209
Lorene::Sym_tensor::get_tilde_B_from_TT_trace
Scalar get_tilde_B_from_TT_trace(const Scalar &tilde_B_tt_in, const Scalar &trace) const
Computes (see Sym_tensor::p_tilde_b ) from its transverse-traceless part and the trace.
Definition sym_tensor_aux.C:531
Lorene::Sym_tensor::mu
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
Definition sym_tensor_aux.C:151
Lorene::Sym_tensor::set_auxiliary
void set_auxiliary(const Scalar &trr, const Scalar &eta_over_r, const Scalar &mu_over_r, const Scalar &www, const Scalar &xxx, const Scalar &ttt)
Assigns the component and the derived members p_eta , p_mu , p_www, p_xxx and p_ttt ,...
Definition sym_tensor_aux.C:266
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
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::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::dec_dzpuis
virtual void dec_dzpuis(Scalar &) const =0
Decreases by 1 the value of dzpuis of a Scalar and changes accordingly its values in the compacti...