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 
26 char 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  */
58 namespace Lorene {
59 Scalar 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
Definition: type_parite.h:180
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene prototypes.
Definition: app_hor.h:64