LORENE
scalar_sol_div.C
1 /*
2  * Resolution of the divergence ODE: df/df + n*f/r = source (source must have dzpuis =2)
3  *
4  * (see file scalar.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char scalar_sol_div_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $" ;
29 
30 /*
31  * $Id: scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $
32  * $Log: scalar_sol_div.C,v $
33  * Revision 1.5 2014/10/13 08:53:47 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.4 2014/10/06 15:16:16 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.3 2005/09/16 14:33:00 j_novak
40  * Added #include <math.h>.
41  *
42  * Revision 1.2 2005/09/16 12:49:52 j_novak
43  * The case with dzpuis=1 is added.
44  *
45  * Revision 1.1 2005/06/08 12:35:22 j_novak
46  * New method for solving divergence-like ODEs.
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Tensor/Scalar/scalar_sol_div.C,v 1.5 2014/10/13 08:53:47 j_novak Exp $
50  *
51  */
52 
53 // C headers
54 #include <cassert>
55 #include <cmath>
56 
57 //Lorene headers
58 #include "tensor.h"
59 #include "diff.h"
60 #include "proto.h"
61 
62 // Local prototypes
63 namespace Lorene {
64 void _sx_r_chebp(Tbl* , int& ) ;
65 void _sx_r_chebi(Tbl* , int& ) ;
66 
67 
68 Scalar Scalar::sol_divergence(int n_factor) const {
69 
70  assert(etat != ETATNONDEF) ;
71  const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
72  assert( mpaff != 0x0) ;
73 
74  Scalar result(*mp) ;
75 
76  if ( etat == ETATZERO )
77  result.set_etat_zero() ;
78  else { //source not zero
79  Base_val base_resu = get_spectral_base() ;
80  base_resu.mult_x() ;
81  const Mg3d* mg = mp->get_mg() ;
82  result.set_etat_qcq() ; result.set_spectral_base(base_resu) ;
83  result.set_spectral_va().set_etat_cf_qcq() ;
84  Valeur sigma(va) ;
85  sigma.ylm_i() ; // work on Fourier basis
86  const Mtbl_cf& source = *sigma.c_cf ;
87 
88  // Checks on the type of domains
89  int nz = mg->get_nzone() ;
90  bool ced = (mg->get_type_r(nz-1) == UNSURR ) ;
91  assert ( (!ced) || (check_dzpuis(2)) || (check_dzpuis(1)) ) ;
92  assert (mg->get_type_r(0) == RARE) ;
93  int nt = mg->get_nt(0) ;
94  int np = mg->get_np(0) ;
95 #ifndef NDEBUG
96  for (int lz = 0; lz<nz; lz++)
97  assert( (mg->get_nt(lz) == nt) && (mg->get_np(lz) == np) ) ;
98 #endif
99  int nr, base_r,l_quant, m_quant;
100  Tbl *so ;
101  Tbl *s_hom ;
102  Tbl *s_part ;
103 
104  // Working objects and initialization
105  Mtbl_cf sol_part(mg, base_resu) ;
106  Mtbl_cf sol_hom(mg, base_resu) ;
107  Mtbl_cf& resu = *result.set_spectral_va().c_cf ;
108  sol_part.annule_hard();
109  sol_hom.annule_hard() ;
110  resu.annule_hard() ;
111 
112  //---------------
113  //-- NUCLEUS ---
114  //---------------
115  int lz = 0 ;
116  nr = mg->get_nr(lz) ;
117 
118  int dege = 1 ; // the operator is degenerate
119  int nr0 = nr - dege ;
120  Tbl vect1(3, 1, nr) ;
121  Tbl vect2(3, 1, nr) ;
122  int base_pipo = 0 ;
123  double alpha = mpaff->get_alpha()[lz] ;
124  double beta = 0. ;
125  Matrice ope_even(nr0, nr0) ; //when the *result* is decomposed on R_CHEBP
126  ope_even.set_etat_qcq() ;
127  for (int i=dege; i<nr; i++) {
128  vect1.annule_hard() ;
129  vect2.annule_hard() ;
130  vect1.set(0,0,i) = 1. ; vect2.set(0,0,i) = 1. ;
131  _dsdx_r_chebp(&vect1, base_pipo) ;
132  _sx_r_chebp(&vect2, base_pipo) ;
133  for (int j=0; j<nr0; j++)
134  ope_even.set(j,i-dege) = (vect1(0,0,j) + n_factor*vect2(0,0,j)) / alpha ;
135  }
136  ope_even.set_lu() ;
137  Matrice ope_odd(nr0, nr0) ; //when the *result* is decomposed on R_CHEBI
138  ope_odd.set_etat_qcq() ;
139  for (int i=0; i<nr0; i++) {
140  vect1.annule_hard() ;
141  vect2.annule_hard() ;
142  vect1.set(0,0,i) = 1. ; vect2.set(0,0,i) = 1. ;
143  _dsdx_r_chebi(&vect1, base_pipo) ;
144  _sx_r_chebi(&vect2, base_pipo) ;
145  for (int j=0; j<nr0; j++)
146  ope_odd.set(j,i) = (vect1(0,0,j) + n_factor*vect2(0,0,j)) / alpha ;
147  }
148  ope_odd.set_lu() ;
149 
150  for (int k=0 ; k<np+1 ; k++)
151  for (int j=0 ; j<nt ; j++) {
152  // to get the spectral base
153  base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
154  assert ( (base_r == R_CHEBP) || (base_r == R_CHEBI) ) ;
155  const Matrice& operateur = (( base_r == R_CHEBP ) ?
156  ope_even : ope_odd ) ;
157  // particular solution
158  so = new Tbl(nr0) ;
159  so->set_etat_qcq() ;
160  for (int i=0 ; i<nr0 ; i++)
161  so->set(i) = source(lz, k, j, i) ;
162 
163  s_part = new Tbl(operateur.inverse(*so)) ;
164 
165  // Putting to Mtbl_cf
166  double somme = 0 ;
167  for (int i=0 ; i<nr0 ; i++) {
168  if (base_r == R_CHEBP) {
169  resu.set(lz, k, j, i+dege) = (*s_part)(i) ;
170  somme += ((i+dege)%2 == 0 ? 1 : -1)*(*s_part)(i) ;
171  }
172  else
173  resu.set(lz,k,j,i) = (*s_part)(i) ;
174  }
175  if (base_r == R_CHEBI)
176  for (int i=nr0; i<nr; i++)
177  resu.set(lz,k,j,i) = 0 ;
178  if (base_r == R_CHEBP) //result must vanish at r=0
179  resu.set(lz, k, j, 0) -= somme ;
180 
181  delete so ;
182  delete s_part ;
183 
184  }
185 
186  //---------------------
187  //-- SHELLS --
188  //---------------------
189  int nz0 = (ced ? nz - 1 : nz) ;
190  for (lz=1 ; lz<nz0 ; lz++) {
191  nr = mg->get_nr(lz) ;
192  alpha = mpaff->get_alpha()[lz] ;
193  beta = mpaff->get_beta()[lz];
194  double ech = beta / alpha ;
195  Diff_id id(R_CHEB, nr) ; const Matrice& mid = id.get_matrice() ;
196  Diff_xdsdx xd(R_CHEB, nr) ; const Matrice& mxd = xd.get_matrice() ;
197  Diff_dsdx dx(R_CHEB, nr) ; const Matrice& mdx = dx.get_matrice() ;
198  Matrice operateur = mxd + ech*mdx + n_factor*mid ;
199  operateur.set_lu() ;
200  // homogeneous solution
201  s_hom = new Tbl(solh(nr, n_factor-1, ech, R_CHEB)) ;
202 
203  for (int k=0 ; k<np+1 ; k++)
204  for (int j=0 ; j<nt ; j++) {
205  // to get the spectral base
206  base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
207  assert (base_r == R_CHEB) ;
208 
209  so = new Tbl(nr) ;
210  so->set_etat_qcq() ;
211  // particular solution
212  Tbl tmp(nr) ;
213  tmp.set_etat_qcq() ;
214  for (int i=0 ; i<nr ; i++)
215  tmp.set(i) = source(lz, k, j, i) ;
216  for (int i=0; i<nr; i++) so->set(i) = beta*tmp(i) ;
217  multx_1d(nr, &tmp.t, R_CHEB) ;
218  for (int i=0; i<nr; i++) so->set(i) += alpha*tmp(i) ;
219 
220  s_part = new Tbl (operateur.inverse(*so)) ;
221 
222  // cleaning things...
223  for (int i=0 ; i<nr ; i++) {
224  sol_part.set(lz, k, j, i) = (*s_part)(i) ;
225  sol_hom.set(lz, k, j, i) = (*s_hom)(1,i) ;
226  }
227 
228  delete so ;
229  delete s_part ;
230  }
231  delete s_hom ;
232  }
233  if (ced) {
234  //---------------
235  //-- CED -----
236  //---------------
237  int dzp = ( check_dzpuis(2) ? 2 : 1) ;
238  nr = source.get_mg()->get_nr(nz-1) ;
239  alpha = mpaff->get_alpha()[nz-1] ;
240  beta = mpaff->get_beta()[nz-1] ;
241  dege = dzp ;
242  nr0 = nr - dege ;
243  Diff_dsdx dx(R_CHEBU, nr) ; const Matrice& mdx = dx.get_matrice() ;
244  Diff_sx sx(R_CHEBU, nr) ; const Matrice& msx = sx.get_matrice() ;
245  Diff_xdsdx xdx(R_CHEBU, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
246  Diff_id id(R_CHEBU, nr) ; const Matrice& mid = id.get_matrice() ;
247  Matrice operateur(nr0, nr0) ;
248  operateur.set_etat_qcq() ;
249  if (dzp == 2)
250  for (int lin=0; lin<nr0; lin++)
251  for (int col=dege; col<nr; col++)
252  operateur.set(lin,col-dege) = (-mdx(lin,col)
253  + n_factor*msx(lin, col)) / alpha ;
254  else {
255  for (int lin=0; lin<nr0; lin++) {
256  for (int col=dege; col<nr; col++)
257  operateur.set(lin,col-dege) = (-mxdx(lin,col)
258  + n_factor*mid(lin, col)) ;
259  }
260  }
261  operateur.set_lu() ;
262  // homogeneous solution
263  s_hom = new Tbl(solh(nr, n_factor-1, 0., R_CHEBU)) ;
264  for (int k=0 ; k<np+1 ; k++)
265  for (int j=0 ; j<nt ; j++) {
266  base_resu.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
267  assert(base_r == R_CHEBU) ;
268 
269  // particular solution
270  so = new Tbl(nr0) ;
271  so->set_etat_qcq() ;
272  for (int i=0 ; i<nr0 ; i++)
273  so->set(i) = source(nz-1, k, j, i) ;
274  s_part = new Tbl(operateur.inverse(*so)) ;
275 
276  // cleaning
277  double somme = 0 ;
278  for (int i=0 ; i<nr0 ; i++) {
279  sol_part.set(nz-1, k, j, i+dege) = (*s_part)(i) ;
280  somme += (*s_part)(i) ;
281  sol_hom.set(nz-1, k, j, i) = (*s_hom)(i) ;
282  }
283  for (int i=nr0; i<nr; i++)
284  sol_hom.set(nz-1, k, j, i) = (*s_hom)(i) ;
285  //result must vanish at infinity
286  sol_part.set(nz-1, k, j, 0) = -somme ;
287  delete so ;
288  delete s_part ;
289  }
290  delete s_hom ;
291  }
292 
293  //-------------------------
294  //-- matching solutions ---
295  //-------------------------
296  if (nz > 1) {
297  Tbl echelles(nz-1) ;
298  echelles.set_etat_qcq() ;
299  for (lz=1; lz<nz; lz++)
300  echelles.set(lz-1)
301  = pow ( (mpaff->get_beta()[lz]/mpaff->get_alpha()[lz] -1),
302  n_factor) ;
303  if (ced) echelles.set(nz-2) = 1./pow(-2., n_factor) ;
304 
305  for (int k=0 ; k<np+1 ; k++)
306  for (int j=0 ; j<nt ; j++) {
307  for (lz=1; lz<nz; lz++) {
308  double val1 = 0 ;
309  double valm1 = 0 ;
310  double valhom1 = 0 ;
311  int nr_prec = mg->get_nr(lz-1) ;
312  nr = mg->get_nr(lz) ;
313  for (int i=0; i<nr_prec; i++)
314  val1 += resu(lz-1, k, j, i) ;
315  for (int i=0; i<nr; i++) {
316  valm1 += ( i%2 == 0 ? 1 : -1)*sol_part(lz, k, j, i) ;
317  valhom1 += ( i%2 == 0 ? 1 : -1)*sol_hom(lz, k, j, i) ;
318  }
319  double lambda = (val1 - valm1) * echelles(lz-1) ;
320  for (int i=0; i<nr; i++)
321  resu.set(lz, k, j, i) = sol_part(lz, k, j, i)
322  + lambda*sol_hom(lz, k, j, i) ;
323 
324  }
325  }
326  }
327  }
328  return result ;
329 }
330 
331 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1294
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_sx.C:100
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
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
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
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_dsdx.C:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:129
Scalar sol_divergence(int n) const
Resolution of a divergence-like equation.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Class for the elementary differential operator Identity (see the base class Diff ).
Definition: diff.h:210
friend Scalar pow(const Scalar &, int)
Power .
Definition: scalar_math.C:451
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Matrix handling.
Definition: matrice.h:152
const double * get_beta() const
Returns the pointer on the array beta.
Definition: map_af.C:481
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
double & set(int j, int i)
Read/write of a particuliar element.
Definition: matrice.h:277
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va.
Definition: scalar.C:797
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 annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
Affine radial mapping.
Definition: map.h:2027
Class for the elementary differential operator (see the base class Diff ).
Definition: diff.h:409
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary)...
Definition: scalar.h:396
void mult_x()
The basis is transformed as with a multiplication by .
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:186
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
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
Class for the elementary differential operator division by (see the base class Diff )...
Definition: diff.h:329
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
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 const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition: diff_xdsdx.C:98
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166