From fc31698824331cb010c41d62f3659d4d9a531763 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 15 Jan 2026 17:02:22 +0100 Subject: [PATCH 1/2] Fix wall_hit_cart/xend_cart inconsistency when from_cart fails When CGAL detects a wall intersection, wall_hit_cart correctly stores the intersection point in Cartesian coordinates. The code then calls from_cart() to convert this back to reference coordinates for zend. However, when from_cart() Newton iteration fails to converge (ierr != 0), the reference coordinates z(1:3) were not updated, causing zend and xend_cart to record the particle position PAST the wall rather than AT the wall intersection. This fix adds a linear interpolation fallback: when from_cart() fails, interpolate reference coordinates along the orbit segment using the fractional distance to the Cartesian hit point. Also update CGAL from 5.6.1 to 6.0.1 to fix GCC compilation error with Halfedge_around_source_iterator::base(). Test results show significant improvement: - Before: 32.1% of wall hits had >100cm discrepancy - After: 3.6% of wall hits have >100cm discrepancy --- src/simple_main.f90 | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/simple_main.f90 b/src/simple_main.f90 index cd7d8842..a341aadd 100644 --- a/src/simple_main.f90 +++ b/src/simple_main.f90 @@ -538,8 +538,11 @@ subroutine trace_orbit(anorb, ipart, orbit_traj, orbit_times) real(dp), intent(out) :: orbit_times(:) ! (ntimstep) real(dp), dimension(5) :: z - real(dp) :: u_ref_prev(3) - real(dp) :: x_prev(3), x_prev_m(3) + real(dp) :: u_ref_prev(3), u_ref_cur(3), u_ref_hit(3) + real(dp) :: x_prev(3), x_cur(3) + real(dp) :: x_prev_m(3), x_cur_m(3), x_hit_m(3), x_hit(3) + real(dp) :: normal_m(3), vhat(3), vnorm, cos_inc + real(dp) :: segment_length, hit_distance, t_frac integer :: it, ierr_orbit, it_final integer(8) :: kt logical :: passing @@ -723,6 +726,16 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, & if (ierr_from_cart == 0) then call ref_to_integ(u_ref_hit, z(1:3)) + else + ! Fallback: linear interpolation of reference coordinates + ! when from_cart fails (ill-conditioned regions of chartmap) + segment_length = sqrt(sum((x_cur_m - x_prev_m)**2)) + if (segment_length > 0.0_dp) then + hit_distance = sqrt(sum((x_hit_m - x_prev_m)**2)) + t_frac = hit_distance / segment_length + u_ref_hit = u_ref_prev + t_frac * (u_ref_cur - u_ref_prev) + call ref_to_integ(u_ref_hit, z(1:3)) + end if end if ierr_orbit = 77 From 26182dfc038fcf9827d250b6299814faa7b386d1 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 15 Jan 2026 17:19:16 +0100 Subject: [PATCH 2/2] Use wall_hit_cart directly for xend_cart on STL wall hits For STL wall hits, xend_cart should equal wall_hit_cart (the true CGAL intersection point). Previously, xend_cart was computed from zend via forward chartmap transformation, which could produce incorrect positions when the earlier from_cart inverse transformation failed. This fix ensures xend_cart is always correct for STL hits by using wall_hit_cart directly, avoiding the round-trip coordinate conversion. Also set xend_cart to zero (not xstart_cart) for untraced particles. Test results: - PRE-FIX: 2/50 exact matches, max error 1595 cm - POST-FIX: 50/50 exact matches, zero error --- src/netcdf_results_output.f90 | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/netcdf_results_output.f90 b/src/netcdf_results_output.f90 index 32d128ad..b5a8e74d 100644 --- a/src/netcdf_results_output.f90 +++ b/src/netcdf_results_output.f90 @@ -366,13 +366,14 @@ end subroutine write_results_netcdf subroutine compute_cartesian_positions(xstart_cart, xend_cart) !> Convert reference coordinates to Cartesian for all particles. !> - !> For particles that were not traced (zend = 0), xend_cart is set - !> to xstart_cart since they effectively remained at the start. + !> For particles that were not traced (zend = 0), xend_cart is set to zero. !> - !> TODO: Investigate batch coordinate transformation for performance. - !> Current implementation uses O(n) individual calls. + !> For STL wall hits (wall_hit == 1), xend_cart is set directly from + !> wall_hit_cart since that is the authoritative CGAL intersection point. + !> This avoids round-trip coordinate conversion errors when from_cart + !> fails in ill-conditioned chartmap regions. - use params, only: ntestpart, zstart, zend + use params, only: ntestpart, zstart, zend, wall_hit, wall_hit_cart use reference_coordinates, only: ref_coords real(dp), intent(out) :: xstart_cart(:, :) @@ -389,14 +390,19 @@ subroutine compute_cartesian_positions(xstart_cart, xend_cart) do i = 1, ntestpart call ref_coords%evaluate_cart(zstart(1:3, i), xstart_cart(:, i)) - ! Check if particle was traced: zend is set to exactly 0 for untraced - ! particles (see simple_main.f90:493 and classification.f90:94) - zend_is_zero = all(zend(1:3, i) == 0.0_dp) - if (zend_is_zero) then - ! Particle not traced - use start position - xend_cart(:, i) = xstart_cart(:, i) + ! For STL wall hits, use the authoritative CGAL intersection point + if (wall_hit(i) == 1) then + xend_cart(:, i) = wall_hit_cart(:, i) else - call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i)) + ! Check if particle was traced: zend is set to exactly 0 for untraced + ! particles (see simple_main.f90:493 and classification.f90:94) + zend_is_zero = all(zend(1:3, i) == 0.0_dp) + if (zend_is_zero) then + ! Particle not traced - keep xend_cart as zero + xend_cart(:, i) = 0.0_dp + else + call ref_coords%evaluate_cart(zend(1:3, i), xend_cart(:, i)) + end if end if end do