LORENE
pois_vect_r0.C
1/*
2 * Solution of the l=0 part of the vector Poisson equation (only r-component)
3 *
4 */
5
6/*
7 * Copyright (c) 2007 Jerome Novak
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26char pois_vect_r0_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pois_vect_r0.C,v 1.2 2014/10/13 08:53:29 j_novak Exp $" ;
27
28/*
29 * $Id: pois_vect_r0.C,v 1.2 2014/10/13 08:53:29 j_novak Exp $
30 * $Log: pois_vect_r0.C,v $
31 * Revision 1.2 2014/10/13 08:53:29 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.1 2007/01/23 17:08:46 j_novak
35 * New function pois_vect_r0.C to solve the l=0 part of the vector Poisson
36 * equation, which involves only the r-component.
37 *
38 *
39 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/pois_vect_r0.C,v 1.2 2014/10/13 08:53:29 j_novak Exp $
40 *
41 */
42
43// Lorene headers
44#include "metric.h"
45#include "proto.h"
46#include "diff.h"
47
48/*
49 * This function solves for the l=0 component of
50 *
51 * d2 f 2 df 2f
52 * ---- + - -- - -- = source
53 * dr2 r dr r2
54 *
55 * and returns the soluton f.
56 * The input Scalar must have dzpuis = 4.
57 */
58namespace Lorene {
59Scalar pois_vect_r0(const Scalar& source) {
60
61 const Map& map0 = source.get_mp() ;
62 const Map_af* map1 = dynamic_cast<const Map_af*>(&map0) ;
63 assert(map1 != 0x0) ;
64 const Map_af& map = *map1 ;
65
66 const Mg3d& mgrid = *map.get_mg() ;
67 int nz = mgrid.get_nzone() ;
68
69 Scalar resu(map) ;
70 if (source.get_etat() == ETATZERO) {
71 resu = 0 ;
72 return resu ;
73 }
74
75 resu.annule_hard() ;
76 resu.std_spectral_base_odd() ;
77 resu.set_spectral_va().ylm() ;
78 Mtbl_cf& sol_coef = (*resu.set_spectral_va().c_cf) ;
79 const Base_val& base = source.get_spectral_base() ;
80 assert(resu.get_spectral_base() == base) ;
81 assert(source.check_dzpuis(4)) ;
82
83 Mtbl_cf sol_part(mgrid, base) ; sol_part.annule_hard() ;
84 Mtbl_cf sol_hom1(mgrid, base) ; sol_hom1.annule_hard() ;
85 Mtbl_cf sol_hom2(mgrid, base) ; sol_hom2.annule_hard() ;
86
87 { int lz = 0 ;
88 int nr = mgrid.get_nr(lz) ;
89 double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
90 assert(mgrid.get_type_r(lz) == RARE) ;
91 int base_r = R_CHEBI ;
92 Matrice ope(nr,nr) ;
93 ope.annule_hard() ;
94 Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
95 Diff_sxdsdx sdx(base_r, nr) ; const Matrice& msdx = sdx.get_matrice() ;
96 Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
97
98 for (int lin=0; lin<nr-1; lin++)
99 for (int col=0; col<nr; col++)
100 ope.set(lin,col) = (mdx2(lin,col) + 2*msdx(lin,col) - 2*ms2(lin,col))/alpha2 ;
101
102 ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
103 for (int i=1; i<nr; i++)
104 ope.set(nr-1, i) = 0 ;
105
106 Tbl rhs(nr) ;
107 rhs.annule_hard() ;
108 for (int i=0; i<nr-1; i++)
109 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
110 rhs.set(nr-1) = 0 ;
111 ope.set_lu() ;
112 Tbl sol = ope.inverse(rhs) ;
113
114 for (int i=0; i<nr; i++)
115 sol_part.set(lz, 0, 0, i) = sol(i) ;
116
117 rhs.annule_hard() ;
118 rhs.set(nr-1) = 1 ;
119 sol = ope.inverse(rhs) ;
120 for (int i=0; i<nr; i++)
121 sol_hom1.set(lz, 0, 0, i) = sol(i) ;
122
123 }
124
125 for (int lz=1; lz<nz-1; lz++) {
126 int nr = mgrid.get_nr(lz) ;
127 double alpha = map.get_alpha()[lz] ;
128 double beta = map.get_beta()[lz] ;
129 double ech = beta / alpha ;
130 assert(mgrid.get_type_r(lz) == FIN) ;
131 int base_r = R_CHEB ;
132
133 Matrice ope(nr,nr) ;
134 ope.annule_hard() ;
135
136 Diff_dsdx dx(base_r, nr) ; const Matrice& mdx = dx.get_matrice() ;
137 Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
138 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
139 Diff_xdsdx xdx(base_r, nr) ; const Matrice& mxdx = xdx.get_matrice() ;
140 Diff_xdsdx2 xdx2(base_r, nr) ; const Matrice& mxdx2 = xdx2.get_matrice() ;
141 Diff_x2dsdx2 x2dx2(base_r, nr) ; const Matrice& mx2dx2 = x2dx2.get_matrice() ;
142
143 for (int lin=0; lin<nr-2; lin++)
144 for (int col=0; col<nr; col++)
145 ope.set(lin, col) = mx2dx2(lin, col) + 2*ech*mxdx2(lin, col) + ech*ech*mdx2(lin, col)
146 + 2*(mxdx(lin, col) + ech*mdx(lin, col)) - 2*mid(lin, col) ;
147
148 ope.set(nr-2, 0) = 0 ;
149 ope.set(nr-2, 1) = 1 ;
150 for (int col=2; col<nr; col++) {
151 ope.set(nr-2, col) = 0 ;
152 }
153 ope.set(nr-1, 0) = 1 ;
154 for (int col=1; col<nr; col++)
155 ope.set(nr-1, col) = 0 ;
156
157 Tbl src(nr) ;
158 src.set_etat_qcq() ;
159 for (int i=0; i<nr; i++)
160 src.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
161 Tbl xsrc = src ; multx_1d(nr, &xsrc.t, base_r) ;
162 Tbl x2src = src ; multx2_1d(nr, &x2src.t, base_r) ;
163 Tbl rhs(nr) ;
164 rhs.set_etat_qcq() ;
165 for (int i=0; i<nr-2; i++)
166 rhs.set(i) = beta*beta*src(i) + 2*beta*alpha*xsrc(i) + alpha*alpha*x2src(i) ;
167 rhs.set(nr-2) = 0 ;
168 rhs.set(nr-1) = 0 ;
169 ope.set_lu() ;
170 Tbl sol = ope.inverse(rhs) ;
171
172 for (int i=0; i<nr; i++)
173 sol_part.set(lz, 0, 0, i) = sol(i) ;
174
175 rhs.annule_hard() ;
176 rhs.set(nr-2) = 1 ;
177 sol = ope.inverse(rhs) ;
178 for (int i=0; i<nr; i++)
179 sol_hom1.set(lz, 0, 0, i) = sol(i) ;
180
181 rhs.set(nr-2) = 0 ;
182 rhs.set(nr-1) = 1 ;
183 sol = ope.inverse(rhs) ;
184 for (int i=0; i<nr; i++)
185 sol_hom2.set(lz, 0, 0, i) = sol(i) ;
186
187 }
188
189 { int lz = nz-1 ;
190 int nr = mgrid.get_nr(lz) ;
191 double alpha2 = map.get_alpha()[lz]*map.get_alpha()[lz] ;
192 assert(mgrid.get_type_r(lz) == UNSURR) ;
193 int base_r = R_CHEBU ;
194
195 Matrice ope(nr,nr) ;
196 ope.annule_hard() ;
197 Diff_dsdx2 dx2(base_r, nr) ; const Matrice& mdx2 = dx2.get_matrice() ;
198 Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
199
200 for (int lin=0; lin<nr-3; lin++)
201 for (int col=0; col<nr; col++)
202 ope.set(lin, col) = (mdx2(lin, col) - 2*ms2(lin, col))/alpha2 ;
203
204 for (int i=0; i<nr; i++) {
205 ope.set(nr-3, i) = i*i ; //for the finite part (derivative = 0 at infty)
206 }
207
208 for (int i=0; i<nr; i++) {
209 ope.set(nr-2, i) = 1 ; //for the limit at inifinity
210 }
211
212 ope.set(nr-1, 0) = 1 ; //for the true homogeneous solution
213 for (int i=1; i<nr; i++)
214 ope.set(nr-1, i) = 0 ;
215
216 Tbl rhs(nr) ;
217 rhs.annule_hard() ;
218 for (int i=0; i<nr-3; i++)
219 rhs.set(i) = (*source.get_spectral_va().c_cf)(lz, 0, 0, i) ;
220 rhs.set(nr-3) = 0 ;
221 rhs.set(nr-2) = 0 ;
222 rhs.set(nr-1) = 0 ;
223 ope.set_lu() ;
224 Tbl sol = ope.inverse(rhs) ;
225 for (int i=0; i<nr; i++)
226 sol_part.set(lz, 0, 0, i) = sol(i) ;
227
228 rhs.annule_hard() ;
229 rhs.set(nr-1) = 1 ;
230 sol = ope.inverse(rhs) ;
231 for (int i=0; i<nr; i++)
232 sol_hom2.set(lz, 0, 0, i) = sol(i) ;
233
234 }
235
236 Mtbl_cf dpart = sol_part ; dpart.dsdx() ;
237 Mtbl_cf dhom1 = sol_hom1 ; dhom1.dsdx() ;
238 Mtbl_cf dhom2 = sol_hom2 ; dhom2.dsdx() ;
239
240 Matrice systeme(2*(nz-1), 2*(nz-1)) ;
241 systeme.annule_hard() ;
242 Tbl rhs(2*(nz-1)) ;
243 rhs.annule_hard() ;
244
245 //Nucleus
246 int lin = 0 ;
247 int col = 0 ;
248 double alpha = map.get_alpha()[0] ;
249 systeme.set(lin,col) = sol_hom1.val_out_bound_jk(0, 0, 0) ;
250 rhs.set(lin) -= sol_part.val_out_bound_jk(0, 0, 0) ;
251 lin++ ;
252 systeme.set(lin,col) = dhom1.val_out_bound_jk(0, 0, 0) / alpha ;
253 rhs.set(lin) -= dpart.val_out_bound_jk(0, 0, 0) / alpha ;
254 col++ ;
255
256 //Shells
257 for (int lz=1; lz<nz-1; lz++) {
258 alpha = map.get_alpha()[lz] ;
259 lin-- ;
260 systeme.set(lin,col) -= sol_hom1.val_in_bound_jk(lz, 0, 0) ;
261 systeme.set(lin,col+1) -= sol_hom2.val_in_bound_jk(lz, 0, 0) ;
262 rhs.set(lin) += sol_part.val_in_bound_jk(lz, 0, 0) ;
263
264 lin++ ;
265 systeme.set(lin,col) -= dhom1.val_in_bound_jk(lz, 0, 0) / alpha ;
266 systeme.set(lin,col+1) -= dhom2.val_in_bound_jk(lz, 0, 0) / alpha ;
267 rhs.set(lin) += dpart.val_in_bound_jk(lz, 0, 0) / alpha;
268
269 lin++ ;
270 systeme.set(lin, col) += sol_hom1.val_out_bound_jk(lz, 0, 0) ;
271 systeme.set(lin, col+1) += sol_hom2.val_out_bound_jk(lz, 0, 0) ;
272 rhs.set(lin) -= sol_part.val_out_bound_jk(lz, 0, 0) ;
273
274 lin++ ;
275 systeme.set(lin, col) += dhom1.val_out_bound_jk(lz, 0, 0) / alpha ;
276 systeme.set(lin, col+1) += dhom2.val_out_bound_jk(lz, 0, 0) / alpha ;
277 rhs.set(lin) -= dpart.val_out_bound_jk(lz, 0, 0) / alpha ;
278 col += 2 ;
279 }
280
281 //CED
282 alpha = map.get_alpha()[nz-1] ;
283 lin-- ;
284 systeme.set(lin,col) -= sol_hom2.val_in_bound_jk(nz-1, 0, 0) ;
285 rhs.set(lin) += sol_part.val_in_bound_jk(nz-1, 0, 0) ;
286
287 lin++ ;
288 systeme.set(lin,col) -= (-4*alpha)*dhom2.val_in_bound_jk(nz-1, 0, 0) ;
289 rhs.set(lin) += (-4*alpha)*dpart.val_in_bound_jk(nz-1, 0, 0) ;
290
291 systeme.set_lu() ;
292 Tbl coef = systeme.inverse(rhs) ;
293 int indice = 0 ;
294
295 for (int i=0; i<mgrid.get_nr(0); i++)
296 sol_coef.set(0, 0, 0, i) = sol_part(0, 0, 0, i)
297 + coef(indice)*sol_hom1(0, 0, 0, i) ;
298 indice++ ;
299 for (int lz=1; lz<nz-1; lz++) {
300 for (int i=0; i<mgrid.get_nr(lz); i++)
301 sol_coef.set(lz, 0, 0, i) = sol_part(lz, 0, 0, i)
302 +coef(indice)*sol_hom1(lz, 0, 0, i)
303 +coef(indice+1)*sol_hom2(lz, 0, 0, i) ;
304 indice += 2 ;
305 }
306 for (int i=0; i<mgrid.get_nr(nz-1); i++)
307 sol_coef.set(nz-1, 0, 0, i) = sol_part(nz-1, 0, 0, i)
308 +coef(indice)*sol_hom2(nz-1, 0, 0, i) ;
309
310 delete resu.set_spectral_va().c ;
311 resu.set_spectral_va().c = 0x0 ;
312 resu.set_dzpuis(0) ;
313 resu.set_spectral_va().ylm_i() ;
314
315 return resu ;
316}
317}
#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)
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64