sym_ttt_poisson.C
/*
* Resolution of the TT tensor Poisson equation
*
* (see file sym_tensor.h for documentation).
*
*/
/*
* Copyright (c) 2003 Eric Gourgoulhon & 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_ttt_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $" ;
/*
* $Id: sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $
* $Log: sym_ttt_poisson.C,v $
* Revision 1.5 2014/10/13 08:53:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2004/12/28 10:37:24 j_novak
* Better way of enforcing zero divergence.
*
* Revision 1.3 2004/12/27 14:33:12 j_novak
* New algorithm for the tensor Poisson eq.
*
* Revision 1.2 2004/03/04 09:50:41 e_gourgoulhon
* Method poisson: use of new methods khi() and set_khi_mu.
*
* Revision 1.1 2003/11/07 16:53:52 e_gourgoulhon
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $
*
*/
// C headers
//#include <>
// Lorene headers
#include "tensor.h"
#include "param_elliptic.h"
namespace Lorene {
Sym_tensor_tt Sym_tensor_tt::poisson(int dzfin) const {
// All this has a meaning only for spherical components...
assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
//## ... and affine mapping, for the moment!
assert(dynamic_cast<const Map_af*>(mp) != 0x0) ;
assert( (dzfin == 0) || (dzfin == 2) ) ;
Sym_tensor_tt resu(*mp, *triad, *met_div) ;
// Solution for the rr-component
// ----------------------------
const Scalar& source_rr = operator()(1,1) ;
Scalar h_rr(*mp) ;
int nz = mp->get_mg()->get_nzone() ;
if (source_rr.get_etat() != ETATZERO) {
//------------------------
// The elliptic operator
//------------------------
Param_elliptic param_hr(source_rr) ;
for (int lz=0; lz<nz; lz++)
param_hr.set_poisson_tens_rr(lz) ;
h_rr = source_rr.sol_elliptic(param_hr) ;
}
else
h_rr.set_etat_zero() ;
h_rr.inc_dzpuis(dzfin) ; //## can we improve here?
resu.set(1,1) = h_rr ;
// Solution for (eta / r)
//-----------------------
// Scalar source_eta = - source_rr ;
// source_eta.mult_r_dzpuis(3) ;
// source_eta.mult_r_dzpuis(2) ;
// h_rr.set_spectral_va().ylm() ;
// Scalar tmp = 2*h_rr + h_rr.lapang() ;
// if (dzfin == 0)
// tmp.inc_dzpuis(2) ;
// source_eta += tmp ;
// source_eta = source_eta.primr() ;
// source_eta.div_r_dzpuis(dzfin) ;
// Scalar etasurr = (h_rr+source_eta).poisson_angu() ;
Scalar source_eta = -resu(1,1).dsdr() ;
source_eta.mult_r_dzpuis(dzfin) ;
source_eta -= 3.*resu(1,1) ;
Scalar etasurr = source_eta.poisson_angu() ;
// Solution for mu
// ---------------
Scalar musurr = mu().poisson() ;
musurr.div_r_dzpuis(dzfin) ;
resu.set(1,1).set_spectral_va().ylm_i() ;
Scalar** rcmp = resu.cmp ;
Itbl idx(2) ;
idx.set(0) = 1 ; // r index
// h^{r theta} :
// ------------
idx.set(1) = 2 ; // theta index
*rcmp[position(idx)] = etasurr.dsdt() - musurr.stdsdp() ;
// h^{r phi} :
// ------------
idx.set(1) = 3 ; // phi index
*rcmp[position(idx)] = etasurr.stdsdp() + musurr.dsdt() ;
// h^{theta phi} and h^{phi phi}
// -----------------------------
//-------------- Computation of T^theta --> taut :
Scalar tautst = resu(1,2).dsdr() ;
// dhrr contains dh^{rt}/dr in all domains but the CED,
// in the CED: r^2 dh^{rt}/dr if dzfin = 0 (1)
// r^3 dh^{rt}/dr if dzfin = 2 (2)
// Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
tautst.mult_r_dzpuis(dzfin) ;
// Addition of the remaining parts :
tautst += 3 * resu(1,2) - resu(1,1).dsdt() ;
tautst.mult_sint() ;
Scalar tmp = resu(1,1) ;
tmp.mult_cost() ; // h^{rr} cos(th)
tautst -= tmp ; // T^th / sin(th)
Scalar taut = tautst ;
taut.mult_sint() ; // T^th
//----------- Computation of T^phi --> taup :
Scalar taupst = - resu(1,3).dsdr() ;
// dhrr contains - dh^{rp}/dr in all domains but the CED,
// in the CED: - r^2 dh^{rp}/dr if dzfin = 0 (3)
// - r^3 dh^{rp}/dr if dzfin = 2 (4)
// Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
taupst.mult_r_dzpuis(dzfin) ;
// Addition of the remaining part :
taupst -= 3 * resu(1,3) ;
taupst.mult_sint() ; // T^ph / sin(th)
Scalar taup = taupst ;
taup.mult_sint() ; // T^ph
//------------------- Computation of F and h^[ph ph}
tmp = tautst ;
tmp.mult_cost() ; // T^th / tan(th)
// dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
tmp = taut.dsdt() + tmp + taup.stdsdp() ;
Scalar tmp2 (*mp) ;
tmp2 = tmp.poisson_angu() ; // F
tmp2.div_sint() ;
tmp2.div_sint() ; // h^{ph ph}
idx.set(0) = 3 ; // phi index
idx.set(1) = 3 ; // phi index
*rcmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
//------------------- Computation of G and h^[th ph}
tmp = taupst ;
tmp.mult_cost() ; // T^ph / tan(th)
// - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
tmp2 = tmp.poisson_angu() ; // G
tmp2.div_sint() ;
tmp2.div_sint() ; // h^{th ph}
idx.set(0) = 2 ; // theta index
idx.set(1) = 3 ; // phi index
*rcmp[position(idx)] = tmp2 ; // h^{th ph} is updated
// h^{th th} (from the trace-free condition)
// ---------
idx.set(1) = 2 ; // theta index
*rcmp[position(idx)] = - resu(1,1) - resu(3,3) ;
return resu ;
}
}
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition param_elliptic.h:135
Lorene::Param_elliptic::set_poisson_tens_rr
void set_poisson_tens_rr(int zone)
Sets the operator to in all domains, for only.
Definition param_elliptic.C:560
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:147
Lorene::Scalar::dsdt
const Scalar & dsdt() const
Returns of *this .
Definition scalar_deriv.C:208
Lorene::Scalar::div_sint
void div_sint()
Division by .
Definition scalar_th_manip.C:98
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
Lorene::Scalar::mult_cost
void mult_cost()
Multiplication by .
Definition scalar_th_manip.C:58
Lorene::Scalar::poisson_angu
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition scalar_pde.C:200
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::stdsdp
const Scalar & stdsdp() const
Returns of *this .
Definition scalar_deriv.C:238
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::mult_sint
void mult_sint()
Multiplication by .
Definition scalar_th_manip.C:84
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Lorene::Scalar::sol_elliptic
Scalar sol_elliptic(Param_elliptic &params) const
Resolution of a general elliptic equation, putting zero at infinity.
Definition scalar_pde.C:234
Lorene::Scalar::mult_r_dzpuis
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:221
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_tt::poisson
Sym_tensor_tt poisson(int dzfin=2) const
Computes the solution of a tensorial TT Poisson equation with *this as a source:
Definition sym_ttt_poisson.C:62
Lorene::Sym_tensor_tt::Sym_tensor_tt
Sym_tensor_tt(const Map &map, const Base_vect &triad_i, const Metric &met)
Standard constructor.
Definition sym_tensor_tt.C:75
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::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
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::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
Lorene::Tensor_sym::position
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor_sym.C:245
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
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