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 }
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
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
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
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
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition map_af.C:477
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
Multi-domain grid.
Definition grilles.h:273
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
friend Scalar cos(const Scalar &)
Cosine.
void smooth_decay(int k, int n)
Performs a matching of the last non-compactified shell with a decaying function where is the spher...
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
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...
friend Scalar pow(const Scalar &, int)
Power .
Valeur va
The numerical value of the Scalar
Definition scalar.h:405
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_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 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 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