LORENE
poisson_falloff.C
1/*
2 * Method of Non_class_members for solving Poisson equations
3 * with a 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_falloff.C,v 1.4 2016/12/05 16:18:10 j_novak Exp $
31 * $Log: poisson_falloff.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:29 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/11/30 20:55:03 k_taniguchi
42 * *** empty log message ***
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/poisson_falloff.C,v 1.4 2016/12/05 16:18:10 j_novak Exp $
46 *
47 */
48
49// Header C :
50#include <cstdlib>
51#include <cmath>
52
53// Headers Lorene :
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_falloff: 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_falloff(const Map_af& mapping, const Mtbl_cf& source, const int k_falloff)
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
103 // donnees sur la zone
104 int nr, nt, np ;
105 int base_r ;
106 double alpha, beta, echelle ;
107 int l_quant, m_quant;
108
109 //Rangement des valeurs intermediaires
110 Tbl *so ;
111 Tbl *sol_hom ;
112 Tbl *sol_part ;
113 Matrice *operateur ;
114 Matrice *nondege ;
115
116
117 // Rangement des solutions, avant raccordement
118 Mtbl_cf solution_part(source.get_mg(), base) ;
119 Mtbl_cf solution_hom_un(source.get_mg(), base) ;
120 Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
121 Mtbl_cf resultat(source.get_mg(), base) ;
122
123 solution_part.set_etat_qcq() ;
124 solution_hom_un.set_etat_qcq() ;
125 solution_hom_deux.set_etat_qcq() ;
126 resultat.set_etat_qcq() ;
127
128 for (int l=0 ; l<nz ; l++) {
129 solution_part.t[l]->set_etat_qcq() ;
130 solution_hom_un.t[l]->set_etat_qcq() ;
131 solution_hom_deux.t[l]->set_etat_qcq() ;
132 resultat.t[l]->set_etat_qcq() ;
133 for (int k=0 ; k<source.get_mg()->get_np(l)+2 ; k++)
134 for (int j=0 ; j<source.get_mg()->get_nt(l) ; j++)
135 for (int i=0 ; i<source.get_mg()->get_nr(l) ; i++)
136 resultat.set(l, k, j, i) = 0 ;
137 }
138
139 // nbre maximum de point en theta et en phi :
140 int np_max, nt_max ;
141
142 //---------------
143 //-- NOYAU -----
144 //---------------
145
146 nr = source.get_mg()->get_nr(0) ;
147 nt = source.get_mg()->get_nt(0) ;
148 np = source.get_mg()->get_np(0) ;
149
150 nt_max = nt ;
151 np_max = np ;
152
153 alpha = mapping.get_alpha()[0] ;
154 beta = mapping.get_beta()[0] ;
155
156 for (int k=0 ; k<np+1 ; k++)
157 for (int j=0 ; j<nt ; j++)
158 if (nullite_plm(j, nt, k, np, base) == 1)
159 {
160 // calcul des nombres quantiques :
161 donne_lm(nz, 0, j, k, base, m_quant, l_quant, base_r) ;
162
163 // Construction operateur
164 operateur = new Matrice(laplacien_mat(nr, l_quant, 0., 0, base_r)) ;
165 (*operateur) = combinaison(*operateur, l_quant, 0., 0, base_r) ;
166
167 //Operateur inversible
168 nondege = new Matrice(prepa_nondege(*operateur, l_quant, 0., 0, base_r)) ;
169
170 // Calcul de la SH
171 sol_hom = new Tbl(solh(nr, l_quant, 0., base_r)) ;
172
173 //Calcul de la SP
174 so = new Tbl(nr) ;
175 so->set_etat_qcq() ;
176 for (int i=0 ; i<nr ; i++)
177 so->set(i) = source(0, k, j, i) ;
178
179 sol_part = new Tbl(solp(*operateur, *nondege, alpha, beta,
180 *so, 0, base_r)) ;
181
182 // Rangement dans les tableaux globaux ;
183
184 for (int i=0 ; i<nr ; i++) {
185 solution_part.set(0, k, j, i) = (*sol_part)(i) ;
186 solution_hom_un.set(0, k, j, i) = (*sol_hom)(i) ;
187 solution_hom_deux.set(0, k, j, i) = 0. ;
188 }
189
190
191
192 delete operateur ;
193 delete nondege ;
194 delete so ;
195 delete sol_hom ;
196 delete sol_part ;
197 }
198
199
200
201 //---------------
202 //- COQUILLES ---
203 //---------------
204
205 for (int zone=1 ; zone<nz ; zone++) {
206 nr = source.get_mg()->get_nr(zone) ;
207 nt = source.get_mg()->get_nt(zone) ;
208 np = source.get_mg()->get_np(zone) ;
209
210 if (np > np_max) np_max = np ;
211 if (nt > nt_max) nt_max = nt ;
212
213 alpha = mapping.get_alpha()[zone] ;
214 beta = mapping.get_beta()[zone] ;
215
216 for (int k=0 ; k<np+1 ; k++)
217 for (int j=0 ; j<nt ; j++)
218 if (nullite_plm(j, nt, k, np, base) == 1)
219 {
220 // calcul des nombres quantiques :
221 donne_lm(nz, zone, j, k, base, m_quant, l_quant, base_r) ;
222
223 // Construction de l'operateur
224 operateur = new Matrice(laplacien_mat
225 (nr, l_quant, beta/alpha, 0, base_r)) ;
226
227 (*operateur) = combinaison(*operateur, l_quant, beta/alpha, 0,
228 base_r) ;
229
230 // Operateur inversible
231 nondege = new Matrice(prepa_nondege(*operateur, l_quant,
232 beta/alpha, 0, base_r)) ;
233
234 // Calcul DES DEUX SH
235 sol_hom = new Tbl(solh(nr, l_quant, beta/alpha, base_r)) ;
236
237 // Calcul de la SP
238 so = new Tbl(nr) ;
239 so->set_etat_qcq() ;
240 for (int i=0 ; i<nr ; i++)
241 so->set(i) = source(zone, k, j, i) ;
242
243 sol_part = new Tbl (solp(*operateur, *nondege, alpha,
244 beta, *so, 0, base_r)) ;
245
246
247 // Rangement
248 for (int i=0 ; i<nr ; i++) {
249 solution_part.set(zone, k, j, i) = (*sol_part)(i) ;
250 solution_hom_un.set(zone, k, j, i) = (*sol_hom)(0, i) ;
251 solution_hom_deux.set(zone, k, j, i) = (*sol_hom)(1, i) ;
252 }
253
254
255 delete operateur ;
256 delete nondege ;
257 delete so ;
258 delete sol_hom ;
259 delete sol_part ;
260 }
261 }
262
263 //-------------------------------------------
264 // On est parti pour le raccord des solutions
265 //-------------------------------------------
266 // Tableau de 0 et 1 sur les zones, indiquant si le Plm considere
267 // intervient dans le developpement de cette zone.
268 int * indic = new int [nz] ;
269 int taille ;
270 Tbl *sec_membre ; // termes constants du systeme
271 Matrice *systeme ; //le systeme a resoudre
272
273 // des compteurs pour le remplissage du systeme
274 int ligne ;
275 int colonne ;
276
277 // compteurs pour les diagonales du systeme :
278 int sup ;
279 int inf ;
280 int sup_new, inf_new ;
281
282 // on boucle sur le plus grand nombre possible de Plm intervenant...
283 for (int k=0 ; k<np_max+1 ; k++)
284 for (int j=0 ; j<nt_max ; j++)
285 if (nullite_plm(j, nt_max, k, np_max, base) == 1) {
286
287 ligne = 0 ;
288 colonne = 0 ;
289 sup = 0 ;
290 inf = 0 ;
291
292 for (int zone=0 ; zone<nz ; zone++)
293 indic[zone] = nullite_plm(j, source.get_mg()->get_nt(zone),
294 k, source.get_mg()->get_np(zone), base);
295
296 // taille du systeme a resoudre pour ce Plm
297 taille = indic[0] ;
298 for (int zone=1 ; zone<nz ; zone++)
299 taille+=2*indic[zone] ;
300
301 // on verifie que la taille est non-nulle.
302 // cas pouvant a priori se produire...
303
304 if (taille != 0) {
305
306 sec_membre = new Tbl(taille) ;
307 sec_membre->set_etat_qcq() ;
308 for (int l=0 ; l<taille ; l++)
309 sec_membre->set(l) = 0 ;
310
311 systeme = new Matrice(taille, taille) ;
312 systeme->set_etat_qcq() ;
313 for (int l=0 ; l<taille ; l++)
314 for (int c=0 ; c<taille ; c++)
315 systeme->set(l, c) = 0 ;
316
317 //Calcul des nombres quantiques
318 //base_r est donne dans le noyau, sa valeur dans les autres
319 //zones etant inutile.
320
321 donne_lm (nz, 0, j, k, base, m_quant, l_quant, base_r) ;
322
323
324 //--------------------------
325 // NOYAU
326 //---------------------------
327
328 if (indic[0] == 1) {
329 nr = source.get_mg()->get_nr(0) ;
330
331 alpha = mapping.get_alpha()[0] ;
332 // valeur de x^l en 1 :
333 systeme->set(ligne, colonne) = 1. ; /* ligne=0, colonne=0 */
334 // coefficient du Plm dans la solp
335 for (int i=0 ; i<nr ; i++)
336 sec_membre->set(ligne) -= solution_part(0, k, j, i) ; /* ligne=0 */
337
338 ligne++ ; /* ligne=1, colonne=0 */
339 // on prend les derivees que si Plm existe
340 //dans la zone suivante
341
342 if (indic[1] == 1) {
343 // derivee de x^l en 1 :
344 systeme->set(ligne, colonne) = 1./alpha*l_quant ; /* ligne=1, colonne=0 */
345
346 // coefficient de la derivee du Plm dans la solp
347 if (base_r == R_CHEBP)
348 // cas en R_CHEBP
349 for (int i=0 ; i<nr ; i++)
350 sec_membre->set(ligne) -=
351 4*i*i/alpha
352 *solution_part(0, k, j, i) ; /* ligne=1 */
353 else
354 // cas en R_CHEBI
355 for (int i=0 ; i<nr ; i++)
356 sec_membre->set(ligne) -=
357 (2*i+1)*(2*i+1)/alpha
358 *solution_part(0, k, j, i) ; /* ligne=1 */
359
360 // on a au moins un diag inferieure dans ce cas ...
361 inf = 1 ;
362 }
363 colonne++ ; /* ligne=1, colonne=1 */
364 }
365
366 //-----------------------------
367 // COQUILLES
368 //------------------------------
369
370 for (int zone=1 ; zone<nz ; zone++)
371 if (indic[zone] == 1) {
372
373 nr = source.get_mg()->get_nr(zone) ;
374 alpha = mapping.get_alpha()[zone] ;
375 echelle = mapping.get_beta()[zone]/alpha ;
376
377 //Frontiere avec la zone precedente :
378 if (indic[zone-1] == 1) ligne -- ; /* ligne=0, colonne=1 */
379
380 //valeur de (x+echelle)^l en -1 :
381 systeme->set(ligne, colonne) =
382 -pow(echelle-1, double(l_quant)) ; /* ligne=0, colonne=1 */
383
384 // valeur de 1/(x+echelle) ^(l+1) en -1 :
385 systeme->set(ligne, colonne+1) =
386 -1/pow(echelle-1, double(l_quant+1)) ; /* ligne=0, colonne=1, colonne+1=2 */
387
388 // la solution particuliere :
389 for (int i=0 ; i<nr ; i++)
390 if (i%2 == 0)
391 sec_membre->set(ligne) += solution_part(zone, k, j, i) ;
392 else sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=0 */
393
394 // les diagonales :
395 sup_new = colonne+1-ligne ; /* ligne=0, colonne=1, colonne+1-ligne=2 */
396 if (sup_new > sup) sup = sup_new ;
397
398 ligne++ ; /* ligne=1 */
399
400 // on prend les derivees si Plm existe dans la zone
401 // precedente :
402
403 if (indic[zone-1] == 1) {
404 // derivee de (x+echelle)^l en -1 :
405 systeme->set(ligne, colonne) =
406 -l_quant*pow(echelle-1, double(l_quant-1))/alpha ; /* ligne=1, colonne=1 */
407 // derivee de 1/(x+echelle)^(l+1) en -1 :
408 systeme->set(ligne, colonne+1) =
409 (l_quant+1)/pow(echelle-1, double(l_quant+2))/alpha ; /* ligne=1, colonne=1, colonne+1=2 */
410
411
412
413 // la solution particuliere :
414 for (int i=0 ; i<nr ; i++)
415 if (i%2 == 0)
416 sec_membre->set(ligne) -=
417 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
418 else
419 sec_membre->set(ligne) +=
420 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=1 */
421
422 // les diagonales :
423 sup_new = colonne+1-ligne ; /* ligne=1, colonne=1, colonne+1-ligne=1 */
424 if (sup_new > sup) sup = sup_new ;
425
426 ligne++ ; /* ligne=2 */
427 }
428
429
430 if(zone < nz-1) {
431
432 // Frontiere avec la zone suivante :
433 //valeur de (x+echelle)^l en 1 :
434 systeme->set(ligne, colonne) =
435 pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
436
437 // valeur de 1/(x+echelle)^(l+1) en 1 :
438 systeme->set(ligne, colonne+1) =
439 1/pow(echelle+1, double(l_quant+1)) ; /* ligne=2, colonne=1, colonne+1=2 */
440
441 // la solution particuliere :
442 for (int i=0 ; i<nr ; i++)
443 sec_membre->set(ligne) -= solution_part(zone, k, j, i) ; /* ligne=2 */
444
445 // les diagonales inf :
446 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
447 if (inf_new > inf) inf = inf_new ;
448
449 ligne ++ ; /* ligne=3 */
450
451 // Utilisation des derivees ssi Plm existe dans la
452 //zone suivante :
453 if (indic[zone+1] == 1) {
454
455 //derivee de (x+echelle)^l en 1 :
456 systeme->set(ligne, colonne) =
457 l_quant*pow(echelle+1, double(l_quant-1))/alpha ; /* ligne=3, colonne=1 */
458
459 //derivee de 1/(echelle+x)^(l+1) en 1 :
460 systeme->set(ligne, colonne+1) =
461 -(l_quant+1)/pow(echelle+1, double(l_quant+2))/alpha ; /* ligne=3, colonne=1, colonne+1=2 */
462
463 // La solution particuliere :
464 for (int i=0 ; i<nr ; i++)
465 sec_membre->set(ligne) -=
466 i*i/alpha*solution_part(zone, k, j, i) ; /* ligne=3 */
467
468 // les diagonales inf :
469 inf_new = ligne-colonne ; /* ligne=3, colonne=1, ligne-colonne=2 */
470 if (inf_new > inf) inf = inf_new ;
471
472 }
473 colonne += 2 ; /* ligne=colonne=3 -> ligne=colonne=1 next block of two */
474 } else {
475
476
477
478 //-------------------------
479 // Falloff outer boundary
480 //---------------------------
481
482 /* ligne=2, colonne=1 */
483
484 // The coefficient for A_n is (k+l)(echelle+1)^l
485 systeme->set(ligne, colonne) =
486 double(k_falloff+l_quant)*pow(echelle+1, double(l_quant)) ; /* ligne=2, colonne=1 */
487
488 // The coefficient for B_n is (k-(l+1))(echelle+1)^{-(l+1)}
489 systeme->set(ligne, colonne+1) =
490 double(k_falloff-l_quant-1)/pow(echelle+1, double(l_quant+1)) ; /* ligne=2, colonne=1, colonne+1=2 */
491
492 // Here we have -(1+echelle)*F'(x=1)-k*F(x=1)
493 // La solution particuliere :
494 for (int i=0 ; i<nr ; i++)
495 sec_membre->set(ligne) -=
496 (k_falloff+(echelle+1)*i*i)*solution_part(zone, k, j, i) ; /* ligne=2 */
497
498 // les diagonales inf :
499 inf_new = ligne-colonne ; /* ligne=2, colonne=1, ligne-colonne=1 */
500 if (inf_new > inf) inf = inf_new ;
501
502 }
503 }
504
505
506 //-------------------------
507 // resolution du systeme
508 //--------------------------
509
510 systeme->set_band(sup, inf) ;
511 systeme->set_lu() ;
512
513 Tbl facteurs(systeme->inverse(*sec_membre)) ;
514 int conte = 0 ;
515
516
517 // rangement dans le noyau :
518
519 if (indic[0] == 1) {
520 nr = source.get_mg()->get_nr(0) ;
521 for (int i=0 ; i<nr ; i++)
522 resultat.set(0, k, j, i) = solution_part(0, k, j, i)
523 +facteurs(conte)*solution_hom_un(0, k, j, i) ;
524 conte++ ;
525 }
526
527 // rangement dans les coquilles :
528 for (int zone=1 ; zone<nz ; zone++)
529 if (indic[zone] == 1) {
530 nr = source.get_mg()->get_nr(zone) ;
531 for (int i=0 ; i<nr ; i++)
532 resultat.set(zone, k, j, i) =
533 solution_part(zone, k, j, i)
534 +facteurs(conte)*solution_hom_un(zone, k, j, i)
535 +facteurs(conte+1)*solution_hom_deux(zone, k, j, i) ;
536 conte+=2 ;
537 }
538
539 delete sec_membre ;
540 delete systeme ;
541 }
542
543 }
544
545 delete [] indic ;
546
547 return resultat;
548}
549}
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