LORENE
map_af_elliptic.C
1/*
2 * Copyright (c) 1999-2001 Eric Gourgoulhon
3 * Copyright (c) 1999-2001 Philippe Grandclement
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char map_af_elliptic_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $" ;
25
26/*
27 * $Id: map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
28 * $Log: map_af_elliptic.C,v $
29 * Revision 1.13 2014/10/13 08:53:02 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.12 2014/10/06 15:13:12 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.11 2007/05/06 10:48:11 p_grandclement
36 * Modification of a few operators for the vorton project
37 *
38 * Revision 1.10 2007/01/16 15:08:07 n_vasset
39 * New constructor, usn Scalar on mono-domain angular grid for boundary,
40 * for function sol_elliptic_boundary()
41 *
42 * Revision 1.9 2005/11/30 11:09:07 p_grandclement
43 * Changes for the Bin_ns_bh project
44 *
45 * Revision 1.8 2005/08/26 14:02:40 p_grandclement
46 * Modification of the elliptic solver that matches with an oscillatory exterior solution
47 * small correction in Poisson tau also...
48 *
49 * Revision 1.7 2005/06/09 07:57:31 f_limousin
50 * Implement a new function sol_elliptic_boundary() and
51 * Vector::poisson_boundary(...) which solve the vectorial poisson
52 * equation (method 6) with an inner boundary condition.
53 *
54 * Revision 1.6 2004/08/24 09:14:42 p_grandclement
55 * Addition of some new operators, like Poisson in 2d... It now requieres the
56 * GSL library to work.
57 *
58 * Also, the way a variable change is stored by a Param_elliptic is changed and
59 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
60 * will requiere some modification. (It should concern only the ones about monopoles)
61 *
62 * Revision 1.5 2004/06/22 08:49:58 p_grandclement
63 * Addition of everything needed for using the logarithmic mapping
64 *
65 * Revision 1.4 2004/03/17 15:58:48 p_grandclement
66 * Slight modification of sol_elliptic_no_zec
67 *
68 * Revision 1.3 2004/02/11 09:47:46 p_grandclement
69 * Addition of a new elliptic solver, matching with the homogeneous solution
70 * at the outer shell and not solving in the external domain (more details
71 * coming soon ; check your local Lorene dealer...)
72 *
73 * Revision 1.2 2004/01/28 16:46:23 p_grandclement
74 * Addition of the sol_elliptic_fixe_der_zero stuff
75 *
76 * Revision 1.1 2003/12/11 14:48:48 p_grandclement
77 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
78 *
79 *
80 * $Header: /cvsroot/Lorene/C++/Source/Map/map_af_elliptic.C,v 1.13 2014/10/13 08:53:02 j_novak Exp $
81 *
82 */
83
84// Header C :
85#include <cstdlib>
86#include <cmath>
87
88// Headers Lorene :
89#include "tbl.h"
90#include "mtbl_cf.h"
91#include "map.h"
92#include "param_elliptic.h"
93#include "proto.h"
94
95 //----------------------------------------------
96 // General elliptic solver
97 //----------------------------------------------
98
99namespace Lorene {
100void Map_af::sol_elliptic(Param_elliptic& ope_var, const Scalar& source,
101 Scalar& pot) const {
102
103 assert(source.get_etat() != ETATNONDEF) ;
104 assert(source.get_mp().get_mg() == mg) ;
105 assert(pot.get_mp().get_mg() == mg) ;
106
107 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
108 source.check_dzpuis(4)) ;
109 // Spherical harmonic expansion of the source
110 // ------------------------------------------
111
112 const Valeur& sourva = source.get_spectral_va() ;
113
114 if (sourva.get_etat() == ETATZERO) {
115 pot.set_etat_zero() ;
116 return ;
117 }
118
119 // Spectral coefficients of the source
120 assert(sourva.get_etat() == ETATQCQ) ;
121
122 Valeur rho(sourva.get_mg()) ;
123 sourva.coef() ;
124 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
125
126 rho.ylm() ; // spherical harmonic transforms
127
128 // On met les bonnes bases dans le changement de variable
129 // de ope_var :
130 ope_var.var_F.set_spectral_va().coef() ;
131 ope_var.var_F.set_spectral_va().ylm() ;
132 ope_var.var_G.set_spectral_va().coef() ;
133 ope_var.var_G.set_spectral_va().ylm() ;
134
135 // Call to the Mtbl_cf version
136 // ---------------------------
137 Mtbl_cf resu = elliptic_solver (ope_var, *(rho.c_cf)) ;
138
139 // Final result returned as a Scalar
140 // ---------------------------------
141
142 pot.set_etat_zero() ; // to call Scalar::del_t().
143
144 pot.set_etat_qcq() ;
145
146 pot.set_spectral_va() = resu ;
147 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
148
149 pot.set_dzpuis(0) ;
150}
151
152
153
154 //-----------------------------------------------
155 // General elliptic solver with boundary as Mtbl-cf
156 //-------------------------------------------------
157
158
160 Scalar& pot, const Mtbl_cf& bound,
161 double fact_dir, double fact_neu ) const {
162
163 assert(source.get_etat() != ETATNONDEF) ;
164 assert(source.get_mp().get_mg() == mg) ;
165 assert(pot.get_mp().get_mg() == mg) ;
166
167 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
168 source.check_dzpuis(4)) ;
169 // Spherical harmonic expansion of the source
170 // ------------------------------------------
171
172 const Valeur& sourva = source.get_spectral_va() ;
173
174 if (sourva.get_etat() == ETATZERO) {
175 pot.set_etat_zero() ;
176 return ;
177 }
178
179 // Spectral coefficients of the source
180 assert(sourva.get_etat() == ETATQCQ) ;
181
182 Valeur rho(sourva.get_mg()) ;
183 sourva.coef() ;
184 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
185
186 rho.ylm() ; // spherical harmonic transforms
187
188 // On met les bonnes bases dans le changement de variable
189 // de ope_var :
190 ope_var.var_F.set_spectral_va().coef() ;
191 ope_var.var_F.set_spectral_va().ylm() ;
192 ope_var.var_G.set_spectral_va().coef() ;
193 ope_var.var_G.set_spectral_va().ylm() ;
194
195 // Call to the Mtbl_cf version
196 // ---------------------------
197 Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound,
198 fact_dir, fact_neu) ;
199
200 // Final result returned as a Scalar
201 // ---------------------------------
202
203 pot.set_etat_zero() ; // to call Scalar::del_t().
204
205 pot.set_etat_qcq() ;
206
207 pot.set_spectral_va() = resu ;
208 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
209
210 pot.set_dzpuis(0) ;
211}
212
213
214
215
216
217 //-----------------------------------------------
218 // General elliptic solver with boundary as Scalar
219 //-------------------------------------------------
220
221
223 Scalar& pot, const Scalar& bound,
224 double fact_dir, double fact_neu ) const {
225
226 assert(source.get_etat() != ETATNONDEF) ;
227 assert(source.get_mp().get_mg() == mg) ;
228 assert(pot.get_mp().get_mg() == mg) ;
229
230 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
231 source.check_dzpuis(4)) ;
232 // Spherical harmonic expansion of the source
233 // ------------------------------------------
234
235 const Valeur& sourva = source.get_spectral_va() ;
236
237 if (sourva.get_etat() == ETATZERO) {
238 pot.set_etat_zero() ;
239 return ;
240 }
241
242 // Spectral coefficients of the source
243 assert(sourva.get_etat() == ETATQCQ) ;
244
245 Valeur rho(sourva.get_mg()) ;
246 sourva.coef() ;
247 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
248
249 rho.ylm() ; // spherical harmonic transforms
250
251 // On met les bonnes bases dans le changement de variable
252 // de ope_var :
253 ope_var.var_F.set_spectral_va().coef() ;
254 ope_var.var_F.set_spectral_va().ylm() ;
255 ope_var.var_G.set_spectral_va().coef() ;
256 ope_var.var_G.set_spectral_va().ylm() ;
257
258 // Call to the Mtbl_cf version
259 // ---------------------------
260
261// REMINDER : The scalar bound must be defined on a mono-domain angular grid corresponding with the full grid of the scalar source (following assert())
262
263 Scalar bbound = bound;
264 bbound.set_spectral_va().ylm() ;
265 const Map& mapp = bbound.get_mp();
266
267 const Mg3d& gri2d = *mapp.get_mg();
268
269 assert(&gri2d == source.get_mp().get_mg()->get_angu_1dom()) ; // Attention à cet assert !!
270
271 Mtbl_cf bound2 (gri2d , bbound.get_spectral_base()) ;
272 bound2.annule_hard() ;
273
274 if (bbound.get_etat() != ETATZERO){
275
276 int nr = gri2d.get_nr(0) ;
277 int nt = gri2d.get_nt(0) ;
278 int np = gri2d.get_np(0) ;
279
280 for(int k=0; k<np+2; k++)
281 for (int j=0; j<=nt-1; j++)
282 for(int xi=0; xi<= nr-1; xi++)
283 {
284
285 bound2.set(0, k , j , xi) = (*bbound.get_spectral_va().c_cf)(0, k, j, xi) ;
286 }
287 }
288 Mtbl_cf resu = elliptic_solver_boundary (ope_var, *(rho.c_cf), bound2,
289 fact_dir, fact_neu) ;
290
291 // Final result returned as a Scalar
292 // ---------------------------------
293
294 pot.set_etat_zero() ; // to call Scalar::del_t().
295
296 pot.set_etat_qcq() ;
297
298 pot.set_spectral_va() = resu ;
299 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
300
301 pot.set_dzpuis(0) ;
302}
303
304
305
306
307
308 //----------------------------------------------
309 // General elliptic solver with no ZEC
310 //----------------------------------------------
311
313 Scalar& pot, double val) const {
314
315 assert(source.get_etat() != ETATNONDEF) ;
316 assert(source.get_mp().get_mg() == mg) ;
317 assert(pot.get_mp().get_mg() == mg) ;
318
319 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
320 source.check_dzpuis(4)) ;
321 // Spherical harmonic expansion of the source
322 // ------------------------------------------
323
324 const Valeur& sourva = source.get_spectral_va() ;
325
326 if (sourva.get_etat() == ETATZERO) {
327 pot.set_etat_zero() ;
328 return ;
329 }
330
331 // Spectral coefficients of the source
332 assert(sourva.get_etat() == ETATQCQ) ;
333
334 Valeur rho(sourva.get_mg()) ;
335 sourva.coef() ;
336 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
337
338 rho.ylm() ; // spherical harmonic transforms
339
340 // On met les bonnes bases dans le changement de variable
341 // de ope_var :
342 ope_var.var_F.set_spectral_va().coef() ;
343 ope_var.var_F.set_spectral_va().ylm() ;
344 ope_var.var_G.set_spectral_va().coef() ;
345 ope_var.var_G.set_spectral_va().ylm() ;
346
347 // Call to the Mtbl_cf version
348 // ---------------------------
349 Mtbl_cf resu = elliptic_solver_no_zec (ope_var, *(rho.c_cf), val) ;
350
351 // Final result returned as a Scalar
352 // ---------------------------------
353
354 pot.set_etat_zero() ; // to call Scalar::del_t().
355
356 pot.set_etat_qcq() ;
357
358 pot.set_spectral_va() = resu ;
359 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
360
361 pot.set_dzpuis(0) ;
362}
363
364 //----------------------------------------------
365 // General elliptic solver only in the ZEC
366 //----------------------------------------------
367
369 const Scalar& source,
370 Scalar& pot, double val) const {
371
372 assert(source.get_etat() != ETATNONDEF) ;
373 assert(source.get_mp().get_mg() == mg) ;
374 assert(pot.get_mp().get_mg() == mg) ;
375
376 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
377 source.check_dzpuis(4)) ;
378 // Spherical harmonic expansion of the source
379 // ------------------------------------------
380
381 const Valeur& sourva = source.get_spectral_va() ;
382
383 if (sourva.get_etat() == ETATZERO) {
384 pot.set_etat_zero() ;
385 return ;
386 }
387
388 // Spectral coefficients of the source
389 assert(sourva.get_etat() == ETATQCQ) ;
390
391 Valeur rho(sourva.get_mg()) ;
392 sourva.coef() ;
393 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
394
395 rho.ylm() ; // spherical harmonic transforms
396
397 // On met les bonnes bases dans le changement de variable
398 // de ope_var :
399 ope_var.var_F.set_spectral_va().coef() ;
400 ope_var.var_F.set_spectral_va().ylm() ;
401 ope_var.var_G.set_spectral_va().coef() ;
402 ope_var.var_G.set_spectral_va().ylm() ;
403
404 // Call to the Mtbl_cf version
405 // ---------------------------
406 Mtbl_cf resu = elliptic_solver_only_zec (ope_var, *(rho.c_cf), val) ;
407
408 // Final result returned as a Scalar
409 // ---------------------------------
410
411 pot.set_etat_zero() ; // to call Scalar::del_t().
412
413 pot.set_etat_qcq() ;
414
415 pot.set_spectral_va() = resu ;
416 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
417
418 pot.set_dzpuis(0) ;
419}
420
421 //----------------------------------------------
422 // General elliptic solver with no ZEC
423 // and a mtaching with sin(r)/r
424 //----------------------------------------------
425
427 const Scalar& source, Scalar& pot, double* amplis, double* phases) const {
428
429 assert(source.get_etat() != ETATNONDEF) ;
430 assert(source.get_mp().get_mg() == mg) ;
431 assert(pot.get_mp().get_mg() == mg) ;
432
433 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
434 source.check_dzpuis(4)) ;
435 // Spherical harmonic expansion of the source
436 // ------------------------------------------
437
438 const Valeur& sourva = source.get_spectral_va() ;
439
440 if (sourva.get_etat() == ETATZERO) {
441 pot.set_etat_zero() ;
442 return ;
443 }
444
445 // Spectral coefficients of the source
446 assert(sourva.get_etat() == ETATQCQ) ;
447
448 Valeur rho(sourva.get_mg()) ;
449 sourva.coef() ;
450 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
451
452 rho.ylm() ; // spherical harmonic transforms
453
454 // On met les bonnes bases dans le changement de variable
455 // de ope_var :
456 ope_var.var_F.set_spectral_va().coef() ;
457 ope_var.var_F.set_spectral_va().ylm() ;
458 ope_var.var_G.set_spectral_va().coef() ;
459 ope_var.var_G.set_spectral_va().ylm() ;
460
461 // Call to the Mtbl_cf version
462 // ---------------------------
463 Mtbl_cf resu = elliptic_solver_sin_zec (ope_var, *(rho.c_cf), amplis, phases) ;
464
465 // Final result returned as a Scalar
466 // ---------------------------------
467
468 pot.set_etat_zero() ; // to call Scalar::del_t().
469
470 pot.set_etat_qcq() ;
471
472 pot.set_spectral_va() = resu ;
473 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
474
475 pot.set_dzpuis(0) ;
476}
477
478
479 //----------------------------------------------
480 // General elliptic solver with no ZEC
481 //----------------------------------------------
482
484 Param_elliptic& ope_var,
485 const Scalar& source,
486 Scalar& pot) const {
487
488 assert(source.get_etat() != ETATNONDEF) ;
489 assert(source.get_mp().get_mg() == mg) ;
490 assert(pot.get_mp().get_mg() == mg) ;
491
492 assert(source.check_dzpuis(2) || source.check_dzpuis(3) ||
493 source.check_dzpuis(4)) ;
494 // Spherical harmonic expansion of the source
495 // ------------------------------------------
496
497 const Valeur& sourva = source.get_spectral_va() ;
498
499 if (sourva.get_etat() == ETATZERO) {
500 pot.set_etat_zero() ;
501 return ;
502 }
503
504 // Spectral coefficients of the source
505 assert(sourva.get_etat() == ETATQCQ) ;
506
507 Valeur rho(sourva.get_mg()) ;
508 sourva.coef() ;
509 rho = *(sourva.c_cf) ; // copy of the coefficients of the source
510
511 rho.ylm() ; // spherical harmonic transforms
512
513 // On met les bonnes bases dans le changement de variable
514 // de ope_var :
515 ope_var.var_F.set_spectral_va().coef() ;
516 ope_var.var_F.set_spectral_va().ylm() ;
517 ope_var.var_G.set_spectral_va().coef() ;
518 ope_var.var_G.set_spectral_va().ylm() ;
519
520 // Call to the Mtbl_cf version
521 // ---------------------------
522 valeur *= alpha[0] ;
523 Mtbl_cf resu = elliptic_solver_fixe_der_zero (valeur, ope_var, *(rho.c_cf)) ;
524
525 // Final result returned as a Scalar
526 // ---------------------------------
527
528 pot.set_etat_zero() ; // to call Scalar::del_t().
529
530 pot.set_etat_qcq() ;
531
532 pot.set_spectral_va() = resu ;
533 pot.set_spectral_va().ylm_i() ; // On repasse en base standard.
534
535 pot.set_dzpuis(0) ;
536}
537
538}
void sol_elliptic_boundary(Param_elliptic &params, const Scalar &so, Scalar &uu, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
General elliptic solver including inner boundary conditions.
void sol_elliptic(Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver.
void sol_elliptic_sin_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double *coefs, double *) const
General elliptic solver.
void sol_elliptic_fixe_der_zero(double val, Param_elliptic &params, const Scalar &so, Scalar &uu) const
General elliptic solver fixing the derivative at the origin and relaxing the continuity of the first ...
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2033
void sol_elliptic_only_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double val) const
General elliptic solver.
void sol_elliptic_no_zec(Param_elliptic &params, const Scalar &so, Scalar &uu, double val) const
General elliptic solver.
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition mg3d.C:494
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
This class contains the parameters needed to call the general elliptic solver.
Scalar var_G
Multiplicative variable change that must be sphericaly symetric !
Scalar var_F
Additive variable change function.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition scalar.C:873
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Values and coefficients of a (real-value) function.
Definition valeur.h:287
int get_etat() const
Returns the logical state.
Definition valeur.h:726
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
const Mg3d * get_mg() const
Returns the Mg3d on which the this is defined.
Definition valeur.h:729
void ylm_i()
Inverse of ylm()
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Lorene prototypes.
Definition app_hor.h:64