Skip to content

Commit 9eca50e

Browse files
krystophnykrystophny
authored
Check physical momentum before accepting symplectic endpoints (#560)
[V] Restore [SIMPLE](https://github.com/itpplasma/SIMPLE)'s existing 0.5% momentum bound using a fresh endpoint field and canonical parallel momentum, including localized ordinary/wall losses. Preserve the carried cache, collision reference and electric work. Main 9a6697e is merged. [V] Golden changes are intended physics corrections; no new production defect was found. Per-marker causes: | Case | Markers | Cause | |---|---|---| | Albert coils | 1, 2, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 19, 21, 22, 24, 25, 26, 29, 31, 32 | #560 rejects the first fresh physical endpoint beyond 0.5%; recovery then changes the trajectory. | | Meiss coils | 1, 2, 6, 7, 9, 10, 11, 12, 13, 15, 16, 17, 19, 21, 22, 24, 25, 26, 29, 31, 32 | Same corrected physical rejection. | | Meiss coils | 8 | #540 recovery-start normalization changes only reported p/p0: 0.999899952176 to 1.0. | | Meiss coils | 18 | Same normalization: 1.000186277524 to 1.0; time, position, pitch and exit status agree exactly. | [V] Albert marker 9 returns from main's p/p0=3.28614485936 to 1.0. Its differences from pre-#540 columns 2, 6, 7, 8, 10 follow recovery after the genuine physical rejection at microstep 87 (physical drift -0.504850499%). The old cache rejected step 1 despite physical drift only +0.003256568%. [V] Independent reconstruction uses v_parallel=(Pphi-Aphi/ro0)/hphi and p²=mu*B+v_parallel²/2; the unchanged bound gives energy interval [0.990025,1.010025]. Four retained canonical/field samples independently violate that interval. The new parent diagnostic preserves the raw table byte-for-byte. For Meiss 8/18, reconstructing boundary fractions 0.331670804921 / 0.580034170284 from loss time and microstep width reproduces both old reported momenta exactly via p_start+fraction*(p_RK-p_start). The #540 recovery start at src/simple_main.f90:2755 makes both interpolation momenta equal to the reference. [V] Refresh the two coil orbit records from the deterministic GCC 12.2 run at 048c02f, libneo 924ba7735366d95c2c6855bf00bd0e46b8c05cf2. Apply them only to exact historical references 9ec87fd / 9a6697e; retain raw reference runs and original field diagnostics. Comparators and tolerances are unchanged. Other raw differences are completed survivor phases (Albert/Meiss 2,13, protected by the existing outcome contract) and pre-#540 collision trajectories (the #558 stream/operator correction, checked by the native analytic collision contract). [V] Native validation with system GCC/GFortran 12.2, deterministic FP64, CGAL and the Python/IMAS interface: endpoint regression 1/1, fast selection 162/162 (sCluster job 3300299 COMPLETED 0:0), and slow selection 45/45 (CTest results from 3300192). Golden job 3300118 COMPLETED 0:0: both 15-test selections passed against historical parent 9ec87fd and main 9a6697e. It executed 6e124d7; production/test sources and all snapshot numerical tokens are identical at final head f1eb9ec. [V] The initial slow-job failure was external setup: --no-build-isolation lacked wheel, so function-generator metadata failed before compilation. The replacement installed its build prerequisites. A separate inherited GOLDEN_RECORD_PERF_GATE=0 conflicted with the comparator's performance-gate test; job 3300277 reproduced this on main and passed all 35 checks on main plus head's registered test with the override removed. The full corrected fast selection then passed; no comparator, tolerance, test or production source was weakened. All 45 slow tests passed before 3300192's wrapper error, which was independently reproduced as replacing an active Bash script; the final fast submission kept its script unchanged. [V] Current-head hosted CI is green, including both build/test shards, changed-line lint, CGAL, documentation and the required Build and Test / CI result checks. Merge through the configured queue. Chris&AI
2 parents ead8d0a + 710fd62 commit 9eca50e

7 files changed

Lines changed: 389 additions & 32 deletions

File tree

‎src/simple_main.f90‎

Lines changed: 73 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@ module simple_main
4141
ORBIT_EXIT_NUMERICAL_DOMAIN, &
4242
ORBIT_EXIT_NUMERICAL_MAXITER, ORBIT_EXIT_NUMERICAL_LINEAR, &
4343
ORBIT_EXIT_NUMERICAL_EVENT, ORBIT_EXIT_NUMERICAL_FULL_ORBIT
44+
use params, only: radial_electric_potential_slope
4445
use params, only: canonical_grid_nr, canonical_grid_ntheta, &
4546
canonical_grid_nphi, canonical_ode_relerr
4647
use diag_counters, only: diag_counters_init, count_event, &
@@ -2529,6 +2530,56 @@ subroutine validate_rk_step_or_reference_boundary(z_start, z_candidate, &
25292530
call validate_rk_state(z_start, z, ierr, expected_momentum)
25302531
end subroutine validate_rk_step_or_reference_boundary
25312532

2533+
subroutine validate_symplectic_endpoint(si_start, anorb, z_start, &
2534+
z_candidate, z, event_fraction, ierr, reference_boundary)
2535+
use field_can_mod, only: get_val
2536+
2537+
type(symplectic_integrator_t), intent(in) :: si_start
2538+
type(tracer_t), intent(in) :: anorb
2539+
real(dp), intent(in) :: z_start(5), z_candidate(5)
2540+
real(dp), intent(inout) :: z(5)
2541+
real(dp), intent(out) :: event_fraction
2542+
integer, intent(out) :: ierr
2543+
logical, intent(out) :: reference_boundary
2544+
type(field_can_t) :: endpoint_field
2545+
real(dp) :: z_physical(5), pphi, momentum_reference
2546+
integer :: momentum_ierr
2547+
2548+
ierr = 0
2549+
call validate_rk_step_or_reference_boundary(z_start, z_candidate, &
2550+
anorb%fper, z=z, event_fraction=event_fraction, ierr=ierr, &
2551+
reference_boundary=reference_boundary)
2552+
if (ierr /= 0) then
2553+
if (.not. reference_boundary .or. ierr /= 1) return
2554+
end if
2555+
2556+
! Euler's carried field is extrapolated before the angular update; it
2557+
! is not the field at si%z. Validate a copy at the accepted endpoint,
2558+
! without changing the next map's carried field or public coordinates.
2559+
endpoint_field = anorb%f
2560+
pphi = anorb%si%z(4)
2561+
if (reference_boundary) pphi = si_start%z(4) + &
2562+
event_fraction*(pphi - si_start%z(4))
2563+
call evaluate_field(endpoint_field, z(1), z(2), z(3), 0)
2564+
call get_val(endpoint_field, pphi)
2565+
z_physical = z
2566+
z_physical(4) = sqrt(endpoint_field%mu*endpoint_field%Bmod + &
2567+
0.5_dp*endpoint_field%vpar**2)
2568+
z_physical(5) = endpoint_field%vpar/(sqrt2*z_physical(4))
2569+
momentum_reference = si_start%pabs
2570+
! Kinetic momentum is invariant only without electric work.
2571+
if (abs(radial_electric_potential_slope) > 0.0_dp) &
2572+
momentum_reference = z_physical(4)
2573+
momentum_ierr = 0
2574+
call validate_rk_state(z_start, z_physical, momentum_ierr, &
2575+
momentum_reference)
2576+
if (momentum_ierr /= 0) then
2577+
z = z_start
2578+
ierr = momentum_ierr
2579+
reference_boundary = .false.
2580+
end if
2581+
end subroutine validate_symplectic_endpoint
2582+
25322583
subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, &
25332584
hold_streak, numerical_hold_any, rk_recovery)
25342585
use alpha_lifetime_sub, only: orbit_timestep_axis
@@ -2546,7 +2597,7 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, &
25462597
logical, intent(out), optional :: numerical_hold_any
25472598
type(rk_recovery_state_t), intent(inout), optional :: rk_recovery
25482599

2549-
integer :: hold_streak_local, ktau, recovery_method
2600+
integer :: hold_streak_local, ktau, recovery_method, symplectic_status
25502601
integer :: ierr_coll, collision_hold_start
25512602
real(dp) :: z_step_start(5), z_step_end(5), loss_fraction
25522603
real(dp) :: u_event(3)
@@ -2659,16 +2710,15 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, &
26592710
f_step_start = anorb%f
26602711
call advance_symplectic_with_retry(anorb%si, anorb%f, &
26612712
orbit_timestep_sympl, ierr_orbit)
2662-
if (ierr_orbit == 0) then
2713+
if (ierr_orbit == 0 .or. &
2714+
ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then
2715+
symplectic_status = ierr_orbit
2716+
ierr_orbit = 0
26632717
call to_standard_z_coordinates(anorb, z)
26642718
z_step_end = z
2665-
! z(4) is reconstructed from H, which the map conserves
2666-
! only to O(dt^order), already within one step; the RK
2667-
! momentum-drift gate does not apply to symplectic steps.
2668-
call validate_rk_step_or_reference_boundary(z_step_start, &
2669-
z_step_end, anorb%fper, z=z, &
2670-
event_fraction=loss_fraction, ierr=ierr_orbit, &
2671-
reference_boundary=reference_boundary)
2719+
call validate_symplectic_endpoint(si_step_start, anorb, &
2720+
z_step_start, z_step_end, z, loss_fraction, ierr_orbit, &
2721+
reference_boundary)
26722722
if (reference_boundary .and. ierr_orbit == 1) then
26732723
call integ_to_ref(z(1:3), u_event)
26742724
anorb%si%last_step_fraction = loss_fraction
@@ -2677,8 +2727,11 @@ subroutine macrostep(anorb, z, kt, ierr_orbit, ntau_local, exit_step, &
26772727
anorb%si%last_event_fraction_width = 1.d0
26782728
ierr_orbit = SYMPLECTIC_STEP_BOUNDARY
26792729
else if (ierr_orbit /= 0) then
2730+
ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED
26802731
anorb%si = si_step_start
26812732
anorb%f = f_step_start
2733+
else if (symplectic_status == SYMPLECTIC_STEP_BOUNDARY) then
2734+
ierr_orbit = SYMPLECTIC_STEP_BOUNDARY
26822735
end if
26832736
end if
26842737
if (ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then
@@ -2903,7 +2956,7 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, &
29032956
logical, intent(out), optional :: numerical_hold_any
29042957
type(rk_recovery_state_t), intent(inout), optional :: rk_recovery
29052958

2906-
integer :: hold_streak_local, ktau, recovery_method
2959+
integer :: hold_streak_local, ktau, recovery_method, symplectic_status
29072960
integer :: ierr_coll, collision_hold_start
29082961
real(dp) :: u_ref_prev(3), u_ref_cur(3), x_cur(3), x_cur_m(3)
29092962
real(dp) :: z_step_start(5), z_step_end(5), segment_duration
@@ -2999,23 +3052,25 @@ subroutine macrostep_with_wall_check(anorb, z, kt, ierr_orbit, ntau_local, &
29993052
f_step_start = anorb%f
30003053
call advance_symplectic_with_retry(anorb%si, anorb%f, &
30013054
orbit_timestep_sympl, ierr_orbit, segment_duration)
3002-
if (ierr_orbit == 0) then
3055+
if (ierr_orbit == 0 .or. &
3056+
ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then
3057+
symplectic_status = ierr_orbit
3058+
ierr_orbit = 0
30033059
call to_standard_z_coordinates(anorb, z)
30043060
z_step_end = z
3005-
! z(4) is reconstructed from H, which the map conserves
3006-
! only to O(dt^order), already within one step; the RK
3007-
! momentum-drift gate does not apply to symplectic steps.
3008-
call validate_rk_step_or_reference_boundary(z_step_start, &
3009-
z_step_end, anorb%fper, z=z, &
3010-
event_fraction=boundary_fraction, ierr=ierr_orbit, &
3011-
reference_boundary=reference_boundary)
3061+
call validate_symplectic_endpoint(si_step_start, anorb, &
3062+
z_step_start, z_step_end, z, boundary_fraction, ierr_orbit, &
3063+
reference_boundary)
30123064
if (reference_boundary .and. ierr_orbit == 1) then
30133065
anorb%si%last_step_fraction = boundary_fraction
30143066
anorb%si%last_event_fraction_width = 1.d0
30153067
ierr_orbit = SYMPLECTIC_STEP_BOUNDARY
30163068
else if (ierr_orbit /= 0) then
3069+
ierr_orbit = SYMPLECTIC_STEP_EVENT_NOT_CONVERGED
30173070
anorb%si = si_step_start
30183071
anorb%f = f_step_start
3072+
else if (symplectic_status == SYMPLECTIC_STEP_BOUNDARY) then
3073+
ierr_orbit = SYMPLECTIC_STEP_BOUNDARY
30193074
end if
30203075
end if
30213076
if (ierr_orbit == SYMPLECTIC_STEP_BOUNDARY) then

‎test/golden_record/baselines/physical_momentum/albert_coils/orbit_exit_code.dat‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
1 1 8.6722971322861804E-006 0.0000000000000000 1.0001000100010002E-008
2+
2 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
3+
3 1 5.7023945654970240E-006 0.0000000000000000 1.0001000100010002E-008
4+
4 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
5+
5 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
6+
6 1 6.2725620305795632E-006 3.9968028886505635E-015 1.0001000100010002E-008
7+
7 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
8+
8 1 2.9988756623702986E-005 3.5527136788005009E-015 1.0001000100010002E-008
9+
9 1 2.8691433990049079E-006 0.0000000000000000 1.0001000100010002E-008
10+
10 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
11+
11 1 1.3496564494001541E-005 4.4408920985006262E-016 1.0001000100010002E-008
12+
12 1 7.5238042189822716E-006 0.0000000000000000 1.0001000100010002E-008
13+
13 1 6.7107312757327157E-006 0.0000000000000000 1.0001000100010002E-008
14+
14 1 8.1559519991643968E-006 0.0000000000000000 1.0001000100010002E-008
15+
15 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
16+
16 1 8.5916975335624126E-006 0.0000000000000000 1.0001000100010002E-008
17+
17 1 2.9888792487099021E-006 0.0000000000000000 1.0001000100010002E-008
18+
18 1 1.4561789102277897E-005 0.0000000000000000 1.0001000100010002E-008
19+
19 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
20+
20 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
21+
21 1 6.9948888470367058E-006 4.4408920985006262E-016 1.0001000100010002E-008
22+
22 1 1.0630958713861306E-005 0.0000000000000000 1.0001000100010002E-008
23+
23 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
24+
24 1 4.2553388549307421E-006 5.7731597280508140E-015 1.0001000100010002E-008
25+
25 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
26+
26 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
27+
27 1 1.0191044132675798E-005 0.0000000000000000 1.0001000100010002E-008
28+
28 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
29+
29 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
30+
30 1 5.7854018273426140E-005 5.4178883601707639E-013 1.0001000100010002E-008
31+
31 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000
32+
32 0 1.0000000000000000E-004 -1.0000000000000000 -1.0000000000000000

‎test/golden_record/baselines/physical_momentum/albert_coils/times_lost.dat‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
1 8.6722971322861804E-006 -4.2014945993859225 0.29999999999999999 2.7118938344299529E-005 1.0000000000000000 4.6534842403114993 0.85396613388812526 1.0000000000000000 -0.71744844123740792
2+
2 1.0000000000000000E-004 -4.6363481730088729 0.29999999999999999 1.1770432502225849E-005 0.17820708030194976 1.6311719733831387 3.6850033344292736 1.0000000000000000 -0.90198694344370234
3+
3 5.7023945654970240E-006 -3.4399485485034105 0.29999999999999999 5.3998317335245138E-005 1.0000000000000000 1.7043680550992477 4.5948232003392953 0.99735089681928579 0.31748450527699351
4+
4 1.0000000000000000E-004 -3.1163840562327016 0.29999999999999999 6.5418785743791920E-005 0.20566240749746117 5.7913124855800673 0.72454833965536025 1.0002798675185098 7.8419491631383467E-002
5+
5 1.0000000000000000E-004 -3.4055110577411241 0.29999999999999999 5.5213816207078333E-005 0.41838573730162332 3.1340204918070795 3.2630965618568508 0.99965184848492972 -0.28325826712458518
6+
6 6.2725620305795632E-006 -3.5913533407714646 0.29999999999999999 4.8654363692911486E-005 1.0000000000000040 4.8042160688140960 4.9943620228879091 1.0000000000000000 -0.39663764231944842
7+
7 1.0000000000000000E-004 -4.4282644039274972 0.29999999999999999 1.9114916114899219E-005 0.41810330529404877 2.2007372279497153 4.0191322250553823 1.0000000000000000 -0.81963867202436202
8+
8 2.9988756623702986E-005 -3.2524984572991502 0.29999999999999999 6.0614518729143487E-005 1.0000000000000036 1.5182350661031161 5.8592337090843358 1.0000000000000000 0.22610774046376622
9+
9 2.8691433990049079E-006 -3.8430613759850543 0.29999999999999999 3.9770126599525060E-005 0.13342231420895095 1.6557238875486129 3.9108287520303167 1.0000000000000000 0.59618938305034364
10+
10 1.0000000000000000E-004 -4.1012354071183772 0.29999999999999999 3.0657666991322872E-005 0.10948897717803756 3.2814987890929115 0.72399520612366075 1.0000000000000000 0.88218112873856991
11+
11 1.3496564494001541E-005 -3.6420115743970425 0.29999999999999999 4.6866340684905832E-005 1.0000000000000004 1.6865156804135983 2.9766362933871724 1.0022752021525394 0.48072103176794401
12+
12 7.5238042189822716E-006 -3.4979862715710159 0.29999999999999999 5.1949829325494603E-005 1.0000000000000000 4.8821276630845452 0.99667471664085472 1.0000000000000000 -0.36989441262359357
13+
13 6.7107312757327157E-006 -4.8250948594608918 0.29999999999999999 5.1084667385668663E-006 1.0000000000000000 1.5909228860835176 3.1525751045286778 1.0000000000000000 0.95754861217480369
14+
14 8.1559519991643968E-006 -3.2090570027439020 0.29999999999999999 6.2147819728545025E-005 1.0000000000000000 1.4014541245863690 4.9374221184456317 1.0004334076364942 6.9424717813672002E-002
15+
15 1.0000000000000000E-004 -4.2047637807166440 0.29999999999999999 2.7003549965780902E-005 0.38182192150000616 2.9463853883534030 4.1897777510386920 1.0000000000000000 -0.76460143267103009
16+
16 8.5916975335624126E-006 -4.2014649905207531 0.29999999999999999 2.7119983412956180E-005 1.0000000000000000 1.8425081680049908 0.43755651582321509 0.97622676869163938 0.72132732608023720
17+
17 2.9888792487099021E-006 -3.8843410627805053 0.29999999999999999 3.8313126927979734E-005 1.0000000000000000 1.3786125848553714 2.9526442137829036 1.0000000000000000 0.62171385436570892
18+
18 1.4561789102277897E-005 -3.2091695779845830 0.29999999999999999 6.2143846295085881E-005 1.0000000000000000 1.2913937884935587 3.7023019274517850 0.99953946804615401 0.26474633090486827
19+
19 1.0000000000000000E-004 -4.3562921369401710 0.29999999999999999 2.1655235034481313E-005 0.28592200759576875 3.7430330712544366 0.96492209087753089 1.0000000000000000 -0.85578206787984434
20+
20 1.0000000000000000E-004 -3.0985573200886440 0.29999999999999999 6.6047994701585652E-005 0.41598018049812407 0.70021264899934899 1.2535617397624268 0.99995075546868184 0.12369152782680505
21+
21 6.9948888470367058E-006 -3.7041892531179879 0.29999999999999999 4.4671729618141369E-005 1.0000000000000004 4.4896324837245034 1.5468153985419990 1.0000000000000000 -0.51122269055724123
22+
22 1.0630958713861306E-005 -3.8546169950433642 0.29999999999999999 3.9362261751715314E-005 1.0000000000000000 4.4985723528833734 5.5799881860474070 1.0000000000000000 -0.58798862054622580
23+
23 1.0000000000000000E-004 -3.5831318346407666 0.29999999999999999 4.8944548349384106E-005 0.38572240229178412 4.2220406886427551 3.6703315477766933 1.0014216569168928 -0.44352704383453179
24+
24 4.2553388549307421E-006 -4.3050872348558826 0.29999999999999999 2.3462553144427330E-005 1.0000000000000058 1.3480031627751119 6.1947736296455620 1.0000000000000000 0.79627847571982879
25+
25 1.0000000000000000E-004 -4.3073159196848021 0.29999999999999999 2.3383889924916852E-005 0.36423608494058668 4.2413693647852506 5.1460818922775813 1.0000000000000000 0.85511686983113144
26+
26 1.0000000000000000E-004 -4.1323492741726531 0.29999999999999999 2.9559478085403629E-005 0.17430221970859322 4.4129270731255543 1.5078653014972447 1.0000000000000000 -0.70418700126114731
27+
27 1.0191044132675798E-005 -3.4723688513570949 0.29999999999999999 5.2854016729397923E-005 1.0000000000000000 1.5017367686682528 2.5109287565801797 0.99781515387667952 0.37638590192188137
28+
28 1.0000000000000000E-004 -3.1968250162142544 0.29999999999999999 6.2579557509026613E-005 0.66900030962483636 5.2942180702637689 6.1520874491008222 1.0008945007122552 -0.23503249930433351
29+
29 1.0000000000000000E-004 -4.7117884444523286 0.29999999999999999 9.1077075823624190E-006 0.35703926288183591 3.7205284901236833 4.3228716661039508 1.0000000000000000 0.92515949169617218
30+
30 5.7854018273426140E-005 -3.0926695318993032 0.29999999999999999 6.6255808910557402E-005 1.0000000000005418 1.2503181602910205 1.5056961802970772 1.0014565977036876 -0.17290621357566818
31+
31 1.0000000000000000E-004 -4.2679326302766505 0.29999999999999999 2.4773954728857561E-005 6.0887962203227335E-002 1.9814823319511152 3.0828863173220880 1.0000000000000000 0.86570225191795491
32+
32 1.0000000000000000E-004 -3.7154489712289691 0.29999999999999999 4.4274308831722006E-005 0.31816538583906917 2.8607282036238129 1.9240043611578264 1.0034293878819840 -0.48883430001461842

0 commit comments

Comments
 (0)