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 char 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 
81 namespace Lorene {
82 Mtbl_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 }
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:721