scalar_sol_div.C
/*
* Resolution of the divergence ODE: df/df + n*f/r = source (source must have dzpuis =2)
*
* (see file scalar.h for documentation).
*
*/
/*
* Copyright (c) 2005 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 scalar_sol_div_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $" ;
/*
* $Id: scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $
* $Log: scalar_sol_div.C,v $
* Revision 1.5 2014/10/13 08:53:47 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2014/10/06 15:16:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.3 2005/09/16 14:33:00 j_novak
* Added #include <math.h>.
*
* Revision 1.2 2005/09/16 12:49:52 j_novak
* The case with dzpuis=1 is added.
*
* Revision 1.1 2005/06/08 12:35:22 j_novak
* New method for solving divergence-like ODEs.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $
*
*/
// C headers
#include <cassert>
#include <cmath>
//Lorene headers
#include "tensor.h"
#include "diff.h"
#include "proto.h"
// Local prototypes
namespace Lorene {
void _sx_r_chebp(Tbl* , int& ) ;
void _sx_r_chebi(Tbl* , int& ) ;
Scalar Scalar::sol_divergence(int n_factor) const {
assert(etat != ETATNONDEF) ;
const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
assert( mpaff != 0x0) ;
Scalar result(*mp) ;
if ( etat == ETATZERO )
result.set_etat_zero() ;
else { //source not zero
Base_val base_resu = get_spectral_base() ;
base_resu.mult_x() ;
const Mg3d* mg = mp->get_mg() ;
result.set_etat_qcq() ; result.set_spectral_base(base_resu) ;
result.set_spectral_va().set_etat_cf_qcq() ;
Valeur sigma(va) ;
sigma.ylm_i() ; // work on Fourier basis
const Mtbl_cf& source = *sigma.c_cf ;
// Checks on the type of domains
int nz = mg->get_nzone() ;
bool ced = (mg->get_type_r(nz-1) == UNSURR ) ;
assert ( (!ced) || (check_dzpuis(2)) || (check_dzpuis(1)) ) ;
assert (mg->get_type_r(0) == RARE) ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
#ifndef NDEBUG
for (int lz = 0; lz<nz; lz++)
assert( (mg->get_nt(lz) == nt) && (mg->get_np(lz) == np) ) ;
#endif
int nr, base_r,l_quant, m_quant;
Tbl *so ;
Tbl *s_hom ;
Tbl *s_part ;
// Working objects and initialization
Mtbl_cf sol_part(mg, base_resu) ;
Mtbl_cf sol_hom(mg, base_resu) ;
Mtbl_cf& resu = *result.set_spectral_va().c_cf ;
sol_part.annule_hard();
sol_hom.annule_hard() ;
resu.annule_hard() ;
//---------------
//-- NUCLEUS ---
//---------------
int lz = 0 ;
nr = mg->get_nr(lz) ;
int dege = 1 ; // the operator is degenerate
int nr0 = nr - dege ;
Tbl vect1(3, 1, nr) ;
Tbl vect2(3, 1, nr) ;
int base_pipo = 0 ;
double alpha = mpaff->get_alpha()[lz] ;
double beta = 0. ;
Matrice ope_even(nr0, nr0) ; //when the *result* is decomposed on R_CHEBP
ope_even.set_etat_qcq() ;
for (int i=dege; i<nr; i++) {
vect1.annule_hard() ;
vect2.annule_hard() ;
vect1.set(0,0,i) = 1. ; vect2.set(0,0,i) = 1. ;
_dsdx_r_chebp(&vect1, base_pipo) ;
_sx_r_chebp(&vect2, base_pipo) ;
for (int j=0; j<nr0; j++)
ope_even.set(j,i-dege) = (vect1(0,0,j) + n_factor*vect2(0,0,j)) / alpha ;
}
ope_even.set_lu() ;
Matrice ope_odd(nr0, nr0) ; //when the *result* is decomposed on R_CHEBI
ope_odd.set_etat_qcq() ;
for (int i=0; i<nr0; i++) {
vect1.annule_hard() ;
vect2.annule_hard() ;
vect1.set(0,0,i) = 1. ; vect2.set(0,0,i) = 1. ;
_dsdx_r_chebi(&vect1, base_pipo) ;
_sx_r_chebi(&vect2, base_pipo) ;
for (int j=0; j<nr0; j++)
ope_odd.set(j,i) = (vect1(0,0,j) + n_factor*vect2(0,0,j)) / alpha ;
}
ope_odd.set_lu() ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
// to get the spectral base
base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
assert ( (base_r == R_CHEBP) || (base_r == R_CHEBI) ) ;
const Matrice& operateur = (( base_r == R_CHEBP ) ?
ope_even : ope_odd ) ;
// particular solution
so = new Tbl(nr0) ;
so->set_etat_qcq() ;
for (int i=0 ; i<nr0 ; i++)
so->set(i) = source(lz, k, j, i) ;
s_part = new Tbl(operateur.inverse(*so)) ;
// Putting to Mtbl_cf
double somme = 0 ;
for (int i=0 ; i<nr0 ; i++) {
if (base_r == R_CHEBP) {
resu.set(lz, k, j, i+dege) = (*s_part)(i) ;
somme += ((i+dege)%2 == 0 ? 1 : -1)*(*s_part)(i) ;
}
else
resu.set(lz,k,j,i) = (*s_part)(i) ;
}
if (base_r == R_CHEBI)
for (int i=nr0; i<nr; i++)
resu.set(lz,k,j,i) = 0 ;
if (base_r == R_CHEBP) //result must vanish at r=0
resu.set(lz, k, j, 0) -= somme ;
delete so ;
delete s_part ;
}
//---------------------
//-- SHELLS --
//---------------------
int nz0 = (ced ? nz - 1 : nz) ;
for (lz=1 ; lz<nz0 ; lz++) {
nr = mg->get_nr(lz) ;
alpha = mpaff->get_alpha()[lz] ;
beta = mpaff->get_beta()[lz];
double ech = beta / alpha ;
Diff_id id(R_CHEB, nr) ; const Matrice& mid = id.get_matrice() ;
Diff_xdsdx xd(R_CHEB, nr) ; const Matrice& mxd = xd.get_matrice() ;
Diff_dsdx dx(R_CHEB, nr) ; const Matrice& mdx = dx.get_matrice() ;
Matrice operateur = mxd + ech*mdx + n_factor*mid ;
operateur.set_lu() ;
// homogeneous solution
s_hom = new Tbl(solh(nr, n_factor-1, ech, R_CHEB)) ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
// to get the spectral base
base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
assert (base_r == R_CHEB) ;
so = new Tbl(nr) ;
so->set_etat_qcq() ;
// particular solution
Tbl tmp(nr) ;
tmp.set_etat_qcq() ;
for (int i=0 ; i<nr ; i++)
tmp.set(i) = source(lz, k, j, i) ;
for (int i=0; i<nr; i++) so->set(i) = beta*tmp(i) ;
multx_1d(nr, &tmp.t, R_CHEB) ;
for (int i=0; i<nr; i++) so->set(i) += alpha*tmp(i) ;
s_part = new Tbl (operateur.inverse(*so)) ;
// cleaning things...
for (int i=0 ; i<nr ; i++) {
sol_part.set(lz, k, j, i) = (*s_part)(i) ;
sol_hom.set(lz, k, j, i) = (*s_hom)(1,i) ;
}
delete so ;
delete s_part ;
}
delete s_hom ;
}
if (ced) {
//---------------
//-- CED -----
//---------------
int dzp = ( check_dzpuis(2) ? 2 : 1) ;
nr = source.get_mg()->get_nr(nz-1) ;
alpha = mpaff->get_alpha()[nz-1] ;
beta = mpaff->get_beta()[nz-1] ;
dege = dzp ;
nr0 = nr - dege ;
Diff_dsdx dx(R_CHEBU, nr) ; const Matrice& mdx = dx.get_matrice() ;
Diff_sx sx(R_CHEBU, nr) ; const Matrice& msx = sx.get_matrice() ;
Diff_xdsdx xdx(R_CHEBU, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
Diff_id id(R_CHEBU, nr) ; const Matrice& mid = id.get_matrice() ;
Matrice operateur(nr0, nr0) ;
operateur.set_etat_qcq() ;
if (dzp == 2)
for (int lin=0; lin<nr0; lin++)
for (int col=dege; col<nr; col++)
operateur.set(lin,col-dege) = (-mdx(lin,col)
+ n_factor*msx(lin, col)) / alpha ;
else {
for (int lin=0; lin<nr0; lin++) {
for (int col=dege; col<nr; col++)
operateur.set(lin,col-dege) = (-mxdx(lin,col)
+ n_factor*mid(lin, col)) ;
}
}
operateur.set_lu() ;
// homogeneous solution
s_hom = new Tbl(solh(nr, n_factor-1, 0., R_CHEBU)) ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
assert(base_r == R_CHEBU) ;
// particular solution
so = new Tbl(nr0) ;
so->set_etat_qcq() ;
for (int i=0 ; i<nr0 ; i++)
so->set(i) = source(nz-1, k, j, i) ;
s_part = new Tbl(operateur.inverse(*so)) ;
// cleaning
double somme = 0 ;
for (int i=0 ; i<nr0 ; i++) {
sol_part.set(nz-1, k, j, i+dege) = (*s_part)(i) ;
somme += (*s_part)(i) ;
sol_hom.set(nz-1, k, j, i) = (*s_hom)(i) ;
}
for (int i=nr0; i<nr; i++)
sol_hom.set(nz-1, k, j, i) = (*s_hom)(i) ;
//result must vanish at infinity
sol_part.set(nz-1, k, j, 0) = -somme ;
delete so ;
delete s_part ;
}
delete s_hom ;
}
//-------------------------
//-- matching solutions ---
//-------------------------
if (nz > 1) {
Tbl echelles(nz-1) ;
echelles.set_etat_qcq() ;
for (lz=1; lz<nz; lz++)
echelles.set(lz-1)
= pow ( (mpaff->get_beta()[lz]/mpaff->get_alpha()[lz] -1),
n_factor) ;
if (ced) echelles.set(nz-2) = 1./pow(-2., n_factor) ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++) {
for (lz=1; lz<nz; lz++) {
double val1 = 0 ;
double valm1 = 0 ;
double valhom1 = 0 ;
int nr_prec = mg->get_nr(lz-1) ;
nr = mg->get_nr(lz) ;
for (int i=0; i<nr_prec; i++)
val1 += resu(lz-1, k, j, i) ;
for (int i=0; i<nr; i++) {
valm1 += ( i%2 == 0 ? 1 : -1)*sol_part(lz, k, j, i) ;
valhom1 += ( i%2 == 0 ? 1 : -1)*sol_hom(lz, k, j, i) ;
}
double lambda = (val1 - valm1) * echelles(lz-1) ;
for (int i=0; i<nr; i++)
resu.set(lz, k, j, i) = sol_part(lz, k, j, i)
+ lambda*sol_hom(lz, k, j, i) ;
}
}
}
}
return result ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::mult_x
void mult_x()
The basis is transformed as with a multiplication by .
Definition base_val_manip.C:157
Lorene::Base_val::give_quant_numbers
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Definition base_val_quantum.C:81
Lorene::Diff_dsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Lorene::Diff_dsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_dsdx.C:94
Lorene::Diff_id
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Lorene::Diff_sx
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
Lorene::Diff_sx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sx.C:100
Lorene::Diff_xdsdx
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Lorene::Diff_xdsdx::get_matrice
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_xdsdx.C:98
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Map_af::get_beta
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
Lorene::Map_af::get_alpha
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Matrice::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:175
Lorene::Matrice::set
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Lorene::Matrice::inverse
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:424
Lorene::Matrice::set_lu
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
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_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::Mg3d::get_type_r
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
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::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:453
Lorene::Mtbl_cf::annule_hard
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
Lorene::Scalar::sol_divergence
Scalar sol_divergence(int n) const
Resolution of a divergence-like equation.
Definition scalar_sol_div.C:68
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::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:204
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
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::va
Valeur va
The numerical value of the Scalar.
Definition scalar.h:405
Lorene::Scalar::pow
friend Scalar pow(const Scalar &, int)
Power .
Definition scalar_math.C:451
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::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::set_etat_cf_qcq
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
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
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition type_parite.h:180
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene
Lorene prototypes.
Definition app_hor.h:64