map_af_elliptic.C
/*
* Copyright (c) 1999-2001 Eric Gourgoulhon
* Copyright (c) 1999-2001 Philippe Grandclement
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char map_af_elliptic_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $" ;
/*
* $Id: map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
* $Log: map_af_elliptic.C,v $
* Revision 1.13 2014/10/13 08:53:02 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/10/06 15:13:12 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.11 2007/05/06 10:48:11 p_grandclement
* Modification of a few operators for the vorton project
*
* Revision 1.10 2007/01/16 15:08:07 n_vasset
* New constructor, usn Scalar on mono-domain angular grid for boundary,
* for function sol_elliptic_boundary()
*
* Revision 1.9 2005/11/30 11:09:07 p_grandclement
* Changes for the Bin_ns_bh project
*
* Revision 1.8 2005/08/26 14:02:40 p_grandclement
* Modification of the elliptic solver that matches with an oscillatory exterior solution
* small correction in Poisson tau also...
*
* Revision 1.7 2005/06/09 07:57:31 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.6 2004/08/24 09:14:42 p_grandclement
* Addition of some new operators, like Poisson in 2d... It now requieres the
* GSL library to work.
*
* Also, the way a variable change is stored by a Param_elliptic is changed and
* no longer uses Change_var but rather 2 Scalars. The codes using that feature
* will requiere some modification. (It should concern only the ones about monopoles)
*
* Revision 1.5 2004/06/22 08:49:58 p_grandclement
* Addition of everything needed for using the logarithmic mapping
*
* Revision 1.4 2004/03/17 15:58:48 p_grandclement
* Slight modification of sol_elliptic_no_zec
*
* Revision 1.3 2004/02/11 09:47:46 p_grandclement
* Addition of a new elliptic solver, matching with the homogeneous solution
* at the outer shell and not solving in the external domain (more details
* coming soon ; check your local Lorene dealer...)
*
* Revision 1.2 2004/01/28 16:46:23 p_grandclement
* Addition of the sol_elliptic_fixe_der_zero stuff
*
* Revision 1.1 2003/12/11 14:48:48 p_grandclement
* Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
*
*
* $Header: /cvsroot/Lorene/C++/Source/Map/map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
*
*/
// Header C :
#include <cstdlib>
#include <cmath>
// Headers Lorene :
#include "tbl.h"
#include "mtbl_cf.h"
#include "map.h"
#include "param_elliptic.h"
#include "proto.h"
//----------------------------------------------
// General elliptic solver
//----------------------------------------------
namespace Lorene {
void Map_af::sol_elliptic(Param_elliptic& ope_var, const Scalar& source,
Scalar& pot) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
Mtbl_cf resu = elliptic_solver (ope_var, *(rho.c_cf)) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//-----------------------------------------------
// General elliptic solver with boundary as Mtbl-cf
//-------------------------------------------------
void Map_af::sol_elliptic_boundary(Param_elliptic& ope_var, const Scalar& source,
Scalar& pot, const Mtbl_cf& bound,
double fact_dir, double fact_neu ) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound,
fact_dir, fact_neu) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//-----------------------------------------------
// General elliptic solver with boundary as Scalar
//-------------------------------------------------
void Map_af::sol_elliptic_boundary(Param_elliptic& ope_var, const Scalar& source,
Scalar& pot, const Scalar& bound,
double fact_dir, double fact_neu ) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
// REMINDER : The scalar bound must be defined on a mono-domain angular grid corresponding with the full grid of the scalar source (following assert())
Scalar bbound = bound;
bbound.set_spectral_va().ylm() ;
const Map& mapp = bbound.get_mp();
const Mg3d& gri2d = *mapp.get_mg();
assert(&gri2d == source.get_mp().get_mg()->get_angu_1dom()) ; // Attention à cet assert !!
Mtbl_cf bound2 (gri2d , bbound.get_spectral_base()) ;
bound2.annule_hard() ;
if (bbound.get_etat() != ETATZERO){
int nr = gri2d.get_nr(0) ;
int nt = gri2d.get_nt(0) ;
int np = gri2d.get_np(0) ;
for(int k=0; k<np+2; k++)
for (int j=0; j<=nt-1; j++)
for(int xi=0; xi<= nr-1; xi++)
{
bound2.set(0, k , j , xi) = (*bbound.get_spectral_va().c_cf)(0, k, j, xi) ;
}
}
Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound2,
fact_dir, fact_neu) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//----------------------------------------------
// General elliptic solver with no ZEC
//----------------------------------------------
void Map_af::sol_elliptic_no_zec(Param_elliptic& ope_var, const Scalar& source,
Scalar& pot, double val) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
Mtbl_cf resu = elliptic_solver_no_zec (ope_var, *(rho.c_cf), val) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//----------------------------------------------
// General elliptic solver only in the ZEC
//----------------------------------------------
void Map_af::sol_elliptic_only_zec(Param_elliptic& ope_var,
const Scalar& source,
Scalar& pot, double val) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
Mtbl_cf resu = elliptic_solver_only_zec (ope_var, *(rho.c_cf), val) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//----------------------------------------------
// General elliptic solver with no ZEC
// and a mtaching with sin(r)/r
//----------------------------------------------
void Map_af::sol_elliptic_sin_zec(Param_elliptic& ope_var,
const Scalar& source, Scalar& pot, double* amplis, double* phases) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
Mtbl_cf resu = elliptic_solver_sin_zec (ope_var, *(rho.c_cf), amplis, phases) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
//----------------------------------------------
// General elliptic solver with no ZEC
//----------------------------------------------
void Map_af::sol_elliptic_fixe_der_zero (double valeur,
Param_elliptic& ope_var,
const Scalar& source,
Scalar& pot) const {
assert(source.get_etat() != ETATNONDEF) ;
assert(source.get_mp().get_mg() == mg) ;
assert(pot.get_mp().get_mg() == mg) ;
assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
source.check_dzpuis(4)) ;
// Spherical harmonic expansion of the source
// ------------------------------------------
const Valeur& sourva = source.get_spectral_va() ;
if (sourva.get_etat() == ETATZERO) {
pot.set_etat_zero() ;
return ;
}
// Spectral coefficients of the source
assert(sourva.get_etat() == ETATQCQ) ;
Valeur rho(sourva.get_mg()) ;
sourva.coef() ;
rho = *(sourva.c_cf) ; // copy of the coefficients of the source
rho.ylm() ; // spherical harmonic transforms
// On met les bonnes bases dans le changement de variable
// de ope_var :
ope_var.var_F.set_spectral_va().coef() ;
ope_var.var_F.set_spectral_va().ylm() ;
ope_var.var_G.set_spectral_va().coef() ;
ope_var.var_G.set_spectral_va().ylm() ;
// Call to the Mtbl_cf version
// ---------------------------
valeur *= alpha[0] ;
Mtbl_cf resu = elliptic_solver_fixe_der_zero (valeur, ope_var, *(rho.c_cf)) ;
// Final result returned as a Scalar
// ---------------------------------
pot.set_etat_zero() ; // to call Scalar::del_t().
pot.set_etat_qcq() ;
pot.set_spectral_va() = resu ;
pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
pot.set_dzpuis(0) ;
}
}
Lorene::Map_af::sol_elliptic_boundary
void sol_elliptic_boundary(Param_elliptic &params, const Scalar &so, Scalar &uu, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
General elliptic solver including inner boundary conditions.
Definition map_af_elliptic.C:159
Lorene::Map_af::sol_elliptic
void sol_elliptic(Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver.
Definition map_af_elliptic.C:100
Lorene::Map_af::sol_elliptic_sin_zec
void sol_elliptic_sin_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double *coefs, double *) const
General elliptic solver.
Definition map_af_elliptic.C:426
Lorene::Map_af::sol_elliptic_fixe_der_zero
void sol_elliptic_fixe_der_zero(double val, Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver fixing the derivative at the origin and relaxing the continuity of the first ...
Definition map_af_elliptic.C:483
Lorene::Map_af::alpha
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
Lorene::Map_af::sol_elliptic_only_zec
void sol_elliptic_only_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double val) const
General elliptic solver.
Definition map_af_elliptic.C:368
Lorene::Map_af::sol_elliptic_no_zec
void sol_elliptic_no_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double val) const
General elliptic solver.
Definition map_af_elliptic.C:312
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
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_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
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Lorene::Mtbl_cf::set
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition param_elliptic.h:135
Lorene::Param_elliptic::var_G
Scalar var_G
Multiplicative variable change that must be sphericaly symetric !
Definition param_elliptic.h:145
Lorene::Param_elliptic::var_F
Scalar var_F
Additive variable change function.
Definition param_elliptic.h:144
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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::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::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
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::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_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::get_etat
int get_etat() const
Returns the logical state.
Definition valeur.h:726
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Valeur::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:131
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139