From 3ac5be2e2d9445682587a3cf2ea8de18c8f127aa Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Mon, 2 Nov 2026 12:23:41 +0100 Subject: [PATCH] Repo - BUGFIX! - Make reruns self-contained, flag vanished harness outputs, and remove the placeholder RLAR and CLAR torque - Rerun: record the kinetic profiles (Input/RawInputs/Kinetic) and a file wall's text (Input/RawInputs/Wall/wall_file_raw) in gpec.h5; a rerun writes them back out as _kinetic.h5 and _wall.dat, and equilibrium_from_h5 reads the profiles from the snapshot, so a replay never re-reads the files the source run named. A stored deck naming either file without its snapshot stops with a re-run-from-the-deck error. - Rerun: refuse an output or side file that is the source under another name (symlink, hardlink, case-only alias), not only the same path. - Regression harness: a quantity that disappears is CHANGED (missing) and a storage-type change is CHANGED (type), and both fail --fail-on-change; a quantity only the newer ref has is NEW. - KineticForces: remove the RLAR method (a linear-amplitude placeholder for the kappa integral) and CLAR (which returned zero); rlar_flag and clar_flag are no longer keys. - Vacuum: build the 3D singular quadrature under a lock and return the data built for the caller's parameters. - Repo: drop the committed __pycache__ bytecode and ignore __pycache__/. Co-Authored-By: Claude Opus 5.5 --- .gitignore | 3 + .../__pycache__/check_subject.cpython-310.pyc | Bin 11359 -> 0 bytes docs/development/hdf5-conventions.md | 4 +- docs/development/regression-harness.md | 7 +- docs/src/workflow.md | 2 +- regression-harness/src/extractor.jl | 15 ++- regression-harness/src/reporter.jl | 26 ++++-- regression-harness/test/runtests_compare.jl | 18 ++++ src/Equilibrium/KineticProfiles.jl | 88 ++++++++++-------- src/ForceFreeStates/InputSnapshot.jl | 13 +++ src/ForceFreeStates/Output.jl | 1 + src/KineticForces/KineticForcesStructs.jl | 9 +- src/KineticForces/Torque.jl | 74 +-------------- src/KineticForces/Utils.jl | 2 +- src/Rerun.jl | 74 ++++++++++++--- src/Vacuum/Kernel3D.jl | 23 +++-- src/Vacuum/PnQuadCache.jl | 2 +- test/runtests_error_fields.jl | 3 + test/runtests_frozen_equilibrium.jl | 15 ++- test/runtests_kinetic.jl | 2 +- test/runtests_rerun_from_h5.jl | 35 +++++++ 21 files changed, 245 insertions(+), 171 deletions(-) delete mode 100644 ci/conventions/__pycache__/check_subject.cpython-310.pyc create mode 100644 regression-harness/test/runtests_compare.jl diff --git a/.gitignore b/.gitignore index 30f8dc867e..5e4d6cc2fe 100644 --- a/.gitignore +++ b/.gitignore @@ -44,3 +44,6 @@ profiling/ # Per-developer package preferences (e.g. the precompile_workload switch) LocalPreferences.toml + +# Python bytecode from the CI convention scripts +__pycache__/ diff --git a/ci/conventions/__pycache__/check_subject.cpython-310.pyc b/ci/conventions/__pycache__/check_subject.cpython-310.pyc deleted file mode 100644 index 6743b0f19cdedc4a53453badf6f885a4271ea18e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 11359 zcmai4TW}=TS?=!Xx%G^+(t5oYt=I0o*xr#flD)Q*uwK@gl}1|GT1nAmZEJTuYRzel zMw;86?%7>+TCipV6c9s#lj7nmNJ)iY;2}`F0C|E?MezbWP(^hG9-v5-9|({DfenfI z{=a8NBY6{MYR;VQ{?F|{|Ktw6-6jZ zicp1CQLEaLrcympi7Z9BM+>)Zc@-le{>A`Xg2-ciLN@hHl};t0y4 z;uy*SaUA7iT%O?aaq$H1nBpYLQ{ptrGvY~{v%YxVlZz8TQ~IcZxaw=&DLCfgMCqGMPkt6T%GmK|9f zSspe=_p-!Xtye1b+f46RRog5{t6H_Bxzcn^OWG!QYTC8cx-8nl%w~<{Y-=9bs#99# zY2%jT3QSenjkbDCx+EEtJo;C!Frczu_MC1WmoDAx;+x|akwg$(H+YVvk2RV zRS8>`Z=oemBV-Lpf+#y1rnE3wi5XipOjn9x9il&EI(5?$SY-**H>6#y-=ex@R_mgP zDK#4;pAIwf7FM;=tXOi>w3=?cih1)ghmH5iNM+O@W!1n~UXL@4Ykd}jW@15nJfB-w zoXZEXON$fZldt1Jer`O7)L4FYEwR1cYU#j#CV*^t?uI0*%>#XHrxg?5fBn?>|)Ra&i>Mn-8> zIu$W0b{E+jq$aP-&dprSU&&7|1o6?Sd~SMiHi!^2Bh&fUg7nPn!sM07rI3#)n#IiA zmE81bK8O(miayE#K&P#^z`I3?#$({PNLJ; zm2e-T`A;lPP6eqhP6UywxzQkY)hafdRm@|2J~%K=>NGBGdmb_Xj%%~|;LyX^80ZtT z`B8ZYi-ZtNO$CW7i&G1evr~|e+yquQFZ<}3cKJnlkSftj`RUQivWM#N{Oh?ZAeH@e zm6*s)<>y|Voc50H8-92KIwd#CwIDJv@ru_wS;N*&S(~;T#&r;zpUSp3b7H6l3GttqR$*K7uO=Ce(7HnB{Al@^JbCdbGAi3y5Qrxob zcu&81_KGb_w&~A7d!=jotK~b*#!26~*|eOslRm6cO`Obl`8UtLJbMA!YdIIM&p(qf zvnJhet?qp*t&9t#t9^%WVRuNd|2bZH5OkXS~C2o<4oud1fy38D!t= z>(1aA@7Xu;cP%@31Ao`?{ba>ju`65Wsg7IJ@hz*;v@Z@1zMSC}gN5s7+gAMJJQrup zx_Ir`?2lR5)m(P`#wR~I^un`Snd_N?cjCgeH!t4!1`3e&Ij;Cf&WuR3Cs+l;wmn{$F{bktOLLtQ??Xe5!!KOTi;P+zpu=L zVmrEn|4hWoO~DeuCz3ZH$IFC$>XBQ+2h1GJ!E);?|3`4zg%MRiV>`Kjbf>_xph_VbK#oGc_Em?%=X21QO zBJ=3z!PrO@rWr$7ZgFj-R;-+R?%c;lZo4C;hFvVcd}lxM?75L*d8Eq=kFeV3w7S_-rZiB6paii28O)zyS5CuqJpJ$sO>zzto41-VfqHnu-xG2-|!OTUT=7^ogs zwp3r;R(2G}D1mt;WlM8`8Aaq%3hsm5W6C`+eoVQM*^2nucEpdY>x5K5r88Az57nbv zx*PTN?U;z}pq+fV-R50I#C&~Se^=q?P=-G<;(RJO8Z$#uGw~MYP znD(--VPoUgO)O}PeOTN7kdp+9THYu}TJkXO24866lOtN}6i^Wf1T?KuUWGQAXl3o| z2ab_EG@%xF=ZbL5SDKYFsOn8vF=vQ73<1)HZ8sbfj?jWD4Lhe!ohHB&dS31}2H4_q zrBJf#Roj)IS~lx|9Q9T1js$ksiw>LPkhts$FkaUp+nyXA2cF2aF-Za|yz~&xg<{Da zw)j$Z0A}jV(i)=^jBCO05bJ2bmZ30<7PO+}%7yJ(<;pPkAumE$_J|TGM+83y+WDQw zPhyc6yacNEv_KnvHc%U^6k#!r`$6(T32yaHLtgw@ zTz-cfoOX9W@tz2$Utr@me1TDub8)ys{X-`uUqKi0RVt>axIzVK;_jew7MDLqw?2Ww z(+AnRl5{&}H5xF;fewIW2NA2b8R(@-eI@9b&5gbar&NGj zntLPA9oMeHXk4bQiE}}erys<^#lWaEVfaER+Bc^l-MwHBjU5w3B9VZS>)L>tRtMBx z^{A@d)qZd=X5d=8tE;jAqCLDIy{8RgT_;cJVuS9fuIg&{)Fapw=1nHTEd!klM$k7s zGhMhcGqyODFXZO(x%rG5E@Cg=f)w0AS%iD?d35cNH1w%*1%E-*t#8y(3{UUb4B`)o^gQ{p z!>E)=TjFhP@zr;xWT$5WN$SQqgfbvq^C&bFd%+DBg^xuMG?9>rj) zC|uy7t(YJ4RgB-@Cn-P1pl_%nY1_75m(Y-im`8%>yGHa7Zlp>~NsP z;jPpU92A!Ow#9^i1JWZh9iT ziaHHM2|=U+Q<`p(g(f*J?@%wK0rDrPml%{2T$Y=35e3r$M|Jw!uVg)Y;N)+C+&POP z0X5TAL(|n5Jd6RwI|YS1toB9HssRN|sE?74(d5sA_(Num0#H`E1yqnf)PNcpF>WE( zC}5ghbw?$2;5yWcP$1VMu6{3k7N=`JtS9&xkR#}l_fuq)e}pYPia|?cOA};ET?JYL z_ft^G_fE|A)=A^kcC*^nMpyeC)d61!RGf?~Viqtn7W(aNu2cb)Tq zO5ORm^QNZ_nO;9`+jokUCam9iLZGJiJZ>$YT!yM49|1qpog!$7;!Arp0tp&#U{5ZG z1`Ck_5_n!3S9QwW&;Z^^!|C;lGl@eeROJZ zEI$^CEL@f%?@>W~%xEllLA*8ges1=Eh;4sHOdUYct7;G}{ey(2{~#UFzn_Te-%G?| z53~>DJJdblzCfRbwb}y-$ZbaPV;LhvSn3`5xuk6|B6i8$E@lMB7+Uad9YDaYvyH2{SY{5xmw+G=w;l|-c zHFB6D0$_L=4ih%VkCMG<#hAT*Ob)s_dWykKzrYkRrikxS_-J!Mx!DUS8D$3@_v3QP zkBcM>$`rRT{P?@teJ~W#MpP-d@r0j{1u(nd#?U_;yHjA=G*hYaEjPZMKzy5IXqxns z@-08PZHOMsrg=dcsQrLw!$Tp5{NgqxZ}8V{+!!U7!I z6f%CuB4?0nP2`%Qj1;QCE4@(>$U`Mc=O`}E-8{*`okpfGro1LM}g1I5Vge3t~| zu#iNuhfBlea*fbcp(`0?J->oeqI1t>3K(w1E?N{m&^ZfjO&ijXg>lSg4c5e@(+Np| zp8P#3>L?IlqT42+1=gHSTn?wf>G;9H4$o=QErdi1r-sXnaU5$Chc^~r75kTvH{k>_ zJ;9+phW~OtH-^-u{36YORF;uSkb($|<}c4ojpgTLgC^E~2G&gcjL6izBOwi|Tm9mjARauTv zImzXmx;X&&69-7rx8sx`+fGo)k!;jWa#WDmP~=B$4#0vNh*y)xC{}@q2?={dw{CJ1PV44HI z0h+!RP1#St^z?iAoCBr9fg?14(QoIRt#ir*m39Eh?K?0A$V3rk3mvo%mWI6pKqIaA zf{)0&V|?^<_S%~x%Qv2R8OiJOBmQ;dUoy)NqD^--`!P7;1Po4{&Yd2p{)cd9C!D4e z0)~LGb~xPy|KRjkGRq*^Cp_{Oaqj@!3`4#{NHvAR#0%_-pF3dhee%B?TgVD(KYT_e zo~{unP49wIflUg*=`S-l#l+^QyI-K<(^TX;#W-q6h&gr$OAoi_Xe!?051(V3)NA1~ z(}R7YqezEN2+6|msb(wAY8{g9Zl~A3TW{^AcdC8R+78Y-c}+UAtURanuSQh#r{q|QhwFrP>5tYt|&Y}|d zjZ)g@i-2+fXaH;|k;r;mqvQlk1v2fFqltP9JY63;aUb0sm+92?BljZQ zf8h3r$K9jr$HWsHD(*=rKBLN*JHX=%+*8|g;=4WK z6vi8XJ2-tmvfI-CQKNjf)$eii@QN7Fl`Uf{g|?4jW@oq3etP@Fj`A7R!MHd3eM6l2 zC@8>>bl1oq^HY33g}q^V%-~7!)b`_qLn%Q;dqeH9qNjbsPaP$fbrwcm5g9*)a`1kf zXnY;Bm=R~NpV-sPGt)p?DPf0)Vs!X@APfqSWI@8uavU$+{USu5k55kJy>oEpyAZWI z-9#BFBv`_iFRD0Dr-KIqKo}>E6Ap}Nn1d^o)1PBGmodS;8ZSN_;cIpFR`!dBMT_x{*(f(K9D!5hZW(ls%Rj-C<=0TaK4X4_G&%A6 zP*3KN&3XzXZ;>bKZTnVAuWYM=qB7K4M?3&Q8SR038ln?gcCzuPQ|?;+B^7&}8y~#9 zNM@G2xOVr9;vK2tWt}ds+dd6i@XU6M{^8JqUwAq%KTjk43JTBs{}LZwZx?T|8QuVe zAmQZL(76RTp$it!j^%-UD36}N4()F zOR6Pc1=`6!AOf{sb%wxB8brgLopwvxFk^JSk5Ty6e~~4iRk;y5wkx062fi zd6RcRM+ZiD11gBDSmI&A66hovPmr{@An$|A?fxdz7^Rg^QoaHIA5i@;fqq=k(Rhpx z$jNcXw8v4>p~7Kys(lZL^no7N{zp$}|E(vr z@9Kv3UwTUWPd%;u2YURD-m86EKcM})-lzSW-mm?uelYeeN;{O97TtQX< diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 99d1887460..403fb7ffc5 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -29,7 +29,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness ## Inputs live only under `Input/` -`Input/gpec_toml_raw` stores the deck that reproduces the run: each result's save adds the sections stating the problem, algorithm and controls it was computed with (`Utilities.record_inputs!`), so a deck run and a script run of the same workflow record the same deck. `Input/RawInputs/` stores the raw equilibrium/forcing/coil data and tolerance file — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer and rerun reader carry cross-reference comments marking the mirrored path pair). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) +`Input/gpec_toml_raw` stores the deck that reproduces the run: each result's save adds the sections stating the problem, algorithm and controls it was computed with (`Utilities.record_inputs!`), so a deck run and a script run of the same workflow record the same deck. `Input/RawInputs/` stores the raw equilibrium/forcing/coil data, the kinetic profiles, the wall file and the tolerance file — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer and rerun reader carry cross-reference comments marking the mirrored path pair). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) ## Schema @@ -38,7 +38,7 @@ Top level (11 groups): | Group | Contents | |---|---| | `Info/` | Run metadata: `git_version`, mode-number ranges (`mpert`, `mlow`, …, `mn_index`), `psilim`, `qlim`, `Runtimes/` (per-stage wall-clock seconds) | -| `Input/` | Rerun snapshot: `gpec_toml_raw`, `RawInputs/{Equilibrium, ForcingTerms, Coils/}` | +| `Input/` | Rerun snapshot: `gpec_toml_raw`, `RawInputs/{Equilibrium, ForcingTerms, Coils/, Kinetic, Wall, ErrorFields}` | | `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) | | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, `gal_basis_output`-gated `Basis/`), `write_el_matrices`-gated `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` | diff --git a/docs/development/regression-harness.md b/docs/development/regression-harness.md index 0dab467f41..f89d7d5ee8 100644 --- a/docs/development/regression-harness.md +++ b/docs/development/regression-harness.md @@ -72,6 +72,11 @@ singular q values [4 elements] [4 elements] 0 Summary: 23 unchanged, 3 missing/N/A ``` +Each row's status is `OK` (within the case's noise threshold), `** CHANGED **` (moved past it), +`** CHANGED (missing) **` (on the earlier ref, gone from the later one), `** CHANGED (type) **` (stored +with a different shape or type), `NEW` (only on the later ref), `N/A` (on neither ref) or `FAILED`. +Every `CHANGED` status counts toward `--fail-on-change`; `NEW` and `N/A` do not. + **Compare your uncommitted working tree against develop:** ```bash regress --cases solovev_n1 --refs develop,local @@ -102,7 +107,7 @@ regress --cases solovev_n1 --ref-range develop~10..develop - `--no-instantiate` — skip `Pkg.instantiate()` (faster if deps are already resolved) - `--no-pin-manifest` — let each ref resolve its own package set (see below) - `--allow-env-mismatch` — reuse cached results produced in a different environment -- `--fail-on-change` — exit non-zero when any tracked quantity changed +- `--fail-on-change` — exit non-zero when any tracked quantity changed, disappeared or changed storage type GPEC subprocesses run with `-t auto` (all cores) so GPEC's threaded kernels are active; set `GPEC_REGRESS_THREADS=1` to force single-threaded runs. Tracked quantities are thread-count independent, and the count each run actually used is recorded in its environment fingerprint (shown in the report's `env:` lines). Thread count is deliberately not part of the cache key, so `Runtime (s)` rows cached from single-threaded runs are not comparable to threaded ones — re-baseline with `--force` if runtime tracking matters. diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 02e227150e..ef6fab0284 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -184,7 +184,7 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The top-leve | Group | Contents | |---|---| | `Info/` | Run metadata: git version, mode-number ranges, ψ limit, `Runtimes/` (per-stage wall-clock seconds) | -| `Input/` | Self-contained rerun snapshot: the recorded deck, raw equilibrium/forcing/coil inputs | +| `Input/` | Self-contained rerun snapshot: the recorded deck, raw equilibrium/forcing/coil inputs, kinetic profiles and wall file | | `Equilibrium/` | Equilibrium scalars (`beta_N`, `q_axis`, `q_95`, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) | | `ForceFreeStates/` | Stability solve: `Solutions/{ForwardIntegration,GalerkinIntegration}`, `EulerLagrangeMatrices/` (opt-in: `write_el_matrices = true`), `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier D_I, resistive interchange D_R, ballooning Δ' profiles | diff --git a/regression-harness/src/extractor.jl b/regression-harness/src/extractor.jl index 33a328d7cb..222312553a 100644 --- a/regression-harness/src/extractor.jl +++ b/regression-harness/src/extractor.jl @@ -155,21 +155,26 @@ function _json_element_abs(x)::Float64 end """ -Compare two extracted values. Returns (abs_diff, rel_diff, status). -For arrays, returns max element-wise absolute diff. +Compare two extracted values, `q1` from the baseline ref and `q2` from the newer one. Returns +(abs_diff, rel_diff, status). For arrays, returns max element-wise absolute diff. A quantity +missing from both refs is `N/A`; one that disappears or changes storage type is a change, and +one that only the newer ref has is `NEW`. """ function compare_values(q1::NamedTuple, q2::NamedTuple) vtype = q1.value_type - # Handle missing values - if vtype == "missing" || q2.value_type == "missing" + if vtype == "missing" && q2.value_type == "missing" return (NaN, NaN, "N/A") + elseif q2.value_type == "missing" + return (NaN, NaN, "CHANGED (missing)") + elseif vtype == "missing" + return (NaN, NaN, "NEW") end # Storage type differs between refs (e.g. a quantity whose schema changed from a # vector to a scalar across the two commits); the values are not directly comparable. if vtype != q2.value_type - return (NaN, NaN, "N/A (type changed)") + return (NaN, NaN, "CHANGED (type)") end threshold = q1.noise_threshold diff --git a/regression-harness/src/reporter.jl b/regression-harness/src/reporter.jl index c4b6b94ac0..d2e6990edc 100644 --- a/regression-harness/src/reporter.jl +++ b/regression-harness/src/reporter.jl @@ -114,11 +114,11 @@ function report_two_ref_comparison(db::SQLite.DB, case_spec::CaseSpec, if info1 === nothing println("ERROR: No results for ref 1 ($(ref1.name))") - return (n_ok=0, n_changed=0, n_missing=0, n_failed=1) + return (n_ok=0, n_changed=0, n_missing=0, n_new=0, n_failed=1) end if info2 === nothing println("ERROR: No results for ref 2 ($(ref2.name))") - return (n_ok=0, n_changed=0, n_missing=0, n_failed=1) + return (n_ok=0, n_changed=0, n_missing=0, n_new=0, n_failed=1) end failed1 = !info1.success @@ -135,6 +135,7 @@ function report_two_ref_comparison(db::SQLite.DB, case_spec::CaseSpec, n_ok = 0 n_changed = 0 n_missing = 0 + n_new = 0 # Helper: pick the value-cell text for one ref/quantity, accounting for # whole-ref failure (which should render as "FAILED" rather than "N/A"). @@ -179,14 +180,16 @@ function report_two_ref_comparison(db::SQLite.DB, case_spec::CaseSpec, q1 = q1_all[qname] q2 = q2_all[qname] abs_diff, rel_diff, status = compare_values(q1, q2) - st = status == "CHANGED" ? "** CHANGED **" : status + st = startswith(status, "CHANGED") ? "** $status **" : status push!(rows, [spec.label, format_value(q1), format_value(q2), format_diff(abs_diff, rel_diff, status, q1.value_type), st]) if status == "OK" n_ok += 1 - elseif contains(status, "CHANGED") + elseif startswith(status, "CHANGED") n_changed += 1 + elseif status == "NEW" + n_new += 1 else n_missing += 1 end @@ -231,10 +234,11 @@ function report_two_ref_comparison(db::SQLite.DB, case_spec::CaseSpec, parts = String[] n_changed > 0 && push!(parts, "$n_changed changed") n_ok > 0 && push!(parts, "$n_ok unchanged") + n_new > 0 && push!(parts, "$n_new new") n_missing > 0 && push!(parts, "$n_missing missing/N/A") println("Summary: ", join(parts, ", ")) println() - return (n_ok=n_ok, n_changed=n_changed, n_missing=n_missing, + return (n_ok=n_ok, n_changed=n_changed, n_missing=n_missing, n_new=n_new, n_failed=count(identity, (failed1, failed2))) end @@ -269,6 +273,7 @@ function report_multi_ref(db::SQLite.DB, case_spec::CaseSpec, n_ok = 0 n_changed = 0 n_missing = 0 + n_new = 0 for spec in case_spec.quantities qname = spec.name @@ -316,13 +321,15 @@ function report_multi_ref(db::SQLite.DB, case_spec::CaseSpec, q_prev = haskey(all_qs[end-1], qname) ? all_qs[end-1][qname] : nothing if q_last !== nothing && q_prev !== nothing abs_diff, rel_diff, status = compare_values(q_prev, q_last) - st = status == "CHANGED" ? "** CHANGED **" : status + st = startswith(status, "CHANGED") ? "** $status **" : status push!(row, format_diff(abs_diff, rel_diff, status, q_prev.value_type)) push!(row, st) if status == "OK" n_ok += 1 - elseif contains(status, "CHANGED") + elseif startswith(status, "CHANGED") n_changed += 1 + elseif status == "NEW" + n_new += 1 else n_missing += 1 end @@ -384,11 +391,12 @@ function report_multi_ref(db::SQLite.DB, case_spec::CaseSpec, parts = String[] n_changed > 0 && push!(parts, "$n_changed changed") n_ok > 0 && push!(parts, "$n_ok unchanged") + n_new > 0 && push!(parts, "$n_new new") n_missing > 0 && push!(parts, "$n_missing missing/N/A") println("Summary (last vs prev): ", join(parts, ", ")) end println() - return (n_ok=n_ok, n_changed=n_changed, n_missing=n_missing, + return (n_ok=n_ok, n_changed=n_changed, n_missing=n_missing, n_new=n_new, n_failed=count(identity, failed_mask)) end @@ -434,7 +442,7 @@ function show_quantity_history(db::SQLite.DB, case_name::String, qty_name::Strin else abs_diff, rel_diff, status = compare_values(prev_q, q) diff_str = format_diff(abs_diff, rel_diff, status, q.value_type) - st = status == "CHANGED" ? "** CHANGED **" : status + st = startswith(status, "CHANGED") ? "** $status **" : status push!(rows, [short, date_str, val_str, diff_str, st]) end prev_q = q diff --git a/regression-harness/test/runtests_compare.jl b/regression-harness/test/runtests_compare.jl new file mode 100644 index 0000000000..3a3ff07742 --- /dev/null +++ b/regression-harness/test/runtests_compare.jl @@ -0,0 +1,18 @@ +using HDF5, JSON, Test + +include(joinpath(@__DIR__, "..", "src", "types.jl")) +include(joinpath(@__DIR__, "..", "src", "extractor.jl")) + +# A stored quantity as get_quantities returns it. +stored(value_type; value_real=nothing) = + (label="q", value_real=value_real, value_int=nothing, value_text=nothing, value_type=value_type, noise_threshold=1e-10) + +@testset "A quantity that disappears or changes type is a change" begin + real1, real2 = stored("real"; value_real=1.0), stored("real"; value_real=2.0) + @test compare_values(real1, real1)[3] == "OK" + @test compare_values(real1, real2)[3] == "CHANGED" + @test compare_values(stored("missing"), stored("missing"))[3] == "N/A" + @test compare_values(real1, stored("missing"))[3] == "CHANGED (missing)" + @test compare_values(stored("missing"), real1)[3] == "NEW" + @test compare_values(real1, stored("integer"))[3] == "CHANGED (type)" +end diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index f7fb5ef392..4ed4a278e9 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -87,25 +87,27 @@ Read the GPEC HDF5 kinetic schema. Only `psi` is required; other standard datase optional and surfaced as-is. Non-standard `n_*` datasets become named per-species densities. """ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") - h5open(path, "r") do f - g = group == "/" ? f : f[group] - haskey(g, "psi") || error("kinetic HDF5 '$path' group '$group': missing required dataset 'psi'") - rd(name) = haskey(g, name) ? Float64.(vec(read(g[name]))) : nothing - prov = "" - ats = attributes(g) - haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) - # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density - # profiles for explicit multi-ion input; other names are reserved for future schema fields. - extras = Dict{String,Vector{Float64}}() - for k in keys(g) - (k == "psi" || any(p -> first(p) == k, _KINETIC_H5_PROFILES) || !startswith(k, "n_")) && continue - v = read(g[k]) - v isa AbstractArray && (extras[k] = Float64.(vec(v))) - end - return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))), - (Symbol(name) => rd(name) for (name, _) in _KINETIC_H5_PROFILES)..., - species_densities=(isempty(extras) ? nothing : extras), provenance=prov) + h5open(f -> _read_kinetic_group(group == "/" ? f : f[group], "kinetic HDF5 '$path' group '$group'"), path, "r") +end + +# The kinetic profiles in the open HDF5 file or group `g`; `label` names it in errors. +function _read_kinetic_group(g, label::AbstractString) + haskey(g, "psi") || error("$label: missing required dataset 'psi'") + rd(name) = haskey(g, name) ? Float64.(vec(read(g[name]))) : nothing + prov = "" + ats = attributes(g) + haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) + # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density + # profiles for explicit multi-ion input; other names are reserved for future schema fields. + extras = Dict{String,Vector{Float64}}() + for k in keys(g) + (k == "psi" || any(p -> first(p) == k, _KINETIC_H5_PROFILES) || !startswith(k, "n_")) && continue + v = read(g[k]) + v isa AbstractArray && (extras[k] = Float64.(vec(v))) end + return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))), + (Symbol(name) => rd(name) for (name, _) in _KINETIC_H5_PROFILES)..., + species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end """ @@ -143,30 +145,38 @@ with `read_kinetic_file`. The group carries `schema_version` and `provenance` at function write_kinetic_h5(path::AbstractString, data::KineticProfileData; group::AbstractString="/", schema_version::AbstractString="1.0", provenance::AbstractString=data.provenance) - h5open(path, "w") do f - g = group == "/" ? f : create_group(f, group) - function put(name, v, units) - v === nothing && return + h5open(f -> write_kinetic_group!(group == "/" ? f : create_group(f, group), data; schema_version, provenance), path, "w") + return path +end + +""" + write_kinetic_group!(g, data; schema_version="1.0", provenance=data.provenance) -> g + +Write the profiles of `data` into the open HDF5 file or group `g` in the GPEC kinetic schema, as +[`write_kinetic_h5`](@ref) does for a new file; the composition and closure fields are not written. +""" +function write_kinetic_group!(g, data::KineticProfileData; schema_version::AbstractString="1.0", provenance::AbstractString=data.provenance) + function put(name, v, units) + v === nothing && return + g[name] = collect(Float64, v) + attributes(g[name])["units"] = units + end + put("psi", data.psi, "normalized poloidal flux") + for (name, units) in _KINETIC_H5_PROFILES + put(name, getfield(data, Symbol(name)), units) + end + # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip. + if data.species_densities !== nothing + for (name, v) in data.species_densities + startswith(name, "n_") || + error("species density dataset `$name` must be named `n_*` to be readable as a per-species profile") g[name] = collect(Float64, v) - attributes(g[name])["units"] = units + attributes(g[name])["units"] = "m^-3" end - put("psi", data.psi, "normalized poloidal flux") - for (name, units) in _KINETIC_H5_PROFILES - put(name, getfield(data, Symbol(name)), units) - end - # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip. - if data.species_densities !== nothing - for (name, v) in data.species_densities - startswith(name, "n_") || - error("species density dataset `$name` must be named `n_*` to be readable as a per-species profile") - g[name] = collect(Float64, v) - attributes(g[name])["units"] = "m^-3" - end - end - attributes(g)["schema_version"] = schema_version - attributes(g)["provenance"] = provenance end - return path + attributes(g)["schema_version"] = schema_version + attributes(g)["provenance"] = provenance + return g end """ diff --git a/src/ForceFreeStates/InputSnapshot.jl b/src/ForceFreeStates/InputSnapshot.jl index 0a5d2c3b3b..1aa1c0c864 100644 --- a/src/ForceFreeStates/InputSnapshot.jl +++ b/src/ForceFreeStates/InputSnapshot.jl @@ -13,3 +13,16 @@ function save_equilibrium_ingest!(out_h5, equil) end return nothing end + +""" + save_input_files!(out_h5, equil, wall_settings) + +Record the external files a run read, so a rerun and `equilibrium_from_h5` replay them without the originals: +the equilibrium's kinetic profiles under `Input/RawInputs/Kinetic/` (the kinetic HDF5 schema; the composition is +restated by the deck) and a file wall's text as `Input/RawInputs/Wall/wall_file_raw`. +""" +function save_input_files!(out_h5, equil, wall_settings) + equil.kinetic === nothing || Equilibrium.write_kinetic_group!(create_group(out_h5, "Input/RawInputs/Kinetic"), equil.kinetic) + wall_settings isa Vacuum.FileWall && (out_h5["Input/RawInputs/Wall/wall_file_raw"] = read(wall_settings.file, String)) + return nothing +end diff --git a/src/ForceFreeStates/Output.jl b/src/ForceFreeStates/Output.jl index dd6b9ed612..ba8480148c 100644 --- a/src/ForceFreeStates/Output.jl +++ b/src/ForceFreeStates/Output.jl @@ -52,6 +52,7 @@ function save_to_hdf5( end save_equilibrium_ingest!(out_h5, equil) + save_input_files!(out_h5, equil, result.wall_settings) Utilities.record_inputs!(out_h5, result) # Write derived run parameters diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 7a91500f09..33cadb66be 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -7,8 +7,7 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTup - `name` — short method identifier used as the HDF5 group key and in `intr.method` - `flag` — the `KineticForcesControl` field symbol that enables the method - `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` - (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the - three special-cased methods) + (`:gar` for the GAR/matrix family, `:fcgl` for the fluid CGL method) - `doc` — one-line description printed in verbose output To add a method: append an entry here and add the matching `*_flag` field @@ -18,8 +17,6 @@ const METHOD_REGISTRY = ( (name="fgar", flag=:fgar_flag, kind=:gar, doc="Full general-aspect-ratio calculation"), (name="tgar", flag=:tgar_flag, kind=:gar, doc="Trapped particle general-aspect-ratio calculation"), (name="pgar", flag=:pgar_flag, kind=:gar, doc="Passing particle general-aspect-ratio calculation"), - (name="rlar", flag=:rlar_flag, kind=:rlar, doc="Trapped particle large-aspect-ratio calculation"), - (name="clar", flag=:clar_flag, kind=:clar, doc="Trapped particle cylindrical large-aspect-ratio calculation"), (name="fcgl", flag=:fcgl_flag, kind=:fcgl, doc="Fluid Chew–Goldberger–Low calculation"), (name="fwmm", flag=:fwmm_flag, kind=:gar, doc="Full energy calculation using MXM Euler–Lagrange matrix"), (name="twmm", flag=:twmm_flag, kind=:gar, doc="Trapped energy calculation using MXM Euler–Lagrange matrix"), @@ -37,7 +34,7 @@ const METHOD_REGISTRY = ( """ method_kind(name) -> Symbol -Return the dispatch routing tag (`:gar`, `:fcgl`, `:rlar`, `:clar`) for the NTV +Return the dispatch routing tag (`:gar` or `:fcgl`) for the NTV method `name`, looked up from [`METHOD_REGISTRY`]. Errors on an unknown method. """ function method_kind(name::AbstractString) @@ -80,8 +77,6 @@ Immutable: vary a field by building a new control. fgar_flag::Bool = true # Full general-aspect-ratio tgar_flag::Bool = false # Trapped particle GAR pgar_flag::Bool = false # Passing particle GAR - rlar_flag::Bool = false # Trapped particle large-aspect-ratio - clar_flag::Bool = false # Trapped cylindrical LAR fcgl_flag::Bool = false # Fluid Chew-Goldberger-Low fwmm_flag::Bool = false # Full energy via MXM E-L matrix twmm_flag::Bool = false # Trapped energy via MXM E-L matrix diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index ff3b38aaf1..b33a0ae343 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -16,7 +16,7 @@ Toroidal torque from nonambipolar transport at one ψ and bounce harmonic `l`, w - `wdfac::Float64`: Drift factor - `divxfac::Float64`: Divergence factor - `electron::Bool`: Calculate quantities for electrons (zi,mi ignored) - - `method::String`: Integration method (RLAR, CLAR, *GAR, *TMM, *WMM, *KMM) + - `method::String`: Integration method (FCGL, *GAR, *TMM, *WMM, *KMM) where * = F,T,P for full,trapped,passing - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines - `intr::KineticForcesInternal`: Internal state with mode indexing and perturbation splines @@ -58,10 +58,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, _torque_surface(psi, zi, mi, electron, method, equil, intr, kinetic_profiles; need_matrices=(!isnothing(op_wmats)), n, wdfac, divxfac) : surface - (; chrg, mass, tspl, bmax, bmin, q, n_s, T_s, welec, wdian, wdiat, wtran, wgyro, nuk, epsr, dVdpsi, dbob_m_f, divx_m_f) = surf - wbhat = (π / 4) * sqrt(epsr / 2) * wtran - wdhat = q^3 * wtran^2 / (4 * epsr * wgyro) * wdfac - nueff = nuk / (2 * epsr) + (; chrg, mass, tspl, bmax, bmin, q, n_s, T_s, welec, wdian, wdiat, nuk, epsr, dbob_m_f, divx_m_f) = surf # Method selection — route on the registry dispatch tag (errors on unknown method) kind = method_kind(method) @@ -69,16 +66,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, tpsi_var[] = calculate_fcgl(psi, n, l, tspl, dbob_m_f, divx_m_f, divxfac, n_s, T_s, equil, intr) - elseif kind == :rlar - tpsi_var[] = calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec, - wdhat, wbhat, nueff, dVdpsi, n_s, T_s, - dbob_m_f, intr.bo, bmin) - - elseif kind == :clar - tpsi_var[] = calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, - nuk, intr.bo, bmax, bmin, n_s, T_s, mass, chrg, - tspl, dbob_m_f, divx_m_f, divxfac, wdfac) - else # :gar — fgar/tgar/pgar + all *mm matrix methods tpsi_var[] = calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, intr.bo, bmax, bmin, n_s, T_s, mass, chrg, @@ -167,63 +154,6 @@ function calculate_fcgl(psi, n, l, tspl, dbob_m_f, divx_m_f, divxfac::Float64, end -""" - calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec, wdhat, wbhat, - nueff, dVdpsi, n_s, T_s, dbob_m_f, bo, bmin)::ComplexF64 - -RLAR (Reduced Large Aspect Ratio) torque: energy-space integration with a pitch-averaged -bounce response. Valid for large aspect ratio (ε ≪ 1); the κ bounce average is a placeholder estimate. - -Reference: [Logan et al., Phys. Plasmas, 2013] -""" -function calculate_rlar(psi, n, l, q, epsr, wdian, wdiat, welec, - wdhat, wbhat, nueff, dVdpsi::Float64, - n_s::Float64, T_s::Float64, - dbob_m_f, bo, bmin=0.5)::ComplexF64 - - # Setup parameters for energy integration - lnq = Float64(l) # Resonant mode for trapped particles - - # Maximum pitch angle parameter (use safe estimate) - lmdamax = min(1.0 / (1 - epsr), bo / bmin) - - # Energy-space quadrature ∫ K(x) dx of the resonance operator K - xint = integrate_energy(wdian, wdiat, welec, wdhat, wbhat, nueff, - l, lnq, n, psi, lmdamax, "rlar") - - # Kappa/bounce averaging (effect of field perturbations) - # Placeholder: simplified estimate - kappaint_val = sqrt(mean(abs.(dbob_m_f) .^ 2)) - - # dV/dpsi gradient term from Clebsch coordinate Jacobian - psi_factor = dVdpsi - - # Normalization: √(ε/(2π³)) * n² * n_s * T_s - norm = sqrt(epsr / (2.0 * π^3)) * Float64(n)^2 * n_s * T_s - - # dT/dψ = -0.5 · dV/dψ · κ_int · x_int · norm - result = psi_factor * kappaint_val * 0.5 * (-xint) * norm - - return result -end - - -""" - calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, - bmax, bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f, - divxfac, wdfac)::ComplexF64 - -CLAR (Circular Large Aspect Ratio) torque. Not implemented: warns once and returns zero. -""" -function calculate_clar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, - bmax, bmin, n_s::Float64, T_s::Float64, mass, chrg, tspl, - dbob_m_f, divx_m_f, divxfac, wdfac)::ComplexF64 - - @warn "CLAR method not yet fully implemented, returning zero" maxlog = 1 - return ComplexF64(0.0, 0.0) -end - - """ calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, bmax, bmin, n_s, T_s, mass, chrg, tspl, dbob_m_f, divx_m_f, diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 0f17f7925d..8902a274b1 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -56,7 +56,7 @@ density peaks (Logan & Park, Phys. Plasmas 20, 122507 (2013), §IV–V). `xeval` peak of the Maxwellian-weighted drive `x^2.5·e^−x`. The ℓ = 0 node is the ω_d-shifted ExB (superbanana-plateau) resonance. -The frequencies are the cheap large-aspect-ratio closed forms of the `rlar` method, +The frequencies are the cheap large-aspect-ratio closed forms, `ω_b = (π/4)·√(ε/2)·wtran` and `ω_d = q·T_s/(2·ε·R₀²·Z·e·B₀)·wdfac`, with ε = ⟨r⟩/⟨R⟩ clamped away from the axis; no bounce averaging, which is enough for panel placement. """ diff --git a/src/Rerun.jl b/src/Rerun.jl index 4363de63f1..0f8de2f56d 100644 --- a/src/Rerun.jl +++ b/src/Rerun.jl @@ -142,7 +142,9 @@ end Decide where the rerun output should go. Default naming rule: take the source basename, drop `.h5`, and append `_rerun.h5`. Refuse to overwrite the source -file if the resolved destination collides. +file if the output or one of its side files (`_RERUN_SIDE_FILES`) is the +source, by path or as the same file under another name (a link, or a case-only +difference on a case-insensitive file system). """ function resolve_rerun_output_path(source_h5::String, output_dir::String, output_name::Union{Nothing,String}) source_abs = abspath(source_h5) @@ -151,13 +153,48 @@ function resolve_rerun_output_path(source_h5::String, output_dir::String, output stem = endswith(lowercase(base), ".h5") ? base[1:end-3] : base output_name = string(stem, "_rerun.h5") end - if abspath(joinpath(output_dir, output_name)) == source_abs - error("Refusing to overwrite source HDF5 at $source_abs. " * - "Pass `--output-dir ` or `--output-name `.") + for name in (output_name, (_rerun_side_file(output_name, s) for s in _RERUN_SIDE_FILES)...) + dest = abspath(joinpath(output_dir, name)) + if dest == source_abs || (isfile(dest) && samefile(dest, source_abs)) + error("Refusing to overwrite source HDF5 at $source_abs with $dest. " * + "Pass `--output-dir ` or `--output-name `.") + end end return (output_dir, output_name) end +# Suffixes of the input files a rerun writes next to its output from the source snapshot. +const _RERUN_SIDE_FILES = (forcing="_forcing.h5", tolerances="_tolerances.toml", kinetic="_kinetic.h5", wall="_wall.dat") + +# The side file a rerun writes next to `output_name` for one of the `_RERUN_SIDE_FILES` suffixes. +_rerun_side_file(output_name::AbstractString, suffix::AbstractString) = splitext(output_name)[1] * suffix + +# Where a run records the external files it read (written by ForceFreeStates.save_input_files!). +const _RECORDED_KINETIC = "Input/RawInputs/Kinetic" +const _RECORDED_WALL = "Input/RawInputs/Wall/wall_file_raw" + +# Point the stored deck's kinetic_file and file wall at copies written from the source snapshot. A deck that names +# either without its snapshot comes from an older GPEC and cannot be replayed faithfully. +function _replay_input_files!(inputs, source_h5, output_dir, output_name, has_kinetic::Bool, wall_raw) + remedy = "re-run from the original deck instead of replaying it" + equil = inputs["Equilibrium"] + if has_kinetic + kinetic_file = _rerun_side_file(output_name, _RERUN_SIDE_FILES.kinetic) + Equilibrium.write_kinetic_h5(joinpath(output_dir, kinetic_file), h5open(f -> Equilibrium._read_kinetic_group(f[_RECORDED_KINETIC], source_h5), source_h5)) + equil["kinetic_file"], equil["kinetic_group"] = kinetic_file, "/" + elseif !isempty(get(equil, "kinetic_file", "")) + error("$source_h5 names a kinetic_file but records no $_RECORDED_KINETIC; $remedy") + end + wall = get(inputs, "Wall", nothing) + if wall !== nothing && get(wall, "shape", "") == "file" + wall_raw === nothing && error("$source_h5 has a file wall but records no $_RECORDED_WALL; $remedy") + wall_file = _rerun_side_file(output_name, _RERUN_SIDE_FILES.wall) + write(joinpath(output_dir, wall_file), wall_raw) + wall["file"] = wall_file + end + return inputs +end + """ build_inputs_from_h5(args::Vector{String}) -> (inputs, eq_config, additional_input, output_dir, current_git, preloaded_coils) @@ -170,6 +207,8 @@ The recorded forcing (`Input/RawInputs/ForcingTerms/`) is written next to the re `_forcing.h5` and `[ForcingTerms]` is pointed at it before overrides apply, so the replay reads and normalizes it against its own equilibrium exactly as a deck run reads its file; a recorded tolerance file is written as `_tolerances.toml` and `[ErrorFields] tolerance_file` points at it. +The recorded kinetic profiles and file wall are written as `_kinetic.h5` and `_wall.dat` the same +way, so a replay never reads the files the source run named. """ function build_inputs_from_h5(args::Vector{String}) cli = parse_rerun_cli(args) @@ -178,7 +217,7 @@ function build_inputs_from_h5(args::Vector{String}) # Stored TOML, equilibrium ingest and the forcing snapshot (present only with [PerturbedEquilibrium]). use_coils = cli.coil_source == "coils" - toml_raw, ingest, source_git, forcing_snapshot, preloaded_coils, tolerance_raw = h5open(source_h5, "r") do in_h5 + toml_raw, ingest, source_git, forcing_snapshot, preloaded_coils, tolerance_raw, has_kinetic, wall_raw = h5open(source_h5, "r") do in_h5 haskey(in_h5, "Input/gpec_toml_raw") || error("Source HDF5 $source_h5 has no Input/gpec_toml_raw — produced by a pre-rerun version of GPEC") forcing_modes = if !haskey(in_h5, "Input/RawInputs/ForcingTerms") @@ -204,7 +243,9 @@ function build_inputs_from_h5(args::Vector{String}) haskey(in_h5, "Info/git_version") ? read(in_h5, "Info/git_version") : "unknown", forcing_modes, coil_sets, - haskey(in_h5, _RECORDED_TOLERANCES) ? read(in_h5, _RECORDED_TOLERANCES) : nothing + haskey(in_h5, _RECORDED_TOLERANCES) ? read(in_h5, _RECORDED_TOLERANCES) : nothing, + haskey(in_h5, _RECORDED_KINETIC), + haskey(in_h5, _RECORDED_WALL) ? read(in_h5, _RECORDED_WALL) : nothing ) end @@ -218,7 +259,7 @@ function build_inputs_from_h5(args::Vector{String}) _stop_on_deprecated_deck_keys(inputs; remedy="$source_h5 was written by an older GPEC version; re-run from the original deck instead of replaying it") if forcing_snapshot !== nothing - forcing_file = splitext(output_name)[1] * "_forcing.h5" + forcing_file = _rerun_side_file(output_name, _RERUN_SIDE_FILES.forcing) modes, normalization = forcing_snapshot h5open(f -> ForcingTerms.save_forcing_to_h5(modes, f; normalization), joinpath(output_dir, forcing_file), "w") ft = get!(inputs, "ForcingTerms", Dict{String,Any}()) @@ -226,10 +267,11 @@ function build_inputs_from_h5(args::Vector{String}) ft["forcing_data_format"] = "hdf5" end if tolerance_raw !== nothing && haskey(inputs, "ErrorFields") - tolerance_file = splitext(output_name)[1] * "_tolerances.toml" + tolerance_file = _rerun_side_file(output_name, _RERUN_SIDE_FILES.tolerances) write(joinpath(output_dir, tolerance_file), tolerance_raw) inputs["ErrorFields"]["tolerance_file"] = tolerance_file end + _replay_input_files!(inputs, source_h5, output_dir, output_name, has_kinetic, wall_raw) if cli.override_file !== nothing isfile(cli.override_file) || error("Override file not found: $(cli.override_file)") file_overrides = TOML.parsefile(cli.override_file) @@ -288,12 +330,14 @@ function rebuild_equilibrium_inputs(inputs::Dict{String,Any}, ingest, output_dir return eq_config, additional_input end -# The stored deck's kinetic profiles, re-read from its file; a file that has since moved leaves the equilibrium without them. -function _replay_kinetic_profiles(config::Equilibrium.EquilibriumConfig) - isempty(config.kinetic_file) && return nothing - isfile(config.kinetic_file) && return Equilibrium.configured_kinetic_profiles(config) - @warn "The stored deck's kinetic_file is not readable; the restored equilibrium carries no kinetic profiles" file = config.kinetic_file - return nothing +# The recorded kinetic profiles with the stored deck's composition; `nothing` when the run had none. +function _replay_kinetic_profiles(f, config::Equilibrium.EquilibriumConfig, h5path) + if !haskey(f, _RECORDED_KINETIC) + isempty(config.kinetic_file) && return nothing + error("$h5path names a kinetic_file but records no $_RECORDED_KINETIC; re-run from the original deck") + end + raw = Equilibrium._read_kinetic_group(f[_RECORDED_KINETIC], "$h5path $_RECORDED_KINETIC") + return Equilibrium.KineticProfileData(raw; Equilibrium.kinetic_composition(config)...) end """ @@ -337,7 +381,7 @@ function equilibrium_from_h5(h5path::AbstractString) eqfun = Equilibrium._eqfun_from_geometry(profiles, rzphi..., params.ro, params.psio, grid2d, opts2d) geometry = Equilibrium.compute_geometry_profiles(xs, ys, rzphi[1], rzphi[2], rzphi[4], params.ro) equil = Equilibrium.PlasmaEquilibrium(config, params, profiles, geometry, xs, ys, rzphi..., eqfun..., params.ro, params.zo, params.psio, - read_equilibrium_ingest(f), analytic_config(inputs), _replay_kinetic_profiles(config)) + read_equilibrium_ingest(f), analytic_config(inputs), _replay_kinetic_profiles(f, config, h5path)) return (; equil, inputs, psilim=Float64(read(f, "Info/psilim"))) end end diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl index c668dbf21b..345671fe3f 100644 --- a/src/Vacuum/Kernel3D.jl +++ b/src/Vacuum/Kernel3D.jl @@ -163,22 +163,27 @@ end # Cached singular quadrature; rebuilt if PATCH_RAD / RAD_DIM / INTERP_ORDER change const SINGULAR_QUAD_CACHE = Ref{Union{Nothing,SingularQuadratureData}}(nothing) +const _SINGULAR_QUAD_LOCK = ReentrantLock() """ get_singular_quadrature(PATCH_RAD, RAD_DIM, INTERP_ORDER) -Return the cached [`SingularQuadratureData`](@ref), rebuilding it if the parameters changed. +Return the cached [`SingularQuadratureData`](@ref), rebuilding it if the parameters changed. Takes a lock, +so concurrent solves on different grids each get the data for their own parameters. """ function get_singular_quadrature(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) - cached = SINGULAR_QUAD_CACHE[] - if !isnothing(cached) && - cached.PATCH_RAD == PATCH_RAD && - cached.RAD_DIM == RAD_DIM && - cached.INTERP_ORDER == INTERP_ORDER - return cached + @lock _SINGULAR_QUAD_LOCK begin + cached = SINGULAR_QUAD_CACHE[] + if !isnothing(cached) && + cached.PATCH_RAD == PATCH_RAD && + cached.RAD_DIM == RAD_DIM && + cached.INTERP_ORDER == INTERP_ORDER + return cached + end + data = SingularQuadratureData(PATCH_RAD, RAD_DIM, INTERP_ORDER) + SINGULAR_QUAD_CACHE[] = data + return data end - SINGULAR_QUAD_CACHE[] = SingularQuadratureData(PATCH_RAD, RAD_DIM, INTERP_ORDER) - return SINGULAR_QUAD_CACHE[] end """ diff --git a/src/Vacuum/PnQuadCache.jl b/src/Vacuum/PnQuadCache.jl index d4a40ad868..0bb58ca23d 100644 --- a/src/Vacuum/PnQuadCache.jl +++ b/src/Vacuum/PnQuadCache.jl @@ -46,7 +46,7 @@ into the package image. """ function reset_caches!() @lock _PN_CACHE_LOCK empty!(_PN_CACHE) - SINGULAR_QUAD_CACHE[] = nothing + @lock _SINGULAR_QUAD_LOCK SINGULAR_QUAD_CACHE[] = nothing return nothing end diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index 8f2f1769e1..8dddc9b61e 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -498,6 +498,9 @@ include("h5_metadata_check.jl") rebuilt = GPEC.equilibrium_from_h5(h5path) @test rebuilt.psilim == ffs.psilim @test rebuilt.equil.params.bt0 ≈ ffs.equil.params.bt0 + # The run records its kinetic profiles, and the rebuilt equilibrium carries them. + @test h5open(f -> haskey(f, "Input/RawInputs/Kinetic"), h5path) + @test rebuilt.equil.kinetic.n_e == ffs.equil.kinetic.n_e post_hoc = EF.compute_coil_sensitivities(h5path, sets[1:1]) @test post_hoc.field_as_designed ≈ sens.field_as_designed[:, 1:1] rtol = 1e-10 @test post_hoc.shift_sensitivity_per_m ≈ sens.shift_sensitivity_per_m[:, :, 1:1] rtol = 1e-10 diff --git a/test/runtests_frozen_equilibrium.jl b/test/runtests_frozen_equilibrium.jl index e6814cb165..6a5061af07 100644 --- a/test/runtests_frozen_equilibrium.jl +++ b/test/runtests_frozen_equilibrium.jl @@ -115,20 +115,19 @@ using GeneralizedPerturbedEquilibrium @test analytic_snapshot.analytic isa Eq.SolovevConfig @test analytic_snapshot.analytic.elongation == 1.8 @test analytic_snapshot.profiles.q_spline.y == profile_values[4] - # The stored deck's kinetic_file resolves against the h5's directory and is re-read; a missing file warns and leaves none. - kinetic_path = joinpath(dir, "kinetic.dat") - cp(joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example", "kinetic.dat"), kinetic_path) + # The kinetic profiles come from the recorded snapshot, never from the deck's kinetic_file (absent here); a deck + # naming one without the snapshot cannot replay. analytic_inputs["Equilibrium"]["kinetic_file"] = "kinetic.dat" h5open(path, "r+") do f delete_object(f, "Input/gpec_toml_raw") f["Input/gpec_toml_raw"] = sprint(io -> TOML.print(io, analytic_inputs)) end + @test_throws r"records no Input/RawInputs/Kinetic" GPEC.equilibrium_from_h5(path) + kinetic = Eq.read_kinetic_file(joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example", "kinetic.dat")) + h5open(f -> Eq.write_kinetic_group!(create_group(f, "Input/RawInputs/Kinetic"), kinetic), path, "r+") with_kinetic = GPEC.equilibrium_from_h5(path).equil - @test with_kinetic.kinetic.n_e == Eq.read_kinetic_file(kinetic_path).n_e - @test with_kinetic.config.kinetic_file == kinetic_path - rm(kinetic_path) - without = @test_logs (:warn, r"kinetic_file is not readable") match_mode = :any GPEC.equilibrium_from_h5(path).equil - @test without.kinetic === nothing + @test with_kinetic.kinetic.n_e == kinetic.n_e + @test with_kinetic.kinetic.omega_E == kinetic.omega_E h5open(path, "r+") do f delete_object(f, "Equilibrium/Profiles/q") end diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 1b03704e1f..b97cbb48ad 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -722,7 +722,7 @@ using TOML for entry in KF.METHOD_REGISTRY @test entry.flag in fieldnames(KF.KineticForcesControl) @test endswith(string(entry.flag), "_flag") - @test entry.kind in (:gar, :fcgl, :rlar, :clar) + @test entry.kind in (:gar, :fcgl) @test KF.method_kind(entry.name) == entry.kind end end diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 952a010058..bb9ce7befb 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -133,6 +133,41 @@ end # Default output is `_rerun.h5`, so clobbering only happens if the user # routes the rerun back onto the source file; exercise that guard. @test_throws ErrorException GeneralizedPerturbedEquilibrium.main([source_h5, "--output-dir", snapshot_dir, "--output-name", "gpec.h5"]) + # The same file under another name, as the output or as one of its side files, is refused too. + mktempdir() do out_dir + symlink(source_h5, joinpath(out_dir, "alias.h5")) + @test_throws r"Refusing to overwrite" GeneralizedPerturbedEquilibrium.main([source_h5, "--output-dir", out_dir, "--output-name", "alias.h5"]) + symlink(source_h5, joinpath(out_dir, "run_forcing.h5")) + @test_throws r"Refusing to overwrite" GeneralizedPerturbedEquilibrium.main([source_h5, "--output-dir", out_dir, "--output-name", "run.h5"]) + end + end + + @testset "kinetic profiles and a file wall replay from the snapshot" begin + mktempdir() do out_dir + doctored = joinpath(out_dir, "source.h5") + cp(source_h5, doctored) + kinetic = GeneralizedPerturbedEquilibrium.Equilibrium.read_kinetic_file(joinpath(@__DIR__, "..", "examples", "TJ_kinetic_NTV_example", "kinetic.dat")) + wall_text = "3\n1.0\nindex R Z\n1 2.0 0.0\n2 1.5 0.5\n3 1.5 -0.5\n" + h5open(doctored, "r+") do h5 + inputs = TOML.parse(read(h5, "Input/gpec_toml_raw")) + # Neither file exists: the replay must take both from the snapshot. + inputs["Equilibrium"]["kinetic_file"] = joinpath(out_dir, "gone.dat") + inputs["Wall"] = Dict{String,Any}("shape" => "file", "file" => joinpath(out_dir, "gone_wall.dat")) + delete_object(h5, "Input/gpec_toml_raw") + h5["Input/gpec_toml_raw"] = sprint(TOML.print, inputs) + end + @test_throws r"records no Input/RawInputs/Kinetic" GeneralizedPerturbedEquilibrium.build_inputs_from_h5([doctored, "--output-dir", out_dir]) + h5open(doctored, "r+") do h5 + GeneralizedPerturbedEquilibrium.Equilibrium.write_kinetic_group!(create_group(h5, "Input/RawInputs/Kinetic"), kinetic) + end + @test_throws r"records no Input/RawInputs/Wall" GeneralizedPerturbedEquilibrium.build_inputs_from_h5([doctored, "--output-dir", out_dir]) + h5open(h5 -> h5["Input/RawInputs/Wall/wall_file_raw"] = wall_text, doctored, "r+") + inputs = first(GeneralizedPerturbedEquilibrium.build_inputs_from_h5([doctored, "--output-dir", out_dir])) + @test inputs["Equilibrium"]["kinetic_file"] == "source_rerun_kinetic.h5" + @test GeneralizedPerturbedEquilibrium.Equilibrium.read_kinetic_file(joinpath(out_dir, "source_rerun_kinetic.h5")).n_e == kinetic.n_e + @test inputs["Wall"]["file"] == "source_rerun_wall.dat" + @test read(joinpath(out_dir, "source_rerun_wall.dat"), String) == wall_text + end end @testset "--coil-source coils errors when source has no coil snapshot" begin