From 9bcad028b18bc6ca55dc57e035221c79fba902a5 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Thu, 14 Jul 2022 00:38:34 +0200 Subject: [PATCH] improve error handling in the gray_main routine MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - rename errocodes → gray_errors - restructure the errors into a `error_spec` type - make the list of errors easily extensible - rewrite the `print_errn`, `print_errhcd` (now `print_err_raytracing`, `print_err_ecrh_cd`) subroutine to handle arbitrary errors - add functions to easily manipulate errors (`raise_error`, `has_error`, `is_critical`) - remove print statements from quadpack - log all errors to stderr using the logger module --- src/dispersion.f90 | 61 +++++++------ src/eccd.f90 | 6 +- src/errcodes.f90 | 83 ------------------ src/gray_core.f90 | 40 ++++----- src/gray_errors.f90 | 208 ++++++++++++++++++++++++++++++++++++++++++++ src/quadpack.f90 | 7 +- 6 files changed, 265 insertions(+), 140 deletions(-) delete mode 100644 src/errcodes.f90 create mode 100644 src/gray_errors.f90 diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 30a1c8e..81ee716 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -225,8 +225,9 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & ! method starting from the cold plasma solution ! ! Reference: https://doi.org/10.13182/FST08-A1660 - use, intrinsic :: ieee_arithmetic, only: ieee_is_finite - + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use gray_errors, only : gray_error, warmdisp_convergence, warmdisp_result, & + raise_error implicit none ! subroutine arguments @@ -246,15 +247,8 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & ! Outputs - ! error code: - ! 0: converged within `max_iters` - ! 1: failed to converge, returned fallback value - ! 2: failed to converge, returned last value - ! 4: converged to invalid value - ! 5: failed to converge, fallback is also invalid - ! 6: failed to converge, last value is invalid - ! Note: invalid means either IEEE 754 NaN or ±Infinity. - integer, intent(out) :: error + ! GRAY error code + integer(kind=gray_error), intent(out) :: error ! orthogonal refractive index: N⊥ complex(wp_), intent(out) :: Npr ! polarization vector: ε (local cartesian frame, z∥B̅) @@ -262,7 +256,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & ! Parameters - ! switch for the dielectic tensor model: + ! switch for the dielectric tensor model: ! 1: weakly relativistic ! 2: fully relativistic (faster variant) ! 3: fully relativistic (slower variant) @@ -315,7 +309,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & Npr2_err = 1 ! Compute the N⊥-independent part of the dielectric tensor - call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) + call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error) iterate_to_solve_disperion: do i = 1, max_iters Npr_prev = sqrt(Npr2_prev) @@ -433,23 +427,23 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & if (i > max_iters) then if (fallback) then ! first try failed, use fallback - error = 1 + error = raise_error(error, warmdisp_convergence, 0) Npr2 = Npr2_fallback else ! first try failed, no retry - error = 2 + error = raise_error(error, warmdisp_convergence, 1) end if end if ! Catch NaN and ±Infinity values if (.not. ieee_is_finite(Npr2%re) .or. .not. ieee_is_finite(Npr2%im)) then Npr2 = 0 - error = error + 4 + error = raise_error(error, warmdisp_result, 0) end if if (Npr2%re < 0 .and. Npr2%im < 0) then Npr2 = sqrt(0.5_wp_) ! Lorenzo, what should we have here?? - error = 99 + error = raise_error(error, warmdisp_result, 1) end if ! Final result @@ -482,11 +476,12 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & end subroutine warmdisp -subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) +subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error) ! Returns the N⊥-independent parts of the dielectric tensor, namely ! ε₃₃⁰ and the 3x3 tensors ε̅^(l) for l=1,nlarmor. ! ! Reference: https://doi.org/10.13182/FST08-A1660 + use gray_errors, only : gray_error implicit none @@ -503,7 +498,7 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) ! Parameters - ! switch for the dielectic tensor model: + ! switch for the dielectric tensor model: ! 1: weakly relativistic ! 2: fully relativistic (faster variant) ! 3: fully relativistic (slower variant) @@ -518,6 +513,8 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) complex(wp_), intent(out) :: e330 ! tensors ε̅_ij^(l), defined in eq. 10 complex(wp_), intent(out) :: epsl(3,3,nlarmor) + ! GRAY error code + integer(kind=gray_error), intent(out) :: error ! local variables integer :: l, n @@ -536,10 +533,10 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) if (model == 1) then ! Compute the (combinations of) Shkarofsky functions F_q(n) ! in which the functions Q⁺_hl(n), Q⁻_hl(n) can be expressed - call fsup(nlarmor, Y, Npl, mu, Fp, Fm) + call fsup(nlarmor, Y, Npl, mu, Fp, Fm, error) else ! Compute the real and imaginary parts of the integrals - ! defining the functions Q_hl(n), Q_hl(n) + ! defining the functions Q_hl(n), Q_hl(-n) block real(wp_) :: cr, ci, cmxw cmxw = 1 + 15/(8*mu) + 105/(128 * mu**2) @@ -552,7 +549,8 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl) call hermitian(rr, Y, mu, Npl, cr, model, nlarmor) else ! slow variant - call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor) + call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor, error) + if (error /= 0) error = 2 end if ! imaginary part call antihermitian(ri, Y, mu, Npl, ci, nlarmor) @@ -868,8 +866,9 @@ end subroutine hermitian ! ! ! -subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm) - use quadpack, only : dqagsmv !dqagimv +subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm,error) + use gray_errors, only : gray_error, dielectric_tensor, raise_error + use quadpack, only : dqagsmv implicit none ! local constants integer,parameter :: lw=5000,liw=lw/4,npar=7 @@ -878,6 +877,7 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm) integer :: lrm,fast real(wp_) :: yg,mu,npl,cr real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr + integer(kind=gray_error), intent(out) :: error ! local variables integer :: n,m,ih,k,n1,nn,llm,neval,ier,last,ihmin integer, dimension(liw) :: iw @@ -919,9 +919,9 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm) if(n.eq.0.and.m.eq.0) ihmin=2 do ih=ihmin,2 apar(7) = real(ih,wp_) -! call dqagimv(fhermit,bound,2,apar,npar,epsa,epsr,resfh, call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, & epp,neval,ier,liw,lw,last,iw,w) + if (ier /= 0) error = raise_error(error, dielectric_tensor, 1) rr(n,ih,m) = resfh end do end do @@ -1216,13 +1216,16 @@ subroutine expinit end subroutine expinit -pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) +pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error) + use gray_errors, only : gray_error, dielectric_tensor, raise_error + implicit none ! subroutine arguments integer, intent(in) :: lrm real(wp_), intent(in) :: yg, npl, mu complex(wp_), intent(out), dimension(0:lrm,0:2) :: cefp, cefm + integer(kind=gray_error), intent(out) :: error ! local constants real(wp_), parameter :: apsicr=0.7_wp_ @@ -1256,7 +1259,9 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) y0=phim end if call zetac (xp,yp,zrp,zip,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) call zetac (xm,ym,zrm,zim,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) ! czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) @@ -1273,6 +1278,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) cphi=phim if(alpha.lt.0) cphi=-im*phim call zetac (x0,y0,zr0,zi0,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) cz0=dcmplx(zr0,zi0) cdz0=2.0_wp_*(one-cphi*cz0) cf32=cdz0 @@ -1319,7 +1325,9 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) y0=phim end if call zetac (xp,yp,zrp,zip,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) call zetac (xm,ym,zrm,zim,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) ! czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) @@ -1336,6 +1344,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm) cphi=phim if(alpha.lt.0) cphi=-im*phim call zetac (x0,y0,zr0,zi0,iflag) + if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0) cz0=dcmplx(zr0,zi0) cdz0=2.0_wp_*(one-cphi*cz0) cf32=cdz0 diff --git a/src/eccd.f90 b/src/eccd.f90 index 930b05c..b578669 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -176,7 +176,7 @@ contains ithn,cst2,fcur,eccdpar,effjcd,iokhawa,ierr) use const_and_precisions, only : pi,qesi=>e_,mesi=>me_, & vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ - use errcodes, only : pcdfp,pcdfc + use gray_errors, only : fpp_integration, fcur_integration, raise_error use quadpack, only : dqagsmv implicit none ! local constants @@ -258,7 +258,7 @@ contains epp,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then - ierr=ibset(ierr,pcdfp) + ierr=raise_error(ierr, fpp_integration) return end if @@ -302,7 +302,7 @@ contains if (abs(resji).lt.1.0e-10_wp_) then resji=0.0_wp_ else - ierr=ibset(ierr,pcdfc+iokhawa+i) + ierr=raise_error(ierr, fcur_integration, iokhawa + i) return end if end if diff --git a/src/errcodes.f90 b/src/errcodes.f90 deleted file mode 100644 index a1a7f7f..0000000 --- a/src/errcodes.f90 +++ /dev/null @@ -1,83 +0,0 @@ -module errcodes - implicit none - integer, parameter :: pnpl = 0, lnpl = 2 ! N// too large (2 thresholds) - integer, parameter :: pconv = pnpl + lnpl, lconv = 2 ! Disp. convergence (disabled/failed) - integer, parameter :: pnprre = pconv + lconv, lnprre = 1 ! Re(Nperp)<0 - integer, parameter :: palph = pnprre+ lnprre, lalph = 1 ! alpha<0 - integer, parameter :: pcdfp = palph + lalph, lcdfp = 1 ! fpp integration - integer, parameter :: pcdfc = pcdfp + lcdfp, lcdfc = 3 ! fcur integration (no trap/trap 1/trap 2) -contains - - subroutine check_err(ierr,istop) - implicit none -! arguments - integer, intent(in) :: ierr - integer, intent(inout) :: istop -! if(ibits(ierr,pnpl, lnpl )>1 .or. & ! N// too large -! ibits(ierr,palph,lalph)==1) then ! alpha < 0 -! istop = 1 -! else -! istop = 0 -! end if - -! if(ibits(ierr,pnpl, lnpl )>1) istop = 1 ! N// too large - if(ibits(ierr,palph,lalph)==1) istop = 1 ! alpha < 0 - - end subroutine check_err - - - subroutine print_errn(ierr,i,anpl) - use const_and_precisions, only : wp_ - implicit none -! arguments - integer, intent(in) :: ierr,i - real(wp_), intent(in) :: anpl -! local variables - integer :: ierrs - - ierrs = ibits(ierr,pnpl,lnpl) - if(ierrs/=0) print*,i,' IERR = ', ierrs*2**pnpl,' N// = ',anpl - end subroutine print_errn - - - subroutine print_errhcd(ierr,i,anprre,anprim,alpha) - use const_and_precisions, only : wp_ - implicit none -! arguments - integer, intent(in) :: ierr,i - real(wp_), intent(in) :: anprre,anprim,alpha -! local variables - integer :: ierrs - - ierrs=ibits(ierr,pconv,lconv) - if(ierrs==1) then - print*,i,' IERR = ', ierrs*2**pconv,' Nwarm = ',anprre,anprim, & - ': convergence disabled.' - else if(ierrs==2) then - print*,i,' IERR = ', ierrs*2**pconv,' Nwarm = ',anprre,anprim, & - ': convergence failed.' - end if - - ierrs=ibits(ierr,pnprre,lnprre) - if(ierrs/=0) & - print*,i,' IERR = ', ierrs*2**pnprre,' Nwarm = ',anprre,anprim, & - ': Re(Nwarm)<0 or Nwarm**2 invalid.' - - ierrs=ibits(ierr,palph,lalph) - if(ierrs/=0) & - print*,i,' IERR = ', ierrs*2**palph,' alpha = ',alpha - - ierrs=ibits(ierr,pcdfp,lcdfp) - if(ierrs/=0) & - print*,i,' IERR = ', ierrs*2**pcdfp,' fpp integration error' - - ierrs=ibits(ierr,pcdfc,lcdfc) - if(ibits(ierrs,0,1)/=0) & - print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (no trapping)' - if(ibits(ierrs,1,1)/=0) & - print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (1st trapping region)' - if(ibits(ierrs,2,1)/=0) & - print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (2nd trapping region)' - end subroutine print_errhcd - -end module errcodes diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 746cfa1..2106d78 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -15,7 +15,7 @@ contains print_parameters use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight, rayi2jk - use errcodes, only : check_err, print_errn, print_errhcd + use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd use magsurf_data, only : flux_average, dealloc_surfvec use beamdata, only : init_btr, dealloc_beam use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & @@ -62,7 +62,7 @@ contains ! i: integration step, jk: global ray index integer :: i, jk - integer :: iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt + integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd,index_rt integer :: ip,ib,iopmin,ipar,iO integer :: igrad_b,istop_pass,nbeam_pass,nlim logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff @@ -271,7 +271,6 @@ contains if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) error = 0 - istop = 0 ! stop flag for current beam iopmin = 10 ! =========== rays loop BEGIN =========== @@ -287,7 +286,7 @@ contains ! update global error code and print message if(ierrn/=0) then error = ior(error,ierrn) - call print_errn(ierrn,i,anpl) + call print_err_raytracing(ierrn, i, anpl) end if ! check entrance/exit plasma/wall @@ -379,7 +378,7 @@ contains ierrn,igrad_b) ! * update derivatives after reflection if(ierrn/=0) then ! * update global error code and print message error = ior(error,ierrn) - call print_errn(ierrn,i,anpl) + call print_err_raytracing(ierrn, i, anpl) end if if(jk == 1 .and. ip == 1) then ! * 1st pass, polarization angles at reflection for central ray @@ -437,11 +436,11 @@ contains alpha, didp, nharm, nhf, iokhawa, ierrhcd) anprre = Npr_warm%re anprim = Npr_warm%im + if(ierrhcd /= 0) then + error = ior(error, ierrhcd) + call print_err_ecrh_cd(ierrhcd, i, Npr_warm, alpha) + end if end block - if(ierrhcd/=0) then - error = ior(error,ierrhcd) - call print_errhcd(ierrhcd,i,anprre,anprim,alpha) - end if else tekev=zero alpha=zero @@ -498,12 +497,10 @@ contains call print_projxyzt(stv,yw,0) end if - ! check for any error code and stop if necessary - call check_err(error,istop) ! test whether further trajectory integration is unnecessary call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling - if(istop == 1) then ! stop propagation for current beam + if(is_critical(error)) then ! stop propagation for current beam istop_pass = istop_pass +1 ! * +1 non propagating beam if(ip < params%raytracing%ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams exit @@ -1029,7 +1026,7 @@ contains bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update - use errcodes, only : pnpl + use gray_errors, only : raise_error, large_npl implicit none ! arguments real(wp_), dimension(3), intent(in) :: xv,anv @@ -1053,13 +1050,12 @@ contains call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, & dery,anpl,anpr,ddr,ddi,dersdst,derdnm) - error=0 - if( abs(anpl) > anplth1) then - if(abs(anpl) > anplth2) then - error=ibset(error,pnpl+1) - else - error=ibset(error,pnpl) - end if + if (abs(anpl) > anplth2) then + error = raise_error(error, large_npl, 1) + else if (abs(anpl) > anplth1) then + error = raise_error(error, large_npl, 0) + else + error = 0 end if end subroutine ywppla_upd @@ -1711,7 +1707,7 @@ contains use equilibrium, only : sgnbphi use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl - use errcodes, only : palph + use gray_errors, only : negative_absorption, raise_error use magsurf_data, only : fluxval implicit none @@ -1814,7 +1810,7 @@ contains end block if (alpha < 0) then - error = ibset(error, palph) + error = raise_error(error, negative_absorption) return end if diff --git a/src/gray_errors.f90 b/src/gray_errors.f90 new file mode 100644 index 0000000..7a02db8 --- /dev/null +++ b/src/gray_errors.f90 @@ -0,0 +1,208 @@ +! This module contains the error codes handled by the gray_main subroutine +! +! An error is actually a 32bit integer bitmask that can encode several +! error cases. The individual errors can be raised using the `raise_error` +! function and combined simply with the intrinsic `ior` function. +module gray_errors + + use logger, only : log_error + use, intrinsic :: iso_fortran_env, only : int32 + + implicit none + + ! The type of a GRAY error is `integer(kind=gray_error)` + integer, parameter :: gray_error = int32 + + ! Specification of a GRAY error + type error_spec + integer :: offset ! bitmask offset + integer :: subcases ! number of error subcases (max: 10) + character(32) :: mod ! module raising the error + character(32) :: proc ! procedure raising the error + character(64) :: msg(0:9) ! error message (one for each subcase) + end type + + ! macros used for defining errors +# define _ , +# define list(x) reshape([character(64) :: x], [10], ['']) +# define after(x) x%offset + x%subcases + + ! All GRAY errors + + type(error_spec), parameter :: large_npl = & + error_spec(offset=0, subcases=2, & + mod='gray_core', proc='gray_main', & + msg=list('N∥ is too large (>0.99)'_ + 'N∥ is too large (>1.05)')) + + type(error_spec), parameter :: dielectric_tensor = & + error_spec(offset=after(large_npl), subcases=2, & + mod='gray_core', proc='gray_main', & + msg=list('ε tensor, overflow in `fsup`'_ + 'ε tensor, integration error in `hermitian_2`')) + + type(error_spec), parameter :: warmdisp_convergence = & + error_spec(offset=after(dielectric_tensor), subcases=2, & + mod='dispersion', proc='warmdisp', & + msg=list('failed to converge, returned fallback value'_ + 'failed to converge, returned last value')) + + type(error_spec), parameter :: warmdisp_result = & + error_spec(offset=after(warmdisp_convergence), subcases=2, & + mod='dispersion', proc='warmdisp', & + msg=list('final N⊥² is NaN or ±Infinity'_ + 'final N⊥² in 3rd quadrant')) + + type(error_spec), parameter :: negative_absorption = & + error_spec(offset=after(warmdisp_result), subcases=1, & + mod='gray_core', proc='alpha_effj', & + msg=list('negative absorption coeff.')) + + type(error_spec), parameter :: fpp_integration = & + error_spec(offset=after(negative_absorption), subcases=1, & + mod='eccd', proc='eccdeff', & + msg=list('fpp integration error')) + + type(error_spec), parameter :: fcur_integration = & + error_spec(offset=after(fpp_integration), subcases=3, & + mod='eccd', proc='eccdeff', & + msg=list('fcur integration error (no trapping)'_ + 'fcur integration error (1st trapping region)'_ + 'fcur integration error (2st trapping region)')) + + ! Errors occuring during raytracing + type(error_spec), parameter :: raytracing_errors(*) = [large_npl] + + ! Errors occuring during absorption and current drive computations + type(error_spec), parameter :: ecrh_cd_errors(*) = & + [dielectric_tensor, & + warmdisp_convergence, & + warmdisp_result, & + negative_absorption, & + fpp_integration, & + fcur_integration] + +contains + + pure function is_critical(error) + ! Checks whether critical errors have occurred + + implicit none + + ! subroutines arguments + integer(kind=gray_error), intent(in) :: error + logical :: is_critical + + is_critical = has_error(error, negative_absorption) + end function is_critical + + + pure function has_error(error, spec) + ! Checks whether the `error` bitmask contains the error given by `spec` + implicit none + + ! function arguments + integer(kind=gray_error), intent(in) :: error + type(error_spec), intent(in) :: spec + logical :: has_error + + has_error = ibits(error, spec%offset, spec%subcases) == 1 + end function has_error + + + pure function raise_error(error, spec, subcase) + ! Raise the bits of error `spec` (with optional `subcase` number) + ! in the `error` bitmask. + implicit none + + ! function arguments + integer(kind=gray_error), intent(in) :: error + type(error_spec), intent(in) :: spec + integer, intent(in), optional :: subcase + integer(kind=gray_error) :: raise_error + + raise_error = ibset(error, spec%offset & + + merge(subcase, 0, present(subcase))) + end function raise_error + + + subroutine print_err_raytracing(error, step, Npl) + ! Pretty prints raytracing errors + ! + ! The error and some context (integration step, N∥) + ! is logged to the stderr using the logger module. + + use const_and_precisions, only : wp_ + implicit none + + ! subroutines arguments + integer, intent(in) :: error, step + real(wp_), intent(in) :: Npl + + ! local variables + integer :: slice ! a slice of the bitmask + character(256) :: line ! formatted log line + type(error_spec) :: spec + integer :: i, j + + ! format specifier of the log line + character(*), parameter :: fmt = & + '(a,": ","error=",g0," N∥=",g0.3," step=",g0)' + + ! iterate on the known errors + do i = 1, size(raytracing_errors) + spec = raytracing_errors(i) + slice = ibits(error, spec%offset, spec%subcases) + if (slice == 0) cycle + + ! iterate on the subcases + do j = 0, spec%subcases - 1 + if (ibits(slice, j, 1) == 0) cycle + write(line, fmt) trim(spec%msg(j)), slice * 2**j, Npl, step + call log_error(line, mod=spec%mod, proc=spec%proc) + end do + end do + end subroutine print_err_raytracing + + + subroutine print_err_ecrh_cd(error, step, Npr, alpha) + ! Pretty prints ECRH & CD errors + ! + ! The error and some context (integration step, N⊥, α) + ! is logged to the stderr using the logger module. + + use const_and_precisions, only : wp_ + implicit none + + ! subroutines arguments + integer, intent(in) :: error, step + complex(wp_), intent(in) :: Npr + real(wp_), intent(in) :: alpha + + ! local variables + integer :: slice ! a slice of the bitmask + character(256) :: line ! formatted log line + type(error_spec) :: spec + integer :: i, j + + ! format specifier of the log line + character(*), parameter :: fmt = & + '(a,": ", "error=",g0, " N⊥=",f0.0,sp,f0.0,"i", " α=",g0, " step=",g0)' + + ! iterate on the known errors + do i = 1, size(ecrh_cd_errors) + spec = ecrh_cd_errors(i) + slice = ibits(error, spec%offset, spec%subcases) + if (slice == 0) cycle + + ! iterate on the subcases + do j = 0, spec%subcases - 1 + if (ibits(slice, j, 1) == 0) cycle + write(line, fmt) trim(spec%msg(j)), slice * 2**j, & + Npr%re, Npr%im, alpha, step + call log_error(line, mod=spec%mod, proc=spec%proc) + end do + end do + end subroutine print_err_ecrh_cd + +end module gray_errors diff --git a/src/quadpack.f90 b/src/quadpack.f90 index 3279453..6c3aab8 100644 --- a/src/quadpack.f90 +++ b/src/quadpack.f90 @@ -195,7 +195,6 @@ contains lvl = 0 end if if(ier==6) lvl = 1 - if(ier/=0) print*,'habnormal return from dqags',ier,lvl end subroutine dqags subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, & @@ -1358,7 +1357,6 @@ contains neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! end if - if(ier/=0) print*,'habnormal return from dqagi' end subroutine dqagi subroutine dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, & @@ -2248,7 +2246,6 @@ contains lvl = 0 end if if(ier.eq.6) lvl = 1 - if(ier.ne.0) print*,'habnormal return from dqaps',ier,lvl end subroutine dqagp subroutine dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result, & @@ -3016,7 +3013,6 @@ contains lvl = 0 end if if(ier==6) lvl = 1 - if(ier/=0) print*,'habnormal return from dqags',ier,lvl end subroutine dqagsmv subroutine dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result,abserr,neval, & @@ -3865,7 +3861,6 @@ contains neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) ! end if - if(ier/=0) print*,'habnormal return from dqagi' end subroutine dqagimv subroutine dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit,result,abserr, & @@ -4538,4 +4533,4 @@ contains ((epmach*0.5e+02_wp_)*resabs,abserr) end subroutine dqk15imv -end module quadpack \ No newline at end of file +end module quadpack