LORENE
map_et_poisson_ylm.C
1/*
2 * Method of the class Map_et for the (iterative) resolution of the scalar
3 * Poisson equation with a multipole falloff condition at the outer boundary
4 *
5 * (see file map.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2004 Joshua A. Faber
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 version 2
16 * as published by the Free Software Foundation.
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
29char map_et_poisson_ylm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson_ylm.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $" ;
30
31/*
32 * $Id: map_et_poisson_ylm.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $
33 * $Log: map_et_poisson_ylm.C,v $
34 * Revision 1.2 2014/10/13 08:53:05 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.1 2004/12/30 15:56:42 k_taniguchi
38 * *** empty log message ***
39 *
40 *
41 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson_ylm.C,v 1.2 2014/10/13 08:53:05 j_novak Exp $
42 *
43 */
44
45// Lorene headers
46#include "map.h"
47#include "cmp.h"
48#include "param.h"
49
50//*****************************************************************************
51
52namespace Lorene {
53
54void Map_et::poisson_ylm(const Cmp& source, Param& par, Cmp& uu, int nylm, double* intvec) const {
55
56 assert(source.get_etat() != ETATNONDEF) ;
57 assert(source.get_mp() == this) ;
58
59 assert(uu.get_mp() == this) ;
60
61
62 int nz = mg->get_nzone() ;
63
64 //-------------------------------
65 // Computation of the prefactor a ---> Cmp apre
66 //-------------------------------
67
68 Mtbl unjj = 1 + srdrdt*srdrdt + srstdrdp*srstdrdp ;
69
70 Mtbl apre1(*mg) ;
71 apre1.set_etat_qcq() ;
72 for (int l=0; l<nz; l++) {
73 *(apre1.t[l]) = alpha[l]*alpha[l] ;
74 }
75
76 apre1 = apre1 * dxdr * dxdr * unjj ;
77
78
79 Cmp apre(*this) ;
80 apre = apre1 ;
81
82 Tbl amax0 = max(apre1) ; // maximum values in each domain
83
84 // The maximum values of a in each domain are put in a Mtbl
85 Mtbl amax1(*mg) ;
86 amax1.set_etat_qcq() ;
87 for (int l=0; l<nz; l++) {
88 *(amax1.t[l]) = amax0(l) ;
89 }
90
91 Cmp amax(*this) ;
92 amax = amax1 ;
93
94 //-------------------
95 // Initializations
96 //-------------------
97
98 int nitermax = par.get_int() ;
99 int& niter = par.get_int_mod() ;
100 double lambda = par.get_double() ;
101 lambda=1.0;
102 double unmlambda = 1. - lambda ;
103 double precis = par.get_double(1) ;
104
105 cout <<"relax in map_et:"<<lambda<<" "<<unmlambda<<endl;
106
107 Cmp& ssj = par.get_cmp_mod() ;
108
109 Cmp ssjm1 = ssj ;
110 Cmp ssjm2 = ssjm1 ;
111
112 Valeur& vuu = uu.va ;
113
114 Valeur vuujm1(*mg) ;
115 if (uu.get_etat() == ETATZERO) {
116 vuujm1 = 1 ; // to take relative differences
117 vuujm1.set_base( vuu.base ) ;
118 }
119 else{
120 vuujm1 = vuu ;
121 }
122
123 // Affine mapping for the Laplacian-tilde
124
125 Map_af mpaff(*this) ;
126 Param par_nul ;
127
128 cout << "Map_et::poisson : relat. diff. u^J <-> u^{J-1} : " << endl ;
129
130//==========================================================================
131//==========================================================================
132// Start of iteration
133//==========================================================================
134//==========================================================================
135
136 Tbl tdiff(nz) ;
137 double diff ;
138 niter = 0 ;
139
140 do {
141
142 //====================================================================
143 // Computation of R(u) (the result is put in uu)
144 //====================================================================
145
146
147 //-----------------------
148 // First operations on uu
149 //-----------------------
150
151 Valeur duudx = (uu.va).dsdx() ; // d/dx
152
153 const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
154
155 const Valeur& std2uudpdx = duudx.stdsdp() ; // 1/sin(theta) d^2/dxdphi
156
157 //------------------
158 // Angular Laplacian
159 //------------------
160
161 Valeur sxlapang = uu.va ;
162
163 sxlapang.ylm() ;
164
165 sxlapang = sxlapang.lapang() ;
166
167 sxlapang = sxlapang.sx() ; // Lap_ang(uu) /x in the nucleus
168 // Lap_ang(uu) in the shells
169 // Lap_ang(uu) /(x-1) in the ZEC
170
171 //---------------------------------------------------------------
172 // Computation of
173 // [ 2 /(dRdx) ( A - 1 ) duu/dx + 1/R (B - 1) Lap_ang(uu) ] / x
174 //
175 // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
176 // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
177 //
178 // The result is put in uu (via vuu)
179 //---------------------------------------------------------------
180
181 Valeur varduudx = duudx ;
182
183 uu.set_etat_qcq() ;
184
185 Base_val sauve_base = varduudx.base ;
186
187 vuu = 2. * dxdr * ( rsxdxdr * unjj - 1.) * varduudx
188 + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
189
190 vuu.set_base(sauve_base) ;
191
192 vuu = vuu.sx() ;
193
194 //---------------------------------------
195 // Computation of R(u)
196 //
197 // The result is put in uu (via vuu)
198 //---------------------------------------
199
200
201 sauve_base = vuu.base ;
202
203 vuu = xsr * vuu
204 + 2. * dxdr * ( sr2drdt * d2uudtdx
205 + sr2stdrdp * std2uudpdx ) ;
206
207 vuu += dxdr * ( lapr_tp + dxdr * (
208 dxdr* unjj * d2rdx2
209 - 2. * ( sr2drdt * d2rdtdx + sr2stdrdp * sstd2rdpdx ) )
210 ) * duudx ;
211
212 vuu.set_base(sauve_base) ;
213
214 // Since the assignment is performed on vuu (uu.va), the treatment
215 // of uu.dzpuis must be performed by hand:
216
217
218 //====================================================================
219 // Computation of the effective source s^J of the ``affine''
220 // Poisson equation
221 //====================================================================
222
223 ssj = lambda * ssjm1 + unmlambda * ssjm2 ;
224
225 ssj = ( source + uu + (amax - apre) * ssj ) / amax ;
226
227 (ssj.va).set_base((source.va).base) ;
228
229 //====================================================================
230 // Resolution of the ``affine'' Poisson equation
231 //====================================================================
232
233 // *****************************************************************
234
235 mpaff.poisson_ylm(ssj, par_nul, uu, nylm, intvec) ;
236
237 // *****************************************************************
238
239 tdiff = diffrel(vuu, vuujm1) ;
240
241 diff = max(tdiff) ;
242
243
244 cout << " iter: " << niter << " : " ;
245 for (int l=0; l<nz; l++) {
246 cout << tdiff(l) << " " ;
247 }
248 cout << endl ;
249
250 //=================================
251 // Updates for the next iteration
252 //=================================
253
254 ssjm2 = ssjm1 ;
255 ssjm1 = ssj ;
256 vuujm1 = vuu ;
257
258 niter++ ;
259
260 } // End of iteration
261 while ( (diff > precis) && (niter < nitermax) ) ;
262
263//==========================================================================
264//==========================================================================
265// End of iteration
266//==========================================================================
267//==========================================================================
268
269
270
271}
272
273
274}
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 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
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