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 $" ;
141 #include "graphique.h"
142 #include "utilitaires.h"
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&,
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) ) ;