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