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