diff --git a/KIM/tests/CMakeLists.txt b/KIM/tests/CMakeLists.txt index 6501935c..83dac984 100644 --- a/KIM/tests/CMakeLists.txt +++ b/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 diff --git a/KIM/tests/test_fp_susceptibility_oracle.f90 b/KIM/tests/test_fp_susceptibility_oracle.f90 new file mode 100644 index 00000000..cd60b010 --- /dev/null +++ b/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 diff --git a/Makefile b/Makefile index 835e9cac..4ecad8d8 100644 --- a/Makefile +++ b/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. diff --git a/verification/FORMULA_INDEX.md b/verification/FORMULA_INDEX.md index 39c3894a..788fb0c9 100644 --- a/verification/FORMULA_INDEX.md +++ b/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,15 +570,16 @@ 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 | @@ -586,8 +587,20 @@ path exists and that each code family has a source-to-check chain. ## 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. diff --git a/verification/mathematica/01_conventions_and_inventory.wl b/verification/mathematica/01_conventions_and_inventory.wl index 87b4c254..6c9786d2 100644 --- a/verification/mathematica/01_conventions_and_inventory.wl +++ b/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", diff --git a/verification/mathematica/02_kinetic_and_susceptibilities.wl b/verification/mathematica/02_kinetic_and_susceptibilities.wl new file mode 100644 index 00000000..b71a7cf3 --- /dev/null +++ b/verification/mathematica/02_kinetic_and_susceptibilities.wl @@ -0,0 +1,209 @@ +scriptDirectory = DirectoryName[ExpandFileName[$InputFileName]]; +If[!MemberQ[$Path, scriptDirectory], AppendTo[$Path, scriptDirectory]]; +Needs["checklib`"]; +resetChecks[]; + +repositoryRoot = ExpandFileName[FileNameJoin[{scriptDirectory, "..", ".."}]]; +oraclePath = FileNameJoin[{repositoryRoot, "verification", "oracles", "fp_susceptibilities.dat"}]; +workingPrecision = 70; +$MaxExtraPrecision = 4000; + +(* Thesis (5.14)-(5.15), defined independently of QL-Balance/W2_arr. The + y-rescaling keeps both small- and large-x1 quadratures well conditioned. *) +generatingIntegrand[x1_, x2_, alpha_, beta_, y_] := Module[ + {scale = 1 + x1^2, tau}, + tau = y/scale; + Exp[(I x2 - x1^2) tau + + (alpha + I x1) (beta + I x1) (Exp[-tau] - 1) + + (alpha + beta)^2/2]/scale +]; + +clearMomentCache[] := Clear[momentMatrix]; +momentMatrix[x1in_?NumericQ, x2in_?NumericQ] := + momentMatrix[x1in, x2in] = Module[ + {x1 = SetPrecision[x1in, workingPrecision], + x2 = SetPrecision[x2in, workingPrecision], expressions}, + expressions = Flatten[Table[ + D[generatingIntegrand[x1, x2, alpha, beta, y], {alpha, k}, {beta, l}] /. + {alpha -> 0, beta -> 0}, + {k, 0, 3}, {l, 0, 3} + ]]; + Partition[ + Quiet[ + NIntegrate[ + Evaluate[expressions], {y, 0, Infinity}, + WorkingPrecision -> workingPrecision, + AccuracyGoal -> 48, PrecisionGoal -> 42, MaxRecursion -> 40, + Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0} + ], + NIntegrate::precw + ], + 4 + ] + ]; + +iDifferential[k_Integer, l_Integer, x1_?NumericQ, x2_?NumericQ] := + momentMatrix[x1, x2][[k + 1, l + 1]]; + +(* Thesis (5.13): rank-one energy-conservation correction to the + Ornstein-Uhlenbeck differential moments. Momentum is intentionally not + corrected in the straight-cylinder model. *) +iConserving[k_Integer, l_Integer, x1_?NumericQ, x2_?NumericQ] := Module[ + {denominator}, + denominator = 1 - iDifferential[0, 0, x1, x2] + + 2 iDifferential[2, 0, x1, x2] - iDifferential[2, 2, x1, x2]; + iDifferential[k, l, x1, x2] + + (iDifferential[k, 0, x1, x2] - iDifferential[k, 2, x1, x2]) * + (iDifferential[l, 0, x1, x2] - iDifferential[l, 2, x1, x2]) / + denominator +]; + +mutation = Environment["KAMEL_VERIFY_NEGATIVE_FIXTURE"]; +regenerate = TrueQ[Environment["KAMEL_REGENERATE_ORACLE"] === "1"]; + +thermalSpeed[temperature_, mass_] := Sqrt[temperature/mass] * + If[mutation === "thermal-speed", Sqrt[2], 1]; +x1FromPhysics[case_] := case["kpar"] thermalSpeed[case["temperature"], case["mass"]]/case["nu"]; +x2FromPhysics[case_] := -( + case["omegaE"] + + If[mutation === "flow-sign", -1, 1] case["kpar"] case["Vpar"] + + case["mphi"] case["omegaC"] - case["omega"] +)/case["nu"]; + +caseRecord[id_, kpar_, temperature_, mass_, nu_, omegaE_, vpar_, mphi_, omegaC_, omega_] := <| + "id" -> id, "kpar" -> SetPrecision[kpar, workingPrecision], + "temperature" -> SetPrecision[temperature, workingPrecision], + "mass" -> SetPrecision[mass, workingPrecision], "nu" -> SetPrecision[nu, workingPrecision], + "omegaE" -> SetPrecision[omegaE, workingPrecision], "Vpar" -> SetPrecision[vpar, workingPrecision], + "mphi" -> mphi, "omegaC" -> SetPrecision[omegaC, workingPrecision], + "omega" -> SetPrecision[omega, workingPrecision] +|>; + +cases = { + caseRecord[1, 3/5, 1, 1, 1, 1/5, 2/5, -1, 11/10, 7/10], + caseRecord[2, 3/5, 1, 1, 1, 1/5, 2/5, 0, 11/10, 7/10], + caseRecord[3, 3/5, 1, 1, 1, 1/5, 2/5, 1, 11/10, 7/10], + caseRecord[4, 3/5, 1, 1, 1, 0, 0, 0, 0, -7/4], + caseRecord[5, 3/4, 1, 1, 1, 0, 0, 0, 0, 1/100000000], + caseRecord[6, 2, 1, 1, 10, 0, 0, 0, 0, 1/2], + caseRecord[7, 3/5, 1, 1, 1/10, 0, 0, 0, 0, -4/5] +}; +consumedPairs = {{0, 0}, {2, 0}, {0, 2}, {0, 1}, {2, 1}, {2, 2}, {1, 1}, {1, 3}}; + +checkNumeric["mphi=-1 x2 mapping", x2FromPhysics[cases[[1]]], 34/25, 10^-55]; +checkNumeric["mphi=0 x2 mapping", x2FromPhysics[cases[[2]]], 13/50, 10^-55]; +checkNumeric["mphi=+1 x2 mapping", x2FromPhysics[cases[[3]]], -21/25, 10^-55]; +check["thermal speed is sqrt(T/m)", thermalSpeed[9, 4] == 3/2]; + +(* Collision-model moment contract: the differential OU operator conserves + particle number only; the rank-one integral term restores energy, not + parallel momentum. Krook would damp all three and is not substituted. *) +ouEigenvalues = <|"particle" -> 0, "momentum" -> -1, "energy" -> -2|>; +correctedEigenvalues = Join[ouEigenvalues, <|"energy" -> 0|>]; +check["OU particle moment conserved", ouEigenvalues["particle"] == 0]; +check["OU momentum intentionally damped", ouEigenvalues["momentum"] == -1]; +check["OU energy requires integral correction", ouEigenvalues["energy"] == -2]; +check["integral term restores energy", correctedEigenvalues["energy"] == 0]; +check["FP model is not non-conserving Krook", Values[ouEigenvalues] =!= {-1, -1, -1}]; +If[MemberQ[{"thermal-speed", "flow-sign"}, mutation], FinishChecks[]]; + +x1Probe = SetPrecision[3/5, workingPrecision]; +x2Probe = SetPrecision[-7/4, workingPrecision]; +recurrenceRight = If[ + mutation === "recurrence-index", + (x2Probe/x1Probe) iConserving[2, 1, x1Probe, x2Probe], + (x2Probe/x1Probe) iConserving[1, 0, x1Probe, x2Probe] +]; +checkNumeric["I11 recurrence (A.4)", iConserving[1, 1, x1Probe, x2Probe], recurrenceRight, 10^-38]; +checkNumeric[ + "I31 recurrence (A.4)", + iConserving[3, 1, x1Probe, x2Probe], + (x2Probe/x1Probe) iConserving[2, 1, x1Probe, x2Probe], + 10^-38 +]; +checkNumeric[ + "susceptibility symmetry", + iConserving[1, 3, x1Probe, x2Probe], + iConserving[3, 1, x1Probe, x2Probe], + 10^-42 +]; +checkNumeric[ + "x1 parity", + iConserving[1, 3, -x1Probe, x2Probe], + iConserving[1, 3, x1Probe, x2Probe], + 10^-38 +]; +branchExpected = (-1)^(1 + 3) Conjugate[iConserving[1, 3, x1Probe, x2Probe]]; +If[mutation === "wrong-branch", branchExpected = -branchExpected]; +checkNumeric[ + "Landau-branch conjugation", + iConserving[1, 3, x1Probe, -x2Probe], + branchExpected, + 10^-38 +]; +If[MemberQ[{"wrong-branch", "recurrence-index"}, mutation], FinishChecks[]]; + +formatNumber[value_] := StringReplace[ + ToString[FortranForm[N[value, 34]], InputForm], + {"*^" -> "E", " " -> ""} +]; + +oracleRows = Flatten[Table[ + Module[{x1 = x1FromPhysics[case], x2 = x2FromPhysics[case], value}, + Table[ + value = iConserving[pair[[1]], pair[[2]], x1, x2]; + { + ToString[case["id"]], formatNumber[case["kpar"]], + formatNumber[case["temperature"]], formatNumber[case["mass"]], + formatNumber[case["nu"]], formatNumber[case["omegaE"]], + formatNumber[case["Vpar"]], ToString[case["mphi"]], + formatNumber[case["omegaC"]], formatNumber[case["omega"]], + ToString[pair[[1]]], ToString[pair[[2]]], formatNumber[x1], formatNumber[x2], + formatNumber[Re[value]], formatNumber[Im[value]] + }, + {pair, consumedPairs} + ] + ], + {case, cases} +], 1]; + +oracleHeader = { + "# KAMEL FP susceptibility oracle; generated from Markl thesis (5.13)-(5.15), not W2_arr.", + "# Wolfram NIntegrate working precision: 70 digits; committed values: 34 digits.", + "# Cases 1-3: finite parallel flow and mphi=-1,0,+1; 4: general complex point;", + "# 5: near-zero resonance numerator; 6: strong collision; 7: collisionless scaling.", + "# case kpar T_erg mass nu omega_E Vpar mphi omega_c omega k l x1 x2 Re(Ikl) Im(Ikl)" +}; + +If[regenerate, + Export[ + oraclePath, + StringRiffle[Join[oracleHeader, StringRiffle[#, " "] & /@ oracleRows], "\n"] <> "\n", + "Text" + ]; + Print["Wrote ", Length[oracleRows], " rows to ", oraclePath], + Module[{lines, parseNumber, committedRows}, + lines = Select[ + StringSplit[Import[oraclePath, "Text"], "\n"], + StringLength[StringTrim[#]] > 0 && !StringStartsQ[StringTrim[#], "#"] & + ]; + parseNumber[token_] := ToExpression[StringReplace[ + token, RegularExpression["[Ee]([+-]?[0-9]+)$"] -> "*^$1" + ]]; + committedRows = (parseNumber /@ StringSplit[#]) & /@ lines; + check["56 committed oracle rows", Length[committedRows] == Length[oracleRows] == 56]; + If[Length[committedRows] == Length[oracleRows], + Do[ + checkNumeric[ + "oracle row " <> ToString[row], + committedRows[[row, -2]] + I committedRows[[row, -1]], + parseNumber[oracleRows[[row, -2]]] + I parseNumber[oracleRows[[row, -1]]], + 10^-30 + ], + {row, 1, Length[oracleRows]} + ] + ] + ] +]; + +FinishChecks[]; diff --git a/verification/mathematica/checklib.wl b/verification/mathematica/checklib.wl index f2be8aa3..f4d80459 100644 --- a/verification/mathematica/checklib.wl +++ b/verification/mathematica/checklib.wl @@ -1,6 +1,7 @@ BeginPackage["checklib`"]; check::usage = "check[name, condition] records one Boolean verification."; +checkNumeric::usage = "checkNumeric[name, actual, expected, tolerance] checks a scaled numerical error."; registerConvention::usage = "registerConvention[name, definition, source] rejects incompatible duplicates."; registeredConventions::usage = "registeredConventions[] returns the convention registry."; resetChecks::usage = "resetChecks[] clears checks and registered conventions."; @@ -23,6 +24,12 @@ check[name_String, condition_] := Module[{passed = TrueQ[condition]}, passed ]; +checkNumeric[name_String, actual_, expected_, tolerance_] := Module[ + {scale = Max[1, Abs[N[expected]]], error}, + error = Abs[N[actual] - N[expected]]/scale; + check[name <> " (scaled error " <> ToString[ScientificForm[error, 3]] <> ")", error <= tolerance] +]; + registerConvention[name_String, definition_String, source_String] := Module[{}, If[StringLength[StringTrim[name]] == 0 || StringLength[StringTrim[definition]] == 0 || diff --git a/verification/oracles/fp_susceptibilities.dat b/verification/oracles/fp_susceptibilities.dat new file mode 100644 index 00000000..05383d14 --- /dev/null +++ b/verification/oracles/fp_susceptibilities.dat @@ -0,0 +1,61 @@ +# KAMEL FP susceptibility oracle; generated from Markl thesis (5.13)-(5.15), not W2_arr. +# Wolfram NIntegrate working precision: 70 digits; committed values: 34 digits. +# Cases 1-3: finite parallel flow and mphi=-1,0,+1; 4: general complex point; +# 5: near-zero resonance numerator; 6: strong collision; 7: collisionless scaling. +# case kpar T_erg mass nu omega_E Vpar mphi omega_c omega k l x1 x2 Re(Ikl) Im(Ikl) +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 0 0 0.6 1.36 0.1398646758452037469735238584320868 0.803447617927279868240472235082642 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 2 0 0.6 1.36 0.4611248446171895266532985801210965 0.8785352035215942782821385929902574 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 0 2 0.6 1.36 0.4611248446171895266532985801210965 0.8785352035215942782821385929902574 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 0 1 0.6 1.36 0.3170265985824618264733207457793968 0.1544812673018343680117370661873218 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 2 1 0.6 1.36 1.045216314465629593747476781607819 0.3246797946489470307728474774445835 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 2 2 0.6 1.36 1.82802398837384869453953490923675 2.477968058646995925997582251554087 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 1 1 0.6 1.36 0.7185936234535801400061936904332995 0.350157539217491234159937350024596 +1 0.6 1. 1. 1. 0.2 0.4 -1 1.1 0.7 1 3 0.6 1.36 2.369156979455427079160947371644389 0.7359408678709466030851209488743892 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 0 0 0.6 0.26 1.208511582165335305212297806950257 2.766757230090777017377503015285048 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 2 0 0.6 0.26 1.006495842029929561428334854932569 0.6701272225475658492497684267787116 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 0 2 0.6 0.26 1.006495842029929561428334854932569 0.6701272225475658492497684267787116 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 0 1 0.6 0.26 0.5236883522716452989253290496784448 -0.4677385336273299591364153600431457 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 2 1 0.6 0.26 0.4361481982129694766189451038041132 -1.376278203562721465325100348395892 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 2 2 0.6 0.26 2.482794558496822548743376792308268 0.1305264421444364927240316887019689 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 1 1 0.6 0.26 0.2269316193177129628676425881939928 -0.2026866979051763156257799893520298 +2 0.6 1. 1. 1. 0.2 0.4 0 1.1 0.7 1 3 0.6 0.26 0.1889975525589534398682095449817824 -0.5963872215438459683075434843048864 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 0 0 0.6 -0.84 0.3535783365377262114073397280376117 -1.084620726581397714204721024215933 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 2 0 0.6 -0.84 0.9400096220351264856584812548943187 -0.6175394093542340131250459062176569 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 0 2 0.6 -0.84 0.9400096220351264856584812548943187 -0.6175394093542340131250459062176569 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 0 1 0.6 -0.84 -0.4950096711528166959702756192526564 -0.1481976494527098667800572327643598 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 2 1 0.6 -0.84 -1.316013470849177079921873756852046 -0.802111493570739048291602397961947 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 2 2 0.6 -0.84 3.17927134847341299237662725619611 -1.070399693749593798928212904273351 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 1 1 0.6 -0.84 0.6930135396139433743583858669537189 0.2074767092337938134920801258701038 +3 0.6 1. 1. 1. 0.2 0.4 1 1.1 0.7 1 3 0.6 -0.84 1.842418859188847911890623259592865 1.122956090999034667608243357146726 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 0 0 0.6 -1.75 0.05286527559367708087321982786766412 -0.6263017022554388431699162747190108 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 2 0 0.6 -1.75 0.182977576551716901741339561351257 -0.7237866540172676353777869169316186 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 0 2 0.6 -1.75 0.182977576551716901741339561351257 -0.7237866540172676353777869169316186 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 0 1 0.6 -1.75 -0.154190387148224819213557831280687 0.1600466315783632925789224679304481 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 2 1 0.6 -1.75 -0.5336845982758409634122403872744995 0.444377740883697269851878507717221 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 2 2 0.6 -1.75 0.8159505101650406935325702833552553 -2.185576074703851976088379626299394 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 1 1 0.6 -1.75 0.4497219625156557227062103412353371 -0.4668026754368929366885238647971403 +4 0.6 1. 1. 1. 0 0 0 0 -1.75 1 3 0.6 -1.75 1.556580078304536143285701129550624 -1.296101744244117037067978980841895 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 0 0 0.75 1.e-8 5.924644544044985648267357299604898 4.085660433046205236750243809563205e-7 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 2 0 0.75 1.e-8 1.777777777777771567651593525854805 8.75492363385775952698273839255135e-8 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 0 2 0.75 1.e-8 1.777777777777771567651593525854805 8.75492363385775952698273839255135e-8 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 0 1 0.75 1.e-8 7.899526058726647531023143066139864e-8 -1.333333333333327885786089271726351 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 2 1 0.75 1.e-8 2.370370370370362090202124701139741e-8 -1.33333333333333216601018215229873 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 2 2 0.75 1.e-8 1.777777777777776537396292252446585 1.382716049382706565589256731788015e-8 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 1 1 0.75 1.e-8 1.053270141163553004136419075485315e-15 -1.777777777777770514381452362301801e-8 +5 0.75 1. 1. 1. 0 0 0 0 1.e-8 1 3 0.75 1.e-8 3.160493827160482786936166268186321e-16 -1.777777777777776221346909536398306e-8 +6 2. 1. 1. 10. 0 0 0 0 0.5 0 0 0.2 0.05 4.599265301175717153746711301920209 15.53410359971119189221208918047826 +6 2. 1. 1. 10. 0 0 0 0 0.5 2 0 0.2 0.05 5.869824581684492456844057980772183 5.469963101451595935446644701180153 +6 2. 1. 1. 10. 0 0 0 0 0.5 0 2 0.2 0.05 5.869824581684492456844057980772183 5.469963101451595935446644701180153 +6 2. 1. 1. 10. 0 0 0 0 0.5 0 1 0.2 0.05 1.149816325293929288436677825480052 -1.116474100072202026946977704880434 +6 2. 1. 1. 10. 0 0 0 0 0.5 2 1 0.2 0.05 1.467456145421123114211014495193046 -3.632509224637101016138338824704962 +6 2. 1. 1. 10. 0 0 0 0 0.5 2 2 0.2 0.05 18.52941015954078585924444774732307 6.429153420946340317020487769788988 +6 2. 1. 1. 10. 0 0 0 0 0.5 1 1 0.2 0.05 0.2874540813234823221091694563700131 -0.2791185250180505067367444262201085 +6 2. 1. 1. 10. 0 0 0 0 0.5 1 3 0.2 0.05 0.3668640363552807785527536237982614 -0.9081273061592752540345847061762405 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 0 0 6. -8. 0.07272348505035489007650962391837135 -0.1180077660340046305094739348686044 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 2 0 6. -8. 0.1308400254152965533561340125507481 -0.003730136293864874256066911748268209 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 0 2 6. -8. 0.1308400254152965533561340125507481 -0.003730136293864874256066911748268209 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 0 1 6. -8. -0.09696464673380652010201283189116179 -0.009322978621327159320701420175194067 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 2 1 6. -8. -0.1744533672203954044748453500676641 -0.1616931516081801676585774510023091 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 2 2 6. -8. 0.2595533482285572339095567085906036 0.1865153076075076561322957096584681 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 1 1 6. -8. 0.1292861956450753601360171091882157 0.01243063816176954576093522690025876 +7 0.6 1. 1. 0.1 0 0 0 0 -0.8 1 3 6. -8. 0.2326044896271938726331271334235521 0.2155908688109068902114366013364121