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 6743b0f19c..0000000000 Binary files a/ci/conventions/__pycache__/check_subject.cpython-310.pyc and /dev/null differ 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 b227f876b6..8cb343f7dd 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -120,25 +120,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 """ @@ -177,30 +179,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/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 3f41526944..1e3347dd1e 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -502,6 +502,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_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