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
148 namespace Lorene {
149 Mtbl_cf sol_poisson_compact(const Mtbl_cf&, double, double, bool) ;
150 Mtbl_cf sol_poisson_compact(const Map_af&, const Mtbl_cf&, const Tbl&,
151  const Tbl&, bool) ;
152 
153 
155  // Single domain version //
157 
158 void 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 
456 void 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 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:87
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:517
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
const Valeur & dsdx() const
Returns of *this.
Definition: valeur_dsdx.C:114
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Valeur & lapang() const
Returns the angular Laplacian of *this.
Definition: valeur_lapang.C:75
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
Multi-domain array.
Definition: mtbl.h:118
const Cmp & srstdsdp() const
Returns of *this .
Definition: cmp_deriv.C:130
Lorene prototypes.
Definition: app_hor.h:67
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
const Cmp & srdsdt() const
Returns of *this .
Definition: cmp_deriv.C:108
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx.C:113
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:215
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:827
const Valeur & d2sdx2() const
Returns of *this.
const Valeur & mult_x() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:292
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
const Map * get_mp() const
Returns pointer on the mapping.
Definition: tenseur.h:702
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
Definition: map.h:1581
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
double * t
The array of double.
Definition: tbl.h:176
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition: param.C:433
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:307
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...
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition: mtbl.h:221
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
void annule_hard()
Sets the Mtbl to zero in a hard way.
Definition: mtbl.C:315
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Definition: map.h:694
Affine radial mapping.
Definition: map.h:2048
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:707
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
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
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
Definition: cmp_deriv.C:245
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Definition: map_af_deriv.C:282
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304