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