LORENE
scalar_raccord_zec.C
1/*
2 * Copyright (c) 2003 Eric Gourgoulhon & Jerome Novak
3 *
4 * Copyright (c) 2000-2001 Philippe Grandclement (for preceding Cmp version)
5 *
6 *
7 * This file is part of LORENE.
8 *
9 * LORENE is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * LORENE is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with LORENE; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 *
23 */
24
25
26char 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 $" ;
27
28/*
29 * $Id: scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $
30 * $Log: scalar_raccord_zec.C,v $
31 * Revision 1.8 2014/10/13 08:53:47 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.7 2014/10/06 15:16:16 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.6 2004/06/04 16:14:18 j_novak
38 * In smooth_decay, the configuration space was not up-to-date.
39 *
40 * Revision 1.5 2004/03/05 15:09:59 e_gourgoulhon
41 * Added method smooth_decay.
42 *
43 * Revision 1.4 2003/10/29 11:04:34 e_gourgoulhon
44 * dec2_dpzuis() replaced by dec_dzpuis(2).
45 * inc2_dpzuis() replaced by inc_dzpuis(2).
46 *
47 * Revision 1.3 2003/10/10 15:57:29 j_novak
48 * Added the state one (ETATUN) to the class Scalar
49 *
50 * Revision 1.2 2003/09/25 09:22:33 j_novak
51 * Added a #include<math.h>
52 *
53 * Revision 1.1 2003/09/25 08:58:10 e_gourgoulhon
54 * First version.
55 *
56 *
57 * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_raccord_zec.C,v 1.8 2014/10/13 08:53:47 j_novak Exp $
58 *
59 */
60
61//standard
62#include <cstdlib>
63#include <cmath>
64
65// LORENE
66#include "matrice.h"
67#include "tensor.h"
68#include "proto.h"
69#include "nbr_spx.h"
70#include "utilitaires.h"
71
72// Fait le raccord C1 dans la zec ...
73namespace Lorene {
74// Suppose (pour le moment, le meme nbre de points sur les angles ...)
75// et que la zone precedente est une coquille
76
77void Scalar::raccord_c1_zec(int puis, int nbre, int lmax) {
78
79 assert (nbre>0) ;
80 assert (etat != ETATNONDEF) ;
81 if ((etat == ETATZERO) || (etat == ETATUN))
82 return ;
83
84 // Le mapping doit etre affine :
85 const Map_af* map = dynamic_cast<const Map_af*>(mp) ;
86 if (map == 0x0) {
87 cout << "Le mapping doit etre affine" << endl ;
88 abort() ;
89 }
90
91 int nz = map->get_mg()->get_nzone() ;
92 int nr = map->get_mg()->get_nr (nz-1) ;
93 int nt = map->get_mg()->get_nt (nz-1) ;
94 int np = map->get_mg()->get_np (nz-1) ;
95
96 double alpha = map->get_alpha()[nz-1] ;
97 double r_cont = -1./2./alpha ; //Rayon de debut de la zec.
98
99 // On calcul les coefficients des puissances de 1./r
100 Tbl coef (nbre+2*lmax, nr) ;
101 coef.set_etat_qcq() ;
102
103 int* deg = new int[3] ;
104 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
105 double* auxi = new double[nr] ;
106
107 for (int conte=0 ; conte<nbre+2*lmax ; conte++) {
108 for (int i=0 ; i<nr ; i++)
109 auxi[i] = pow(-1-cos(M_PI*i/(nr-1)), (conte+puis)) ;
110
111 cfrcheb(deg, deg, auxi, deg, auxi) ;
112 for (int i=0 ; i<nr ; i++)
113 coef.set(conte, i) = auxi[i]*pow (alpha, conte+puis) ;
114 }
115
116 delete[] deg ;
117 // Maintenant on va calculer les valeurs de la ieme derivee :
118 Tbl valeurs (nbre, nt, np+1) ;
119 valeurs.set_etat_qcq() ;
120
121 Scalar courant (*this) ;
122 double* res_val = new double[1] ;
123
124 for (int conte=0 ; conte<nbre ; conte++) {
125
126 courant.va.coef() ;
127 courant.va.ylm() ;
128 courant.va.c_cf->t[nz-1]->annule_hard() ;
129
130 int base_r = courant.va.base.get_base_r(nz-2) ;
131 for (int k=0 ; k<np+1 ; k++)
132 for (int j=0 ; j<nt ; j++)
133 if (nullite_plm(j, nt, k, np, courant.va.base) == 1) {
134
135 for (int i=0 ; i<nr ; i++)
136 auxi[i] = (*courant.va.c_cf)(nz-2, k, j, i) ;
137
138 switch (base_r) {
139 case R_CHEB :
140 som_r_cheb (auxi, nr, 1, 1, 1, res_val) ;
141 break ;
142 default :
143 cout << "Cas non prevu dans raccord_zec" << endl ;
144 abort() ;
145 break ;
146 }
147 valeurs.set(conte, k, j) = res_val[0] ;
148 }
149 Scalar copie (courant) ;
150 copie.dec_dzpuis(2) ;
151 courant = copie.dsdr() ;
152 }
153
154 delete [] auxi ;
155 delete [] res_val ;
156
157 // On boucle sur les harmoniques : construction de la matrice
158 // et du second membre
159 va.coef() ;
160 va.ylm() ;
161 va.c_cf->t[nz-1]->annule_hard() ;
163
164 const Base_val& base = va.base ;
165 int base_r, l_quant, m_quant ;
166 for (int k=0 ; k<np+1 ; k++)
167 for (int j=0 ; j<nt ; j++)
168 if (nullite_plm(j, nt, k, np, va.base) == 1) {
169
170 donne_lm (nz, nz-1, j, k, base, m_quant, l_quant, base_r) ;
171
172 if (l_quant<=lmax) {
173
174 Matrice systeme (nbre, nbre) ;
175 systeme.set_etat_qcq() ;
176
177 for (int col=0 ; col<nbre ; col++)
178 for (int lig=0 ; lig<nbre ; lig++) {
179
180 int facteur = (lig%2==0) ? 1 : -1 ;
181 for (int conte=0 ; conte<lig ; conte++)
182 facteur *= puis+col+conte+2*l_quant ;
183 systeme.set(lig, col) = facteur/pow(r_cont, puis+col+lig+2*l_quant) ;
184 }
185
186 systeme.set_band(nbre, nbre) ;
187 systeme.set_lu() ;
188
189 Tbl sec_membre (nbre) ;
190 sec_membre.set_etat_qcq() ;
191 for (int conte=0 ; conte<nbre ; conte++)
192 sec_membre.set(conte) = valeurs(conte, k, j) ;
193
194 Tbl inv (systeme.inverse(sec_membre)) ;
195
196 for (int conte=0 ; conte<nbre ; conte++)
197 for (int i=0 ; i<nr ; i++)
198 va.c_cf->set(nz-1, k, j, i)+=
199 inv(conte)*coef(conte+2*l_quant, i) ;
200 }
201 else for (int i=0 ; i<nr ; i++)
202 va.c_cf->set(nz-1, k, j, i)
203 = 0 ;
204 }
205
206 va.ylm_i() ;
207 set_dzpuis (0) ;
208}
209
210
211//***************************************************************************
212
213void Scalar::smooth_decay(int kk, int nn) {
214
215 assert(kk >= 0) ;
216
217 if (etat == ETATZERO) return ;
218 if (va.get_etat() == ETATZERO) return ;
219
220 const Mg3d& mgrid = *(mp->get_mg()) ;
221
222 int nz = mgrid.get_nzone() ;
223 int nzm1 = nz - 1 ;
224 int np = mgrid.get_np(nzm1) ;
225 int nt = mgrid.get_nt(nzm1) ;
226 int nr_ced = mgrid.get_nr(nzm1) ;
227 int nr_shell = mgrid.get_nr(nzm1-1) ;
228
229 assert(mgrid.get_np(nzm1-1) == np) ;
230 assert(mgrid.get_nt(nzm1-1) == nt) ;
231
232 assert(mgrid.get_type_r(nzm1) == UNSURR) ;
233
234 // In the present version, the mapping must be affine :
235 const Map_af* mapaf = dynamic_cast<const Map_af*>(mp) ;
236 if (mapaf == 0x0) {
237 cout << "Scalar::smooth_decay: present version supports only \n"
238 << " affine mappings !" << endl ;
239 abort() ;
240 }
241
242
243 assert(va.get_etat() == ETATQCQ) ;
244
245
246 // Spherical harmonics are required
247 va.ylm() ;
248
249 // Array for the spectral coefficients in the CED:
250 assert( va.c_cf->t[nzm1] != 0x0) ;
251 Tbl& coefresu = *( va.c_cf->t[nzm1] ) ;
252 coefresu.set_etat_qcq() ;
253
254 // 1-D grid
255 //----------
256 int nbr1[] = {nr_shell, nr_ced} ;
257 int nbt1[] = {1, 1} ;
258 int nbp1[] = {1, 1} ;
259 int typr1[] = {mgrid.get_type_r(nzm1-1), UNSURR} ;
260 Mg3d grid1d(2, nbr1, typr1, nbt1, SYM, nbp1, SYM) ;
261
262 // 1-D mapping
263 //------------
264 double rbound = mapaf->get_alpha()[nzm1-1]
265 + mapaf->get_beta()[nzm1-1] ;
266 double rmin = - mapaf->get_alpha()[nzm1-1]
267 + mapaf->get_beta()[nzm1-1] ;
268 double bound1[] = {rmin, rbound, __infinity} ;
269
270 Map_af map1d(grid1d, bound1) ;
271
272 // 1-D scalars
273 // -----------
274 Scalar uu(map1d) ;
275 uu.std_spectral_base() ;
276 Scalar duu(map1d) ;
277 Scalar vv(map1d) ;
278 Scalar tmp(map1d) ;
279
280 const Base_val& base = va.get_base() ;
281
282 // Loop on the spherical harmonics
283 // -------------------------------
284 for (int k=0; k<=np; k++) {
285 for (int j=0; j<nt; j++) {
286
287 if (nullite_plm(j, nt, k, np, base) != 1) {
288 for (int i=0 ; i<nr_ced ; i++) {
289 coefresu.set(k, j, i) = 0 ;
290 }
291 }
292 else {
293 int base_r_ced, base_r_shell , l_quant, m_quant ;
294 donne_lm(nz, nzm1, j, k, base, m_quant, l_quant, base_r_ced) ;
295 donne_lm(nz, nzm1-1, j, k, base, m_quant, l_quant, base_r_shell) ;
296
297 int nn0 = l_quant + nn ; // lowerst power of decay
298
299 uu.set_etat_qcq() ;
300 uu.va.set_etat_cf_qcq() ;
301 uu.va.c_cf->set_etat_qcq() ;
302 uu.va.c_cf->t[0]->set_etat_qcq() ;
303 uu.va.c_cf->t[1]->set_etat_qcq() ;
304 for (int i=0; i<nr_shell; i++) {
305 uu.va.c_cf->t[0]->set(0, 0, i) =
306 va.c_cf->operator()(nzm1-1, k, j, i) ;
307 uu.va.c_cf->t[0]->set(1, 0, i) = 0. ;
308 uu.va.c_cf->t[0]->set(2, 0, i) = 0. ;
309 }
310
311 uu.va.c_cf->t[1]->set_etat_zero() ;
312
313 uu.va.set_base_r(0, base_r_shell) ;
314 uu.va.set_base_r(1, base_r_ced) ;
315
316 // Computation of the derivatives up to order k at the outer
317 // of the last shell:
318 // ----------------------------------------------------------
319 Tbl derivb(kk+1) ;
320 derivb.set_etat_qcq() ;
321 duu = uu ;
322 derivb.set(0) = uu.val_grid_point(0, 0, 0, nr_shell-1) ;
323 for (int p=1; p<=kk; p++) {
324 tmp = duu.dsdr() ;
325 duu = tmp ;
326 derivb.set(p) = duu.val_grid_point(0, 0, 0, nr_shell-1) ;
327 }
328
329 // Matrix of the linear system to be solved
330 // ----------------------------------------
331
332 Matrice mat(kk+1,kk+1) ;
333 mat.set_etat_qcq() ;
334
335 for (int im=0; im<=kk; im++) {
336
337 double fact = (im%2 == 0) ? 1. : -1. ;
338 fact /= pow(rbound, nn0 + im) ;
339
340 for (int jm=0; jm<=kk; jm++) {
341
342 double prod = 1 ;
343 for (int ir=0; ir<im; ir++) {
344 prod *= nn0 + jm + ir ;
345 }
346
347 mat.set(im, jm) = fact * prod / pow(rbound, jm) ;
348
349 }
350 }
351
352 // Resolution of the linear system to get the coefficients
353 // of the 1/r^i functions
354 //---------------------------------------------------------
355 mat.set_band(kk+1, kk+1) ;
356 mat.set_lu() ;
357 Tbl aa = mat.inverse( derivb ) ;
358
359 // Decaying function
360 // -----------------
361 vv = 0 ;
362 const Coord& r = map1d.r ;
363 for (int p=0; p<=kk; p++) {
364 tmp = aa(p) / pow(r, nn0 + p) ;
365 vv += tmp ;
366 }
367 vv.va.set_base( uu.va.get_base() ) ;
368
369 // The coefficients of the decaying function
370 // are set to the result
371 //-------------------------------------------
372 vv.va.coef() ;
373
374 if (vv.get_etat() == ETATZERO) {
375 for (int i=0; i<nr_ced; i++) {
376 coefresu.set(k,j,i) = 0 ;
377 }
378 }
379 else {
380 assert( vv.va.c_cf != 0x0 ) ;
381 for (int i=0; i<nr_ced; i++) {
382 coefresu.set(k,j,i) = vv.va.c_cf->operator()(1,0,0,i) ;
383 }
384 }
385
386
387 }
388
389
390 }
391 }
392
393 if (va.c != 0x0) {
394 delete va.c ;
395 va.c = 0x0 ;
396 }
397 va.ylm_i() ;
398
399 // Test of the computation
400 // -----------------------
401
402 Scalar tmp1(*this) ;
403 Scalar tmp2(*mp) ;
404 for (int p=0; p<=kk; p++) {
405 double rd = pow(rbound, tmp1.get_dzpuis()) ;
406 double err = 0 ;
407 for (int k=0; k<np; k++) {
408 for (int j=0; j<nt; j++) {
409 double diff = fabs( tmp1.val_grid_point(nzm1, k, j, 0) / rd -
410 tmp1.val_grid_point(nzm1-1, k, j, nr_shell-1) ) ;
411 if (diff > err) err = diff ;
412 }
413 }
414 cout <<
415 "Scalar::smooth_decay: Max error matching of " << p
416 << "-th derivative: " << err << endl ;
417 tmp2 = tmp1.dsdr() ;
418 tmp1 = tmp2 ;
419 }
420
421
422}
423}
Bases of the spectral expansions.
Definition base_val.h:322
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
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Affine radial mapping.
Definition map.h:2027
const double * get_beta() const
Returns the pointer on the array beta.
Definition map_af.C:481
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
Coord r
r coordinate centered on the grid
Definition map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition matrice.C:424
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition matrice.C:364
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition matrice.C:392
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
Definition mtbl_cf.h:205
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl_cf.C:300
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
void smooth_decay(int k, int n)
Performs a matching of the last non-compactified shell with a decaying function where is the spher...
friend Scalar cos(const Scalar &)
Cosine.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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
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
const Scalar & dsdr() const
Returns of *this .
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
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...
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
friend Scalar pow(const Scalar &, int)
Power .
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...
Basic array class.
Definition tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:347
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition valeur.h:480
int get_etat() const
Returns the logical state.
Definition valeur.h:726
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition valeur.C:836
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
#define R_CHEB
base de Chebychev ordinaire (fin)
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64