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 $" ;
141#include "graphique.h"
142#include "utilitaires.h"
146Mtbl_cf sol_poisson_compact(
const Mtbl_cf&,
double,
double,
bool) ;
147Mtbl_cf sol_poisson_compact(
const Map_af&,
const Mtbl_cf&,
const Tbl&,
163 assert(source.
get_etat() != ETATNONDEF) ;
164 assert(aa.
get_etat() != ETATNONDEF) ;
165 assert(bb.
get_etat() != ETATNONDEF) ;
181 double unmrelax = 1. - relax ;
186 if ( source.
get_etat() == ETATZERO ) {
196 double dxdr_c = tmp(0, 0, 0, 0) ;
198 double alph = dxdr_c * dxdr_c * aa(0, 0, 0, 0) ;
200 const Valeur& b_r = bb(0).va ;
204 double bet =
min(vatmp)(0) ;
206 cout <<
"Map_radial::poisson_compact : alph, bet : " << alph <<
" "
220 Cmp b_grad_psi(
this) ;
244 b_grad_psi = bb(0) % psi.
dsdr() +
259 aux_psi = 2.*dpsi + (vpsi.
lapang()).sx() ;
264 lap_xi_psi = d2psi + aux_psi.
sx() ;
269 + alph * ( lap_xi_psi - (lap_xi_psi.
mult_x()).mult_x() )
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 ;
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 ;
297 bool reamorce = (niter == 0) ;
299 assert(sour_j.
c_cf != 0x0) ;
303 vpsi = sol_poisson_compact( *(sour_j.
c_cf), alph, bet, reamorce ) ;
381 aux_psi = vpsi.
dsdx() ;
383 aux_psi = 2.*aux_psi + (vpsi.
lapang()).sx() ;
385 lap_xi_psi = vpsi.
d2sdx2() + aux_psi.
sx() ;
387 oper_psi = alph * ( lap_xi_psi - (lap_xi_psi.
mult_x()).mult_x() )
388 + bet * (vpsi.
dsdx()).mult_x() ;
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 ;
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 ;
400 maxc = fabs(
max(diff_opsou)(0) ) ;
401 minc = fabs(
min(diff_opsou)(0) ) ;
402 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
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 ;
419 psi = relax * psi + unmrelax * psi_jm1 ;
421 tdiff =
diffrel(psi, psi_jm1) ;
426 " Relative difference psi^J <-> psi^{J-1} : "
427 << tdiff(0) << endl ;
438 while ( (diff > precis) && (niter < nitermax) ) ;
466 assert(source.
get_etat() != ETATNONDEF) ;
467 assert(aa.
get_etat() != ETATNONDEF) ;
468 assert(bb.
get_etat() != ETATNONDEF) ;
481 if ( source.
get_etat() == ETATZERO ) {
492 double unmrelax = 1. - relax ;
509 ac.
set(0,0) = ap(0,0,0,0) ;
510 ac.
set(0,2) = ap(0,0,0,nrm1) - ac(0,0) ;
512 bc.
set(0,1) = bp(0,0,0,nrm1) ;
515 for (
int lz=1; lz<nzet-1; lz++) {
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) ) ;
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) ) ;
525 int lext = nzet - 1 ;
527 ac.
set(lext,0) = 0.5 * ap(lext,0,0,0) ;
528 ac.
set(lext,1) = - ac(lext,0) ;
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) ) ;
533 cout <<
"ac : " << ac << endl ;
534 cout <<
"bc : " << bc << endl ;
543 for (
int lz=0; lz<nzet; lz++) {
545 double* tta = ta.
set(lz).
t ;
546 double* ttb = tb.
set(lz).
t ;
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]) ;
590 Cmp b_grad_psi(
this) ;
621 Cmp lap_zeta(mpaff) ;
624 Cmp grad_zeta(mpaff) ;
625 mpaff.
dsdr(psi, grad_zeta) ;
629 + tb * grad_zeta.
va - b_grad_psi.
va ;
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 ;
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 ;
655 bool reamorce = (niter == 0) ;
657 assert(sour_j.
c_cf != 0x0) ;
662 vpsi = sol_poisson_compact(mpaff, *(sour_j.
c_cf), ac, bc, reamorce) ;
668 mpaff.
dsdr(psi, grad_zeta) ;
670 oper_psi = ta * lap_zeta.
va + tb * grad_zeta.
va ;
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 ;
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 ;
689 maxc = fabs(
max(diff_opsou)(lz) ) ;
690 minc = fabs(
min(diff_opsou)(lz) ) ;
691 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
693 cout <<
" lz = " << lz <<
" : max |psi| |sou| |oper(psi)-sou|: "
694 << max_abs_psi <<
" " << max_abs_sou <<
" "
695 << max_abs_diff << endl ;
703 psi = relax * psi + unmrelax * psi_jm1 ;
705 tdiff =
diffrel(psi, psi_jm1) ;
709 cout <<
" Relative difference psi^J <-> psi^{J-1} : "
719 while ( (diff > precis) && (niter < nitermax) ) ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
int get_etat() const
Returns the logical state.
Valeur va
The numerical value of the Cmp
const Cmp & srstdsdp() const
Returns of *this .
void annule(int l)
Sets the Cmp to zero in a given domain.
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
const Cmp & srdsdt() const
Returns of *this .
const Map * get_mp() const
Returns the mapping.
const Cmp & dsdr() const
Returns of *this .
double * x
Array of values of at the nr collocation points.
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
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.
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Coefficients storage for the multi-domain spectral method.
Tbl & set(int l)
Read/write of the Tbl in a given domain.
void annule_hard()
Sets the Mtbl to zero in a hard way.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
void annule_hard()
Sets the Tbl to zero in a hard way.
double & set(int i)
Read/write of a particular element (index i) (1D case)
double * t
The array of double.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
const Map * get_mp() const
Returns pointer on the mapping.
int get_etat() const
Returns the logical state.
Values and coefficients of a (real-value) function.
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
void ylm()
Computes the coefficients of *this.
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.
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.
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).
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.