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
28char 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
63namespace Lorene {
64void _sx_r_chebp(Tbl* , int& ) ;
65void _sx_r_chebi(Tbl* , int& ) ;
66
67
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
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
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() ;
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() ;
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}
Bases of the spectral expansions.
Definition base_val.h:322
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sx.C:100
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
Matrix handling.
Definition matrice.h:152
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
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Scalar sol_divergence(int n) const
Resolution of a divergence-like equation.
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
int etat
The logical state ETATNONDEF (undefined), ETATZERO (null), ETATUN (one), or ETATQCQ (ordinary).
Definition scalar.h:396
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
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
Values and coefficients of a (real-value) function.
Definition valeur.h:287
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
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