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
23char 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
112namespace Lorene {
113Mtbl_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}
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
const Mg3d * get_mg() const
Returns the Mg3d on which the Mtbl_cf is defined.
Definition mtbl_cf.h:453
#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:64