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 char 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
145 namespace Lorene {
146 Mtbl_cf sol_poisson_compact(const Mtbl_cf&, double, double, bool) ;
147 Mtbl_cf sol_poisson_compact(const Map_af&, const Mtbl_cf&, const Tbl&,
148  const Tbl&, bool) ;
149 
150 
152  // Single domain version //
154 
155 void 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 
453 void 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
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
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 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.
Definition: map_af_deriv.C:279
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
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
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_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 Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:699
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
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.
Definition: valeur_dsdx.C:111
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.
Definition: valeur_coef.C:148
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
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.
Definition: valeur_lapang.C:72
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