scalar_raccord_zec.C
/*
* Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
*
* Copyright (c) 2000-2001 Philippe Grandclement (for preceding Cmp version)
*
*
* 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 scalar_raccord_zec_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $" ;
/*
* $Id: scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $
* $Log: scalar_raccord_zec.C,v $
* Revision 1.8 2014/10/13 08:53:47 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.7 2014/10/06 15:16:16 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.6 2004/06/04 16:14:18 j_novak
* In smooth_decay, the configuration space was not up-to-date.
*
* Revision 1.5 2004/03/05 15:09:59 e_gourgoulhon
* Added method smooth_decay.
*
* Revision 1.4 2003/10/29 11:04:34 e_gourgoulhon
* dec2_dpzuis() replaced by dec_dzpuis(2).
* inc2_dpzuis() replaced by inc_dzpuis(2).
*
* Revision 1.3 2003/10/10 15:57:29 j_novak
* Added the state one (ETATUN) to the class Scalar
*
* Revision 1.2 2003/09/25 09:22:33 j_novak
* Added a #include<math.h>
*
* Revision 1.1 2003/09/25 08:58:10 e_gourgoulhon
* First version.
*
*
* $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $
*
*/
//standard
#include <cstdlib>
#include <cmath>
// LORENE
#include "matrice.h"
#include "tensor.h"
#include "proto.h"
#include "nbr_spx.h"
#include "utilitaires.h"
// Fait le raccord C1 dans la zec ...
namespace Lorene {
// Suppose (pour le moment, le meme nbre de points sur les angles ...)
// et que la zone precedente est une coquille
void Scalar::raccord_c1_zec(int puis, int nbre, int lmax) {
assert (nbre>0) ;
assert (etat != ETATNONDEF) ;
if ((etat == ETATZERO) || (etat == ETATUN))
return ;
// Le mapping doit etre affine :
const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
if (map == 0x0) {
cout << "Le mapping doit etre affine" << endl ;
abort() ;
}
int nz = map->get_mg()->get_nzone() ;
int nr = map->get_mg()->get_nr (nz-1) ;
int nt = map->get_mg()->get_nt (nz-1) ;
int np = map->get_mg()->get_np (nz-1) ;
double alpha = map->get_alpha()[nz-1] ;
double r_cont = -1./2./alpha ; //Rayon de debut de la zec.
// On calcul les coefficients des puissances de 1./r
Tbl coef (nbre+2*lmax, nr) ;
coef.set_etat_qcq() ;
int* deg = new int[3] ;
deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
double* auxi = new double[nr] ;
for (int conte=0 ; conte<nbre+2*lmax ; conte++) {
for (int i=0 ; i<nr ; i++)
auxi[i] = pow(-1-cos(M_PI*i/(nr-1)), (conte+puis)) ;
cfrcheb(deg, deg, auxi, deg, auxi) ;
for (int i=0 ; i<nr ; i++)
coef.set(conte, i) = auxi[i]*pow (alpha, conte+puis) ;
}
delete[] deg ;
// Maintenant on va calculer les valeurs de la ieme derivee :
Tbl valeurs (nbre, nt, np+1) ;
valeurs.set_etat_qcq() ;
Scalar courant (*this) ;
double* res_val = new double[1] ;
for (int conte=0 ; conte<nbre ; conte++) {
courant.va.coef() ;
courant.va.ylm() ;
courant.va.c_cf->t[nz-1]->annule_hard() ;
int base_r = courant.va.base.get_base_r(nz-2) ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++)
if (nullite_plm(j, nt, k, np, courant.va.base) == 1) {
for (int i=0 ; i<nr ; i++)
auxi[i] = (*courant.va.c_cf)(nz-2, k, j, i) ;
switch (base_r) {
case R_CHEB :
som_r_cheb (auxi, nr, 1, 1, 1, res_val) ;
break ;
default :
cout << "Cas non prevu dans raccord_zec" << endl ;
abort() ;
break ;
}
valeurs.set(conte, k, j) = res_val[0] ;
}
Scalar copie (courant) ;
copie.dec_dzpuis(2) ;
courant = copie.dsdr() ;
}
delete [] auxi ;
delete [] res_val ;
// On boucle sur les harmoniques : construction de la matrice
// et du second membre
va.coef() ;
va.ylm() ;
va.c_cf->t[nz-1]->annule_hard() ;
va.set_etat_cf_qcq() ;
const Base_val& base = va.base ;
int base_r, l_quant, m_quant ;
for (int k=0 ; k<np+1 ; k++)
for (int j=0 ; j<nt ; j++)
if (nullite_plm(j, nt, k, np, va.base) == 1) {
donne_lm (nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
if (l_quant<=lmax) {
Matrice systeme (nbre, nbre) ;
systeme.set_etat_qcq() ;
for (int col=0 ; col<nbre ; col++)
for (int lig=0 ; lig<nbre ; lig++) {
int facteur = (lig%2==0) ? 1 : -1 ;
for (int conte=0 ; conte<lig ; conte++)
facteur *= puis+col+conte+2*l_quant ;
systeme.set(lig, col) = facteur/pow(r_cont, puis+col+lig+2*l_quant) ;
}
systeme.set_band(nbre, nbre) ;
systeme.set_lu() ;
Tbl sec_membre (nbre) ;
sec_membre.set_etat_qcq() ;
for (int conte=0 ; conte<nbre ; conte++)
sec_membre.set(conte) = valeurs(conte, k, j) ;
Tbl inv (systeme.inverse(sec_membre)) ;
for (int conte=0 ; conte<nbre ; conte++)
for (int i=0 ; i<nr ; i++)
va.c_cf->set(nz-1, k, j, i)+=
inv(conte)*coef(conte+2*l_quant, i) ;
}
else for (int i=0 ; i<nr ; i++)
va.c_cf->set(nz-1, k, j, i)
= 0 ;
}
va.ylm_i() ;
set_dzpuis (0) ;
}
//***************************************************************************
void Scalar::smooth_decay(int kk, int nn) {
assert(kk >= 0) ;
if (etat == ETATZERO) return ;
if (va.get_etat() == ETATZERO) return ;
const Mg3d& mgrid = *(mp->get_mg()) ;
int nz = mgrid.get_nzone() ;
int nzm1 = nz - 1 ;
int np = mgrid.get_np(nzm1) ;
int nt = mgrid.get_nt(nzm1) ;
int nr_ced = mgrid.get_nr(nzm1) ;
int nr_shell = mgrid.get_nr(nzm1-1) ;
assert(mgrid.get_np(nzm1-1) == np) ;
assert(mgrid.get_nt(nzm1-1) == nt) ;
assert(mgrid.get_type_r(nzm1) == UNSURR) ;
// In the present version, the mapping must be affine :
const Map_af* mapaf = dynamic_cast<const Map_af*>(mp) ;
if (mapaf == 0x0) {
cout << "Scalar::smooth_decay: present version supports only \n"
<< " affine mappings !" << endl ;
abort() ;
}
assert(va.get_etat() == ETATQCQ) ;
// Spherical harmonics are required
va.ylm() ;
// Array for the spectral coefficients in the CED:
assert( va.c_cf->t[nzm1] != 0x0) ;
Tbl& coefresu = *( va.c_cf->t[nzm1] ) ;
coefresu.set_etat_qcq() ;
// 1-D grid
//----------
int nbr1[] = {nr_shell, nr_ced} ;
int nbt1[] = {1, 1} ;
int nbp1[] = {1, 1} ;
int typr1[] = {mgrid.get_type_r(nzm1-1), UNSURR} ;
Mg3d grid1d(2, nbr1, typr1, nbt1, SYM, nbp1, SYM) ;
// 1-D mapping
//------------
double rbound = mapaf->get_alpha()[nzm1-1]
+ mapaf->get_beta()[nzm1-1] ;
double rmin = - mapaf->get_alpha()[nzm1-1]
+ mapaf->get_beta()[nzm1-1] ;
double bound1[] = {rmin, rbound, __infinity} ;
Map_af map1d(grid1d, bound1) ;
// 1-D scalars
// -----------
Scalar uu(map1d) ;
uu.std_spectral_base() ;
Scalar duu(map1d) ;
Scalar vv(map1d) ;
Scalar tmp(map1d) ;
const Base_val& base = va.get_base() ;
// Loop on the spherical harmonics
// -------------------------------
for (int k=0; k<=np; k++) {
for (int j=0; j<nt; j++) {
if (nullite_plm(j, nt, k, np, base) != 1) {
for (int i=0 ; i<nr_ced ; i++) {
coefresu.set(k, j, i) = 0 ;
}
}
else {
int base_r_ced, base_r_shell , l_quant, m_quant ;
donne_lm(nz, nzm1, j, k, base, m_quant, l_quant, base_r_ced) ;
donne_lm(nz, nzm1-1, j, k, base, m_quant, l_quant, base_r_shell) ;
int nn0 = l_quant + nn ; // lowerst power of decay
uu.set_etat_qcq() ;
uu.va.set_etat_cf_qcq() ;
uu.va.c_cf->set_etat_qcq() ;
uu.va.c_cf->t[0]->set_etat_qcq() ;
uu.va.c_cf->t[1]->set_etat_qcq() ;
for (int i=0; i<nr_shell; i++) {
uu.va.c_cf->t[0]->set(0, 0, i) =
va.c_cf->operator()(nzm1-1, k, j, i) ;
uu.va.c_cf->t[0]->set(1, 0, i) = 0. ;
uu.va.c_cf->t[0]->set(2, 0, i) = 0. ;
}
uu.va.c_cf->t[1]->set_etat_zero() ;
uu.va.set_base_r(0, base_r_shell) ;
uu.va.set_base_r(1, base_r_ced) ;
// Computation of the derivatives up to order k at the outer
// of the last shell:
// ----------------------------------------------------------
Tbl derivb(kk+1) ;
derivb.set_etat_qcq() ;
duu = uu ;
derivb.set(0) = uu.val_grid_point(0, 0, 0, nr_shell-1) ;
for (int p=1; p<=kk; p++) {
tmp = duu.dsdr() ;
duu = tmp ;
derivb.set(p) = duu.val_grid_point(0, 0, 0, nr_shell-1) ;
}
// Matrix of the linear system to be solved
// ----------------------------------------
Matrice mat(kk+1,kk+1) ;
mat.set_etat_qcq() ;
for (int im=0; im<=kk; im++) {
double fact = (im%2 == 0) ? 1. : -1. ;
fact /= pow(rbound, nn0 + im) ;
for (int jm=0; jm<=kk; jm++) {
double prod = 1 ;
for (int ir=0; ir<im; ir++) {
prod *= nn0 + jm + ir ;
}
mat.set(im, jm) = fact * prod / pow(rbound, jm) ;
}
}
// Resolution of the linear system to get the coefficients
// of the 1/r^i functions
//---------------------------------------------------------
mat.set_band(kk+1, kk+1) ;
mat.set_lu() ;
Tbl aa = mat.inverse( derivb ) ;
// Decaying function
// -----------------
vv = 0 ;
const Coord& r = map1d.r ;
for (int p=0; p<=kk; p++) {
tmp = aa(p) / pow(r, nn0 + p) ;
vv += tmp ;
}
vv.va.set_base( uu.va.get_base() ) ;
// The coefficients of the decaying function
// are set to the result
//-------------------------------------------
vv.va.coef() ;
if (vv.get_etat() == ETATZERO) {
for (int i=0; i<nr_ced; i++) {
coefresu.set(k,j,i) = 0 ;
}
}
else {
assert( vv.va.c_cf != 0x0 ) ;
for (int i=0; i<nr_ced; i++) {
coefresu.set(k,j,i) = vv.va.c_cf->operator()(1,0,0,i) ;
}
}
}
}
}
if (va.c != 0x0) {
delete va.c ;
va.c = 0x0 ;
}
va.ylm_i() ;
// Test of the computation
// -----------------------
Scalar tmp1(*this) ;
Scalar tmp2(*mp) ;
for (int p=0; p<=kk; p++) {
double rd = pow(rbound, tmp1.get_dzpuis()) ;
double err = 0 ;
for (int k=0; k<np; k++) {
for (int j=0; j<nt; j++) {
double diff = fabs( tmp1.val_grid_point(nzm1, k, j, 0) / rd -
tmp1.val_grid_point(nzm1-1, k, j, nr_shell-1) ) ;
if (diff > err) err = diff ;
}
}
cout <<
"Scalar::smooth_decay: Max error matching of " << p
<< "-th derivative: " << err << endl ;
tmp2 = tmp1.dsdr() ;
tmp1 = tmp2 ;
}
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Base_val::get_base_r
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition base_val.h:400
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
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_band
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition matrice.C:364
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::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
Lorene::Mtbl_cf::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl_cf.C:300
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
Lorene::Scalar::smooth_decay
void smooth_decay(int k, int n)
Performs a matching of the last non-compactified shell with a decaying function where is the spher...
Definition scalar_raccord_zec.C:213
Lorene::Scalar::cos
friend Scalar cos(const Scalar &)
Cosine.
Definition scalar_math.C:104
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::Scalar
Scalar(const Map &mpi)
Constructor from mapping.
Definition scalar.C:204
Lorene::Scalar::etat
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
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::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Scalar::raccord_c1_zec
void raccord_c1_zec(int puis, int nbre, int lmax)
Performs the matching of the external domain with respect to the last shell using function like wit...
Definition scalar_raccord_zec.C:77
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::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:418
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_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
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::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::get_base
const Base_val & get_base() const
Return the bases for spectral expansions (member base ).
Definition valeur.h:480
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:810
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::set_base_r
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition valeur.C:836
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
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
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:718