LORENE
dalembert.C
1 /*
2  * Copyright (c) 2000-2001 Jerome Novak
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char dalembert_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/dalembert.C,v 1.13 2014/10/13 08:53:28 j_novak Exp $" ;
24 
25 /*
26  * $Id: dalembert.C,v 1.13 2014/10/13 08:53:28 j_novak Exp $
27  * $Log: dalembert.C,v $
28  * Revision 1.13 2014/10/13 08:53:28 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.12 2014/10/06 15:16:08 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.11 2013/06/05 15:10:43 j_novak
35  * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
36  * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
37  *
38  * Revision 1.10 2008/08/27 08:51:15 jl_cornou
39  * Added Jacobi(0,2) polynomials
40  *
41  * Revision 1.9 2006/08/31 08:56:40 j_novak
42  * Added the possibility to have a shift in the quantum number l in the operator.
43  *
44  * Revision 1.8 2004/10/05 15:44:21 j_novak
45  * Minor speed enhancements.
46  *
47  * Revision 1.7 2004/03/01 09:57:03 j_novak
48  * the wave equation is solved with Scalars. It now accepts a grid with a
49  * compactified external domain, which the solver ignores and where it copies
50  * the values of the field from one time-step to the next.
51  *
52  * Revision 1.6 2003/12/19 16:21:49 j_novak
53  * Shadow hunt
54  *
55  * Revision 1.5 2003/07/25 08:31:20 j_novak
56  * Error corrected in the case of only nucleus
57  *
58  * Revision 1.4 2003/06/18 08:45:27 j_novak
59  * In class Mg3d: added the member get_radial, returning only a radial grid
60  * For dAlembert solver: the way the coefficients of the operator are defined has been changed.
61  *
62  * Revision 1.3 2002/01/03 13:18:41 j_novak
63  * Optimization: the members set(i,j) and operator(i,j) of class Matrice are
64  * now defined inline. Matrice is a friend class of Tbl.
65  *
66  * Revision 1.2 2002/01/02 14:07:57 j_novak
67  * Dalembert equation is now solved in the shells. However, the number of
68  * points in theta and phi must be the same in each domain. The solver is not
69  * completely tested (beta version!).
70  *
71  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
72  * LORENE
73  *
74  * Revision 1.1 2000/12/04 14:24:15 novak
75  * Initial revision
76  *
77  *
78  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/dalembert.C,v 1.13 2014/10/13 08:53:28 j_novak Exp $
79  *
80  */
81 
82 
83 // Header C :
84 #include <cmath>
85 
86 // Headers Lorene :
87 #include "param.h"
88 #include "matrice.h"
89 #include "map.h"
90 #include "base_val.h"
91 #include "proto.h"
92 
93 
94  //----------------------------------------------
95  // Version Mtbl_cf
96  //----------------------------------------------
97 
98 /*
99  *
100  * Solution de l'equation de d'Alembert
101  *
102  * Entree : mapping : le mapping affine
103  * source : les coefficients de la source
104  * La base de decomposition doit etre Ylm
105  * Sortie : renvoie les coefficients de la solution dans la meme base de
106  * decomposition (a savoir Ylm)
107  *
108  */
109 
110 
111 
112 namespace Lorene {
113 Mtbl_cf sol_dalembert(Param& par, const Map_af& mapping, const Mtbl_cf& source)
114 {
115 
116  // Verifications d'usage sur les zones
117  int nz = source.get_mg()->get_nzone() ;
118  bool ced = (source.get_mg()->get_type_r(nz-1) == UNSURR ) ;
119  int nz0 = (ced ? nz - 1 : nz ) ;
120  assert ((source.get_mg()->get_type_r(0) == RARE)||(source.get_mg()->get_type_r(0) == FIN)) ;
121  for (int l=1 ; l<nz0 ; l++) {
122  assert(source.get_mg()->get_type_r(l) == FIN) ;
123  assert(source.get_mg()->get_nt(l) == source.get_mg()->get_nt(0)) ;
124  assert(source.get_mg()->get_np(l) == source.get_mg()->get_np(0)) ;
125  } // Same number of points in theta and phi in all domains...
126  assert (par.get_n_double() > 0) ;
127  assert (par.get_n_tbl_mod() > 1) ;
128 
129  //Is there a shift in the quantum number l?
130  int dl = 0 ; //value of the shift
131  int l_min = 0 ; //the wave equation is solved only for l+dl >= l_min
132  if (par.get_n_int() > 1) {
133  dl = -1 ;
134  l_min = par.get_int(1) ;
135  }
136 
137  // Bases spectrales
138  const Base_val& base = source.base ;
139 
140  // donnees sur la zone
141  int nr, nt, np ;
142  int base_r, type_dal ;
143  double alpha, beta ;
144  int l_quant, m_quant;
145  nt = source.get_mg()->get_nt(0) ;
146  np = source.get_mg()->get_np(0) ;
147 
148  //Rangement des valeurs intermediaires
149  Tbl *so ;
150  Tbl *sol_hom ;
151  Tbl *sol_hom2 ;
152  Tbl *sol_part ;
153 
154  // Rangement des solutions, avant raccordement
155  Mtbl_cf solution_part(source.get_mg(), base) ;
156  Mtbl_cf solution_hom_un(source.get_mg(), base) ;
157  Mtbl_cf solution_hom_deux(source.get_mg(), base) ;
158  Mtbl_cf resultat(source.get_mg(), base) ;
159 
160  solution_part.set_etat_qcq() ;
161  solution_hom_un.set_etat_qcq() ;
162  solution_hom_deux.set_etat_qcq() ;
163  resultat.annule_hard() ;
164 
165  // Tbls for the boundary condition
166  double* bc1 = &par.get_double_mod(1) ;
167  double* bc2 = &par.get_double_mod(2) ;
168  Tbl* tbc3 = &par.get_tbl_mod(1) ;
169 
170  for (int l=0 ; l<nz ; l++) {
171  solution_part.t[l]->annule_hard() ;
172  solution_hom_un.t[l]->annule_hard() ;
173  solution_hom_deux.t[l]->annule_hard() ;
174  }
175 
176  //---------------
177  //-- NUCLEUS ---
178  //---------------
179  int lz = 0 ;
180  nr = source.get_mg()->get_nr(lz) ;
181  so = new Tbl(nr) ;
182 
183  alpha = mapping.get_alpha()[lz] ;
184 
185  for (int k=0 ; k<np+1 ; k++) {
186  for (int j=0 ; j<nt ; j++) {
187  // quantic numbers and spectral bases
188  base.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
189  assert( (source.get_mg()->get_type_r(0) == RARE) ||
190  (base_r == R_JACO02) ) ;
191  l_quant += dl ;
192  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_quant >=l_min) )
193  {
194  //Calculation of the coefficients of the operator
195  par.get_tbl_mod().set(4,lz) = 2*par.get_tbl_mod()(2,lz) ;
196  par.get_tbl_mod().set(5,lz) = 2*par.get_tbl_mod()(3,lz) ;
197  par.get_tbl_mod().set(6,lz) = 2*par.get_tbl_mod()(1,lz) ;
198  par.get_tbl_mod().set(7,lz) =
199  -l_quant*(l_quant+1)*par.get_tbl_mod()(3,lz) ;
200  par.get_tbl_mod().set(8,lz) =
201  -l_quant*(l_quant+1)*par.get_tbl_mod()(2,lz) ;
202  par.get_tbl_mod().set(9,lz) =
203  -l_quant*(l_quant+1)*par.get_tbl_mod()(1,lz) ;
204 
205  Matrice operateur(nr,nr) ;
206 
207  get_operateur_dal(par, lz, base_r, type_dal, operateur) ;
208 
209  // Getting the particular solution
210  so->set_etat_qcq() ;
211  for (int i=0 ; i<nr ; i++)
212  so->set(i) = source(lz, k, j, i) ;
213  if ((type_dal == ORDRE1_LARGE) || (type_dal == O2DEGE_LARGE)
214  || (type_dal == O2NOND_LARGE))
215  so->set(nr-1) = 0 ;
216  sol_part = new Tbl(dal_inverse(base_r, type_dal, operateur,
217  *so, true)) ;
218 
219  // Getting the homogeneous solution
220  sol_hom = new Tbl(dal_inverse(base_r, type_dal, operateur,
221  *so, false)) ;
222 
223  // Putting to Mtbl_cf
224  for (int i=0 ; i<nr ; i++) {
225  solution_part.set(lz, k, j, i) = (*sol_part)(i) ;
226  solution_hom_un.set(lz, k, j, i) = (*sol_hom)(i) ;
227  solution_hom_deux.set(lz, k, j, i) = 0. ;
228  }
229 
230  // If only one zone, the BC is set
231  if (nz0 == 1) {
232 
233  int base_pipo = 0 ;
234  double part, dpart, hom, dhom;
235  Tbl der_part(3,1,nr) ;
236  der_part.set_etat_qcq() ;
237  for (int i=0; i<nr; i++)
238  der_part.set(0,0,i) = (*sol_part)(i) ;
239  Tbl der_hom(3,1,nr) ;
240  der_hom.set_etat_qcq() ;
241  for (int i=0; i<nr; i++)
242  der_hom.set(0,0,i) = (*sol_hom)(i) ;
243 
244  if (base_r == R_CHEBP) {
245  som_r_chebp(sol_part->t, nr, 1, 1, 1., &part) ;
246  _dsdx_r_chebp(&der_part, base_pipo) ;
247  som_r_chebi(der_part.t, nr, 1, 1, 1., &dpart) ;
248  som_r_chebp(sol_hom->t, nr, 1, 1, 1., &hom) ;
249  _dsdx_r_chebp(&der_hom, base_pipo) ;
250  som_r_chebi(der_hom.t, nr, 1, 1, 1., &dhom) ;
251  }
252  else {
253  som_r_chebi(sol_part->t, nr, 1, 1, 1., &part) ;
254  _dsdx_r_chebi(&der_part, base_pipo) ;
255  som_r_chebp(der_part.t, nr, 1, 1, 1., &dpart) ;
256  som_r_chebi(sol_hom->t, nr, 1, 1, 1., &hom) ;
257  _dsdx_r_chebi(&der_hom, base_pipo) ;
258  som_r_chebp(der_hom.t, nr, 1, 1, 1., &dhom) ;
259  }
260 
261  part = part*(*bc1) + dpart*(*bc2)/alpha ;
262  hom = hom*(*bc1) + dhom*(*bc2)/alpha ;
263  double lambda = ((*tbc3)(k,j) - part) / hom ;
264  for (int i=0 ; i<nr ; i++)
265  resultat.set(lz, k, j, i) =
266  solution_part(lz, k, j, i)
267  +lambda*solution_hom_un(lz, k, j, i) ;
268  }
269 
270  delete sol_hom ;
271  delete sol_part ;
272  } // nullite_plm
273  } // theta loop
274  } // phi loop
275  delete so ;
276 
277  //---------------------
278  //-- SHELLS --
279  //---------------------
280  for (lz=1 ; lz<nz0 ; lz++) {
281  nr = source.get_mg()->get_nr(lz) ;
282  so = new Tbl(nr) ;
283  alpha = mapping.get_alpha()[lz] ;
284  beta = mapping.get_beta()[lz] ;
285 
286  for (int k=0 ; k<np+1 ; k++)
287  for (int j=0 ; j<nt ; j++) {
288  // quantic numbers and spectral bases
289  base.give_quant_numbers(lz, k, j, m_quant, l_quant, base_r) ;
290  l_quant += dl ;
291  if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_quant >=l_min) )
292  {
293  //Calculation of the coefficients of the operator
294  par.get_tbl_mod().set(4,lz) = 2*par.get_tbl_mod()(2,lz) ;
295  par.get_tbl_mod().set(5,lz) = 2*par.get_tbl_mod()(3,lz) ;
296  par.get_tbl_mod().set(6,lz) = 2*par.get_tbl_mod()(1,lz) ;
297  par.get_tbl_mod().set(7,lz) =
298  -l_quant*(l_quant+1)*par.get_tbl_mod()(3,lz) ;
299  par.get_tbl_mod().set(8,lz) =
300  -l_quant*(l_quant+1)*par.get_tbl_mod()(2,lz) ;
301  par.get_tbl_mod().set(9,lz) =
302  -l_quant*(l_quant+1)*par.get_tbl_mod()(1,lz) ;
303 
304  Matrice operateur(nr,nr) ;
305 
306  get_operateur_dal(par, lz, base_r, type_dal, operateur) ;
307 
308  // Calcul DES DEUX SH
309  so->set_etat_qcq() ;
310  for (int i=0; i<nr; i++) so->set(i) = 0. ;
311  so->set(nr-2) = 1. ;
312  sol_hom = new Tbl(dal_inverse(base_r, type_dal, operateur, *so,
313  false)) ;
314  so->set(nr-2) = 0. ;
315  so->set(nr-1) = 1. ;
316  sol_hom2 = new Tbl(dal_inverse(base_r, type_dal, operateur, *so,
317  false)) ;
318 
319  // Calcul de la SP
320  double *tmp = new double[nr] ;
321  for (int i=0 ; i<nr ; i++)
322  tmp[i] = source(lz, k, j, i) ;
323  if ((type_dal == O2DEGE_SMALL) || (type_dal == O2DEGE_LARGE)) {
324  for (int i=0; i<nr; i++) so->set(i) = beta*tmp[i] ;
325  multx_1d(nr, &tmp, R_CHEB) ;
326  for (int i=0; i<nr; i++) so->set(i) += alpha*tmp[i] ;
327  }
328  else {
329  for (int i=0; i<nr; i++) so->set(i) = beta*beta*tmp[i] ;
330  multx_1d(nr, &tmp, R_CHEB) ;
331  for (int i=0; i<nr; i++) so->set(i) += 2*alpha*beta*tmp[i] ;
332  multx_1d(nr, &tmp, R_CHEB) ;
333  for (int i=0; i<nr; i++) so->set(i) += alpha*alpha*tmp[i] ;
334  }
335  so->set(nr-2) = 0. ;
336  so->set(nr-1) = 0. ;
337 
338  sol_part = new Tbl (dal_inverse(base_r, type_dal, operateur,
339  *so, true)) ;
340  // Rangement
341  for (int i=0 ; i<nr ; i++) {
342  solution_part.set(lz, k, j, i) = (*sol_part)(i) ;
343  solution_hom_un.set(lz, k, j, i) = (*sol_hom)(i) ;
344  solution_hom_deux.set(lz, k, j, i) = (*sol_hom2)(i) ;
345  }
346 
347  delete [] tmp ;
348  delete sol_hom ;
349  delete sol_hom2 ;
350  delete sol_part ;
351  }
352  } // theta loop
353  delete so ;
354  } // domain loop
355  if (nz0 > 1) {
356  //--------------------------------------------------------------------
357  //
358  // Combinations of particular and homogeneous solutions
359  // to verify continuity and boundary conditions
360  //
361  //--------------------------------------------------------------------
362  int taille = 2*nz0 - 1 ;
363  Tbl deuz(taille) ;
364  deuz.set_etat_qcq() ;
365  Matrice systeme(taille,taille) ;
366  systeme.set_etat_qcq() ;
367  int sup = 2;
368  int inf = (nz0>2) ? 2 : 1 ;
369  for (int k=0; k<np+1; k++) {
370  for (int j=0; j<nt; j++) {
371  // To get the r basis in the nucleus
372  base.give_quant_numbers(0, k, j, m_quant, l_quant, base_r) ;
373  if ( (nullite_plm(j, nt, k, np, base)) && (l_quant + dl >= l_min) ) {
374  assert ((base_r == R_CHEBP)||(base_r == R_CHEBI)) ;
375  int parite = (base_r == R_CHEBP) ? 0 : 1 ;
376  int l, c ;
377  double xx = 0.;
378  for (l=0; l<taille; l++)
379  for (c=0; c<taille; c++) systeme.set(l,c) = xx ;
380  for (l=0; l<taille; l++) deuz.set(l) = xx ;
381 
382  //---------
383  // Nucleus
384  //---------
385  nr = source.get_mg()->get_nr(0) ;
386  alpha = mapping.get_alpha()[0] ;
387  l=0 ; c=0 ;
388  for (int i=0; i<nr; i++)
389  systeme.set(l,c) += solution_hom_un(0, k, j, i) ;
390  for (int i=0; i<nr; i++) deuz.set(l) -= solution_part(0, k, j, i) ;
391 
392  l++ ;
393  xx = 0. ;
394  for (int i=0; i<nr; i++)
395  xx +=(2*i+parite)*(2*i+parite)
396  *solution_hom_un(0, k, j, i) ;
397  systeme.set(l,c) += xx/alpha ;
398  xx = 0. ;
399  for (int i=0; i<nr; i++) xx -= (2*i+parite)*
400  (2*i+parite)*solution_part(0, k, j, i) ;
401  deuz.set(l) += xx/alpha ;
402 
403  //----------
404  // Shells
405  //----------
406  for (lz=1; lz<nz0; lz++) {
407  nr = source.get_mg()->get_nr(lz) ;
408  alpha = mapping.get_alpha()[lz] ;
409  l-- ;
410  c = l+1 ;
411  for (int i=0; i<nr; i++)
412  if (i%2 == 0)
413  systeme.set(l,c) -= solution_hom_un(lz, k, j, i) ;
414  else
415  systeme.set(l,c) += solution_hom_un(lz, k, j, i) ;
416  c++ ;
417  for (int i=0; i<nr; i++)
418  if (i%2 == 0)
419  systeme.set(l,c) -= solution_hom_deux(lz, k, j, i) ;
420  else
421  systeme.set(l,c) += solution_hom_deux(lz, k, j, i) ;
422  for (int i=0; i<nr; i++)
423  if (i%2 == 0) deuz.set(l) += solution_part(lz, k, j, i) ;
424  else deuz.set(l) -= solution_part(lz, k, j, i) ;
425 
426  l++ ; c-- ;
427  xx = 0. ;
428  for (int i=0; i<nr; i++)
429  if (i%2 == 0)
430  xx += i*i*solution_hom_un(lz, k, j, i) ;
431  else
432  xx -= i*i*solution_hom_un(lz, k, j, i) ;
433  systeme.set(l,c) += xx/alpha ;
434  c++ ;
435  xx = 0. ;
436  for (int i=0; i<nr; i++)
437  if (i%2 == 0)
438  xx += i*i*solution_hom_deux(lz, k, j, i) ;
439  else
440  xx -= i*i*solution_hom_deux(lz, k, j, i) ;
441  systeme.set(l,c) += xx/alpha ;
442  xx = 0. ;
443  for (int i=0; i<nr; i++)
444  if (i%2 == 0) xx -= i*i*solution_part(lz, k, j, i) ;
445  else xx += i*i*solution_part(lz, k, j, i) ;
446  deuz.set(l) += xx/alpha ;
447 
448  l++ ; c--;
449  if (lz == nz0-1) { // Last domain, the outer BC is set
450  for (int i=0; i<nr; i++)
451  systeme.set(l,c) +=
452  ((*bc1)+(*bc2)*i*i/alpha)*solution_hom_un(lz, k, j, i) ;
453  c++ ;
454  for (int i=0; i<nr; i++)
455  systeme.set(l,c) +=
456  ((*bc1)+(*bc2)*i*i/alpha)*solution_hom_deux(lz, k, j, i) ;
457  for (int i=0; i<nr; i++)
458  deuz.set(l) -=
459  ((*bc1)+(*bc2)*i*i/alpha)*solution_part(lz, k, j, i) ;
460  deuz.set(l) += (*tbc3)(k,j) ;
461  }
462  else { // At least one more shell
463  for (int i=0; i<nr; i++)
464  systeme.set(l,c) += solution_hom_un(lz, k, j, i) ;
465  c++ ;
466  for (int i=0; i<nr; i++)
467  systeme.set(l,c) += solution_hom_deux(lz, k, j, i) ;
468  for (int i=0; i<nr; i++)
469  deuz.set(l) -= solution_part(lz, k, j, i) ;
470  l++ ; c-- ;
471  xx = 0. ;
472  for (int i=0; i<nr; i++) xx += i*i*solution_hom_un(lz, k, j, i) ;
473  systeme.set(l,c) += xx/alpha ;
474  c++ ;
475  xx = 0. ;
476  for (int i=0; i<nr; i++)
477  xx += i*i*solution_hom_deux(lz, k, j, i) ;
478  systeme.set(l,c) += xx/alpha ;
479  xx = 0. ;
480  for (int i=0; i<nr; i++)
481  xx -= i*i*solution_part(lz, k, j, i) ;
482  deuz.set(l) += xx/alpha ;
483  }
484  }
485 
486  //--------------------------------------
487  // Solution of the linear system
488  //--------------------------------------
489 
490  systeme.set_band(sup, inf) ;
491  systeme.set_lu() ;
492  Tbl facteur(systeme.inverse(deuz)) ;
493 
494  //Linear Combination in the nucleus
495  nr = source.get_mg()->get_nr(0) ;
496  for (int i=0; i<nr; i++)
497  resultat.set(0, k, j, i) = solution_part(0, k, j, i)
498  + facteur(0)*solution_hom_un(0, k, j, i) ;
499 
500  //Linear combination in the shells
501  for (lz=1; lz<nz0; lz++) {
502  nr = source.get_mg()->get_nr(lz) ;
503  for (int i=0; i<nr; i++)
504  resultat.set(lz, k, j, i) = solution_part(lz, k, j, i)
505  + facteur(2*lz-1)*solution_hom_un(lz, k, j, i)
506  + facteur(2*lz)*solution_hom_deux(lz, k, j, i) ;
507  }
508  }
509  } //End of j/theta loop
510  } //End of k/phi loop
511  } //End of case nz0>1
512 
513  return resultat ;
514 
515 }
516 }
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:282
Lorene prototypes.
Definition: app_hor.h:64
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
Definition: type_parite.h:286
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition: mtbl_cf.h:453
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#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
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
Definition: type_parite.h:280
#define ORDRE1_LARGE
Operateur du premier ordre .
Definition: type_parite.h:278
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166