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 
84 namespace Lorene {
85 Mtbl_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
Lorene prototypes.
Definition: app_hor.h:67
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
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724