[math-02] Derive kinetic response and Fokker-Planck susceptibilities

KIM/tests/CMakeLists.txt

@@ -147,6 +147,22 @@ add_test(NAME test_fp_parallel_flow
 set_tests_properties(test_fp_parallel_flow PROPERTIES
                      WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/tests/)
 
+# Independent arbitrary-precision Mathematica oracle for every FP
+# susceptibility consumed by KIM.
+add_executable(test_fp_susceptibility_oracle
+               ${CMAKE_SOURCE_DIR}/KIM/tests/test_fp_susceptibility_oracle.f90)
+set_target_properties(test_fp_susceptibility_oracle PROPERTIES
+                        OUTPUT_NAME test_fp_susceptibility_oracle.x
+                        RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/")
+target_link_libraries(test_fp_susceptibility_oracle
+                      KIM_lib kilca_lib lapack cerf ddeabm slatec)
+configure_file(${CMAKE_SOURCE_DIR}/verification/oracles/fp_susceptibilities.dat
+               ${CMAKE_BINARY_DIR}/tests/fp_susceptibilities.dat COPYONLY)
+add_test(NAME test_fp_susceptibility_oracle
+         COMMAND ${CMAKE_BINARY_DIR}/tests/test_fp_susceptibility_oracle.x)
+set_tests_properties(test_fp_susceptibility_oracle PROPERTIES
+                     WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/tests/)
+
 # Assemble the dense periodic Fourier matrices K_{m,m'} (Phase-2.4).
 add_executable(test_periodic_assembly ${CMAKE_SOURCE_DIR}/KIM/tests/test_periodic_assembly.f90)
 set_target_properties(test_periodic_assembly PROPERTIES

KIM/tests/test_fp_susceptibility_oracle.f90

@@ -0,0 +1,111 @@
+program test_fp_susceptibility_oracle
+    use KIM_kinds_m, only: dp
+
+    implicit none
+
+    integer, parameter :: nmmax = 3
+    real(dp), parameter :: map_tol = 5.0e-14_dp
+    real(dp), parameter :: oracle_tol = 1.0e-12_dp
+    complex(dp) :: symbI(0:nmmax, 0:nmmax), expected
+    real(dp) :: kpar, temperature_erg, mass, nu, omega_E, Vpar
+    real(dp) :: omega_c, omega, x1_ref, x2_ref, expected_re, expected_im
+    real(dp) :: vT, x1, x2, error_scale, scaled_error, max_scaled_error
+    integer :: case_id, mphi, k, l, iunit, ios, nrows
+    character(len=1024) :: line
+    logical :: seen_pair(0:nmmax, 0:nmmax), seen_harmonic(-1:1)
+
+    interface
+        subroutine getIfunc(x1, x2, symbI)
+            double precision, intent(in) :: x1, x2
+            double complex, intent(out) :: symbI(0:3, 0:3)
+        end subroutine getIfunc
+    end interface
+
+    seen_pair = .false.
+    seen_harmonic = .false.
+    nrows = 0
+    max_scaled_error = 0.0_dp
+    open(newunit=iunit, file='fp_susceptibilities.dat', status='old', &
+         action='read', iostat=ios)
+    if (ios /= 0) error stop 'cannot open fp_susceptibilities.dat'
+
+    do
+        read(iunit, '(A)', iostat=ios) line
+        if (ios < 0) exit
+        if (ios > 0) error stop 'cannot read susceptibility oracle line'
+        if (len_trim(line) == 0 .or. line(1:1) == '#') cycle
+        read(line, *, iostat=ios) case_id, kpar, temperature_erg, mass, nu, &
+            omega_E, Vpar, mphi, omega_c, omega, k, l, x1_ref, x2_ref, &
+            expected_re, expected_im
+        if (ios /= 0) error stop 'malformed susceptibility oracle row'
+        if (k < 0 .or. k > nmmax .or. l < 0 .or. l > nmmax) then
+            error stop 'oracle susceptibility index out of range'
+        end if
+        if (mphi < -1 .or. mphi > 1) error stop 'oracle harmonic out of range'
+
+        vT = sqrt(temperature_erg/mass)
+        x1 = kpar*vT/nu
+        x2 = -(omega_E + kpar*Vpar + real(mphi, dp)*omega_c - omega)/nu
+        call assert_close('x1 mapping', x1, x1_ref, map_tol, case_id, k, l)
+        call assert_close('x2 mapping', x2, x2_ref, map_tol, case_id, k, l)
+
+        call getIfunc(x1, x2, symbI)
+        expected = cmplx(expected_re, expected_im, kind=dp)
+        error_scale = max(1.0_dp, abs(expected))
+        scaled_error = abs(symbI(k, l) - expected)/error_scale
+        max_scaled_error = max(max_scaled_error, scaled_error)
+        if (scaled_error > oracle_tol) then
+            print *, 'FAIL: susceptibility oracle mismatch'
+            print *, ' case, (k,l), mphi = ', case_id, k, l, mphi
+            print *, ' x1, x2 = ', x1, x2
+            print *, ' actual   = ', symbI(k, l)
+            print *, ' expected = ', expected
+            print *, ' scaled error = ', scaled_error
+            error stop
+        end if
+
+        nrows = nrows + 1
+        seen_pair(k, l) = .true.
+        seen_harmonic(mphi) = .true.
+    end do
+    close(iunit)
+
+    if (nrows < 56) error stop 'susceptibility oracle is missing regime coverage'
+    call require_pair(0, 0)
+    call require_pair(2, 0)
+    call require_pair(0, 2)
+    call require_pair(0, 1)
+    call require_pair(2, 1)
+    call require_pair(2, 2)
+    call require_pair(1, 1)
+    call require_pair(1, 3)
+    if (.not. all(seen_harmonic)) error stop 'oracle lacks mphi=-1,0,+1'
+
+    print *, 'PASS: ', nrows, ' independent FP susceptibility oracle rows'
+    print *, 'Maximum scaled oracle error: ', max_scaled_error
+    print *, 'All tests PASSED'
+    stop 0
+
+contains
+
+    subroutine assert_close(label, actual, reference, tolerance, icase, ik, il)
+        character(len=*), intent(in) :: label
+        real(dp), intent(in) :: actual, reference, tolerance
+        integer, intent(in) :: icase, ik, il
+
+        if (abs(actual - reference) > tolerance*max(1.0_dp, abs(reference))) then
+            print *, 'FAIL: ', trim(label), ' case/(k,l) = ', icase, ik, il
+            print *, ' actual/reference = ', actual, reference
+            error stop
+        end if
+    end subroutine assert_close
+
+    subroutine require_pair(ik, il)
+        integer, intent(in) :: ik, il
+        if (.not. seen_pair(ik, il)) then
+            print *, 'FAIL: missing susceptibility pair ', ik, il
+            error stop
+        end if
+    end subroutine require_pair
+
+end program test_fp_susceptibility_oracle

Makefile

@@ -14,7 +14,9 @@ ifeq ($(origin LIBNEO_PATH),command line)
 endif
 
 .PHONY: all ninja test clean KIM KiLCA QL-Balance PreProc install install-kim ctest golden pytest \
-	verify-mathematica-inventory verify-mathematica-inventory-negative
+	verify-mathematica-inventory verify-mathematica-inventory-negative \
+	verify-mathematica-susceptibilities verify-mathematica-susceptibilities-negative \
+	regenerate-fp-susceptibility-oracle
 
 all: ninja
 
@@ -57,6 +59,23 @@ verify-mathematica-inventory-negative:
 		fi; \
 	done
 
+verify-mathematica-susceptibilities:
+	@test -n "$(WOLFRAM_KERNEL)" || { echo "WolframKernel not found; set WOLFRAM_KERNEL"; exit 1; }
+	$(WOLFRAM_KERNEL) -noinit -noprompt -script verification/mathematica/02_kinetic_and_susceptibilities.wl
+
+verify-mathematica-susceptibilities-negative:
+	@for fixture in wrong-branch thermal-speed flow-sign recurrence-index; do \
+		if KAMEL_VERIFY_NEGATIVE_FIXTURE=$$fixture $(MAKE) --no-print-directory verify-mathematica-susceptibilities >/dev/null 2>&1; then \
+			echo "negative fixture unexpectedly passed: $$fixture"; exit 1; \
+		else \
+			echo "negative fixture rejected: $$fixture"; \
+		fi; \
+	done
+
+regenerate-fp-susceptibility-oracle:
+	@test -n "$(WOLFRAM_KERNEL)" || { echo "WolframKernel not found; set WOLFRAM_KERNEL"; exit 1; }
+	KAMEL_REGENERATE_ORACLE=1 $(WOLFRAM_KERNEL) -noinit -noprompt -script verification/mathematica/02_kinetic_and_susceptibilities.wl
+
 # Golden-record regression now lives in test/golden/ and runs in the dedicated
 # GitHub Actions job, NOT in `make test`/ctest. This target is for manual local
 # runs only; it does the A/B double build and needs the golden-baseline tag.

verification/FORMULA_INDEX.md

@@ -559,7 +559,7 @@ path exists and that each code family has a source-to-check chain.
 | Site | Code symbol | Class | Thesis mapping | Verification |
 | --- | --- | --- | --- | --- |
 | KIM-01 | `KIM/src/background_equilibrium/species_mod.f90::calculate_plasma_backs` | vT, gyrofrequency, Larmor/Debye lengths, collision scale | (4.14)–(4.16), (5.8) | conventions script; #198 pending |
-| KIM-02 | `KIM/src/background_equilibrium/species_mod.f90::calculate_thermodynamic_forces_and_susc` | A1/A2, x1/x2, FP susceptibilities | (5.8), (5.12)–(5.20) | flow test; #194 pending |
+| KIM-02 | `KIM/src/background_equilibrium/species_mod.f90::calculate_thermodynamic_forces_and_susc` | A1/A2, x1/x2, FP susceptibilities | (5.8), (5.12)–(5.20) | flow test; independent #194 oracle |
 | KIM-03 | `KIM/src/background_equilibrium/calculate_equil.f90::calculate_equil` | h_theta/h_z, k_s, k_parallel, omega_E | (3.5), (5.18)–(5.20) | conventions script; #198 pending |
 | KIM-04 | `KIM/src/asymptotics/flr2_fourier_kernel.f90` | b_plus, b_cross, phases, charge/current kernels | Ch. 14 | #196 pending |
 | KIM-05 | `KIM/src/kernels/FP_kernel_plasma_prefacs.f90` | finite-radius FP kernel prefactors | Ch. 14 | #196 pending |
@@ -570,24 +570,37 @@ path exists and that each code family has a source-to-check chain.
 | KIM-10 | `KIM/src/grid/prepare_resonances.f90` | q_res=abs(m/n) | (3.5)–(3.6) | conventions script |
 | KIM-11 | `KIM/src/background_equilibrium/profile_input_m.f90` | radial force balance and Vz-to-Vparallel projection | (8.8)–(8.9) | profile/flow tests |
 | KIM-12 | `KIM/src/background_equilibrium/periodic_background.f90` | periodic primitive state and derivative reconstruction | Ch. 3, 8 | #198 pending |
-| KILCA-01 | `KiLCA/flre/conductivity/calc_I_array.f90::calc_Imn_array` | susceptibility-function definition | (5.12)–(5.17), App. A | #194 pending |
-| KILCA-02 | `KiLCA/solver/VER_5_STABLE/wave_stuff.f90` | plasma/cyclotron frequencies and dielectric response | Ch. 4–5 | #194/#197 pending |
-| KILCA-03 | `KiLCA/flre/conductivity/calc_I_array_drift_serg.f90` | drift/FP conductivity integrals | Ch. 5, App. A | #194 pending |
+| KILCA-01 | `KiLCA/flre/conductivity/calc_I_array.f90::calc_Imn_array` | susceptibility-function definition | (5.12)–(5.17), App. A | not exercised by KIM #194 oracle |
+| KILCA-02 | `KiLCA/solver/VER_5_STABLE/wave_stuff.f90` | plasma/cyclotron frequencies and dielectric response | Ch. 4–5 | not exercised by KIM #194 oracle; #197 pending |
+| KILCA-03 | `KiLCA/flre/conductivity/calc_I_array_drift_serg.f90` | drift/FP conductivity integrals | Ch. 5, App. A | not exercised by KIM #194 oracle |
 | KILCA-04 | `PreProc/fourier/src/rhs_flt.f90` | exp(-i m iota phi) field-line Fourier phase | (4.3) | conventions script |
-| QLB-01 | `QL-Balance/src/base/getIfunc.f90` | susceptibility functions | App. A | #194 pending |
+| QLB-01 | `QL-Balance/src/base/getIfunc.f90` | susceptibility functions | App. A | independent #194 oracle |
 | QLB-02 | `QL-Balance/src/base/get_dql.f90` | quasilinear transport coefficients | (6.11)–(6.16) | broader #199 |
 | QLB-03 | `QL-Balance/src/base/rhs_balance_m.f90::compute_particle_fluxes` | particle fluxes | (6.2), (6.11)–(6.16) | unit tests; broader #199 |
 | QLB-04 | `QL-Balance/src/base/rhs_balance_m.f90::compute_total_heat_fluxes` | heat fluxes | (6.9), (6.11)–(6.16) | unit tests; broader #199 |
 | QLB-05 | `QL-Balance/src/base/calc_current_densities.f90` | current moments and thermal speeds | (4.8), (5.8) | broader #199 |
+| QLB-06 | `QL-Balance/src/base/W2_arr.f90::W2_arr` | production differential FP moments consumed by KIM | (5.14)–(5.15), (A.1)–(A.3) | independent #194 oracle |
 | EQ-01 | `common/equil/mag_wrapper.f90::magfie` | cylindrical metric and field unit vector | Ch. 3 | #198 pending |
 | EQ-02 | `common/equil/equil_profiles.f90` | flux-surface equilibrium profiles | Ch. 3, 8 | #198 pending |
 | PRE-01 | `PreProc/fourier/src/fouriermodes.f90` | straight-field-line angle and equilibrium Fourier modes | (4.3), Ch. 7–8 | #198 pending |
 <!-- CODE-SITES-END -->
 
 ## Coverage policy
 
+### Fokker–Planck susceptibility oracle (#194)
+
+`verification/mathematica/02_kinetic_and_susceptibilities.wl` independently
+evaluates the generating integral (5.14)–(5.15), applies the energy-conserving
+rank-one correction (5.13), and traces `x1/x2` through (4.1), (4.26),
+(5.16)–(5.20), (2.52)–(2.54), and Appendix (A.1)–(A.29). The generated
+`verification/oracles/fp_susceptibilities.dat` covers every `I00`, `I20`,
+`I02`, `I01`, `I21`, `I22`, `I11`, and `I13` consumed by KIM, including
+finite flow, `m_phi=-1,0,+1`, near-resonance, strong-collision, and
+collisionless-scaled points. `test_fp_susceptibility_oracle` reads this file
+without requiring Mathematica and compares it to the production implementation.
+
 The inventory script requires all 494 thesis labels, all 19 convention rows,
-and all 24 semantic code sites. It rejects a missing equation, missing code
+and all 25 semantic code sites. It rejects a missing equation, missing code
 site/path, incompatible duplicate convention, or unclassified convention.
 Subsequent oracle issues replace “pending” statuses with deterministic
 script/test links; copying a formula here never upgrades its status.

verification/mathematica/01_conventions_and_inventory.wl

@@ -42,7 +42,7 @@ check["formula index exists", FileExistsQ[indexPath]];
 check["494 unique thesis equations indexed", Length[DeleteDuplicates[indexedEquations]] == 494];
 check["complete thesis equation labels", Sort[indexedEquations] === Sort[expectedEquations]];
 check["no duplicate thesis equation rows", DuplicateFreeQ[indexedEquations]];
-check["24 semantic code formula sites classified", Length[codeSiteRows] == 24];
+check["25 semantic code formula sites classified", Length[codeSiteRows] == 25];
 check["19 conventions classified", Length[conventionRows] == 19];
 
 codePaths = {
@@ -62,6 +62,7 @@ codePaths = {
   "KiLCA/flre/conductivity/calc_I_array_drift_serg.f90",
   "PreProc/fourier/src/rhs_flt.f90",
   "QL-Balance/src/base/getIfunc.f90",
+  "QL-Balance/src/base/W2_arr.f90",
   "QL-Balance/src/base/get_dql.f90",
   "QL-Balance/src/base/rhs_balance_m.f90",
   "QL-Balance/src/base/calc_current_densities.f90",