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
118namespace Lorene {
119Mtbl_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}
Bases of the spectral expansions.
Definition base_val.h:325
Affine radial mapping.
Definition map.h:2042
Matrix handling.
Definition matrice.h:152
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:463
Parameter storage.
Definition param.h:125
Basic array class.
Definition tbl.h:161
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define O2DEGE_SMALL
Operateur du deuxieme ordre degenere .
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define ORDRE1_LARGE
Operateur du premier ordre .
#define O2NOND_LARGE
Operateur du deuxieme ordre non degenere .
#define O2DEGE_LARGE
Operateur du deuxieme ordre degenere .
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:67