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 char map_et_poisson_regu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson_regu.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $" ;
32 
33 /*
34  * $Id: map_et_poisson_regu.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $
35  * $Log: map_et_poisson_regu.C,v $
36  * Revision 1.2 2014/10/13 08:53:05 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
40  * LORENE
41  *
42  * Revision 2.8 2000/09/27 14:07:14 keisuke
43  * Traitement des bases spectrales de d_logn_auto_div.
44  *
45  * Revision 2.7 2000/09/26 15:41:20 keisuke
46  * Correction erreur: la triade de duu_div doit etre celle de *this et
47  * non celle de l'objet temporaire mpaff.
48  *
49  * Revision 2.6 2000/09/25 15:03:34 keisuke
50  * Correct the derivative duu_div.
51  *
52  * Revision 2.5 2000/09/11 14:00:20 keisuke
53  * Suppress "uu = uu_regu + uu_div" because of double setting (in poisson_regular).
54  *
55  * Revision 2.4 2000/09/07 15:51:29 keisuke
56  * Minor change.
57  *
58  * Revision 2.3 2000/09/07 15:30:07 keisuke
59  * Add a new argument Cmp& uu.
60  *
61  * Revision 2.2 2000/09/04 15:56:15 keisuke
62  * Change the argumant Cmp& duu_div_r into Tenseur& duu_div.
63  *
64  * Revision 2.1 2000/09/04 14:52:17 keisuke
65  * Change the scheme of code into that of map_et_poisson.C.
66  *
67  * Revision 2.0 2000/09/01 09:55:33 keisuke
68  * *** empty log message ***
69  *
70  *
71  * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson_regu.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $
72  *
73  */
74 
75 // Header Lorene:
76 #include "map.h"
77 #include "cmp.h"
78 #include "tenseur.h"
79 #include "param.h"
80 
81 //*****************************************************************************
82 
83 namespace Lorene {
84 
85 void Map_et::poisson_regular(const Cmp& source, int k_div, int nzet,
86  double unsgam1, Param& par, Cmp& uu,
87  Cmp& uu_regu, Cmp& uu_div, Tenseur& duu_div,
88  Cmp& source_regu, Cmp& source_div) const {
89 
90 
91  assert(source.get_etat() != ETATNONDEF) ;
92  assert(source.get_mp() == this) ;
93 
94  assert( source.check_dzpuis(2) || source.check_dzpuis(4)
95  || source.check_dzpuis(3)) ;
96 
97  assert(uu.get_mp() == this) ;
98  assert(uu.check_dzpuis(0)) ;
99 
100  int nz = mg->get_nzone() ;
101  int nzm1 = nz - 1 ;
102 
103  // Indicator of existence of a compactified external domain
104  bool zec = false ;
105  if (mg->get_type_r(nzm1) == UNSURR) {
106  zec = true ;
107  }
108 
109  //-------------------------------
110  // Computation of the prefactor a ---> Cmp apre
111  //-------------------------------
112 
113  Mtbl unjj = 1 + srdrdt*srdrdt + srstdrdp*srstdrdp ;
114 
115  Mtbl apre1(*mg) ;
116  apre1.set_etat_qcq() ;
117  for (int l=0; l<nz; l++) {
118  *(apre1.t[l]) = alpha[l]*alpha[l] ;
119  }
120 
121  apre1 = apre1 * dxdr * dxdr * unjj ;
122 
123  Cmp apre(*this) ;
124  apre = apre1 ;
125 
126  Tbl amax0 = max(apre1) ; // maximum values in each domain
127 
128  // The maximum values of a in each domain are put in a Mtbl
129  Mtbl amax1(*mg) ;
130  amax1.set_etat_qcq() ;
131  for (int l=0; l<nz; l++) {
132  *(amax1.t[l]) = amax0(l) ;
133  }
134 
135  Cmp amax(*this) ;
136  amax = amax1 ;
137 
138  //-------------------
139  // Initializations
140  //-------------------
141 
142  int nitermax = par.get_int() ;
143  int& niter = par.get_int_mod() ;
144  double lambda = par.get_double() ;
145  double unmlambda = 1. - lambda ;
146  double precis = par.get_double(1) ;
147 
148  Cmp& ssj = par.get_cmp_mod() ;
149 
150  Cmp ssjm1 = ssj ;
151  Cmp ssjm2 = ssjm1 ;
152 
153  Valeur& vuu = uu.va ;
154 
155  Valeur vuujm1(*mg) ;
156  if (uu.get_etat() == ETATZERO) {
157  vuujm1 = 1 ; // to take relative differences
158  vuujm1.set_base( vuu.base ) ;
159  }
160  else{
161  vuujm1 = vuu ;
162  }
163 
164  // Affine mapping for the Laplacian-tilde
165 
166  Map_af mpaff(*this) ;
167  Param par_nul ;
168 
169  cout << "Map_et::poisson_regular : relat. diff. u^J <-> u^{J-1} : "
170  << endl ;
171 
172 //==========================================================================
173 //==========================================================================
174 // Start of iteration
175 //==========================================================================
176 //==========================================================================
177 
178  Tbl tdiff(nz) ;
179  double diff ;
180  niter = 0 ;
181 
182  do {
183 
184  //====================================================================
185  // Computation of R(u) (the result is put in uu)
186  //====================================================================
187 
188 
189  //------------------------
190  // First operations on uu
191  //------------------------
192 
193  Valeur duudx = (uu.va).dsdx() ; // d/dx
194 
195  const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
196 
197  const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
198 
199 
200  //------------------
201  // Angular Laplacian
202  //------------------
203 
204  Valeur sxlapang = uu.va ;
205 
206  sxlapang.ylm() ;
207 
208  sxlapang = sxlapang.lapang() ;
209 
210  sxlapang = sxlapang.sx() ; // Lap_ang(uu) /x in the nucleus
211  // Lap_ang(uu) in the shells
212  // Lap_ang(uu) /(x-1) in the ZEC
213 
214  //------------------------------------------------------------------
215  // Computation of
216  // [ 2 /(dRdx) ( A - 1 ) duu/dx + 1/R (B - 1) Lap_ang(uu) ] / x
217  //
218  // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
219  // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
220  //
221  // The result is put in uu (via vuu)
222  //------------------------------------------------------------------
223 
224  Valeur varduudx = duudx ;
225 
226  if (zec) {
227  varduudx.annule(nzm1) ; // term in d/dx set to zero in the ZEC
228  }
229 
230  uu.set_etat_qcq() ;
231 
232  Base_val sauve_base = varduudx.base ;
233 
234  vuu = 2. * dxdr * ( rsxdxdr * unjj - 1.) * varduudx
235  + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
236 
237  vuu.set_base(sauve_base) ;
238 
239  vuu = vuu.sx() ;
240 
241  //----------------------------------------
242  // Computation of R(u)
243  //
244  // The result is put in uu (via vuu)
245  //----------------------------------------
246 
247  sauve_base = vuu.base ;
248 
249  vuu = xsr * vuu
250  + 2. * dxdr * ( sr2drdt * d2uudtdx
251  + sr2stdrdp * std2uudpdx ) ;
252 
253  vuu += dxdr * ( lapr_tp + dxdr * (
254  dxdr* unjj * d2rdx2
255  - 2. * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
256  ) * duudx ;
257 
258  vuu.set_base(sauve_base) ;
259 
260  // Since the assignment is performed on vuu (uu.va), the treatment
261  // of uu.dzpuis must be performed by hand:
262 
263  uu.set_dzpuis(4) ;
264 
265  if (source.get_dzpuis() == 2) {
266  uu.dec2_dzpuis() ; // uu.dzpuis: 4 -> 2
267  }
268 
269  if (source.get_dzpuis() == 3) {
270  uu.dec_dzpuis() ; //uu.dzpuis 4 -> 3
271  }
272 
273  //====================================================================
274  // Computation of the effective source s^J of the ``affine''
275  // Poisson equation
276  //====================================================================
277 
278  ssj = lambda * ssjm1 + unmlambda * ssjm2 ;
279 
280  ssj = ( source + uu + (amax - apre) * ssj ) / amax ;
281 
282  (ssj.va).set_base((source.va).base) ;
283 
284  //====================================================================
285  // Resolution of the ``affine'' Poisson equation
286  //====================================================================
287 
288  if ( source.get_dzpuis() == 0 ){
289  ssj.set_dzpuis( 4 ) ;
290  }
291  else {
292  ssj.set_dzpuis( source.get_dzpuis() ) ;
293  // Choice of the resolution
294  // dzpuis = 2, 3 or 4
295  }
296 
297  assert( uu.check_dzpuis( ssj.get_dzpuis() ) ) ;
298 
299  mpaff.poisson_regular(ssj, k_div, nzet, unsgam1, par_nul, uu,
300  uu_regu, uu_div, duu_div,
301  source_regu, source_div) ;
302 
303  //======================================
304  // Gradient of the diverging part (from that computed on the Map_af)
305  //======================================
306 
307  Valeur& dr_uu_div = duu_div.set(0).va ;
308  Valeur& dt_uu_div = duu_div.set(1).va ;
309  Valeur& dp_uu_div = duu_div.set(2).va ;
310 
311  Base_val bv = dr_uu_div.base ;
312  dr_uu_div = alpha[0] * dr_uu_div * dxdr ;
313  dr_uu_div.set_base( bv ) ;
314 
315  bv = dt_uu_div.base ;
316  dt_uu_div = alpha[0] * dt_uu_div * xsr - srdrdt * dr_uu_div ;
317  dt_uu_div.set_base( bv ) ;
318 
319  bv = dp_uu_div.base ;
320  dp_uu_div = alpha[0] * dp_uu_div * xsr - srstdrdp * dr_uu_div ;
321  dp_uu_div.set_base( bv ) ;
322 
323  duu_div.set_triad( this->get_bvect_spher() ) ;
324 
325 
326  //========================================
327  // Relative difference with previous step
328  //========================================
329 
330  tdiff = diffrel(vuu, vuujm1) ;
331 
332  diff = max(tdiff) ;
333 
334  cout << " step " << niter << " : " ;
335  for (int l=0; l<nz; l++) {
336  cout << tdiff(l) << " " ;
337  }
338  cout << endl ;
339 
340  //=================================
341  // Updates for the next iteration
342  //=================================
343 
344  ssjm2 = ssjm1 ;
345  ssjm1 = ssj ;
346  vuujm1 = vuu ;
347 
348  niter++ ;
349 
350  } // End of iteration
351  while ( (diff > precis) && (niter < nitermax) ) ;
352 
353 //==========================================================================
354 //==========================================================================
355 // End of iteration
356 //==========================================================================
357 //==========================================================================
358 
359 
360 }
361 }
Bases of the spectral expansions.
Definition: base_val.h:322
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
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:154
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
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:715
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:180
Affine radial mapping.
Definition: map.h:2027
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.
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 rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
Definition: map.h:2834
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2758
Coord d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1619
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1600
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1592
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1640
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1648
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1631
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1608
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1584
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1549
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1560
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition: map.h:676
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:474
Multi-domain array.
Definition: mtbl.h:118
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition: mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Parameter storage.
Definition: param.h:125
Cmp & get_cmp_mod(int position=0) const
Returns the reference of a modifiable Cmp stored in the list.
Definition: param.C:1049
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:361
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:430
Basic array class.
Definition: tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition: valeur_sx.C:110
const Valeur & stdsdp() const
Returns of *this.
Definition: valeur_stdsdp.C:60
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:112
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition: valeur.C:744
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:305
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:72
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene prototypes.
Definition: app_hor.h:64