LORENE
vector_df_poisson.C
1/*
2 * Resolution of the divergence-free vector Poisson equation
3 *
4 * (see file vector.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2003 Eric Gourgoulhon & 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_df_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $" ;
29
30/*
31 * $Id: vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $
32 * $Log: vector_df_poisson.C,v $
33 * Revision 1.15 2014/10/13 08:53:44 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.14 2014/10/06 15:13:20 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.13 2005/02/15 09:45:22 j_novak
40 * Correction of an error in the matching.
41 *
42 * Revision 1.12 2005/02/09 16:53:11 j_novak
43 * Now V^r and eta are matched across domains, but not any of their derivatives.
44 *
45 * Revision 1.11 2005/02/09 14:52:01 j_novak
46 * Better solution in the shells.
47 *
48 * Revision 1.10 2005/02/09 13:20:27 j_novak
49 * Completely new way of solving the vector Poisson equation in the Div_free
50 * case: the system is inverted "as a whole" for V^r and eta. This works only
51 * with Map_af...
52 *
53 *
54 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_df_poisson.C,v 1.15 2014/10/13 08:53:44 j_novak Exp $
55 *
56 */
57
58// C headers
59#include <cassert>
60#include <cmath>
61
62// Lorene headers
63#include "tensor.h"
64#include "diff.h"
65#include "proto.h"
66#include "param.h"
67
68namespace Lorene {
70
71 // All this has a meaning only for spherical components:
72#ifndef NDEBUG
73 const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
74 assert(bvs != 0x0) ;
75#endif
76
77 int nitermax = par.get_int() ;
78 int& niter = par.get_int_mod() ;
79 double relax = par.get_double() ;
80 double precis = par.get_double(1) ;
81 Cmp& ss_khi = par.get_cmp_mod(0) ;
82 Cmp& ss_mu = par.get_cmp_mod(1) ;
83
84 // Solution for the r-component
85 // ----------------------------
86
87 Scalar source_r = *(cmp[0]) ;
88 source_r.mult_r() ;
89
91 par_khi.add_int(nitermax, 0) ;
92 par_khi.add_int_mod(niter, 0) ;
93 par_khi.add_double(relax, 0) ;
94 par_khi.add_double(precis, 1) ;
95 par_khi.add_cmp_mod(ss_khi, 0) ;
96
97 Scalar khi (*mp) ;
98 khi.set_etat_zero() ;
99
100 source_r.poisson(par_khi, khi) ;
101 khi.div_r() ; // khi now contains V^r
102
103 // Solution for mu
104 // ---------------
105
106 Param par_mu ;
107 par_mu.add_int(nitermax, 0) ;
108 par_mu.add_int_mod(niter, 0) ;
109 par_mu.add_double(relax, 0) ;
110 par_mu.add_double(precis, 1) ;
111 par_mu.add_cmp_mod(ss_mu, 0) ;
112
113 Scalar mu_resu (*mp) ;
114 mu_resu.set_etat_zero() ;
115
116 mu().poisson(par_mu, mu_resu) ;
117
118 // Final result
119 // ------------
120
122
123 resu.set_vr_mu(khi, mu_resu) ;
124
125 return resu ;
126
127}
128
129/*
130 * In the case without parameters, first is solved the equation for mu and then
131 * the system of equations for (eta, V^r) is inverted as a whole:
132 * d2 eta / dr2 + 2/r d eta / dr - 1/r d V^r / dr = S(eta)
133 * d V^r / dr + 2/r V^r - l(l+1)/r eta = 0 (div free condition)
134 *
135 * There is no l=0 contribution (divergence free in all space!).
136 * In the nucleus and the CED the system is inverted for h(r) and v(r) ,
137 * such that eta = r^2 h and V^r = r^2 v in the nucleus,
138 * in the compactified domain one has eta = u^2 h and V^r = u^2 v (where u=1/r);
139 * In the shells, both equations are multiplied by r.
140 * These methods are used only to get particular solutions.
141 *
142 * Homogeneous solutions are known analitically: r^(l-1) and/or 1/r^(l+2)
143 * It is then only possible to match eta and V^r, but not their derivatives,
144 * due to the div-free condition.
145 */
147
148 const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
149#ifndef NDEBUG
150 for (int i=0; i<3; i++)
151 assert(cmp[i]->check_dzpuis(4)) ;
152 // All this has a meaning only for spherical components:
153 const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
154 assert(bvs != 0x0) ;
155 //## ... and affine mapping, for the moment!
156 assert( mpaff != 0x0) ;
157#endif
158
159 // Final result
160 // ------------
162
163 // Solution for mu
164 // ---------------
165 Scalar mu_resu = mu().poisson() ;
166
167 Scalar f_r(*mpaff) ;
168 if (cmp[0]->get_etat() == ETATZERO) { // no work needed ...
169 f_r.set_etat_zero() ;
170 resu.set_vr_mu(f_r, mu_resu) ;
171 return resu ;
172 }
173
174 // Some working objects
175 //---------------------
176 const Mg3d& mg = *mpaff->get_mg() ;
177 int nz = mg.get_nzone() ; int nzm1 = nz - 1;
178 assert(mg.get_type_r(0) == RARE) ;
179 assert(mg.get_type_r(nzm1) == UNSURR) ;
180 Scalar S_r = *cmp[0] ;
181 Scalar S_eta = eta() ;
182 S_r.set_spectral_va().ylm() ;
183 S_eta.set_spectral_va().ylm() ;
184 const Base_val& base = S_eta.get_spectral_va().base ;
185 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
186 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
187 Mtbl_cf solution_hom_un(mg, base) ; solution_hom_un.annule_hard() ;
188 Mtbl_cf solution_hom_deux(mg, base) ; solution_hom_deux.annule_hard() ;
189
190 // Build-up & inversion of the system for (eta, V^r) in each domain
191 //-----------------------------------------------------------------
192
193 // Nucleus
194 //--------
195 int nr = mg.get_nr(0) ;
196 int nt = mg.get_nt(0) ;
197 int np = mg.get_np(0) ;
198 double alpha = mpaff->get_alpha()[0] ;
199 double beta = mpaff->get_beta()[0] ;
200 int l_q = 0 ; int m_q = 0; int base_r = 0 ;
201 int nr0 = nr - 1 ; //one degree of freedom less because of division by r^2
202
203 // Loop on l and m
204 //----------------
205 for (int k=0 ; k<np+1 ; k++) {
206 for (int j=0 ; j<nt ; j++) {
207 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
208 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
209 int dege1 = (l_q <3 ? 0 : 1) ; //degeneracy of eq.1
210 int dege2 = 0 ; //degeneracy of eq.2
211 int nr_eq1 = nr0 - dege1 ; //Eq.1 is for h (eta/r^2)
212 int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
213 int nrtot = nr_eq1 + nr_eq2 ;
214 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
215 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
216 Diff_x2dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
217 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
218 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
219
220 // Building the operator
221 //----------------------
222 for (int lin=0; lin<nr_eq1; lin++) { //eq.1
223 for (int col=dege1; col<nr0; col++)
224 oper.set(lin,col-dege1)
225 = md2(lin,col) + 6*mxd(lin,col) + 6*mid(lin,col) ;
226 for (int col=dege2; col<nr0; col++)
227 oper.set(lin,col-dege2+nr_eq1) = -mxd(lin,col) -2*mid(lin,col) ;
228 }
229 for (int lin=0; lin<nr0; lin++) { //eq.2
230 for (int col=dege1; col<nr0; col++)
231 oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
232 for (int col=dege2; col<nr0; col++)
233 oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) + 4*mid(lin,col) ;
234 }
235 oper.set_lu() ;
236
237 // Filling the r.h.s
238 //------------------
239 for (int i=0; i<nr_eq1; i++) //eq.1
240 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(0, k, j, i) ;
241 for (int i=0; i<nr0; i++) //eq.2
242 sec_membre.set(i+nr_eq1) = 0 ;
243
244 // Inversion of the "big" operator
245 //--------------------------------
246 Tbl big_res = oper.inverse(sec_membre) ;
247 Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
248 Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
249
250 // Putting coefficients of h and v to individual arrays
251 //-----------------------------------------------------
252 for (int i=0; i<dege1; i++)
253 res_eta.set(i) = 0 ;
254 for (int i=dege1; i<nr0; i++)
255 res_eta.set(i) = big_res(i-dege1) ;
256 res_eta.set(nr0) = 0 ;
257 for (int i=0; i<dege2; i++)
258 res_vr.set(i) = 0 ;
259 for (int i=dege2; i<nr0; i++)
260 res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
261 res_vr.set(nr0) = 0 ;
262
263 // Multiplication by xi^2 (the alpha^2 factor comes soon)
264 //-------------------------------------------------------
265 multx2_1d(nr, &res_eta.t, base_r) ;
266 multx2_1d(nr, &res_vr.t, base_r) ;
267
268 // Homogeneous solution (only r^(l-1) in the nucleus)
269 Tbl sol_hom = solh(nr, l_q-1, 0., base_r) ;
270 for (int i=0 ; i<nr ; i++) {
271 sol_part_eta.set(0, k, j, i) = alpha*alpha*res_eta(i) ;
272 sol_part_vr.set(0, k, j, i) = alpha*alpha*res_vr(i) ;
273 solution_hom_un.set(0, k, j, i) = sol_hom(i) ;
274 solution_hom_deux.set(0, k, j, i) = 0. ;
275 }
276 }
277 }
278 }
279
280
281 // Shells
282 //-------
283 for (int zone=1 ; zone<nzm1 ; zone++) {
284 nr = mg.get_nr(zone) ;
285 assert (nr > 5) ;
286 assert(nt == mg.get_nt(zone)) ;
287 assert(np == mg.get_np(zone)) ;
288 alpha = mpaff->get_alpha()[zone] ;
289 beta = mpaff->get_beta()[zone] ;
290 double ech = beta / alpha ;
291
292 // Loop on l and m
293 //----------------
294 for (int k=0 ; k<np+1 ; k++) {
295 for (int j=0 ; j<nt ; j++) {
296 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
297 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
298 int dege1 = 2 ; //degeneracy of eq.1
299 int dege2 = 0 ; //degeneracy of eq.2
300 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
301 int nr_eq2 = nr - dege2 ; //Eq.2 is the div-free condition
302 int nrtot = nr_eq1 + nr_eq2 + 1;
303 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
304 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
305 Diff_x2dsdx2 x2d2(base_r, nr+1); const Matrice& m2d2 = x2d2.get_matrice() ;
306 Diff_xdsdx2 xd2(base_r, nr+1) ; const Matrice& mxd2 = xd2.get_matrice() ;
307 Diff_dsdx2 d2(base_r, nr+1) ; const Matrice& md2 = d2.get_matrice() ;
308 Diff_xdsdx xd(base_r, nr+1) ; const Matrice& mxd = xd.get_matrice() ;
309 Diff_dsdx d1(base_r, nr+1) ; const Matrice& md = d1.get_matrice() ;
310 Diff_id id(base_r, nr+1) ; const Matrice& mid = id.get_matrice() ;
311
312 // Building the operator
313 //----------------------
314 for (int lin=0; lin<nr_eq1; lin++) {
315 for (int col=dege1; col<nr; col++)
316 oper.set(lin,col-dege1)
317 = mxd2(lin,col) + ech*md2(lin,col) + 2*md(lin,col) ;
318 for (int col=dege2; col<nr+1; col++)
319 oper.set(lin,col-dege2+nr_eq1) = -md(lin,col) ;
320 }
321 for (int lin=0; lin<nr_eq2; lin++) {
322 for (int col=dege1; col<nr; col++)
323 oper.set(lin+nr_eq1,col-dege1) = -l_q*(l_q+1)*mid(lin,col) ;
324 for (int col=dege2; col<nr+1; col++)
325 oper.set(lin+nr_eq1, col-dege2+nr_eq1) =
326 mxd(lin,col) + ech*md(lin,col) + 2*mid(lin,col) ;
327 }
328 //Additional line to avoid spurious homogeneous solutions
329 //this line is the first one of the V^r eq.
330 for (int col=dege1; col<nr; col++)
331 oper.set(nrtot-1, col-dege1) = 0 ;
332 for (int col=dege2; col<nr+1; col++)
333 oper.set(nrtot-1, col-dege2+nr_eq1) =
334 m2d2(0,col) + ech*(2*mxd2(0,col) + ech*md2(0,col))
335 +4*(mxd(0,col) +ech*md(0,col))
336 +(2 - l_q*(l_q+1))*mid(0,col) ;
337 oper.set_lu() ;
338
339 // Filling the r.h.s
340 //------------------
341 Tbl sr(5) ; sr.set_etat_qcq() ;
342 Tbl seta(nr) ; seta.set_etat_qcq() ;
343 for (int i=0; i<5; i++) {
344 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
345 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
346 }
347 for (int i=5; i<nr; i++)
348 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
349 Tbl xsr= sr ; Tbl x2sr= sr ;
350 Tbl xseta= seta ;
351 multx2_1d(5, &x2sr.t, base_r) ; multx_1d(5, &xsr.t, base_r) ;
352 multx_1d(nr, &xseta.t, base_r) ;
353
354 for (int i=0; i<nr_eq1; i++)
355 sec_membre.set(i) = alpha*(alpha*xseta(i) + beta*seta(i)) ;
356 for (int i=0; i<nr_eq2; i++)
357 sec_membre.set(i+nr_eq1) = 0 ;
358 sec_membre.set(nr_eq1+nr_eq2) = alpha*alpha*x2sr(0) + 2*alpha*beta*xsr(0)
359 + beta*beta*sr(0) ;
360
361 // Inversion of the "big" operator
362 //--------------------------------
363 Tbl big_res = oper.inverse(sec_membre) ;
364 Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
365 Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
366
367 // Putting coefficients of h and v to individual arrays
368 //-----------------------------------------------------
369 for (int i=0; i<dege1; i++)
370 res_eta.set(i) = 0 ;
371 for (int i=dege1; i<nr; i++)
372 res_eta.set(i) = big_res(i-dege1) ;
373 for (int i=0; i<dege2; i++)
374 res_vr.set(i) = 0 ;
375 for (int i=dege2; i<nr; i++)
376 res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
377
378 //homogeneous solutions
379 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
380 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
381 for (int i=0 ; i<nr ; i++) {
382 sol_part_eta.set(zone, k, j, i) = res_eta(i) ;
383 sol_part_vr.set(zone, k, j, i) = res_vr(i) ;
384 solution_hom_un.set(zone, k, j, i) = sol_hom1(0,i) ;
385 solution_hom_deux.set(zone, k, j, i) = sol_hom2(1,i) ;
386 }
387 }
388 }
389 }
390 }
391
392 // Compactified external domain
393 //-----------------------------
394 nr = mg.get_nr(nzm1) ;
395 assert(nt == mg.get_nt(nzm1)) ;
396 assert(np == mg.get_np(nzm1)) ;
397 alpha = mpaff->get_alpha()[nzm1] ;
398 assert (nr > 4) ;
399 nr0 = nr - 2 ; //two degrees of freedom less because of division by r^2
400
401 // Loop on l and m
402 //----------------
403 for (int k=0 ; k<np+1 ; k++) {
404 for (int j=0 ; j<nt ; j++) {
405 base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
406 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
407 int dege1 = 0; //degeneracy of eq.1
408 int dege2 = 1; //degeneracy of eq.2
409 int nr_eq1 = nr0 - dege1 ; //Eq.1 is for eta
410 int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
411 int nrtot = nr_eq1 + nr_eq2 ;
412 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
413 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
414 Diff_x2dsdx2 x2d2(base_r, nr) ; const Matrice& m2d2 = x2d2.get_matrice() ;
415 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
416 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
417
418 // Building the operator
419 //----------------------
420 for (int lin=0; lin<nr_eq1; lin++) {
421 for (int col=dege1; col<nr0; col++)
422 oper.set(lin,col-dege1)
423 = m2d2(lin,col) + 4*mxd(lin,col) + 2*mid(lin,col) ;
424 for (int col=dege2; col<nr0; col++)
425 oper.set(lin,col-dege2+nr_eq1) =
426 mxd(lin,col) + 2*mid(lin,col) ;
427 }
428 for (int lin=0; lin<nr_eq2; lin++) {
429 for (int col=dege1; col<nr0; col++)
430 oper.set(lin+nr_eq1,col-dege1) = l_q*(l_q+1)*mid(lin,col) ;
431 for (int col=dege2; col<nr0; col++)
432 oper.set(lin+nr_eq1, col-dege2+nr_eq1) = mxd(lin,col) ;
433 }
434 oper.set_lu() ;
435
436 // Filling the r.h.s
437 //------------------
438 for (int i=0; i<nr_eq1; i++)
439 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
440 for (int i=0; i<nr_eq2; i++)
441 sec_membre.set(i+nr_eq1) = 0 ;
442 Tbl big_res = oper.inverse(sec_membre) ;
443 Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
444 Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
445
446 // Putting coefficients of h and v to individual arrays
447 //-----------------------------------------------------
448 for (int i=0; i<dege1; i++)
449 res_eta.set(i) = 0 ;
450 for (int i=dege1; i<nr0; i++)
451 res_eta.set(i) = big_res(i-dege1) ;
452 res_eta.set(nr0) = 0 ;
453 res_eta.set(nr0+1) = 0 ;
454 for (int i=0; i<dege2; i++)
455 res_vr.set(i) = 0 ;
456 for (int i=dege2; i<nr0; i++)
457 res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
458 res_vr.set(nr0) = 0 ;
459 res_vr.set(nr0+1) = 0 ;
460
461 // Multiplication by r^2
462 //-----------------------
463 Tbl res_v2(nr) ; res_v2.set_etat_qcq() ;
464 Tbl res_e2(nr) ; res_e2.set_etat_qcq() ;
465 mult2_xm1_1d_cheb(nr, res_eta.t, res_e2.t) ;
466 mult2_xm1_1d_cheb(nr, res_vr.t, res_v2.t) ;
467
468 // Homogeneous solution (only 1/r^(l+2) in the CED)
469 Tbl sol_hom = solh(nr, l_q+1, 0., base_r) ;
470 for (int i=0 ; i<nr ; i++) {
471 sol_part_eta.set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
472 sol_part_vr.set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
473 solution_hom_un.set(nzm1, k, j, i) = sol_hom(i) ;
474 solution_hom_deux.set(nzm1, k, j, i) = 0. ;
475 }
476 }
477 }
478 }
479
480 // Now let's match everything ...
481 //-------------------------------
482
483 // Resulting V^r & eta
484 Scalar vr(*mpaff) ; vr.set_etat_qcq() ;
485 vr.set_spectral_base(base) ;
486 vr.set_spectral_va().set_etat_cf_qcq() ;
487 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
488 cf_vr.annule_hard() ;
489 Scalar het(*mpaff) ; het.set_etat_qcq() ;
490 het.set_spectral_base(base) ;
491 het.set_spectral_va().set_etat_cf_qcq() ;
492 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
493 cf_eta.annule_hard() ;
494 int taille = 2*nzm1 ;
495 Tbl sec_membre(taille) ;
496 Matrice systeme(taille, taille) ;
497 systeme.set_etat_qcq() ;
498 int ligne ; int colonne ;
499
500 // Loop on l and m
501 //----------------
502 for (int k=0 ; k<np+1 ; k++)
503 for (int j=0 ; j<nt ; j++) {
504 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
505 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
506
507 ligne = 0 ;
508 colonne = 0 ;
509 sec_membre.annule_hard() ;
510 for (int l=0; l<taille; l++)
511 for (int c=0; c<taille; c++)
512 systeme.set(l,c) = 0 ;
513 //Nucleus
514 nr = mg.get_nr(0) ;
515 alpha = mpaff->get_alpha()[0] ;
516 // value of x^(l-1) at 1 ...
517 systeme.set(ligne, colonne) = 1. ;
518 for (int i=0 ; i<nr ; i++)
519 sec_membre.set(ligne) -= sol_part_eta(0, k, j, i) ;
520 ligne++ ;
521 // ... and of its couterpart for V^r
522 systeme.set(ligne, colonne) = l_q;
523 for (int i=0; i<nr; i++)
524 sec_membre.set(ligne) -= sol_part_vr(0,k,j,i) ;
525 colonne++ ;
526 //shells
527 for (int zone=1 ; zone<nzm1 ; zone++) {
528 nr = mg.get_nr(zone) ;
529 alpha = mpaff->get_alpha()[zone] ;
530 double echelle = mpaff->get_beta()[zone]/alpha ;
531 ligne -- ;
532 //value of (x+echelle)^(l-1) at -1
533 systeme.set(ligne, colonne) = -pow(echelle-1., double(l_q-1)) ;
534 // value of 1/(x+echelle) ^(l+2) at -1
535 systeme.set(ligne, colonne+1) = -1/pow(echelle-1., double(l_q+2)) ;
536 for (int i=0 ; i<nr ; i++)
537 if (i%2 == 0)
538 sec_membre.set(ligne) += sol_part_eta(zone, k, j, i) ;
539 else sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
540 ligne++ ;
541 // ... and their couterparts for V^r
542 systeme.set(ligne, colonne) = -l_q*pow(echelle-1., double(l_q-1)) ;
543 systeme.set(ligne, colonne+1) = (l_q+1)/pow(echelle-1., double(l_q+2));
544 for (int i=0 ; i<nr ; i++)
545 if (i%2 == 0)
546 sec_membre.set(ligne) += sol_part_vr(zone, k, j, i) ;
547 else sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
548 ligne++ ;
549
550 //value of (x+echelle)^(l-1) at 1 :
551 systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
552 // value of 1/(x+echelle)^(l+2) at 1 :
553 systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
554 for (int i=0 ; i<nr ; i++)
555 sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
556 ligne ++ ;
557 //... and their couterparts for V^r
558 systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
559 systeme.set(ligne, colonne+1) = -double(l_q+1)
560 / pow(echelle+1., double(l_q+2)) ;
561 for (int i=0 ; i<nr ; i++)
562 sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i);
563 colonne += 2 ;
564 }
565 //Compactified external domain
566 nr = mg.get_nr(nzm1) ;
567
568 alpha = mpaff->get_alpha()[nzm1] ;
569 ligne -- ;
570 //value of (x-1)^(l+2) at -1 :
571 systeme.set(ligne, colonne) = -pow(-2, double(l_q+2)) ;
572 for (int i=0 ; i<nr ; i++)
573 if (i%2 == 0) sec_membre.set(ligne) += sol_part_eta(nzm1, k, j, i) ;
574 else sec_membre.set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
575 //... and of its couterpart for V^r
576 systeme.set(ligne+1, colonne) = double(l_q+1)*pow(-2, double(l_q+2)) ;
577 for (int i=0 ; i<nr ; i++)
578 if (i%2 == 0) sec_membre.set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
579 else sec_membre.set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
580
581 // Solution of the system giving the coefficients for the homogeneous
582 // solutions
583 //-------------------------------------------------------------------
584 if (taille > 2) systeme.set_band(2,2) ;
585 else systeme.set_band(1,1) ;
586 systeme.set_lu() ;
587 Tbl facteurs(systeme.inverse(sec_membre)) ;
588 int conte = 0 ;
589
590 // everything is put to the right place, the same combination of hom.
591 // solutions (with some l or -(l+1) factors) must be used for V^r
592 //-------------------------------------------------------------------
593 nr = mg.get_nr(0) ; //nucleus
594 for (int i=0 ; i<nr ; i++) {
595 cf_eta.set(0, k, j, i) = sol_part_eta(0, k, j, i)
597 cf_vr.set(0, k, j, i) = sol_part_vr(0, k, j, i)
599 }
600 conte++ ;
601 for (int zone=1 ; zone<nzm1 ; zone++) { //shells
602 nr = mg.get_nr(zone) ;
603 for (int i=0 ; i<nr ; i++) {
604 cf_eta.set(zone, k, j, i) =
605 sol_part_eta(zone, k, j, i)
608 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
611 }
612 conte+=2 ;
613 }
614 nr = mg.get_nr(nz-1) ; //compactified external domain
615 for (int i=0 ; i<nr ; i++) {
616 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
618 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
620 }
621 } // End of nullite_plm
622 } //End of loop on theta
623
624 vr.set_spectral_va().ylm_i() ;
625 het.set_spectral_va().ylm_i() ;
626
627 resu.set_vr_eta_mu(vr, het, mu_resu) ;
628
629 return resu ;
630
631}
632
633}
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
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
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 (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
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void div_r()
Division by r everywhere; dzpuis is not changed.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Basic array class.
Definition tbl.h:161
double * t
The array of double.
Definition tbl.h:173
Divergence-free vectors.
Definition vector.h:724
Vector_divfree poisson() const
Computes the solution of a vectorial Poisson equation with *this as a source:
const Metric *const met_div
Metric with respect to which the divergence is defined.
Definition vector.h:730
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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