LORENE
map_radial_poisson_cpt.C
1/*
2 * Method of the class Map_radial for resolution of the equation
3 *
4 * a \Delta\psi + {\bf b} \cdot \nabla \psi = \sigma
5 *
6 * Computes the unique solution such that psi(0) = 0.
7 * bb must be given in spherical coordinates.
8 *
9 * (see file map.h for documentation)
10 *
11 */
12
13/*
14 * Copyright (c) 2000-2007 Eric Gourgoulhon
15 * Copyright (c) 2007 Michal Bejger
16 * Copyright (c) 2000-2001 Philippe Grandclement
17 * Copyright (c) 2001 Keisuke Taniguchi
18 *
19 * This file is part of LORENE.
20 *
21 * LORENE is free software; you can redistribute it and/or modify
22 * it under the terms of the GNU General Public License as published by
23 * the Free Software Foundation; either version 2 of the License, or
24 * (at your option) any later version.
25 *
26 * LORENE is distributed in the hope that it will be useful,
27 * but WITHOUT ANY WARRANTY; without even the implied warranty of
28 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29 * GNU General Public License for more details.
30 *
31 * You should have received a copy of the GNU General Public License
32 * along with LORENE; if not, write to the Free Software
33 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
34 *
35 */
36
37
38
39
40/*
41 * $Id: map_radial_poisson_cpt.C,v 1.8 2016/12/05 16:17:58 j_novak Exp $
42 * $Log: map_radial_poisson_cpt.C,v $
43 * Revision 1.8 2016/12/05 16:17:58 j_novak
44 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
45 *
46 * Revision 1.7 2014/10/13 08:53:06 j_novak
47 * Lorene classes and functions now belong to the namespace Lorene.
48 *
49 * Revision 1.6 2014/10/06 15:13:14 j_novak
50 * Modified #include directives to use c++ syntax.
51 *
52 * Revision 1.5 2007/10/18 08:19:32 e_gourgoulhon
53 * Suppression of the abort for nzet > 2 : the function should be able
54 * to treat an arbitrary number of domains inside the star.
55 * NB: tested only for nzet = 2.
56 *
57 * Revision 1.4 2007/10/16 21:53:13 e_gourgoulhon
58 * Added new function poisson_compact (multi-domain version)
59 *
60 * Revision 1.3 2003/10/03 15:58:48 j_novak
61 * Cleaning of some headers
62 *
63 * Revision 1.2 2002/12/11 13:17:07 k_taniguchi
64 * Change the multiplication "*" to "%".
65 *
66 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
67 * LORENE
68 *
69 * Revision 2.17 2001/08/28 15:08:06 keisuke
70 * Uncomment "sour_j.std_base_scal()" and "sour_j.ylm()".
71 *
72 * Revision 2.16 2001/08/28 14:45:06 keisuke
73 * Uncomment "sour_j.coef()".
74 *
75 * Revision 2.15 2001/08/28 14:10:42 keisuke
76 * Comment out "sour_j.std_base_scal()", "sour_j.coef()", and "sour_j.ylm()".
77 *
78 * Revision 2.14 2000/03/30 09:18:53 eric
79 * Modifs affichage.
80 *
81 * Revision 2.13 2000/03/17 15:24:01 phil
82 * suppression du nettoyage brutal
83 *
84 * Revision 2.12 2000/03/16 16:26:07 phil
85 * Ajout du nettoyage des hautes frequences
86 *
87 * Revision 2.11 2000/03/10 09:18:36 eric
88 * Annulation du terme l=0 de la source effective.
89 *
90 * Revision 2.10 2000/03/09 13:52:19 phil
91 * *** empty log message ***
92 *
93 * Revision 2.9 2000/02/28 14:34:31 eric
94 * *** empty log message ***
95 *
96 * Revision 2.8 2000/02/28 14:31:32 eric
97 * Suppression de mpaff: remplacement de pre_lap = (1-r^2/R^2) par
98 * Id - .mult_x().mult_x().
99 * Introduction de dpsi.
100 * Modif affichages.
101 *
102 * Revision 2.7 2000/02/25 13:48:00 eric
103 * Suppression des appels a Mtbl_cf::nettoie().
104 *
105 * Revision 2.6 2000/02/22 11:43:10 eric
106 * Appel de ylm() sur d2psi.
107 * Modif affichage.
108 *
109 * Revision 2.5 2000/02/21 12:58:12 eric
110 * Modif affichage.
111 *
112 * Revision 2.4 2000/01/27 15:58:36 eric
113 * Utilisation du nouveau constructeur Map_af::Map_af(const Map&)
114 * pour la construction du mapping auxiliaire mpaff.
115 * Suppression des affichages intermediaires.
116 *
117 * Revision 2.3 2000/01/14 17:33:44 eric
118 * Amelioration du calcul (le Cmp intermediaire psi0 est supprime).
119 *
120 * Revision 2.2 2000/01/13 16:43:59 eric
121 * Premiere version testee : OK !
122 *
123 * Revision 2.1 2000/01/12 16:03:23 eric
124 * Premiere version complete.
125 *
126 * Revision 2.0 2000/01/10 09:14:52 eric
127 * *** empty log message ***
128 *
129 *
130 * $Header: /cvsroot/Lorene/C++/Source/Map/map_radial_poisson_cpt.C,v 1.8 2016/12/05 16:17:58 j_novak Exp $
131 *
132 */
133
134// Headers C++
135
136// Headers C
137#include <cstdlib>
138#include <cmath>
139
140// Headers Lorene
141#include "tenseur.h"
142#include "param.h"
143#include "proto.h"
144#include "graphique.h"
145#include "utilitaires.h"
146
147// Local prototypes
148namespace Lorene {
149Mtbl_cf sol_poisson_compact(const Mtbl_cf&, double, double, bool) ;
150Mtbl_cf sol_poisson_compact(const Map_af&, const Mtbl_cf&, const Tbl&,
151 const Tbl&, bool) ;
152
153
155 // Single domain version //
157
158void Map_radial::poisson_compact(const Cmp& source, const Cmp& aa,
159 const Tenseur& bb, const Param& par,
160 Cmp& psi) const {
161
162
163 // Protections
164 // -----------
165
166 assert(source.get_etat() != ETATNONDEF) ;
167 assert(aa.get_etat() != ETATNONDEF) ;
168 assert(bb.get_etat() != ETATNONDEF) ;
169 assert(aa.get_mp() == source.get_mp()) ;
170 assert(bb.get_mp() == source.get_mp()) ;
171 assert(psi.get_mp() == source.get_mp()) ;
172
173
174 // The components of vector b must be given in the spherical basis
175 // associated with the mapping :
176 assert(*(bb.get_triad()) == bvect_spher) ;
177
178 // Computation parameters
179 // ----------------------
180 int nitermax = par.get_int() ;
181 int& niter = par.get_int_mod() ;
182 double precis = par.get_double(0) ;
183 double relax = par.get_double(1) ;
184 double unmrelax = 1. - relax ;
185
186 // Maybe nothing to do ?
187 // ---------------------
188
189 if ( source.get_etat() == ETATZERO ) {
190 psi.set_etat_zero() ;
191 return ;
192 }
193
194 // Factors alpha ( in front of (1-xi^2) Lap_xi(psi) )
195 // and beta ( in front of xi dpsi/dxi )
196 // --------------------------------------------------
197
198 Mtbl tmp = dxdr ;
199 double dxdr_c = tmp(0, 0, 0, 0) ; // central value of 1/(dR/dxi)
200
201 double alph = dxdr_c * dxdr_c * aa(0, 0, 0, 0) ;
202
203 const Valeur& b_r = bb(0).va ; // radial component b^r of vector b
204
205 Valeur vatmp = dxdr*dxdr*b_r ;
206
207 double bet = min(vatmp)(0) ; // Minimal value in domain no. 0
208
209 cout << "Map_radial::poisson_compact : alph, bet : " << alph << " "
210 << bet << endl ;
211
212
213 // Everything is set to zero except in the innermost domain (nucleus) :
214 // ------------------------------------------------------------------
215
216 int nz = mg->get_nzone() ;
217
218 psi.annule(1, nz-1) ;
219
220 // Auxilary quantities
221 // -------------------
222 Cmp psi_jm1 = psi ;
223 Cmp b_grad_psi(this) ;
224 Valeur sour_j(*mg) ;
225 Valeur aux_psi(*mg) ;
226 Valeur lap_xi_psi(*mg) ;
227 Valeur oper_psi(*mg) ;
228 Valeur dpsi(*mg) ;
229 Valeur d2psi(*mg) ;
230
231 Valeur& vpsi = psi.va ;
232
233//==========================================================================
234// Start of iteration
235//==========================================================================
236
237 Tbl tdiff(nz) ;
238 double diff ;
239 niter = 0 ;
240
241
242 do {
243
244 // Computation of the source for sol_poisson_compact
245 // -------------------------------------------------
246
247 b_grad_psi = bb(0) % psi.dsdr() +
248 bb(1) % psi.srdsdt() +
249 bb(2) % psi.srstdsdp() ;
250
251
252 vpsi.ylm() ; // Expansion of psi onto spherical harmonics
253
254 // Lap_xi(psi) :
255
256 dpsi = vpsi.dsdx() ;
257 dpsi.ylm() ; // necessary because vpsi.p_dsdx
258 // has been already computed (in
259 // non-spherical harmonics bases) by
260 // the call psi.dsdr() above
261
262 aux_psi = 2.*dpsi + (vpsi.lapang()).sx() ;
263
264 d2psi = vpsi.d2sdx2() ;
265 d2psi.ylm() ;
266
267 lap_xi_psi = d2psi + aux_psi.sx() ;
268
269 // Effective source :
270
271 sour_j = source.va
272 + alph * ( lap_xi_psi - (lap_xi_psi.mult_x()).mult_x() )
273 - aa.va * (psi.laplacien()).va
274 + bet * dpsi.mult_x()
275 - b_grad_psi.va ;
276
277 sour_j.std_base_scal() ;
278 sour_j.coef() ;
279 sour_j.ylm() ; // Spherical harmonics expansion
280
281 // The term l=0 of the effective source is set to zero :
282 // ---------------------------------------------------
283 double somlzero = 0 ;
284 for (int i=0; i<mg->get_nr(0); i++) {
285 somlzero += fabs( (*(sour_j.c_cf))(0, 0, 0, i) ) ;
286 (sour_j.c_cf)->set(0, 0, 0, i) = 0 ;
287 }
288
289 if (somlzero > 1.e-10) {
290 cout << "### WARNING : Map_radial::poisson_compact : " << endl
291 << " the l=0 part of the effective source is > 1.e-10 : "
292 << somlzero << endl ;
293 }
294
295
296 // Resolution of the equation
297 // a (1-xi^2) Lap_xi(psi) + b xi dpsi/dxi = sour_j
298 //---------------------------------------------------
299
300 bool reamorce = (niter == 0) ;
301
302 assert(sour_j.c_cf != 0x0) ;
303
304 psi.set_etat_zero() ; // to call Cmp::del_deriv().
305 psi.set_etat_qcq() ;
306 vpsi = sol_poisson_compact( *(sour_j.c_cf), alph, bet, reamorce ) ;
307
308
309 // Test: multiplication by the operator matrix
310 // -------------------------------------------
311
312// Mtbl_cf cvpsi = *(vpsi.c_cf) ;
313// Mtbl_cf csour = *(sour_j.c_cf) ;
314//
315// int np = mg->get_np(0) ;
316// int nt = mg->get_nt(0) ;
317// int nr = mg->get_nr(0) ;
318//
319// for (int k=0 ; k<np+1 ; k++) {
320// for (int j=0 ; j<nt ; j++) {
321// if (nullite_plm(j, nt, k, np, cvpsi.base) == 1) {
322// int l_quant, m_quant, base_r ;
323// donne_lm(nz, 0, j, k, cvpsi.base, m_quant, l_quant, base_r) ;
324//
325// cout << "### k, j, l, m : " << k << " " << j << " "
326// << l_quant << " " << m_quant << endl ;
327// Matrice operateur = alph * laplacien_nul_mat(nr, l_quant, base_r)
328// + bet * xdsdx_mat(nr, l_quant, base_r) ;
329// Tbl coef(nr) ;
330// operateur = combinaison_nul(operateur, l_quant, coef, base_r) ;
331//
332// Tbl so(nr) ;
333// so.set_etat_qcq() ;
334// for (int i=0 ; i<nr ; i++)
335// so.set(i) = csour(0, k, j, i) ;
336// so = combinaison_nul(so, l_quant, coef, base_r) ;
337//
338// Tbl psi1(nr) ;
339// Tbl opi1(nr) ;
340// psi1.set_etat_qcq() ;
341// opi1.set_etat_qcq() ;
342//
343// bool discrep = false ;
344//
345// for (int i=0; i<nr; i++) {
346//
347// double op = 0 ;
348// for (int i1=0; i1<nr; i1++) {
349// op += operateur(i, i1) * cvpsi(0, k, j, i1) ;
350// }
351//
352// psi1.set(i) = cvpsi(0, k, j, i) ;
353// opi1.set(i) = op ;
354//
355// double x1 = so(i) ;
356// double x2 = op - so(i) ;
357// double seuil = 1e-11 ;
358// if ( fabs(x2) > seuil )
359// {
360// discrep = true ;
361// cout << "i : op , sou, diff : " << i << " : " << op << " "
362// << x1 << " " << x2 << endl ;
363// }
364//
365// }
366//
367// if ( discrep ) {
368// cout << "Matrice pour k, j = " << k << " " << j << endl ;
369// cout << operateur << endl ;
370// cout << "psi : " << psi1 << endl ;
371// cout << "op(psi) : " << opi1 << endl ;
372// cout << "so : " << so << endl ;
373// }
374//
375// } // fin du cas ou m<=l
376// } // fin de boucle sur j
377// } // fin de boucle sur k
378
379
380 // Test: has the equation been correctly solved ?
381 // ---------------------------------------------
382
383 // Lap_xi(psi) :
384 aux_psi = vpsi.dsdx() ;
385
386 aux_psi = 2.*aux_psi + (vpsi.lapang()).sx() ;
387
388 lap_xi_psi = vpsi.d2sdx2() + aux_psi.sx() ;
389
390 oper_psi = alph * ( lap_xi_psi - (lap_xi_psi.mult_x()).mult_x() )
391 + bet * (vpsi.dsdx()).mult_x() ;
392
393
394 double maxc = fabs( max(*(vpsi.c_cf))(0) ) ;
395 double minc = fabs( min(*(vpsi.c_cf))(0) ) ;
396 double max_abs_psi = ( maxc > minc ) ? maxc : minc ;
397
398 maxc = fabs( max(*(sour_j.c_cf))(0) ) ;
399 minc = fabs( min(*(sour_j.c_cf))(0) ) ;
400 double max_abs_sou = ( maxc > minc ) ? maxc : minc ;
401
402 Mtbl_cf diff_opsou = *oper_psi.c_cf - *sour_j.c_cf ;
403 maxc = fabs( max(diff_opsou)(0) ) ;
404 minc = fabs( min(diff_opsou)(0) ) ;
405 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
406
407
408// cout << " Coef of oper_psi - sour_j : " << endl ;
409// diff_opsou.affiche_seuil(cout, 4, 1.e-11) ;
410
411 cout << " Step " << niter
412 << " : test (1-xi^2) Lap_xi(psi) + b xi dpsi/dxi : " << endl ;
413 cout << " max |psi| |sou| |oper(psi)-sou|: "
414 << max_abs_psi << " " << max_abs_sou << " "
415 << max_abs_diff << endl ;
416
417 // Relaxation :
418 // -----------
419
420 vpsi.ylm_i() ; // Inverse spherical harmonics transform
421
422 psi = relax * psi + unmrelax * psi_jm1 ;
423
424 tdiff = diffrel(psi, psi_jm1) ;
425
426 diff = max(tdiff) ;
427
428 cout <<
429 " Relative difference psi^J <-> psi^{J-1} : "
430 << tdiff(0) << endl ;
431
432// arrete() ;
433
434 // Updates for the next iteration
435 // -------------------------------
436
437 psi_jm1 = psi ;
438 niter++ ;
439
440 } // End of iteration
441 while ( (diff > precis) && (niter < nitermax) ) ;
442
443//==========================================================================
444// End of iteration
445//==========================================================================
446
447}
448
449
450
452 // Multi domain version //
454
455
456void Map_radial::poisson_compact(int nzet, const Cmp& source, const Cmp& aa,
457 const Tenseur& bb, const Param& par,
458 Cmp& psi) const {
459
460 if (nzet == 1) {
461 poisson_compact(source, aa, bb, par, psi) ;
462 return ;
463 }
464
465
466 // Protections
467 // -----------
468
469 assert(source.get_etat() != ETATNONDEF) ;
470 assert(aa.get_etat() != ETATNONDEF) ;
471 assert(bb.get_etat() != ETATNONDEF) ;
472 assert(aa.get_mp() == source.get_mp()) ;
473 assert(bb.get_mp() == source.get_mp()) ;
474 assert(psi.get_mp() == source.get_mp()) ;
475
476
477 // The components of vector b must be given in the spherical basis
478 // associated with the mapping :
479 assert(*(bb.get_triad()) == bvect_spher) ;
480
481 // Maybe nothing to do ?
482 // ---------------------
483
484 if ( source.get_etat() == ETATZERO ) {
485 psi.set_etat_zero() ;
486 return ;
487 }
488
489 // Computation parameters
490 // ----------------------
491 int nitermax = par.get_int() ;
492 int& niter = par.get_int_mod() ;
493 double precis = par.get_double(0) ;
494 double relax = par.get_double(1) ;
495 double unmrelax = 1. - relax ;
496
497 // Auxiliary affine mapping
498 Map_af mpaff(*this) ;
499
500 // Coefficients to fit the profiles of aa and bb in each domain
501 // ------------------------------------------------------------
502 Tbl ac(nzet,3) ;
503 ac.annule_hard() ; // initialization to zero
504 Tbl bc(nzet,3) ;
505 bc.annule_hard() ; // initialization to zero
506
507 Valeur ap = aa.va ;
508 Valeur bp = bb(0).va ;
509
510 // Coefficients in the nucleus
511 int nrm1 = mg->get_nr(0) - 1 ;
512 ac.set(0,0) = ap(0,0,0,0) ;
513 ac.set(0,2) = ap(0,0,0,nrm1) - ac(0,0) ;
514
515 bc.set(0,1) = bp(0,0,0,nrm1) ;
516
517 // Coefficients in the intermediate shells
518 for (int lz=1; lz<nzet-1; lz++) {
519 nrm1 = mg->get_nr(lz) - 1 ;
520 ac.set(lz,0) = 0.5 * ( ap(lz,0,0,nrm1) + ap(lz,0,0,0) ) ;
521 ac.set(lz,1) = 0.5 * ( ap(lz,0,0,nrm1) - ap(lz,0,0,0) ) ;
522
523 bc.set(lz,0) = 0.5 * ( bp(lz,0,0,nrm1) + bp(lz,0,0,0) ) ;
524 bc.set(lz,1) = 0.5 * ( bp(lz,0,0,nrm1) - bp(lz,0,0,0) ) ;
525 }
526
527 // Coefficients in the external shell
528 int lext = nzet - 1 ;
529 nrm1 = mg->get_nr(lext) - 1 ;
530 ac.set(lext,0) = 0.5 * ap(lext,0,0,0) ;
531 ac.set(lext,1) = - ac(lext,0) ;
532
533 bc.set(lext,0) = 0.5 * ( bp(lext,0,0,nrm1) + bp(lext,0,0,0) ) ;
534 bc.set(lext,1) = 0.5 * ( bp(lext,0,0,nrm1) - bp(lext,0,0,0) ) ;
535
536 cout << "ac : " << ac << endl ;
537 cout << "bc : " << bc << endl ;
538
539 // Prefactor of Lap_xi(Psi) and dPsi/dxi
540 // -------------------------------------
541
542 Mtbl ta(mg) ;
543 Mtbl tb(mg) ;
544 ta.annule_hard() ;
545 tb.annule_hard() ;
546 for (int lz=0; lz<nzet; lz++) {
547 const double* xi = mg->get_grille3d(lz)->x ;
548 double* tta = ta.set(lz).t ;
549 double* ttb = tb.set(lz).t ;
550 int np = mg->get_np(lz) ;
551 int nt = mg->get_nt(lz) ;
552 int nr = mg->get_nr(lz) ;
553 int pt = 0 ;
554 for (int k=0; k<np; k++) {
555 for (int j=0; j<nt; j++) {
556 for (int i=0; i<nr; i++) {
557 tta[pt] = ac(lz,0) + xi[i] * (ac(lz,1) + ac(lz,2) * xi[i]) ;
558 ttb[pt] = bc(lz,0) + xi[i] * (bc(lz,1) + bc(lz,2) * xi[i]) ;
559 pt++ ;
560 }
561 }
562 }
563 }
564
565// Verification
566// cout << "Map :" << *(aa.get_mp()) << endl ;
567// Cmp tverif(*this) ;
568// tverif = ap ;
569// tverif.std_base_scal() ;
570// des_profile(tverif,0., 4., 0., 0.) ;
571// tverif = ta ;
572// tverif.std_base_scal() ;
573// des_profile(tverif,0., 4., 0., 0.) ;
574//
575// tverif = bp ;
576// tverif.std_base_scal() ;
577// des_profile(tverif,0., 4., 0., 0.) ;
578// tverif = tb ;
579// tverif.std_base_scal() ;
580// des_profile(tverif,0., 4., 0., 0.) ;
581
582
583 // Everything is set to zero except inside the star
584 // -------------------------------------------------
585
586 int nz = mg->get_nzone() ;
587
588 psi.annule(nzet, nz-1) ;
589
590 // Auxilary quantities
591 // -------------------
592 Cmp psi_jm1 = psi ;
593 Cmp b_grad_psi(this) ;
594 Valeur sour_j(*mg) ;
595 Valeur aux_psi(*mg) ;
596 Valeur lap_xi_psi(*mg) ;
597 Valeur oper_psi(*mg) ;
598 Valeur dpsi(*mg) ;
599 Valeur d2psi(*mg) ;
600
601 Valeur& vpsi = psi.va ;
602
603//==========================================================================
604// Start of iteration
605//==========================================================================
606
607 Tbl tdiff(nz) ;
608 double diff ;
609 niter = 0 ;
610
611
612 do {
613
614 // Computation of the source for sol_poisson_compact
615 // -------------------------------------------------
616
617 b_grad_psi = bb(0) % psi.dsdr() + bb(1) % psi.srdsdt() + bb(2) % psi.srstdsdp() ;
618
619//??
620 vpsi.ylm() ; // Expansion of psi onto spherical harmonics
621
622 // Effective source :
623
624 Cmp lap_zeta(mpaff) ;
625 mpaff.laplacien(psi, 0, lap_zeta) ;
626
627 Cmp grad_zeta(mpaff) ;
628 mpaff.dsdr(psi, grad_zeta) ;
629
630 sour_j = source.va
631 + ta * lap_zeta.va - aa.va * (psi.laplacien()).va
632 + tb * grad_zeta.va - b_grad_psi.va ;
633
634 sour_j.std_base_scal() ;
635 sour_j.coef() ;
636 sour_j.ylm() ; // Spherical harmonics expansion
637
638
639 // The term l=0 of the effective source is set to zero :
640 // ---------------------------------------------------
641
642 for (int lz=0; lz<nzet; lz++) {
643 double somlzero = 0 ;
644 for (int i=0; i<mg->get_nr(lz); i++) {
645 somlzero += fabs( (*(sour_j.c_cf))(lz, 0, 0, i) ) ;
646 (sour_j.c_cf)->set(lz, 0, 0, i) = 0 ;
647 }
648 if (somlzero > 1.e-10) {
649 cout << "### WARNING : Map_radial::poisson_compact : " << endl
650 << " domain no. " << lz << " : the l=0 part of the effective source is > 1.e-10 : "
651 << somlzero << endl ;
652 }
653 }
654
655 // Resolution of the equation
656 //----------------------------
657
658 bool reamorce = (niter == 0) ;
659
660 assert(sour_j.c_cf != 0x0) ;
661
662 psi.set_etat_zero() ; // to call Cmp::del_deriv().
663 psi.set_etat_qcq() ;
664
665 vpsi = sol_poisson_compact(mpaff, *(sour_j.c_cf), ac, bc, reamorce) ;
666
667 // Test: has the equation been correctly solved ?
668 // ---------------------------------------------
669
670 mpaff.laplacien(psi, 0, lap_zeta) ;
671 mpaff.dsdr(psi, grad_zeta) ;
672
673 oper_psi = ta * lap_zeta.va + tb * grad_zeta.va ;
674 oper_psi.std_base_scal() ;
675 oper_psi.coef() ;
676 oper_psi.ylm() ;
677
678 Mtbl_cf diff_opsou = *oper_psi.c_cf - *sour_j.c_cf ;
679 //cout << " Coef of oper_psi - sour_j : " << endl ;
680 // diff_opsou.affiche_seuil(cout, 4, 1.e-11) ;
681
682 cout << "poisson_compact: step " << niter << " : " << endl ;
683 for (int lz=0; lz<nzet; lz++) {
684 double maxc = fabs( max(*(vpsi.c_cf))(lz) ) ;
685 double minc = fabs( min(*(vpsi.c_cf))(lz) ) ;
686 double max_abs_psi = ( maxc > minc ) ? maxc : minc ;
687
688 maxc = fabs( max(*(sour_j.c_cf))(lz) ) ;
689 minc = fabs( min(*(sour_j.c_cf))(lz) ) ;
690 double max_abs_sou = ( maxc > minc ) ? maxc : minc ;
691
692 maxc = fabs( max(diff_opsou)(lz) ) ;
693 minc = fabs( min(diff_opsou)(lz) ) ;
694 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
695
696 cout << " lz = " << lz << " : max |psi| |sou| |oper(psi)-sou|: "
697 << max_abs_psi << " " << max_abs_sou << " "
698 << max_abs_diff << endl ;
699 }
700
701 // Relaxation :
702 // -----------
703
704 vpsi.ylm_i() ; // Inverse spherical harmonics transform
705
706 psi = relax * psi + unmrelax * psi_jm1 ;
707
708 tdiff = diffrel(psi, psi_jm1) ;
709
710 diff = max(tdiff) ;
711
712 cout << " Relative difference psi^J <-> psi^{J-1} : "
713 << tdiff << endl ;
714
715 // Updates for the next iteration
716 // -------------------------------
717
718 psi_jm1 = psi ;
719 niter++ ;
720
721 } // End of iteration
722 while ( (diff > precis) && (niter < nitermax) ) ;
723
724//==========================================================================
725// End of iteration
726//==========================================================================
727
728 psi.annule(nzet, nz-1) ;
729
730}
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::srstdsdp
const Cmp & srstdsdp() const
Returns of *this .
Definition cmp_deriv.C:130
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:351
Lorene::Cmp::laplacien
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
Definition cmp_deriv.C:245
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition cmp.C:307
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:292
Lorene::Cmp::srdsdt
const Cmp & srdsdt() const
Returns of *this .
Definition cmp_deriv.C:108
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Cmp::dsdr
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:87
Lorene::Map_af
Affine radial mapping.
Definition map.h:2042
Lorene::Map_af::laplacien
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
Definition map_af_lap.C:182
Lorene::Map_af::dsdr
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition map_af_deriv.C:282
Lorene::Map_radial::poisson_compact
virtual void poisson_compact(const Cmp &source, const Cmp &aa, const Tenseur &bb, const Param &par, Cmp &psi) const
Resolution of the elliptic equation in the case where the stellar interior is covered by a single do...
Definition map_radial_poisson_cpt.C:158
Lorene::Map_radial::dxdr
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1575
Lorene::Mtbl_cf
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:196
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::set
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition mtbl.h:221
Lorene::Mtbl::annule_hard
void annule_hard()
Sets the Mtbl to zero in a hard way.
Definition mtbl.C:315
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:295
Lorene::Param::get_double
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:364
Lorene::Param::get_int_mod
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition param.C:433
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::t
double * t
The array of double.
Definition tbl.h:173
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:304
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:707
Lorene::Tenseur::get_mp
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:702
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition tenseur.h:710
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Lorene::Valeur::sx
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ).
Definition valeur_sx.C:113
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::dsdx
const Valeur & dsdx() const
Returns of *this.
Definition valeur_dsdx.C:114
Lorene::Valeur::mult_x
const Valeur & mult_x() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ).
Definition valeur_mult_x.C:108
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:312
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:151
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:134
Lorene::Valeur::d2sdx2
const Valeur & d2sdx2() const
Returns of *this.
Definition valeur_d2sdx2.C:117
Lorene::Valeur::std_base_scal
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:827
Lorene::Valeur::lapang
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition valeur_lapang.C:75
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:461
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:438
Lorene::bvect_spher
Base_vect_spher bvect_spher
Base class for coordinate mappings.
Definition map.h:701
Lorene
Lorene prototypes.
Definition app_hor.h:67