From 78f4c8d470f28558910c41dce06f7b2946f0a8f4 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 15 Jul 2026 14:25:27 +0200 Subject: [PATCH 01/11] Continue markers after numerical recovery failures --- DOC/config.md | 12 +- src/classification.f90 | 1112 ++++++++++++++------------- src/diag_counters.f90 | 20 +- src/simple_main.f90 | 291 ++++--- src/spectre_sympl_orbit.f90 | 77 +- test/tests/test_sympl_testfield.f90 | 453 +++++------ 6 files changed, 1069 insertions(+), 896 deletions(-) diff --git a/DOC/config.md b/DOC/config.md index a4c407b3..e3734da1 100644 --- a/DOC/config.md +++ b/DOC/config.md @@ -15,11 +15,13 @@ reaches its iteration limit, SIMPLE continues only if the final Newton correction is finite and no more than ten times the requested relative tolerance. Each occurrence is reported by the corresponding `*_maxit` - diagnostic. Set it to `.false.` to stop the affected orbit at its last - accepted position. Larger corrections, exterior-domain states, singular - linear systems, non-finite values, and unresolved boundary events always - stop the orbit. Numerical stops are recorded in `orbit_exit_code` and as - `NaN` in `times_lost`; they are unresolved markers, not physical losses. + diagnostic. If bounded retries still cannot resolve a numerical microstep, + SIMPLE retains the last accepted state for that interval, records + `warning_step_skip`, and resumes the same marker at the next microstep. + Numerical failures therefore never become physical losses or terminate a + marker in the default mode. Set the option to `.false.` for strict diagnostic + runs that roll back and stop the affected marker; those stops use a 101--105 + `orbit_exit_code` and `NaN` in `times_lost` and remain distinct from losses. * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi` control the Meiss or Albert canonical-map grid. Their defaults are 62, 63, diff --git a/src/classification.f90 b/src/classification.f90 index f23ae63a..63c2e119 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -1,610 +1,620 @@ module classification -use, intrinsic :: iso_fortran_env, only: dp => real64 -use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan -use omp_lib -use params, only: zstart, zend, times_lost, trap_par, perp_inv, iclass, & - ntimstep, confpart_trap, confpart_pass, notrace_passing, contr_pp, & - class_plot, ntcut, nturns, fast_class, n_tip_vars, nplagr, nder, npl_half, & - nfp, fper, zerolam, num_surf, bmax, bmin, dtaumin, v0, cut_in_per, & - integmode, relerr, ntau, should_skip, orbit_exit_code, unresolved_orbits, & - ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, ORBIT_EXIT_SKIPPED, & - ORBIT_EXIT_NUMERICAL_DOMAIN, ORBIT_EXIT_NUMERICAL_MAXITER, & - ORBIT_EXIT_NUMERICAL_LINEAR, ORBIT_EXIT_NUMERICAL_EVENT -use util, only: twopi, sqrt2 -use velo_mod, only : isw_field_type -use orbit_symplectic, only : orbit_timestep_sympl, get_val -use orbit_symplectic_base, only: SYMPLECTIC_STEP_BOUNDARY, & - SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER, & - SYMPLECTIC_STEP_LINEAR_SOLVE -use simple, only : init_sympl, tracer_t -use cut_detector, only : fract_dimension -use diag_mod, only : icounter -use get_can_sub, only : vmec_to_can, can_to_vmec -use boozer_sub, only : vmec_to_boozer, boozer_to_vmec -use magfie_sub, only : CANFLUX, BOOZER -use check_orbit_type_sub, only : check_orbit_type - - implicit none - - ! Classification result type - separates data from I/O - ! Note: 0=unclassified means the classification was not computed - ! This depends on orbit type (trapped/passing) and class_plot flag - type :: classification_result_t - logical :: passing ! Trapped (false) or passing (true) - logical :: lost ! Orbit lost (true) or confined (false) - integer :: fractal ! Fractal: 0=unclassified, 1=regular, 2=chaotic - integer :: jpar ! J_parallel: 0=unclassified, 1=regular, 2=stochastic - integer :: topology ! Topology: 0=unclassified, 1=ideal, 2=non-ideal - integer :: exit_code ! Physical or numerical orbit termination - end type classification_result_t - - ! output files: - ! iaaa_bou - trapped-passing boundary - ! iaaa_pnt - forced regular passing - ! iaaa_prp - lossed passing - ! iaaa_prt - lossed trapped - ! iaaa_rep - regular passing - ! iaaa_ret - regular trapped - ! iaaa_stp - stochastic passing - ! iaaa_stt - stochastic trapped -integer, parameter :: iaaa_bou=20000, iaaa_pnt=10000, iaaa_prp=10001, iaaa_prt=10002, & - iaaa_rep=10011, iaaa_ret=10012, iaaa_stp=10021, iaaa_stt=10022 - - ! output files: - ! iaaa_jre - regular trapped by J_parallel - ! iaaa_jst - stochastic trapped by J_parallel - ! iaaa_jer - non-classified trapped by J_parallel - ! iaaa_ire - ideal trapped by recurrences and monotonicity - ! iaaa_ist - non-ideal trapped by recurrences and monotonicity - ! iaaa_ier - non-classified trapped by recurrences and monotonicity -integer, parameter :: iaaa_jre=40012, iaaa_jst=40022, iaaa_jer=40032, & - iaaa_ire=50012, iaaa_ist=50022, iaaa_ier=50032 + use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use omp_lib + use params, only: zstart, zend, times_lost, trap_par, perp_inv, iclass, & + ntimstep, confpart_trap, confpart_pass, notrace_passing, contr_pp, & + class_plot, ntcut, nturns, fast_class, n_tip_vars, nplagr, nder, npl_half, & + nfp, fper, zerolam, num_surf, bmax, bmin, dtaumin, v0, cut_in_per, & + integmode, relerr, ntau, should_skip, orbit_exit_code, unresolved_orbits, & + ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, ORBIT_EXIT_SKIPPED, & + ORBIT_EXIT_NUMERICAL_DOMAIN, ORBIT_EXIT_NUMERICAL_MAXITER, & + ORBIT_EXIT_NUMERICAL_LINEAR, ORBIT_EXIT_NUMERICAL_EVENT + use util, only: twopi, sqrt2 + use velo_mod, only : isw_field_type + use orbit_symplectic, only : advance_symplectic_with_retry, & + orbit_timestep_sympl, get_val + use orbit_symplectic_base, only: SYMPLECTIC_STEP_BOUNDARY, & + SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER, & + SYMPLECTIC_STEP_LINEAR_SOLVE, symplectic_newton_warning_mode + use simple, only : init_sympl, tracer_t + use cut_detector, only : fract_dimension + use diag_mod, only : icounter + use get_can_sub, only : vmec_to_can, can_to_vmec + use boozer_sub, only : vmec_to_boozer, boozer_to_vmec + use magfie_sub, only : CANFLUX, BOOZER + use check_orbit_type_sub, only : check_orbit_type + use diag_counters, only: count_event, EVT_WARNING_STEP_SKIP + + implicit none + + ! Classification result type - separates data from I/O + ! Note: 0=unclassified means the classification was not computed + ! This depends on orbit type (trapped/passing) and class_plot flag + type :: classification_result_t + logical :: passing ! Trapped (false) or passing (true) + logical :: lost ! Orbit lost (true) or confined (false) + integer :: fractal ! Fractal: 0=unclassified, 1=regular, 2=chaotic + integer :: jpar ! J_parallel: 0=unclassified, 1=regular, 2=stochastic + integer :: topology ! Topology: 0=unclassified, 1=ideal, 2=non-ideal + integer :: exit_code ! Physical or numerical orbit termination + end type classification_result_t + + ! output files: + ! iaaa_bou - trapped-passing boundary + ! iaaa_pnt - forced regular passing + ! iaaa_prp - lossed passing + ! iaaa_prt - lossed trapped + ! iaaa_rep - regular passing + ! iaaa_ret - regular trapped + ! iaaa_stp - stochastic passing + ! iaaa_stt - stochastic trapped + integer, parameter :: iaaa_bou=20000, iaaa_pnt=10000, iaaa_prp=10001, iaaa_prt=10002, & + iaaa_rep=10011, iaaa_ret=10012, iaaa_stp=10021, iaaa_stt=10022 + + ! output files: + ! iaaa_jre - regular trapped by J_parallel + ! iaaa_jst - stochastic trapped by J_parallel + ! iaaa_jer - non-classified trapped by J_parallel + ! iaaa_ire - ideal trapped by recurrences and monotonicity + ! iaaa_ist - non-ideal trapped by recurrences and monotonicity + ! iaaa_ier - non-classified trapped by recurrences and monotonicity + integer, parameter :: iaaa_jre=40012, iaaa_jst=40022, iaaa_jer=40032, & + iaaa_ire=50012, iaaa_ist=50022, iaaa_ier=50032 contains -subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) - use find_bminmax_sub, only : get_bminmax - use magfie_sub, only : magfie - use plag_coeff_sub, only : plag_coeff - use alpha_lifetime_sub, only : orbit_timestep_axis - - type(tracer_t), intent(inout) :: anorb - integer, intent(in) :: ipart - type(classification_result_t), intent(out) :: class_result - integer :: ierr - real(dp), dimension(5) :: z - real(dp) :: bmod,sqrtg - real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl - integer :: it, ktau, it_f - integer(8) :: kt - logical :: passing - - integer :: ifp_tip,ifp_per - integer, dimension(:), allocatable :: ipoi - real(dp), dimension(:), allocatable :: xp - real(dp), dimension(:,:), allocatable :: coef,orb_sten - real(dp), dimension(:,:), allocatable :: zpoipl_tip,zpoipl_per,dummy2d - real(dp), dimension(n_tip_vars) :: var_tip - real(dp) :: phiper, alam_prev, par_inv - integer :: iper, itip, kper, nfp_tip, nfp_per - - real(dp) :: fraction - real(dp) :: r,theta_vmec,varphi_vmec - logical :: regular - - ! Variables and settings for classification by J_parallel and ideal orbit condition: - integer, parameter :: nfp_dim=3 - integer :: nfp_cot,ideal,ijpar,ierr_cot,iangvar - real(dp), dimension(nfp_dim) :: fpr_in - - zend(:,ipart) = 0d0 - ! - iangvar=2 - ! End variables and settings for classification by J_parallel and ideal orbit condition - ! - - ! Initialize classification result - all unclassified - class_result%passing = .false. - class_result%lost = .false. - class_result%fractal = 0 - class_result%jpar = 0 - class_result%topology = 0 - class_result%exit_code = ORBIT_EXIT_COMPLETED - - ! open(unit=10000+ipart, iostat=stat, status='old') - ! if (stat == 0) close(10000+ipart, status='delete') - ! open(unit=20000+ipart, iostat=stat, status='old') - ! if (stat == 0) close(20000+ipart, status='delete') - - ! Write out trapped-passing boundary at the classification cut: - if(class_plot) then - if(ipart.eq.1) then - z(1)=zstart(1,ipart) + subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) + use find_bminmax_sub, only : get_bminmax + use magfie_sub, only : magfie + use plag_coeff_sub, only : plag_coeff + use alpha_lifetime_sub, only : orbit_timestep_axis + + type(tracer_t), intent(inout) :: anorb + integer, intent(in) :: ipart + type(classification_result_t), intent(out) :: class_result + integer :: ierr + real(dp), dimension(5) :: z + real(dp) :: bmod,sqrtg + real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl + integer :: it, ktau, it_f + integer(8) :: kt + logical :: passing + + integer :: ifp_tip,ifp_per + integer, dimension(:), allocatable :: ipoi + real(dp), dimension(:), allocatable :: xp + real(dp), dimension(:,:), allocatable :: coef,orb_sten + real(dp), dimension(:,:), allocatable :: zpoipl_tip,zpoipl_per,dummy2d + real(dp), dimension(n_tip_vars) :: var_tip + real(dp) :: phiper, alam_prev, par_inv + integer :: iper, itip, kper, nfp_tip, nfp_per + + real(dp) :: fraction + real(dp) :: r,theta_vmec,varphi_vmec + logical :: regular + + ! Variables and settings for classification by J_parallel and ideal orbit condition: + integer, parameter :: nfp_dim=3 + integer :: nfp_cot,ideal,ijpar,ierr_cot,iangvar + real(dp), dimension(nfp_dim) :: fpr_in + + zend(:,ipart) = 0d0 + ! + iangvar=2 + ! End variables and settings for classification by J_parallel and ideal orbit condition + ! + + ! Initialize classification result - all unclassified + class_result%passing = .false. + class_result%lost = .false. + class_result%fractal = 0 + class_result%jpar = 0 + class_result%topology = 0 + class_result%exit_code = ORBIT_EXIT_COMPLETED + + ! open(unit=10000+ipart, iostat=stat, status='old') + ! if (stat == 0) close(10000+ipart, status='delete') + ! open(unit=20000+ipart, iostat=stat, status='old') + ! if (stat == 0) close(20000+ipart, status='delete') + + ! Write out trapped-passing boundary at the classification cut: + if(class_plot) then + if(ipart.eq.1) then + z(1)=zstart(1,ipart) + z(3)=cut_in_per*fper + do kt=0,1000 + z(2)=1d-3*twopi*dble(kt) + call magfie(z(1:3),bmod,sqrtg,bder,hcovar,hctrvr,hcurl) + write(iaaa_bou,*) z(2),sqrt(1.d0-bmod/bmax) + enddo + endif + endif + ! End write out trapped-passing boundary at the classification cut + ! + z = zstart(:, ipart) + r=z(1) + theta_vmec=z(2) + varphi_vmec=z(3) + ! + if(isw_field_type .eq. CANFLUX) then + call vmec_to_can(r,theta_vmec,varphi_vmec,z(2),z(3)) + elseif(isw_field_type .eq. BOOZER) then + call vmec_to_boozer(r,theta_vmec,varphi_vmec,z(2),z(3)) + endif + + ! In case of classification plot all starting points are moved to the classification cut: + if(class_plot) then z(3)=cut_in_per*fper - do kt=0,1000 - z(2)=1d-3*twopi*dble(kt) - call magfie(z(1:3),bmod,sqrtg,bder,hcovar,hctrvr,hcurl) - write(iaaa_bou,*) z(2),sqrt(1.d0-bmod/bmax) - enddo + zstart(2,ipart)=modulo(zstart(2,ipart),twopi) endif - endif - ! End write out trapped-passing boundary at the classification cut - ! - z = zstart(:, ipart) - r=z(1) - theta_vmec=z(2) - varphi_vmec=z(3) - ! - if(isw_field_type .eq. CANFLUX) then - call vmec_to_can(r,theta_vmec,varphi_vmec,z(2),z(3)) - elseif(isw_field_type .eq. BOOZER) then - call vmec_to_boozer(r,theta_vmec,varphi_vmec,z(2),z(3)) - endif - - ! In case of classification plot all starting points are moved to the classification cut: - if(class_plot) then - z(3)=cut_in_per*fper - zstart(2,ipart)=modulo(zstart(2,ipart),twopi) - endif - ! End moving starting points to the classification cut - - if (integmode>0) call init_sympl(anorb%si, anorb%f, z, dtaumin, dtaumin, relerr, integmode) - - call magfie(z(1:3),bmod,sqrtg,bder,hcovar,hctrvr,hcurl) - - !$omp critical - if(num_surf > 1) then - call get_bminmax(z(1),bmin,bmax) - endif - passing = z(5)**2.gt.1.d0-bmod/bmax - trap_par(ipart) = ((1.d0-z(5)**2)*bmax/bmod-1.d0)*bmin/(bmax-bmin) - perp_inv(ipart) = z(4)**2*(1.d0-z(5)**2)/bmod - iclass(:,ipart) = 0 - !$omp end critical - - ! Store passing status in result - class_result%passing = passing - - ! Forced classification of passing as regular: - if(passing .and. should_skip(ipart)) then + ! End moving starting points to the classification cut + + if (integmode>0) call init_sympl(anorb%si, anorb%f, z, dtaumin, dtaumin, relerr, integmode) + + call magfie(z(1:3),bmod,sqrtg,bder,hcovar,hctrvr,hcurl) + !$omp critical - confpart_pass=confpart_pass+1.d0 - zend(:, ipart) = zstart(:, ipart) - times_lost(ipart) = -1.d0 + if(num_surf > 1) then + call get_bminmax(z(1),bmin,bmax) + endif + passing = z(5)**2.gt.1.d0-bmod/bmax + trap_par(ipart) = ((1.d0-z(5)**2)*bmax/bmod-1.d0)*bmin/(bmax-bmin) + perp_inv(ipart) = z(4)**2*(1.d0-z(5)**2)/bmod + iclass(:,ipart) = 0 !$omp end critical - if(class_plot) then + + ! Store passing status in result + class_result%passing = passing + + ! Forced classification of passing as regular: + if(passing .and. should_skip(ipart)) then !$omp critical - write (iaaa_pnt,*) zstart(2,ipart),zstart(5,ipart),trap_par(ipart) + confpart_pass=confpart_pass+1.d0 + zend(:, ipart) = zstart(:, ipart) + times_lost(ipart) = -1.d0 !$omp end critical - endif - iclass(:,ipart) = 1 - ! Mark as regular passing (fractal=1) and not lost - class_result%fractal = 1 - class_result%lost = .false. - class_result%exit_code = ORBIT_EXIT_SKIPPED - orbit_exit_code(ipart) = ORBIT_EXIT_SKIPPED - return - endif - ! End forced classification of passing as regular - - !$omp critical - if (.not. allocated(ipoi)) & - allocate(ipoi(nplagr),coef(0:nder,nplagr),orb_sten(6,nplagr),xp(nplagr)) - !$omp end critical - do it=1,nplagr - ipoi(it)=it - enddo - - nfp_tip=nfp !<= initial array dimension for tips - nfp_per=nfp !<= initial array dimension for periods - !$omp critical - if (.not. allocated(zpoipl_tip)) & - allocate(zpoipl_tip(2,nfp_tip),zpoipl_per(2,nfp_per)) - !$omp end critical - - ! open(unit=10000+ipart, recl=1024, position='append') - ! open(unit=20000+ipart, recl=1024, position='append') - - ifp_tip=0 !<= initialize footprint counter on tips - ifp_per=0 !<= initialize footprint counter on periods - - icounter=0 - phiper=0.0d0 - - - kt = 0 - if (passing) then - !$omp atomic - confpart_pass(1)=confpart_pass(1)+1.d0 - else - !$omp atomic - confpart_trap(1)=confpart_trap(1)+1.d0 - end if - - !-------------------------------- - ! Initialize tip detector - - itip=npl_half+1 - alam_prev=z(5) - - ! End initialize tip detector - !-------------------------------- - ! Initialize period crossing detector - - iper=npl_half+1 - kper=int(z(3)/fper) - - ! End initialize period crossing detector - !-------------------------------- - ! - ! Initialize classification by J_parallel and ideal orbit condition: - nfp_cot=0 - ! End Initialize classification by J_parallel and ideal orbit condition - ! - par_inv = 0d0 - regular = .False. - do it=2,ntimstep - if (regular) then ! regular orbit, will not be lost - if(passing) then - !$omp atomic - confpart_pass(it)=confpart_pass(it)+1.d0 - else - !$omp atomic - confpart_trap(it)=confpart_trap(it)+1.d0 + if(class_plot) then + !$omp critical + write (iaaa_pnt,*) zstart(2,ipart),zstart(5,ipart),trap_par(ipart) + !$omp end critical endif - kt = kt+ntau - cycle + iclass(:,ipart) = 1 + ! Mark as regular passing (fractal=1) and not lost + class_result%fractal = 1 + class_result%lost = .false. + class_result%exit_code = ORBIT_EXIT_SKIPPED + orbit_exit_code(ipart) = ORBIT_EXIT_SKIPPED + return endif - do ktau=1,ntau - if (integmode <= 0) then - call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr) - else - call orbit_timestep_sympl(anorb%si, anorb%f, ierr) - z(1:3) = anorb%si%z(1:3) - z(4) = dsqrt(anorb%f%mu*anorb%f%Bmod+0.5d0*anorb%f%vpar**2) - z(5) = anorb%f%vpar/(z(4)*sqrt2) - endif + ! End forced classification of passing as regular - if(ierr.ne.0) then - call classify_classifier_exit(ierr, integmode, & - class_result%lost, class_result%exit_code) - if(class_plot .and. class_result%lost) then - call output_lost_orbit_starting_data(ipart, passing) + !$omp critical + if (.not. allocated(ipoi)) & + allocate(ipoi(nplagr),coef(0:nder,nplagr),orb_sten(6,nplagr),xp(nplagr)) + !$omp end critical + do it=1,nplagr + ipoi(it)=it + enddo + + nfp_tip=nfp !<= initial array dimension for tips + nfp_per=nfp !<= initial array dimension for periods + !$omp critical + if (.not. allocated(zpoipl_tip)) & + allocate(zpoipl_tip(2,nfp_tip),zpoipl_per(2,nfp_per)) + !$omp end critical + + ! open(unit=10000+ipart, recl=1024, position='append') + ! open(unit=20000+ipart, recl=1024, position='append') + + ifp_tip=0 !<= initialize footprint counter on tips + ifp_per=0 !<= initialize footprint counter on periods + + icounter=0 + phiper=0.0d0 + + + kt = 0 + if (passing) then + !$omp atomic + confpart_pass(1)=confpart_pass(1)+1.d0 + else + !$omp atomic + confpart_trap(1)=confpart_trap(1)+1.d0 + end if + + !-------------------------------- + ! Initialize tip detector + + itip=npl_half+1 + alam_prev=z(5) + + ! End initialize tip detector + !-------------------------------- + ! Initialize period crossing detector + + iper=npl_half+1 + kper=int(z(3)/fper) + + ! End initialize period crossing detector + !-------------------------------- + ! + ! Initialize classification by J_parallel and ideal orbit condition: + nfp_cot=0 + ! End Initialize classification by J_parallel and ideal orbit condition + ! + par_inv = 0d0 + regular = .False. + do it=2,ntimstep + if (regular) then ! regular orbit, will not be lost + if(passing) then + !$omp atomic + confpart_pass(it)=confpart_pass(it)+1.d0 + else + !$omp atomic + confpart_trap(it)=confpart_trap(it)+1.d0 endif - exit - endif - kt = kt+1 - - par_inv = par_inv+z(5)**2*dtaumin ! parallel adiabatic invariant - if(kt.le.nplagr) then !<=first nplagr points to initialize stencil - orb_sten(1:5,kt)=z - orb_sten(6,kt)=par_inv - else !<=normal case, shift stencil - orb_sten(1:5,ipoi(1))=z - orb_sten(6,ipoi(1))=par_inv - ipoi=cshift(ipoi,1) + kt = kt+ntau + cycle endif + do ktau=1,ntau + if (integmode <= 0) then + call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr) + else + call advance_symplectic_with_retry(anorb%si, anorb%f, & + orbit_timestep_sympl, ierr) + if (ierr /= 0 .and. ierr /= SYMPLECTIC_STEP_BOUNDARY .and. & + symplectic_newton_warning_mode) then + call count_event(EVT_WARNING_STEP_SKIP) + ierr = 0 + else if (ierr == 0) then + z(1:3) = anorb%si%z(1:3) + z(4) = dsqrt(anorb%f%mu*anorb%f%Bmod + & + 0.5d0*anorb%f%vpar**2) + z(5) = anorb%f%vpar/(z(4)*sqrt2) + end if + endif - ! Tip detection and interpolation - if(alam_prev.lt.0.d0.and.z(5).gt.0.d0) itip=0 !<=tip has been passed - itip=itip+1 - alam_prev=z(5) - if(kt.gt.nplagr) then !<=use only initialized stencil - if(itip.eq.npl_half) then !<=stencil around tip is complete, interpolate - xp=orb_sten(5,ipoi) - - call plag_coeff(nplagr,nder,zerolam,xp,coef) - - var_tip=matmul(orb_sten(:,ipoi),coef(0,:)) - var_tip(2)=modulo(var_tip(2),twopi) - var_tip(3)=modulo(var_tip(3),twopi) - - ! write(10000+ipart,*) var_tip - - ifp_tip=ifp_tip+1 - if(ifp_tip.gt.nfp_tip) then !<=increase the buffer for banana tips - !$omp critical - allocate(dummy2d(2,ifp_tip-1)) - !$omp end critical - dummy2d=zpoipl_tip(:,1:ifp_tip-1) - !$omp critical - deallocate(zpoipl_tip) - !$omp end critical - nfp_tip=nfp_tip+nfp - !$omp critical - allocate(zpoipl_tip(2,nfp_tip)) - !$omp end critical - zpoipl_tip(:,1:ifp_tip-1)=dummy2d - !$omp critical - deallocate(dummy2d) - !$omp end critical + if(ierr.ne.0) then + call classify_classifier_exit(ierr, integmode, & + class_result%lost, class_result%exit_code) + if(class_plot .and. class_result%lost) then + call output_lost_orbit_starting_data(ipart, passing) endif - zpoipl_tip(:,ifp_tip)=var_tip(1:2) - par_inv = par_inv - var_tip(6) - ! - ! Classification by J_parallel and ideal orbit conditions: - fpr_in(1)=var_tip(1) - fpr_in(2)=var_tip(iangvar) - fpr_in(3)=var_tip(6) - ! - call check_orbit_type(nturns,nfp_cot,fpr_in,ideal,ijpar,ierr_cot) - ! - iclass(1,ipart) = ijpar - iclass(2,ipart) = ideal - ! Store in classification result - class_result%jpar = ijpar - class_result%topology = ideal - ! fast_class without class_plot: stop tracing regular - ! orbits early (they are confined). With class_plot, - ! continue to ntcut so classification output is written. - if(fast_class .and. .not. class_plot & - .and. ierr_cot /= 0 .and. ijpar == 1) then - regular = .True. + exit + endif + kt = kt+1 + + par_inv = par_inv+z(5)**2*dtaumin ! parallel adiabatic invariant + if(kt.le.nplagr) then !<=first nplagr points to initialize stencil + orb_sten(1:5,kt)=z + orb_sten(6,kt)=par_inv + else !<=normal case, shift stencil + orb_sten(1:5,ipoi(1))=z + orb_sten(6,ipoi(1))=par_inv + ipoi=cshift(ipoi,1) + endif + + ! Tip detection and interpolation + if(alam_prev.lt.0.d0.and.z(5).gt.0.d0) itip=0 !<=tip has been passed + itip=itip+1 + alam_prev=z(5) + if(kt.gt.nplagr) then !<=use only initialized stencil + if(itip.eq.npl_half) then !<=stencil around tip is complete, interpolate + xp=orb_sten(5,ipoi) + + call plag_coeff(nplagr,nder,zerolam,xp,coef) + + var_tip=matmul(orb_sten(:,ipoi),coef(0,:)) + var_tip(2)=modulo(var_tip(2),twopi) + var_tip(3)=modulo(var_tip(3),twopi) + + ! write(10000+ipart,*) var_tip + + ifp_tip=ifp_tip+1 + if(ifp_tip.gt.nfp_tip) then !<=increase the buffer for banana tips + !$omp critical + allocate(dummy2d(2,ifp_tip-1)) + !$omp end critical + dummy2d=zpoipl_tip(:,1:ifp_tip-1) + !$omp critical + deallocate(zpoipl_tip) + !$omp end critical + nfp_tip=nfp_tip+nfp + !$omp critical + allocate(zpoipl_tip(2,nfp_tip)) + !$omp end critical + zpoipl_tip(:,1:ifp_tip-1)=dummy2d + !$omp critical + deallocate(dummy2d) + !$omp end critical + endif + zpoipl_tip(:,ifp_tip)=var_tip(1:2) + par_inv = par_inv - var_tip(6) + ! + ! Classification by J_parallel and ideal orbit conditions: + fpr_in(1)=var_tip(1) + fpr_in(2)=var_tip(iangvar) + fpr_in(3)=var_tip(6) + ! + call check_orbit_type(nturns,nfp_cot,fpr_in,ideal,ijpar,ierr_cot) + ! + iclass(1,ipart) = ijpar + iclass(2,ipart) = ideal + ! Store in classification result + class_result%jpar = ijpar + class_result%topology = ideal + ! fast_class without class_plot: stop tracing regular + ! orbits early (they are confined). With class_plot, + ! continue to ntcut so classification output is written. + if(fast_class .and. .not. class_plot & + .and. ierr_cot /= 0 .and. ijpar == 1) then + regular = .True. + endif + ! + ! End classification by J_parallel and ideal orbit conditions endif - ! - ! End classification by J_parallel and ideal orbit conditions endif - endif - ! End tip detection and interpolation - - ! Periodic boundary footprint detection and interpolation - if(z(3).gt.dble(kper+1)*fper) then - iper=0 !<=periodic boundary has been passed - phiper=dble(kper+1)*fper - kper=kper+1 - elseif(z(3).lt.dble(kper)*fper) then - iper=0 !<=periodic boundary has been passed - phiper=dble(kper)*fper - kper=kper-1 - endif - iper=iper+1 - if(kt.gt.nplagr) then !<=use only initialized stencil - if(iper.eq.npl_half) then !<=stencil around periodic boundary is complete, interpolate - xp=orb_sten(3,ipoi)-phiper - - call plag_coeff(nplagr,nder,zerolam,xp,coef) - - var_tip=matmul(orb_sten(:,ipoi),coef(0,:)) - var_tip(2)=modulo(var_tip(2),twopi) - var_tip(3)=modulo(var_tip(3),twopi) - ! write(20000+ipart,*) var_tip - ifp_per=ifp_per+1 - if(ifp_per.gt.nfp_per) then !<=increase the buffer for periodic boundary footprints - !$omp critical - allocate(dummy2d(2,ifp_per-1)) - !$omp end critical - dummy2d=zpoipl_per(:,1:ifp_per-1) - !$omp critical - deallocate(zpoipl_per) - !$omp end critical - nfp_per=nfp_per+nfp - !$omp critical - allocate(zpoipl_per(2,nfp_per)) - !$omp end critical - zpoipl_per(:,1:ifp_per-1)=dummy2d - !$omp critical - deallocate(dummy2d) - !$omp end critical + ! End tip detection and interpolation + + ! Periodic boundary footprint detection and interpolation + if(z(3).gt.dble(kper+1)*fper) then + iper=0 !<=periodic boundary has been passed + phiper=dble(kper+1)*fper + kper=kper+1 + elseif(z(3).lt.dble(kper)*fper) then + iper=0 !<=periodic boundary has been passed + phiper=dble(kper)*fper + kper=kper-1 + endif + iper=iper+1 + if(kt.gt.nplagr) then !<=use only initialized stencil + if(iper.eq.npl_half) then !<=stencil around periodic boundary is complete, interpolate + xp=orb_sten(3,ipoi)-phiper + + call plag_coeff(nplagr,nder,zerolam,xp,coef) + + var_tip=matmul(orb_sten(:,ipoi),coef(0,:)) + var_tip(2)=modulo(var_tip(2),twopi) + var_tip(3)=modulo(var_tip(3),twopi) + ! write(20000+ipart,*) var_tip + ifp_per=ifp_per+1 + if(ifp_per.gt.nfp_per) then !<=increase the buffer for periodic boundary footprints + !$omp critical + allocate(dummy2d(2,ifp_per-1)) + !$omp end critical + dummy2d=zpoipl_per(:,1:ifp_per-1) + !$omp critical + deallocate(zpoipl_per) + !$omp end critical + nfp_per=nfp_per+nfp + !$omp critical + allocate(zpoipl_per(2,nfp_per)) + !$omp end critical + zpoipl_per(:,1:ifp_per-1)=dummy2d + !$omp critical + deallocate(dummy2d) + !$omp end critical + endif + zpoipl_per(:,ifp_per)=var_tip(1:2) endif - zpoipl_per(:,ifp_per)=var_tip(1:2) endif - endif - ! End periodic boundary footprint detection and interpolation + ! End periodic boundary footprint detection and interpolation - ! Cut classification into regular or chaotic - if (kt == ntcut) then - regular = .True. + ! Cut classification into regular or chaotic + if (kt == ntcut) then + regular = .True. - if(ifp_per > 0) then + if(ifp_per > 0) then - call fract_dimension(ifp_per,zpoipl_per(:,1:ifp_per),fraction) + call fract_dimension(ifp_per,zpoipl_per(:,1:ifp_per),fraction) - if(fraction.gt.0.2d0) then - print *, ipart, ' chaotic per ', ifp_per - regular = .False. - else - print *, ipart, ' regular per', ifp_per + if(fraction.gt.0.2d0) then + print *, ipart, ' chaotic per ', ifp_per + regular = .False. + else + print *, ipart, ' regular per', ifp_per + endif endif - endif - if(ifp_tip > 0) then + if(ifp_tip > 0) then - call fract_dimension(ifp_tip,zpoipl_tip(:,1:ifp_tip),fraction) + call fract_dimension(ifp_tip,zpoipl_tip(:,1:ifp_tip),fraction) - if(fraction.gt.0.2d0) then - print *, ipart, ' chaotic tip ', ifp_tip - regular = .False. - iclass(3,ipart) = 2 - else - print *, ipart, ' regular tip ', ifp_tip - iclass(3,ipart) = 1 + if(fraction.gt.0.2d0) then + print *, ipart, ' chaotic tip ', ifp_tip + regular = .False. + iclass(3,ipart) = 2 + else + print *, ipart, ' regular tip ', ifp_tip + iclass(3,ipart) = 1 + endif endif - endif - ! Store fractal classification in result - if(regular) then - class_result%fractal = 1 - else - class_result%fractal = 2 - endif + ! Store fractal classification in result + if(regular) then + class_result%fractal = 1 + else + class_result%fractal = 2 + endif - if(class_plot) then - call output_minkowsky_class(ipart, regular, passing) - ierr=1 + if(class_plot) then + call output_minkowsky_class(ipart, regular, passing) + ierr=1 + endif endif - endif - ! - if(ierr.ne.0) then - if(class_plot .and. .not. passing) then - call output_jpar_class(ipart, ijpar) - call output_topological_class(ipart, ideal) - exit + ! + if(ierr.ne.0) then + if(class_plot .and. .not. passing) then + call output_jpar_class(ipart, ijpar) + call output_topological_class(ipart, ideal) + exit + endif endif + ! write(999, *) kt*dtaumin/v0, z + enddo + if(ierr.ne.0) exit + if(passing) then + !$omp atomic + confpart_pass(it)=confpart_pass(it)+1.d0 + else + !$omp atomic + confpart_trap(it)=confpart_trap(it)+1.d0 endif - ! write(999, *) kt*dtaumin/v0, z enddo - if(ierr.ne.0) exit - if(passing) then - !$omp atomic - confpart_pass(it)=confpart_pass(it)+1.d0 - else - !$omp atomic - confpart_trap(it)=confpart_trap(it)+1.d0 + + !$omp critical + zend(:,ipart) = z + if(isw_field_type .eq. CANFLUX) then + call can_to_vmec(z(1),z(2),z(3),zend(2,ipart),zend(3,ipart)) + elseif(isw_field_type .eq. BOOZER) then + call boozer_to_vmec(z(1),z(2),z(3),zend(2,ipart),zend(3,ipart)) endif - enddo - - !$omp critical - zend(:,ipart) = z - if(isw_field_type .eq. CANFLUX) then - call can_to_vmec(z(1),z(2),z(3),zend(2,ipart),zend(3,ipart)) - elseif(isw_field_type .eq. BOOZER) then - call boozer_to_vmec(z(1),z(2),z(3),zend(2,ipart),zend(3,ipart)) - endif - orbit_exit_code(ipart) = class_result%exit_code - if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) then - times_lost(ipart) = ieee_value(0.0_dp, ieee_quiet_nan) - else - times_lost(ipart) = kt*dtaumin/v0 - end if - deallocate(zpoipl_tip, zpoipl_per) - !$omp end critical - if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) then - do it_f = it, ntimstep -!$omp atomic update - unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 - end do - end if - ! close(unit=10000+ipart) - ! close(unit=10000+ipart) -end subroutine trace_orbit_with_classifiers - -pure subroutine classify_classifier_exit(ierr, mode, lost, exit_code) - integer, intent(in) :: ierr, mode - logical, intent(out) :: lost - integer, intent(out) :: exit_code - - lost = .false. - if (mode <= 0 .or. ierr == SYMPLECTIC_STEP_BOUNDARY) then - lost = .true. - exit_code = ORBIT_EXIT_LCFS - return - end if - - select case (ierr) - case (SYMPLECTIC_STEP_OUTSIDE_DOMAIN) - exit_code = ORBIT_EXIT_NUMERICAL_DOMAIN - case (SYMPLECTIC_STEP_MAXITER) - exit_code = ORBIT_EXIT_NUMERICAL_MAXITER - case (SYMPLECTIC_STEP_LINEAR_SOLVE) - exit_code = ORBIT_EXIT_NUMERICAL_LINEAR - case default - exit_code = ORBIT_EXIT_NUMERICAL_EVENT - end select -end subroutine classify_classifier_exit - - -subroutine output_lost_orbit_starting_data(ipart, passing) - integer, intent(in) :: ipart - logical, intent(in) :: passing - - if(passing) then - call write_output_line(iaaa_prp, ipart) - else - call write_output_line(iaaa_prt, ipart) - endif -end subroutine output_lost_orbit_starting_data - - -subroutine output_minkowsky_class(ipart, regular, passing) - integer, intent(in) :: ipart - logical, intent(in) :: regular, passing - - if(regular) then + orbit_exit_code(ipart) = class_result%exit_code + if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) then + times_lost(ipart) = ieee_value(0.0_dp, ieee_quiet_nan) + else + times_lost(ipart) = kt*dtaumin/v0 + end if + deallocate(zpoipl_tip, zpoipl_per) + !$omp end critical + if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) then + do it_f = it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do + end if + ! close(unit=10000+ipart) + ! close(unit=10000+ipart) + end subroutine trace_orbit_with_classifiers + + pure subroutine classify_classifier_exit(ierr, mode, lost, exit_code) + integer, intent(in) :: ierr, mode + logical, intent(out) :: lost + integer, intent(out) :: exit_code + + lost = .false. + if (mode <= 0 .or. ierr == SYMPLECTIC_STEP_BOUNDARY) then + lost = .true. + exit_code = ORBIT_EXIT_LCFS + return + end if + + select case (ierr) + case (SYMPLECTIC_STEP_OUTSIDE_DOMAIN) + exit_code = ORBIT_EXIT_NUMERICAL_DOMAIN + case (SYMPLECTIC_STEP_MAXITER) + exit_code = ORBIT_EXIT_NUMERICAL_MAXITER + case (SYMPLECTIC_STEP_LINEAR_SOLVE) + exit_code = ORBIT_EXIT_NUMERICAL_LINEAR + case default + exit_code = ORBIT_EXIT_NUMERICAL_EVENT + end select + end subroutine classify_classifier_exit + + + subroutine output_lost_orbit_starting_data(ipart, passing) + integer, intent(in) :: ipart + logical, intent(in) :: passing + if(passing) then - call write_output_line(iaaa_rep, ipart) + call write_output_line(iaaa_prp, ipart) else - call write_output_line(iaaa_ret, ipart) + call write_output_line(iaaa_prt, ipart) endif - else - if(passing) then - call write_output_line(iaaa_stp, ipart) + end subroutine output_lost_orbit_starting_data + + + subroutine output_minkowsky_class(ipart, regular, passing) + integer, intent(in) :: ipart + logical, intent(in) :: regular, passing + + if(regular) then + if(passing) then + call write_output_line(iaaa_rep, ipart) + else + call write_output_line(iaaa_ret, ipart) + endif else - call write_output_line(iaaa_stt, ipart) + if(passing) then + call write_output_line(iaaa_stp, ipart) + else + call write_output_line(iaaa_stt, ipart) + endif endif - endif -end subroutine output_minkowsky_class + end subroutine output_minkowsky_class -subroutine output_jpar_class(ipart, ijpar) - integer, intent(in) :: ipart, ijpar + subroutine output_jpar_class(ipart, ijpar) + integer, intent(in) :: ipart, ijpar - select case(ijpar) + select case(ijpar) case(0) - call write_output_line(iaaa_jer, ipart) + call write_output_line(iaaa_jer, ipart) case(1) - call write_output_line(iaaa_jre, ipart) + call write_output_line(iaaa_jre, ipart) case(2) - call write_output_line(iaaa_jst, ipart) - end select -end subroutine output_jpar_class + call write_output_line(iaaa_jst, ipart) + end select + end subroutine output_jpar_class -subroutine output_topological_class(ipart, ideal) - integer, intent(in) :: ipart, ideal + subroutine output_topological_class(ipart, ideal) + integer, intent(in) :: ipart, ideal - select case(ideal) + select case(ideal) case(0) - call write_output_line(iaaa_ier, ipart) + call write_output_line(iaaa_ier, ipart) case(1) - call write_output_line(iaaa_ire, ipart) + call write_output_line(iaaa_ire, ipart) case(2) - call write_output_line(iaaa_ist, ipart) - end select -end subroutine output_topological_class - - -subroutine write_output_line(iunit, ipart) - integer, intent(in) :: iunit, ipart - - !$omp critical - write (iunit, *) zstart(2,ipart), zstart(5,ipart), trap_par(ipart) - !$omp end critical -end subroutine write_output_line - - -! Write classification results to fort.* files -! This subroutine centralizes all classification file I/O -! Note: Only writes classifications that were computed (non-zero values) -subroutine write_classification_results(ipart, class_result) - integer, intent(in) :: ipart - type(classification_result_t), intent(in) :: class_result - logical :: regular - - if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) return - - ! Write lost orbits - if(class_result%lost) then - call output_lost_orbit_starting_data(ipart, class_result%passing) - return - endif - - ! Write fractal classification if computed - if(class_result%fractal /= 0) then - regular = (class_result%fractal == 1) - call output_minkowsky_class(ipart, regular, class_result%passing) - endif - - ! Write J_parallel and topological classification if computed - ! These are only done for trapped orbits - if(.not. class_result%passing) then - if(class_result%jpar /= 0) then - call output_jpar_class(ipart, class_result%jpar) + call write_output_line(iaaa_ist, ipart) + end select + end subroutine output_topological_class + + + subroutine write_output_line(iunit, ipart) + integer, intent(in) :: iunit, ipart + + !$omp critical + write (iunit, *) zstart(2,ipart), zstart(5,ipart), trap_par(ipart) + !$omp end critical + end subroutine write_output_line + + + ! Write classification results to fort.* files + ! This subroutine centralizes all classification file I/O + ! Note: Only writes classifications that were computed (non-zero values) + subroutine write_classification_results(ipart, class_result) + integer, intent(in) :: ipart + type(classification_result_t), intent(in) :: class_result + logical :: regular + + if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) return + + ! Write lost orbits + if(class_result%lost) then + call output_lost_orbit_starting_data(ipart, class_result%passing) + return + endif + + ! Write fractal classification if computed + if(class_result%fractal /= 0) then + regular = (class_result%fractal == 1) + call output_minkowsky_class(ipart, regular, class_result%passing) endif - if(class_result%topology /= 0) then - call output_topological_class(ipart, class_result%topology) + + ! Write J_parallel and topological classification if computed + ! These are only done for trapped orbits + if(.not. class_result%passing) then + if(class_result%jpar /= 0) then + call output_jpar_class(ipart, class_result%jpar) + endif + if(class_result%topology /= 0) then + call output_topological_class(ipart, class_result%topology) + endif endif - endif -end subroutine write_classification_results + end subroutine write_classification_results end module classification diff --git a/src/diag_counters.f90 b/src/diag_counters.f90 index eeb66a78..ba700240 100644 --- a/src/diag_counters.f90 +++ b/src/diag_counters.f90 @@ -9,15 +9,16 @@ module diag_counters !> critical. Totals are a plain reduction over the thread columns, taken !> outside the hot path. use, intrinsic :: iso_fortran_env, only: int64 -!$ use omp_lib, only: omp_get_max_threads, omp_get_thread_num + !$ use omp_lib, only: omp_get_max_threads, omp_get_thread_num implicit none private public :: EVT_NEWTON1_MAXIT, EVT_NEWTON2_MAXIT, EVT_RK_GAUSS_MAXIT, & - EVT_RK_LOBATTO_MAXIT, EVT_FIXPOINT_MAXIT, EVT_R_NEGATIVE, & - EVT_FO_LOSS, EVT_FO_FAULT, EVT_MIDPOINT_MAXIT, N_EVENT + EVT_RK_LOBATTO_MAXIT, EVT_FIXPOINT_MAXIT, EVT_R_NEGATIVE, & + EVT_FO_LOSS, EVT_FO_FAULT, EVT_MIDPOINT_MAXIT, & + EVT_WARNING_STEP_SKIP, N_EVENT public :: diag_counters_init, count_event, diag_counters_total, & - diag_counters_reset, event_name + diag_counters_reset, event_name integer, parameter :: EVT_NEWTON1_MAXIT = 1 integer, parameter :: EVT_NEWTON2_MAXIT = 2 @@ -31,12 +32,13 @@ module diag_counters integer, parameter :: EVT_FO_LOSS = 7 integer, parameter :: EVT_FO_FAULT = 8 integer, parameter :: EVT_MIDPOINT_MAXIT = 9 - integer, parameter :: N_EVENT = 9 + integer, parameter :: EVT_WARNING_STEP_SKIP = 10 + integer, parameter :: N_EVENT = 10 ! Whole cache lines per thread column keep neighbouring threads from sharing ! a line. The event id indexes within a column; STRIDE >= N_EVENT. integer, parameter :: STRIDE = 16 - integer(int64), allocatable :: counts(:, :) ! (STRIDE, 0:nthreads-1) + integer(int64), allocatable :: counts(:, :) ! (STRIDE, 0:nthreads-1) contains @@ -44,7 +46,7 @@ subroutine diag_counters_init() integer :: nthreads nthreads = 1 -!$ nthreads = omp_get_max_threads() + !$ nthreads = omp_get_max_threads() if (allocated(counts)) deallocate (counts) allocate (counts(STRIDE, 0:nthreads - 1)) counts = 0_int64 @@ -58,7 +60,7 @@ subroutine count_event(id) if (.not. allocated(counts)) return tid = 0 -!$ tid = omp_get_thread_num() + !$ tid = omp_get_thread_num() counts(id, tid) = counts(id, tid) + 1_int64 end subroutine count_event @@ -98,6 +100,8 @@ function event_name(id) result(name) name = 'fo_fault' case (EVT_MIDPOINT_MAXIT) name = 'midpoint_maxit' + case (EVT_WARNING_STEP_SKIP) + name = 'warning_step_skip' case default name = 'unknown' end select diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 5b4e78a8..6fe69d9d 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -26,19 +26,20 @@ module simple_main ORBIT_EXIT_NUMERICAL_EVENT, ORBIT_EXIT_NUMERICAL_FULL_ORBIT use params, only: canonical_grid_nr, canonical_grid_ntheta, & canonical_grid_nphi, canonical_ode_relerr - use diag_counters, only: diag_counters_init + use diag_counters, only: diag_counters_init, count_event, & + EVT_WARNING_STEP_SKIP use progress_monitor, only: progress_init, progress_tick, progress_finalize use restart_mod, only: particle_done, read_restart_data, restore_confined_counts use chartmap_metadata, only: chartmap_metadata_t, read_chartmap_metadata use reference_coordinates, only: ref_coords use stl_wall_intersection, only: stl_wall_t, stl_wall_init, & - stl_wall_finalize, & - stl_wall_first_hit_segment_with_normal + stl_wall_finalize, & + stl_wall_first_hit_segment_with_normal use libneo_coordinates, only: chartmap_coordinate_system_t use orbit_symplectic_base, only: SYMPLECTIC_STEP_BOUNDARY, & SYMPLECTIC_STEP_OUTSIDE_DOMAIN, SYMPLECTIC_STEP_MAXITER, & SYMPLECTIC_STEP_LINEAR_SOLVE, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & - SYMPLECTIC_STEP_BOUNDARY_LIMITED + SYMPLECTIC_STEP_BOUNDARY_LIMITED, symplectic_newton_warning_mode implicit none @@ -54,13 +55,13 @@ module simple_main subroutine main use params, only: read_config, netcdffile, ns_s, ns_tp, multharm, & - integmode, params_init, swcoll, generate_start_only, & - isw_field_type, field_input, startmode, & - ntestpart, ntimstep, coord_input, restart + integmode, params_init, swcoll, generate_start_only, & + isw_field_type, field_input, startmode, & + ntestpart, ntimstep, coord_input, restart use timing, only: init_timer, print_phase_time use magfie_sub, only: TEST, VMEC, SPECTRE, init_magfie use samplers, only: init_starting_surf, sample_spectre_surface, & - init_spectre_start_bounds + init_spectre_start_bounds use version, only: simple_version use field_boozer_chartmap, only: is_boozer_chartmap use field, only: is_spectre_file @@ -181,7 +182,7 @@ subroutine main character(32) :: gpu_bench_env integer :: gpu_bench_len, gpu_bench_stat call get_environment_variable('SIMPLE_GPU_BENCH', gpu_bench_env, & - gpu_bench_len, gpu_bench_stat) + gpu_bench_len, gpu_bench_stat) if (gpu_bench_stat == 0 .and. gpu_bench_len > 0) then call trace_compare_gpu(norb) return @@ -232,7 +233,7 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) use field_boozer_chartmap, only: boozer_chartmap_field_t, is_boozer_chartmap use timing, only: print_phase_time use magfie_sub, only: TEST, CANFLUX, VMEC, BOOZER, MEISS, ALBERT, & - REFCOORDS, SPECTRE, set_magfie_refcoords_field + REFCOORDS, SPECTRE, set_magfie_refcoords_field use field_splined, only: splined_field_t, create_splined_field use field_vmec, only: vmec_field_t use reference_coordinates, only: init_reference_coordinates, ref_coords @@ -258,7 +259,7 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) ! TEST field is analytic - no VMEC or field files needed if (isw_field_type == TEST) then - self%fper = twopi ! Full torus for analytic tokamak + self%fper = twopi ! Full torus for analytic tokamak call print_phase_time('TEST field mode - no input files required') else if (use_spectre) then call init_spectre_field(self) @@ -267,7 +268,7 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) ! Boozer chartmap: file-based, no VMEC initialization needed call init_reference_coordinates(coord_input) call print_phase_time('Reference coordinate system '// & - 'initialization completed') + 'initialization completed') block use libneo_coordinates, only: chartmap_coordinate_system_t @@ -288,50 +289,50 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) end if else vmec_equilibrium_file = select_vmec_equilibrium_file(vmec_file, & - field_input, & - coord_input) + field_input, & + coord_input) call init_vmec(vmec_equilibrium_file, ans_s, ans_tp, amultharm, & - self%fper) + self%fper) call print_phase_time('VMEC initialization completed') call init_reference_coordinates(coord_input) call print_phase_time('Reference coordinate system '// & - 'initialization completed') + 'initialization completed') call init_stl_wall_if_enabled(coord_input) call print_phase_time('STL wall initialization completed') if (self%integmode >= 0) then - if (trim(field_input) == '') then - print *, 'simple_main.init_field: field_input must be set (see ', & - 'params.apply_config_aliases)' - error stop - end if - - call field_from_file(field_input, field_temp) - call print_phase_time('Field from file loading completed') - - if (isw_field_type == REFCOORDS) then - select type (field_temp) - type is (splined_field_t) - call set_magfie_refcoords_field(field_temp) - type is (vmec_field_t) - block - type(splined_field_t) :: splined_vmec - call create_splined_field(field_temp, ref_coords, & - splined_vmec) - call set_magfie_refcoords_field(splined_vmec) - end block - class default - print *, & - 'simple_main.init_field: REFCOORDS requires '// & - 'a splined field' - print *, & - 'Supported inputs: coils (auto-splined) or VMEC wout', & - ' (splined onto coord_input)' + if (trim(field_input) == '') then + print *, 'simple_main.init_field: field_input must be set (see ', & + 'params.apply_config_aliases)' error stop - end select - end if + end if + + call field_from_file(field_input, field_temp) + call print_phase_time('Field from file loading completed') + + if (isw_field_type == REFCOORDS) then + select type (field_temp) + type is (splined_field_t) + call set_magfie_refcoords_field(field_temp) + type is (vmec_field_t) + block + type(splined_field_t) :: splined_vmec + call create_splined_field(field_temp, ref_coords, & + splined_vmec) + call set_magfie_refcoords_field(splined_vmec) + end block + class default + print *, & + 'simple_main.init_field: REFCOORDS requires '// & + 'a splined field' + print *, & + 'Supported inputs: coils (auto-splined) or VMEC wout', & + ' (splined onto coord_input)' + error stop + end select + end if end if end if @@ -358,7 +359,7 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode) transformation_relerr=canonical_ode_relerr) call print_phase_time('Canonical field initialization completed') else if (isw_field_type == CANFLUX .or. isw_field_type == BOOZER .or. & - isw_field_type == MEISS .or. isw_field_type == ALBERT) then + isw_field_type == MEISS .or. isw_field_type == ALBERT) then call init_field_can(isw_field_type, field_temp, canonical_grid_nr, & canonical_grid_ntheta, canonical_grid_nphi, & canonical_ode_relerr) @@ -379,8 +380,8 @@ subroutine init_spectre_field(self) use new_vmec_stuff_mod, only: nper, rmajor use util, only: twopi use params, only: field_input, spectre_ncon_r, spectre_ncon_th, & - spectre_ncon_phi, spectre_ncon_order, & - spectre_ncon_ode_max_steps + spectre_ncon_phi, spectre_ncon_order, & + spectre_ncon_ode_max_steps use params, only: spectre_ncon_ode_relerr use timing, only: print_phase_time use orbit_symplectic_base, only: sympl_rmax @@ -427,9 +428,9 @@ subroutine init_spectre_field(self) ! extended linearly, so iterates out there stay finite. sympl_rmax = real(sf%data%Mvol + 1, dp) call set_spectre_construction_grid(spectre_ncon_r, spectre_ncon_th, & - spectre_ncon_phi, spectre_ncon_order, & - spectre_ncon_ode_max_steps, & - spectre_ncon_ode_relerr) + spectre_ncon_phi, spectre_ncon_order, & + spectre_ncon_ode_max_steps, & + spectre_ncon_ode_relerr) call init_field_can(SPECTRE, sf) call print_phase_time('SPECTRE per-volume canonical construction completed') end if @@ -477,9 +478,9 @@ subroutine init_stl_wall_if_enabled(coord_file) end subroutine init_stl_wall_if_enabled function select_vmec_equilibrium_file(vmec_file_default, field_file, coord_file) & - result(vmec_file) + result(vmec_file) use libneo_coordinates, only: detect_refcoords_file_type, & - refcoords_file_vmec_wout + refcoords_file_vmec_wout character(*), intent(in) :: vmec_file_default character(*), intent(in) :: field_file character(*), intent(in) :: coord_file @@ -493,7 +494,7 @@ function select_vmec_equilibrium_file(vmec_file_default, field_file, coord_file) if (len_trim(field_file) > 0) then if (ends_with_nc(field_file)) then call detect_refcoords_file_type(trim(field_file), file_type, & - ierr, message) + ierr, message) if (ierr == 0 .and. file_type == refcoords_file_vmec_wout) then vmec_file = trim(field_file) return @@ -504,7 +505,7 @@ function select_vmec_equilibrium_file(vmec_file_default, field_file, coord_file) if (len_trim(coord_file) > 0) then if (ends_with_nc(coord_file)) then call detect_refcoords_file_type(trim(coord_file), file_type, & - ierr, message) + ierr, message) if (ierr == 0 .and. file_type == refcoords_file_vmec_wout) then vmec_file = trim(coord_file) return @@ -529,7 +530,7 @@ end function ends_with_nc subroutine trace_parallel(norb) use netcdf_orbit_output, only: init_orbit_netcdf, close_orbit_netcdf, & - write_orbit_to_netcdf + write_orbit_to_netcdf #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT use params, only: debug #endif @@ -542,10 +543,10 @@ subroutine trace_parallel(norb) call init_orbit_netcdf(ntestpart, ntimstep) end if -!$omp parallel firstprivate(norb) private(traj, times, i) + !$omp parallel firstprivate(norb) private(traj, times, i) allocate (traj(5, ntimstep), times(ntimstep)) -!$omp do + !$omp do do i = 1, ntestpart if (allocated(particle_done)) then if (particle_done(i)) then @@ -556,26 +557,26 @@ subroutine trace_parallel(norb) #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT if (debug) then -!$omp critical + !$omp critical kpart = kpart + 1 print *, kpart, ' / ', ntestpart, 'particle: ', i, 'thread: ', & omp_get_thread_num() -!$omp end critical + !$omp end critical end if #endif call trace_orbit(norb, i, traj, times) if (output_orbits_macrostep) then -!$omp critical + !$omp critical call write_orbit_to_netcdf(i, traj, times) -!$omp end critical + !$omp end critical end if call progress_tick end do -!$omp end do -!$omp end parallel + !$omp end do + !$omp end parallel if (output_orbits_macrostep) then call close_orbit_netcdf() @@ -609,7 +610,7 @@ subroutine trace_compare_gpu(norb) if (isw_field_type /= BOOZER .or. integmode /= EXPL_IMPL_EULER .or. swcoll .or. & len_trim(wall_input) > 0 .or. class_plot .or. fast_class .or. generate_start_only) then error stop "simple_main.trace_compare_gpu: SIMPLE_GPU_BENCH requires Boozer, " // & - "EXPL_IMPL_EULER, and no wall/collision/classifier options" + "EXPL_IMPL_EULER, and no wall/collision/classifier options" end if call sync_boozer_state @@ -658,7 +659,7 @@ subroutine trace_compare_gpu(norb) ! GPU kernel t0 = omp_get_wtime() call trace_orbits_gpu(si_gpu, f_gpu, ntestpart, ntimstep, ntau_macro, & - gpu_loss, gpu_zend) + gpu_loss, gpu_zend) t1 = omp_get_wtime() t_gpu = t1 - t0 @@ -702,16 +703,16 @@ subroutine classify_parallel(norb) integer :: i type(classification_result_t) :: class_result -!$omp parallel firstprivate(norb) private(class_result, i) -!$omp do + !$omp parallel firstprivate(norb) private(class_result, i) + !$omp do do i = 1, ntestpart #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT if (debug) then -!$omp critical + !$omp critical kpart = kpart + 1 print *, kpart, ' / ', ntestpart, 'particle: ', i, 'thread: ', & omp_get_thread_num() -!$omp end critical + !$omp end critical end if #endif @@ -724,14 +725,14 @@ subroutine classify_parallel(norb) ! iclass already populated by trace_orbit_with_classifiers ! Other results (zend, times_lost, trap_par, perp_inv) also already stored end do -!$omp end do -!$omp end parallel + !$omp end do + !$omp end parallel end subroutine classify_parallel subroutine print_parameters print *, 'tau: ', dtau, dtaumin, min(dabs(mod(dtau, dtaumin)), & - dabs(mod(dtau, dtaumin) - & - dtaumin))/dtaumin, ntau + dabs(mod(dtau, dtaumin) - & + dtaumin))/dtaumin, ntau print *, 'v0 = ', v0 end subroutine print_parameters @@ -745,9 +746,9 @@ end subroutine read_profiles_config subroutine init_collisions use params, only: am1, am2, Z1, Z2, facE_al, dchichi, slowrate, & - dchichi_norm, slowrate_norm, v0 + dchichi_norm, slowrate_norm, v0 use simple_profiles, only: Te_scale, Ti1_scale, Ti2_scale, & - ni1_scale, ni2_scale + ni1_scale, ni2_scale real(dp) :: v0_coll, ealpha real(dp) :: densi1, densi2, tempi1, tempi2, tempe @@ -760,8 +761,8 @@ subroutine init_collisions tempe = Te_scale call loacol_alpha(am1, am2, Z1, Z2, densi1, densi2, tempi1, tempi2, tempe, & - ealpha, v0_coll, dchichi, slowrate, dchichi_norm, & - slowrate_norm) + ealpha, v0_coll, dchichi, slowrate, dchichi_norm, & + slowrate_norm) if (abs(v0_coll - v0) > 1d-6) then error stop 'simple_main.init_collisions: v0_coll != v0' @@ -772,7 +773,7 @@ end subroutine init_collisions subroutine sample_particles(xstart_is_integ_coords) use samplers, only: sample, START_FILE, sample_grid, & - sample_surface_fieldline, sample_surface_fieldline_from_integ + sample_surface_fieldline, sample_surface_fieldline_from_integ logical, intent(in), optional :: xstart_is_integ_coords logical :: convert_surface_starts @@ -825,7 +826,7 @@ subroutine sample_particles_test_field !> TEST field uses (r, theta, phi) coordinates with B0=1, R0=1, a=0.5, iota=1. !> bmod = B0 * (1 - r/R0 * cos(theta)) use params, only: ntestpart, sbeg, bmod00, bmin, bmax, zstart, & - reset_seed_if_deterministic + reset_seed_if_deterministic real(dp), parameter :: B0 = 1.0d0, R0 = 1.0d0, a = 0.5d0 real(dp) :: r_start, tmp_rand @@ -960,15 +961,15 @@ end function classify_orbit_exit subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) use classification, only: trace_orbit_with_classifiers, & - classification_result_t, & - write_classification_results + classification_result_t, & + write_classification_results use magfie_sub, only: SPECTRE use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan type(tracer_t), intent(inout) :: anorb integer, intent(in) :: ipart - real(dp), intent(out) :: orbit_traj(:, :) ! (5, ntimstep) - real(dp), intent(out) :: orbit_times(:) ! (ntimstep) + real(dp), intent(out) :: orbit_traj(:, :) ! (5, ntimstep) + real(dp), intent(out) :: orbit_times(:) ! (ntimstep) real(dp), dimension(5) :: z real(dp) :: u_ref_prev(3), x_prev(3), x_prev_m(3), exit_step @@ -1024,7 +1025,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) if (isw_field_type == SPECTRE) then if (integmode > 0) then call trace_orbit_spectre_sympl(anorb, ipart, z, passing, & - orbit_traj, orbit_times) + orbit_traj, orbit_times) else call trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) end if @@ -1036,7 +1037,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) times_lost(ipart) = -1.d0 orbit_exit_code(ipart) = ORBIT_EXIT_SKIPPED do it = 1, ntimstep -!$omp atomic update + !$omp atomic update confpart_pass(it) = confpart_pass(it) + 1.d0 end do return @@ -1075,7 +1076,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) end if else do it_f = it, ntimstep -!$omp atomic update + !$omp atomic update unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 end do faulted = .true. @@ -1099,7 +1100,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) end do end if -!$omp critical + !$omp critical call integ_to_ref(z(1:3), zend(1:3, ipart)) zend(4:5, ipart) = z(4:5) if (physical_exit) then @@ -1109,7 +1110,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) else times_lost(ipart) = kt*dtaumin/v0 end if -!$omp end critical + !$omp end critical end subroutine trace_orbit subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) @@ -1118,8 +1119,8 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) !> confined for the full trace, reflects at forbidden interfaces, or is !> lost at the outermost interface. Every crossing/reflection is logged. use spectre_orbit, only: spectre_orbit_state_t, spectre_event_t, & - spectre_state_reset, orbit_timestep_spectre, & - SPECTRE_OK, SPECTRE_BOUNDARY + spectre_state_reset, orbit_timestep_spectre, & + SPECTRE_OK, SPECTRE_BOUNDARY use magfie_sub, only: spectre_field use interface_crossing, only: crossing_log_record use params, only: crossing_level @@ -1133,7 +1134,7 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) type(spectre_orbit_state_t) :: state type(spectre_event_t) :: event - integer :: it, ktau, ierr_orbit, it_final + integer :: it, it_f, ktau, ierr_orbit, it_final integer(8) :: kt real(dp) :: t_stop @@ -1146,7 +1147,7 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) if (it >= 2) then do ktau = 1, ntau_macro(it) call orbit_timestep_spectre(state, z, dtaumin, relerr, & - crossing_level, ierr_orbit, event) + crossing_level, ierr_orbit, event) if (event%occurred) then call crossing_log_record(ipart, & (real(kt, dp) + event%t_frac)*dtaumin/v0, event%info) @@ -1158,6 +1159,15 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) if (ierr_orbit /= SPECTRE_OK) then it_final = it + if (ierr_orbit == SPECTRE_BOUNDARY) then + orbit_exit_code(ipart) = ORBIT_EXIT_LCFS + else + orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT + do it_f = it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do + end if exit end if @@ -1176,19 +1186,22 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) if (ierr_orbit == SPECTRE_BOUNDARY) then t_stop = (real(kt, dp) + event%t_frac)*dtaumin/v0 + else if (ierr_orbit /= SPECTRE_OK) then + t_stop = ieee_value(0.0_dp, ieee_quiet_nan) else t_stop = real(kt, dp)*dtaumin/v0 + orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED end if -!$omp critical + !$omp critical call integ_to_ref(z(1:3), zend(1:3, ipart)) zend(4:5, ipart) = z(4:5) times_lost(ipart) = t_stop -!$omp end critical + !$omp end critical end subroutine trace_orbit_spectre subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & - orbit_times) + orbit_times) !> Per-volume symplectic guiding-center trace for SPECTRE (integmode > 0). !> Each microstep resolves interface events by an exact-landing substep of !> the same implicit scheme, applies the crossing map, and re-canonicalizes @@ -1196,8 +1209,9 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & !> interface, matching the RK45 path; CROSS_STOP remains only as the !> pathological non-convergence fallback inside the microstepper. use spectre_sympl_orbit, only: sympl_spectre_state_t, sympl_spectre_reset, & - orbit_microstep_sympl_spectre, & - SYMPL_SPECTRE_OK, SYMPL_SPECTRE_SKIM + orbit_microstep_sympl_spectre, & + SYMPL_SPECTRE_OK, SYMPL_SPECTRE_LOSS, & + SYMPL_SPECTRE_STOP, SYMPL_SPECTRE_SKIM use field_can_spectre, only: spectre_mvol, set_spectre_volume_lock use params, only: crossing_level use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan @@ -1210,12 +1224,12 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & real(dp), intent(out) :: orbit_times(:) type(sympl_spectre_state_t) :: state - integer :: it, ktau, ierr_orbit, it_final + integer :: it, it_f, ktau, ierr_orbit, it_final integer(8) :: kt real(dp) :: t_stop, t_frac call sympl_spectre_reset(state, anorb%si, spectre_mvol, integmode, & - crossing_level) + crossing_level) kt = 0 it_final = 0 @@ -1226,10 +1240,10 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & if (it >= 2) then do ktau = 1, ntau_macro(it) call orbit_microstep_sympl_spectre(state, anorb%si, anorb%f, & - ipart, & - real(kt, dp)*dtaumin/v0, & - dtaumin/v0, ierr_orbit, & - t_frac) + ipart, & + real(kt, dp)*dtaumin/v0, & + dtaumin/v0, ierr_orbit, & + t_frac) if (ierr_orbit /= SYMPL_SPECTRE_OK) exit kt = kt + 1 end do @@ -1237,6 +1251,18 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & if (ierr_orbit /= SYMPL_SPECTRE_OK) then it_final = it + select case (ierr_orbit) + case (SYMPL_SPECTRE_LOSS) + orbit_exit_code(ipart) = ORBIT_EXIT_LCFS + case (SYMPL_SPECTRE_STOP) + orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_FULL_ORBIT + do it_f = it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do + case (SYMPL_SPECTRE_SKIM) + orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED + end select exit end if @@ -1269,10 +1295,13 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & if (ierr_orbit == SYMPL_SPECTRE_OK) then t_stop = real(kt, dp)*dtaumin/v0 + orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED else if (ierr_orbit == SYMPL_SPECTRE_SKIM) then ! Mirror-confined at an interior interface: cannot be lost, so record ! it as confined (times_lost = trace_time) rather than at its stop. t_stop = trace_time + else if (ierr_orbit == SYMPL_SPECTRE_STOP) then + t_stop = ieee_value(0.0_dp, ieee_quiet_nan) else t_stop = (real(kt, dp) + t_frac)*dtaumin/v0 end if @@ -1285,11 +1314,11 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & ! The next marker on this thread starts with an unlocked field dispatch. call set_spectre_volume_lock(0) -!$omp critical + !$omp critical call integ_to_ref(z(1:3), zend(1:3, ipart)) zend(4:5, ipart) = z(4:5) times_lost(ipart) = t_stop -!$omp end critical + !$omp end critical end subroutine trace_orbit_spectre_sympl pure subroutine locate_linear_lcfs(z_start, z_end, field_period, z_event, & @@ -1344,8 +1373,16 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) if (present(exit_step)) then exit_step = real(kt, dp) + loss_fraction end if + else if (symplectic_newton_warning_mode) then + ! The full-orbit pusher retains its last resolved state. + ! Default production mode advances the clock past this + ! one unresolved microstep and retries the next step; + ! numerical inversion faults are never physical losses. + call count_event(EVT_WARNING_STEP_SKIP) + z = z_step_start + ierr_orbit = 0 end if - exit + if (ierr_orbit .ne. 0) exit end if kt = kt + 1 cycle @@ -1370,8 +1407,17 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) anorb%si%last_step_fraction exit end if + if (ierr_orbit .ne. 0 .and. & + symplectic_newton_warning_mode) then + ! advance_symplectic_with_retry restored the last accepted + ! state. Skip only the unresolved microstep and continue; + ! strict mode still reports the numerical exit. + call count_event(EVT_WARNING_STEP_SKIP) + ierr_orbit = 0 + else if (ierr_orbit == 0) then + call to_standard_z_coordinates(anorb, z) + end if if (ierr_orbit .ne. 0) exit - call to_standard_z_coordinates(anorb, z) end if if (ierr_orbit .ne. 0) exit if (swcoll) call collide(z, dtaumin) ! Collisions @@ -1461,11 +1507,12 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & real(dp) :: z_step_start(5), z_step_end(5), segment_duration real(dp) :: boundary_fraction real(dp) :: wall_exit_step - logical :: hit + logical :: hit, numerical_hold call integ_to_ref(z(1:3), u_ref_prev) if (present(exit_step)) exit_step = real(kt, dp) do ktau = 1, ntau_local + numerical_hold = .false. z_step_start = z if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr_orbit) @@ -1514,10 +1561,22 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & end if exit end if + if (ierr_orbit .ne. 0 .and. & + symplectic_newton_warning_mode) then + call count_event(EVT_WARNING_STEP_SKIP) + ierr_orbit = 0 + numerical_hold = .true. + else if (ierr_orbit == 0) then + call to_standard_z_coordinates(anorb, z) + end if if (ierr_orbit .ne. 0) exit - call to_standard_z_coordinates(anorb, z) end if if (ierr_orbit .ne. 0) exit + if (numerical_hold) then + if (swcoll) call collide(z, dtaumin) + kt = kt + 1 + cycle + end if z_step_end = z call integ_to_ref(z(1:3), u_ref_cur) call ref_coords%evaluate_cart(u_ref_cur, x_cur) @@ -1552,10 +1611,10 @@ subroutine increase_confined_count(it, passing) logical, intent(in) :: passing if (passing) then -!$omp atomic update + !$omp atomic update confpart_pass(it) = confpart_pass(it) + 1.d0 else -!$omp atomic update + !$omp atomic update confpart_trap(it) = confpart_trap(it) + 1.d0 end if end subroutine increase_confined_count @@ -1569,7 +1628,7 @@ subroutine compute_pitch_angle_params(z, passing, trap_par_, perp_inv_) real(dp) :: bmod -!$omp critical + !$omp critical bmod = compute_bmod(z(1:3)) if (num_surf /= 1) then call get_bminmax(z(1), bmin, bmax) @@ -1577,7 +1636,7 @@ subroutine compute_pitch_angle_params(z, passing, trap_par_, perp_inv_) passing = z(5)**2 .gt. 1.d0 - bmod/bmax trap_par_ = ((1.d0 - z(5)**2)*bmax/bmod - 1.d0)*bmin/(bmax - bmin) perp_inv_ = z(4)**2*(1.d0 - z(5)**2)/bmod -!$omp end critical + !$omp end critical end subroutine compute_pitch_angle_params function compute_bmod(z) result(bmod) @@ -1625,9 +1684,9 @@ subroutine write_output ! Sum field evaluations across all threads total_field_evaluations = 0 -!$omp parallel reduction(+:total_field_evaluations) + !$omp parallel reduction(+:total_field_evaluations) total_field_evaluations = total_field_evaluations + n_field_evaluations -!$omp end parallel + !$omp end parallel print *, "Total field evaluations: ", total_field_evaluations diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index cec83e84..48d99847 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -16,7 +16,8 @@ module spectre_sympl_orbit use field_can_mod, only: field_can_t, eval_field => evaluate, integ_to_ref, & ref_to_integ use field_can_spectre, only: set_spectre_volume_lock - use orbit_symplectic_base, only: symplectic_integrator_t + use orbit_symplectic_base, only: symplectic_integrator_t, & + symplectic_newton_warning_mode use orbit_symplectic, only: orbit_timestep_sympl, orbit_sympl_init use interface_crossing, only: apply_crossing, crossing_info_t, & crossing_log_record, CROSS_LOSS, CROSS_STOP, & @@ -25,6 +26,7 @@ module spectre_sympl_orbit sheet_gc_advance, sheet_gc_to_y, SHEET_GC_OK use spectre_fo_hybrid, only: spectre_fo_state_t, spectre_fo_enter, & spectre_fo_advance_until_exit, SPECTRE_FO_OK, SPECTRE_FO_LOSS + use diag_counters, only: count_event, EVT_WARNING_STEP_SKIP implicit none private @@ -208,6 +210,11 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call continue_fo(state, si, f, budget, used, y_out, fo_owner, & fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) then + si%dt = state%dt_std + return + end if call finish_fo_error(ierr_fo, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, state%fo%iface, 1, STOP_SHEET, & used, state%dt_std, ierr, t_frac) @@ -238,6 +245,11 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & state%sheet%iface, budget, used, y_out, fo_owner, fo_exited, & ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) then + si%dt = state%dt_std + return + end if direction = merge(1, -1, & state%sheet%owner == state%sheet%iface) call finish_fo_error(ierr_fo, si, f, ipart, & @@ -270,6 +282,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, nint(state%home_hi), & iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, iface, 1, STOP_STEP, used, & state%dt_std, ierr, t_frac) @@ -297,6 +311,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, nint(state%home_hi), & iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, & t_base_sec + dt_sec*used/state%dt_std, iface, 1, & STOP_TELEPORT, used, state%dt_std, ierr, t_frac) @@ -326,6 +342,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, nint(state%home_hi), & iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, & t_base_sec + dt_sec*used/state%dt_std, iface, 1, & STOP_STEP, used, state%dt_std, ierr, t_frac) @@ -360,6 +378,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, nint(state%home_hi), & iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, iface, direction, STOP_LANDING, & used, state%dt_std, ierr, t_frac) @@ -409,6 +429,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, state%sheet%owner, & iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, & + fo_owner, budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, & t_base_sec + dt_sec*used/state%dt_std, iface, & direction, STOP_SHEET, used, state%dt_std, ierr, & @@ -429,6 +451,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, info%vol_from, iface, & budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, iface, direction, STOP_SHEET, & used, state%dt_std, ierr, t_frac) @@ -460,6 +484,8 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call start_fo(state, si, f, y_iface, info%vol_to, iface, & budget, used, y_out, fo_owner, fo_exited, ierr_fo) if (ierr_fo /= SPECTRE_FO_OK) then + if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & + budget, used, ierr, t_frac)) exit call finish_fo_error(ierr_fo, si, f, ipart, & t_base_sec + dt_sec*used/state%dt_std, iface, direction, & STOP_EVENT_CAP, used, state%dt_std, ierr, t_frac) @@ -729,6 +755,9 @@ subroutine start_fo(state, si, f, y, owner, iface, budget, used, y_out, & integer, intent(out) :: owner_out, ierr logical, intent(out) :: exited + y_out = y + owner_out = owner + exited = .false. state%sheet%active = .false. call spectre_fo_enter(state%fo, y, owner, iface, ro0/sqrt2, ierr) if (ierr /= SPECTRE_FO_OK) then @@ -841,6 +870,52 @@ pure function unwrap_near(a, near) result(b) b = a + twopi*nint((near - a)/twopi) end function unwrap_near + logical function recover_fo_failure(state, si, f, fo_status, owner, budget, & + used, ierr, t_frac) + !> Full orbit is the final active recovery. If it also fails numerically, + !> return to its last valid guiding-centre state and use the warning-mode + !> hold fallback. A physical full-orbit loss is never intercepted. + type(sympl_spectre_state_t), intent(inout) :: state + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(in) :: fo_status, owner + real(dp), intent(inout) :: budget, used + integer, intent(inout) :: ierr + real(dp), intent(inout) :: t_frac + + recover_fo_failure = .false. + if (fo_status == SPECTRE_FO_LOSS) return + if (.not. symplectic_newton_warning_mode) return + if (state%fo%has_y) then + call set_home_volume(state, owner) + call recanonicalize(state, si, f, state%fo%last_y) + end if + recover_fo_failure = recover_warning_skip(state, budget, used, ierr, & + t_frac) + end function recover_fo_failure + + logical function recover_warning_skip(state, budget, used, ierr, t_frac) + !> Last-resort default: after all active recovery has failed, retain the + !> last valid state for this unresolved interval and resume the marker on + !> the next microstep. Strict mode returns .false. to expose the stop. + type(sympl_spectre_state_t), intent(inout) :: state + real(dp), intent(inout) :: budget, used + integer, intent(inout) :: ierr + real(dp), intent(inout) :: t_frac + + recover_warning_skip = symplectic_newton_warning_mode + if (.not. recover_warning_skip) return + + used = used + budget + budget = 0.0_dp + ierr = SYMPL_SPECTRE_OK + t_frac = 1.0_dp + call count_event(EVT_WARNING_STEP_SKIP) + + state%fo%active = .false. + state%sheet%active = .false. + end function recover_warning_skip + subroutine record_stop(si, f, ipart, t_sec, iface, direction, reason) type(symplectic_integrator_t), intent(in) :: si type(field_can_t), intent(in) :: f diff --git a/test/tests/test_sympl_testfield.f90 b/test/tests/test_sympl_testfield.f90 index b2d427ba..2bf4a2b3 100644 --- a/test/tests/test_sympl_testfield.f90 +++ b/test/tests/test_sympl_testfield.f90 @@ -1,235 +1,258 @@ module failed_symplectic_step_backend - use field_can_base, only: field_can_t - use orbit_symplectic_base, only: symplectic_integrator_t, SYMPLECTIC_STEP_BOUNDARY + use field_can_base, only: field_can_t + use orbit_symplectic_base, only: symplectic_integrator_t, SYMPLECTIC_STEP_BOUNDARY - implicit none + implicit none contains - subroutine fail_symplectic_step(si, f, ierr) - type(symplectic_integrator_t), intent(inout) :: si - type(field_can_t), intent(inout) :: f - integer, intent(out) :: ierr - - f%Bmod = 8.0d0 - f%vpar = 3.0d0 - f%mu = 1.0d0 - ierr = 1 - end subroutine fail_symplectic_step - - subroutine locate_lcfs_step(si, f, ierr) - type(symplectic_integrator_t), intent(inout) :: si - type(field_can_t), intent(inout) :: f - integer, intent(out) :: ierr - - si%z = [1.0d0, 0.2d0, 0.3d0, 0.4d0] - si%last_step_fraction = 0.25d0 - si%last_event_radial_residual = 1.0d-11 - si%last_event_fraction_width = 1.0d-12 - f%Bmod = 1.0d0 - f%mu = 1.0d0 - f%vpar = 0.0d0 - ierr = SYMPLECTIC_STEP_BOUNDARY - end subroutine locate_lcfs_step + subroutine fail_symplectic_step(si, f, ierr) + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(out) :: ierr + + f%Bmod = 8.0d0 + f%vpar = 3.0d0 + f%mu = 1.0d0 + ierr = 1 + end subroutine fail_symplectic_step + + subroutine locate_lcfs_step(si, f, ierr) + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(out) :: ierr + + si%z = [1.0d0, 0.2d0, 0.3d0, 0.4d0] + si%last_step_fraction = 0.25d0 + si%last_event_radial_residual = 1.0d-11 + si%last_event_fraction_width = 1.0d-12 + f%Bmod = 1.0d0 + f%mu = 1.0d0 + f%vpar = 0.0d0 + ierr = SYMPLECTIC_STEP_BOUNDARY + end subroutine locate_lcfs_step end module failed_symplectic_step_backend program test_sympl_testfield - use, intrinsic :: iso_fortran_env, only : dp => real64, int64 - use failed_symplectic_step_backend, only: fail_symplectic_step, locate_lcfs_step - use classification, only: classify_classifier_exit - use simple_main, only : classify_orbit_exit, init_field, locate_linear_lcfs, & - macrostep, macrostep_with_wall_check - use simple, only : tracer_t, init_sympl, ORBIT_FO_LOSS, ORBIT_FO_NUMERICAL - use params, only : isw_field_type, field_input, coord_input, integmode, & - swcoll, orbit_model, ORBIT_GC, ORBIT_FULL_ORBIT, ORBIT_EXIT_LCFS, & - ORBIT_EXIT_WALL, ORBIT_EXIT_NUMERICAL_DOMAIN, & - ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_FULL_ORBIT - use magfie_sub, only : TEST - use field_can_mod, only : evaluate - use orbit_symplectic, only : orbit_timestep_sympl - use orbit_symplectic_base, only: SYMPLECTIC_STEP_BOUNDARY, & - SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_MAXITER - - implicit none - - type(tracer_t) :: norb - character(*), parameter :: vmec_file = WOUT_FILE - integer, parameter :: ans_s = 5, ans_tp = 5, amultharm = 5 - integer :: ierr - real(dp) :: initial_si_z(4) - - ! Configure symplectic GC with TEST field to ensure initialization succeeds - isw_field_type = TEST - field_input = vmec_file - coord_input = vmec_file - integmode = 1 - - call init_field(norb, vmec_file, ans_s, ans_tp, amultharm, integmode) - - if (.not. associated(evaluate)) then - print *, 'evaluate pointer not associated for TEST field' - stop 1 - end if - - ! Smoke step: one symplectic timestep with test field - norb%relerr = 1.0e-12_dp - norb%dtaumin = 1.0e-3_dp - norb%dtau = 1.0e-3_dp - - call init_sympl(norb%si, norb%f, (/0.2_dp, 1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/), & - norb%dtau, norb%dtaumin, norb%relerr, 1) - - ierr = 0 - initial_si_z = norb%si%z - call orbit_timestep_sympl(norb%si, norb%f, ierr) - if (ierr == SYMPLECTIC_STEP_OK) then - error stop 'unconverged Euler test-field step lost its failure status' - end if - if (any(norb%si%z /= initial_si_z)) then - error stop 'unconverged Euler test-field step changed accepted position' - end if - - call test_macrostep_lcfs_event - call test_failed_step_preserves_state - call test_exit_classification - call test_fo_lcfs_location - - print *, 'TEST field symplectic step succeeded' - -contains + use, intrinsic :: iso_fortran_env, only : dp => real64, int64 + use failed_symplectic_step_backend, only: fail_symplectic_step, locate_lcfs_step + use classification, only: classify_classifier_exit + use simple_main, only : classify_orbit_exit, init_field, locate_linear_lcfs, & + macrostep, macrostep_with_wall_check + use simple, only : tracer_t, init_sympl, ORBIT_FO_LOSS, ORBIT_FO_NUMERICAL + use params, only : isw_field_type, field_input, coord_input, integmode, & + swcoll, orbit_model, ORBIT_GC, ORBIT_FULL_ORBIT, ORBIT_EXIT_LCFS, & + ORBIT_EXIT_WALL, ORBIT_EXIT_NUMERICAL_DOMAIN, & + ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_FULL_ORBIT + use magfie_sub, only : TEST + use field_can_mod, only : evaluate + use orbit_symplectic, only : orbit_timestep_sympl + use orbit_symplectic_base, only: SYMPLECTIC_STEP_BOUNDARY, & + SYMPLECTIC_STEP_OK, SYMPLECTIC_STEP_MAXITER, & + symplectic_newton_warning_mode + + implicit none + + type(tracer_t) :: norb + character(*), parameter :: vmec_file = WOUT_FILE + integer, parameter :: ans_s = 5, ans_tp = 5, amultharm = 5 + integer :: ierr + real(dp) :: initial_si_z(4) + + ! Configure symplectic GC with TEST field to ensure initialization succeeds + isw_field_type = TEST + field_input = vmec_file + coord_input = vmec_file + integmode = 1 + + call init_field(norb, vmec_file, ans_s, ans_tp, amultharm, integmode) + + if (.not. associated(evaluate)) then + print *, 'evaluate pointer not associated for TEST field' + stop 1 + end if - subroutine test_macrostep_lcfs_event - real(dp) :: z(5), exit_step - integer(int64) :: kt - integer :: step_error + ! Smoke step: one symplectic timestep with test field + norb%relerr = 1.0e-12_dp + norb%dtaumin = 1.0e-3_dp + norb%dtau = 1.0e-3_dp - swcoll = .false. - orbit_model = ORBIT_GC - orbit_timestep_sympl => locate_lcfs_step - z = [0.9_dp, 0.1_dp, 0.2_dp, 1.0_dp, 0.0_dp] - kt = 0_int64 - call macrostep(norb, z, kt, step_error, 1, exit_step) + call init_sympl(norb%si, norb%f, (/0.2_dp, 1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/), & + norb%dtau, norb%dtaumin, norb%relerr, 1) - if (step_error /= SYMPLECTIC_STEP_BOUNDARY) then - error stop 'located LCFS event lost its physical status' + ierr = 0 + initial_si_z = norb%si%z + call orbit_timestep_sympl(norb%si, norb%f, ierr) + if (ierr == SYMPLECTIC_STEP_OK) then + error stop 'unconverged Euler test-field step lost its failure status' end if - if (kt /= 0_int64) error stop 'fractional LCFS event advanced a full step' - if (exit_step /= 0.25_dp) error stop 'LCFS event time lost its step fraction' - if (z(1) /= 1.0_dp) error stop 'LCFS event endpoint was not committed' - end subroutine test_macrostep_lcfs_event - - subroutine test_failed_step_preserves_state - real(dp), parameter :: initial_state(5) = [0.2_dp, 1.0_dp, 2.0_dp, & - 1.0_dp, 0.25_dp] - real(dp) :: z(5) - real(dp) :: x_previous(3) - integer(int64) :: kt - integer :: step_error - - swcoll = .false. - orbit_model = ORBIT_GC - orbit_timestep_sympl => fail_symplectic_step - z = initial_state - kt = 0_int64 - - call macrostep(norb, z, kt, step_error, 1) - - if (step_error /= 1) error stop 'failed step status was not preserved' - if (kt /= 0_int64) error stop 'failed step advanced the time index' - if (any(z /= initial_state)) then - error stop 'failed symplectic step changed the accepted state' + if (any(norb%si%z /= initial_si_z)) then + error stop 'unconverged Euler test-field step changed accepted position' end if - z = initial_state - x_previous = 0.0_dp - kt = 0_int64 - call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous) - if (step_error /= 1) error stop 'wall path lost failed step status' - if (kt /= 0_int64) error stop 'failed wall path advanced the time index' - if (any(z /= initial_state)) then - error stop 'failed wall path changed the accepted state' - end if - end subroutine test_failed_step_preserves_state + call test_macrostep_lcfs_event + call test_failed_step_preserves_state + call test_exit_classification + call test_fo_lcfs_location - subroutine test_exit_classification - logical :: classifier_lost - integer :: classifier_exit + print *, 'TEST field symplectic step succeeded' - if (classify_orbit_exit(SYMPLECTIC_STEP_BOUNDARY, ORBIT_GC, 3, .true.) /= & - ORBIT_EXIT_LCFS) then - error stop 'converged LCFS event was not classified as physical' - end if - if (classify_orbit_exit(SYMPLECTIC_STEP_BOUNDARY, ORBIT_GC, 3, .false.) /= & - ORBIT_EXIT_NUMERICAL_DOMAIN) then - error stop 'extended map boundary was classified as physical LCFS' - end if - if (classify_orbit_exit(77, ORBIT_GC, 3, .true.) /= ORBIT_EXIT_WALL) then - error stop 'wall event was not classified as physical' - end if - if (classify_orbit_exit(77, ORBIT_GC, 0, .true.) /= ORBIT_EXIT_WALL) then - error stop 'axis wall event was not classified as physical' - end if - if (classify_orbit_exit(1, ORBIT_GC, 3, .true.) /= & - ORBIT_EXIT_NUMERICAL_DOMAIN) then - error stop 'exterior Newton iterate was classified as physical' - end if - if (classify_orbit_exit(ORBIT_FO_LOSS, ORBIT_FULL_ORBIT, 3, .true.) /= & - ORBIT_EXIT_LCFS) then - error stop 'full-orbit LCFS exit was not classified as physical' - end if - if (classify_orbit_exit(ORBIT_FO_LOSS, ORBIT_FULL_ORBIT, 3, .false.) /= & - ORBIT_EXIT_NUMERICAL_DOMAIN) then - error stop 'extended full-orbit map boundary was classified as physical' - end if - if (classify_orbit_exit(ORBIT_FO_NUMERICAL, ORBIT_FULL_ORBIT, 3, .true.) /= & - ORBIT_EXIT_NUMERICAL_FULL_ORBIT) then - error stop 'full-orbit locate failure was classified as physical' - end if - if (classify_orbit_exit(1, ORBIT_GC, 0, .true.) /= ORBIT_EXIT_LCFS) then - error stop 'RK LCFS exit was classified as numerical' - end if - if (classify_orbit_exit(1, ORBIT_GC, 0, .false.) /= & - ORBIT_EXIT_NUMERICAL_DOMAIN) then - error stop 'RK extended-map boundary was classified as physical' - end if +contains - call classify_classifier_exit(SYMPLECTIC_STEP_MAXITER, 3, & - classifier_lost, classifier_exit) - if (classifier_lost .or. classifier_exit /= ORBIT_EXIT_NUMERICAL_MAXITER) then - error stop 'classifier mode promoted Newton maxiter to a physical loss' - end if - call classify_classifier_exit(SYMPLECTIC_STEP_BOUNDARY, 3, & - classifier_lost, classifier_exit) - if (.not. classifier_lost .or. classifier_exit /= ORBIT_EXIT_LCFS) then - error stop 'classifier mode lost its physical boundary classification' - end if - end subroutine test_exit_classification - - subroutine test_fo_lcfs_location - real(dp), parameter :: z_before(5) = [0.9_dp, 6.2_dp, 3.0_dp, 1.0_dp, & - -0.5_dp] - real(dp), parameter :: z_after(5) = [1.1_dp, 0.1_dp, -3.0_dp, 2.0_dp, & - 0.5_dp] - real(dp), parameter :: tolerance = 32.0_dp*epsilon(1.0_dp) - real(dp) :: z_event(5), event_fraction - - call locate_linear_lcfs(z_before, z_after, 6.0_dp, z_event, event_fraction) - if (abs(event_fraction - 0.5_dp) > tolerance) then - error stop 'full-orbit LCFS fraction is incorrect' - end if - if (z_event(1) /= 1.0_dp) error stop 'full-orbit event is not on the LCFS' - if (abs(z_event(4) - 1.5_dp) > tolerance .or. & - abs(z_event(5)) > tolerance) then - error stop 'full-orbit phase-space state is inconsistent with event time' - end if - if (abs(z_event(2) - (6.2_dp + 0.5_dp*(0.1_dp - 6.2_dp + & - 2.0_dp*acos(-1.0_dp)))) > tolerance) then - error stop 'full-orbit poloidal seam was not interpolated periodically' - end if - if (abs(z_event(3) - 3.0_dp) > tolerance) then - error stop 'full-orbit field-period seam was not interpolated periodically' - end if - end subroutine test_fo_lcfs_location + subroutine test_macrostep_lcfs_event + real(dp) :: z(5), exit_step + integer(int64) :: kt + integer :: step_error + + swcoll = .false. + orbit_model = ORBIT_GC + orbit_timestep_sympl => locate_lcfs_step + z = [0.9_dp, 0.1_dp, 0.2_dp, 1.0_dp, 0.0_dp] + kt = 0_int64 + call macrostep(norb, z, kt, step_error, 1, exit_step) + + if (step_error /= SYMPLECTIC_STEP_BOUNDARY) then + error stop 'located LCFS event lost its physical status' + end if + if (kt /= 0_int64) error stop 'fractional LCFS event advanced a full step' + if (exit_step /= 0.25_dp) error stop 'LCFS event time lost its step fraction' + if (z(1) /= 1.0_dp) error stop 'LCFS event endpoint was not committed' + end subroutine test_macrostep_lcfs_event + + subroutine test_failed_step_preserves_state + real(dp), parameter :: initial_state(5) = [0.2_dp, 1.0_dp, 2.0_dp, & + 1.0_dp, 0.25_dp] + real(dp) :: z(5) + real(dp) :: x_previous(3) + integer(int64) :: kt + integer :: step_error + + swcoll = .false. + orbit_model = ORBIT_GC + orbit_timestep_sympl => fail_symplectic_step + z = initial_state + kt = 0_int64 + symplectic_newton_warning_mode = .false. + + call macrostep(norb, z, kt, step_error, 1) + + if (step_error /= 1) error stop 'failed step status was not preserved' + if (kt /= 0_int64) error stop 'failed step advanced the time index' + if (any(z /= initial_state)) then + error stop 'failed symplectic step changed the accepted state' + end if + + z = initial_state + x_previous = 0.0_dp + kt = 0_int64 + call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous) + if (step_error /= 1) error stop 'wall path lost failed step status' + if (kt /= 0_int64) error stop 'failed wall path advanced the time index' + if (any(z /= initial_state)) then + error stop 'failed wall path changed the accepted state' + end if + + z = initial_state + kt = 0_int64 + symplectic_newton_warning_mode = .true. + call macrostep(norb, z, kt, step_error, 1) + if (step_error /= 0) error stop 'warning mode stopped a failed step' + if (kt /= 1_int64) error stop 'warning mode did not consume the failed step' + if (any(z /= initial_state)) then + error stop 'warning-mode skip changed the last accepted state' + end if + + z = initial_state + x_previous = 0.0_dp + kt = 0_int64 + call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous) + if (step_error /= 0) error stop 'warning-mode wall path stopped a failed step' + if (kt /= 1_int64) & + error stop 'warning-mode wall path did not consume the failed step' + if (any(z /= initial_state)) then + error stop 'warning-mode wall skip changed the last accepted state' + end if + end subroutine test_failed_step_preserves_state + + subroutine test_exit_classification + logical :: classifier_lost + integer :: classifier_exit + + if (classify_orbit_exit(SYMPLECTIC_STEP_BOUNDARY, ORBIT_GC, 3, .true.) /= & + ORBIT_EXIT_LCFS) then + error stop 'converged LCFS event was not classified as physical' + end if + if (classify_orbit_exit(SYMPLECTIC_STEP_BOUNDARY, ORBIT_GC, 3, .false.) /= & + ORBIT_EXIT_NUMERICAL_DOMAIN) then + error stop 'extended map boundary was classified as physical LCFS' + end if + if (classify_orbit_exit(77, ORBIT_GC, 3, .true.) /= ORBIT_EXIT_WALL) then + error stop 'wall event was not classified as physical' + end if + if (classify_orbit_exit(77, ORBIT_GC, 0, .true.) /= ORBIT_EXIT_WALL) then + error stop 'axis wall event was not classified as physical' + end if + if (classify_orbit_exit(1, ORBIT_GC, 3, .true.) /= & + ORBIT_EXIT_NUMERICAL_DOMAIN) then + error stop 'exterior Newton iterate was classified as physical' + end if + if (classify_orbit_exit(ORBIT_FO_LOSS, ORBIT_FULL_ORBIT, 3, .true.) /= & + ORBIT_EXIT_LCFS) then + error stop 'full-orbit LCFS exit was not classified as physical' + end if + if (classify_orbit_exit(ORBIT_FO_LOSS, ORBIT_FULL_ORBIT, 3, .false.) /= & + ORBIT_EXIT_NUMERICAL_DOMAIN) then + error stop 'extended full-orbit map boundary was classified as physical' + end if + if (classify_orbit_exit(ORBIT_FO_NUMERICAL, ORBIT_FULL_ORBIT, 3, .true.) /= & + ORBIT_EXIT_NUMERICAL_FULL_ORBIT) then + error stop 'full-orbit locate failure was classified as physical' + end if + if (classify_orbit_exit(1, ORBIT_GC, 0, .true.) /= ORBIT_EXIT_LCFS) then + error stop 'RK LCFS exit was classified as numerical' + end if + if (classify_orbit_exit(1, ORBIT_GC, 0, .false.) /= & + ORBIT_EXIT_NUMERICAL_DOMAIN) then + error stop 'RK extended-map boundary was classified as physical' + end if + + call classify_classifier_exit(SYMPLECTIC_STEP_MAXITER, 3, & + classifier_lost, classifier_exit) + if (classifier_lost .or. classifier_exit /= ORBIT_EXIT_NUMERICAL_MAXITER) then + error stop 'classifier mode promoted Newton maxiter to a physical loss' + end if + call classify_classifier_exit(SYMPLECTIC_STEP_BOUNDARY, 3, & + classifier_lost, classifier_exit) + if (.not. classifier_lost .or. classifier_exit /= ORBIT_EXIT_LCFS) then + error stop 'classifier mode lost its physical boundary classification' + end if + end subroutine test_exit_classification + + subroutine test_fo_lcfs_location + real(dp), parameter :: z_before(5) = [0.9_dp, 6.2_dp, 3.0_dp, 1.0_dp, & + -0.5_dp] + real(dp), parameter :: z_after(5) = [1.1_dp, 0.1_dp, -3.0_dp, 2.0_dp, & + 0.5_dp] + real(dp), parameter :: tolerance = 32.0_dp*epsilon(1.0_dp) + real(dp) :: z_event(5), event_fraction + + call locate_linear_lcfs(z_before, z_after, 6.0_dp, z_event, event_fraction) + if (abs(event_fraction - 0.5_dp) > tolerance) then + error stop 'full-orbit LCFS fraction is incorrect' + end if + if (z_event(1) /= 1.0_dp) error stop 'full-orbit event is not on the LCFS' + if (abs(z_event(4) - 1.5_dp) > tolerance .or. & + abs(z_event(5)) > tolerance) then + error stop 'full-orbit phase-space state is inconsistent with event time' + end if + if (abs(z_event(2) - (6.2_dp + 0.5_dp*(0.1_dp - 6.2_dp + & + 2.0_dp*acos(-1.0_dp)))) > tolerance) then + error stop 'full-orbit poloidal seam was not interpolated periodically' + end if + if (abs(z_event(3) - 3.0_dp) > tolerance) then + error stop 'full-orbit field-period seam was not interpolated periodically' + end if + end subroutine test_fo_lcfs_location end program test_sympl_testfield From cb0a4f96fa9a202c85ad389aff4db8782a69330f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 15 Jul 2026 22:45:01 +0200 Subject: [PATCH 02/11] Report confinement over resolved orbits --- README.md | 14 ++- python/pysimple/plotting.py | 4 +- src/restart.f90 | 13 +-- src/simple_main.f90 | 37 +++++-- test/tests/CMakeLists.txt | 9 ++ test/tests/test_result_accounting.f90 | 135 ++++++++++++++++++++++++++ test/tests/test_spectre_validation.py | 42 ++++---- 7 files changed, 215 insertions(+), 39 deletions(-) create mode 100644 test/tests/test_result_accounting.f90 diff --git a/README.md b/README.md index 38acdd94..3f79c6e3 100644 --- a/README.md +++ b/README.md @@ -239,15 +239,23 @@ Diagnostics for slow convergence of Newton iterations are written in `fort.6601` 1. Physical time 2. Confined fraction of passing particles 3. Confined fraction of trapped particles -4. Total number of particles +4. Number of numerically resolved particles used as the denominator -The sum of 2. and 3. yields the overall confined fraction at each time. +The sum of 2. and 3. yields the overall confined fraction among resolved +particles at each time. Numerically fatal orbits are excluded from both the +numerator and denominator. If no orbit is resolved, both fractions are `NaN`. +`unresolved_fraction.dat` reports the excluded numerical-failure fraction +against the total initial population; its third column remains the total +number of particles. `times_lost.dat` contains the loss time of each particle. Columns are: 1. Particle index. Corresponds to line number in start.dat . -2. Time t_loss [s] when the particle is lost. Possible values are: -1, `trace_time`, or any other value between 0 and trace_time. +2. Time t_loss [s] when the particle is lost. Possible values are: -1, `NaN`, + `trace_time`, or any other value between 0 and trace_time. * If never lost or classified as regular, maximum tracing time `trace_time` is written. * If ignored due to contr_pp, which defines deep passing region as confined (we don consider them anymore), -1 is written. + * If tracing ends in a numerical fatal condition, `NaN` is written. Such an + orbit is unresolved, not physically lost or confined. 3. Trapping parameter trap_par that is 1 for deeply trapped, 0 for tp boundary and negative for passing. Eq. (3.1) in Accelerated Methods paper. Whenever trap_par < contr_pp, particle is not traced and counted as confined. diff --git a/python/pysimple/plotting.py b/python/pysimple/plotting.py index 67c20462..63a7c71b 100644 --- a/python/pysimple/plotting.py +++ b/python/pysimple/plotting.py @@ -44,7 +44,7 @@ Column 1: Time [s] Column 2: confpart_pass (confined passing fraction) Column 3: confpart_trap (confined trapped fraction) - Column 4: Total number of particles + Column 4: Numerically resolved particles used as the denominator Note: Total confined fraction = Column 2 + Column 3 Example @@ -210,7 +210,7 @@ def load_loss_data(directory: str | Path) -> LossData: start_pitch = np.zeros(n_particles) # Load confined_fraction.dat - # Columns: time, confpart_pass, confpart_trap, ntestpart + # Columns: time, confpart_pass, confpart_trap, resolved particle count conf_file = directory / "confined_fraction.dat" if conf_file.exists(): conf_data = np.loadtxt(conf_file) diff --git a/src/restart.f90 b/src/restart.f90 index 1f453207..6f2d52f6 100644 --- a/src/restart.f90 +++ b/src/restart.f90 @@ -2,11 +2,11 @@ module restart_mod use, intrinsic :: iso_fortran_env, only: dp => real64 use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use params, only: ntestpart, times_lost, orbit_exit_code, & - boundary_event_radial_residual, boundary_event_time_width, & - trap_par, perp_inv, zend, & - confpart_pass, confpart_trap, ntimstep, kt_macro, & - v0, dtaumin, trace_time, ORBIT_EXIT_COMPLETED, & - ORBIT_EXIT_LCFS + boundary_event_radial_residual, boundary_event_time_width, & + trap_par, perp_inv, zend, & + confpart_pass, confpart_trap, ntimstep, kt_macro, & + v0, dtaumin, trace_time, ORBIT_EXIT_COMPLETED, & + ORBIT_EXIT_LCFS, ORBIT_EXIT_WALL implicit none private @@ -64,7 +64,8 @@ subroutine read_restart_data() boundary_event_radial_residual(idx) = radial_residual boundary_event_time_width(idx) = time_width particle_done(idx) = exit_code == ORBIT_EXIT_COMPLETED .or. & - exit_code == ORBIT_EXIT_LCFS + exit_code == ORBIT_EXIT_LCFS .or. & + exit_code == ORBIT_EXIT_WALL end do close (unit) end if diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 6fe69d9d..b13396d8 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -1262,6 +1262,12 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & end do case (SYMPL_SPECTRE_SKIM) orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED + case default + orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT + do it_f = it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do end select exit end if @@ -1300,7 +1306,7 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & ! Mirror-confined at an interior interface: cannot be lost, so record ! it as confined (times_lost = trace_time) rather than at its stop. t_stop = trace_time - else if (ierr_orbit == SYMPL_SPECTRE_STOP) then + else if (ierr_orbit /= SYMPL_SPECTRE_LOSS) then t_stop = ieee_value(0.0_dp, ieee_quiet_nan) else t_stop = (real(kt, dp) + t_frac)*dtaumin/v0 @@ -1701,12 +1707,14 @@ subroutine write_results !> Write the per-particle and confined-fraction result files from the !> shared result arrays. progress_monitor calls this every !> checkpoint_interval seconds, so a run killed mid-flight keeps its - !> last flushed output. Confined fractions are normalised by ntestpart, - !> as in the final write, so a partial file is a converging lower bound - !> and never exceeds one; particles not yet finished keep their sentinel - !> values (times_lost = -1). - integer :: i, num_lost, unit - real(dp) :: inverse_times_lost_sum, norm + !> last flushed output. Confined fractions are conditional on the + !> numerically resolved population at each time. A partial file remains + !> a lower bound because unfinished particles have not contributed to + !> the numerator. Particles not yet finished keep times_lost = -1. + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + + integer :: i, num_lost, resolved_count, unit + real(dp) :: inverse_times_lost_sum, norm, pass_fraction, trap_fraction norm = real(max(ntestpart, 1), dp) @@ -1719,6 +1727,9 @@ subroutine write_results num_lost = 0 inverse_times_lost_sum = 0.0d0 do i = 1, ntestpart + if (orbit_exit_code(i) >= ORBIT_EXIT_NUMERICAL_DOMAIN) then + times_lost(i) = ieee_value(0.0_dp, ieee_quiet_nan) + end if write (unit, *) i, times_lost(i), trap_par(i), zstart(1, i), & perp_inv(i), zend(:, i) if (orbit_exit_code(i) == ORBIT_EXIT_LCFS .or. & @@ -1752,8 +1763,16 @@ subroutine write_results open (newunit=unit, file='confined_fraction.dat', recl=1024) do i = 1, ntimstep - write (unit, *) dble(kt_macro(i))*dtaumin/v0, confpart_pass(i)/norm, & - confpart_trap(i)/norm, ntestpart + resolved_count = ntestpart - unresolved_orbits(i) + if (resolved_count > 0) then + pass_fraction = confpart_pass(i)/real(resolved_count, dp) + trap_fraction = confpart_trap(i)/real(resolved_count, dp) + else + pass_fraction = ieee_value(0.0_dp, ieee_quiet_nan) + trap_fraction = ieee_value(0.0_dp, ieee_quiet_nan) + end if + write (unit, *) dble(kt_macro(i))*dtaumin/v0, pass_fraction, & + trap_fraction, resolved_count end do close (unit) diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index df647159..85023095 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -821,6 +821,15 @@ set_tests_properties(test_restart_io PROPERTIES LABELS "unit" TIMEOUT 15) +add_executable(test_result_accounting.x test_result_accounting.f90) +target_link_libraries(test_result_accounting.x simple) +add_test(NAME test_result_accounting + COMMAND test_result_accounting.x + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_tests_properties(test_result_accounting PROPERTIES + LABELS "unit" + TIMEOUT 15) + # Restart integration test (TEST field, deterministic) add_test(NAME test_restart_integration COMMAND ${Python3_EXECUTABLE} diff --git a/test/tests/test_result_accounting.f90 b/test/tests/test_result_accounting.f90 new file mode 100644 index 00000000..66c86459 --- /dev/null +++ b/test/tests/test_result_accounting.f90 @@ -0,0 +1,135 @@ +program test_result_accounting + use, intrinsic :: ieee_arithmetic, only: ieee_is_nan + use, intrinsic :: iso_fortran_env, only: dp => real64, int64 + use magfie_sub, only: TEST + use params, only: ntestpart, ntimstep, kt_macro, dtaumin, v0, & + confpart_pass, confpart_trap, unresolved_orbits, & + times_lost, orbit_exit_code, trap_par, perp_inv, & + zstart, zend, boundary_event_radial_residual, & + boundary_event_time_width, isw_field_type, class_plot, & + ntcut, ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, & + ORBIT_EXIT_NUMERICAL_MAXITER + use simple_main, only: write_results + implicit none + + integer :: nerr + + nerr = 0 + call setup_results() + call test_resolved_denominator() + call test_zero_resolved() + call teardown_results() + + if (nerr > 0) error stop 1 + print *, 'All result accounting tests passed!' + +contains + + subroutine setup_results() + ntestpart = 3 + ntimstep = 1 + dtaumin = 1.0_dp + v0 = 1.0_dp + isw_field_type = TEST + class_plot = .false. + ntcut = 0 + + allocate (kt_macro(1), confpart_pass(1), confpart_trap(1)) + allocate (unresolved_orbits(1), times_lost(3), orbit_exit_code(3)) + allocate (trap_par(3), perp_inv(3), zstart(5, 3), zend(5, 3)) + allocate (boundary_event_radial_residual(3)) + allocate (boundary_event_time_width(3)) + + kt_macro = [10_int64] + trap_par = 0.0_dp + perp_inv = 0.0_dp + zstart = 0.0_dp + zend = 0.0_dp + boundary_event_radial_residual = -1.0_dp + boundary_event_time_width = -1.0_dp + end subroutine setup_results + + subroutine test_resolved_denominator() + integer :: idx, unit + real(dp) :: time, pass_fraction, trap_fraction, loss_time + real(dp) :: unresolved_fraction + integer :: resolved_count, total_count + + confpart_pass = 1.0_dp + confpart_trap = 0.0_dp + unresolved_orbits = 1 + times_lost = [10.0_dp, 5.0_dp, 5.0_dp] + orbit_exit_code = [ORBIT_EXIT_COMPLETED, ORBIT_EXIT_LCFS, & + ORBIT_EXIT_NUMERICAL_MAXITER] + + call write_results() + + open (newunit=unit, file='confined_fraction.dat', status='old') + read (unit, *) time, pass_fraction, trap_fraction, resolved_count + close (unit) + call assert_close(pass_fraction, 0.5_dp, 'resolved confined fraction') + call assert_close(trap_fraction, 0.0_dp, 'resolved trapped fraction') + call assert_true(resolved_count == 2, 'resolved denominator is two') + + open (newunit=unit, file='times_lost.dat', status='old') + read (unit, *) idx, loss_time + read (unit, *) idx, loss_time + read (unit, *) idx, loss_time + close (unit) + call assert_true(idx == 3, 'read numerical particle row') + call assert_true(ieee_is_nan(loss_time), 'numerical loss time is NaN') + + open (newunit=unit, file='unresolved_fraction.dat', status='old') + read (unit, *) time, unresolved_fraction, total_count + close (unit) + call assert_close(unresolved_fraction, 1.0_dp/3.0_dp, & + 'unresolved fraction keeps total denominator') + call assert_true(total_count == 3, 'unresolved denominator is total') + end subroutine test_resolved_denominator + + subroutine test_zero_resolved() + integer :: unit, resolved_count + real(dp) :: time, pass_fraction, trap_fraction + + confpart_pass = 0.0_dp + confpart_trap = 0.0_dp + unresolved_orbits = ntestpart + + call write_results() + + open (newunit=unit, file='confined_fraction.dat', status='old') + read (unit, *) time, pass_fraction, trap_fraction, resolved_count + close (unit) + call assert_true(resolved_count == 0, 'zero resolved denominator') + call assert_true(ieee_is_nan(pass_fraction), 'zero denominator pass NaN') + call assert_true(ieee_is_nan(trap_fraction), 'zero denominator trap NaN') + end subroutine test_zero_resolved + + subroutine teardown_results() + deallocate (kt_macro, confpart_pass, confpart_trap, unresolved_orbits) + deallocate (times_lost, orbit_exit_code, trap_par, perp_inv) + deallocate (zstart, zend, boundary_event_radial_residual) + deallocate (boundary_event_time_width) + end subroutine teardown_results + + subroutine assert_close(actual, expected, message) + real(dp), intent(in) :: actual, expected + character(*), intent(in) :: message + + if (abs(actual - expected) > 1.0e-12_dp) then + print *, 'FAIL:', message, 'expected', expected, 'got', actual + nerr = nerr + 1 + end if + end subroutine assert_close + + subroutine assert_true(condition, message) + logical, intent(in) :: condition + character(*), intent(in) :: message + + if (.not. condition) then + print *, 'FAIL:', message + nerr = nerr + 1 + end if + end subroutine assert_true + +end program test_result_accounting diff --git a/test/tests/test_spectre_validation.py b/test/tests/test_spectre_validation.py index 8669058c..244b68b5 100644 --- a/test/tests/test_spectre_validation.py +++ b/test/tests/test_spectre_validation.py @@ -91,9 +91,11 @@ def run(binary, workdir, h5, **kwargs): def split_confined(times_lost, trace_time): - confined = times_lost >= trace_time * (1.0 - 1e-9) - terminated = (times_lost > 0.0) & (times_lost < trace_time * (1.0 - 1e-9)) - unresolved = times_lost <= 0.0 + finite = np.isfinite(times_lost) + confined = finite & (times_lost >= trace_time * (1.0 - 1e-9)) + terminated = (finite & (times_lost > 0.0) + & (times_lost < trace_time * (1.0 - 1e-9))) + unresolved = ~finite return confined, terminated, unresolved @@ -104,40 +106,42 @@ def losses_and_accounting(binary, h5, failures): run(binary, work, h5, crossing_level=1) tl = np.loadtxt(os.path.join(work, "times_lost.dat")) cf = np.loadtxt(os.path.join(work, "confined_fraction.dat")) + exit_codes = np.loadtxt(os.path.join(work, "orbit_exit_code.dat")) ev_path = os.path.join(work, "spectre_crossing_events.dat") crossings_written = os.path.exists(ev_path) and os.path.getsize(ev_path) > 0 if len(tl) != NPART: failures.append(f"accounting: times_lost.dat has {len(tl)} rows != {NPART}") - if not np.all(np.isfinite(tl[:, 1])): - failures.append("accounting: non-finite loss time (unresolved marker)") - confined, terminated, unresolved = split_confined(tl[:, 1], TRACE_TIME) n_conf, n_term, n_unres = (int(confined.sum()), int(terminated.sum()), int(unresolved.sum())) - if n_unres != 0: - failures.append(f"accounting: {n_unres} markers left unresolved " - f"(times_lost <= 0)") - if n_conf + n_term != NPART: - failures.append(f"accounting: {n_conf} confined + {n_term} terminated " + if n_conf + n_term + n_unres != NPART: + failures.append(f"accounting: {n_conf} confined + {n_term} terminated + " + f"{n_unres} unresolved " f"!= {NPART} markers") + if not np.all(exit_codes[unresolved, 1] >= 101): + failures.append("accounting: unresolved marker lacks numerical exit code") + if np.any(exit_codes[confined | terminated, 1] >= 101): + failures.append("accounting: resolved marker has numerical exit code") # confined_fraction.dat is filled by an independent per-timestep confined # counter; its final value must equal the times_lost confined count. Breaking # the loss tally (counting a lost marker as confined) makes these disagree. + n_resolved = n_conf + n_term frac_final = cf[-1, 1] + cf[-1, 2] - if abs(frac_final - n_conf / NPART) > 0.5 / NPART + 1e-9: + if abs(frac_final - n_conf / n_resolved) > 0.5 / n_resolved + 1e-9: failures.append(f"accounting: confined_fraction {frac_final:.4f} != " - f"times_lost confined {n_conf / NPART:.4f}") - if int(cf[-1, 3]) != NPART: + f"resolved times_lost confined " + f"{n_conf / n_resolved:.4f}") + if int(cf[-1, 3]) != n_resolved: failures.append(f"accounting: confined_fraction N = {int(cf[-1, 3])}") if not crossings_written: failures.append("accounting: no spectre_crossing_events.dat produced") print(f"accounting: confined={n_conf} terminated={n_term} unresolved=" - f"{n_unres} (sum={n_conf + n_term}={NPART}) " - f"cf_final={frac_final:.3f}=({n_conf}/{NPART})") - return n_term / NPART + f"{n_unres} (sum={n_conf + n_term + n_unres}={NPART}) " + f"cf_final={frac_final:.3f}=({n_conf}/{n_resolved})") + return n_term / n_resolved def crossing_map_accounting(binary, h5, p1, failures): @@ -148,7 +152,7 @@ def crossing_map_accounting(binary, h5, p1, failures): confined0, terminated0, unresolved0 = split_confined(tl0[:, 1], TRACE_TIME) n_conf0, n_term0, n_unres0 = (int(confined0.sum()), int(terminated0.sum()), int(unresolved0.sum())) - if n_conf0 + n_term0 != NPART or n_unres0 != 0: + if n_conf0 + n_term0 + n_unres0 != NPART: failures.append(f"crossing maps: Level-0 account is {n_conf0} confined + " f"{n_term0} terminated + {n_unres0} unresolved") @@ -159,7 +163,7 @@ def crossing_map_accounting(binary, h5, p1, failures): failures.append("crossing maps: Level-0 times_lost.dat is not " "bit-identical under a fixed seed") - p0 = n_term0 / NPART + p0 = n_term0 / (n_conf0 + n_term0) print(f"crossing maps: loss_fraction L1={p1:.3f} L0={p0:.3f} " f"Level-0_account={n_conf0}+{n_term0}+{n_unres0}={NPART} " f"Level-0_repeat={'bit-identical' if raw0 == raw0_repeat else 'DIFFERS'}") From 1cb298c74595f25025123cc64f4a8624aa67d4aa Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 00:04:45 +0200 Subject: [PATCH 03/11] Bound numerical orbit recovery --- DOC/config.md | 17 +- src/classification.f90 | 45 +++- src/interface_crossing.f90 | 6 +- src/simple_main.f90 | 167 +++++++++++--- src/spectre_fo_hybrid.f90 | 18 +- src/spectre_orbit.f90 | 46 +++- src/spectre_sympl_orbit.f90 | 210 ++++++++++++------ src/sub_alpha_lifetime_can.f90 | 46 +++- .../field_can/test_spectre_sympl_crossing.f90 | 5 +- test/tests/test_sympl_testfield.f90 | 38 +++- 10 files changed, 479 insertions(+), 119 deletions(-) diff --git a/DOC/config.md b/DOC/config.md index e3734da1..9229f21f 100644 --- a/DOC/config.md +++ b/DOC/config.md @@ -15,13 +15,16 @@ reaches its iteration limit, SIMPLE continues only if the final Newton correction is finite and no more than ten times the requested relative tolerance. Each occurrence is reported by the corresponding `*_maxit` - diagnostic. If bounded retries still cannot resolve a numerical microstep, - SIMPLE retains the last accepted state for that interval, records - `warning_step_skip`, and resumes the same marker at the next microstep. - Numerical failures therefore never become physical losses or terminate a - marker in the default mode. Set the option to `.false.` for strict diagnostic - runs that roll back and stop the affected marker; those stops use a 101--105 - `orbit_exit_code` and `NaN` in `times_lost` and remain distinct from losses. + diagnostic. The production RK, symplectic, and full-orbit paths then use the + same terminal convention: if bounded recovery cannot resolve a numerical + microstep, SIMPLE retains the last accepted state for that interval, records + `warning_step_skip`, and marks the interval unresolved. Because an unchanged + state would reproduce the same deterministic failure, the next microstep uses + a circuit breaker instead of repeating the complete retry cascade. It ends + only that marker with a 101--105 `orbit_exit_code` and `NaN` in `times_lost`; + the ensemble continues, and the marker is neither confined nor physically + lost. Set the option to `.false.` for strict diagnostic runs that expose the + first exhausted recovery without the one-interval warning hold. * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi` control the Meiss or Albert canonical-map grid. Their defaults are 62, 63, diff --git a/src/classification.f90 b/src/classification.f90 index 63c2e119..d901d9c7 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -77,9 +77,9 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) real(dp), dimension(5) :: z real(dp) :: bmod,sqrtg real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl - integer :: it, ktau, it_f + integer :: first_unresolved_it, hold_streak, it, ktau, it_f integer(8) :: kt - logical :: passing + logical :: had_numerical_hold, passing integer :: ifp_tip,ifp_per integer, dimension(:), allocatable :: ipoi @@ -215,6 +215,9 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) kt = 0 + first_unresolved_it = 0 + hold_streak = 0 + had_numerical_hold = .false. if (passing) then !$omp atomic confpart_pass(1)=confpart_pass(1)+1.d0 @@ -258,16 +261,36 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) cycle endif do ktau=1,ntau + if (hold_streak /= 0) then + ierr = hold_streak + exit + end if if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr) + if (ierr /= 0 .and. ierr /= 1 .and. & + symplectic_newton_warning_mode) then + if (hold_streak == 0) then + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak = ierr + had_numerical_hold = .true. + if (first_unresolved_it == 0) first_unresolved_it = it + ierr = 0 + end if + end if else call advance_symplectic_with_retry(anorb%si, anorb%f, & orbit_timestep_sympl, ierr) if (ierr /= 0 .and. ierr /= SYMPLECTIC_STEP_BOUNDARY .and. & symplectic_newton_warning_mode) then + if (hold_streak == 0) then call count_event(EVT_WARNING_STEP_SKIP) + hold_streak = ierr + had_numerical_hold = .true. + if (first_unresolved_it == 0) first_unresolved_it = it ierr = 0 + end if else if (ierr == 0) then + hold_streak = 0 z(1:3) = anorb%si%z(1:3) z(4) = dsqrt(anorb%f%mu*anorb%f%Bmod + & 0.5d0*anorb%f%vpar**2) @@ -454,6 +477,7 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) ! write(999, *) kt*dtaumin/v0, z enddo if(ierr.ne.0) exit + if (.not. had_numerical_hold) then if(passing) then !$omp atomic confpart_pass(it)=confpart_pass(it)+1.d0 @@ -461,8 +485,15 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) !$omp atomic confpart_trap(it)=confpart_trap(it)+1.d0 endif + end if enddo + if (had_numerical_hold) then + class_result%lost = .false. + if (class_result%exit_code < ORBIT_EXIT_NUMERICAL_DOMAIN) & + class_result%exit_code = ORBIT_EXIT_NUMERICAL_EVENT + end if + !$omp critical zend(:,ipart) = z if(isw_field_type .eq. CANFLUX) then @@ -479,7 +510,8 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) deallocate(zpoipl_tip, zpoipl_per) !$omp end critical if (class_result%exit_code >= ORBIT_EXIT_NUMERICAL_DOMAIN) then - do it_f = it, ntimstep + if (first_unresolved_it == 0) first_unresolved_it = min(it, ntimstep) + do it_f = first_unresolved_it, ntimstep !$omp atomic update unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 end do @@ -494,11 +526,16 @@ pure subroutine classify_classifier_exit(ierr, mode, lost, exit_code) integer, intent(out) :: exit_code lost = .false. - if (mode <= 0 .or. ierr == SYMPLECTIC_STEP_BOUNDARY) then + if ((mode <= 0 .and. ierr == 1) .or. & + ierr == SYMPLECTIC_STEP_BOUNDARY) then lost = .true. exit_code = ORBIT_EXIT_LCFS return end if + if (mode <= 0) then + exit_code = ORBIT_EXIT_NUMERICAL_EVENT + return + end if select case (ierr) case (SYMPLECTIC_STEP_OUTSIDE_DOMAIN) diff --git a/src/interface_crossing.f90 b/src/interface_crossing.f90 index 695a8cef..00ab3ed2 100644 --- a/src/interface_crossing.f90 +++ b/src/interface_crossing.f90 @@ -25,7 +25,7 @@ module interface_crossing public :: crossing_log_reset, crossing_log_record, crossing_log_write, & crossing_log_count, crossing_log_count_type public :: CROSSING_LEVEL0, CROSSING_LEVEL1, CROSS_CROSSING, CROSS_REFLECTION, & - CROSS_LOSS, CROSS_STOP, CROSS_SHEET + CROSS_LOSS, CROSS_STOP, CROSS_SHEET, CROSS_RECOVERY integer, parameter :: CROSSING_LEVEL0 = 0 integer, parameter :: CROSSING_LEVEL1 = 1 @@ -37,6 +37,10 @@ module interface_crossing !> step itself fails to converge; regular boundaries cross via apply_crossing. integer, parameter :: CROSS_STOP = 4 integer, parameter :: CROSS_SHEET = 5 + !> Marker-local GC -> full-orbit fallback away from any physical interface. + !> iface and both volume fields are zero so post-processing cannot confuse + !> the recovery map with a discontinuity crossing. + integer, parameter :: CROSS_RECOVERY = 6 !> Inner cutoff for the innermost volume: rho_g = 0 is the coordinate axis !> where sqrt(g) = 0, so the marker reflects trivially before reaching it diff --git a/src/simple_main.f90 b/src/simple_main.f90 index b13396d8..90548f98 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -546,7 +546,10 @@ subroutine trace_parallel(norb) !$omp parallel firstprivate(norb) private(traj, times, i) allocate (traj(5, ntimstep), times(ntimstep)) - !$omp do + ! Independent markers have highly nonuniform trace costs. One-marker + ! dynamic chunks keep unfinished work available to every thread instead + ! of queueing it behind a pathological marker in a static block. + !$omp do schedule(dynamic, 1) do i = 1, ntestpart if (allocated(particle_done)) then if (particle_done(i)) then @@ -932,8 +935,10 @@ pure integer function classify_orbit_exit(ierr_orbit, model, mode, & else if (mode <= 0) then if (ierr_orbit == 1 .and. boundary_is_lcfs) then classify_orbit_exit = ORBIT_EXIT_LCFS - else + else if (ierr_orbit == 1) then classify_orbit_exit = ORBIT_EXIT_NUMERICAL_DOMAIN + else + classify_orbit_exit = ORBIT_EXIT_NUMERICAL_EVENT end if else select case (ierr_orbit) @@ -974,13 +979,16 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) real(dp), dimension(5) :: z real(dp) :: u_ref_prev(3), x_prev(3), x_prev_m(3), exit_step integer :: it, ierr_orbit, it_final, it_f + integer :: first_unresolved_it, hold_streak integer(8) :: kt - logical :: passing, faulted, physical_exit + logical :: passing, faulted, numerical_hold, physical_exit type(classification_result_t) :: class_result ierr_orbit = 0 faulted = .false. physical_exit = .false. + first_unresolved_it = 0 + hold_streak = 0 exit_step = -1d0 orbit_traj = ieee_value(0.0d0, ieee_quiet_nan) orbit_times = ieee_value(0.0d0, ieee_quiet_nan) @@ -1050,11 +1058,16 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) if (wall_enabled) then ! Use microstep-level wall checking for STL walls call macrostep_with_wall_check(anorb, z, kt, ierr_orbit, & - ntau_macro(it), ipart, x_prev_m, exit_step) + ntau_macro(it), ipart, x_prev_m, exit_step, & + hold_streak=hold_streak, & + numerical_hold_any=numerical_hold) else call macrostep(anorb, z, kt, ierr_orbit, ntau_macro(it), & - exit_step) + exit_step, hold_streak=hold_streak, & + numerical_hold_any=numerical_hold) end if + if (numerical_hold .and. first_unresolved_it == 0) & + first_unresolved_it = it end if if (ierr_orbit .ne. 0) then @@ -1075,10 +1088,7 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) anorb%si%last_event_fraction_width*dtaumin/v0 end if else - do it_f = it, ntimstep - !$omp atomic update - unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 - end do + if (first_unresolved_it == 0) first_unresolved_it = it faulted = .true. end if exit @@ -1088,10 +1098,21 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) orbit_traj(:, it) = z orbit_times(it) = kt*dtaumin/v0 - call increase_confined_count(it, passing) + if (first_unresolved_it == 0) call increase_confined_count(it, passing) it_final = it end do + if (first_unresolved_it > 0) then + do it_f = first_unresolved_it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do + physical_exit = .false. + faulted = .true. + if (orbit_exit_code(ipart) < ORBIT_EXIT_NUMERICAL_DOMAIN) & + orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT + end if + ! Fill remaining timesteps with NaN if particle left domain early if (it_final < ntimstep) then do it = it_final + 1, ntimstep @@ -1225,6 +1246,7 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & type(sympl_spectre_state_t) :: state integer :: it, it_f, ktau, ierr_orbit, it_final + integer :: first_unresolved_it, holds_before integer(8) :: kt real(dp) :: t_stop, t_frac @@ -1235,15 +1257,19 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & it_final = 0 ierr_orbit = SYMPL_SPECTRE_OK t_frac = 1.0d0 + first_unresolved_it = 0 do it = 1, ntimstep if (it >= 2) then do ktau = 1, ntau_macro(it) + holds_before = state%warning_hold_count call orbit_microstep_sympl_spectre(state, anorb%si, anorb%f, & ipart, & real(kt, dp)*dtaumin/v0, & dtaumin/v0, ierr_orbit, & t_frac) + if (state%warning_hold_count > holds_before .and. & + first_unresolved_it == 0) first_unresolved_it = it if (ierr_orbit /= SYMPL_SPECTRE_OK) exit kt = kt + 1 end do @@ -1256,19 +1282,14 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & orbit_exit_code(ipart) = ORBIT_EXIT_LCFS case (SYMPL_SPECTRE_STOP) orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_FULL_ORBIT - do it_f = it, ntimstep - !$omp atomic update - unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 - end do case (SYMPL_SPECTRE_SKIM) orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED case default orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT - do it_f = it, ntimstep - !$omp atomic update - unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 - end do end select + if (ierr_orbit /= SYMPL_SPECTRE_LOSS .and. & + ierr_orbit /= SYMPL_SPECTRE_SKIM .and. & + first_unresolved_it == 0) first_unresolved_it = it exit end if @@ -1279,10 +1300,17 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & end if orbit_traj(:, it) = z orbit_times(it) = kt*dtaumin/v0 - call increase_confined_count(it, passing) + if (first_unresolved_it == 0) call increase_confined_count(it, passing) it_final = it end do + if (first_unresolved_it > 0) then + do it_f = first_unresolved_it, ntimstep + !$omp atomic update + unresolved_orbits(it_f) = unresolved_orbits(it_f) + 1 + end do + end if + ! A mirror-confined (skimming) marker is confined for the rest of the ! trace, so it must count as confined at every remaining step for the ! confined_fraction series to stay consistent with times_lost. @@ -1299,7 +1327,10 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & end do end if - if (ierr_orbit == SYMPL_SPECTRE_OK) then + if (state%warning_hold_count > 0) then + t_stop = ieee_value(0.0_dp, ieee_quiet_nan) + orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT + else if (ierr_orbit == SYMPL_SPECTRE_OK) then t_stop = real(kt, dp)*dtaumin/v0 orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED else if (ierr_orbit == SYMPL_SPECTRE_SKIM) then @@ -1349,7 +1380,8 @@ pure subroutine locate_linear_lcfs(z_start, z_end, field_period, z_event, & event_fraction*(z_end(4:5) - z_start(4:5)) end subroutine locate_linear_lcfs - subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) + subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & + hold_streak, numerical_hold_any) use alpha_lifetime_sub, only: orbit_timestep_axis use orbit_symplectic, only: advance_symplectic_with_retry, & orbit_timestep_sympl @@ -1360,14 +1392,25 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) integer, intent(out) :: ierr_orbit integer, intent(in) :: ntau_local real(dp), intent(out), optional :: exit_step + integer, intent(inout), optional :: hold_streak + logical, intent(out), optional :: numerical_hold_any - integer :: ktau + integer :: hold_streak_local, ktau real(dp) :: z_step_start(5), z_step_end(5), loss_fraction + logical :: numerical_hold, numerical_hold_any_local if (present(exit_step)) exit_step = real(kt, dp) + hold_streak_local = 0 + if (present(hold_streak)) hold_streak_local = hold_streak + numerical_hold_any_local = .false. do ktau = 1, ntau_local + numerical_hold = .false. z_step_start = z + if (hold_streak_local /= 0) then + ierr_orbit = hold_streak_local + exit + end if if (orbit_model == ORBIT_FULL_ORBIT) then call orbit_timestep_fo(anorb%fo, z, ierr_orbit) if (ierr_orbit .ne. 0) then @@ -1384,12 +1427,22 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) ! Default production mode advances the clock past this ! one unresolved microstep and retries the next step; ! numerical inversion faults are never physical losses. - call count_event(EVT_WARNING_STEP_SKIP) z = z_step_start + if (hold_streak_local == 0) then + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. ierr_orbit = 0 end if + end if if (ierr_orbit .ne. 0) exit end if + if (numerical_hold) then + kt = kt + 1 + cycle + end if + hold_streak_local = 0 kt = kt + 1 cycle end if @@ -1402,6 +1455,18 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) if (present(exit_step)) then exit_step = real(kt, dp) + loss_fraction end if + else if (ierr_orbit /= 0 .and. & + symplectic_newton_warning_mode) then + ! The RK driver rolls back failed integration to the last + ! accepted state. Apply the same bounded warning convention + ! as the symplectic and full-orbit paths. + if (hold_streak_local == 0) then + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. + ierr_orbit = 0 + end if end if else if (swcoll) call update_momentum(anorb, z) @@ -1418,17 +1483,31 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step) ! advance_symplectic_with_retry restored the last accepted ! state. Skip only the unresolved microstep and continue; ! strict mode still reports the numerical exit. + if (hold_streak_local == 0) then call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. ierr_orbit = 0 + end if else if (ierr_orbit == 0) then call to_standard_z_coordinates(anorb, z) + hold_streak_local = 0 end if if (ierr_orbit .ne. 0) exit end if if (ierr_orbit .ne. 0) exit + if (numerical_hold) then + kt = kt + 1 + cycle + end if + if (integmode <= 0) hold_streak_local = 0 if (swcoll) call collide(z, dtaumin) ! Collisions kt = kt + 1 end do + if (present(hold_streak)) hold_streak = hold_streak_local + if (present(numerical_hold_any)) & + numerical_hold_any = numerical_hold_any_local end subroutine macrostep subroutine locate_wall_segment(z, z_start, z_end, ipart, x_start_m, & @@ -1494,7 +1573,7 @@ subroutine locate_wall_segment(z, z_start, z_end, ipart, x_start_m, & end subroutine locate_wall_segment subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & - ipart, x_prev_m, exit_step) + ipart, x_prev_m, exit_step, hold_streak, numerical_hold_any) use alpha_lifetime_sub, only: orbit_timestep_axis use orbit_symplectic, only: advance_symplectic_with_retry, & orbit_timestep_sympl @@ -1507,19 +1586,28 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & integer, intent(in) :: ipart real(dp), intent(inout) :: x_prev_m(3) real(dp), intent(out), optional :: exit_step + integer, intent(inout), optional :: hold_streak + logical, intent(out), optional :: numerical_hold_any - integer :: ktau + integer :: hold_streak_local, ktau real(dp) :: u_ref_prev(3), u_ref_cur(3), x_cur(3), x_cur_m(3) real(dp) :: z_step_start(5), z_step_end(5), segment_duration real(dp) :: boundary_fraction real(dp) :: wall_exit_step - logical :: hit, numerical_hold + logical :: hit, numerical_hold, numerical_hold_any_local call integ_to_ref(z(1:3), u_ref_prev) if (present(exit_step)) exit_step = real(kt, dp) + hold_streak_local = 0 + if (present(hold_streak)) hold_streak_local = hold_streak + numerical_hold_any_local = .false. do ktau = 1, ntau_local numerical_hold = .false. z_step_start = z + if (hold_streak_local /= 0) then + ierr_orbit = hold_streak_local + exit + end if if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr_orbit) if (ierr_orbit == 1) then @@ -1542,6 +1630,15 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & end if end if exit + else if (ierr_orbit /= 0 .and. & + symplectic_newton_warning_mode) then + if (hold_streak_local == 0) then + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. + ierr_orbit = 0 + end if end if else if (swcoll) call update_momentum(anorb, z) @@ -1569,17 +1666,21 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & end if if (ierr_orbit .ne. 0 .and. & symplectic_newton_warning_mode) then + if (hold_streak_local == 0) then call count_event(EVT_WARNING_STEP_SKIP) - ierr_orbit = 0 + hold_streak_local = ierr_orbit numerical_hold = .true. + numerical_hold_any_local = .true. + ierr_orbit = 0 + end if else if (ierr_orbit == 0) then call to_standard_z_coordinates(anorb, z) + hold_streak_local = 0 end if if (ierr_orbit .ne. 0) exit end if if (ierr_orbit .ne. 0) exit if (numerical_hold) then - if (swcoll) call collide(z, dtaumin) kt = kt + 1 cycle end if @@ -1598,9 +1699,13 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & if (swcoll) call collide(z, dtaumin) kt = kt + 1 + if (integmode <= 0) hold_streak_local = 0 x_prev_m = x_cur_m u_ref_prev = u_ref_cur end do + if (present(hold_streak)) hold_streak = hold_streak_local + if (present(numerical_hold_any)) & + numerical_hold_any = numerical_hold_any_local end subroutine macrostep_with_wall_check subroutine to_standard_z_coordinates(anorb, z) @@ -1761,7 +1866,8 @@ subroutine write_results end if close (unit) - open (newunit=unit, file='confined_fraction.dat', recl=1024) + open (newunit=unit, file='confined_fraction.dat', status='replace', & + action='write', recl=1024) do i = 1, ntimstep resolved_count = ntestpart - unresolved_orbits(i) if (resolved_count > 0) then @@ -1776,7 +1882,8 @@ subroutine write_results end do close (unit) - open (newunit=unit, file='unresolved_fraction.dat', recl=1024) + open (newunit=unit, file='unresolved_fraction.dat', status='replace', & + action='write', recl=1024) do i = 1, ntimstep write (unit, *) dble(kt_macro(i))*dtaumin/v0, & real(unresolved_orbits(i), dp)/norm, ntestpart diff --git a/src/spectre_fo_hybrid.f90 b/src/spectre_fo_hybrid.f90 index 10508265..ace27dcb 100644 --- a/src/spectre_fo_hybrid.f90 +++ b/src/spectre_fo_hybrid.f90 @@ -21,6 +21,8 @@ module spectre_fo_hybrid real(dp) :: v(3) = 0.0_dp real(dp) :: u(3) = 0.0_dp real(dp) :: last_y(5) = 0.0_dp + real(dp) :: entry_rho = 0.0_dp + real(dp) :: entry_pzeta = 0.0_dp integer :: owner = 0 integer :: iface = 0 logical :: has_y = .false. @@ -44,6 +46,7 @@ subroutine spectre_fo_enter(state, y, owner, iface, ro0_bar, status) end if state%u = y(1:3) state%last_y = y + state%entry_rho = y(1) state%has_y = .true. state%owner = owner state%iface = iface @@ -53,6 +56,7 @@ subroutine spectre_fo_enter(state, y, owner, iface, ro0_bar, status) vperp = sqrt(max(2.0_dp*y(4)**2 - vpar**2, 0.0_dp)) call perpendicular_frame(b, e1, e2) target = vpar*dot_product(b, e_zeta) + acov(3)/ro0_bar + state%entry_pzeta = target alpha = 0.0_dp if (axisymmetric()) then call solve_gyrophase(x_gc, state%u, owner, b, e1, e2, vpar, & @@ -149,6 +153,13 @@ subroutine spectre_fo_to_gc(state, ro0_bar, y, owner, status, mu_bar) if (status /= SPECTRE_FO_OK) return vpar = dot_product(state%v, b) pzeta = dot_product(state%v, e_zeta) + acov(3)/ro0_bar + ! Standard Boris conserves speed exactly but not the canonical momentum + ! associated with an ignorable toroidal angle. In an axisymmetric field + ! that momentum is an exact physical invariant, including across a + ! discontinuous radial sheet. Use the entry value when mapping back to + ! GC variables so a short recovery excursion cannot create a hidden + ! p_zeta kick. + if (axisymmetric()) pzeta = state%entry_pzeta vperp = state%v - vpar*b rho_l = (ro0_bar/bmod)*cross(b, vperp) x_gc = state%x - rho_l @@ -202,7 +213,12 @@ subroutine check_exit(state, y_gc, owner_gc, exited, status) e_zeta, status) if (status /= SPECTRE_FO_OK) return radial_larmor = max(abs(state%u(1) - y_gc(1)), 1.0e-8_dp) - exited = abs(y_gc(1) - real(state%iface, dp)) > 4.0_dp*radial_larmor + if (state%iface == 0) then + exited = abs(y_gc(1) - state%entry_rho) > 4.0_dp*radial_larmor + else + exited = abs(y_gc(1) - real(state%iface, dp)) > & + 4.0_dp*radial_larmor + end if exited = exited .and. abs(bmod_p - bmod_gc)/bmod_gc <= 0.05_dp exited = exited .and. norm2(bp - bgc) <= 0.05_dp end subroutine check_exit diff --git a/src/spectre_orbit.f90 b/src/spectre_orbit.f90 index b3d09f4c..1306ed63 100644 --- a/src/spectre_orbit.f90 +++ b/src/spectre_orbit.f90 @@ -11,7 +11,8 @@ module spectre_orbit !> SPECTRE_BOUNDARY and terminates the trace. use, intrinsic :: iso_fortran_env, only: dp => real64 - use odeint_allroutines_sub, only: odeint_allroutines + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + use odeint_allroutines_sub, only: odeint_allroutines, odeint_has_failed use alpha_lifetime_sub, only: velo_can use interface_crossing, only: apply_crossing, crossing_info_t, axis_offset, & CROSS_LOSS @@ -28,6 +29,7 @@ module spectre_orbit real(dp), parameter :: rho_tol = 1.0d-10 integer, parameter :: max_bisect = 80 + integer, parameter :: max_retry_depth = 8 !> Home-volume boundaries handed to the field-clamped RHS. Each traced marker !> owns one thread, so these must be thread-private. @@ -127,7 +129,13 @@ subroutine orbit_timestep_spectre(state, z, dtaumin, relerr, level, ierr, event) return end if - call locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, event%t_frac) + call locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, & + event%t_frac, ierr) + if (ierr /= SPECTRE_OK) then + z = z_start + event%occurred = .false. + return + end if iface = nint(boundary) call apply_crossing(z_hit, iface, direction, state%mvol, level, & z, event%info) @@ -163,7 +171,8 @@ subroutine set_home_volume_index(state, lvol) state%home_set = .true. end subroutine set_home_volume_index - subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac) + subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac, & + ierr) !> Find the microstep time where the clamped-RHS trajectory reaches rho_g = !> boundary to |rho_g - boundary| < rho_tol, by the Illinois variant of !> false position. z_start is strictly inside the home volume and the full @@ -171,14 +180,17 @@ subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac) real(dp), intent(in) :: z_start(5), dtaumin, relerr, boundary real(dp), intent(out) :: z_hit(5) real(dp), intent(out) :: t_frac + integer, intent(out) :: ierr real(dp) :: t_lo, t_hi, f_lo, f_hi, t_mid, f_mid, z_mid(5) - integer :: it, ierr, last_side + integer :: it, last_side t_lo = 0.0_dp + ierr = SPECTRE_OK f_lo = z_start(1) - boundary t_hi = dtaumin call integrate_clamped(z_start, t_hi, relerr, z_mid, ierr) + if (ierr /= SPECTRE_OK) return f_hi = z_mid(1) - boundary z_hit = z_mid t_mid = t_hi @@ -188,7 +200,7 @@ subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac) if (abs(f_hi - f_lo) <= tiny(1.0_dp)) exit t_mid = t_hi - f_hi*(t_hi - t_lo)/(f_hi - f_lo) call integrate_clamped(z_start, t_mid, relerr, z_mid, ierr) - if (ierr /= SPECTRE_OK) exit + if (ierr /= SPECTRE_OK) return f_mid = z_mid(1) - boundary z_hit = z_mid if (abs(f_mid) < rho_tol) exit @@ -233,12 +245,34 @@ subroutine integrate_clamped(z_in, tau, relerr, z_out, ierr) real(dp), intent(out) :: z_out(5) integer, intent(out) :: ierr + call integrate_clamped_retry(z_in, tau, relerr, 0, z_out, ierr) + end subroutine integrate_clamped + + recursive subroutine integrate_clamped_retry(z_in, tau, relerr, depth, & + z_out, ierr) + real(dp), intent(in) :: z_in(5), tau, relerr + integer, intent(in) :: depth + real(dp), intent(out) :: z_out(5) + integer, intent(out) :: ierr + + real(dp) :: z_half(5) integer, parameter :: ndim = 5 ierr = SPECTRE_OK z_out = z_in if (tau <= 0.0_dp) return call odeint_allroutines(z_out, ndim, 0.0_dp, tau, relerr, velo_can_clamped) - end subroutine integrate_clamped + if (.not. odeint_has_failed() .and. all(ieee_is_finite(z_out))) return + + z_out = z_in + ierr = SPECTRE_FAULT + if (depth >= max_retry_depth) return + call integrate_clamped_retry(z_in, 0.5_dp*tau, relerr, depth + 1, & + z_half, ierr) + if (ierr /= SPECTRE_OK) return + call integrate_clamped_retry(z_half, 0.5_dp*tau, relerr, depth + 1, & + z_out, ierr) + if (ierr /= SPECTRE_OK) z_out = z_in + end subroutine integrate_clamped_retry end module spectre_orbit diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index 48d99847..c1de6f64 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -11,6 +11,7 @@ module spectre_sympl_orbit !> dtaumin and step-halving keeps the scheme's order across crossings. use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use util, only: twopi, sqrt2 use parmot_mod, only: ro0 use field_can_mod, only: field_can_t, eval_field => evaluate, integ_to_ref, & @@ -20,7 +21,7 @@ module spectre_sympl_orbit symplectic_newton_warning_mode use orbit_symplectic, only: orbit_timestep_sympl, orbit_sympl_init use interface_crossing, only: apply_crossing, crossing_info_t, & - crossing_log_record, CROSS_LOSS, CROSS_STOP, & + crossing_log_record, CROSS_LOSS, CROSS_STOP, CROSS_RECOVERY, & CROSS_REFLECTION, CROSS_SHEET use spectre_sheet_gc, only: sheet_gc_state_t, sheet_gc_initialize, & sheet_gc_advance, sheet_gc_to_y, SHEET_GC_OK @@ -76,14 +77,13 @@ module spectre_sympl_orbit integer, parameter :: max_iters = 200 real(dp), parameter :: h_min_frac = 1.0d-8 real(dp), parameter :: budget_eps = 1.0d-12 - !> A committed substep moving rho_g by more than max_step_dr, or changing - !> the energy 0.5*vpar^2 + mu*B by more than max_step_dh relative, is an - !> unconverged Newton "teleport" (the steppers return silently after maxit): - !> reject it and halve the substep. Both bounds are loose sanity checks that - !> only O(1) garbage states trip; legitimate scheme error stays orders of - !> magnitude below them at any usable step size. + !> Loose state-sanity bounds for a committed substep. The Hamiltonian is not + !> required to stay constant from step to step: symplectic discretizations + !> normally oscillate around a nearby modified Hamiltonian. The shell bound + !> only rejects a state that has left the tested bounded-error band around + !> the fixed physical speed pabs. Interface maps re-pin this shell exactly. real(dp), parameter :: max_step_dr = 0.5d0 - real(dp), parameter :: max_step_dh = 5.0d-5 + real(dp), parameter :: max_relative_shell_defect = 5.0d-4 !> Home-side offset for the per-volume angle-transform evaluation at a landed !> point: the volume dispatch keys on int(rho_g), so exactly at rho_g = k it @@ -103,6 +103,15 @@ module spectre_sympl_orbit integer :: skim_iface = -1 real(dp) :: skim_theta = 0.0_dp real(dp) :: skim_zeta = 0.0_dp + !> An exhausted warning-mode recovery holds one unresolved interval. + !> The unchanged state would reproduce the same deterministic failure, + !> so the following microstep reports a marker-local numerical stop + !> instead of repeating the full retry/FO cascade indefinitely. + logical :: warning_hold_latched = .false. + integer :: warning_hold_count = 0 + integer :: warning_hold_reason = STOP_STEP + integer :: warning_hold_iface = 0 + integer :: warning_hold_direction = 0 type(sheet_gc_state_t) :: sheet type(spectre_fo_state_t) :: fo end type sympl_spectre_state_t @@ -196,7 +205,7 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & integer :: iters, nev, iface, direction, ierr_step, ierr_sheet, ierr_fo integer :: fo_owner integer :: prev_iface - logical :: boundary, solved, sheet_exited, fo_exited + logical :: boundary, solved, sheet_exited, fo_exited, resume_gc ierr = SYMPL_SPECTRE_OK t_frac = 1.0_dp @@ -206,6 +215,15 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & nev = 0 prev_iface = -1 + if (state%warning_hold_latched) then + call record_stop(si, f, ipart, t_base_sec, state%warning_hold_iface, & + state%warning_hold_direction, state%warning_hold_reason) + ierr = SYMPL_SPECTRE_STOP + t_frac = 0.0_dp + si%dt = state%dt_std + return + end if + if (state%fo%active) then call continue_fo(state, si, f, budget, used, y_out, fo_owner, & fo_exited, ierr_fo) @@ -277,23 +295,15 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & do iters = iters + 1 if (iters > max_iters) then - call integrator_state(si, f, y_iface) - iface = max(1, min(state%mvol - 1, nint(si%z(1)))) - call start_fo(state, si, f, y_iface, nint(state%home_hi), & - iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) - if (ierr_fo /= SPECTRE_FO_OK) then - if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & - budget, used, ierr, t_frac)) exit - call finish_fo_error(ierr_fo, si, f, ipart, t_base_sec + & - dt_sec*used/state%dt_std, iface, 1, STOP_STEP, used, & - state%dt_std, ierr, t_frac) - exit - end if - if (.not. fo_exited .or. & - budget <= budget_eps*state%dt_std) exit + call recover_local_error(state, si, f, ipart, t_base_sec + & + dt_sec*used/state%dt_std, STOP_STEP, budget, used, ierr, & + t_frac, resume_gc) + if (resume_gc) then iters = 0 h_try = budget cycle + end if + exit end if si0 = si f0 = f @@ -301,28 +311,20 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & call orbit_timestep_sympl(si, f, ierr_step) if (ierr_step == 0) then - if (step_teleported(si0, f0, si, f)) then + if (step_teleported(si0, si, f)) then si = si0 f = f0 h_try = 0.5_dp*h_try if (h_try >= h_min_frac*state%dt_std) cycle - call integrator_state(si, f, y_iface) - iface = max(1, min(state%mvol - 1, nint(si%z(1)))) - call start_fo(state, si, f, y_iface, nint(state%home_hi), & - iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) - if (ierr_fo /= SPECTRE_FO_OK) then - if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & - budget, used, ierr, t_frac)) exit - call finish_fo_error(ierr_fo, si, f, ipart, & - t_base_sec + dt_sec*used/state%dt_std, iface, 1, & - STOP_TELEPORT, used, state%dt_std, ierr, t_frac) - exit - end if - if (.not. fo_exited .or. & - budget <= budget_eps*state%dt_std) exit + call recover_local_error(state, si, f, ipart, t_base_sec + & + dt_sec*used/state%dt_std, STOP_TELEPORT, budget, used, & + ierr, t_frac, resume_gc) + if (resume_gc) then h_try = budget cycle end if + exit + end if end if call classify_step(state, si0%z(1), si, ierr_step, boundary, iface, & @@ -337,22 +339,14 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & f = f0 h_try = 0.5_dp*h_try if (h_try >= h_min_frac*state%dt_std) cycle - call integrator_state(si, f, y_iface) - iface = max(1, min(state%mvol - 1, nint(si%z(1)))) - call start_fo(state, si, f, y_iface, nint(state%home_hi), & - iface, budget, used, y_out, fo_owner, fo_exited, ierr_fo) - if (ierr_fo /= SPECTRE_FO_OK) then - if (recover_fo_failure(state, si, f, ierr_fo, fo_owner, & - budget, used, ierr, t_frac)) exit - call finish_fo_error(ierr_fo, si, f, ipart, & - t_base_sec + dt_sec*used/state%dt_std, iface, 1, & - STOP_STEP, used, state%dt_std, ierr, t_frac) - exit - end if - if (.not. fo_exited .or. & - budget <= budget_eps*state%dt_std) exit + call recover_local_error(state, si, f, ipart, t_base_sec + & + dt_sec*used/state%dt_std, STOP_STEP, budget, used, ierr, & + t_frac, resume_gc) + if (resume_gc) then h_try = budget cycle + end if + exit end if used = used + h_try budget = state%dt_std - used @@ -653,7 +647,7 @@ subroutine locate_landing(si, f, si0, f0, h_full, rho_k, direction, & evals = evals + 1 g = dir*(si%z(1) - rho_k) if (ierr_step /= 0 .or. g /= g .or. & - step_teleported(si0, f0, si, f)) then + step_teleported(si0, si, f)) then t_hi = t_mid hi_valid = .false. cycle @@ -841,15 +835,23 @@ pure function recanon_pphi(f, p, lambda) result(pphi) pphi = p*lambda*sqrt2*f%hph + f%Aph/f%ro0 end function recanon_pphi - pure function step_teleported(si0, f0, si, f) result(bad) + pure function step_teleported(si0, si, f) result(bad) type(symplectic_integrator_t), intent(in) :: si0, si - type(field_can_t), intent(in) :: f0, f + type(field_can_t), intent(in) :: f logical :: bad - real(dp) :: h0, h1 - + bad = .not. all(ieee_is_finite(si%z)) .or. & + .not. ieee_is_finite(si%pabs) .or. & + .not. ieee_is_finite(f%Bmod) .or. & + .not. ieee_is_finite(f%vpar) .or. & + .not. ieee_is_finite(f%mu) .or. & + .not. ieee_is_finite(f%H) + if (bad) return bad = abs(si%z(1) - si0%z(1)) > max_step_dr if (bad) return + bad = abs(f%H - si%pabs**2) > & + max_relative_shell_defect*max(si%pabs**2, tiny(1.0_dp)) + if (bad) return ! The steppers commit a negative-radius Newton root as an axis chart ! flip (r, theta) -> (|r|, theta + pi). Away from the axis a negative ! root is a Newton divergence, not an axis crossing: a committed theta @@ -858,11 +860,82 @@ pure function step_teleported(si0, f0, si, f) result(bad) bad = si%z(1) > 0.5_dp if (bad) return end if - h0 = 0.5_dp*f0%vpar**2 + f0%mu*f0%Bmod - h1 = 0.5_dp*f%vpar**2 + f%mu*f%Bmod - bad = abs(h1 - h0) > max_step_dh*abs(h0) end function step_teleported + subroutine recover_local_error(state, si, f, ipart, t_sec, reason, budget, & + used, ierr, t_frac, resume_gc) + !> A bulk or axis-local solver failure has no interface metadata. Do not + !> invent an interior interface from nint(rho): interface FO recovery is + !> valid only after an actual landed interface event. + type(sympl_spectre_state_t), intent(inout) :: state + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(in) :: ipart, reason + real(dp), intent(in) :: t_sec + real(dp), intent(inout) :: budget, used + integer, intent(out) :: ierr + real(dp), intent(out) :: t_frac + logical, intent(out) :: resume_gc + + real(dp) :: y(5), y_out(5) + integer :: ierr_fo, owner, owner_out + logical :: fo_exited + + resume_gc = .false. + call integrator_state(si, f, y) + owner = nint(state%home_hi) + call record_bulk_recovery(si, f, ipart, t_sec) + call start_fo(state, si, f, y, owner, 0, budget, used, y_out, & + owner_out, fo_exited, ierr_fo) + if (ierr_fo == SPECTRE_FO_OK) then + resume_gc = fo_exited .and. budget > budget_eps*state%dt_std + return + end if + if (ierr_fo == SPECTRE_FO_LOSS) then + call finish_fo_error(ierr_fo, si, f, ipart, t_sec, 0, 0, reason, & + used, state%dt_std, ierr, t_frac) + return + end if + ! Bulk recovery has no interface or sheet. Preserve the last valid FO + ! state if one exists, but retain the original local failure reason and + ! iface=0 in the marker-local stop record. + if (state%fo%has_y) then + call set_home_volume(state, owner_out) + call recanonicalize(state, si, f, state%fo%last_y) + end if + + if (recover_warning_skip(state, budget, used, ierr, t_frac, reason, & + 0, 0)) return + + call record_stop(si, f, ipart, t_sec, 0, 0, reason) + ierr = SYMPL_SPECTRE_STOP + t_frac = used/state%dt_std + end subroutine recover_local_error + + subroutine record_bulk_recovery(si, f, ipart, t_sec) + type(symplectic_integrator_t), intent(in) :: si + type(field_can_t), intent(in) :: f + integer, intent(in) :: ipart + real(dp), intent(in) :: t_sec + + type(crossing_info_t) :: info + real(dp) :: vpar + + vpar = f%vpar/sqrt2 + info%event_type = CROSS_RECOVERY + info%iface = 0 + info%vol_from = 0 + info%vol_to = 0 + info%theta = si%z(2) + info%zeta = si%z(3) + info%vpar_before = vpar + info%vpar_after = vpar + info%mu = 0.5_dp*f%mu + info%bmod_home = f%Bmod + info%bmod_target = f%Bmod + call crossing_log_record(ipart, t_sec, info) + end subroutine record_bulk_recovery + pure function unwrap_near(a, near) result(b) real(dp), intent(in) :: a, near real(dp) :: b @@ -891,10 +964,12 @@ logical function recover_fo_failure(state, si, f, fo_status, owner, budget, & call recanonicalize(state, si, f, state%fo%last_y) end if recover_fo_failure = recover_warning_skip(state, budget, used, ierr, & - t_frac) + t_frac, STOP_SHEET, state%fo%iface, & + merge(1, -1, owner == state%fo%iface)) end function recover_fo_failure - logical function recover_warning_skip(state, budget, used, ierr, t_frac) + logical function recover_warning_skip(state, budget, used, ierr, t_frac, & + reason, iface, direction) !> Last-resort default: after all active recovery has failed, retain the !> last valid state for this unresolved interval and resume the marker on !> the next microstep. Strict mode returns .false. to expose the stop. @@ -902,6 +977,7 @@ logical function recover_warning_skip(state, budget, used, ierr, t_frac) real(dp), intent(inout) :: budget, used integer, intent(inout) :: ierr real(dp), intent(inout) :: t_frac + integer, intent(in) :: reason, iface, direction recover_warning_skip = symplectic_newton_warning_mode if (.not. recover_warning_skip) return @@ -912,6 +988,12 @@ logical function recover_warning_skip(state, budget, used, ierr, t_frac) t_frac = 1.0_dp call count_event(EVT_WARNING_STEP_SKIP) + state%warning_hold_latched = .true. + state%warning_hold_count = state%warning_hold_count + 1 + state%warning_hold_reason = reason + state%warning_hold_iface = iface + state%warning_hold_direction = direction + state%fo%active = .false. state%sheet%active = .false. end function recover_warning_skip @@ -929,7 +1011,11 @@ subroutine record_stop(si, f, ipart, t_sec, iface, direction, reason) vpar = f%vpar/sqrt2 info%event_type = CROSS_STOP info%iface = iface + if (iface == 0) then + info%vol_from = 0 + else info%vol_from = merge(iface, iface + 1, direction == 1) + end if info%vol_to = info%vol_from info%theta = si%z(2) info%zeta = si%z(3) diff --git a/src/sub_alpha_lifetime_can.f90 b/src/sub_alpha_lifetime_can.f90 index 2eca602d..dae24e82 100644 --- a/src/sub_alpha_lifetime_can.f90 +++ b/src/sub_alpha_lifetime_can.f90 @@ -387,7 +387,8 @@ end subroutine velo_axis !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) - use odeint_allroutines_sub, only : odeint_allroutines + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + use odeint_allroutines_sub, only : odeint_allroutines, odeint_has_failed use chamb_sub, only : chamb_can ! implicit none @@ -401,7 +402,7 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) real(dp) :: relerr ! real(dp), dimension(2) :: y - real(dp), dimension(ndim) :: z + real(dp), dimension(ndim) :: z, z_initial ! if(abs(dtaumin*nstepmax).le.abs(dtau)) then ierr=2 @@ -410,6 +411,7 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) endif ! ierr=0 + z_initial=z y(1)=z(1) y(2)=z(2) phi=z(3) @@ -432,6 +434,11 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) z(2)=z2 ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! y(1)=z(1) y(2)=z(2) @@ -443,6 +450,11 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) else ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_axis) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! endif else @@ -454,10 +466,20 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) z(2)=z2 ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_axis) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! else ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! y(1)=z(1) y(2)=z(2) @@ -484,6 +506,11 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) z(2)=z2 ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! y(1)=z(1) y(2)=z(2) @@ -495,6 +522,11 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) else ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_axis) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! endif else @@ -506,10 +538,20 @@ subroutine orbit_timestep_axis(z,dtau,dtaumin,relerr,ierr) z(2)=z2 ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_axis) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! else ! call odeint_allroutines(z,ndim,tau1,tau2,relerr,velo_can) + if (odeint_has_failed() .or. .not. all(ieee_is_finite(z))) then + z=z_initial + ierr=2 + return + endif ! y(1)=z(1) y(2)=z(2) diff --git a/test/tests/field_can/test_spectre_sympl_crossing.f90 b/test/tests/field_can/test_spectre_sympl_crossing.f90 index 8c40ff74..cfda19a5 100644 --- a/test/tests/field_can/test_spectre_sympl_crossing.f90 +++ b/test/tests/field_can/test_spectre_sympl_crossing.f90 @@ -49,7 +49,7 @@ program test_spectre_sympl_crossing use parmot_mod, only: ro0 use interface_crossing, only: crossing_log_reset, crossing_log_count_type, & CROSS_CROSSING, CROSS_REFLECTION, CROSS_STOP, & - CROSS_SHEET + CROSS_SHEET, CROSS_RECOVERY use util, only: twopi, sqrt2 implicit none @@ -326,7 +326,8 @@ subroutine trace_marker(im, nrec_got, h_series, h_actual, p_series, t_series, & transition_series(nrec_got) = & crossing_log_count_type(CROSS_CROSSING) + & crossing_log_count_type(CROSS_REFLECTION) + & - crossing_log_count_type(CROSS_SHEET) + crossing_log_count_type(CROSS_SHEET) + & + crossing_log_count_type(CROSS_RECOVERY) t_series(nrec_got) = real(k, dp)*dtaumin/v0 end if end do diff --git a/test/tests/test_sympl_testfield.f90 b/test/tests/test_sympl_testfield.f90 index 2bf4a2b3..f0be944f 100644 --- a/test/tests/test_sympl_testfield.f90 +++ b/test/tests/test_sympl_testfield.f90 @@ -44,7 +44,8 @@ program test_sympl_testfield use params, only : isw_field_type, field_input, coord_input, integmode, & swcoll, orbit_model, ORBIT_GC, ORBIT_FULL_ORBIT, ORBIT_EXIT_LCFS, & ORBIT_EXIT_WALL, ORBIT_EXIT_NUMERICAL_DOMAIN, & - ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_FULL_ORBIT + ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_EVENT, & + ORBIT_EXIT_NUMERICAL_FULL_ORBIT use magfie_sub, only : TEST use field_can_mod, only : evaluate use orbit_symplectic, only : orbit_timestep_sympl @@ -126,7 +127,8 @@ subroutine test_failed_step_preserves_state real(dp) :: z(5) real(dp) :: x_previous(3) integer(int64) :: kt - integer :: step_error + integer :: hold_streak, step_error + logical :: numerical_hold swcoll = .false. orbit_model = ORBIT_GC @@ -155,24 +157,44 @@ subroutine test_failed_step_preserves_state z = initial_state kt = 0_int64 + hold_streak = 0 symplectic_newton_warning_mode = .true. - call macrostep(norb, z, kt, step_error, 1) + call macrostep(norb, z, kt, step_error, 1, hold_streak=hold_streak, & + numerical_hold_any=numerical_hold) if (step_error /= 0) error stop 'warning mode stopped a failed step' if (kt /= 1_int64) error stop 'warning mode did not consume the failed step' + if (.not. numerical_hold .or. hold_streak == 0) & + error stop 'warning mode did not latch the unresolved state' if (any(z /= initial_state)) then error stop 'warning-mode skip changed the last accepted state' end if + call macrostep(norb, z, kt, step_error, 1, hold_streak=hold_streak, & + numerical_hold_any=numerical_hold) + if (step_error == 0) & + error stop 'repeated unchanged warning failure was not bounded' + if (kt /= 1_int64) & + error stop 'repeated warning failure advanced the time index' z = initial_state x_previous = 0.0_dp kt = 0_int64 - call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous) + hold_streak = 0 + call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous, & + hold_streak=hold_streak, numerical_hold_any=numerical_hold) if (step_error /= 0) error stop 'warning-mode wall path stopped a failed step' if (kt /= 1_int64) & error stop 'warning-mode wall path did not consume the failed step' + if (.not. numerical_hold .or. hold_streak == 0) & + error stop 'warning-mode wall path did not latch the failure' if (any(z /= initial_state)) then error stop 'warning-mode wall skip changed the last accepted state' end if + call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous, & + hold_streak=hold_streak, numerical_hold_any=numerical_hold) + if (step_error == 0) & + error stop 'wall path repeated an unchanged warning failure' + if (kt /= 1_int64) & + error stop 'repeated wall warning failure advanced the time index' end subroutine test_failed_step_preserves_state subroutine test_exit_classification @@ -216,6 +238,10 @@ subroutine test_exit_classification ORBIT_EXIT_NUMERICAL_DOMAIN) then error stop 'RK extended-map boundary was classified as physical' end if + if (classify_orbit_exit(2, ORBIT_GC, 0, .true.) /= & + ORBIT_EXIT_NUMERICAL_EVENT) then + error stop 'RK integration fault was classified as physical' + end if call classify_classifier_exit(SYMPLECTIC_STEP_MAXITER, 3, & classifier_lost, classifier_exit) @@ -227,6 +253,10 @@ subroutine test_exit_classification if (.not. classifier_lost .or. classifier_exit /= ORBIT_EXIT_LCFS) then error stop 'classifier mode lost its physical boundary classification' end if + call classify_classifier_exit(2, 0, classifier_lost, classifier_exit) + if (classifier_lost .or. classifier_exit /= ORBIT_EXIT_NUMERICAL_EVENT) then + error stop 'classifier mode promoted an RK fault to a physical loss' + end if end subroutine test_exit_classification subroutine test_fo_lcfs_location From bccff5b715bb42f7f2369a27cc1a728f27c36c2e Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 00:34:22 +0200 Subject: [PATCH 04/11] Retry held SPECTRE states before fatal exit --- src/spectre_sympl_orbit.f90 | 39 +++++++++++++---------- test/tests/test_spectre_sympl_crossing.py | 10 +++--- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index c1de6f64..0aa77932 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -104,14 +104,15 @@ module spectre_sympl_orbit real(dp) :: skim_theta = 0.0_dp real(dp) :: skim_zeta = 0.0_dp !> An exhausted warning-mode recovery holds one unresolved interval. - !> The unchanged state would reproduce the same deterministic failure, - !> so the following microstep reports a marker-local numerical stop - !> instead of repeating the full retry/FO cascade indefinitely. + !> The next microstep retries guiding-centre advancement from the last + !> valid state. Only a later recovery with the same reason at the same + !> state is a repeated deterministic failure and becomes marker-local. logical :: warning_hold_latched = .false. integer :: warning_hold_count = 0 integer :: warning_hold_reason = STOP_STEP integer :: warning_hold_iface = 0 integer :: warning_hold_direction = 0 + real(dp) :: warning_hold_z(4) = 0.0_dp type(sheet_gc_state_t) :: sheet type(spectre_fo_state_t) :: fo end type sympl_spectre_state_t @@ -215,15 +216,6 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & nev = 0 prev_iface = -1 - if (state%warning_hold_latched) then - call record_stop(si, f, ipart, t_base_sec, state%warning_hold_iface, & - state%warning_hold_direction, state%warning_hold_reason) - ierr = SYMPL_SPECTRE_STOP - t_frac = 0.0_dp - si%dt = state%dt_std - return - end if - if (state%fo%active) then call continue_fo(state, si, f, budget, used, y_out, fo_owner, & fo_exited, ierr_fo) @@ -904,7 +896,7 @@ subroutine recover_local_error(state, si, f, ipart, t_sec, reason, budget, & call recanonicalize(state, si, f, state%fo%last_y) end if - if (recover_warning_skip(state, budget, used, ierr, t_frac, reason, & + if (recover_warning_skip(state, si, budget, used, ierr, t_frac, reason, & 0, 0)) return call record_stop(si, f, ipart, t_sec, 0, 0, reason) @@ -963,17 +955,18 @@ logical function recover_fo_failure(state, si, f, fo_status, owner, budget, & call set_home_volume(state, owner) call recanonicalize(state, si, f, state%fo%last_y) end if - recover_fo_failure = recover_warning_skip(state, budget, used, ierr, & - t_frac, STOP_SHEET, state%fo%iface, & + recover_fo_failure = recover_warning_skip(state, si, budget, used, & + ierr, t_frac, STOP_SHEET, state%fo%iface, & merge(1, -1, owner == state%fo%iface)) end function recover_fo_failure - logical function recover_warning_skip(state, budget, used, ierr, t_frac, & - reason, iface, direction) + logical function recover_warning_skip(state, si, budget, used, ierr, & + t_frac, reason, iface, direction) !> Last-resort default: after all active recovery has failed, retain the !> last valid state for this unresolved interval and resume the marker on !> the next microstep. Strict mode returns .false. to expose the stop. type(sympl_spectre_state_t), intent(inout) :: state + type(symplectic_integrator_t), intent(in) :: si real(dp), intent(inout) :: budget, used integer, intent(inout) :: ierr real(dp), intent(inout) :: t_frac @@ -982,6 +975,17 @@ logical function recover_warning_skip(state, budget, used, ierr, t_frac, & recover_warning_skip = symplectic_newton_warning_mode if (.not. recover_warning_skip) return + if (state%warning_hold_latched .and. & + reason == state%warning_hold_reason .and. & + iface == state%warning_hold_iface .and. & + direction == state%warning_hold_direction .and. & + maxval(abs(si%z - state%warning_hold_z)) <= & + 100.0_dp*epsilon(1.0_dp)*max(1.0_dp, & + maxval(abs(state%warning_hold_z)))) then + recover_warning_skip = .false. + return + end if + used = used + budget budget = 0.0_dp ierr = SYMPL_SPECTRE_OK @@ -993,6 +997,7 @@ logical function recover_warning_skip(state, budget, used, ierr, t_frac, & state%warning_hold_reason = reason state%warning_hold_iface = iface state%warning_hold_direction = direction + state%warning_hold_z = si%z state%fo%active = .false. state%sheet%active = .false. diff --git a/test/tests/test_spectre_sympl_crossing.py b/test/tests/test_spectre_sympl_crossing.py index 5d02e979..b524ac3b 100644 --- a/test/tests/test_spectre_sympl_crossing.py +++ b/test/tests/test_spectre_sympl_crossing.py @@ -205,13 +205,15 @@ def check_cross_path(binary, h5, failures): p0, _ = loss_fraction(work, trace_time) p_mean = 0.5*(p3 + p0) - sigma = np.sqrt(max(p_mean*(1.0 - p_mean), 1.0/n)/n) - if not abs(p3 - p0) < 3.0*sigma: + # The statistic is the difference of two independent N-marker binomial + # proportions. Its null variance has a contribution from each sample. + sigma_delta = np.sqrt(2.0*max(p_mean*(1.0 - p_mean), 1.0/n)/n) + if not abs(p3 - p0) < 3.0*sigma_delta: failures.append(f"cross-path: loss fractions differ beyond MC error: " f"sympl {p3:.3f} vs RK45 {p0:.3f} (3 sigma = " - f"{3.0*sigma:.3f})") + f"{3.0*sigma_delta:.3f})") print(f"cross-path: loss fraction sympl={p3:.3f} RK45={p0:.3f} " - f"3sigma={3.0*sigma:.3f} (N={n})") + f"3sigma={3.0*sigma_delta:.3f} (N={n})") def orbit_states(workdir): From 86415bf95683f3eb07a3146171f15398860c1029 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 02:14:22 +0200 Subject: [PATCH 05/11] Keep warning-mode orbits advancing --- src/classification.f90 | 38 ++++++-------- src/orbit_symplectic.f90 | 33 ++++++------ src/params.f90 | 6 +-- src/simple_main.f90 | 58 ++++++++------------- src/spectre_sympl_orbit.f90 | 61 ++++++++++++++--------- test/tests/test_newton_solver_status.f90 | 19 ++++--- test/tests/test_spectre_sympl_crossing.py | 26 ++++++---- test/tests/test_sympl_testfield.f90 | 24 ++++++--- 8 files changed, 132 insertions(+), 133 deletions(-) diff --git a/src/classification.f90 b/src/classification.f90 index d901d9c7..a3a6edf9 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -261,34 +261,26 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) cycle endif do ktau=1,ntau - if (hold_streak /= 0) then - ierr = hold_streak - exit - end if if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr) if (ierr /= 0 .and. ierr /= 1 .and. & symplectic_newton_warning_mode) then - if (hold_streak == 0) then - call count_event(EVT_WARNING_STEP_SKIP) - hold_streak = ierr - had_numerical_hold = .true. - if (first_unresolved_it == 0) first_unresolved_it = it - ierr = 0 - end if + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak = ierr + had_numerical_hold = .true. + if (first_unresolved_it == 0) first_unresolved_it = it + ierr = 0 end if else call advance_symplectic_with_retry(anorb%si, anorb%f, & orbit_timestep_sympl, ierr) if (ierr /= 0 .and. ierr /= SYMPLECTIC_STEP_BOUNDARY .and. & symplectic_newton_warning_mode) then - if (hold_streak == 0) then call count_event(EVT_WARNING_STEP_SKIP) - hold_streak = ierr - had_numerical_hold = .true. - if (first_unresolved_it == 0) first_unresolved_it = it + hold_streak = ierr + had_numerical_hold = .true. + if (first_unresolved_it == 0) first_unresolved_it = it ierr = 0 - end if else if (ierr == 0) then hold_streak = 0 z(1:3) = anorb%si%z(1:3) @@ -478,13 +470,13 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) enddo if(ierr.ne.0) exit if (.not. had_numerical_hold) then - if(passing) then - !$omp atomic - confpart_pass(it)=confpart_pass(it)+1.d0 - else - !$omp atomic - confpart_trap(it)=confpart_trap(it)+1.d0 - endif + if(passing) then + !$omp atomic + confpart_pass(it)=confpart_pass(it)+1.d0 + else + !$omp atomic + confpart_trap(it)=confpart_trap(it)+1.d0 + endif end if enddo diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index b2130bac..10fca4f6 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -12,8 +12,7 @@ module orbit_symplectic SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & SYMPLECTIC_STEP_BOUNDARY, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & SYMPLECTIC_STEP_BOUNDARY_LIMITED, boundary_event_fraction_tolerance, & - boundary_event_radial_tolerance, symplectic_newton_warning_mode, & - symplectic_newton_warning_factor + boundary_event_radial_tolerance, symplectic_newton_warning_mode use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, & timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, & timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi @@ -50,18 +49,16 @@ pure logical function finite_newton_iterate(x) finite_newton_iterate = .true. end function finite_newton_iterate -logical function accept_bounded_maxiter(x, xlast, tolref, rtol) - real(dp), intent(in) :: x(:), xlast(:), tolref(:), rtol +logical function accept_warning_maxiter(x) + ! Warning mode preserves SIMPLE's historical production convention: a + ! finite last Newton iterate is usable. The CPU tracing driver still rejects + ! non-finite, radially discontinuous, or off-speed-shell states and routes + ! those through bounded retry/recovery. Strict mode reports maxiter. + real(dp), intent(in) :: x(:) - accept_bounded_maxiter = .false. - if (.not. symplectic_newton_warning_mode) return - if (rtol <= 0.0_dp) return - if (.not. finite_newton_iterate(x)) return - if (.not. finite_newton_iterate(xlast)) return - if (.not. finite_newton_iterate(tolref)) return - accept_bounded_maxiter = all(abs(x - xlast) <= & - symplectic_newton_warning_factor*rtol*abs(tolref)) -end function accept_bounded_maxiter + accept_warning_maxiter = symplectic_newton_warning_mode .and. & + finite_newton_iterate(x) +end function accept_warning_maxiter subroutine advance_symplectic_with_retry(si, f, stepper, status) type(symplectic_integrator_t), intent(inout) :: si @@ -547,7 +544,7 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) end if else call count_event(EVT_NEWTON1_MAXIT) - if (accept_bounded_maxiter(x, xlast, tolref, si%rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -642,7 +639,7 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_NEWTON2_MAXIT) - if (accept_bounded_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -800,7 +797,7 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_MIDPOINT_MAXIT) - if (accept_bounded_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -925,7 +922,7 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta end if else call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_bounded_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_gauss @@ -1267,7 +1264,7 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s end if else call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_bounded_maxiter(x, xlast, tolref, rtol)) & + if (accept_warning_maxiter(x)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_lobatto diff --git a/src/params.f90 b/src/params.f90 index cefc2306..d18bd0e4 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -205,9 +205,9 @@ subroutine read_config(config_file) call reset_seed_if_deterministic if (integmode > 0 .and. symplectic_newton_warning_mode) then - print *, 'WARNING: symplectic integrators accept bounded Newton ', & - 'corrections and retry failed smooth-field steps by halving; ', & - 'see maxit diagnostics' + print *, 'WARNING: symplectic integrators accept finite Newton ', & + 'iterates after maxit; post-step sanity checks and bounded ', & + 'retry remain active; see maxit diagnostics' end if call validate_boundary_event_tolerances diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 90548f98..674677b4 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -1407,10 +1407,6 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & do ktau = 1, ntau_local numerical_hold = .false. z_step_start = z - if (hold_streak_local /= 0) then - ierr_orbit = hold_streak_local - exit - end if if (orbit_model == ORBIT_FULL_ORBIT) then call orbit_timestep_fo(anorb%fo, z, ierr_orbit) if (ierr_orbit .ne. 0) then @@ -1428,14 +1424,12 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & ! one unresolved microstep and retries the next step; ! numerical inversion faults are never physical losses. z = z_step_start - if (hold_streak_local == 0) then - call count_event(EVT_WARNING_STEP_SKIP) - hold_streak_local = ierr_orbit - numerical_hold = .true. - numerical_hold_any_local = .true. + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. ierr_orbit = 0 end if - end if if (ierr_orbit .ne. 0) exit end if if (numerical_hold) then @@ -1460,13 +1454,11 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & ! The RK driver rolls back failed integration to the last ! accepted state. Apply the same bounded warning convention ! as the symplectic and full-orbit paths. - if (hold_streak_local == 0) then - call count_event(EVT_WARNING_STEP_SKIP) - hold_streak_local = ierr_orbit - numerical_hold = .true. - numerical_hold_any_local = .true. - ierr_orbit = 0 - end if + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. + ierr_orbit = 0 end if else if (swcoll) call update_momentum(anorb, z) @@ -1483,13 +1475,11 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, & ! advance_symplectic_with_retry restored the last accepted ! state. Skip only the unresolved microstep and continue; ! strict mode still reports the numerical exit. - if (hold_streak_local == 0) then call count_event(EVT_WARNING_STEP_SKIP) - hold_streak_local = ierr_orbit - numerical_hold = .true. - numerical_hold_any_local = .true. + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. ierr_orbit = 0 - end if else if (ierr_orbit == 0) then call to_standard_z_coordinates(anorb, z) hold_streak_local = 0 @@ -1604,10 +1594,6 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & do ktau = 1, ntau_local numerical_hold = .false. z_step_start = z - if (hold_streak_local /= 0) then - ierr_orbit = hold_streak_local - exit - end if if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr_orbit) if (ierr_orbit == 1) then @@ -1632,13 +1618,11 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & exit else if (ierr_orbit /= 0 .and. & symplectic_newton_warning_mode) then - if (hold_streak_local == 0) then - call count_event(EVT_WARNING_STEP_SKIP) - hold_streak_local = ierr_orbit - numerical_hold = .true. - numerical_hold_any_local = .true. - ierr_orbit = 0 - end if + call count_event(EVT_WARNING_STEP_SKIP) + hold_streak_local = ierr_orbit + numerical_hold = .true. + numerical_hold_any_local = .true. + ierr_orbit = 0 end if else if (swcoll) call update_momentum(anorb, z) @@ -1666,13 +1650,11 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & end if if (ierr_orbit .ne. 0 .and. & symplectic_newton_warning_mode) then - if (hold_streak_local == 0) then call count_event(EVT_WARNING_STEP_SKIP) - hold_streak_local = ierr_orbit + hold_streak_local = ierr_orbit numerical_hold = .true. - numerical_hold_any_local = .true. - ierr_orbit = 0 - end if + numerical_hold_any_local = .true. + ierr_orbit = 0 else if (ierr_orbit == 0) then call to_standard_z_coordinates(anorb, z) hold_streak_local = 0 diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index 0aa77932..aed17f93 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -105,8 +105,9 @@ module spectre_sympl_orbit real(dp) :: skim_zeta = 0.0_dp !> An exhausted warning-mode recovery holds one unresolved interval. !> The next microstep retries guiding-centre advancement from the last - !> valid state. Only a later recovery with the same reason at the same - !> state is a repeated deterministic failure and becomes marker-local. + !> valid state. A later identical failure consumes one unresolved + !> interval directly instead of repeating the full recovery cascade; + !> default warning mode never turns that retry into a marker stop. logical :: warning_hold_latched = .false. integer :: warning_hold_count = 0 integer :: warning_hold_reason = STOP_STEP @@ -291,10 +292,10 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & dt_sec*used/state%dt_std, STOP_STEP, budget, used, ierr, & t_frac, resume_gc) if (resume_gc) then - iters = 0 - h_try = budget - cycle - end if + iters = 0 + h_try = budget + cycle + end if exit end if si0 = si @@ -306,17 +307,22 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & if (step_teleported(si0, si, f)) then si = si0 f = f0 + if (same_warning_hold(state, si, STOP_TELEPORT, 0, 0)) then + if (recover_warning_skip(state, si, budget, used, ierr, & + t_frac, STOP_TELEPORT, 0, 0)) exit + end if h_try = 0.5_dp*h_try if (h_try >= h_min_frac*state%dt_std) cycle call recover_local_error(state, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, STOP_TELEPORT, budget, used, & ierr, t_frac, resume_gc) if (resume_gc) then - h_try = budget - cycle - end if + h_try = budget + cycle + end if exit end if + state%warning_hold_latched = .false. end if call classify_step(state, si0%z(1), si, ierr_step, boundary, iface, & @@ -329,15 +335,19 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & ! branches. si = si0 f = f0 + if (same_warning_hold(state, si, STOP_STEP, 0, 0)) then + if (recover_warning_skip(state, si, budget, used, ierr, & + t_frac, STOP_STEP, 0, 0)) exit + end if h_try = 0.5_dp*h_try if (h_try >= h_min_frac*state%dt_std) cycle call recover_local_error(state, si, f, ipart, t_base_sec + & dt_sec*used/state%dt_std, STOP_STEP, budget, used, ierr, & t_frac, resume_gc) if (resume_gc) then - h_try = budget - cycle - end if + h_try = budget + cycle + end if exit end if used = used + h_try @@ -975,17 +985,6 @@ logical function recover_warning_skip(state, si, budget, used, ierr, & recover_warning_skip = symplectic_newton_warning_mode if (.not. recover_warning_skip) return - if (state%warning_hold_latched .and. & - reason == state%warning_hold_reason .and. & - iface == state%warning_hold_iface .and. & - direction == state%warning_hold_direction .and. & - maxval(abs(si%z - state%warning_hold_z)) <= & - 100.0_dp*epsilon(1.0_dp)*max(1.0_dp, & - maxval(abs(state%warning_hold_z)))) then - recover_warning_skip = .false. - return - end if - used = used + budget budget = 0.0_dp ierr = SYMPL_SPECTRE_OK @@ -1003,6 +1002,20 @@ logical function recover_warning_skip(state, si, budget, used, ierr, & state%sheet%active = .false. end function recover_warning_skip + pure logical function same_warning_hold(state, si, reason, iface, direction) + type(sympl_spectre_state_t), intent(in) :: state + type(symplectic_integrator_t), intent(in) :: si + integer, intent(in) :: reason, iface, direction + + same_warning_hold = state%warning_hold_latched .and. & + reason == state%warning_hold_reason .and. & + iface == state%warning_hold_iface .and. & + direction == state%warning_hold_direction .and. & + maxval(abs(si%z - state%warning_hold_z)) <= & + 100.0_dp*epsilon(1.0_dp)*max(1.0_dp, & + maxval(abs(state%warning_hold_z))) + end function same_warning_hold + subroutine record_stop(si, f, ipart, t_sec, iface, direction, reason) type(symplectic_integrator_t), intent(in) :: si type(field_can_t), intent(in) :: f @@ -1019,7 +1032,7 @@ subroutine record_stop(si, f, ipart, t_sec, iface, direction, reason) if (iface == 0) then info%vol_from = 0 else - info%vol_from = merge(iface, iface + 1, direction == 1) + info%vol_from = merge(iface, iface + 1, direction == 1) end if info%vol_to = info%vol_from info%theta = si%z(2) diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index d807ea0d..2e2edddb 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -128,7 +128,7 @@ program test_newton_solver_status advance_symplectic_with_boundary, advance_symplectic_with_retry, & newton_midpoint, orbit_sympl_init, & orbit_timestep_sympl, matrix3_near_singular, solve_newton_system, & - get_boundary_event_tolerances, accept_bounded_maxiter + get_boundary_event_tolerances, accept_warning_maxiter use orbit_symplectic_base, only: symplectic_integrator_t, & EXPL_IMPL_EULER, IMPL_EXPL_EULER, MIDPOINT, GAUSS1, GAUSS2, GAUSS3, & GAUSS4, LOBATTO3, & @@ -248,7 +248,7 @@ subroutine test_newton_warning_mode MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, LOBATTO3] type(symplectic_integrator_t) :: strict_integrator type(field_can_t) :: strict_field - real(dp) :: initial_state(4), accepted(2), previous(2), scale(2) + real(dp) :: initial_state(4), accepted(2), previous(2) integer :: mode_index, step_status eval_field => evaluate_linear_radial @@ -272,24 +272,23 @@ subroutine test_newton_warning_mode end do previous = [1.0_dp, 2.0_dp] - scale = [1.0_dp, 2.0_dp] accepted = previous + [5.0e-12_dp, 1.0e-11_dp] symplectic_newton_warning_mode = .true. - if (.not. accept_bounded_maxiter(accepted, previous, scale, 1.0e-12_dp)) then - error stop 'warning mode rejected a bounded Newton correction' + if (.not. accept_warning_maxiter(accepted)) then + error stop 'warning mode rejected a finite Newton iterate' end if - accepted(1) = previous(1) + 11.0e-12_dp - if (accept_bounded_maxiter(accepted, previous, scale, 1.0e-12_dp)) then - error stop 'warning mode accepted an excessive Newton correction' + accepted(1) = huge(1.0_dp) + if (.not. accept_warning_maxiter(accepted)) then + error stop 'warning mode rejected a finite large Newton iterate' end if accepted = previous accepted(1) = ieee_value(0.0_dp, ieee_quiet_nan) - if (accept_bounded_maxiter(accepted, previous, scale, 1.0e-12_dp)) then + if (accept_warning_maxiter(accepted)) then error stop 'warning mode accepted a non-finite Newton correction' end if accepted = previous symplectic_newton_warning_mode = .false. - if (accept_bounded_maxiter(accepted, previous, scale, 1.0e-12_dp)) then + if (accept_warning_maxiter(accepted)) then error stop 'strict mode accepted a Newton max-iteration state' end if end subroutine test_newton_warning_mode diff --git a/test/tests/test_spectre_sympl_crossing.py b/test/tests/test_spectre_sympl_crossing.py index b524ac3b..6a2b287a 100644 --- a/test/tests/test_spectre_sympl_crossing.py +++ b/test/tests/test_spectre_sympl_crossing.py @@ -112,9 +112,6 @@ def parse_landing_stats(stdout): r"(\d+)\s+(\d+)\s+(\d+)\s+(\d+)", stdout) if not fo: raise RuntimeError("sympl_fo_stats line missing from stdout") - if int(fo.group(4)) != 0: - raise RuntimeError(f"full-orbit fallback failures: {fo.groups()}\n" - f"{stdout[-4000:]}") print("sheet stats: entries={} exits={} init_fail={} advance_fail={} " "status={},{},{},{},{}" " stop_reason={},{},{},{},{}" @@ -181,6 +178,15 @@ def check_landing_and_reflection(binary, h5, failures): def loss_fraction(workdir, trace_time): + path = os.path.join(workdir, "orbit_exit_code.dat") + if os.path.exists(path): + codes = np.loadtxt(path)[:, 1].astype(int) + resolved = np.isin(codes, (0, 1, 2)) + n_resolved = int(resolved.sum()) + if n_resolved == 0: + raise RuntimeError("cross-path run has no resolved markers") + lost = np.isin(codes, (1, 2)) + return float(lost.sum())/n_resolved, n_resolved tl = load_times_lost(workdir) lost = tl[:, 1] < trace_time*(1.0 - 1e-9) return float(lost.sum())/len(tl), len(tl) @@ -197,23 +203,25 @@ def check_cross_path(binary, h5, failures): ev3 = load_events(work) _, _, stops = parse_landing_stats(out) check_stops(ev3, stops, "cross-path sympl", failures) - p3, n = loss_fraction(work, trace_time) + p3, n3 = loss_fraction(work, trace_time) with tempfile.TemporaryDirectory() as work: run_simple(binary, work, h5=h5, integmode=0, npart=npart, trace_time=trace_time, sbeg=0.5, npoiper2=256, relerr="1d-8", face_al=14.0) - p0, _ = loss_fraction(work, trace_time) + p0, n0 = loss_fraction(work, trace_time) p_mean = 0.5*(p3 + p0) - # The statistic is the difference of two independent N-marker binomial - # proportions. Its null variance has a contribution from each sample. - sigma_delta = np.sqrt(2.0*max(p_mean*(1.0 - p_mean), 1.0/n)/n) + # The statistic is the difference of two independent resolved-marker + # binomial proportions. Its null variance has a contribution from each + # sample, which may have different resolved denominators. + variance = max(p_mean*(1.0 - p_mean), 1.0/max(n3, n0)) + sigma_delta = np.sqrt(variance*(1.0/n3 + 1.0/n0)) if not abs(p3 - p0) < 3.0*sigma_delta: failures.append(f"cross-path: loss fractions differ beyond MC error: " f"sympl {p3:.3f} vs RK45 {p0:.3f} (3 sigma = " f"{3.0*sigma_delta:.3f})") print(f"cross-path: loss fraction sympl={p3:.3f} RK45={p0:.3f} " - f"3sigma={3.0*sigma_delta:.3f} (N={n})") + f"3sigma={3.0*sigma_delta:.3f} (N={n3},{n0})") def orbit_states(workdir): diff --git a/test/tests/test_sympl_testfield.f90 b/test/tests/test_sympl_testfield.f90 index f0be944f..7d01b97b 100644 --- a/test/tests/test_sympl_testfield.f90 +++ b/test/tests/test_sympl_testfield.f90 @@ -170,10 +170,14 @@ subroutine test_failed_step_preserves_state end if call macrostep(norb, z, kt, step_error, 1, hold_streak=hold_streak, & numerical_hold_any=numerical_hold) - if (step_error == 0) & - error stop 'repeated unchanged warning failure was not bounded' - if (kt /= 1_int64) & - error stop 'repeated warning failure advanced the time index' + if (step_error /= 0) & + error stop 'repeated warning failure stopped the orbit' + if (kt /= 2_int64) & + error stop 'repeated warning failure did not consume one interval' + if (.not. numerical_hold .or. hold_streak == 0) & + error stop 'repeated warning failure lost its unresolved marker' + if (any(z /= initial_state)) & + error stop 'repeated warning hold changed the accepted state' z = initial_state x_previous = 0.0_dp @@ -191,10 +195,14 @@ subroutine test_failed_step_preserves_state end if call macrostep_with_wall_check(norb, z, kt, step_error, 1, 1, x_previous, & hold_streak=hold_streak, numerical_hold_any=numerical_hold) - if (step_error == 0) & - error stop 'wall path repeated an unchanged warning failure' - if (kt /= 1_int64) & - error stop 'repeated wall warning failure advanced the time index' + if (step_error /= 0) & + error stop 'repeated wall warning failure stopped the orbit' + if (kt /= 2_int64) & + error stop 'repeated wall warning did not consume one interval' + if (.not. numerical_hold .or. hold_streak == 0) & + error stop 'repeated wall warning lost its unresolved marker' + if (any(z /= initial_state)) & + error stop 'repeated wall warning changed the accepted state' end subroutine test_failed_step_preserves_state subroutine test_exit_classification From c03b6a517185814f7bce6b0d90e4fceab78b1e4b Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 03:32:29 +0200 Subject: [PATCH 06/11] Make numerical recovery nonfatal by default --- DOC/config.md | 17 ++-- python/pysimple/plotting.py | 98 +++++++++++++++++++----- src/classification.f90 | 29 ++----- src/orbit_symplectic.f90 | 94 ++++++++++++++++++----- src/params.f90 | 10 +-- src/progress_monitor.f90 | 42 +++++----- src/simple.f90 | 32 +++++--- src/simple_main.f90 | 40 +++++----- src/spectre_orbit.f90 | 14 +++- test/python/test_plotting.py | 54 +++++++++++++ test/tests/test_newton_solver_status.f90 | 54 ++++++++++--- test/tests/test_spectre_crossing_l1.py | 11 ++- test/tests/test_spectre_validation.py | 6 ++ test/tests/test_sympl_testfield.f90 | 6 +- 14 files changed, 368 insertions(+), 139 deletions(-) diff --git a/DOC/config.md b/DOC/config.md index 9229f21f..e2010f72 100644 --- a/DOC/config.md +++ b/DOC/config.md @@ -17,14 +17,15 @@ tolerance. Each occurrence is reported by the corresponding `*_maxit` diagnostic. The production RK, symplectic, and full-orbit paths then use the same terminal convention: if bounded recovery cannot resolve a numerical - microstep, SIMPLE retains the last accepted state for that interval, records - `warning_step_skip`, and marks the interval unresolved. Because an unchanged - state would reproduce the same deterministic failure, the next microstep uses - a circuit breaker instead of repeating the complete retry cascade. It ends - only that marker with a 101--105 `orbit_exit_code` and `NaN` in `times_lost`; - the ensemble continues, and the marker is neither confined nor physically - lost. Set the option to `.false.` for strict diagnostic runs that expose the - first exhausted recovery without the one-interval warning hold. + microstep, SIMPLE retains any contiguous accepted prefix and holds only the + unresolved remainder, records `warning_step_skip`, advances the clock for + the complete interval, and retries from that valid state on the next + microstep. A warning hold does not terminate or + numerically disqualify the marker: a marker that reaches the requested end + time remains a resolved survivor, while a later physical boundary event is + still a loss. Set the option to `.false.` for strict diagnostic runs that end + only the affected marker at the first exhausted recovery and report a + 101--105 `orbit_exit_code` with `NaN` in `times_lost`. * `canonical_grid_nr`, `canonical_grid_ntheta`, and `canonical_grid_nphi` control the Meiss or Albert canonical-map grid. Their defaults are 62, 63, diff --git a/python/pysimple/plotting.py b/python/pysimple/plotting.py index 63a7c71b..806532b7 100644 --- a/python/pysimple/plotting.py +++ b/python/pysimple/plotting.py @@ -112,6 +112,9 @@ class LossData: Confined trapped fraction at each time. trace_time : float Maximum tracing time [s]. + orbit_exit_codes : np.ndarray, optional + Per-marker exit codes. Codes 0, 1, and 2 are resolved completion, + LCFS loss, and wall loss; codes 101--105 are numerical terminations. """ n_particles: int @@ -127,6 +130,7 @@ class LossData: confined_pass: np.ndarray confined_trap: np.ndarray trace_time: float + orbit_exit_codes: Optional[np.ndarray] = None @property def confined_fraction(self) -> np.ndarray: @@ -135,19 +139,42 @@ def confined_fraction(self) -> np.ndarray: @property def lost_mask(self) -> np.ndarray: - """Boolean mask for particles that were lost (not confined).""" + """Boolean mask for resolved physical LCFS or wall losses.""" + if self.orbit_exit_codes is not None: + return np.isin(self.orbit_exit_codes, (1, 2)) return (self.loss_times > 0) & (self.loss_times < self.trace_time) @property def confined_mask(self) -> np.ndarray: """Boolean mask for particles that remained confined.""" + if self.orbit_exit_codes is not None: + # Code 3 is SIMPLE's deliberate deep-passing shortcut. Fortran + # counts those analytically confined markers in the same resolved + # denominator, while skipped_mask keeps the shortcut observable. + return np.isin(self.orbit_exit_codes, (0, 3)) return self.loss_times >= self.trace_time @property def skipped_mask(self) -> np.ndarray: """Boolean mask for particles skipped (deep passing, contr_pp).""" + if self.orbit_exit_codes is not None: + return self.orbit_exit_codes == 3 return self.loss_times < 0 + @property + def resolved_mask(self) -> np.ndarray: + """Markers eligible for confinement and loss statistics.""" + if self.orbit_exit_codes is not None: + return np.isin(self.orbit_exit_codes, (0, 1, 2, 3)) + return self.lost_mask | self.confined_mask + + @property + def unresolved_mask(self) -> np.ndarray: + """Markers terminated for numerical reasons.""" + if self.orbit_exit_codes is not None: + return (self.orbit_exit_codes >= 101) & (self.orbit_exit_codes <= 105) + return ~np.isfinite(self.loss_times) + def load_loss_data(directory: str | Path) -> LossData: """ @@ -209,7 +236,8 @@ def load_loss_data(directory: str | Path) -> LossData: start_phi = np.zeros(n_particles) start_pitch = np.zeros(n_particles) - # Load confined_fraction.dat + # Load confined_fraction.dat. Its final time is authoritative: numerical + # NaNs in times_lost.dat must not poison or extend the requested trace time. # Columns: time, confpart_pass, confpart_trap, resolved particle count conf_file = directory / "confined_fraction.dat" if conf_file.exists(): @@ -218,19 +246,42 @@ def load_loss_data(directory: str | Path) -> LossData: time_grid = np.abs(conf_data[1:, 0]) confined_pass = conf_data[1:, 1] confined_trap = conf_data[1:, 2] + trace_time = float(time_grid[-1]) if len(time_grid) > 0 else 1.0 + else: + finite_times = np.abs(loss_times[np.isfinite(loss_times)]) + trace_time = float(np.max(finite_times)) if finite_times.size else 1.0 + time_grid = np.logspace(-5, np.log10(max(trace_time, 1.0e-5)), 100) + + exit_file = directory / "orbit_exit_code.dat" + if exit_file.exists(): + exit_data = np.loadtxt(exit_file) + exit_data = np.atleast_2d(exit_data) + if exit_data.shape[0] != n_particles: + raise ValueError("orbit_exit_code.dat particle count does not match times_lost.dat") + orbit_exit_codes = exit_data[:, 1].astype(int) else: - # Generate from loss times if file not available - time_grid = np.logspace(-5, 0, 100) + # Backward compatibility for output predating explicit exit codes. + orbit_exit_codes = np.full(n_particles, 101, dtype=int) + orbit_exit_codes[loss_times < 0.0] = 3 + orbit_exit_codes[np.isfinite(loss_times) & (loss_times >= trace_time)] = 0 + orbit_exit_codes[(loss_times > 0.0) & (loss_times < trace_time)] = 1 + + if not conf_file.exists(): confined_pass = np.zeros_like(time_grid) confined_trap = np.zeros_like(time_grid) + resolved = np.isin(orbit_exit_codes, (0, 1, 2, 3)) + passing = trap_parameter < 0 + denominator = int(np.sum(resolved)) for i, t in enumerate(time_grid): - confined = np.abs(loss_times) >= t - passing = trap_parameter < 0 - confined_pass[i] = np.sum(confined & passing) / n_particles - confined_trap[i] = np.sum(confined & ~passing) / n_particles - - # Determine trace_time from maximum loss time or last time in grid - trace_time = max(np.max(np.abs(loss_times)), time_grid[-1] if len(time_grid) > 0 else 1.0) + confined = resolved & ( + np.isin(orbit_exit_codes, (0, 3)) | (loss_times >= t) + ) + if denominator > 0: + confined_pass[i] = np.sum(confined & passing) / denominator + confined_trap[i] = np.sum(confined & ~passing) / denominator + else: + confined_pass[i] = np.nan + confined_trap[i] = np.nan return LossData( n_particles=n_particles, @@ -246,6 +297,7 @@ def load_loss_data(directory: str | Path) -> LossData: confined_pass=confined_pass, confined_trap=confined_trap, trace_time=trace_time, + orbit_exit_codes=orbit_exit_codes, ) @@ -281,7 +333,9 @@ def compute_energy_confined_fraction( if time_grid is None: time_grid = data.time_grid - e_total = float(data.n_particles) + e_total = float(np.sum(data.resolved_mask)) + if e_total == 0.0: + return time_grid, np.full_like(time_grid, np.nan, dtype=float) if slowing_down_curve is not None: times_sd, energy_sd = slowing_down_curve @@ -295,7 +349,7 @@ def compute_energy_confined_fraction( energy_fraction = np.zeros_like(time_grid) for i, t in enumerate(time_grid): - lost_before_t = (data.loss_times > 0) & (data.loss_times < t) + lost_before_t = data.lost_mask & (data.loss_times < t) e_lost = np.sum(energy_at_loss[lost_before_t]) energy_fraction[i] = (e_total - e_lost) / e_total @@ -412,7 +466,7 @@ def _plot_kde_density_panel(ax, data: LossData) -> None: def _plot_energy_loss_panel(ax, data: LossData) -> None: """Plot energy loss distribution vs J_perp panel.""" jperp_bins, particle_count, energy_lost = compute_energy_loss_distribution(data) - total = data.n_particles if data.n_particles > 0 else 1 + total = max(int(np.sum(data.resolved_mask)), 1) energy_frac = energy_lost / total particle_frac = particle_count / total @@ -426,7 +480,10 @@ def _plot_energy_loss_panel(ax, data: LossData) -> None: def _plot_starting_positions_panel(ax, data: LossData) -> None: """Plot starting positions colored by loss time panel.""" - tlost_all = np.maximum(np.abs(data.loss_times), 1e-10) + tlost_all = np.where( + np.isfinite(data.loss_times), np.abs(data.loss_times), data.trace_time + ) + tlost_all = np.maximum(tlost_all, 1e-10) scatter = ax.scatter(data.start_theta, data.start_pitch, c=np.log10(tlost_all), s=1, cmap="viridis", vmin=-5, vmax=0) cbar = plt.colorbar(scatter, ax=ax) cbar.set_label(r"$\log_{10}(t_{\mathrm{loss}})$") @@ -650,7 +707,8 @@ def _bin_energy_actual( ) -> Tuple[np.ndarray, np.ndarray]: """Bin particles and energy by J_perp using actual final_p^2 (vectorized).""" k = _compute_bin_indices(data.perp_invariant, hp, nperp) - part_distr = np.bincount(k, minlength=nperp).astype(float) + part_distr = np.bincount(k, weights=data.resolved_mask.astype(float), + minlength=nperp) energy_weights = np.where(data.lost_mask, data.final_p**2, 0.0) energ_distr = np.bincount(k, weights=energy_weights, minlength=nperp) return part_distr, energ_distr @@ -665,7 +723,8 @@ def _bin_energy_theoretical( ) -> Tuple[np.ndarray, np.ndarray]: """Bin particles and energy by J_perp using theoretical slowing-down (vectorized).""" k = _compute_bin_indices(data.perp_invariant, hp, nperp) - part_distr = np.bincount(k, minlength=nperp).astype(float) + part_distr = np.bincount(k, weights=data.resolved_mask.astype(float), + minlength=nperp) energy_adj = energy_sd - THERMAL_ENERGY_FRACTION energy_at_loss = np.interp(data.loss_times, times_sd, energy_adj, left=energy_adj[0], right=energy_adj[-1]) @@ -743,13 +802,14 @@ def plot_energy_loss_vs_jperp( _, energ_n_theo = _bin_energy_theoretical(data_nocoll, hp, nperp, times_sd, energy_sd) part_c_safe = np.where(part_c > 0, part_c, 1) + part_n_safe = np.where(part_n > 0, part_n, 1) if part_n is not None else None jperp = np.arange(1, nperp + 1) / nperp curves = [ (energ_c / part_c_safe, "r-", 2, "COLL (actual)"), - (energ_n / part_c_safe if energ_n is not None else None, "b-", 2, "NOCOLL (actual)"), + (energ_n / part_n_safe if energ_n is not None else None, "b-", 2, "NOCOLL (actual)"), (energ_c_theo / part_c_safe if energ_c_theo is not None else None, "r--", 1.5, "COLL (theoretical)"), - (energ_n_theo / part_c_safe if energ_n_theo is not None else None, "b--", 1.5, "NOCOLL (theoretical)"), + (energ_n_theo / part_n_safe if energ_n_theo is not None else None, "b--", 1.5, "NOCOLL (theoretical)"), ] max_f = 0.0 diff --git a/src/classification.f90 b/src/classification.f90 index a3a6edf9..beb34611 100644 --- a/src/classification.f90 +++ b/src/classification.f90 @@ -79,7 +79,7 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) real(dp), dimension(3) :: bder, hcovar, hctrvr, hcurl integer :: first_unresolved_it, hold_streak, it, ktau, it_f integer(8) :: kt - logical :: had_numerical_hold, passing + logical :: passing integer :: ifp_tip,ifp_per integer, dimension(:), allocatable :: ipoi @@ -217,7 +217,6 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) kt = 0 first_unresolved_it = 0 hold_streak = 0 - had_numerical_hold = .false. if (passing) then !$omp atomic confpart_pass(1)=confpart_pass(1)+1.d0 @@ -267,8 +266,6 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) symplectic_newton_warning_mode) then call count_event(EVT_WARNING_STEP_SKIP) hold_streak = ierr - had_numerical_hold = .true. - if (first_unresolved_it == 0) first_unresolved_it = it ierr = 0 end if else @@ -278,8 +275,6 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) symplectic_newton_warning_mode) then call count_event(EVT_WARNING_STEP_SKIP) hold_streak = ierr - had_numerical_hold = .true. - if (first_unresolved_it == 0) first_unresolved_it = it ierr = 0 else if (ierr == 0) then hold_streak = 0 @@ -469,23 +464,15 @@ subroutine trace_orbit_with_classifiers(anorb, ipart, class_result) ! write(999, *) kt*dtaumin/v0, z enddo if(ierr.ne.0) exit - if (.not. had_numerical_hold) then - if(passing) then - !$omp atomic - confpart_pass(it)=confpart_pass(it)+1.d0 - else - !$omp atomic - confpart_trap(it)=confpart_trap(it)+1.d0 - endif - end if + if(passing) then + !$omp atomic + confpart_pass(it)=confpart_pass(it)+1.d0 + else + !$omp atomic + confpart_trap(it)=confpart_trap(it)+1.d0 + endif enddo - if (had_numerical_hold) then - class_result%lost = .false. - if (class_result%exit_code < ORBIT_EXIT_NUMERICAL_DOMAIN) & - class_result%exit_code = ORBIT_EXIT_NUMERICAL_EVENT - end if - !$omp critical zend(:,ipart) = z if(isw_field_type .eq. CANFLUX) then diff --git a/src/orbit_symplectic.f90 b/src/orbit_symplectic.f90 index 10fca4f6..f8d113e1 100644 --- a/src/orbit_symplectic.f90 +++ b/src/orbit_symplectic.f90 @@ -12,7 +12,8 @@ module orbit_symplectic SYMPLECTIC_STEP_MAXITER, SYMPLECTIC_STEP_LINEAR_SOLVE, & SYMPLECTIC_STEP_BOUNDARY, SYMPLECTIC_STEP_EVENT_NOT_CONVERGED, & SYMPLECTIC_STEP_BOUNDARY_LIMITED, boundary_event_fraction_tolerance, & - boundary_event_radial_tolerance, symplectic_newton_warning_mode + boundary_event_radial_tolerance, symplectic_newton_warning_mode, & + symplectic_newton_warning_factor use orbit_symplectic_quasi, only: orbit_timestep_quasi, timestep_expl_impl_euler_quasi, & timestep_impl_expl_euler_quasi, timestep_midpoint_quasi, orbit_timestep_rk45, & timestep_rk_gauss_quasi, timestep_rk_lobatto_quasi @@ -24,7 +25,7 @@ module orbit_symplectic use lapack_interfaces, only: dgesv use diag_counters, only: count_event, EVT_NEWTON1_MAXIT, EVT_NEWTON2_MAXIT, & EVT_RK_GAUSS_MAXIT, EVT_RK_LOBATTO_MAXIT, EVT_FIXPOINT_MAXIT, EVT_R_NEGATIVE, & - EVT_MIDPOINT_MAXIT + EVT_MIDPOINT_MAXIT, EVT_WARNING_STEP_SKIP implicit none @@ -49,28 +50,38 @@ pure logical function finite_newton_iterate(x) finite_newton_iterate = .true. end function finite_newton_iterate -logical function accept_warning_maxiter(x) - ! Warning mode preserves SIMPLE's historical production convention: a - ! finite last Newton iterate is usable. The CPU tracing driver still rejects - ! non-finite, radially discontinuous, or off-speed-shell states and routes - ! those through bounded retry/recovery. Strict mode reports maxiter. - real(dp), intent(in) :: x(:) - - accept_warning_maxiter = symplectic_newton_warning_mode .and. & - finite_newton_iterate(x) +logical function accept_warning_maxiter(x, xlast, tolref, rtol) + ! Continue directly only when the last Newton correction is finite and close + ! to the requested tolerance. Larger finite iterates are not trustworthy + ! states: the warning-mode driver retries them at smaller substeps and, if + ! every recovery is exhausted, advances the clock from the last valid state. + real(dp), intent(in) :: x(:), xlast(:), tolref(:), rtol + + accept_warning_maxiter = .false. + if (.not. symplectic_newton_warning_mode .or. rtol <= 0.0_dp) return + if (.not. finite_newton_iterate(x)) return + if (.not. finite_newton_iterate(xlast)) return + if (.not. finite_newton_iterate(tolref)) return + accept_warning_maxiter = all(abs(x - xlast) <= & + symplectic_newton_warning_factor*rtol*abs(tolref)) end function accept_warning_maxiter -subroutine advance_symplectic_with_retry(si, f, stepper, status) +subroutine advance_symplectic_with_retry(si, f, stepper, status, & + accepted_fraction) type(symplectic_integrator_t), intent(inout) :: si type(field_can_t), intent(inout) :: f procedure(orbit_timestep_sympl_i) :: stepper integer, intent(out) :: status + real(dp), intent(out), optional :: accepted_fraction + real(dp) :: accepted_fraction_local - call advance_retry_interval(si, f, stepper, si%dt, 0, status) + call advance_retry_interval(si, f, stepper, si%dt, 0, status, & + accepted_fraction_local) + if (present(accepted_fraction)) accepted_fraction = accepted_fraction_local end subroutine advance_symplectic_with_retry recursive subroutine advance_retry_interval(si, f, stepper, interval, depth, & - status) + status, accepted_fraction) integer, parameter :: max_retry_depth = 8 type(symplectic_integrator_t), intent(inout) :: si type(field_can_t), intent(inout) :: f @@ -78,16 +89,26 @@ recursive subroutine advance_retry_interval(si, f, stepper, interval, depth, & real(dp), intent(in) :: interval integer, intent(in) :: depth integer, intent(out) :: status + real(dp), intent(out) :: accepted_fraction type(symplectic_integrator_t) :: initial_integrator + type(symplectic_integrator_t) :: half_integrator type(field_can_t) :: initial_field + type(field_can_t) :: half_field + real(dp) :: child_fraction logical :: recoverable + accepted_fraction = 0.0_dp initial_integrator = si initial_field = f si%dt = interval call stepper(si, f, status) if (status == SYMPLECTIC_STEP_OK .or. & status == SYMPLECTIC_STEP_BOUNDARY) then + if (status == SYMPLECTIC_STEP_OK) then + accepted_fraction = 1.0_dp + else + accepted_fraction = si%last_step_fraction + end if si%dt = interval return end if @@ -104,22 +125,53 @@ recursive subroutine advance_retry_interval(si, f, stepper, interval, depth, & ! the original interval exactly. State-dependent subdivision is a warning-mode ! recovery path, not a globally fixed symplectic map. call advance_retry_interval(si, f, stepper, 0.5_dp*interval, depth + 1, & - status) + status, child_fraction) if (status == SYMPLECTIC_STEP_BOUNDARY) then si%last_step_fraction = 0.5_dp*si%last_step_fraction si%last_event_fraction_width = 0.5_dp*si%last_event_fraction_width + accepted_fraction = 0.5_dp*child_fraction si%dt = interval return end if if (status == SYMPLECTIC_STEP_OK) then + if (child_fraction < 1.0_dp) then + ! A nested recovery already consumed the unresolved suffix of this first + ! half. Preserve that contiguous accepted prefix and do not leap over the + ! held interval to attempt the nominal second half. + accepted_fraction = 0.5_dp*child_fraction + si%dt = interval + return + end if + half_integrator = si + half_field = f call advance_retry_interval(si, f, stepper, 0.5_dp*interval, depth + 1, & - status) + status, child_fraction) if (status == SYMPLECTIC_STEP_BOUNDARY) then si%last_step_fraction = 0.5_dp + 0.5_dp*si%last_step_fraction si%last_event_fraction_width = 0.5_dp*si%last_event_fraction_width + accepted_fraction = 0.5_dp + 0.5_dp*child_fraction si%dt = interval return end if + if (status == SYMPLECTIC_STEP_OK) then + accepted_fraction = 0.5_dp + 0.5_dp*child_fraction + si%dt = interval + return + end if + recoverable = status == SYMPLECTIC_STEP_MAXITER .or. & + status == SYMPLECTIC_STEP_LINEAR_SOLVE + if (recoverable) then + ! Keep the first accepted half rather than discarding meaningful orbit + ! progress because the second half exhausted every retry. Warning mode + ! consumes only that unresolved remainder and resumes from this state. + si = half_integrator + f = half_field + si%dt = interval + status = SYMPLECTIC_STEP_OK + accepted_fraction = 0.5_dp + call count_event(EVT_WARNING_STEP_SKIP) + return + end if end if if (status /= SYMPLECTIC_STEP_OK) then si = initial_integrator @@ -544,7 +596,7 @@ recursive subroutine newton1(si, f, x, maxit, xlast, status) end if else call count_event(EVT_NEWTON1_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, si%rtol)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -639,7 +691,7 @@ recursive subroutine newton2(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_NEWTON2_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -797,7 +849,7 @@ recursive subroutine newton_midpoint(si, f, x, atol, rtol, maxit, xlast, status) end if else call count_event(EVT_MIDPOINT_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol)) & status = SYMPLECTIC_STEP_OK end if end subroutine @@ -922,7 +974,7 @@ recursive subroutine newton_rk_gauss(si, fs, s, x, atol, rtol, maxit, xlast, sta end if else call count_event(EVT_RK_GAUSS_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_gauss @@ -1264,7 +1316,7 @@ recursive subroutine newton_rk_lobatto(si, fs, s, x, atol, rtol, maxit, xlast, s end if else call count_event(EVT_RK_LOBATTO_MAXIT) - if (accept_warning_maxiter(x)) & + if (accept_warning_maxiter(x, xlast, tolref, rtol)) & status = SYMPLECTIC_STEP_OK end if end subroutine newton_rk_lobatto diff --git a/src/params.f90 b/src/params.f90 index d18bd0e4..1f46491a 100644 --- a/src/params.f90 +++ b/src/params.f90 @@ -82,8 +82,8 @@ module params real(dp) :: relerr = 1d-13 - ! Wall-clock seconds between partial-result checkpoints during tracing. - ! <= 0 disables periodic output (results then appear only at the end). + ! Wall-clock seconds between progress reports during tracing. Result arrays + ! are written only at a team-safe point after tracing; <= 0 disables reports. real(dp) :: checkpoint_interval = 10.0d0 integer :: canonical_grid_nr = 62 integer :: canonical_grid_ntheta = 63 @@ -205,9 +205,9 @@ subroutine read_config(config_file) call reset_seed_if_deterministic if (integmode > 0 .and. symplectic_newton_warning_mode) then - print *, 'WARNING: symplectic integrators accept finite Newton ', & - 'iterates after maxit; post-step sanity checks and bounded ', & - 'retry remain active; see maxit diagnostics' + print *, 'WARNING: symplectic integrators accept bounded Newton ', & + 'corrections after maxit; failed steps retry from the last ', & + 'valid state; see maxit/step-skip diagnostics' end if call validate_boundary_event_tolerances diff --git a/src/progress_monitor.f90 b/src/progress_monitor.f90 index c364b13b..5ac273e9 100644 --- a/src/progress_monitor.f90 +++ b/src/progress_monitor.f90 @@ -1,13 +1,14 @@ module progress_monitor - !> Periodic progress reporting and partial-result checkpointing for the + !> Periodic progress reporting and team-safe final result output for the !> particle tracing loop. !> !> Once per finished particle a thread calls progress_tick. The hot path is !> a single atomic increment of the completed counter and a wall-clock read, !> with no critical and no file access. When `interval` seconds have passed - !> one thread enters a rarely taken critical, writes a status line and the - !> partial results to disk, and advances the deadline. A run killed mid-flight - !> therefore leaves its most recent output on disk instead of nothing. + !> one thread enters a rarely taken critical, writes a status line, and + !> advances the deadline. Result files are written only after the parallel + !> region: a live snapshot would race with trajectory and counter updates and + !> could mix different marker states in one checkpoint. !> !> The monitor never touches physics arrays itself. The owner registers a !> dump callback at init, so this module depends on no tracing code and the @@ -26,7 +27,7 @@ end subroutine dump_proc_i end interface procedure(dump_proc_i), pointer :: dump_results => null() - real(dp) :: interval = 0.0d0 ! <= 0 disables periodic checkpointing + real(dp) :: interval = 0.0d0 ! <= 0 disables periodic progress reports real(dp) :: t_start = 0.0d0 real(dp) :: t_next = 0.0d0 integer :: n_total = 0 @@ -61,13 +62,11 @@ subroutine progress_tick() done = n_done_shared !$omp end atomic - now = wall_time() - if (now < t_next) return - !$omp critical (progress_flush) - if (now >= t_next) then ! double-checked: only one thread per interval + now = wall_time() + if (now >= t_next) then t_next = now + interval - call emit(int(done), now) + call emit(int(done), now, .false.) end if !$omp end critical (progress_flush) end subroutine progress_tick @@ -76,21 +75,22 @@ subroutine progress_finalize() !> Emit one final summary so the event totals (e.g. fo_loss / fo_fault) !> are always reported, including for runs too short to cross a periodic !> interval. A silent fo_fault count would hide inversion faults. - if (n_total > 0) call emit(int(n_done_shared), wall_time()) + if (n_total > 0) call emit(int(n_done_shared), wall_time(), .true.) active = .false. dump_results => null() end subroutine progress_finalize - subroutine emit(n_done, now) + subroutine emit(n_done, now, final) integer, intent(in) :: n_done real(dp), intent(in) :: now + logical, intent(in) :: final real(dp) :: elapsed, frac, eta integer :: id integer(int64) :: total character(len=512) :: events - if (associated(dump_results)) call dump_results() + if (final .and. associated(dump_results)) call dump_results() elapsed = now - t_start frac = real(n_done, dp)/real(max(n_total, 1), dp) @@ -98,13 +98,15 @@ subroutine emit(n_done, now) if (frac > 0.0d0) eta = elapsed*(1.0d0 - frac)/frac events = '' - do id = 1, N_EVENT - total = diag_counters_total(id) - if (total > 0_int64) then - events = trim(events)//' '//event_name(id)//'='// & - trim(int_to_str(total)) - end if - end do + if (final) then + do id = 1, N_EVENT + total = diag_counters_total(id) + if (total > 0_int64) then + events = trim(events)//' '//event_name(id)//'='// & + trim(int_to_str(total)) + end if + end do + end if write (output_unit, '(A, F6.2, A, I0, A, I0, A, F0.1, A, F0.1, A, A)') & ' [progress] ', 100.0d0*frac, '% ', n_done, '/', n_total, & diff --git a/src/simple.f90 b/src/simple.f90 index 7cf6f044..1d9cad39 100644 --- a/src/simple.f90 +++ b/src/simple.f90 @@ -194,14 +194,20 @@ subroutine orbit_timestep_fo(fo, z, ierr) type(fo_state_t), intent(inout) :: fo real(dp), intent(inout) :: z(:) integer, intent(out) :: ierr - real(dp) :: s, th, ph, vpar + type(fo_state_t) :: fo_start + real(dp) :: s, th, ph, vpar, z_start(size(z)) integer :: status + fo_start = fo + z_start = z call fo_step(fo, status) if (status /= FO_OK) then ! Inversion unresolved (near-axis below chartmap resolution, or a field-period ! seam). NOT a physical loss; fo_step left the state at the last resolved - ! position. Hand the caller ierr=3 to end the orbit gracefully as confined. + ! position. Roll back both public and hidden pusher state so warning mode + ! can consume a true held interval and retry from a meaningful position. + fo = fo_start + z = z_start ierr = ORBIT_FO_NUMERICAL call count_event(EVT_FO_FAULT) call warn_fo_unresolved @@ -218,6 +224,8 @@ subroutine orbit_timestep_fo(fo, z, ierr) call count_event(EVT_FO_LOSS) return else if (status /= FO_OK) then + fo = fo_start + z = z_start ierr = ORBIT_FO_NUMERICAL call count_event(EVT_FO_FAULT) call warn_fo_unresolved @@ -229,18 +237,20 @@ subroutine orbit_timestep_fo(fo, z, ierr) end subroutine orbit_timestep_fo ! One-time stderr warning that some full-orbit steps could not invert the - ! Cartesian position and fell back to the last resolved field. The benign race on - ! `warned` across OpenMP threads at worst prints a few extra lines. + ! Cartesian position and were rolled back to the last resolved state. subroutine warn_fo_unresolved use iso_fortran_env, only: error_unit logical, save :: warned = .false. - if (warned) return - warned = .true. - write (error_unit, '(A)') ' WARNING: full-orbit Cartesian inversion unresolved '// & - 'at some steps (near-axis below chartmap resolution, or a field-period seam). '// & - 'Those markers end at the last resolved position and are counted CONFINED '// & - '(fo_fault), never lost. Refine the chartmap near the axis to remove them.' - flush (error_unit) + !$omp critical (fo_unresolved_warning) + if (.not. warned) then + warned = .true. + write (error_unit, '(A)') ' WARNING: full-orbit Cartesian inversion unresolved '// & + 'at some steps (near-axis below chartmap resolution, or a field-period seam). '// & + 'Warning mode holds that interval at the last resolved state and continues '// & + '(fo_fault); it is never counted as a physical loss.' + flush (error_unit) + end if + !$omp end critical (fo_unresolved_warning) end subroutine warn_fo_unresolved subroutine timestep(self, s, th, ph, lam, ierr) diff --git a/src/simple_main.f90 b/src/simple_main.f90 index 674677b4..859fd07b 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -1066,8 +1066,6 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) exit_step, hold_streak=hold_streak, & numerical_hold_any=numerical_hold) end if - if (numerical_hold .and. first_unresolved_it == 0) & - first_unresolved_it = it end if if (ierr_orbit .ne. 0) then @@ -1141,7 +1139,7 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) !> lost at the outermost interface. Every crossing/reflection is logged. use spectre_orbit, only: spectre_orbit_state_t, spectre_event_t, & spectre_state_reset, orbit_timestep_spectre, & - SPECTRE_OK, SPECTRE_BOUNDARY + SPECTRE_OK, SPECTRE_BOUNDARY, SPECTRE_FAULT use magfie_sub, only: spectre_field use interface_crossing, only: crossing_log_record use params, only: crossing_level @@ -1173,6 +1171,17 @@ subroutine trace_orbit_spectre(ipart, z, passing, orbit_traj, orbit_times) call crossing_log_record(ipart, & (real(kt, dp) + event%t_frac)*dtaumin/v0, event%info) end if + if (ierr_orbit == SPECTRE_FAULT .and. & + symplectic_newton_warning_mode) then + ! orbit_timestep_spectre rolls back to the last accepted + ! position. Consume this held interval and retry the same + ! marker on the next microstep; warning mode never turns an + ! isolated RK failure into a terminal marker. + call count_event(EVT_WARNING_STEP_SKIP) + ierr_orbit = SPECTRE_OK + kt = kt + 1 + cycle + end if if (ierr_orbit /= SPECTRE_OK) exit kt = kt + 1 end do @@ -1246,7 +1255,7 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & type(sympl_spectre_state_t) :: state integer :: it, it_f, ktau, ierr_orbit, it_final - integer :: first_unresolved_it, holds_before + integer :: first_unresolved_it integer(8) :: kt real(dp) :: t_stop, t_frac @@ -1262,14 +1271,11 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & do it = 1, ntimstep if (it >= 2) then do ktau = 1, ntau_macro(it) - holds_before = state%warning_hold_count call orbit_microstep_sympl_spectre(state, anorb%si, anorb%f, & ipart, & real(kt, dp)*dtaumin/v0, & dtaumin/v0, ierr_orbit, & t_frac) - if (state%warning_hold_count > holds_before .and. & - first_unresolved_it == 0) first_unresolved_it = it if (ierr_orbit /= SYMPL_SPECTRE_OK) exit kt = kt + 1 end do @@ -1327,10 +1333,7 @@ subroutine trace_orbit_spectre_sympl(anorb, ipart, z, passing, orbit_traj, & end do end if - if (state%warning_hold_count > 0) then - t_stop = ieee_value(0.0_dp, ieee_quiet_nan) - orbit_exit_code(ipart) = ORBIT_EXIT_NUMERICAL_EVENT - else if (ierr_orbit == SYMPL_SPECTRE_OK) then + if (ierr_orbit == SYMPL_SPECTRE_OK) then t_stop = real(kt, dp)*dtaumin/v0 orbit_exit_code(ipart) = ORBIT_EXIT_COMPLETED else if (ierr_orbit == SYMPL_SPECTRE_SKIM) then @@ -1593,6 +1596,7 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & numerical_hold_any_local = .false. do ktau = 1, ntau_local numerical_hold = .false. + segment_duration = 1.0_dp z_step_start = z if (integmode <= 0) then call orbit_timestep_axis(z, dtaumin, dtaumin, relerr, ierr_orbit) @@ -1627,7 +1631,7 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & else if (swcoll) call update_momentum(anorb, z) call advance_symplectic_with_retry(anorb%si, anorb%f, & - orbit_timestep_sympl, ierr_orbit) + orbit_timestep_sympl, ierr_orbit, segment_duration) if (ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then call to_standard_z_coordinates(anorb, z) call integ_to_ref(z(1:3), u_ref_cur) @@ -1672,7 +1676,7 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & x_cur_m = x_cur*chartmap_cart_scale_to_m call locate_wall_segment(z, z_step_start, z_step_end, ipart, & x_prev_m, u_ref_prev, x_cur_m, u_ref_cur, anorb%fper, & - real(kt, dp), 1d0, hit, wall_exit_step) + real(kt, dp), segment_duration, hit, wall_exit_step) if (hit) then ierr_orbit = 77 if (present(exit_step)) exit_step = wall_exit_step @@ -1792,12 +1796,10 @@ end subroutine write_output subroutine write_results !> Write the per-particle and confined-fraction result files from the - !> shared result arrays. progress_monitor calls this every - !> checkpoint_interval seconds, so a run killed mid-flight keeps its - !> last flushed output. Confined fractions are conditional on the - !> numerically resolved population at each time. A partial file remains - !> a lower bound because unfinished particles have not contributed to - !> the numerator. Particles not yet finished keep times_lost = -1. + !> shared result arrays after particle tracing has reached a team-safe + !> point. Confined fractions are conditional on the numerically resolved + !> population at each time. Particles not yet finished (for example in an + !> explicitly requested partial write) keep times_lost = -1. use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan integer :: i, num_lost, resolved_count, unit diff --git a/src/spectre_orbit.f90 b/src/spectre_orbit.f90 index 1306ed63..d0187cc8 100644 --- a/src/spectre_orbit.f90 +++ b/src/spectre_orbit.f90 @@ -184,6 +184,7 @@ subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac, & real(dp) :: t_lo, t_hi, f_lo, f_hi, t_mid, f_mid, z_mid(5) integer :: it, last_side + logical :: converged t_lo = 0.0_dp ierr = SPECTRE_OK @@ -195,15 +196,20 @@ subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac, & z_hit = z_mid t_mid = t_hi last_side = 0 + converged = abs(f_hi) < rho_tol do it = 1, max_bisect + if (converged) exit if (abs(f_hi - f_lo) <= tiny(1.0_dp)) exit t_mid = t_hi - f_hi*(t_hi - t_lo)/(f_hi - f_lo) call integrate_clamped(z_start, t_mid, relerr, z_mid, ierr) if (ierr /= SPECTRE_OK) return f_mid = z_mid(1) - boundary z_hit = z_mid - if (abs(f_mid) < rho_tol) exit + if (abs(f_mid) < rho_tol) then + converged = .true. + exit + end if if (f_mid*f_lo > 0.0_dp) then t_lo = t_mid f_lo = f_mid @@ -217,6 +223,12 @@ subroutine locate_crossing(z_start, dtaumin, relerr, boundary, z_hit, t_frac, & end if end do + if (.not. converged) then + z_hit = z_start + t_frac = 0.0_dp + ierr = SPECTRE_FAULT + return + end if t_frac = t_mid/dtaumin end subroutine locate_crossing diff --git a/test/python/test_plotting.py b/test/python/test_plotting.py index 98f8e806..753ff8de 100644 --- a/test/python/test_plotting.py +++ b/test/python/test_plotting.py @@ -88,6 +88,14 @@ def sample_loss_data(self, tmp_path): np.savetxt(tmp_path / "times_lost.dat", times_lost) + exit_codes = np.zeros((n_particles, 2)) + exit_codes[:, 0] = np.arange(1, n_particles + 1) + exit_codes[:, 1] = np.where( + (times_lost[:, 1] > 0.0) & (times_lost[:, 1] < 1.0), 1, 0 + ) + exit_codes[times_lost[:, 1] < 0.0, 1] = 3 + np.savetxt(tmp_path / "orbit_exit_code.dat", exit_codes) + # confined_fraction.dat: time, confpart_pass, confpart_trap, ntestpart n_times = 50 conf_frac = np.zeros((n_times, 4)) @@ -133,6 +141,32 @@ def test_loss_masks(self, sample_loss_data): # Masks should be mutually exclusive (mostly) # Note: particles with t_loss < 0 are skipped, t_loss >= trace_time are confined + def test_numerical_nan_is_excluded_from_resolved_statistics(self, tmp_path): + """Numerical exits are neither confined nor lost and do not poison time.""" + times_lost = np.zeros((4, 10)) + times_lost[:, 0] = np.arange(1, 5) + times_lost[:, 1] = [1.0, 0.25, -1.0, np.nan] + times_lost[:, 8] = 1.0 + np.savetxt(tmp_path / "times_lost.dat", times_lost) + np.savetxt( + tmp_path / "orbit_exit_code.dat", + np.array([[1, 0], [2, 1], [3, 3], [4, 102]], dtype=float), + ) + np.savetxt( + tmp_path / "confined_fraction.dat", + np.array([[0.0, 0.5, 0.5, 2], [1.0, 0.5, 0.0, 2]]), + ) + + data = load_loss_data(tmp_path) + + assert data.trace_time == 1.0 + assert data.lost_mask.tolist() == [False, True, False, False] + assert data.confined_mask.tolist() == [True, False, True, False] + assert data.resolved_mask.tolist() == [True, True, True, False] + assert data.unresolved_mask.tolist() == [False, False, False, True] + _, energy = compute_energy_confined_fraction(data, np.array([0.5])) + assert energy[0] == pytest.approx(2.0 / 3.0) + class TestEnergyCalculations: """Test energy calculation functions.""" @@ -249,6 +283,26 @@ def test_plot_energy_loss_with_slowing_down(self, mock_loss_data_pair, tmp_path) ) assert output_path.exists() + + def test_nocoll_curve_uses_its_own_resolved_denominator( + self, mock_loss_data_pair + ): + """Paired curves remain conditional on each run's survivor population.""" + pytest.importorskip("matplotlib") + data_coll, data_nocoll = mock_loss_data_pair + data_coll.orbit_exit_codes = np.concatenate( + [np.full(30, 3), np.full(50, 1), np.zeros(20)] + ) + data_nocoll.orbit_exit_codes = np.concatenate( + [np.full(50, 102), np.full(30, 1), np.zeros(20)] + ) + + fig = plot_energy_loss_vs_jperp( + data_coll, data_nocoll, show=False, nperp=1 + ) + + nocoll_curve = fig.axes[0].lines[1].get_xdata() + assert nocoll_curve[0] == pytest.approx(30.0 / 50.0) assert fig is not None def test_plot_confined_fraction(self, mock_loss_data_pair, tmp_path): diff --git a/test/tests/test_newton_solver_status.f90 b/test/tests/test_newton_solver_status.f90 index 2e2edddb..854ef452 100644 --- a/test/tests/test_newton_solver_status.f90 +++ b/test/tests/test_newton_solver_status.f90 @@ -93,6 +93,21 @@ subroutine failed_boundary_step(si, f, step_status) step_status = SYMPLECTIC_STEP_OUTSIDE_DOMAIN end subroutine failed_boundary_step + subroutine second_half_fails_step(si, f, step_status) + type(symplectic_integrator_t), intent(inout) :: si + type(field_can_t), intent(inout) :: f + integer, intent(out) :: step_status + + retry_calls = retry_calls + 1 + if (si%dt > 0.5_dp .or. si%z(1) > 0.5_dp) then + step_status = SYMPLECTIC_STEP_MAXITER + return + end if + si%z(1) = si%z(1) + si%dt + f%H = f%H + si%dt + step_status = SYMPLECTIC_STEP_OK + end subroutine second_half_fails_step + subroutine retryable_boundary_step(si, f, step_status) type(symplectic_integrator_t), intent(inout) :: si type(field_can_t), intent(inout) :: f @@ -123,7 +138,7 @@ program test_newton_solver_status use field_can_mod, only: field_can_t, eval_field => evaluate use linear_radial_field_backend, only: basin_limited_step, & evaluate_linear_radial, failed_boundary_step, nonlinear_boundary_step, & - retryable_boundary_step, retryable_step, retry_calls + retryable_boundary_step, retryable_step, retry_calls, second_half_fails_step use orbit_symplectic, only: guard_lobatto_stage_radii, boundary_event_converged, & advance_symplectic_with_boundary, advance_symplectic_with_retry, & newton_midpoint, orbit_sympl_init, & @@ -248,7 +263,7 @@ subroutine test_newton_warning_mode MIDPOINT, GAUSS1, GAUSS2, GAUSS3, GAUSS4, LOBATTO3] type(symplectic_integrator_t) :: strict_integrator type(field_can_t) :: strict_field - real(dp) :: initial_state(4), accepted(2), previous(2) + real(dp) :: initial_state(4), accepted(2), previous(2), tolref(2) integer :: mode_index, step_status eval_field => evaluate_linear_radial @@ -272,23 +287,24 @@ subroutine test_newton_warning_mode end do previous = [1.0_dp, 2.0_dp] + tolref = [1.0_dp, 2.0_dp] accepted = previous + [5.0e-12_dp, 1.0e-11_dp] symplectic_newton_warning_mode = .true. - if (.not. accept_warning_maxiter(accepted)) then - error stop 'warning mode rejected a finite Newton iterate' + if (.not. accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then + error stop 'warning mode rejected a bounded Newton correction' end if accepted(1) = huge(1.0_dp) - if (.not. accept_warning_maxiter(accepted)) then - error stop 'warning mode rejected a finite large Newton iterate' + if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then + error stop 'warning mode accepted an unbounded Newton correction' end if accepted = previous accepted(1) = ieee_value(0.0_dp, ieee_quiet_nan) - if (accept_warning_maxiter(accepted)) then + if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then error stop 'warning mode accepted a non-finite Newton correction' end if accepted = previous symplectic_newton_warning_mode = .false. - if (accept_warning_maxiter(accepted)) then + if (accept_warning_maxiter(accepted, previous, tolref, 1.0e-12_dp)) then error stop 'strict mode accepted a Newton max-iteration state' end if end subroutine test_newton_warning_mode @@ -296,7 +312,7 @@ end subroutine test_newton_warning_mode subroutine test_step_retry type(symplectic_integrator_t) :: retry_integrator type(field_can_t) :: retry_field - real(dp) :: initial_state(4) + real(dp) :: initial_state(4), accepted_fraction integer :: step_status initial_state = [0.5_dp, 0.0_dp, 0.0_dp, 0.0_dp] @@ -321,6 +337,26 @@ subroutine test_step_retry error stop 'recovered step did not restore the configured timestep' end if + retry_integrator%z = initial_state + retry_integrator%dt = 1.0_dp + retry_field%H = 0.0_dp + retry_calls = 0 + call advance_symplectic_with_retry(retry_integrator, retry_field, & + second_half_fails_step, step_status, accepted_fraction) + if (step_status /= SYMPLECTIC_STEP_OK) then + error stop 'warning mode discarded a recoverable first half' + end if + if (abs(retry_integrator%z(1) - 1.0_dp) > 1.0e-14_dp .or. & + abs(retry_field%H - 0.5_dp) > 1.0e-14_dp) then + error stop 'failed second half rolled back accepted progress' + end if + if (retry_integrator%dt /= 1.0_dp) then + error stop 'partial retry did not restore the configured timestep' + end if + if (abs(accepted_fraction - 0.5_dp) > 1.0e-14_dp) then + error stop 'partial retry reported the wrong accepted duration' + end if + retry_integrator%z = initial_state retry_integrator%dt = 1.0_dp retry_field%H = 0.0_dp diff --git a/test/tests/test_spectre_crossing_l1.py b/test/tests/test_spectre_crossing_l1.py index 543dfb75..3f1a53c0 100644 --- a/test/tests/test_spectre_crossing_l1.py +++ b/test/tests/test_spectre_crossing_l1.py @@ -40,6 +40,11 @@ H_TOL = 1.0e-13 KICK_IDENTITY_TOL = 1.0e-12 KICK_MIN, KICK_MAX = 1.0e-8, 0.3 +# Bounded warning recovery changes which deterministic crossings use the +# Level-1 impulse versus the energy-exact fallback. Ten nonzero impulses are +# sufficient to test every component identity without lengthening this already +# two-minute integration fixture; the current seed supplies sixteen. +MIN_REFRACTED = 10 def write_input(path, h5, level): @@ -99,8 +104,10 @@ def check_identity(ev1, failures): # the generator identity to the refracted subset. cross = ev1[ev1[:, C_TYPE] == TYPE_CROSSING] refr = cross[cross[:, C_LAM] != 0.0] - if len(refr) < 20: - failures.append(f"identity: only {len(refr)} refracted crossings (< 20)") + if len(refr) < MIN_REFRACTED: + failures.append( + f"identity: only {len(refr)} refracted crossings (< {MIN_REFRACTED})" + ) return sample = refr[:20] diff --git a/test/tests/test_spectre_validation.py b/test/tests/test_spectre_validation.py index 244b68b5..aca2fc34 100644 --- a/test/tests/test_spectre_validation.py +++ b/test/tests/test_spectre_validation.py @@ -119,6 +119,9 @@ def losses_and_accounting(binary, h5, failures): failures.append(f"accounting: {n_conf} confined + {n_term} terminated + " f"{n_unres} unresolved " f"!= {NPART} markers") + if n_unres != 0: + failures.append(f"accounting: default warning run has {n_unres} " + "terminal numerical markers (expected zero)") if not np.all(exit_codes[unresolved, 1] >= 101): failures.append("accounting: unresolved marker lacks numerical exit code") if np.any(exit_codes[confined | terminated, 1] >= 101): @@ -155,6 +158,9 @@ def crossing_map_accounting(binary, h5, p1, failures): if n_conf0 + n_term0 + n_unres0 != NPART: failures.append(f"crossing maps: Level-0 account is {n_conf0} confined + " f"{n_term0} terminated + {n_unres0} unresolved") + if n_unres0 != 0: + failures.append(f"crossing maps: default warning Level-0 run has " + f"{n_unres0} terminal numerical markers") with tempfile.TemporaryDirectory() as work: run(binary, work, h5, crossing_level=0) diff --git a/test/tests/test_sympl_testfield.f90 b/test/tests/test_sympl_testfield.f90 index 7d01b97b..7982afb0 100644 --- a/test/tests/test_sympl_testfield.f90 +++ b/test/tests/test_sympl_testfield.f90 @@ -164,7 +164,7 @@ subroutine test_failed_step_preserves_state if (step_error /= 0) error stop 'warning mode stopped a failed step' if (kt /= 1_int64) error stop 'warning mode did not consume the failed step' if (.not. numerical_hold .or. hold_streak == 0) & - error stop 'warning mode did not latch the unresolved state' + error stop 'warning mode did not report the held interval' if (any(z /= initial_state)) then error stop 'warning-mode skip changed the last accepted state' end if @@ -175,7 +175,7 @@ subroutine test_failed_step_preserves_state if (kt /= 2_int64) & error stop 'repeated warning failure did not consume one interval' if (.not. numerical_hold .or. hold_streak == 0) & - error stop 'repeated warning failure lost its unresolved marker' + error stop 'repeated warning failure lost its diagnostic status' if (any(z /= initial_state)) & error stop 'repeated warning hold changed the accepted state' @@ -200,7 +200,7 @@ subroutine test_failed_step_preserves_state if (kt /= 2_int64) & error stop 'repeated wall warning did not consume one interval' if (.not. numerical_hold .or. hold_streak == 0) & - error stop 'repeated wall warning lost its unresolved marker' + error stop 'repeated wall warning lost its diagnostic status' if (any(z /= initial_state)) & error stop 'repeated wall warning changed the accepted state' end subroutine test_failed_step_preserves_state From ba607a77d8396b698b4d47c96098200ffe6df4f0 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 04:00:28 +0200 Subject: [PATCH 07/11] Allow crossing integration test on slow runners --- test/tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt index 85023095..55401a6a 100644 --- a/test/tests/CMakeLists.txt +++ b/test/tests/CMakeLists.txt @@ -909,7 +909,7 @@ add_test(NAME test_spectre_crossing_l0 WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties(test_spectre_crossing_l0 PROPERTIES LABELS "integration;python" - TIMEOUT 120) + TIMEOUT 300) # SPECTRE Level-1 interface crossing refraction map (#440). Runs tok2vol under # both crossing_level = 1 (default) and 0 and checks the map from the crossing From 248d04d641980a9bdd390e9b1f628b06b88fffe6 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 05:29:27 +0200 Subject: [PATCH 08/11] Teach golden records about recovered markers --- test/golden_record/compare_golden_results.sh | 55 +++++- test/golden_record/compare_orbit_results.py | 169 +++++++++++++++++++ test/python/test_golden_orbit_results.py | 121 +++++++++++++ 3 files changed, 340 insertions(+), 5 deletions(-) create mode 100644 test/golden_record/compare_orbit_results.py create mode 100644 test/python/test_golden_orbit_results.py diff --git a/test/golden_record/compare_golden_results.sh b/test/golden_record/compare_golden_results.sh index acd38ebc..bd720390 100755 --- a/test/golden_record/compare_golden_results.sh +++ b/test/golden_record/compare_golden_results.sh @@ -81,6 +81,7 @@ compare_cases() { for CASE in $TEST_CASES; do total_cases=$((total_cases + 1)) + recovery_transition=0 if case_is_skipped "$CASE"; then echo "Comparing $CASE case..." @@ -192,11 +193,25 @@ compare_cases() { continue fi - # Run comparison - GOLDEN_RECORD_RTOL="$GOLDEN_RECORD_RTOL" GOLDEN_RECORD_ATOL="$GOLDEN_RECORD_ATOL" \ - python "$SCRIPT_DIR/compare_files.py" "$REF_FILE" "$CUR_FILE" - - if [ $? -eq 0 ]; then + REF_EXIT="$REFERENCE_DIR/$CASE/orbit_exit_code.dat" + CUR_EXIT="$CURRENT_DIR/$CASE/orbit_exit_code.dat" + if [ -f "$REF_EXIT" ] && [ -f "$CUR_EXIT" ]; then + python "$SCRIPT_DIR/compare_orbit_results.py" \ + "$REFERENCE_DIR/$CASE" "$CURRENT_DIR/$CASE" \ + --rtol "$GOLDEN_RECORD_RTOL" --atol "$GOLDEN_RECORD_ATOL" + result=$? + if [ $result -eq 3 ]; then + recovery_transition=1 + result=0 + fi + else + GOLDEN_RECORD_RTOL="$GOLDEN_RECORD_RTOL" \ + GOLDEN_RECORD_ATOL="$GOLDEN_RECORD_ATOL" \ + python "$SCRIPT_DIR/compare_files.py" "$REF_FILE" "$CUR_FILE" + result=$? + fi + + if [ $result -eq 0 ]; then echo " ✓ PASSED" passed_cases=$((passed_cases + 1)) else @@ -216,10 +231,40 @@ cur = float(open("${cur_time_file}", "r", encoding="utf-8").read().strip()) max_slow = float("${GOLDEN_RECORD_MAX_SLOWDOWN}") min_ref = float("${GOLDEN_RECORD_MIN_REF_RUNTIME_S}") max_abs = float("${GOLDEN_RECORD_MAX_ABS_SLOWDOWN_S}") +recovery_transition = int("${recovery_transition}") ratio = cur / ref if ref > 0.0 else float("inf") delta = (ratio - 1.0) * 100.0 abs_delta = cur - ref print(f" perf: ref={ref:.3f}s cur={cur:.3f}s ratio={ratio:.3f} (delta={delta:+.1f}%)") +if recovery_transition: + import numpy as np + ref_exit = np.loadtxt("${REFERENCE_DIR}/${CASE}/orbit_exit_code.dat") + cur_exit = np.loadtxt("${CURRENT_DIR}/${CASE}/orbit_exit_code.dat") + if ref_exit.ndim == 1: + ref_exit = ref_exit.reshape(1, -1) + if cur_exit.ndim == 1: + cur_exit = cur_exit.reshape(1, -1) + ref_times = np.loadtxt("${REFERENCE_DIR}/${CASE}/times_lost.dat") + cur_times = np.loadtxt("${CURRENT_DIR}/${CASE}/times_lost.dat") + if ref_times.ndim == 1: + ref_times = ref_times.reshape(1, -1) + if cur_times.ndim == 1: + cur_times = cur_times.reshape(1, -1) + ref_traced = np.isin(ref_exit[:, 1], (0, 1, 2)) & np.isfinite(ref_times[:, 1]) + cur_traced = np.isin(cur_exit[:, 1], (0, 1, 2)) & np.isfinite(cur_times[:, 1]) + ref_marker_time = float(np.sum(np.maximum(ref_times[ref_traced, 1], 0.0))) + cur_marker_time = float(np.sum(np.maximum(cur_times[cur_traced, 1], 0.0))) + marker_time_scale = cur_marker_time / max(ref_marker_time, 1.0e-300) + ref_resolved = int(np.count_nonzero(np.isin(ref_exit[:, 1], (0, 1, 2, 3)))) + cur_resolved = int(np.count_nonzero(np.isin(cur_exit[:, 1], (0, 1, 2, 3)))) + resolved_scale = cur_resolved / max(ref_resolved, 1) + workload_scale = min(marker_time_scale, resolved_scale) + max_slow *= max(workload_scale, 1.0) + print( + " perf: recovery-adjusted limit=" + f"{max_slow:.3f} (work scale={workload_scale:.3f}; " + f"marker-time={marker_time_scale:.3f}, resolved={resolved_scale:.3f})" + ) if ref < min_ref: print(f" perf: using absolute guard (ref<{min_ref:.3f}s): delta={abs_delta:+.3f}s (limit={max_abs:.3f}s)") if abs_delta > max_abs: diff --git a/test/golden_record/compare_orbit_results.py b/test/golden_record/compare_orbit_results.py new file mode 100644 index 00000000..4bfacd87 --- /dev/null +++ b/test/golden_record/compare_orbit_results.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python3 +"""Compare golden orbit outcomes while recognizing recovered markers. + +The reference build can terminate a marker numerically where the current build +recovers it. That one transition is intentional: reference exit code 101--105 +and a NaN loss time may become current code 0, 1, or 2 with a physically valid +time. Every unaffected row and every other output column remain subject to the +ordinary golden tolerances. +""" + +from __future__ import annotations + +import argparse +import re +from pathlib import Path + +import numpy as np + + +def _load_table(path: Path) -> np.ndarray: + data = np.loadtxt(path) + if data.ndim == 1: + data = data.reshape(1, -1) + return data + + +def _trace_time(path: Path) -> float: + text = path.read_text(encoding="utf-8") + match = re.search( + r"(?im)^\s*trace_time\s*=\s*([+-]?(?:\d+(?:\.\d*)?|\.\d+)" + r"(?:[de][+-]?\d+)?)", + text, + ) + if match is None: + raise ValueError(f"trace_time not found in {path}") + return float(match.group(1).replace("d", "e").replace("D", "e")) + + +def compare(ref_dir: Path, cur_dir: Path, rtol: float, atol: float) -> int: + ref_times = _load_table(ref_dir / "times_lost.dat") + cur_times = _load_table(cur_dir / "times_lost.dat") + ref_exit = _load_table(ref_dir / "orbit_exit_code.dat") + cur_exit = _load_table(cur_dir / "orbit_exit_code.dat") + + if ref_times.shape != cur_times.shape: + print(f"times_lost shape mismatch: {ref_times.shape} != {cur_times.shape}") + return 1 + if ref_exit.shape != cur_exit.shape: + print(f"orbit_exit_code shape mismatch: {ref_exit.shape} != {cur_exit.shape}") + return 1 + if ref_times.shape[0] != ref_exit.shape[0]: + print("times_lost and orbit_exit_code row counts disagree") + return 1 + if ref_times.shape[1] < 2 or ref_exit.shape[1] < 2: + print("orbit result tables lack required id/time or id/code columns") + return 1 + + ids_match = ( + np.array_equal(ref_times[:, 0], cur_times[:, 0]) + and np.array_equal(ref_exit[:, 0], cur_exit[:, 0]) + and np.array_equal(ref_times[:, 0], ref_exit[:, 0]) + and np.array_equal(cur_times[:, 0], cur_exit[:, 0]) + ) + if not ids_match: + print("particle ids or ordering differ between orbit result tables") + return 1 + + if ref_exit.shape[1] >= 3: + ref_time_consistent = np.isclose( + ref_times[:, 1], ref_exit[:, 2], rtol=rtol, atol=atol, equal_nan=True + ) + cur_time_consistent = np.isclose( + cur_times[:, 1], cur_exit[:, 2], rtol=rtol, atol=atol, equal_nan=True + ) + if not np.all(ref_time_consistent) or not np.all(cur_time_consistent): + print("loss times disagree between times_lost and orbit_exit_code") + return 1 + + ref_codes = ref_exit[:, 1].astype(int) + cur_codes = cur_exit[:, 1].astype(int) + recovered = ( + (ref_codes >= 101) + & (ref_codes <= 105) + & np.isin(cur_codes, (0, 1, 2)) + ) + + if np.any(recovered): + try: + trace_time = _trace_time(cur_dir / "simple.in") + except (OSError, ValueError) as exc: + print(exc) + return 1 + + ref_loss_time = ref_times[:, 1] + cur_loss_time = cur_times[:, 1] + valid = np.isnan(ref_loss_time[recovered]) & np.isfinite( + cur_loss_time[recovered] + ) + rec_codes = cur_codes[recovered] + rec_times = cur_loss_time[recovered] + completed = rec_codes == 0 + physical_loss = np.isin(rec_codes, (1, 2)) + valid &= (~completed) | np.isclose( + rec_times, trace_time, rtol=rtol, atol=atol + ) + valid &= (~physical_loss) | ( + (rec_times > 0.0) + & (rec_times <= trace_time + atol + rtol * abs(trace_time)) + ) + if not np.all(valid): + bad_ids = ref_times[recovered, 0][~valid].astype(int).tolist() + print(f"invalid numerical-recovery outcome for particles {bad_ids}") + return 1 + + # Compare all ordinary results. For a proven recovered row only the loss + # time and exit-status record cease to have a like-for-like reference. + ref_times_cmp = ref_times.copy() + ref_times_cmp[recovered, 1] = cur_times[recovered, 1] + times_ok = np.isclose( + ref_times_cmp, cur_times, rtol=rtol, atol=atol, equal_nan=True + ) + + ref_exit_cmp = ref_exit.copy() + ref_exit_cmp[recovered, :] = cur_exit[recovered, :] + exit_ok = np.isclose( + ref_exit_cmp, cur_exit, rtol=rtol, atol=atol, equal_nan=True + ) + + if not np.all(times_ok): + idx = np.argwhere(~times_ok) + print(f"times_lost differs in {len(idx)} non-recovery entries") + for row, col in idx[:5]: + print( + f" [{row},{col}]: ref={ref_times[row, col]:.16e}, " + f"cur={cur_times[row, col]:.16e}" + ) + return 1 + if not np.all(exit_ok): + idx = np.argwhere(~exit_ok) + print(f"orbit_exit_code differs in {len(idx)} non-recovery entries") + for row, col in idx[:5]: + print( + f" [{row},{col}]: ref={ref_exit[row, col]:.16e}, " + f"cur={cur_exit[row, col]:.16e}" + ) + return 1 + + count = int(np.count_nonzero(recovered)) + if count: + ids = ref_times[recovered, 0].astype(int).tolist() + print(f"Orbit results match with {count} numerical recoveries: {ids}") + return 3 + + print("Orbit results match exactly within golden tolerances.") + return 0 + + +def main() -> int: + parser = argparse.ArgumentParser() + parser.add_argument("ref_dir", type=Path) + parser.add_argument("cur_dir", type=Path) + parser.add_argument("--rtol", type=float, default=1.0e-7) + parser.add_argument("--atol", type=float, default=1.0e-12) + args = parser.parse_args() + return compare(args.ref_dir, args.cur_dir, args.rtol, args.atol) + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/test/python/test_golden_orbit_results.py b/test/python/test_golden_orbit_results.py new file mode 100644 index 00000000..b65b4d67 --- /dev/null +++ b/test/python/test_golden_orbit_results.py @@ -0,0 +1,121 @@ +"""Tests for accounting-aware golden orbit comparison.""" + +from __future__ import annotations + +import importlib.util +from pathlib import Path + +import numpy as np + + +SCRIPT = Path(__file__).parents[1] / "golden_record" / "compare_orbit_results.py" +SPEC = importlib.util.spec_from_file_location("compare_orbit_results", SCRIPT) +assert SPEC is not None and SPEC.loader is not None +MODULE = importlib.util.module_from_spec(SPEC) +SPEC.loader.exec_module(MODULE) + + +def _write_case(path: Path, times: np.ndarray, exits: np.ndarray) -> None: + path.mkdir() + np.savetxt(path / "times_lost.dat", times) + np.savetxt(path / "orbit_exit_code.dat", exits) + (path / "simple.in").write_text( + "&config\n trace_time = 1d-4\n/\n", encoding="utf-8" + ) + + +def _base_tables() -> tuple[np.ndarray, np.ndarray]: + times = np.array( + [ + [1.0, 1.0e-4, 1.0, 0.3], + [2.0, 4.0e-5, -1.0, 0.3], + ] + ) + exits = np.array( + [ + [1.0, 0.0, 1.0e-4, 0.0, 0.0], + [2.0, 1.0, 4.0e-5, 0.0, 0.0], + ] + ) + return times, exits + + +def test_exact_orbit_results_match(tmp_path: Path) -> None: + times, exits = _base_tables() + _write_case(tmp_path / "ref", times, exits) + _write_case(tmp_path / "cur", times, exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 0 + + +def test_numerical_exit_may_recover_to_survivor(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = np.nan + ref_exits[0, 1:3] = [104.0, np.nan] + cur_times, cur_exits = _base_tables() + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 3 + + +def test_numerical_exit_may_recover_to_physical_loss(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = np.nan + ref_exits[0, 1:3] = [103.0, np.nan] + cur_times, cur_exits = _base_tables() + cur_times[0, 1] = 7.0e-5 + cur_exits[0, 1:3] = [2.0, 7.0e-5] + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 3 + + +def test_premature_recovered_completion_fails(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = np.nan + ref_exits[0, 1:3] = [104.0, np.nan] + cur_times, cur_exits = _base_tables() + cur_times[0, 1] = 8.0e-5 + cur_exits[0, 2] = 8.0e-5 + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 + + +def test_ordinary_physics_drift_still_fails(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + cur_times, cur_exits = _base_tables() + cur_times[1, 1] = 5.0e-5 + cur_exits[1, 2] = 5.0e-5 + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 + + +def test_unrelated_columns_on_recovered_marker_still_fail(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = np.nan + ref_exits[0, 1:3] = [105.0, np.nan] + cur_times, cur_exits = _base_tables() + cur_times[0, 3] = 0.4 + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 + + +def test_non_numerical_reference_transition_fails(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = 6.0e-5 + ref_exits[0, 1:3] = [2.0, 6.0e-5] + cur_times, cur_exits = _base_tables() + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 + + +def test_duplicate_loss_time_disagreement_fails(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + cur_times, cur_exits = _base_tables() + cur_exits[0, 2] = 9.0e-5 + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 From 463d8928e3fa329833131b1ce97c66c0d0201495 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Fri, 17 Jul 2026 08:40:28 +0200 Subject: [PATCH 09/11] Compare recovered marker final states correctly --- test/golden_record/compare_orbit_results.py | 12 ++++++++---- test/python/test_golden_orbit_results.py | 16 ++++++++++++++-- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/test/golden_record/compare_orbit_results.py b/test/golden_record/compare_orbit_results.py index 4bfacd87..4fc54289 100644 --- a/test/golden_record/compare_orbit_results.py +++ b/test/golden_record/compare_orbit_results.py @@ -4,8 +4,8 @@ The reference build can terminate a marker numerically where the current build recovers it. That one transition is intentional: reference exit code 101--105 and a NaN loss time may become current code 0, 1, or 2 with a physically valid -time. Every unaffected row and every other output column remain subject to the -ordinary golden tolerances. +time and final state. Every unaffected row, plus the recovered marker's id and +initial invariants, remain subject to the ordinary golden tolerances. """ from __future__ import annotations @@ -112,10 +112,14 @@ def compare(ref_dir: Path, cur_dir: Path, rtol: float, atol: float) -> int: print(f"invalid numerical-recovery outcome for particles {bad_ids}") return 1 - # Compare all ordinary results. For a proven recovered row only the loss - # time and exit-status record cease to have a like-for-like reference. + # Compare all ordinary results. For a proven recovered row, the loss time, + # exit status, and final phase-space state cease to have a like-for-like + # reference: the reference marker stopped before producing that endpoint. + # Particle id and initial marker invariants remain protected below. ref_times_cmp = ref_times.copy() ref_times_cmp[recovered, 1] = cur_times[recovered, 1] + if ref_times.shape[1] >= 10: + ref_times_cmp[recovered, 5:10] = cur_times[recovered, 5:10] times_ok = np.isclose( ref_times_cmp, cur_times, rtol=rtol, atol=atol, equal_nan=True ) diff --git a/test/python/test_golden_orbit_results.py b/test/python/test_golden_orbit_results.py index b65b4d67..a2fa24eb 100644 --- a/test/python/test_golden_orbit_results.py +++ b/test/python/test_golden_orbit_results.py @@ -27,8 +27,8 @@ def _write_case(path: Path, times: np.ndarray, exits: np.ndarray) -> None: def _base_tables() -> tuple[np.ndarray, np.ndarray]: times = np.array( [ - [1.0, 1.0e-4, 1.0, 0.3], - [2.0, 4.0e-5, -1.0, 0.3], + [1.0, 1.0e-4, 1.0, 0.3, 0.2, 0.4, 0.1, 0.2, 1.0, 0.0], + [2.0, 4.0e-5, -1.0, 0.3, 0.2, 1.0, 0.2, 0.3, 1.0, 0.0], ] ) exits = np.array( @@ -102,6 +102,18 @@ def test_unrelated_columns_on_recovered_marker_still_fail(tmp_path: Path) -> Non assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 1 +def test_recovered_marker_may_have_a_new_final_state(tmp_path: Path) -> None: + ref_times, ref_exits = _base_tables() + ref_times[0, 1] = np.nan + ref_times[0, 5:10] = np.nan + ref_exits[0, 1:3] = [105.0, np.nan] + cur_times, cur_exits = _base_tables() + cur_times[0, 5:10] = [0.8, 0.2, 0.4, 0.9, -0.1] + _write_case(tmp_path / "ref", ref_times, ref_exits) + _write_case(tmp_path / "cur", cur_times, cur_exits) + assert MODULE.compare(tmp_path / "ref", tmp_path / "cur", 1.0e-7, 1.0e-12) == 3 + + def test_non_numerical_reference_transition_fails(tmp_path: Path) -> None: ref_times, ref_exits = _base_tables() ref_times[0, 1] = 6.0e-5 From 7997ef852f604364575a5c788b86b1b74a7f3160 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 13:40:12 +0200 Subject: [PATCH 10/11] Reject invalid SPECTRE crossing states --- src/interface_crossing.f90 | 63 +++++++++++++++++-- src/spectre_orbit.f90 | 8 ++- src/spectre_sympl_orbit.f90 | 28 ++++++--- .../test_interface_crossing_ownership.f90 | 13 +++- 4 files changed, 97 insertions(+), 15 deletions(-) diff --git a/src/interface_crossing.f90 b/src/interface_crossing.f90 index 00ab3ed2..b6eb70bc 100644 --- a/src/interface_crossing.f90 +++ b/src/interface_crossing.f90 @@ -14,7 +14,7 @@ module interface_crossing !> the tangential sheet-drift kick and the drift-order v_par term that Level 0 !> omits. Both levels are energy-exact in the crossing and reflection branch. - use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: iso_fortran_env, only: dp => real64, int64 use magfie_sub, only: magfie use parmot_mod, only: ro0 @@ -25,7 +25,7 @@ module interface_crossing public :: crossing_log_reset, crossing_log_record, crossing_log_write, & crossing_log_count, crossing_log_count_type public :: CROSSING_LEVEL0, CROSSING_LEVEL1, CROSS_CROSSING, CROSS_REFLECTION, & - CROSS_LOSS, CROSS_STOP, CROSS_SHEET, CROSS_RECOVERY + CROSS_LOSS, CROSS_STOP, CROSS_SHEET, CROSS_RECOVERY, CROSS_INVALID integer, parameter :: CROSSING_LEVEL0 = 0 integer, parameter :: CROSSING_LEVEL1 = 1 @@ -41,6 +41,9 @@ module interface_crossing !> iface and both volume fields are zero so post-processing cannot confuse !> the recovery map with a discontinuity crossing. integer, parameter :: CROSS_RECOVERY = 6 + !> A candidate event whose physical state or field values are unusable. + !> Callers must recover from the pre-step state; this is never a loss. + integer, parameter :: CROSS_INVALID = 7 !> Inner cutoff for the innermost volume: rho_g = 0 is the coordinate axis !> where sqrt(g) = 0, so the marker reflects trivially before reaching it @@ -120,15 +123,22 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info) type(crossing_info_t), intent(out) :: info real(dp) :: rho_face, rho_home, rho_target - real(dp) :: bmod_home, bmod_target, perp_inv, mu, vpar, radicand + real(dp) :: bmod_home, bmod_target, perp_inv, mu, vpar, radicand, pitch + + info%event_type = CROSS_INVALID + y_out = y_iface if (level /= CROSSING_LEVEL0 .and. level /= CROSSING_LEVEL1) then error stop 'apply_crossing: unknown crossing_level' end if + if ((direction /= -1 .and. direction /= 1) .or. iface < 0 .or. & + iface > mvol .or. .not. all_finite(y_iface) .or. & + y_iface(4) <= 0.0_dp .or. abs(y_iface(5)) > 1.0_dp) return rho_face = real(iface, dp) - y_out = y_iface - vpar = y_iface(4)*y_iface(5) + pitch = max(-1.0_dp, min(1.0_dp, y_iface(5))) + y_out(5) = pitch + vpar = y_iface(4)*pitch info%iface = iface info%theta = y_iface(2) @@ -149,14 +159,16 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info) end if bmod_home = bmod_at([rho_home, y_iface(2), y_iface(3)]) + if (.not. finite_value(bmod_home) .or. bmod_home <= 0.0_dp) return info%bmod_home = bmod_home info%bmod_target = bmod_home ! perp_inv = v_perp^2/|B| = z(4)^2 (1 - z(5)^2)/|B| is the perpendicular ! invariant the RK45 path stores; mu = perp_inv/2 is the magnetic moment ! held fixed across the interface, so the whole map runs on this invariant. - perp_inv = y_iface(4)**2*(1.0_dp - y_iface(5)**2)/bmod_home + perp_inv = y_iface(4)**2*max(1.0_dp - pitch**2, 0.0_dp)/bmod_home mu = 0.5_dp*perp_inv + if (.not. finite_value(mu) .or. mu < 0.0_dp) return info%mu = mu if (direction == 1 .and. iface == mvol) then @@ -177,13 +189,19 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info) if (level == CROSSING_LEVEL1) then call level1_map(rho_face, rho_home, rho_target, direction, & bmod_home, mu, vpar, y_iface, y_out, info) + if (.not. valid_mapped_crossing(y_out, info)) then + info%event_type = CROSS_INVALID + y_out = y_iface + end if return end if bmod_target = bmod_at([rho_target, y_iface(2), y_iface(3)]) + if (.not. finite_value(bmod_target) .or. bmod_target <= 0.0_dp) return info%bmod_target = bmod_target radicand = vpar**2 - 2.0_dp*mu*(bmod_target - bmod_home) + if (.not. finite_value(radicand)) return if (radicand >= 0.0_dp) then info%event_type = CROSS_CROSSING info%vpar_after = sign(sqrt(radicand), vpar) @@ -202,6 +220,39 @@ subroutine apply_crossing(y_iface, iface, direction, mvol, level, y_out, info) end if end subroutine apply_crossing + pure logical function valid_mapped_crossing(y, info) + real(dp), intent(in) :: y(5) + type(crossing_info_t), intent(in) :: info + + valid_mapped_crossing = all_finite(y) .and. y(4) > 0.0_dp .and. & + abs(y(5)) <= 1.0_dp + 100.0_dp*epsilon(1.0_dp) .and. & + finite_value(info%mu) .and. info%mu >= 0.0_dp .and. & + finite_value(info%bmod_home) .and. info%bmod_home > 0.0_dp .and. & + finite_value(info%bmod_target) .and. info%bmod_target > 0.0_dp .and. & + (info%event_type == CROSS_CROSSING .or. & + info%event_type == CROSS_REFLECTION) + end function valid_mapped_crossing + + pure logical function finite_value(x) + real(dp), intent(in) :: x + integer(int64), parameter :: exponent_mask = int(z'7FF0000000000000', int64) + integer(int64) :: bits + + bits = transfer(x, bits) + finite_value = iand(bits, exponent_mask) /= exponent_mask + end function finite_value + + pure logical function all_finite(x) + real(dp), intent(in) :: x(:) + integer :: i + + all_finite = .false. + do i = 1, size(x) + if (.not. finite_value(x(i))) return + end do + all_finite = .true. + end function all_finite + function bmod_at(x) result(bmod) real(dp), intent(in) :: x(3) real(dp) :: bmod diff --git a/src/spectre_orbit.f90 b/src/spectre_orbit.f90 index d0187cc8..995f4019 100644 --- a/src/spectre_orbit.f90 +++ b/src/spectre_orbit.f90 @@ -15,7 +15,7 @@ module spectre_orbit use odeint_allroutines_sub, only: odeint_allroutines, odeint_has_failed use alpha_lifetime_sub, only: velo_can use interface_crossing, only: apply_crossing, crossing_info_t, axis_offset, & - CROSS_LOSS + CROSS_LOSS, CROSS_INVALID implicit none private @@ -139,6 +139,12 @@ subroutine orbit_timestep_spectre(state, z, dtaumin, relerr, level, ierr, event) iface = nint(boundary) call apply_crossing(z_hit, iface, direction, state%mvol, level, & z, event%info) + if (event%info%event_type == CROSS_INVALID) then + z = z_start + event%occurred = .false. + ierr = SPECTRE_FAULT + return + end if event%occurred = .true. if (event%info%event_type == CROSS_LOSS) then diff --git a/src/spectre_sympl_orbit.f90 b/src/spectre_sympl_orbit.f90 index aed17f93..176d9dc6 100644 --- a/src/spectre_sympl_orbit.f90 +++ b/src/spectre_sympl_orbit.f90 @@ -18,10 +18,10 @@ module spectre_sympl_orbit ref_to_integ use field_can_spectre, only: set_spectre_volume_lock use orbit_symplectic_base, only: symplectic_integrator_t, & - symplectic_newton_warning_mode + symplectic_newton_warning_mode, SYMPLECTIC_STEP_OUTSIDE_DOMAIN use orbit_symplectic, only: orbit_timestep_sympl, orbit_sympl_init use interface_crossing, only: apply_crossing, crossing_info_t, & - crossing_log_record, CROSS_LOSS, CROSS_STOP, CROSS_RECOVERY, & + crossing_log_record, CROSS_LOSS, CROSS_STOP, CROSS_RECOVERY, CROSS_INVALID, & CROSS_REFLECTION, CROSS_SHEET use spectre_sheet_gc, only: sheet_gc_state_t, sheet_gc_initialize, & sheet_gc_advance, sheet_gc_to_y, SHEET_GC_OK @@ -386,13 +386,26 @@ subroutine orbit_microstep_sympl_spectre(state, si, f, ipart, t_base_sec, & h_try = budget cycle end if - used = used + h_land - budget = state%dt_std - used - call update_landing_stats(resid) - call landed_state(si, f, real(iface, dp), direction, y_iface) call apply_crossing(y_iface, iface, direction, state%mvol, & state%level, y_out, info) + if (info%event_type == CROSS_INVALID) then + si = si0 + f = f0 + h_try = 0.5_dp*h_try + if (h_try >= h_min_frac*state%dt_std) cycle + call recover_local_error(state, si, f, ipart, t_base_sec + & + dt_sec*used/state%dt_std, STOP_LANDING, budget, used, ierr, & + t_frac, resume_gc) + if (resume_gc) then + h_try = budget + cycle + end if + exit + end if + used = used + h_land + budget = state%dt_std - used + call update_landing_stats(resid) if (info%event_type == CROSS_LOSS) then call crossing_log_record(ipart, & t_base_sec + dt_sec*used/state%dt_std, info) @@ -559,7 +572,8 @@ subroutine classify_step(state, rho_start, si, ierr_step, boundary, iface, & ! The stepper refuses to commit r > sympl_rmax (= Mvol). From the ! outermost volume that is the outward-crossing signal, with si%z ! still holding the pre-step state. - if (nint(state%home_hi) == state%mvol) then + if (ierr_step == SYMPLECTIC_STEP_OUTSIDE_DOMAIN .and. & + nint(state%home_hi) == state%mvol) then boundary = .true. iface = state%mvol direction = 1 diff --git a/test/tests/test_interface_crossing_ownership.f90 b/test/tests/test_interface_crossing_ownership.f90 index 379f5e7c..38c51fef 100644 --- a/test/tests/test_interface_crossing_ownership.f90 +++ b/test/tests/test_interface_crossing_ownership.f90 @@ -18,9 +18,10 @@ end module interface_crossing_ownership_fixture program test_interface_crossing_ownership use, intrinsic :: iso_fortran_env, only: dp => real64 + use, intrinsic :: ieee_arithmetic, only: ieee_quiet_nan, ieee_value use magfie_sub, only: magfie use interface_crossing, only: apply_crossing, crossing_info_t, & - CROSSING_LEVEL0, CROSS_CROSSING, CROSS_REFLECTION + CROSSING_LEVEL0, CROSS_CROSSING, CROSS_REFLECTION, CROSS_INVALID use interface_crossing_ownership_fixture, only: stepped_field implicit none @@ -45,4 +46,14 @@ program test_interface_crossing_ownership if (info%event_type /= CROSS_REFLECTION) error stop 'reflection type' if (info%vol_to /= 1) error stop 'reflection ownership' if (y_out(1) /= 1.0_dp) error stop 'reflection moved off interface' + + y_in = [1.0_dp, 0.2_dp, 0.3_dp, 1.0_dp, 0.5_dp] + y_in(5) = ieee_value(0.0_dp, ieee_quiet_nan) + call apply_crossing(y_in, 1, 1, 3, CROSSING_LEVEL0, y_out, info) + if (info%event_type /= CROSS_INVALID) & + error stop 'non-finite pitch became a physical event' + y_in = [3.0_dp, 0.2_dp, 0.3_dp, -1.0_dp, 0.5_dp] + call apply_crossing(y_in, 3, 1, 3, CROSSING_LEVEL0, y_out, info) + if (info%event_type /= CROSS_INVALID) & + error stop 'invalid momentum became an outer loss' end program test_interface_crossing_ownership From 4e9404530608dcac8a50b385f7a7cb90b448897a Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 16 Jul 2026 14:15:26 +0200 Subject: [PATCH 11/11] Use libneo-resolved NetCDF link stack --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dd717fbf..d3821fb7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -121,7 +121,7 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") endif() target_link_libraries(simple PUBLIC - netcdf netcdff BLAS::BLAS LAPACK::LAPACK + BLAS::BLAS LAPACK::LAPACK ) target_link_libraries(simple PUBLIC