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 char 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 
99 namespace Lorene {
100 void 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 
312 void Map_af::sol_elliptic_no_zec(Param_elliptic& ope_var, const Scalar& source,
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 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1294
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:148
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2033
Lorene prototypes.
Definition: app_hor.h:64
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:138
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Base class for coordinate mappings.
Definition: map.h:670
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
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:452
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
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition: mtbl_cf.C:312
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:186
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:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
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
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:494
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:601