144 #include "graphique.h" 145 #include "utilitaires.h" 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&,
166 assert(source.
get_etat() != ETATNONDEF) ;
167 assert(aa.
get_etat() != ETATNONDEF) ;
168 assert(bb.
get_etat() != ETATNONDEF) ;
184 double unmrelax = 1. - relax ;
189 if ( source.
get_etat() == ETATZERO ) {
199 double dxdr_c = tmp(0, 0, 0, 0) ;
201 double alph = dxdr_c * dxdr_c * aa(0, 0, 0, 0) ;
203 const Valeur& b_r = bb(0).va ;
207 double bet =
min(vatmp)(0) ;
209 cout <<
"Map_radial::poisson_compact : alph, bet : " << alph <<
" " 223 Cmp b_grad_psi(
this) ;
247 b_grad_psi = bb(0) % psi.
dsdr() +
262 aux_psi = 2.*dpsi + (vpsi.
lapang()).sx() ;
267 lap_xi_psi = d2psi + aux_psi.
sx() ;
272 + alph * ( lap_xi_psi - (lap_xi_psi.
mult_x()).mult_x() )
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 ;
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 ;
300 bool reamorce = (niter == 0) ;
302 assert(sour_j.
c_cf != 0x0) ;
306 vpsi = sol_poisson_compact( *(sour_j.
c_cf), alph, bet, reamorce ) ;
384 aux_psi = vpsi.
dsdx() ;
386 aux_psi = 2.*aux_psi + (vpsi.
lapang()).sx() ;
388 lap_xi_psi = vpsi.
d2sdx2() + aux_psi.
sx() ;
390 oper_psi = alph * ( lap_xi_psi - (lap_xi_psi.
mult_x()).mult_x() )
391 + bet * (vpsi.
dsdx()).mult_x() ;
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 ;
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 ;
403 maxc = fabs(
max(diff_opsou)(0) ) ;
404 minc = fabs(
min(diff_opsou)(0) ) ;
405 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
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 ;
422 psi = relax * psi + unmrelax * psi_jm1 ;
424 tdiff =
diffrel(psi, psi_jm1) ;
429 " Relative difference psi^J <-> psi^{J-1} : " 430 << tdiff(0) << endl ;
441 while ( (diff > precis) && (niter < nitermax) ) ;
469 assert(source.
get_etat() != ETATNONDEF) ;
470 assert(aa.
get_etat() != ETATNONDEF) ;
471 assert(bb.
get_etat() != ETATNONDEF) ;
484 if ( source.
get_etat() == ETATZERO ) {
495 double unmrelax = 1. - relax ;
512 ac.
set(0,0) = ap(0,0,0,0) ;
513 ac.
set(0,2) = ap(0,0,0,nrm1) - ac(0,0) ;
515 bc.
set(0,1) = bp(0,0,0,nrm1) ;
518 for (
int lz=1; lz<nzet-1; lz++) {
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) ) ;
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) ) ;
528 int lext = nzet - 1 ;
530 ac.
set(lext,0) = 0.5 * ap(lext,0,0,0) ;
531 ac.
set(lext,1) = - ac(lext,0) ;
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) ) ;
536 cout <<
"ac : " << ac << endl ;
537 cout <<
"bc : " << bc << endl ;
546 for (
int lz=0; lz<nzet; lz++) {
548 double* tta = ta.
set(lz).
t ;
549 double* ttb = tb.
set(lz).
t ;
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]) ;
593 Cmp b_grad_psi(
this) ;
624 Cmp lap_zeta(mpaff) ;
627 Cmp grad_zeta(mpaff) ;
628 mpaff.
dsdr(psi, grad_zeta) ;
632 + tb * grad_zeta.
va - b_grad_psi.
va ;
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 ;
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 ;
658 bool reamorce = (niter == 0) ;
660 assert(sour_j.
c_cf != 0x0) ;
665 vpsi = sol_poisson_compact(mpaff, *(sour_j.
c_cf), ac, bc, reamorce) ;
671 mpaff.
dsdr(psi, grad_zeta) ;
673 oper_psi = ta * lap_zeta.
va + tb * grad_zeta.
va ;
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 ;
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 ;
692 maxc = fabs(
max(diff_opsou)(lz) ) ;
693 minc = fabs(
min(diff_opsou)(lz) ) ;
694 double max_abs_diff = ( maxc > minc ) ? maxc : minc ;
696 cout <<
" lz = " << lz <<
" : max |psi| |sou| |oper(psi)-sou|: " 697 << max_abs_psi <<
" " << max_abs_sou <<
" " 698 << max_abs_diff << endl ;
706 psi = relax * psi + unmrelax * psi_jm1 ;
708 tdiff =
diffrel(psi, psi_jm1) ;
712 cout <<
" Relative difference psi^J <-> psi^{J-1} : " 722 while ( (diff > precis) && (niter < nitermax) ) ;
const Cmp & dsdr() const
Returns of *this .
const Map * get_mp() const
Returns the mapping.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
const Valeur & dsdx() const
Returns of *this.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Valeur & lapang() const
Returns the angular Laplacian of *this.
void ylm_i()
Inverse of ylm()
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
void coef() const
Computes the coeffcients of *this.
void annule(int l)
Sets the Cmp to zero in a given domain.
const Cmp & srstdsdp() const
Returns of *this .
void ylm()
Computes the coefficients of *this.
int get_etat() const
Returns the logical state.
double & set(int i)
Read/write of a particular element (index i) (1D case)
const Cmp & srdsdt() const
Returns of *this .
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Values and coefficients of a (real-value) function.
double * x
Array of values of at the nr collocation points.
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
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).
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const Map * get_mp() const
Returns pointer on the mapping.
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
double * t
The array of double.
int get_nzone() const
Returns the number of domains.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
int get_etat() const
Returns the logical state.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
void annule_hard()
Sets the Mtbl to zero in a hard way.
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined.
Coefficients storage for the multi-domain spectral method.
Base_vect_spher bvect_spher
Orthonormal vectorial basis associated with the coordinates of the mapping.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
virtual void laplacien(const Scalar &uu, int zec_mult_r, Scalar &lap) const
Computes the Laplacian of a scalar field.
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
void annule_hard()
Sets the Tbl to zero in a hard way.
Valeur va
The numerical value of the Cmp.
virtual void dsdr(const Cmp &ci, Cmp &resu) const
Computes of a Cmp.
Tensor handling *** DEPRECATED : use class Tensor instead ***.