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 
29 char 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 
52 namespace Lorene {
53 
54 void 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 d2rdx2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1619
Coord sr2stdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1608
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition: map.h:2758
Lorene prototypes.
Definition: app_hor.h:64
Coord sr2drdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1600
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Coord srstdrdp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1592
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1560
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Coord sstd2rdpdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1648
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1549
Coord rsxdxdr
in the nucleus; \ in the shells; \ in the outermost compactified domain.
Definition: map.h:2834
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:676
Coord lapr_tp
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1631
Coord srdrdt
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1584
Coord d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1640