LORENE
poisson_ylm.C
1/*
2 * Method of Non_class_members for solving Poisson equations
3 * with a multipole falloff condition at the outer boundary
4 *
5 */
6
7/*
8 * Copyright (c) 2004 Joshua A. Faber
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28
29/*
30 * $Id: poisson_ylm.C,v 1.4 2016/12/05 16:18:10 j_novak Exp $
31 * $Log: poisson_ylm.C,v $
32 * Revision 1.4 2016/12/05 16:18:10 j_novak
33 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34 *
35 * Revision 1.3 2014/10/13 08:53:30 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.2 2014/10/06 15:16:09 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.1 2004/12/29 16:32:02 k_taniguchi
42 * *** empty log message ***
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_ylm.C,v 1.4 2016/12/05 16:18:10 j_novak Exp $
46 *
47 */
48
49// C headers
50#include <cstdlib>
51#include <cmath>
52
53// Lorene headers
54#include "matrice.h"
55#include "tbl.h"
56#include "mtbl_cf.h"
57#include "map.h"
58#include "base_val.h"
59#include "proto.h"
60#include "type_parite.h"
61
62
63
64 //----------------------------------------------
65 // Version Mtbl_cf
66 //----------------------------------------------
67
68/*
69 *
70 * Solution de l'equation de poisson
71 *
72 * Entree : mapping : le mapping affine
73 * source : les coefficients de la source qui a ete multipliee par
74 * r^4 ou r^2 dans la ZEC.
75 * La base de decomposition doit etre Ylm
76 * k_ylm: exponent in radial dependence of field: phi \propto r^{-k}
77 * Sortie : renvoie les coefficients de la solution dans la meme base de
78 * decomposition (a savoir Ylm)
79 *
80 */
81
82
83
84namespace Lorene {
85Mtbl_cf sol_poisson_ylm(const Map_af& mapping, const Mtbl_cf& source, const int nylm, const double* intvec)
86{
87
88 // Verifications d'usage sur les zones
89 int nz = source.get_mg()->get_nzone() ;
90 assert (nz>1) ;
91 assert (source.get_mg()->get_type_r(0) == RARE) ;
92 // assert (source.get_mg()->get_type_r(nz-1) == UNSURR) ;
93 for (int l=1 ; l<nz ; l++)
94 assert(source.get_mg()->get_type_r(l) == FIN) ;
95
96 // assert ((dzpuis==4) || (dzpuis==2) || (dzpuis==3)) ;
97
98
99 // Bases spectrales
100 const Base_val& base = source.base ;
101
102 // donnees sur la zone
103 int nr, nt, np ;
104 int base_r ;
105 double alpha, beta, echelle ;
106 int l_quant, m_quant;
107
108 //Rangement des valeurs intermediaires
109 Tbl *so ;
110 Tbl *sol_hom ;
111 Tbl *sol_part ;
112 Matrice *operateur ;
113 Matrice *nondege ;
114
115
116 // Rangement des solutions, avant raccordement
117 Mtbl_cf solution_part(source.get_mg(), base) ;
118 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
119 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
120 Mtbl_cf resultat(source.get_mg(), base) ;
121
122 solution_part.set_etat_qcq() ;
123 solution_hom_un.set_etat_qcq() ;
124 solution_hom_deux.set_etat_qcq() ;
125 resultat.set_etat_qcq() ;
126
127 for (int l=0 ; l<nz ; l++) {
128 solution_part.t[l]->set_etat_qcq() ;
129 solution_hom_un.t[l]->set_etat_qcq() ;
130 solution_hom_deux.t[l]->set_etat_qcq() ;
131 resultat.t[l]->set_etat_qcq() ;
132 for (int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
133 for (int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
134 for (int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
135 resultat.set(l, k, j, i) = 0 ;
136 }
137
138 // nbre maximum de point en theta et en phi :
139 int np_max, nt_max ;
140
141 //---------------
142 //-- NOYAU -----
143 //---------------
144
145 nr = source.get_mg()->get_nr(0) ;
146 nt = source.get_mg()->get_nt(0) ;
147 np = source.get_mg()->get_np(0) ;
148
149 nt_max = nt ;
150 np_max = np ;
151
152 alpha = mapping.get_alpha()[0] ;
153 beta = mapping.get_beta()[0] ;
154
155 for (int k=0 ; k<np+1 ; k++)
156 for (int j=0 ; j<nt ; j++)
157 if (nullite_plm(j, nt, k, np, base) == 1)
158 {
159 // calcul des nombres quantiques :
160 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
161
162 // Construction operateur
163 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
164 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
165
166 //Operateur inversible
167 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
168
169 // Calcul de la SH
170 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
171
172 //Calcul de la SP
173 so = new Tbl(nr) ;
174 so->set_etat_qcq() ;
175 for (int i=0 ; i<nr ; i++)
176 so->set(i) = source(0, k, j, i) ;
177
178 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
179 *so, 0, base_r)) ;
180
181 // Rangement dans les tableaux globaux ;
182
183 for (int i=0 ; i<nr ; i++) {
184 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
185 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
186 solution_hom_deux.set(0, k, j, i) = 0. ;
187 }
188
189
190
191 delete operateur ;
192 delete nondege ;
193 delete so ;
194 delete sol_hom ;
195 delete sol_part ;
196 }
197
198
199
200 //---------------
201 //- COQUILLES ---
202 //---------------
203
204 for (int zone=1 ; zone<nz ; zone++) {
205 nr = source.get_mg()->get_nr(zone) ;
206 nt = source.get_mg()->get_nt(zone) ;
207 np = source.get_mg()->get_np(zone) ;
208
209 if (np > np_max) np_max = np ;
210 if (nt > nt_max) nt_max = nt ;
211
212 alpha = mapping.get_alpha()[zone] ;
213 beta = mapping.get_beta()[zone] ;
214
215 for (int k=0 ; k<np+1 ; k++)
216 for (int j=0 ; j<nt ; j++)
217 if (nullite_plm(j, nt, k, np, base) == 1)
218 {
219 // calcul des nombres quantiques :
220 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
221
222 // Construction de l'operateur
223 operateur = new Matrice(laplacien_mat
224 (nr, l_quant, beta/alpha, 0, base_r)) ;
225
226 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
227 base_r) ;
228
229 // Operateur inversible
230 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
231 beta/alpha, 0, base_r)) ;
232
233 // Calcul DES DEUX SH
234 sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
235
236 // Calcul de la SP
237 so = new Tbl(nr) ;
238 so->set_etat_qcq() ;
239 for (int i=0 ; i<nr ; i++)
240 so->set(i) = source(zone, k, j, i) ;
241
242 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
243 beta, *so, 0, base_r)) ;
244
245
246 // Rangement
247 for (int i=0 ; i<nr ; i++) {
248 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
249 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
250 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
251 }
252
253
254 delete operateur ;
255 delete nondege ;
256 delete so ;
257 delete sol_hom ;
258 delete sol_part ;
259 }
260 }
261
262 //-------------------------------------------
263 // On est parti pour le raccord des solutions
264 //-------------------------------------------
265 // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
266 // intervient dans le developpement de cette zone.
267 int * indic = new int [nz] ;
268 int taille ;
269 Tbl *sec_membre ; // termes constants du systeme
270 Matrice *systeme ; //le systeme a resoudre
271
272 // des compteurs pour le remplissage du systeme
273 int ligne ;
274 int colonne ;
275
276 // compteurs pour les diagonales du systeme :
277 int sup ;
278 int inf ;
279 int sup_new, inf_new ;
280
281 // on boucle sur le plus grand nombre possible de Plm intervenant...
282 for (int k=0 ; k<np_max+1 ; k++)
283 for (int j=0 ; j<nt_max ; j++)
284 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
285
286 ligne = 0 ;
287 colonne = 0 ;
288 sup = 0 ;
289 inf = 0 ;
290
291 for (int zone=0 ; zone<nz ; zone++)
292 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
293 k, source.get_mg()->get_np(zone), base);
294
295 // taille du systeme a resoudre pour ce Plm
296 taille = indic[0] ;
297 for (int zone=1 ; zone<nz ; zone++)
298 taille+=2*indic[zone] ;
299
300 // on verifie que la taille est non-nulle.
301 // cas pouvant a priori se produire...
302
303 if (taille != 0) {
304
305 sec_membre = new Tbl(taille) ;
306 sec_membre->set_etat_qcq() ;
307 for (int l=0 ; l<taille ; l++)
308 sec_membre->set(l) = 0 ;
309
310 systeme = new Matrice(taille, taille) ;
311 systeme->set_etat_qcq() ;
312 for (int l=0 ; l<taille ; l++)
313 for (int c=0 ; c<taille ; c++)
314 systeme->set(l, c) = 0 ;
315
316 //Calcul des nombres quantiques
317 //base_r est donne dans le noyau, sa valeur dans les autres
318 //zones etant inutile.
319
320 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
321
322
323 //--------------------------
324 // NOYAU
325 //---------------------------
326
327 if (indic[0] == 1) {
328 nr = source.get_mg()->get_nr(0) ;
329
330 alpha = mapping.get_alpha()[0] ;
331 // valeur de x^l en 1 :
332 systeme->set(ligne, colonne) = 1. ; /* ligne=0, colonne=0 */
333 // coefficient du Plm dans la solp
334 for (int i=0 ; i<nr ; i++)
335 sec_membre->set(ligne) -= solution_part(0, k, j, i) ; /* ligne=0 */
336
337 ligne++ ; /* ligne=1, colonne=0 */
338 // on prend les derivees que si Plm existe
339 //dans la zone suivante
340
341 if (indic[1] == 1) {
342 // derivee de x^l en 1 :
343 systeme->set(ligne, colonne) = 1./alpha*l_quant ; /* ligne=1, colonne=0 */
344
345 // coefficient de la derivee du Plm dans la solp
346 if (base_r == R_CHEBP)
347 // cas en R_CHEBP
348 for (int i=0 ; i<nr ; i++)
349 sec_membre->set(ligne) -=
350 4*i*i/alpha
351 *solution_part(0, k, j, i) ; /* ligne=1 */
352 else
353 // cas en R_CHEBI
354 for (int i=0 ; i<nr ; i++)
355 sec_membre->set(ligne) -=
356 (2*i+1)*(2*i+1)/alpha
357 *solution_part(0, k, j, i) ; /* ligne=1 */
358
359 // on a au moins un diag inferieure dans ce cas ...
360 inf = 1 ;
361 }
362 colonne++ ; /* ligne=1, colonne=1 */
363 }
364
365 //-----------------------------
366 // COQUILLES
367 //------------------------------
368
369 for (int zone=1 ; zone<nz ; zone++)
370 if (indic[zone] == 1) {
371
372 nr = source.get_mg()->get_nr(zone) ;
373 alpha = mapping.get_alpha()[zone] ;
374 echelle = mapping.get_beta()[zone]/alpha ;
375
376 //Frontiere avec la zone precedente :
377 if (indic[zone-1] == 1) ligne -- ; /* ligne=0, colonne=1 */
378
379 //valeur de (x+echelle)^l en -1 :
380 systeme->set(ligne, colonne) =
381 -pow(echelle-1, double(l_quant)) ; /* ligne=0, colonne=1 */
382
383 // valeur de 1/(x+echelle) ^(l+1) en -1 :
384 systeme->set(ligne, colonne+1) =
385 -1/pow(echelle-1, double(l_quant+1)) ; /* ligne=0, colonne=1, colonne+1=2 */
386
387 // la solution particuliere :
388 for (int i=0 ; i<nr ; i++)
389 if (i%2 == 0)
390 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
391 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=0 */
392
393 // les diagonales :
394 sup_new = colonne+1-ligne ; /* ligne=0, colonne=1, colonne+1-ligne=2 */
395 if (sup_new > sup) sup = sup_new ;
396
397 ligne++ ; /* ligne=1 */
398
399 // on prend les derivees si Plm existe dans la zone
400 // precedente :
401
402 if (indic[zone-1] == 1) {
403 // derivee de (x+echelle)^l en -1 :
404 systeme->set(ligne, colonne) =
405 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ; /* ligne=1, colonne=1 */
406 // derivee de 1/(x+echelle)^(l+1) en -1 :
407 systeme->set(ligne, colonne+1) =
408 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ; /* ligne=1, colonne=1, colonne+1=2 */
409
410
411
412 // la solution particuliere :
413 for (int i=0 ; i<nr ; i++)
414 if (i%2 == 0)
415 sec_membre->set(ligne) -=
416 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
417 else
418 sec_membre->set(ligne) +=
419 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
420
421 // les diagonales :
422 sup_new = colonne+1-ligne ; /* ligne=1, colonne=1, colonne+1-ligne=1 */
423 if (sup_new > sup) sup = sup_new ;
424
425 ligne++ ; /* ligne=2 */
426 }
427
428
429 if(zone < nz-1) {
430
431 // Frontiere avec la zone suivante :
432 //valeur de (x+echelle)^l en 1 :
433 systeme->set(ligne, colonne) =
434 pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
435
436 // valeur de 1/(x+echelle)^(l+1) en 1 :
437 systeme->set(ligne, colonne+1) =
438 1/pow(echelle+1, double(l_quant+1)) ; /* ligne=2, colonne=1, colonne+1=2 */
439
440 // la solution particuliere :
441 for (int i=0 ; i<nr ; i++)
442 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=2 */
443
444 // les diagonales inf :
445 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
446 if (inf_new > inf) inf = inf_new ;
447
448 ligne ++ ; /* ligne=3 */
449
450 // Utilisation des derivees ssi Plm existe dans la
451 //zone suivante :
452 if (indic[zone+1] == 1) {
453
454 //derivee de (x+echelle)^l en 1 :
455 systeme->set(ligne, colonne) =
456 l_quant*pow(echelle+1, double(l_quant-1))/alpha ; /* ligne=3, colonne=1 */
457
458 //derivee de 1/(echelle+x)^(l+1) en 1 :
459 systeme->set(ligne, colonne+1) =
460 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ; /* ligne=3, colonne=1, colonne+1=2 */
461
462 // La solution particuliere :
463 for (int i=0 ; i<nr ; i++)
464 sec_membre->set(ligne) -=
465 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=3 */
466
467 // les diagonales inf :
468 inf_new = ligne-colonne ; /* ligne=3, colonne=1, ligne-colonne=2 */
469 if (inf_new > inf) inf = inf_new ;
470
471 }
472 colonne += 2 ; /* ligne=colonne=3 -> ligne=colonne=1 next block of two */
473 } else {
474
475
476
477 //-------------------------
478 // Ylm outer boundary
479 //---------------------------
480
481 /* ligne=2, colonne=1 */
482
483 // The coefficient for A_n is (k+l)(echelle+1)^l
484 systeme->set(ligne, colonne) =
485 pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
486
487 // The coefficient for B_n is (k-(l+1))(echelle+1)^{-(l+1)}
488 systeme->set(ligne, colonne+1) =
489 pow(echelle+1, double(-1-l_quant)) ; /* ligne=2, colonne=1, colonne+1=2 */
490
491 // Here we have -f - multipole integral
492 //(pos source->neg field!!)
493
494 int indylm;
495 double intterm=0.;
496 if(l_quant*l_quant < nylm && m_quant<=l_quant) {
497 indylm=l_quant*l_quant+2*m_quant;
498 if(k%2==0 && k>=2)indylm-=1;
499 intterm=intvec[indylm];
500 cout <<"l,m:"<<l_quant<<" "<<m_quant<<" "<<j<<" "<<k<<" "<<indylm<<" "<<intterm<<endl;
501 }
502 // This is just for testing!!!
503 //if(l_quant==1 && m_quant==1&& k%2==0) {
504 // intterm=1.0 ;
505 //} else {
506 // intterm=0.;
507 //}
508
509 sec_membre->set(ligne) -= intterm;
510
511 // La solution particuliere :
512 for (int i=0 ; i<nr ; i++)
513 sec_membre->set(ligne) -=
514 solution_part(zone, k, j, i) ; /* ligne=2 */
515
516 // les diagonales inf :
517 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
518 if (inf_new > inf) inf = inf_new ;
519
520 }
521 }
522
523
524 //-------------------------
525 // resolution du systeme
526 //--------------------------
527
528 systeme->set_band(sup, inf) ;
529 systeme->set_lu() ;
530
531 Tbl facteurs(systeme->inverse(*sec_membre)) ;
532 int conte = 0 ;
533
534
535 // rangement dans le noyau :
536
537 if (indic[0] == 1) {
538 nr = source.get_mg()->get_nr(0) ;
539 for (int i=0 ; i<nr ; i++)
540 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
541 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
542 conte++ ;
543 }
544
545 // rangement dans les coquilles :
546 for (int zone=1 ; zone<nz ; zone++)
547 if (indic[zone] == 1) {
548 nr = source.get_mg()->get_nr(zone) ;
549 for (int i=0 ; i<nr ; i++)
550 resultat.set(zone, k, j, i) =
551 solution_part(zone, k, j, i)
552 +facteurs(conte)*solution_hom_un(zone, k, j, i)
553 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
554 conte+=2 ;
555 }
556
557 delete sec_membre ;
558 delete systeme ;
559 }
560
561 }
562
563 delete [] indic ;
564
565 return resultat;
566}
567}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:325
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Lorene::Map_af
Affine radial mapping.
Definition map.h:2042
Lorene::Matrice
Matrix handling.
Definition matrice.h:152
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
Lorene::Mtbl_cf::get_mg
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:463
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67