LORENE
vector_poisson_boundary.C
1/*
2 * Method for vector Poisson equation inverting eqs. for V^r and eta as a block
3 * (with a boundary condition on the excised sphere).
4 *
5 * (see file vector.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2005 Jerome Novak
11 * Francois Limousin
12 * Jose Luis Jaramillo
13 *
14 * This file is part of LORENE.
15 *
16 * LORENE is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License version 2
18 * as published by the Free Software Foundation.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31char vector_poisson_boundary_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary.C,v 1.3 2014/10/13 08:53:45 j_novak Exp $" ;
32
33/*
34 * $Id: vector_poisson_boundary.C,v 1.3 2014/10/13 08:53:45 j_novak Exp $
35 * $Log: vector_poisson_boundary.C,v $
36 * Revision 1.3 2014/10/13 08:53:45 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.2 2014/10/06 15:13:21 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.1 2005/06/09 08:00:09 f_limousin
43 * Implement a new function sol_elliptic_boundary() and
44 * Vector::poisson_boundary(...) which solve the vectorial poisson
45 * equation (method 6) with an inner boundary condition.
46 *
47 *
48 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_poisson_boundary.C,v 1.3 2014/10/13 08:53:45 j_novak Exp $
49 *
50 */
51
52// C headers
53#include <cassert>
54#include <cstdlib>
55#include <cmath>
56
57// Lorene headers
58#include "metric.h"
59#include "diff.h"
60#include "param_elliptic.h"
61#include "proto.h"
62#include "utilitaires.h"
63
64namespace Lorene {
66 const Mtbl_cf& bound_eta, const Mtbl_cf& bound_mu,
67 int num_front, double fact_dir, double fact_neu,
68 Vector& resu) const {
69
70 const Map_af* mpaff = dynamic_cast<const Map_af*>(mp) ;
71#ifndef NDEBUG
72 for (int i=0; i<3; i++)
73 assert(cmp[i]->check_dzpuis(4)) ;
74 // All this has a meaning only for spherical components:
75 const Base_vect_spher* bvs = dynamic_cast<const Base_vect_spher*>(triad) ;
76 assert(bvs != 0x0) ;
77 //## ... and affine mapping, for the moment!
78 assert( mpaff != 0x0) ;
79#endif
80
81 if (fabs(lam + 1.) < 1.e-8) {
82 cout << "Not implemented yet !!" << endl ;
83 abort() ;
84 /*
85 const Metric_flat& mets = mp->flat_met_spher() ;
86 Vector_divfree sou(*mp, *triad, mets) ;
87 for (int i=1; i<=3; i++) sou.set(i) = *cmp[i-1] ;
88 resu = sou.poisson() ;
89 return ;
90 */
91 }
92
93 // Some working objects
94 //---------------------
95 const Mg3d& mg = *mpaff->get_mg() ;
96 int nz = mg.get_nzone() ; int nzm1 = nz - 1;
97 assert(mg.get_type_r(nzm1) == UNSURR) ;
98 Scalar S_r = *cmp[0] ;
99 if (S_r.get_etat() == ETATZERO) S_r.annule_hard() ;
100 Scalar S_eta = eta() ;
101 if (S_eta.get_etat() == ETATZERO) S_eta.annule_hard() ;
102 S_r.set_spectral_va().ylm() ;
103 S_eta.set_spectral_va().ylm() ;
104 const Base_val& base = S_eta.get_spectral_va().base ;
105 Mtbl_cf sol_part_eta(mg, base) ; sol_part_eta.annule_hard() ;
106 Mtbl_cf sol_part_vr(mg, base) ; sol_part_vr.annule_hard() ;
107 Mtbl_cf solution_hom_un(mg, base) ; solution_hom_un.annule_hard() ;
108 Mtbl_cf solution_hom_deux(mg, base) ; solution_hom_deux.annule_hard() ;
109 Mtbl_cf solution_hom_trois(mg, base) ; solution_hom_trois.annule_hard() ;
110 Mtbl_cf solution_hom_quatre(mg, base) ; solution_hom_quatre.annule_hard() ;
111
112
113 // The l_0 component is solved independently // Understand this step
114 //------------------------------------------
115 Scalar sou_l0 = (*cmp[0]) / (1. + lam) ;
117 for (int l=0; l<nz; l++)
118 param_l0.set_poisson_vect_r(l, true) ;
119
120 // Scalar vrl0 = sou_l0.sol_elliptic(param_l0) ;
121 Scalar vrl0 = sou_l0.sol_elliptic_boundary(param_l0, bound_vr, fact_dir, fact_neu) ;
122
123 // Build-up & inversion of the system for (eta, V^r) in each domain
124 //-----------------------------------------------------------------
125
126 // Shells
127 //-------
128
129 int nr ;
130 int nt = mg.get_nt(0) ;
131 int np = mg.get_np(0) ;
132 int l_q = 0 ; int m_q = 0; int base_r = 0 ;
133 double alpha, beta, ech ;
134
135 assert(num_front+1 < nzm1) ; // Minimum one shell
136 for (int zone=num_front+1 ; zone<nzm1 ; zone++) { //begins the loop on zones
137 nr = mg.get_nr(zone) ;
138 alpha = mpaff->get_alpha()[zone] ;
139 beta = mpaff->get_beta()[zone] ;
140 ech = beta / alpha ;
141 assert (nr > 5) ;
142 assert(nt == mg.get_nt(zone)) ;
143 assert(np == mg.get_np(zone)) ;
144
145 // Loop on l and m
146 //----------------
147 for (int k=0 ; k<np+1 ; k++) {
148 for (int j=0 ; j<nt ; j++) {
149 base.give_quant_numbers(zone, k, j, m_q, l_q, base_r) ;
150 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
151 int dege1 = 2 ; //degeneracy of eq.1
152 int dege2 = 2 ; //degeneracy of eq.2
153 int nr_eq1 = nr - dege1 ; //Eq.1 is for eta
154 int nr_eq2 = nr - dege2 ; //Eq.2 is for V^r
155 int nrtot = nr_eq1 + nr_eq2 ;
156 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
157 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
158 Diff_x2dsdx2 x2d2(base_r, nr); const Matrice& m2d2 = x2d2.get_matrice() ;
159 Diff_xdsdx2 xd2(base_r, nr) ; const Matrice& mxd2 = xd2.get_matrice() ;
160 Diff_dsdx2 d2(base_r, nr) ; const Matrice& md2 = d2.get_matrice() ;
161 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
162 Diff_dsdx d1(base_r, nr) ; const Matrice& md = d1.get_matrice() ;
163 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
164
165 // Building the operator // Which is the eq. from the notes that it is actually implemented?
166 //----------------------
167 for (int lin=0; lin<nr_eq1; lin++) {
168 for (int col=dege1; col<nr; col++)
169 oper.set(lin,col-dege1)
170 = m2d2(lin,col) + 2*ech*mxd2(lin,col) + ech*ech*md2(lin,col)
171 + 2*(mxd(lin,col) + ech*md(lin,col))
172 - (lam+1)*l_q*(l_q+1)*mid(lin,col) ;
173 for (int col=dege2; col<nr; col++)
174 oper.set(lin,col-dege2+nr_eq1)
175 = lam*(mxd(lin,col) + ech*md(lin,col)) + 2*(1+lam)*mid(lin,col) ;
176 }
177 for (int lin=0; lin<nr_eq2; lin++) {
178 for (int col=dege1; col<nr; col++)
179 oper.set(lin+nr_eq1,col-dege1)
180 = -l_q*(l_q+1)*(lam*(mxd(lin,col) + ech*md(lin,col))
181 - (lam+2)*mid(lin,col)) ;
182 for (int col=dege2; col<nr; col++)
184 = (lam+1)*(m2d2(lin,col) + 2*ech*mxd2(lin,col)
185 + ech*ech*md2(lin,col)
186 + 2*(mxd(lin,col) + ech*md(lin,col)))
187 -(2*(lam+1)+l_q*(l_q+1))*mid(lin,col) ;
188 }
189 oper.set_lu() ;
190
191 // Filling the r.h.s
192 //------------------
193 Tbl sr(nr) ; sr.set_etat_qcq() ;
194 Tbl seta(nr) ; seta.set_etat_qcq() ;
195 for (int i=0; i<nr; i++) {
196 sr.set(i) = (*S_r.get_spectral_va().c_cf)(zone, k, j, i);
197 seta.set(i) = (*S_eta.get_spectral_va().c_cf)(zone, k, j, i) ;
198 }
199 Tbl xsr= sr ; Tbl x2sr= sr ;
200 Tbl xseta= seta ; Tbl x2seta = seta ;
201 multx2_1d(nr, &x2sr.t, base_r) ; multx_1d(nr, &xsr.t, base_r) ;
202 multx2_1d(nr, &x2seta.t, base_r) ; multx_1d(nr, &xseta.t, base_r) ;
203
204 for (int i=0; i<nr_eq1; i++)
205 sec_membre.set(i) = alpha*(alpha*x2seta(i) + 2*beta*xseta(i))
206 + beta*beta*seta(i);
207 for (int i=0; i<nr_eq2; i++)
208 sec_membre.set(i+nr_eq1) = beta*beta*sr(i)
209 + alpha*(alpha*x2sr(i) + 2*beta*xsr(i)) ;
210
211 // Inversion of the "big" operator
212 //--------------------------------
213 Tbl big_res = oper.inverse(sec_membre) ;
214 Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
215 Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
216
217 // Putting coefficients of h and v to individual arrays
218 //-----------------------------------------------------
219 for (int i=0; i<dege1; i++)
220 res_eta.set(i) = 0 ;
221 for (int i=dege1; i<nr; i++)
222 res_eta.set(i) = big_res(i-dege1) ;
223 for (int i=0; i<dege2; i++)
224 res_vr.set(i) = 0 ;
225 for (int i=dege2; i<nr; i++)
226 res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
227
228 //homogeneous solutions //I do not understand!!!
229 Tbl sol_hom1 = solh(nr, l_q-1, ech, base_r) ;
230 Tbl sol_hom2 = solh(nr, l_q+1, ech, base_r) ;
231 for (int i=0 ; i<nr ; i++) {
232 sol_part_eta.set(zone, k, j, i) = res_eta(i) ;
233 sol_part_vr.set(zone, k, j, i) = res_vr(i) ;
234 solution_hom_un.set(zone, k, j, i) = sol_hom1(0,i) ;
235 solution_hom_deux.set(zone, k, j, i) = sol_hom2(1,i) ;
236 solution_hom_trois.set(zone, k, j, i) = sol_hom2(0,i) ;
237 solution_hom_quatre.set(zone, k, j, i) = sol_hom1(1,i) ;
238 }
239 }
240 }
241 }
242 }
243
244 // Compactified external domain
245 //-----------------------------
246 nr = mg.get_nr(nzm1) ;
247 assert(nt == mg.get_nt(nzm1)) ;
248 assert(np == mg.get_np(nzm1)) ;
249 alpha = mpaff->get_alpha()[nzm1] ;
250 assert (nr > 4) ;
251 int nr0 = nr - 2 ; //two degrees of freedom less because of division by u^2
252
253 // Loop on l and m
254 //----------------
255 for (int k=0 ; k<np+1 ; k++) {
256 for (int j=0 ; j<nt ; j++) {
257 base.give_quant_numbers(nzm1, k, j, m_q, l_q, base_r) ;
258 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q != 0) ) {
259 int dege1 = 1; //degeneracy of eq.1
260 int dege2 = 0; //degeneracy of eq.2
261 int nr_eq1 = nr0 - dege1 ; //Eq.1 is for eta
262 int nr_eq2 = nr0 - dege2 ; //Eq.2 is the div-free condition
263 int nrtot = nr_eq1 + nr_eq2 ;
264 Matrice oper(nrtot, nrtot) ; oper.set_etat_qcq() ;
265 Tbl sec_membre(nrtot) ; sec_membre.set_etat_qcq() ;
266 Diff_x2dsdx2 x2d2(base_r, nr) ; const Matrice& m2d2 = x2d2.get_matrice() ;
267 Diff_xdsdx xd(base_r, nr) ; const Matrice& mxd = xd.get_matrice() ;
268 Diff_id id(base_r, nr) ; const Matrice& mid = id.get_matrice() ;
269
270 // Building the operator
271 //----------------------
272 for (int lin=0; lin<nr_eq1; lin++) {
273 for (int col=dege1; col<nr0; col++)
274 oper.set(lin,col-dege1)
275 = m2d2(lin,col) + 4*mxd(lin,col)
276 + (2-(lam+1)*l_q*(l_q+1))*mid(lin,col) ;
277 for (int col=dege2; col<nr0; col++)
278 oper.set(lin,col-dege2+nr_eq1) =
279 -lam*mxd(lin,col) + 2*mid(lin,col) ;
280 }
281 for (int lin=0; lin<nr_eq2; lin++) {
282 for (int col=dege1; col<nr0; col++)
283 oper.set(lin+nr_eq1,col-dege1)
284 = l_q*(l_q+1)*(lam*mxd(lin,col) + (3*lam+2)*mid(lin,col)) ;
285 for (int col=dege2; col<nr0; col++)
287 = (lam+1)*(m2d2(lin,col) + 4*mxd(lin,col))
288 - l_q*(l_q+1)*mid(lin,col) ;
289 }
290 oper.set_lu() ;
291
292 // Filling the r.h.s
293 //------------------
294 for (int i=0; i<nr_eq1; i++)
295 sec_membre.set(i) = (*S_eta.get_spectral_va().c_cf)(nzm1, k, j, i) ;
296 for (int i=0; i<nr_eq2; i++)
297 sec_membre.set(i+nr_eq1) =(*S_r.get_spectral_va().c_cf)(nzm1, k, j, i);
298 Tbl big_res = oper.inverse(sec_membre) ;
299 Tbl res_eta(nr) ; res_eta.set_etat_qcq() ;
300 Tbl res_vr(nr) ; res_vr.set_etat_qcq() ;
301
302 // Putting coefficients of h and v to individual arrays
303 //-----------------------------------------------------
304 for (int i=0; i<dege1; i++)
305 res_eta.set(i) = 0 ;
306 for (int i=dege1; i<nr0; i++)
307 res_eta.set(i) = big_res(i-dege1) ;
308 res_eta.set(nr0) = 0 ;
309 res_eta.set(nr0+1) = 0 ;
310 for (int i=0; i<dege2; i++)
311 res_vr.set(i) = 0 ;
312 for (int i=dege2; i<nr0; i++)
313 res_vr.set(i) = big_res(i-dege2+nr_eq1) ;
314 res_vr.set(nr0) = 0 ;
315 res_vr.set(nr0+1) = 0 ;
316
317 // Multiplication by u^2
318 //-----------------------
319 Tbl res_v2(nr) ; res_v2.set_etat_qcq() ;
320 Tbl res_e2(nr) ; res_e2.set_etat_qcq() ;
321 mult2_xm1_1d_cheb(nr, res_eta.t, res_e2.t) ;
322 mult2_xm1_1d_cheb(nr, res_vr.t, res_v2.t) ;
323
324 // Homogeneous solution (only 1/r^(l+2) and 1/r^l in the CED)
325 Tbl sol_hom1 = solh(nr, l_q-1, 0., base_r) ;
326 Tbl sol_hom2 = solh(nr, l_q+1, 0., base_r) ;
327 for (int i=0 ; i<nr ; i++) {
328 sol_part_eta.set(nzm1, k, j, i) = alpha*alpha*res_e2(i) ;
329 sol_part_vr.set(nzm1, k, j, i) = alpha*alpha*res_v2(i) ;
330 solution_hom_un.set(nzm1, k, j, i) = 0. ;
331 solution_hom_deux.set(nzm1, k, j, i) = sol_hom2(i) ;
332 solution_hom_trois.set(nzm1, k, j, i) = 0. ;
333 solution_hom_quatre.set(nzm1, k, j, i) = sol_hom1(i) ;
334 }
335 }
336 }
337 }
338
339 // Now let's match everything ...
340 //-------------------------------
341
342 // Resulting V^r & eta
343 Scalar vr(*mpaff) ; vr.set_etat_qcq() ;
344 vr.set_spectral_base(base) ;
345 vr.set_spectral_va().set_etat_cf_qcq() ;
346 Mtbl_cf& cf_vr = *vr.set_spectral_va().c_cf ;
347 cf_vr.annule_hard() ;
348 Scalar het(*mpaff) ; het.set_etat_qcq() ;
349 het.set_spectral_base(base) ;
350 het.set_spectral_va().set_etat_cf_qcq() ;
351 Mtbl_cf& cf_eta = *het.set_spectral_va().c_cf ;
352 cf_eta.annule_hard() ;
353 int taille = 4*(nzm1-num_front)-2 ; //## a verifier
354 Tbl sec_membre(taille) ;
355 Matrice systeme(taille, taille) ;
356 systeme.set_etat_qcq() ;
357 int ligne ; int colonne ;
358
359 // Loop on l and m
360 //----------------
361 for (int k=0 ; k<np+1 ; k++)
362 for (int j=0 ; j<nt ; j++) {
363 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
364 if ((nullite_plm(j, nt, k, np, base) == 1)&&(l_q != 0)) {
365
366 double f3_eta = lam*l_q + 3*lam + 2 ;
367 double f4_eta = 2 + 2*lam - lam*l_q ;
368 double f3_vr = (l_q+1)*(lam*l_q - 2) ;
369 double f4_vr = l_q*(lam*l_q + lam + 2) ;
370 ligne = 0 ;
371 colonne = 0 ;
372 sec_membre.annule_hard() ;
373 for (int l=0; l<taille; l++)
374 for (int c=0; c<taille; c++)
375 systeme.set(l,c) = 0 ;
376
377 // First shell
378 nr = mg.get_nr(num_front+1) ;
379 alpha = mpaff->get_alpha()[num_front+1] ;
380 double echelle = mpaff->get_beta()[num_front+1]/alpha ;
381 // Conditions on eta (configuration space)
382 //value and derivative of (x+echelle)^(l-1) at -1
383 systeme.set(ligne, colonne) = pow(echelle-1., double(l_q-1)) ;
384
385 // value and derivative of 1/(x+echelle) ^(l+2) at -1
386 systeme.set(ligne, colonne+1) = 1/pow(echelle-1., double(l_q+2)) ;
387
388 //value and derivative of (x+echelle)^(l+1) at -1
389 systeme.set(ligne, colonne+2) = f3_eta*pow(echelle-1., double(l_q+1)) ;
390 // value and derivative of 1/(x+echelle) ^l at -1
391 systeme.set(ligne, colonne+3) = f4_eta/pow(echelle-1., double(l_q)) ;
392 for (int i=0 ; i<nr ; i++)
393 if (i%2 == 0)
394 sec_membre.set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
395 else sec_membre.set(ligne) += sol_part_eta(num_front+1, k, j, i) ;
396 sec_membre.set(ligne) += bound_eta(num_front+1, k, j, 0) ;
397 ligne++ ;
398
399 // ... and their couterparts for V^r
400 systeme.set(ligne, colonne) = fact_dir*l_q*pow(echelle-1., double(l_q-1)) + fact_neu*l_q*(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
401 systeme.set(ligne, colonne+1) = -fact_dir*(l_q+1)/pow(echelle-1., double(l_q+2)) + fact_neu*(l_q+1)*(l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
402 systeme.set(ligne, colonne+2) = fact_dir*f3_vr*pow(echelle-1., double(l_q+1)) + fact_neu*f3_vr*(l_q+1)*pow(echelle-1., double(l_q))/alpha ;
403 systeme.set(ligne, colonne+3) = fact_dir*f4_vr/pow(echelle-1., double(l_q)) - fact_neu*(f4_vr*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
404 for (int i=0 ; i<nr ; i++)
405 if (i%2 == 0)
407 else sec_membre.set(ligne) += fact_dir*sol_part_vr(num_front+1, k, j, i) - fact_neu*i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
408 sec_membre.set(ligne) += bound_vr(num_front+1, k, j, 0) ;
409
410 ligne++ ;
411
412
413 // Values at 1
414 // eta
415 //value of (x+echelle)^(l-1) at 1
416 systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
417 // value of 1/(x+echelle) ^(l+2) at 1
418 systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
419 //value of (x+echelle)^(l+1) at 1
420 systeme.set(ligne, colonne+2) = f3_eta*pow(echelle+1., double(l_q+1));
421 // value of 1/(x+echelle) ^l at 1
422 systeme.set(ligne, colonne+3) = f4_eta/pow(echelle+1., double(l_q)) ;
423 for (int i=0 ; i<nr ; i++)
424 sec_membre.set(ligne) -= sol_part_eta(num_front+1, k, j, i) ;
425 ligne++ ;
426 // ... and their couterparts for V^r
427 systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
428 systeme.set(ligne, colonne+1)
429 = -double(l_q+1) / pow(echelle+1., double(l_q+2));
430 systeme.set(ligne, colonne+2) = f3_vr*pow(echelle+1., double(l_q+1)) ;
431 systeme.set(ligne, colonne+3) = f4_vr/pow(echelle+1., double(l_q));
432 for (int i=0 ; i<nr ; i++)
433 sec_membre.set(ligne) -= sol_part_vr(num_front+1, k, j, i) ;
434 ligne++ ;
435
436 //Derivatives at 1
437 // eta
438 //derivative of (x+echelle)^(l-1) at 1
439 systeme.set(ligne, colonne)
440 = (l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
441 // derivative of 1/(x+echelle) ^(l+2) at 1
442 systeme.set(ligne, colonne+1)
443 = -(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
444 // derivative of (x+echelle)^(l+1) at 1
445 systeme.set(ligne, colonne+2)
446 = f3_eta*(l_q+1) * pow(echelle+1., double(l_q))/alpha;
447 // derivative of 1/(x+echelle) ^l at 1
448 systeme.set(ligne, colonne+3)
449 = -f4_eta*l_q / pow(echelle+1., double(l_q+1))/alpha ;
450 for (int i=0 ; i<nr ; i++)
451 sec_membre.set(ligne) -= i*i/alpha*sol_part_eta(num_front+1, k, j, i) ;
452 ligne++ ;
453 // ... and their couterparts for V^r
454 systeme.set(ligne, colonne)
455 = l_q*(l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
456 systeme.set(ligne, colonne+1)
457 = (l_q+1)*(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
458 systeme.set(ligne, colonne+2)
459 = f3_vr*(l_q+1) * pow(echelle+1., double(l_q))/alpha ;
460 systeme.set(ligne, colonne+3)
461 = -f4_vr*l_q / pow(echelle+1., double(l_q+1))/alpha ;
462 for (int i=0 ; i<nr ; i++)
463 sec_membre.set(ligne) -= i*i/alpha*sol_part_vr(num_front+1, k, j, i) ;
464
465 colonne += 4 ; // We pass to the next domain
466
467
468 // Next shells
469 if (num_front+2<nzm1){
470 for (int zone=num_front+2 ; zone<nzm1 ; zone++) {
471 nr = mg.get_nr(zone) ;
472 alpha = mpaff->get_alpha()[zone] ;
473 echelle = mpaff->get_beta()[zone]/alpha ;
474 ligne -= 3 ;
475 //value of (x+echelle)^(l-1) at -1
476 systeme.set(ligne, colonne) = -pow(echelle-1., double(l_q-1)) ;
477 // value of 1/(x+echelle) ^(l+2) at -1
478 systeme.set(ligne, colonne+1) = -1/pow(echelle-1., double(l_q+2)) ;
479 //value of (x+echelle)^(l+1) at -1
480 systeme.set(ligne, colonne+2) = -f3_eta*pow(echelle-1., double(l_q+1));
481 // value of 1/(x+echelle) ^l at -1
482 systeme.set(ligne, colonne+3) = -f4_eta/pow(echelle-1., double(l_q)) ;
483 for (int i=0 ; i<nr ; i++)
484 if (i%2 == 0)
485 sec_membre.set(ligne) += sol_part_eta(zone, k, j, i) ;
486 else sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
487 ligne++ ;
488 // ... and their couterparts for V^r
489 systeme.set(ligne, colonne) = -l_q*pow(echelle-1., double(l_q-1)) ;
490 systeme.set(ligne, colonne+1) = (l_q+1)/pow(echelle-1., double(l_q+2));
491 systeme.set(ligne, colonne+2) = -f3_vr*pow(echelle-1., double(l_q+1)) ;
492 systeme.set(ligne, colonne+3) = -f4_vr/pow(echelle-1., double(l_q));
493 for (int i=0 ; i<nr ; i++)
494 if (i%2 == 0)
495 sec_membre.set(ligne) += sol_part_vr(zone, k, j, i) ;
496 else sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
497 ligne++ ;
498
499 //derivative of (x+echelle)^(l-1) at -1
500 systeme.set(ligne, colonne)
501 = -(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
502 // derivative of 1/(x+echelle) ^(l+2) at -1
503 systeme.set(ligne, colonne+1)
504 = (l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
505 // derivative of (x+echelle)^(l+1) at -1
506 systeme.set(ligne, colonne+2)
507 = -f3_eta*(l_q+1)*pow(echelle-1., double(l_q))/alpha;
508 // derivative of 1/(x+echelle) ^l at -1
509 systeme.set(ligne, colonne+3)
510 = (f4_eta*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
511 for (int i=0 ; i<nr ; i++)
512 if (i%2 == 0) sec_membre.set(ligne)
513 -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
514 else sec_membre.set(ligne) +=
515 i*i/alpha*sol_part_eta(zone, k, j, i) ;
516 ligne++ ;
517 // ... and their couterparts for V^r
518 systeme.set(ligne, colonne)
519 = -l_q*(l_q-1)*pow(echelle-1., double(l_q-2))/alpha ;
520 systeme.set(ligne, colonne+1)
521 = -(l_q+1)*(l_q+2)/pow(echelle-1., double(l_q+3))/alpha ;
522 systeme.set(ligne, colonne+2)
523 = -f3_vr*(l_q+1)*pow(echelle-1., double(l_q))/alpha ;
524 systeme.set(ligne, colonne+3)
525 = (f4_vr*l_q/pow(echelle-1., double(l_q+1)))/alpha ;
526 for (int i=0 ; i<nr ; i++)
527 if (i%2 == 0) sec_membre.set(ligne)
528 -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
529 else sec_membre.set(ligne) +=
530 i*i/alpha*sol_part_vr(zone, k, j, i) ;
531 ligne++ ;
532
533 //value of (x+echelle)^(l-1) at 1
534 systeme.set(ligne, colonne) = pow(echelle+1., double(l_q-1)) ;
535 // value of 1/(x+echelle) ^(l+2) at 1
536 systeme.set(ligne, colonne+1) = 1./pow(echelle+1., double(l_q+2)) ;
537 //value of (x+echelle)^(l+1) at 1
538 systeme.set(ligne, colonne+2) = f3_eta*pow(echelle+1., double(l_q+1));
539 // value of 1/(x+echelle) ^l at 1
540 systeme.set(ligne, colonne+3) = f4_eta/pow(echelle+1., double(l_q)) ;
541 for (int i=0 ; i<nr ; i++)
542 sec_membre.set(ligne) -= sol_part_eta(zone, k, j, i) ;
543 ligne++ ;
544 // ... and their couterparts for V^r
545 systeme.set(ligne, colonne) = l_q*pow(echelle+1., double(l_q-1)) ;
546 systeme.set(ligne, colonne+1)
547 = -double(l_q+1) / pow(echelle+1., double(l_q+2));
548 systeme.set(ligne, colonne+2) = f3_vr*pow(echelle+1., double(l_q+1)) ;
549 systeme.set(ligne, colonne+3) = f4_vr/pow(echelle+1., double(l_q));
550 for (int i=0 ; i<nr ; i++)
551 sec_membre.set(ligne) -= sol_part_vr(zone, k, j, i) ;
552 ligne++ ;
553
554 //derivative of (x+echelle)^(l-1) at 1
555 systeme.set(ligne, colonne)
556 = (l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
557 // derivative of 1/(x+echelle) ^(l+2) at 1
558 systeme.set(ligne, colonne+1)
559 = -(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
560 // derivative of (x+echelle)^(l+1) at 1
561 systeme.set(ligne, colonne+2)
562 = f3_eta*(l_q+1) * pow(echelle+1., double(l_q))/alpha;
563 // derivative of 1/(x+echelle) ^l at 1
564 systeme.set(ligne, colonne+3)
565 = -f4_eta*l_q / pow(echelle+1., double(l_q+1))/alpha ;
566 for (int i=0 ; i<nr ; i++)
567 sec_membre.set(ligne) -= i*i/alpha*sol_part_eta(zone, k, j, i) ;
568 ligne++ ;
569 // ... and their couterparts for V^r
570 systeme.set(ligne, colonne)
571 = l_q*(l_q-1) * pow(echelle+1., double(l_q-2))/alpha ;
572 systeme.set(ligne, colonne+1)
573 = (l_q+1)*(l_q+2) / pow(echelle+1., double(l_q+3))/alpha ;
574 systeme.set(ligne, colonne+2)
575 = f3_vr*(l_q+1) * pow(echelle+1., double(l_q))/alpha ;
576 systeme.set(ligne, colonne+3)
577 = -f4_vr*l_q / pow(echelle+1., double(l_q+1))/alpha ;
578 for (int i=0 ; i<nr ; i++)
579 sec_membre.set(ligne) -= i*i/alpha*sol_part_vr(zone, k, j, i) ;
580
581 colonne += 4 ;
582 }
583 }
584 //Compactified external domain
585 nr = mg.get_nr(nzm1) ;
586
587 alpha = mpaff->get_alpha()[nzm1] ;
588 ligne -= 3 ;
589 //value of (x-1)^(l+2) at -1 :
590 systeme.set(ligne, colonne) = -pow(-2, double(l_q+2)) ;
591 //value of (x-1)^l at -1 :
592 systeme.set(ligne, colonne+1) = -f4_eta*pow(-2, double(l_q)) ;
593 for (int i=0 ; i<nr ; i++)
594 if (i%2 == 0) sec_membre.set(ligne) += sol_part_eta(nzm1, k, j, i) ;
595 else sec_membre.set(ligne) -= sol_part_eta(nzm1, k, j, i) ;
596 //... and of its couterpart for V^r
597 systeme.set(ligne+1, colonne) = double(l_q+1)*pow(-2, double(l_q+2)) ;
598 systeme.set(ligne+1, colonne+1) = -f4_vr*pow(-2, double(l_q)) ;
599 for (int i=0 ; i<nr ; i++)
600 if (i%2 == 0) sec_membre.set(ligne+1) += sol_part_vr(nzm1, k, j, i) ;
601 else sec_membre.set(ligne+1) -= sol_part_vr(nzm1, k, j, i) ;
602
603 ligne += 2 ;
604 //derivative of (x-1)^(l+2) at -1 :
605 systeme.set(ligne, colonne) = alpha*(l_q+2)*pow(-2, double(l_q+3)) ;
606 //derivative of (x-1)^l at -1 :
607 systeme.set(ligne, colonne+1) = alpha*l_q*f4_eta*pow(-2, double(l_q+1)) ;
608 for (int i=0 ; i<nr ; i++)
609 if (i%2 == 0) sec_membre.set(ligne)
610 -= -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
611 else sec_membre.set(ligne)
612 += -4*alpha*i*i*sol_part_eta(nzm1, k, j, i) ;
613 //... and of its couterpart for V^r
614 systeme.set(ligne+1, colonne)
615 = -alpha*double((l_q+1)*(l_q+2))*pow(-2, double(l_q+3)) ;
616 systeme.set(ligne+1, colonne+1)
617 = alpha*double(l_q)*f4_vr*pow(-2, double(l_q+1)) ;
618 for (int i=0 ; i<nr ; i++)
619 if (i%2 == 0) sec_membre.set(ligne+1)
620 -= -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
621 else sec_membre.set(ligne+1)
622 += -4*alpha*i*i*sol_part_vr(nzm1, k, j, i) ;
623
624 // Solution of the system giving the coefficients for the homogeneous
625 // solutions
626 //-------------------------------------------------------------------
627 if (taille > 2) systeme.set_band(5,5) ;
628 else systeme.set_band(1,1) ;
629 systeme.set_lu() ;
630 Tbl facteurs(systeme.inverse(sec_membre)) ;
631 int conte = 0 ;
632
633 // everything is put to the right place, the same combination of hom.
634 // solutions (with some l or -(l+1) factors) must be used for V^r
635 //-------------------------------------------------------------------
636
637 for (int zone=1 ; zone<nzm1 ; zone++) { //shells
638 nr = mg.get_nr(zone) ;
639 for (int i=0 ; i<nr ; i++) {
640 cf_eta.set(zone, k, j, i) =
641 sol_part_eta(zone, k, j, i)
646 cf_vr.set(zone, k, j, i) = sol_part_vr(zone, k, j, i)
648 -double(l_q+1)*facteurs(conte+1)*solution_hom_deux(zone, k, j, i) // What is the origin of these factors?!
651 }
652 conte+=4 ;
653 }
654 nr = mg.get_nr(nz-1) ; //compactified external domain
655 for (int i=0 ; i<nr ; i++) {
656 cf_eta.set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
659 cf_vr.set(nzm1, k, j, i) = sol_part_vr(nzm1, k, j, i)
662
663 }
664 } // End of nullite_plm
665 } //End of loop on theta
666 vr.set_spectral_va().ylm_i() ;
667 vr += vrl0 ;
668 het.set_spectral_va().ylm_i() ;
669
670 Valeur temp_mu(mg.get_angu()) ;
671 temp_mu = bound_mu ;
672 const Valeur& limit_mu (temp_mu) ;
673
674 resu.set_vr_eta_mu(vr, 0*het, mu().poisson_dirichlet(limit_mu, num_front)) ;
675
676 return ;
677}
678
679
680Vector Vector::poisson_dirichlet(double lam, const Valeur& bound_vr,
681 const Valeur& bound_vt, const Valeur& bound_vp,
682 int num_front) const {
683
684 Vector resu(*mp, CON, triad) ;
686
687 return resu ;
688
689}
690
692 const Valeur& bound_vt, const Valeur& bound_vp,
693 int num_front) const {
694
695 Vector resu(*mp, CON, triad) ;
697
698 return resu ;
699
700}
701
703 const Valeur& bound_vt, const Valeur& bound_vp,
704 double fact_dir, double fact_neu,
705 int num_front) const {
706
707
708 // Boundary condition for V^r //Construction of a Mtbl_cf from Valeur with Ylm coefficients
710
711 limit_vr.coef() ;
712 limit_vr.ylm() ; // Spherical harmonics transform.
713 Mtbl_cf lim_vr (*(limit_vr.c_cf)) ;
714
715 // bound_vt and bound_vp are only known at the boundary --> we fill
716 // all the zones extending the values at the boundary before calling to poisson_angu.
717 Scalar temp_vt (*mp) ;
718 Scalar temp_vp (*mp) ;
719 temp_vt.annule_hard() ;
720 temp_vp.annule_hard() ;
721 int nz = mp->get_mg()->get_nzone() ;
722 for (int l=0; l<nz; l++)
723 for (int j=0; j<mp->get_mg()->get_nt(l); j++)
724 for (int k=0; k<mp->get_mg()->get_np(l); k++) {
725 temp_vt.set_grid_point(l, k, j, 0) = bound_vt(1, k, j, 0) ;
726 temp_vp.set_grid_point(l, k, j, 0) = bound_vp(1, k, j, 0) ;
727 }
728 temp_vt.set_spectral_va().set_base(bound_vt.base) ; // We set the basis
729 temp_vp.set_spectral_va().set_base(bound_vp.base) ;
730
731 cout << "temp_vp" << endl << norme(temp_vp) << endl ;
732
733 //Source for eta
736 vtstant.div_tant() ;
737 source_eta = temp_vt.dsdt() + vtstant + temp_vp.stdsdp() ;
738
739 //Source for mu
742 vpstant.div_tant() ;
743 source_mu = temp_vp.dsdt() + vpstant - temp_vt.stdsdp() ; //There was a wrong sign here
744
745 Scalar temp_mu (source_mu.poisson_angu()) ;
746 Scalar temp_eta (source_eta.poisson_angu()) ;
747
748 // Boundary condition for mu
749 Valeur limit_mu ((*mp).get_mg()->get_angu() ) ;
750 int nnp = (*mp).get_mg()->get_np(1) ;
751 int nnt = (*mp).get_mg()->get_nt(1) ;
752 limit_mu= 1 ;
753 for (int k=0 ; k<nnp ; k++)
754 for (int j=0 ; j<nnt ; j++)
755 limit_mu.set(1, k, j, 0) = temp_mu.val_grid_point(1, k, j, 0) ;
756 limit_mu.set_base(temp_mu.get_spectral_va().get_base()) ;
757
758 limit_mu.coef() ;
759 limit_mu.ylm() ; // Spherical harmonics transform.
760 Mtbl_cf lim_mu (*(limit_mu.c_cf)) ;
761
762 // Boundary condition for eta
763 Valeur limit_eta ((*mp).get_mg()->get_angu() ) ;
764 limit_eta = 1 ;
765 for (int k=0 ; k<nnp ; k++)
766 for (int j=0 ; j<nnt ; j++)
767 limit_eta.set(1, k, j, 0) = temp_eta.val_grid_point(1, k, j, 0) ;
768 limit_eta.set_base(temp_eta.get_spectral_va().get_base()) ;
769
770 limit_eta.coef() ;
771 limit_eta.ylm() ; // Spherical harmonics transform.
772 Mtbl_cf lim_eta (*(limit_eta.c_cf)) ;
773
774
775 // Call to poisson_boundary(...)
776 bool nullite = true ;
777 for (int i=0; i<3; i++) {
778 assert(cmp[i]->check_dzpuis(4)) ;
779 if (cmp[i]->get_etat() != ETATZERO || bound_vr.get_etat() != ETATZERO ||
780 bound_vt.get_etat() != ETATZERO || bound_vp.get_etat() != ETATZERO)
781 nullite = false ;
782 }
783
784 Vector resu(*mp, CON, triad) ;
785 if (nullite)
786 resu.set_etat_zero() ;
787 else
789 fact_neu, resu) ;
790
791 return resu ;
792
793}
794}
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 (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
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tensor field of valence 1.
Definition vector.h:188
void poisson_boundary(double lambda, const Mtbl_cf &limit_vr, const Mtbl_cf &limit_eta, const Mtbl_cf &limit_mu, int num_front, double fact_dir, double fact_neu, Vector &resu) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
Vector poisson_neumann(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
Vector poisson_robin(double lambda, const Valeur &limit_vr, const Valeur &limit_vt, const Valeur &limit_vp, double fact_dir, double fact_neu, int num_front) const
Solves the vector Poisson equation with *this as a source with a boundary condition on the excised sp...
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
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