From cff201eb521ee9fa6959baa3630bc210958d9a8f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 15 Jul 2026 14:25:27 +0200 Subject: [PATCH] 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