LORENE
vector_poisson_boundary2.C
1/*
2 * Method for vector Poisson equation inverting eqs. for V^r and eta as a block.
3 *
4 * (see file vector.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 vector_poisson_boundary2_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $" ;
29
30/*
31 * $Id: vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $
32 * $Log: vector_poisson_boundary2.C,v $
33 * Revision 1.9 2014/10/13 08:53:45 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.8 2014/10/06 15:13:21 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.7 2008/08/20 15:07:36 n_vasset
40 * Cleaning up the code...
41 *
42 * Revision 1.5 2007/09/05 12:35:18 j_novak
43 * Homogeneous solutions are no longer obtained through the analytic formula, but
44 * by solving (again) the operator system with almost zero as r.h.s. This is to
45 * lower the condition number of the matching system.
46 *
47 * Revision 1.4 2007/01/23 17:08:46 j_novak
48 * New function pois_vect_r0.C to solve the l=0 part of the vector Poisson
49 * equation, which involves only the r-component.
50 *
51 * Revision 1.3 2005/12/30 13:39:38 j_novak
52 * Changed the Galerkin base in the nucleus to (hopefully) stabilise the solver
53 * when used in an iteration. Similar changes in the CED too.
54 *
55 * Revision 1.2 2005/02/15 15:43:18 j_novak
56 * First version of the block inversion for the vector Poisson equation (method 6).
57 *
58 *
59 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary2.C,v 1.9 2014/10/13 08:53:45 j_novak Exp $
60 *
61 */
62
63// C headers
64#include <cassert>
65#include <cstdlib>
66#include <cmath>
67
68// Lorene headers
69#include "metric.h"
70#include "diff.h"
71#include "param_elliptic.h"
72#include "proto.h"
73#include "graphique.h"
74#include "map.h"
75#include "utilitaires.h"
76
77
78
79namespace Lorene {
80void Vector::poisson_boundary2(double lam, Vector& resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const {
81
82 const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
83#ifndef NDEBUG
84 for (int i=0; i<3; i++)
85 assert(cmp[i]->check_dzpuis(4)) ;
86 // All this has a meaning only for spherical components:
87 const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
88 assert(bvs != 0x0) ;
89 //## ... and affine mapping, for the moment!
90 assert( mpaff != 0x0) ;
91#endif
92
93 if (fabs(lam + 1.) < 1.e-8) {
94
95 cout <<" lambda = -1 is not supported by this code!" << endl;
96 abort();
97 }
98
99
100
101 // Extraction of Mtbl_cf objects from boundary informations;
103 bound_vr2.set_spectral_va().ylm();
105 bound_eta2.set_spectral_va().ylm();
106 Scalar bound_mu2 = boundmu; bound_mu2.set_spectral_va().ylm();
107
108 const Mtbl_cf *bound_vr = bound_vr2.get_spectral_va().c_cf;
109 const Mtbl_cf *bound_mu = bound_mu2.get_spectral_va().c_cf;
110 const Mtbl_cf *bound_eta = bound_eta2.get_spectral_va().c_cf;
111
112
113
114
115 // Some working objects
116 //---------------------
117 const Mg3d& mg = *mpaff->get_mg() ;
118 int nz = mg.get_nzone() ; int nzm1 = nz - 1;
119 assert(mg.get_type_r(0) == RARE) ;
120 assert(mg.get_type_r(nzm1) == UNSURR) ;
121 Scalar S_r = *cmp[0] ;
122 Scalar S_eta = eta() ;
123 Scalar het(*mpaff) ;
124 Scalar vr(*mpaff) ;
125
126 if (S_r.get_etat() == ETATZERO) {
127 if (S_eta.get_etat() == ETATZERO) {
128
129 S_r.annule_hard() ;
130 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
131 S_eta.annule_hard() ;
132 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
133 }
134 else {
135 S_r.annule_hard() ;
136 S_r.set_spectral_base(S_eta.get_spectral_base()) ;
137 }
138 }
139 if (S_eta.get_etat() == ETATZERO) {
140 S_eta.annule_hard() ;
141 S_eta.set_spectral_base(S_r.get_spectral_base()) ;
142 }
143
144 S_r.set_spectral_va().ylm() ;
145 S_eta.set_spectral_va().ylm() ;
146 const Base_val& base = S_eta.get_spectral_va().base ;
147 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
148 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
149 Mtbl_cf sol_hom_un_eta(mg, base) ; sol_hom_un_eta.annule_hard() ;
150 Mtbl_cf sol_hom_deux_eta(mg, base) ; sol_hom_deux_eta.annule_hard() ;
151 Mtbl_cf sol_hom_trois_eta(mg, base) ; sol_hom_trois_eta.annule_hard() ;
152 Mtbl_cf sol_hom_quatre_eta(mg, base) ; sol_hom_quatre_eta.annule_hard() ;
153 Mtbl_cf sol_hom_un_vr(mg, base) ; sol_hom_un_vr.annule_hard() ;
154 Mtbl_cf sol_hom_deux_vr(mg, base) ; sol_hom_deux_vr.annule_hard() ;
155 Mtbl_cf sol_hom_trois_vr(mg, base) ; sol_hom_trois_vr.annule_hard() ;
156 Mtbl_cf sol_hom_quatre_vr(mg, base) ; sol_hom_quatre_vr.annule_hard() ;
157
158 // Solution of the l=0 part for Vr
159 //---------------------------------
160 Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
161 sou_l0.set_spectral_va().ylm() ;
162
164 for (int l=0; l<nz; l++){
165 param_l0.set_poisson_vect_r(l, true) ;
166 }
167
168
169 Scalar vrl0 = sou_l0.sol_elliptic_boundary(param_l0, *bound_vr, dir_vr, neum_vr) ;
170
171
172 // Build-up & inversion of the system for (eta, V^r) in each domain (except for the nucleus!)
173 //-----------------------------------------------------------------
174
175
176 // Shells
177 //-------
178
179 int nr = mg.get_nr(1) ;
180 int nt = mg.get_nt(1) ;
181 int np = mg.get_np(1) ;
182 double alpha = mpaff->get_alpha()[1] ; double alp2 = alpha*alpha ;
183 double beta = mpaff->get_beta()[1] ;
184
185 int l_q = 0 ; int m_q = 0; int base_r = 0 ;
186 for (int zone=1 ; zone<nzm1 ; zone++) {
187 nr = mg.get_nr(zone) ;
188 assert (nr > 5) ;
189 assert(nt == mg.get_nt(zone)) ;
190 assert(np == mg.get_np(zone)) ;
191 alpha = mpaff->get_alpha()[zone] ;
192 beta = mpaff->get_beta()[zone] ;
193 double ech = beta / alpha ;
194
195 // Loop on l and m
196 //----------------
197 for (int k=0 ; k<np+1 ; k++) {
198 for (int j=0 ; j<nt ; j++) {
199 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
200 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
201 int dege1 = 2 ; //degeneracy of eq.1
202 int dege2 = 2 ; //degeneracy of eq.2
203 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
204 int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
205 int nrtot = 2*nr ;
206 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
207 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
208 Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
209 Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
210 Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
211 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
212 Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
213 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
214
215 // Building the operator
216 //----------------------
217 for (int lin=0; lin<nr_eq1; lin++) {
218 for (int col=0; col<nr; col++)
219 oper.set(lin,col)
220 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
221 + 2*(mxd(lin,col) + ech*md(lin,col))
222 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
223 for (int col=0; col<nr; col++)
224 oper.set(lin,col+nr)
225 = lam*(mxd(lin,col) + ech*md(lin,col))
226 + 2*(1+lam)*mid(lin,col) ;
227 }
228 oper.set(nr_eq1, 0) = 1 ;
229 for (int col=1; col<2*nr; col++)
230 oper.set(nr_eq1, col) = 0 ;
231 oper.set(nr_eq1+1, 0) = 0 ;
232 oper.set(nr_eq1+1, 1) = 1 ;
233 for (int col=2; col<2*nr; col++)
234 oper.set(nr_eq1+1, col) = 0 ;
235 for (int lin=0; lin<nr_eq2; lin++) {
236 for (int col=0; col<nr; col++)
237 oper.set(lin+nr,col)
238 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
239 - (lam+2)*mid(lin,col)) ;
240 for (int col=0; col<nr; col++)
241 oper.set(lin+nr, col+nr)
242 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
243 + ech*ech*md2(lin,col)
244 + 2*(mxd(lin,col) + ech*md(lin,col)))
245 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
246 }
247 for (int col=0; col<nr; col++)
248 oper.set(nr+nr_eq2, col) = 0 ;
249 oper.set(nr+nr_eq2, nr) = 1 ;
250 for (int col=nr+1; col<2*nr; col++)
251 oper.set(nr+nr_eq2, col) = 0 ;
252 for (int col=0; col<nr+1; col++)
253 oper.set(nr+nr_eq2+1, col) = 0 ;
254 oper.set(nr+nr_eq2+1, nr+1) = 1 ;
255 for (int col=nr+2; col<2*nr; col++)
256 oper.set(nr+nr_eq2+1, col) = 0 ;
257
258 oper.set_lu() ;
259
260 // Filling the r.h.s
261 //------------------
262 Tbl sr(nr) ; sr.set_etat_qcq() ;
263 Tbl seta(nr) ; seta.set_etat_qcq() ;
264 for (int i=0; i<nr; i++) {
265 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
266 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
267 }
268 Tbl xsr= sr ; Tbl x2sr= sr ;
269 Tbl xseta= seta ; Tbl x2seta = seta ;
270 multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
271 multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
272
273 for (int i=0; i<nr_eq1; i++)
274 sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
275 + beta*beta*seta(i);
276 sec_membre.set(nr_eq1) = 0 ;
277 sec_membre.set(nr_eq1+1) = 0 ;
278 for (int i=0; i<nr_eq2; i++)
279 sec_membre.set(i+nr) = beta*beta*sr(i)
280 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
281 sec_membre.set(nr+nr_eq2) = 0 ;
282 sec_membre.set(nr+nr_eq2+1) = 0 ;
283
284 // Inversion of the "big" operator
285 //--------------------------------
286 Tbl big_res = oper.inverse(sec_membre) ;
287 for (int i=0; i<nr; i++) {
288 sol_part_eta.set(zone, k, j, i) = big_res(i) ;
289 sol_part_vr.set(zone, k, j, i) = big_res(nr+i) ;
290 }
291
292 // Getting the homogeneous solutions
293 //----------------------------------
294 sec_membre.annule_hard() ;
295 sec_membre.set(nr_eq1) = 1 ;
296 big_res = oper.inverse(sec_membre) ;
297 for (int i=0 ; i<nr ; i++) {
298 sol_hom_un_eta.set(zone, k, j, i) = big_res(i) ;
299 sol_hom_un_vr.set(zone, k, j, i) = big_res(nr+i) ;
300 }
301 sec_membre.set(nr_eq1) = 0 ;
302 sec_membre.set(nr_eq1+1) = 1 ;
303 big_res = oper.inverse(sec_membre) ;
304 for (int i=0 ; i<nr ; i++) {
305 sol_hom_deux_eta.set(zone, k, j, i) = big_res(i) ;
306 sol_hom_deux_vr.set(zone, k, j, i) = big_res(nr+i) ;
307 }
308 sec_membre.set(nr_eq1+1) = 0 ;
309 sec_membre.set(nr+nr_eq2) = 1 ;
310 big_res = oper.inverse(sec_membre) ;
311 for (int i=0 ; i<nr ; i++) {
312 sol_hom_trois_eta.set(zone, k, j, i) = big_res(i) ;
313 sol_hom_trois_vr.set(zone, k, j, i) = big_res(nr+i) ;
314 }
315 sec_membre.set(nr+nr_eq2) = 0 ;
316 sec_membre.set(nr+nr_eq2+1) = 1 ;
317 big_res = oper.inverse(sec_membre) ;
318 for (int i=0 ; i<nr ; i++) {
319 sol_hom_quatre_eta.set(zone, k, j, i) = big_res(i) ;
320 sol_hom_quatre_vr.set(zone, k, j, i) = big_res(nr+i) ;
321 }
322
323 }
324 }
325 }
326 }
327
328 // Compactified external domain
329 //-----------------------------
330 nr = mg.get_nr(nzm1) ;
331 assert(nt == mg.get_nt(nzm1)) ;
332 assert(np == mg.get_np(nzm1)) ;
333 alpha = mpaff->get_alpha()[nzm1] ; alp2 = alpha*alpha ;
334 assert (nr > 4) ;
335
336 // Loop on l and m
337 //----------------
338 for (int k=0 ; k<np+1 ; k++) {
339 for (int j=0 ; j<nt ; j++) {
340 base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
341 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
342 int dege1 = 3; //degeneracy of eq.1
343 int dege2 = (l_q == 1 ? 2 : 3); //degeneracy of eq.2
344 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
345 int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
346 int nrtot = 2*nr ;
347 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
348 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
349 Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
350 Diff_sxdsdx sxd(base_r, nr) ; const Matrice& mxd = sxd.get_matrice() ;
351 Diff_sx2 sx2(base_r, nr) ; const Matrice& ms2 = sx2.get_matrice() ;
352
353 // Building the operator
354 //----------------------
355 for (int lin=0; lin<nr_eq1; lin++) {
356 for (int col=0; col<nr; col++)
357 oper.set(lin,col)
358 = (md2(lin,col) - (lam+1)*l_q*(l_q+1)*ms2(lin,col))/alp2 ;
359 for (int col=0; col<nr; col++)
360 oper.set(lin,col+nr) =
361 (-lam*mxd(lin,col) + 2*(1+lam)*ms2(lin,col)) / alp2 ;
362 }
363 for (int col=0; col<nr; col++)
364 oper.set(nr_eq1, col) = 1 ;
365 for (int col=nr; col<nrtot; col++)
366 oper.set(nr_eq1, col) = 0 ;
367 int pari = -1 ;
368 for (int col=0; col<nr; col++) {
369 oper.set(nr_eq1+1, col) = pari*col*col ;
370 pari = pari ;
371 }
372 for (int col=nr; col<nrtot; col++)
373 oper.set(nr_eq1+1, col) = 0 ;
374 oper.set(nr_eq1+2, 0) = 1 ;
375 for (int col=1; col<nrtot; col++)
376 oper.set(nr_eq1+2, col) = 0 ;
377 for (int lin=0; lin<nr_eq2; lin++) {
378 for (int col=0; col<nr; col++)
379 oper.set(lin+nr,col) = (l_q*(l_q+1)*(lam*mxd(lin,col)
380 + (lam+2)*ms2(lin,col))) / alp2 ;
381 for (int col=0; col<nr; col++)
382 oper.set(lin+nr, col+nr) = ((lam+1)*md2(lin,col)
383 - (2*(lam+1) + l_q*(l_q+1))*ms2(lin,col)) / alp2 ;
384 }
385 for (int col=0; col<nr; col++)
386 oper.set(nr+nr_eq2, col) = 0 ;
387 for (int col=nr; col<nrtot; col++)
388 oper.set(nr+nr_eq2, col) = 1 ;
389 for (int col=0; col<nr; col++)
390 oper.set(nr+nr_eq2+1, col) = 0 ;
391 pari = -1 ;
392 for (int col=0; col<nr; col++) {
393 oper.set(nr+nr_eq2+1, nr+col) = pari*col*col ;
394 pari = pari ;
395 }
396 if (l_q > 1) {
397 for (int col=0; col<nr; col++)
398 oper.set(nr+nr_eq2+2, col) = 0 ;
399 oper.set(nr+nr_eq2+2, nr) = 1 ;
400 for (int col=nr+1; col<nrtot; col++)
401 oper.set(nr+nr_eq2+2, col) = 0 ;
402 }
403 oper.set_lu() ;
404
405 // Filling the r.h.s
406 //------------------
407 for (int i=0; i<nr_eq1; i++)
408 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
409 for (int i=nr_eq1; i<nr; i++)
410 sec_membre.set(i) = 0 ;
411 for (int i=0; i<nr_eq2; i++)
412 sec_membre.set(i+nr) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
413 for (int i=nr_eq2; i<nr; i++)
414 sec_membre.set(nr+i) = 0 ;
415 Tbl big_res = oper.inverse(sec_membre) ;
416
417 // Putting coefficients of h and v to individual arrays
418 //-----------------------------------------------------
419 for (int i=0; i<nr; i++) {
420 sol_part_eta.set(nzm1, k, j, i) = big_res(i) ;
421 sol_part_vr.set(nzm1, k, j, i) = big_res(i+nr) ;
422 }
423
424 // Homogeneous solutions (only 1/r^(l+2) and 1/r^l in the CED)
425 //------------------------------------------------------------
426 if (l_q == 1) {
427 Tbl sol_hom1 = solh(nr, 0, 0., base_r) ;
428 Tbl sol_hom2 = solh(nr, 2, 0., base_r) ;
429 double fac_eta = lam + 2. ;
430 double fac_vr = 2*lam + 2. ;
431 for (int i=0 ; i<nr ; i++) {
432 sol_hom_deux_eta.set(nzm1, k, j, i) = sol_hom2(i) ;
434 sol_hom_deux_vr.set(nzm1, k, j, i) = -2*sol_hom2(i) ;
436 }
437 }
438 else {
439 sec_membre.annule_hard() ;
440 sec_membre.set(nr-1) = 1 ;
441 big_res = oper.inverse(sec_membre) ;
442
443 for (int i=0; i<nr; i++) {
444 sol_hom_deux_eta.set(nzm1, k, j, i) = big_res(i) ;
445 sol_hom_deux_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
446 }
447 sec_membre.set(nr-1) = 0 ;
448 sec_membre.set(2*nr-1) = 1 ;
449 big_res = oper.inverse(sec_membre) ;
450 for (int i=0; i<nr; i++) {
451 sol_hom_quatre_eta.set(nzm1, k, j, i) = big_res(i) ;
452 sol_hom_quatre_vr.set(nzm1, k, j, i) = big_res(nr+i) ;
453 }
454 }
455 }
456 }
457 }
458
459 // Now let's match everything ...
460 //-------------------------------
461
462 // Resulting V^r & eta
463 vr.set_etat_qcq() ;
464 vr.set_spectral_base(base) ;
465 vr.set_spectral_va().set_etat_cf_qcq() ;
466 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
467 cf_vr.annule_hard() ;
468 het.set_etat_qcq() ;
469 het.set_spectral_base(base) ;
470 het.set_spectral_va().set_etat_cf_qcq() ;
471 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
472 cf_eta.annule_hard() ;
473 int taille = 4*nzm1 -2 ; // Two less than R3 case
474 Tbl sec_membre(taille) ;
475 Matrice systeme(taille, taille) ;
476 systeme.set_etat_qcq() ;
477 int ligne ; int colonne ;
478
479 // Loop on l and m
480 //----------------
481 for (int k=0 ; k<np+1 ; k++)
482 for (int j=0 ; j<nt ; j++) {
483 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
484 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
485
486 ligne = 0 ;
487 colonne = 0 ;
488 systeme.annule_hard() ;
489 sec_membre.annule_hard() ;
490
491
492 //shell 1 (boundary condition)
493 int zone1 = 1;
494 nr = mg.get_nr(zone1) ;
495 alpha = mpaff->get_alpha()[zone1] ;
496
497 // Here we prepare for a Robyn-type boundary condition on eta.
498 // Parameters are dir_eta and neum_eta
499 systeme.set(ligne, colonne)
500 = -dir_eta*sol_hom_un_eta.val_in_bound_jk(zone1, j, k) ;
501 systeme.set(ligne, colonne+1)
502 = -dir_eta * sol_hom_deux_eta.val_in_bound_jk(zone1, j, k) ;
503 systeme.set(ligne, colonne+2)
504 = -dir_eta*sol_hom_trois_eta.val_in_bound_jk(zone1, j, k) ;
505 systeme.set(ligne, colonne+3)
506 = -dir_eta*sol_hom_quatre_eta.val_in_bound_jk(zone1, j, k) ;
507
508
509
510
511 sec_membre.set(ligne) += dir_eta* sol_part_eta.val_in_bound_jk(zone1, j, k) ;
512
513
514
515 int pari = -1 ;
516 for (int i=0; i<nr; i++) {
517 systeme.set(ligne, colonne)
518 -= neum_eta*pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
519 systeme.set(ligne, colonne+1)
520 -= neum_eta*pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
521 systeme.set(ligne, colonne+2)
522 -= neum_eta*pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
523 systeme.set(ligne, colonne+3)
524 -= neum_eta*pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
525 sec_membre.set(ligne)
526 += neum_eta*pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
527 pari = -pari ;
528 }
529
530 sec_membre.set(ligne) -= (*bound_eta).val_in_bound_jk(1,j,k);
531 ligne++ ;
532
533
534
535 // ... and their couterparts for V^r
536 systeme.set(ligne, colonne)
537 = - dir_vr*sol_hom_un_vr.val_in_bound_jk(zone1, j, k) ;
538 systeme.set(ligne, colonne+1)
539 = - dir_vr*sol_hom_deux_vr.val_in_bound_jk(zone1, j, k) ;
540 systeme.set(ligne, colonne+2)
541 = -dir_vr*sol_hom_trois_vr.val_in_bound_jk(zone1, j, k) ;
542 systeme.set(ligne, colonne+3)
543 = -dir_vr*sol_hom_quatre_vr.val_in_bound_jk(zone1, j, k) ;
544 sec_membre.set(ligne) += dir_vr*sol_part_vr.val_in_bound_jk(zone1, j, k) ;
545
546
547
548
549 pari = -1 ;
550 for (int i=0; i<nr; i++) {
551 systeme.set(ligne, colonne)
552 -= neum_vr*pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
553 systeme.set(ligne, colonne+1)
554 -= neum_vr*pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
555 systeme.set(ligne, colonne+2)
556 -= neum_vr*pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
557 systeme.set(ligne, colonne+3)
558 -= neum_vr*pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
559 sec_membre.set(ligne)
560 += neum_vr*pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
561 pari = -pari ;
562 }
563
564 sec_membre.set(ligne) -= (*bound_vr).val_in_bound_jk(1,j,k);
565 ligne++ ;
566
567
568
569
570 // Values in 1: beginning with solution in eta...
571
572 systeme.set(ligne, colonne)
573 += sol_hom_un_eta.val_out_bound_jk(zone1, j, k) ;
574 systeme.set(ligne, colonne+1)
575 += sol_hom_deux_eta.val_out_bound_jk(zone1, j, k) ;
576 systeme.set(ligne, colonne+2)
577 += sol_hom_trois_eta.val_out_bound_jk(zone1, j, k) ;
578 systeme.set(ligne, colonne+3)
579 += sol_hom_quatre_eta.val_out_bound_jk(zone1, j, k) ;
580 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone1, j, k) ;
581 ligne++ ;
582 // ... and their counterparts for V^r
583 systeme.set(ligne, colonne)
584 += sol_hom_un_vr.val_out_bound_jk(zone1, j, k) ;
585 systeme.set(ligne, colonne+1)
586 += sol_hom_deux_vr.val_out_bound_jk(zone1, j, k) ;
587 systeme.set(ligne, colonne+2)
588 += sol_hom_trois_vr.val_out_bound_jk(zone1, j, k) ;
589 systeme.set(ligne, colonne+3)
590 += sol_hom_quatre_vr.val_out_bound_jk(zone1, j, k) ;
591 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone1, j, k) ;
592 ligne++ ;
593
594 //derivative at 1
595 pari = 1 ;
596 for (int i=0; i<nr; i++) {
597 systeme.set(ligne, colonne)
598 += pari*i*i*sol_hom_un_eta(zone1, k, j, i)/alpha ;
599 systeme.set(ligne, colonne+1)
600 += pari*i*i*sol_hom_deux_eta(zone1, k, j, i)/alpha ;
601 systeme.set(ligne, colonne+2)
602 += pari*i*i*sol_hom_trois_eta(zone1, k, j, i)/alpha ;
603 systeme.set(ligne, colonne+3)
604 += pari*i*i*sol_hom_quatre_eta(zone1, k, j, i)/alpha ;
605 sec_membre.set(ligne)
606 -= pari*i*i* sol_part_eta(zone1, k, j, i)/alpha ;
607 pari = pari ;
608 }
609 ligne++ ;
610 // ... and their counterparts for V^r
611 pari = 1 ;
612 for (int i=0; i<nr; i++) {
613 systeme.set(ligne, colonne)
614 += pari*i*i*sol_hom_un_vr(zone1, k, j, i)/alpha ;
615 systeme.set(ligne, colonne+1)
616 += pari*i*i*sol_hom_deux_vr(zone1, k, j, i)/alpha ;
617 systeme.set(ligne, colonne+2)
618 += pari*i*i*sol_hom_trois_vr(zone1, k, j, i)/alpha ;
619 systeme.set(ligne, colonne+3)
620 += pari*i*i*sol_hom_quatre_vr(zone1, k, j, i)/alpha ;
621 sec_membre.set(ligne)
622 -= pari*i*i* sol_part_vr(zone1, k, j, i)/alpha ;
623 pari = pari ;
624 }
625 colonne += 4 ;
626
627 // Other shells
628 if ( nz <= 2){
629 cout <<"WARNING!! Mapping must contain at least 2 shells!!" << endl;}
630 for (int zone=2 ; zone<nzm1 ; zone++) {
631
632 nr = mg.get_nr(zone) ;
633 alpha = mpaff->get_alpha()[zone] ;
634 ligne -= 3 ;
635 systeme.set(ligne, colonne)
636 = -sol_hom_un_eta.val_in_bound_jk(zone, j, k) ;
637 systeme.set(ligne, colonne+1)
638 = -sol_hom_deux_eta.val_in_bound_jk(zone, j, k) ;
639 systeme.set(ligne, colonne+2)
640 = -sol_hom_trois_eta.val_in_bound_jk(zone, j, k) ;
641 systeme.set(ligne, colonne+3)
642 = -sol_hom_quatre_eta.val_in_bound_jk(zone, j, k) ;
643 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
644 ligne++ ;
645 // ... and their counterparts for V^r
646 systeme.set(ligne, colonne)
647 = -sol_hom_un_vr.val_in_bound_jk(zone, j, k) ;
648 systeme.set(ligne, colonne+1)
649 = -sol_hom_deux_vr.val_in_bound_jk(zone, j, k) ;
650 systeme.set(ligne, colonne+2)
651 = -sol_hom_trois_vr.val_in_bound_jk(zone, j, k) ;
652 systeme.set(ligne, colonne+3)
653 = -sol_hom_quatre_vr.val_in_bound_jk(zone, j, k) ;
654 sec_membre.set(ligne) += sol_part_vr.val_in_bound_jk(zone, j, k) ;
655 ligne++ ;
656
657 //derivative of (x+echelle)^(l-1) at -1
658 pari = -1 ;
659 for (int i=0; i<nr; i++) {
660 systeme.set(ligne, colonne)
661 -= pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
662 systeme.set(ligne, colonne+1)
663 -= pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
664 systeme.set(ligne, colonne+2)
665 -= pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
666 systeme.set(ligne, colonne+3)
667 -= pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
668 sec_membre.set(ligne)
669 += pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
670 pari = -pari ;
671 }
672 ligne++ ;
673 // ... and their counterparts for V^r
674 pari = -1 ;
675 for (int i=0; i<nr; i++) {
676 systeme.set(ligne, colonne)
677 -= pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
678 systeme.set(ligne, colonne+1)
679 -= pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
680 systeme.set(ligne, colonne+2)
681 -= pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
682 systeme.set(ligne, colonne+3)
683 -= pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
684 sec_membre.set(ligne)
685 += pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
686 pari = -pari ;
687 }
688 ligne++ ;
689
690 systeme.set(ligne, colonne)
691 += sol_hom_un_eta.val_out_bound_jk(zone, j, k) ;
692 systeme.set(ligne, colonne+1)
693 += sol_hom_deux_eta.val_out_bound_jk(zone, j, k) ;
694 systeme.set(ligne, colonne+2)
695 += sol_hom_trois_eta.val_out_bound_jk(zone, j, k) ;
696 systeme.set(ligne, colonne+3)
697 += sol_hom_quatre_eta.val_out_bound_jk(zone, j, k) ;
698 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
699 ligne++ ;
700 // ... and their counterparts for V^r
701 systeme.set(ligne, colonne)
702 += sol_hom_un_vr.val_out_bound_jk(zone, j, k) ;
703 systeme.set(ligne, colonne+1)
704 += sol_hom_deux_vr.val_out_bound_jk(zone, j, k) ;
705 systeme.set(ligne, colonne+2)
706 += sol_hom_trois_vr.val_out_bound_jk(zone, j, k) ;
707 systeme.set(ligne, colonne+3)
708 += sol_hom_quatre_vr.val_out_bound_jk(zone, j, k) ;
709 sec_membre.set(ligne) -= sol_part_vr.val_out_bound_jk(zone, j, k) ;
710 ligne++ ;
711
712 //derivative at 1
713 pari = 1 ;
714 for (int i=0; i<nr; i++) {
715 systeme.set(ligne, colonne)
716 += pari*i*i*sol_hom_un_eta(zone, k, j, i)/alpha ;
717 systeme.set(ligne, colonne+1)
718 += pari*i*i*sol_hom_deux_eta(zone, k, j, i)/alpha ;
719 systeme.set(ligne, colonne+2)
720 += pari*i*i*sol_hom_trois_eta(zone, k, j, i)/alpha ;
721 systeme.set(ligne, colonne+3)
722 += pari*i*i*sol_hom_quatre_eta(zone, k, j, i)/alpha ;
723 sec_membre.set(ligne)
724 -= pari*i*i* sol_part_eta(zone, k, j, i)/alpha ;
725 pari = pari ;
726 }
727 ligne++ ;
728 // ... and their couterparts for V^r
729 pari = 1 ;
730 for (int i=0; i<nr; i++) {
731 systeme.set(ligne, colonne)
732 += pari*i*i*sol_hom_un_vr(zone, k, j, i)/alpha ;
733 systeme.set(ligne, colonne+1)
734 += pari*i*i*sol_hom_deux_vr(zone, k, j, i)/alpha ;
735 systeme.set(ligne, colonne+2)
736 += pari*i*i*sol_hom_trois_vr(zone, k, j, i)/alpha ;
737 systeme.set(ligne, colonne+3)
738 += pari*i*i*sol_hom_quatre_vr(zone, k, j, i)/alpha ;
739 sec_membre.set(ligne)
740 -= pari*i*i* sol_part_vr(zone, k, j, i)/alpha ;
741 pari = pari ;
742 }
743 colonne += 4 ;
744 }
745
746 //Compactified external domain
747 nr = mg.get_nr(nzm1) ;
748
749 alpha = mpaff->get_alpha()[nzm1] ;
750 ligne -= 3 ;
751 //value of (x-1)^(l+2) at -1 :
752 systeme.set(ligne, colonne)
753 -= sol_hom_deux_eta.val_in_bound_jk(nzm1, j, k) ;
754 systeme.set(ligne, colonne+1)
755 -= sol_hom_quatre_eta.val_in_bound_jk(nzm1, j, k) ;
756 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nzm1, j, k) ;
757 //... and of its counterpart for V^r
758 systeme.set(ligne+1, colonne)
759 -= sol_hom_deux_vr.val_in_bound_jk(nzm1, j, k) ;
760 systeme.set(ligne+1, colonne+1)
761 -= sol_hom_quatre_vr.val_in_bound_jk(nzm1, j, k) ;
762 sec_membre.set(ligne+1) += sol_part_vr.val_in_bound_jk(nzm1, j, k) ;
763
764 ligne += 2 ;
765 //derivative of (x-1)^(l+2) at -1 :
766 pari = 1 ;
767 for (int i=0; i<nr; i++) {
768 systeme.set(ligne, colonne)
769 -= 4*alpha*pari*i*i*sol_hom_deux_eta(nzm1, k, j, i) ;
770 systeme.set(ligne, colonne+1)
771 -= 4*alpha*pari*i*i*sol_hom_quatre_eta(nzm1, k, j, i) ;
772 sec_membre.set(ligne)
773 += 4*alpha*pari*i*i* sol_part_eta(nzm1, k, j, i) ;
774 pari = -pari ;
775 }
776 ligne++ ;
777 // ... and their counterparts for V^r
778 pari = 1 ;
779 for (int i=0; i<nr; i++) {
780 systeme.set(ligne, colonne)
781 -= 4*alpha*pari*i*i*sol_hom_deux_vr(nzm1, k, j, i) ;
782 systeme.set(ligne, colonne+1)
783 -= 4*alpha*pari*i*i*sol_hom_quatre_vr(nzm1, k, j, i) ;
784 sec_membre.set(ligne)
785 += 4*alpha*pari*i*i* sol_part_vr(nzm1, k, j, i) ;
786 pari = -pari ;
787 }
788
789 // Solution of the system giving the coefficients for the homogeneous
790 // solutions
791 //-------------------------------------------------------------------
792 systeme.set_lu() ;
793 Tbl facteurs(systeme.inverse(sec_membre)) ;
794 int conte = 0 ;
795
796 // everything is put to the right place, the same combination of hom.
797 // solutions (with some l or -(l+1) factors) must be used for V^r
798 //-------------------------------------------------------------------
799 nr = mg.get_nr(0) ; //nucleus
800 for (int i=0 ; i<nr ; i++) {
801 cf_eta.set(0, k, j, i) = 0.;
802 cf_vr.set(0, k, j, i) = 0.;
803 }
804
805 for (int zone=1 ; zone<nzm1 ; zone++) { //shells
806 nr = mg.get_nr(zone) ;
807 for (int i=0 ; i<nr ; i++) {
808 cf_eta.set(zone, k, j, i) =
809 sol_part_eta(zone, k, j, i)
814 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
819 }
820 conte+=4 ;
821 }
822 nr = mg.get_nr(nz-1) ; //compactified external domain
823 for (int i=0 ; i<nr ; i++) {
824 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
827 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
830
831 }
832 } // End of nullite_plm
833 } //End of loop on theta
834 vr.set_spectral_va().ylm_i() ;
835 vr += vrl0 ;
836 het.set_spectral_va().ylm_i() ;
837
838
839
840 Scalar pipo(*mpaff);
841 pipo.annule_hard();
842 pipo.std_spectral_base();
843 pipo.set_dzpuis(3);
845
846
847
848
849 resu.set_vr_eta_mu(vr, het, mu().sol_elliptic_boundary(param_mu, (*bound_mu), dir_mu , neum_mu)) ;
850
851 return ;
852
853}
854}
Bases of the spectral expansions.
Definition base_val.h:322
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:172
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:369
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Definition diff_sx2.C:98
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:450
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:490
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:531
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
This class contains the parameters needed to call the general elliptic solver.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Basic array class.
Definition tbl.h:161
double * t
The array of double.
Definition tbl.h:173
Tensor field of valence 1.
Definition vector.h:188
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64