LORENE
map_et_poisson2d.C
1/*
2 * Method of the class Map_et for the resolution of the 2-D Poisson
3 * equation.
4 *
5 * (see file map.h for documentation).
6 */
7
8/*
9 * Copyright (c) 2000-2001 Eric Gourgoulhon
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
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
30char map_et_poisson2d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson2d.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $" ;
31
32/*
33 * $Id: map_et_poisson2d.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $
34 * $Log: map_et_poisson2d.C,v $
35 * Revision 1.4 2014/10/13 08:53:05 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.3 2014/10/06 15:13:13 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.2 2002/02/07 14:55:58 e_gourgoulhon
42 * Corrected a bug when the source is known only in the coefficient
43 * space.
44 *
45 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
46 * LORENE
47 *
48 * Revision 2.4 2000/11/07 14:21:03 eric
49 * Correction d'une erreur dans le cas T_SIN_I (calcul de R(u)).
50 *
51 * Revision 2.3 2000/10/26 15:58:00 eric
52 * Correction cas T_COS_P : l'import de saff_q se fait par copie du Tbl.
53 *
54 * Revision 2.2 2000/10/12 15:37:43 eric
55 * Traitement des bases spectrales dans le cas T_COS_P.
56 *
57 * Revision 2.1 2000/10/11 15:15:43 eric
58 * 1ere version operationnelle.
59 *
60 * Revision 2.0 2000/10/09 13:47:17 eric
61 * *** empty log message ***
62 *
63 *
64 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_poisson2d.C,v 1.4 2014/10/13 08:53:05 j_novak Exp $
65 *
66 */
67
68// Headers C
69#include <cmath>
70
71// Headers Lorene:
72#include "map.h"
73#include "cmp.h"
74#include "param.h"
75
76//*****************************************************************************
77
78namespace Lorene {
79
80void Map_et::poisson2d(const Cmp& source_mat, const Cmp& source_quad,
81 Param& par, Cmp& uu) const {
82
83 assert(source_mat.get_etat() != ETATNONDEF) ;
84 assert(source_quad.get_etat() != ETATNONDEF) ;
85 assert(source_mat.get_mp()->get_mg() == mg) ;
86 assert(source_quad.get_mp()->get_mg() == mg) ;
87 assert(uu.get_mp()->get_mg() == mg) ;
88
89 assert( source_quad.check_dzpuis(4) ) ;
90
91 double& lambda = par.get_double_mod(0) ;
92 int nz = mg->get_nzone() ;
93 int nzm1 = nz-1 ;
94
95 // Special case of a vanishing source
96 // ----------------------------------
97
98 if ( (source_mat.get_etat() == ETATZERO)
99 && (source_quad.get_etat() == ETATZERO) ) {
100
101 uu = 0 ;
102 lambda = 1 ;
103 return ;
104 }
105
106 int base_t = ((source_mat.va).base).get_base_t(0) ;
107
108 switch (base_t) {
109
110 //==================================================================
111 // case T_COS_P
112 //==================================================================
113
114 case T_COS_P : {
115
116 // Construction of a Map_af which coincides with *this on the equator
117 // ------------------------------------------------------------------
118
119 double theta0 = M_PI / 2 ; // Equator
120 double phi0 = 0 ;
121
122 Map_af mpaff(*this) ;
123
124 for (int l=0 ; l<nz ; l++) {
125 double rmax = val_r(l, double(1), theta0, phi0) ;
126 switch ( mg->get_type_r(l) ) {
127 case RARE: {
128 double rmin = val_r(l, double(0), theta0, phi0) ;
129 mpaff.set_alpha(rmax - rmin, l) ;
130 mpaff.set_beta(rmin, l) ;
131 break ;
132 }
133
134 case FIN: {
135 double rmin = val_r(l, double(-1), theta0, phi0) ;
136 mpaff.set_alpha( double(.5) * (rmax - rmin), l ) ;
137 mpaff.set_beta( double(.5) * (rmax + rmin), l) ;
138 break ;
139 }
140
141 case UNSURR: {
142 double rmin = val_r(l, double(-1), theta0, phi0) ;
143 double umax = double(1) / rmin ;
144 double umin = double(1) / rmax ;
145 mpaff.set_alpha( double(.5) * (umin - umax), l) ;
146 mpaff.set_beta( double(.5) * (umin + umax), l) ;
147 break ;
148 }
149
150 default: {
151 cout << "Map_et::poisson2d: unknown type_r ! " << endl ;
152 abort () ;
153 break ;
154 }
155
156 }
157 }
158
159 // Importation of source_mat and source_quad of the affine mapping
160 // ---------------------------------------------------------------
161 Cmp saff_m(mpaff) ;
162 saff_m.import( nzm1, source_mat ) ;
163 (saff_m.va).set_base( (source_mat.va).base ) ;
164
165 Cmp saff_q(mpaff) ;
166
167 // In order to use Cmp::import with dzpuis != 0 :
168 Cmp set_q = source_quad ;
169 set_q.set_dzpuis(0) ; // dzpuis artificially set to 0
170
171 saff_q.import( nzm1, set_q ) ;
172 (saff_q.va).set_base( (set_q.va).base ) ;
173
174 // Copy in the external domain :
175 if ( (set_q.va).get_etat() == ETATQCQ) {
176 (set_q.va).coef_i() ; // the values in configuration space are required
177 assert( (set_q.va).c->get_etat() == ETATQCQ ) ;
178 assert( (saff_q.va).c->get_etat() == ETATQCQ ) ;
179 *( (saff_q.va).c->t[nzm1] ) = *( (set_q.va).c->t[nzm1] ) ;
180 }
181
182 // the true dzpuis is restored :
183 saff_q.set_dzpuis( source_quad.get_dzpuis() ) ;
184
185 // Resolution of the 2-D Poisson equation on the spherical domains
186 // ---------------------------------------------------------------
187
188 Cmp uaff(mpaff) ;
189
190 mpaff.poisson2d(saff_m, saff_q, par, uaff) ;
191
192 // Importation of the solution on the Map_et mapping *this
193 // -------------------------------------------------------
194
195 uu.import(uaff) ;
196
197 uu.va.set_base( uaff.va.base ) ; // same spectral bases
198
199 break ;
200 }
201
202 //==================================================================
203 // case T_SIN_I
204 //==================================================================
205
206 case T_SIN_I : {
207
208 //-------------------------------
209 // Computation of the prefactor a ---> Cmp apre
210 //-------------------------------
211
212 Mtbl unjj = 1 + srdrdt*srdrdt ;
213
214 Mtbl apre1(*mg) ;
215 apre1.set_etat_qcq() ;
216 for (int l=0; l<nz; l++) {
217 *(apre1.t[l]) = alpha[l]*alpha[l] ;
218 }
219
220 apre1 = apre1 * dxdr * dxdr * unjj ;
221
222 Cmp apre(*this) ;
223 apre = apre1 ;
224
225 Tbl amax0 = max(apre1) ; // maximum values in each domain
226
227 // The maximum values of a in each domain are put in a Mtbl
228 Mtbl amax1(*mg) ;
229 amax1.set_etat_qcq() ;
230 for (int l=0; l<nz; l++) {
231 *(amax1.t[l]) = amax0(l) ;
232 }
233
234 Cmp amax(*this) ;
235 amax = amax1 ;
236
237
238 //-------------------
239 // Initializations
240 //-------------------
241
242 int nitermax = par.get_int() ;
243 int& niter = par.get_int_mod() ;
244 double lambda_relax = par.get_double() ;
245 double unmlambda_relax = 1. - lambda_relax ;
246 double precis = par.get_double(1) ;
247
248 Cmp& ssj = par.get_cmp_mod() ;
249
250 Cmp ssjm1 = ssj ;
251 Cmp ssjm2 = ssjm1 ;
252
253 Cmp ssj_q(*this) ;
254 ssj_q = 0 ;
255
256 Valeur& vuu = uu.va ;
257
258 Valeur vuujm1(*mg) ;
259 if (uu.get_etat() == ETATZERO) {
260 vuujm1 = 1 ; // to take relative differences
261 vuujm1.set_base( vuu.base ) ;
262 }
263 else{
264 vuujm1 = vuu ;
265 }
266
267 // Affine mapping for the Laplacian-tilde
268
269 Map_af mpaff(*this) ;
270
271 cout << "Map_et::poisson2d : relat. diff. u^J <-> u^{J-1} : " << endl ;
272
273//==========================================================================
274//==========================================================================
275// Start of iteration
276//==========================================================================
277//==========================================================================
278
279 Tbl tdiff(nz) ;
280 double diff ;
281 niter = 0 ;
282
283 do {
284
285 //====================================================================
286 // Computation of R(u) (the result is put in uu)
287 //====================================================================
288
289
290 //-----------------------
291 // First operations on uu
292 //-----------------------
293
294 Valeur duudx = (uu.va).dsdx() ; // d/dx
295
296 const Valeur& d2uudtdx = duudx.dsdt() ; // d^2/dxdtheta
297
298 //-------------------
299 // 1/x d^2uu/dtheta^2
300 //-------------------
301
302 Valeur sxlapang = uu.va ;
303
304 sxlapang = sxlapang.d2sdt2() ;
305
306 sxlapang = sxlapang.sx() ; // d^2(uu)/dth^2 /x in the nucleus
307 // d^2(uu)/dth^2 in the shells
308 // d^2(uu)/dth^2 /(x-1) in the ZEC
309
310 //---------------------------------------------------------------
311 // Computation of
312 // [ (dR/dx)^{-1} ( A - 1 ) duu/dx + 1/R (B - 1) d^2uu/dth^2 ] / x
313 //
314 // with A = 1/(dRdx) R/(x+beta_l/alpha_l) unjj
315 // B = [1/(dRdx) R/(x+beta_l/alpha_l)]^2 unjj
316 //
317 // The result is put in uu (via vuu)
318 //---------------------------------------------------------------
319
320 // Intermediate quantity jac which value is
321 // (dR/dx)^{-1} in the nucleus and the shells
322 // +(dU/dx)^{-1} in the ZEC
323
324 Mtbl jac = dxdr ;
325 if (mg->get_type_r(nzm1) == UNSURR) {
326 jac.annule(nzm1, nzm1) ;
327 Mtbl jac_ext = dxdr ;
328 jac_ext.annule(0, nzm1-1) ;
329 jac_ext = - jac_ext ;
330 jac = jac + jac_ext ;
331 }
332
333 uu.set_etat_qcq() ;
334
335 Base_val sauve_base = duudx.base ;
336
337 vuu = jac * ( rsxdxdr * unjj - 1.) * duudx
338 + ( rsxdxdr*rsxdxdr * unjj - 1.) * xsr * sxlapang ;
339
340 vuu.set_base(sauve_base) ;
341
342 vuu = vuu.sx() ;
343
344 //---------------------------------------
345 // Computation of R(u)
346 //
347 // The result is put in uu (via vuu)
348 //---------------------------------------
349
350
351 sauve_base = vuu.base ;
352
353 vuu = xsr * vuu
354 + 2. * dxdr * sr2drdt * d2uudtdx ;
355
356 vuu += dxdr * ( sr2d2rdt2 + dxdr * (
357 dxdr* unjj * d2rdx2
358 - 2. * sr2drdt * d2rdtdx )
359 ) * duudx ;
360
361 vuu.set_base(sauve_base) ;
362
363 // Since the assignment is performed on vuu (uu.va), the treatment
364 // of uu.dzpuis must be performed by hand:
365
366 uu.set_dzpuis(4) ;
367
368 //====================================================================
369 // Computation of the effective source s^J of the ``affine''
370 // Poisson equation
371 //====================================================================
372
373 ssj = lambda_relax * ssjm1 + unmlambda_relax * ssjm2 ;
374
375 ssj = ( source_mat + source_quad + uu + (amax - apre) * ssj ) / amax ;
376
377 (ssj.va).set_base((source_mat.va).base) ;
378
379 //====================================================================
380 // Resolution of the ``affine'' Poisson equation
381 //====================================================================
382
383 assert( uu.check_dzpuis( ssj.get_dzpuis() ) ) ;
384
385 mpaff.poisson2d(ssj, ssj_q, par, uu) ;
386
387 tdiff = diffrel(vuu, vuujm1) ;
388
389 diff = max(tdiff) ;
390
391
392 cout << " step " << niter << " : " ;
393 for (int l=0; l<nz; l++) {
394 cout << tdiff(l) << " " ;
395 }
396 cout << endl ;
397
398 //=================================
399 // Updates for the next iteration
400 //=================================
401
402 ssjm2 = ssjm1 ;
403 ssjm1 = ssj ;
404 vuujm1 = vuu ;
405
406 niter++ ;
407
408 } // End of iteration
409 while ( (diff > precis) && (niter < nitermax) ) ;
410
411//==========================================================================
412//==========================================================================
413// End of iteration
414//==========================================================================
415//==========================================================================
416
417
418 break ;
419 }
420
421 default : {
422 cout << "Map_et::poisson2d : unkown theta basis !" << endl ;
423 cout << " basis : " << hex << base_t << endl ;
424 abort() ;
425 break ;
426 }
427 } // End of switch on base_t
428}
429
430
431}
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 import(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
Definition cmp_import.C:73
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
Affine radial mapping.
Definition map.h:2027
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
Definition map_af.C:641
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const
Computes the solution of a 2-D Poisson equation.
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Definition map_af.C:630
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const
Computes the solution of a 2-D 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 d2rdtdx
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1640
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
Coord sr2d2rdt2
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1657
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
void annule(int l_min, int l_max)
Sets the Mtbl to zero in some domains.
Definition mtbl.C:329
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
double & get_double_mod(int position=0) const
Returns the reference of a stored modifiable double .
Definition param.C:498
Basic array class.
Definition tbl.h:161
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 & d2sdt2() const
Returns of *this.
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
const Valeur & dsdt() const
Returns of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
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
#define T_COS_P
dev. cos seulement, harmoniques paires
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Lorene prototypes.
Definition app_hor.h:64