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 
26 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 $" ;
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 ...
73 namespace Lorene {
74 // Suppose (pour le moment, le meme nbre de points sur les angles ...)
75 // et que la zone precedente est une coquille
76 
77 void 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() ;
162  va.set_etat_cf_qcq() ;
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 
213 void 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 ** 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
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
friend Scalar cos(const Scalar &)
Cosine.
Definition: scalar_math.C:104
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...
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 .
Definition: scalar_deriv.C:113
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...
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:451
Valeur va
The numerical value of the Scalar
Definition: scalar.h:405
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
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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
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.
Definition: valeur_coef.C:148
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()
Definition: valeur_ylm_i.C:131
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
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