LORENE
map_et_poisson_regu.C
1 /*
2  * Method of the class Map_et for the (iterative) resolution of the scalar
3  * Poisson equation by using regularized source.
4  *
5  * (see file map.h for the documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2000-2001 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: map_et_poisson_regu.C,v 1.3 2016/12/05 16:17:58 j_novak Exp $
35  * $Log: map_et_poisson_regu.C,v $
36  * Revision 1.3 2016/12/05 16:17:58 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.2 2014/10/13 08:53:05 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
43  * LORENE
44  *
45  * Revision 2.8 2000/09/27 14:07:14 keisuke
46  * Traitement des bases spectrales de d_logn_auto_div.
47  *
48  * Revision 2.7 2000/09/26 15:41:20 keisuke
49  * Correction erreur: la triade de duu_div doit etre celle de *this et
50  * non celle de l'objet temporaire mpaff.
51  *
52  * Revision 2.6 2000/09/25 15:03:34 keisuke
53  * Correct the derivative duu_div.
54  *
55  * Revision 2.5 2000/09/11 14:00:20 keisuke
56  * Suppress "uu = uu_regu + uu_div" because of double setting (in poisson_regular).
57  *
58  * Revision 2.4 2000/09/07 15:51:29 keisuke
59  * Minor change.
60  *
61  * Revision 2.3 2000/09/07 15:30:07 keisuke
62  * Add a new argument Cmp& uu.
63  *
64  * Revision 2.2 2000/09/04 15:56:15 keisuke
65  * Change the argumant Cmp& duu_div_r into Tenseur& duu_div.
66  *
67  * Revision 2.1 2000/09/04 14:52:17 keisuke
68  * Change the scheme of code into that of map_et_poisson.C.
69  *
70  * Revision 2.0 2000/09/01 09:55:33 keisuke
71  * *** empty log message ***
72  *
73  *
74  * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson_regu.C,v 1.3 2016/12/05 16:17:58 j_novak Exp $
75  *
76  */
77 
78 // Header Lorene:
79 #include "map.h"
80 #include "cmp.h"
81 #include "tenseur.h"
82 #include "param.h"
83 
84 //*****************************************************************************
85 
86 namespace Lorene {
87 
88 void Map_et::poisson_regular(const Cmp& source, int k_div, int nzet,
89  double unsgam1, Param& par, Cmp& uu,
90  Cmp& uu_regu, Cmp& uu_div, Tenseur& duu_div,
91  Cmp& source_regu, Cmp& source_div) const {
92 
93 
94  assert(source.get_etat() != ETATNONDEF) ;
95  assert(source.get_mp() == this) ;
96 
97  assert( source.check_dzpuis(2) || source.check_dzpuis(4)
98  || source.check_dzpuis(3)) ;
99 
100  assert(uu.get_mp() == this) ;
101  assert(uu.check_dzpuis(0)) ;
102 
103  int nz = mg->get_nzone() ;
104  int nzm1 = nz - 1 ;
105 
106  // Indicator of existence of a compactified external domain
107  bool zec = false ;
108  if (mg->get_type_r(nzm1) == UNSURR) {
109  zec = true ;
110  }
111 
112  //-------------------------------
113  // Computation of the prefactor a ---> Cmp apre
114  //-------------------------------
115 
116  Mtbl unjj = 1 + srdrdt*srdrdt + srstdrdp*srstdrdp ;
117 
118  Mtbl apre1(*mg) ;
119  apre1.set_etat_qcq() ;
120  for (int l=0; l<nz; l++) {
121  *(apre1.t[l]) = alpha[l]*alpha[l] ;
122  }
123 
124  apre1 = apre1 * dxdr * dxdr * unjj ;
125 
126  Cmp apre(*this) ;
127  apre = apre1 ;
128 
129  Tbl amax0 = max(apre1) ; // maximum values in each domain
130 
131  // The maximum values of a in each domain are put in a Mtbl
132  Mtbl amax1(*mg) ;
133  amax1.set_etat_qcq() ;
134  for (int l=0; l<nz; l++) {
135  *(amax1.t[l]) = amax0(l) ;
136  }
137 
138  Cmp amax(*this) ;
139  amax = amax1 ;
140 
141  //-------------------
142  // Initializations
143  //-------------------
144 
145  int nitermax = par.get_int() ;
146  int& niter = par.get_int_mod() ;
147  double lambda = par.get_double() ;
148  double unmlambda = 1. - lambda ;
149  double precis = par.get_double(1) ;
150 
151  Cmp& ssj = par.get_cmp_mod() ;
152 
153  Cmp ssjm1 = ssj ;
154  Cmp ssjm2 = ssjm1 ;
155 
156  Valeur& vuu = uu.va ;
157 
158  Valeur vuujm1(*mg) ;
159  if (uu.get_etat() == ETATZERO) {
160  vuujm1 = 1 ; // to take relative differences
161  vuujm1.set_base( vuu.base ) ;
162  }
163  else{
164  vuujm1 = vuu ;
165  }
166 
167  // Affine mapping for the Laplacian-tilde
168 
169  Map_af mpaff(*this) ;
170  Param par_nul ;
171 
172  cout << "Map_et::poisson_regular : relat. diff. u^J <-> u^{J-1} : "
173  << endl ;
174 
175 //==========================================================================
176 //==========================================================================
177 // Start of iteration
178 //==========================================================================
179 //==========================================================================
180 
181  Tbl tdiff(nz) ;
182  double diff ;
183  niter = 0 ;
184 
185  do {
186 
187  //====================================================================
188  // Computation of R(u) (the result is put in uu)
189  //====================================================================
190 
191 
192  //------------------------
193  // First operations on uu
194  //------------------------
195 
196  Valeur duudx = (uu.va).dsdx() ; // d/dx
197 
198  const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
199 
200  const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
201 
202 
203  //------------------
204  // Angular Laplacian
205  //------------------
206 
207  Valeur sxlapang = uu.va ;
208 
209  sxlapang.ylm() ;
210 
211  sxlapang = sxlapang.lapang() ;
212 
213  sxlapang = sxlapang.sx() ; // Lap_ang(uu) /x in the nucleus
214  // Lap_ang(uu) in the shells
215  // Lap_ang(uu) /(x-1) in the ZEC
216 
217  //------------------------------------------------------------------
218  // Computation of
219  // [ 2 /(dRdx) ( A - 1 ) duu/dx + 1/R (B - 1) Lap_ang(uu) ] / x
220  //
221  // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
222  // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
223  //
224  // The result is put in uu (via vuu)
225  //------------------------------------------------------------------
226 
227  Valeur varduudx = duudx ;
228 
229  if (zec) {
230  varduudx.annule(nzm1) ; // term in d/dx set to zero in the ZEC
231  }
232 
233  uu.set_etat_qcq() ;
234 
235  Base_val sauve_base = varduudx.base ;
236 
237  vuu = 2. * dxdr * ( rsxdxdr * unjj - 1.) * varduudx
238  + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
239 
240  vuu.set_base(sauve_base) ;
241 
242  vuu = vuu.sx() ;
243 
244  //----------------------------------------
245  // Computation of R(u)
246  //
247  // The result is put in uu (via vuu)
248  //----------------------------------------
249 
250  sauve_base = vuu.base ;
251 
252  vuu = xsr * vuu
253  + 2. * dxdr * ( sr2drdt * d2uudtdx
254  + sr2stdrdp * std2uudpdx ) ;
255 
256  vuu += dxdr * ( lapr_tp + dxdr * (
257  dxdr* unjj * d2rdx2
258  - 2. * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
259  ) * duudx ;
260 
261  vuu.set_base(sauve_base) ;
262 
263  // Since the assignment is performed on vuu (uu.va), the treatment
264  // of uu.dzpuis must be performed by hand:
265 
266  uu.set_dzpuis(4) ;
267 
268  if (source.get_dzpuis() == 2) {
269  uu.dec2_dzpuis() ; // uu.dzpuis: 4 -> 2
270  }
271 
272  if (source.get_dzpuis() == 3) {
273  uu.dec_dzpuis() ; //uu.dzpuis 4 -> 3
274  }
275 
276  //====================================================================
277  // Computation of the effective source s^J of the ``affine''
278  // Poisson equation
279  //====================================================================
280 
281  ssj = lambda * ssjm1 + unmlambda * ssjm2 ;
282 
283  ssj = ( source + uu + (amax - apre) * ssj ) / amax ;
284 
285  (ssj.va).set_base((source.va).base) ;
286 
287  //====================================================================
288  // Resolution of the ``affine'' Poisson equation
289  //====================================================================
290 
291  if ( source.get_dzpuis() == 0 ){
292  ssj.set_dzpuis( 4 ) ;
293  }
294  else {
295  ssj.set_dzpuis( source.get_dzpuis() ) ;
296  // Choice of the resolution
297  // dzpuis = 2, 3 or 4
298  }
299 
300  assert( uu.check_dzpuis( ssj.get_dzpuis() ) ) ;
301 
302  mpaff.poisson_regular(ssj, k_div, nzet, unsgam1, par_nul, uu,
303  uu_regu, uu_div, duu_div,
304  source_regu, source_div) ;
305 
306  //======================================
307  // Gradient of the diverging part (from that computed on the Map_af)
308  //======================================
309 
310  Valeur& dr_uu_div = duu_div.set(0).va ;
311  Valeur& dt_uu_div = duu_div.set(1).va ;
312  Valeur& dp_uu_div = duu_div.set(2).va ;
313 
314  Base_val bv = dr_uu_div.base ;
315  dr_uu_div = alpha[0] * dr_uu_div * dxdr ;
316  dr_uu_div.set_base( bv ) ;
317 
318  bv = dt_uu_div.base ;
319  dt_uu_div = alpha[0] * dt_uu_div * xsr - srdrdt * dr_uu_div ;
320  dt_uu_div.set_base( bv ) ;
321 
322  bv = dp_uu_div.base ;
323  dp_uu_div = alpha[0] * dp_uu_div * xsr - srstdrdp * dr_uu_div ;
324  dp_uu_div.set_base( bv ) ;
325 
326  duu_div.set_triad( this->get_bvect_spher() ) ;
327 
328 
329  //========================================
330  // Relative difference with previous step
331  //========================================
332 
333  tdiff = diffrel(vuu, vuujm1) ;
334 
335  diff = max(tdiff) ;
336 
337  cout << " step " << niter << " : " ;
338  for (int l=0; l<nz; l++) {
339  cout << tdiff(l) << " " ;
340  }
341  cout << endl ;
342 
343  //=================================
344  // Updates for the next iteration
345  //=================================
346 
347  ssjm2 = ssjm1 ;
348  ssjm1 = ssj ;
349  vuujm1 = vuu ;
350 
351  niter++ ;
352 
353  } // End of iteration
354  while ( (diff > precis) && (niter < nitermax) ) ;
355 
356 //==========================================================================
357 //==========================================================================
358 // End of iteration
359 //==========================================================================
360 //==========================================================================
361 
362 
363 }
364 }
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
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:115
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1629
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:75
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:690
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:157
Multi-domain array.
Definition: mtbl.h:118
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2789
Lorene prototypes.
Definition: app_hor.h:67
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1621
virtual void poisson_regular(const Cmp &source, int k_div, int nzet, double unsgam1, Param &par, Cmp &uu, Cmp &uu_regu, Cmp &uu_div, Tenseur &duu_div, Cmp &source_regu, Cmp &source_div) const
Computes the solution of a scalar Poisson equation.
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 annule(int l)
Sets the Valeur to zero in a given domain.
Definition: valeur.C:747
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
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1613
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1052
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
virtual void poisson_regular(const Cmp &source, int k_div, int nzet, double unsgam1, Param &par, Cmp &uu, Cmp &uu_regu, Cmp &uu_div, Tenseur &duu_div, Cmp &source_regu, Cmp &source_div) const
Computes the solution of a scalar Poisson equation.
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
const Valeur & stdsdp() const
Returns of *this.
Definition: valeur_stdsdp.C:63
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
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1669
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:183
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1570
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
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1652
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
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
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1661