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