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