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