LORENE
map_et_poisson2d.C
1 /*
2  * Method of the class Map_et for the resolution of the 2-D Poisson
3  * equation.
4  *
5  * (see file map.h for documentation).
6  */
7 
8 /*
9  * Copyright (c) 2000-2001 Eric Gourgoulhon
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: map_et_poisson2d.C,v 1.5 2016/12/05 16:17:58 j_novak Exp $
34  * $Log: map_et_poisson2d.C,v $
35  * Revision 1.5 2016/12/05 16:17:58 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.4 2014/10/13 08:53:05 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.3 2014/10/06 15:13:13 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.2 2002/02/07 14:55:58 e_gourgoulhon
45  * Corrected a bug when the source is known only in the coefficient
46  * space.
47  *
48  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
49  * LORENE
50  *
51  * Revision 2.4 2000/11/07 14:21:03 eric
52  * Correction d'une erreur dans le cas T_SIN_I (calcul de R(u)).
53  *
54  * Revision 2.3 2000/10/26 15:58:00 eric
55  * Correction cas T_COS_P : l'import de saff_q se fait par copie du Tbl.
56  *
57  * Revision 2.2 2000/10/12 15:37:43 eric
58  * Traitement des bases spectrales dans le cas T_COS_P.
59  *
60  * Revision 2.1 2000/10/11 15:15:43 eric
61  * 1ere version operationnelle.
62  *
63  * Revision 2.0 2000/10/09 13:47:17 eric
64  * *** empty log message ***
65  *
66  *
67  * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson2d.C,v 1.5 2016/12/05 16:17:58 j_novak Exp $
68  *
69  */
70 
71 // Headers C
72 #include <cmath>
73 
74 // Headers Lorene:
75 #include "map.h"
76 #include "cmp.h"
77 #include "param.h"
78 
79 //*****************************************************************************
80 
81 namespace Lorene {
82 
83 void Map_et::poisson2d(const Cmp& source_mat, const Cmp& source_quad,
84  Param& par, Cmp& uu) const {
85 
86  assert(source_mat.get_etat() != ETATNONDEF) ;
87  assert(source_quad.get_etat() != ETATNONDEF) ;
88  assert(source_mat.get_mp()->get_mg() == mg) ;
89  assert(source_quad.get_mp()->get_mg() == mg) ;
90  assert(uu.get_mp()->get_mg() == mg) ;
91 
92  assert( source_quad.check_dzpuis(4) ) ;
93 
94  double& lambda = par.get_double_mod(0) ;
95  int nz = mg->get_nzone() ;
96  int nzm1 = nz-1 ;
97 
98  // Special case of a vanishing source
99  // ----------------------------------
100 
101  if ( (source_mat.get_etat() == ETATZERO)
102  && (source_quad.get_etat() == ETATZERO) ) {
103 
104  uu = 0 ;
105  lambda = 1 ;
106  return ;
107  }
108 
109  int base_t = ((source_mat.va).base).get_base_t(0) ;
110 
111  switch (base_t) {
112 
113  //==================================================================
114  // case T_COS_P
115  //==================================================================
116 
117  case T_COS_P : {
118 
119  // Construction of a Map_af which coincides with *this on the equator
120  // ------------------------------------------------------------------
121 
122  double theta0 = M_PI / 2 ; // Equator
123  double phi0 = 0 ;
124 
125  Map_af mpaff(*this) ;
126 
127  for (int l=0 ; l<nz ; l++) {
128  double rmax = val_r(l, double(1), theta0, phi0) ;
129  switch ( mg->get_type_r(l) ) {
130  case RARE: {
131  double rmin = val_r(l, double(0), theta0, phi0) ;
132  mpaff.set_alpha(rmax - rmin, l) ;
133  mpaff.set_beta(rmin, l) ;
134  break ;
135  }
136 
137  case FIN: {
138  double rmin = val_r(l, double(-1), theta0, phi0) ;
139  mpaff.set_alpha( double(.5) * (rmax - rmin), l ) ;
140  mpaff.set_beta( double(.5) * (rmax + rmin), l) ;
141  break ;
142  }
143 
144  case UNSURR: {
145  double rmin = val_r(l, double(-1), theta0, phi0) ;
146  double umax = double(1) / rmin ;
147  double umin = double(1) / rmax ;
148  mpaff.set_alpha( double(.5) * (umin - umax), l) ;
149  mpaff.set_beta( double(.5) * (umin + umax), l) ;
150  break ;
151  }
152 
153  default: {
154  cout << "Map_et::poisson2d: unknown type_r ! " << endl ;
155  abort () ;
156  break ;
157  }
158 
159  }
160  }
161 
162  // Importation of source_mat and source_quad of the affine mapping
163  // ---------------------------------------------------------------
164  Cmp saff_m(mpaff) ;
165  saff_m.import( nzm1, source_mat ) ;
166  (saff_m.va).set_base( (source_mat.va).base ) ;
167 
168  Cmp saff_q(mpaff) ;
169 
170  // In order to use Cmp::import with dzpuis != 0 :
171  Cmp set_q = source_quad ;
172  set_q.set_dzpuis(0) ; // dzpuis artificially set to 0
173 
174  saff_q.import( nzm1, set_q ) ;
175  (saff_q.va).set_base( (set_q.va).base ) ;
176 
177  // Copy in the external domain :
178  if ( (set_q.va).get_etat() == ETATQCQ) {
179  (set_q.va).coef_i() ; // the values in configuration space are required
180  assert( (set_q.va).c->get_etat() == ETATQCQ ) ;
181  assert( (saff_q.va).c->get_etat() == ETATQCQ ) ;
182  *( (saff_q.va).c->t[nzm1] ) = *( (set_q.va).c->t[nzm1] ) ;
183  }
184 
185  // the true dzpuis is restored :
186  saff_q.set_dzpuis( source_quad.get_dzpuis() ) ;
187 
188  // Resolution of the 2-D Poisson equation on the spherical domains
189  // ---------------------------------------------------------------
190 
191  Cmp uaff(mpaff) ;
192 
193  mpaff.poisson2d(saff_m, saff_q, par, uaff) ;
194 
195  // Importation of the solution on the Map_et mapping *this
196  // -------------------------------------------------------
197 
198  uu.import(uaff) ;
199 
200  uu.va.set_base( uaff.va.base ) ; // same spectral bases
201 
202  break ;
203  }
204 
205  //==================================================================
206  // case T_SIN_I
207  //==================================================================
208 
209  case T_SIN_I : {
210 
211  //-------------------------------
212  // Computation of the prefactor a ---> Cmp apre
213  //-------------------------------
214 
215  Mtbl unjj = 1 + srdrdt*srdrdt ;
216 
217  Mtbl apre1(*mg) ;
218  apre1.set_etat_qcq() ;
219  for (int l=0; l<nz; l++) {
220  *(apre1.t[l]) = alpha[l]*alpha[l] ;
221  }
222 
223  apre1 = apre1 * dxdr * dxdr * unjj ;
224 
225  Cmp apre(*this) ;
226  apre = apre1 ;
227 
228  Tbl amax0 = max(apre1) ; // maximum values in each domain
229 
230  // The maximum values of a in each domain are put in a Mtbl
231  Mtbl amax1(*mg) ;
232  amax1.set_etat_qcq() ;
233  for (int l=0; l<nz; l++) {
234  *(amax1.t[l]) = amax0(l) ;
235  }
236 
237  Cmp amax(*this) ;
238  amax = amax1 ;
239 
240 
241  //-------------------
242  // Initializations
243  //-------------------
244 
245  int nitermax = par.get_int() ;
246  int& niter = par.get_int_mod() ;
247  double lambda_relax = par.get_double() ;
248  double unmlambda_relax = 1. - lambda_relax ;
249  double precis = par.get_double(1) ;
250 
251  Cmp& ssj = par.get_cmp_mod() ;
252 
253  Cmp ssjm1 = ssj ;
254  Cmp ssjm2 = ssjm1 ;
255 
256  Cmp ssj_q(*this) ;
257  ssj_q = 0 ;
258 
259  Valeur& vuu = uu.va ;
260 
261  Valeur vuujm1(*mg) ;
262  if (uu.get_etat() == ETATZERO) {
263  vuujm1 = 1 ; // to take relative differences
264  vuujm1.set_base( vuu.base ) ;
265  }
266  else{
267  vuujm1 = vuu ;
268  }
269 
270  // Affine mapping for the Laplacian-tilde
271 
272  Map_af mpaff(*this) ;
273 
274  cout << "Map_et::poisson2d : relat. diff. u^J <-> u^{J-1} : " << endl ;
275 
276 //==========================================================================
277 //==========================================================================
278 // Start of iteration
279 //==========================================================================
280 //==========================================================================
281 
282  Tbl tdiff(nz) ;
283  double diff ;
284  niter = 0 ;
285 
286  do {
287 
288  //====================================================================
289  // Computation of R(u) (the result is put in uu)
290  //====================================================================
291 
292 
293  //-----------------------
294  // First operations on uu
295  //-----------------------
296 
297  Valeur duudx = (uu.va).dsdx() ; // d/dx
298 
299  const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
300 
301  //-------------------
302  // 1/x d^2uu/dtheta^2
303  //-------------------
304 
305  Valeur sxlapang = uu.va ;
306 
307  sxlapang = sxlapang.d2sdt2() ;
308 
309  sxlapang = sxlapang.sx() ; // d^2(uu)/dth^2 /x in the nucleus
310  // d^2(uu)/dth^2 in the shells
311  // d^2(uu)/dth^2 /(x-1) in the ZEC
312 
313  //---------------------------------------------------------------
314  // Computation of
315  // [ (dR/dx)^{-1} ( A - 1 ) duu/dx + 1/R (B - 1) d^2uu/dth^2 ] / x
316  //
317  // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
318  // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
319  //
320  // The result is put in uu (via vuu)
321  //---------------------------------------------------------------
322 
323  // Intermediate quantity jac which value is
324  // (dR/dx)^{-1} in the nucleus and the shells
325  // +(dU/dx)^{-1} in the ZEC
326 
327  Mtbl jac = dxdr ;
328  if (mg->get_type_r(nzm1) == UNSURR) {
329  jac.annule(nzm1, nzm1) ;
330  Mtbl jac_ext = dxdr ;
331  jac_ext.annule(0, nzm1-1) ;
332  jac_ext = - jac_ext ;
333  jac = jac + jac_ext ;
334  }
335 
336  uu.set_etat_qcq() ;
337 
338  Base_val sauve_base = duudx.base ;
339 
340  vuu = jac * ( rsxdxdr * unjj - 1.) * duudx
341  + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
342 
343  vuu.set_base(sauve_base) ;
344 
345  vuu = vuu.sx() ;
346 
347  //---------------------------------------
348  // Computation of R(u)
349  //
350  // The result is put in uu (via vuu)
351  //---------------------------------------
352 
353 
354  sauve_base = vuu.base ;
355 
356  vuu = xsr * vuu
357  + 2. * dxdr * sr2drdt * d2uudtdx ;
358 
359  vuu += dxdr * ( sr2d2rdt2 + dxdr * (
360  dxdr* unjj * d2rdx2
361  - 2. * sr2drdt * d2rdtdx )
362  ) * duudx ;
363 
364  vuu.set_base(sauve_base) ;
365 
366  // Since the assignment is performed on vuu (uu.va), the treatment
367  // of uu.dzpuis must be performed by hand:
368 
369  uu.set_dzpuis(4) ;
370 
371  //====================================================================
372  // Computation of the effective source s^J of the ``affine''
373  // Poisson equation
374  //====================================================================
375 
376  ssj = lambda_relax * ssjm1 + unmlambda_relax * ssjm2 ;
377 
378  ssj = ( source_mat + source_quad + uu + (amax - apre) * ssj ) / amax ;
379 
380  (ssj.va).set_base((source_mat.va).base) ;
381 
382  //====================================================================
383  // Resolution of the ``affine'' Poisson equation
384  //====================================================================
385 
386  assert( uu.check_dzpuis( ssj.get_dzpuis() ) ) ;
387 
388  mpaff.poisson2d(ssj, ssj_q, par, uu) ;
389 
390  tdiff = diffrel(vuu, vuujm1) ;
391 
392  diff = max(tdiff) ;
393 
394 
395  cout << " step " << niter << " : " ;
396  for (int l=0; l<nz; l++) {
397  cout << tdiff(l) << " " ;
398  }
399  cout << endl ;
400 
401  //=================================
402  // Updates for the next iteration
403  //=================================
404 
405  ssjm2 = ssjm1 ;
406  ssjm1 = ssj ;
407  vuujm1 = vuu ;
408 
409  niter++ ;
410 
411  } // End of iteration
412  while ( (diff > precis) && (niter < nitermax) ) ;
413 
414 //==========================================================================
415 //==========================================================================
416 // End of iteration
417 //==========================================================================
418 //==========================================================================
419 
420 
421  break ;
422  }
423 
424  default : {
425  cout << "Map_et::poisson2d : unkown theta basis !" << endl ;
426  cout << " basis : " << hex << base_t << endl ;
427  abort() ;
428  break ;
429  }
430  } // End of switch on base_t
431 }
432 
433 
434 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1640
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1678
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:115
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition: param.C:501
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Multi-domain array.
Definition: mtbl.h:118
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2789
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const
Computes the solution of a 2-D Poisson equation.
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Definition: map_et_radius.C:95
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1621
void annule(int l_min, int l_max)
Sets the Mtbl to zero in some domains.
Definition: mtbl.C:332
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx.C:113
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:813
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1052
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:771
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1581
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:433
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition: map_af.C:760
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1570
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const
Computes the solution of a 2-D Poisson equation.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
Coord rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
Definition: map.h:2865
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:694
Bases of the spectral expansions.
Definition: base_val.h:325
Affine radial mapping.
Definition: map.h:2048
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
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: cmp.C:718
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:657
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
Basic array class.
Definition: tbl.h:164
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
void import(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
Definition: cmp_import.C:76
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1605
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
const Valeur & d2sdt2() const
Returns of *this.
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1661