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