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
31char 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
83namespace Lorene {
84
85void 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...
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
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
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
void dec2_dzpuis()
Decreases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
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.
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.
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.
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