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 }
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
const double * get_alpha() const
Returns the pointer on the array alpha.
Definition: map_af.C:477
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Definition: matrice.C:392
Lorene prototypes.
Definition: app_hor.h:64
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition: mtbl_cf.h:294
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Definition: matrice.C:424
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
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
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int get_etat() const
Returns the logical state.
Definition: valeur.h:726
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:451
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.
Definition: scalar_math.C:104
Matrix handling.
Definition: matrice.h:152
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Valeur va
The numerical value of the Scalar.
Definition: scalar.h:405
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl_cf.C:300
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
void set_band(int up, int low) const
Calculate the band storage of *std.
Definition: matrice.C:364
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
Bases of the spectral expansions.
Definition: base_val.h:322
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition: valeur.C:836
Affine radial mapping.
Definition: map.h:2027
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...
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Definition: scalar.h:396
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: matrice.C:175
Basic array class.
Definition: tbl.h:161
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_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
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...
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s which contain the spectral coefficients in each domain...
Definition: mtbl_cf.h:205
Coord r
r coordinate centered on the grid
Definition: map.h:718
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:480
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166