diff --git a/Project.toml b/Project.toml index 4c1c7c78..b084c60c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,39 +1,69 @@ name = "JPEC" uuid = "462872dd-e066-4d2e-b993-6468b5239634" license = "MIT" -authors = ["Matthew Pharr ", "Rithik Banerjee ", "Jaebeom Cho ", "Jacob Halpern "] version = "0.1.0" +authors = ["Matthew Pharr ", "Rithik Banerjee ", "Jaebeom Cho ", "Jacob Halpern "] [deps] +AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b" +DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" +EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +FastInterpolations = "9ea80cae-fc13-4c00-8066-6eaedb12f34b" +FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" +IMASdd = "c5a45a97-b3f9-491c-b9a7-aa88c3bc0067" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +AbstractTrees = "0.4.5" +BenchmarkTools = "1.6.3" +DataFrames = "1.8.1" DiffEqCallbacks = "4.9.0" +Distributions = "0.25.123" Documenter = "1.14.1" -FastInterpolations = ">=0.2.8" +DynamicalSystems = "3.6.7" +EFIT = "1.2.5" FFTW = "1.9.0" +FastInterpolations = ">=0.2.8" +FileIO = "1.18.0" HDF5 = "0.17.2" +IJulia = "1.34.2" +IMASdd = ">=8" JLD2 = "0.6.3" +LaTeXStrings = "1.4.0" LinearAlgebra = "1.11.0" +Measurements = "2.14.1" OrdinaryDiffEq = "6.102.0" Pkg = "1.11.0" Plots = "1.40.15" Printf = "1.11.0" +PyPlot = "2.11.6" SpecialFunctions = "2.5.1" StaticArrays = "1.9.15" TOML = "1.0.3" Test = "1.11.0" + +[extras] +EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" + +[targets] +test = ["EFIT"] diff --git a/src/DCON/DCON.jl b/src/DCON/DCON.jl index d32b4e19..058c58e4 100644 --- a/src/DCON/DCON.jl +++ b/src/DCON/DCON.jl @@ -16,6 +16,7 @@ import ..Util import ..Vacuum using Printf import StaticArrays: @MVector, @MMatrix +import IMASdd # Include all necessary files include("DconStructs.jl") @@ -28,6 +29,7 @@ include("Fourfit.jl") include("FixedBoundaryStability.jl") include("Utils.jl") include("Free.jl") +include("WriteImas.jl") # These are used for various small tolerances and root finders throughout DCON global eps = 1e-10 diff --git a/src/DCON/Main.jl b/src/DCON/Main.jl index 4b1513ad..cd57786d 100644 --- a/src/DCON/Main.jl +++ b/src/DCON/Main.jl @@ -1,4 +1,22 @@ +""" + Main(path::String="./") + +Run DCON from a directory containing `dcon.toml` and `equil.toml`. +""" function Main(path::String="./") + return Main(path, nothing) +end + +""" + Main(path::String, dd::IMASdd.dd) + +Run DCON using IMAS data dictionary `dd` for the equilibrium instead of an +`equil.toml` file. The directory `path` must still contain `dcon.toml` for +DCON control parameters. Equilibrium resolution settings (mpsi, mtheta, etc.) +use their defaults; override by constructing your own `EquilibriumConfig` and +calling `Equilibrium.setup_equilibrium` directly if needed. +""" +function Main(path::String, dd::Union{IMASdd.dd, Nothing}) println("DCON START") println("----------------------------------") @@ -9,7 +27,21 @@ function Main(path::String="./") # TODO: leaving DCON_CONTROL as a part of the toml file, eventually can combine equil, gpec, etc. into one input file? inputs = TOML.parsefile(joinpath(intr.dir_path, "dcon.toml")) ctrl = DconControl(; (Symbol(k) => v for (k, v) in inputs["DCON_CONTROL"])...) - equil = Equilibrium.setup_equilibrium(joinpath(intr.dir_path, "equil.toml")) + + if dd === nothing + # Standard file-based equilibrium: read equil.toml from the run directory + equil = Equilibrium.setup_equilibrium(joinpath(intr.dir_path, "equil.toml")) + else + # IMAS-based equilibrium: build config with defaults; data comes from dd. + # eq_filename is a placeholder — its directory is used for any diagnostic outputs. + eq_config = Equilibrium.EquilibriumConfig(; + control = Equilibrium.EquilibriumControl(; + eq_type = "imas", + eq_filename = joinpath(path, "imas_equilibrium"), + ) + ) + equil = Equilibrium.setup_equilibrium(eq_config, dd) + end if "WALL" in keys(inputs) intr.wall_settings = Vacuum.WallShapeSettings(; (Symbol(k) => v for (k, v) in inputs["WALL"])...) else @@ -178,6 +210,17 @@ function Main(path::String="./") write_outputs_to_HDF5(ctrl, equil, intr, odet, ctrl.vac_flag ? vac_data : nothing) end + result = (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, vac_data=ctrl.vac_flag ? vac_data : nothing) + + # Write linear stability results back to the IMAS data dictionary (Task 5). + # Only done when dd was supplied by the caller (IMAS workflow). + if dd !== nothing + write_imas(dd, result) + if ctrl.verbose + println("Stability results written to dd.mhd_linear ✓") + end + end + end_time = time() - start_time println("----------------------------------") println("Run time: $(@sprintf("%.3e", end_time)) seconds") @@ -185,7 +228,7 @@ function Main(path::String="./") # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found - return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet, vac_data=ctrl.vac_flag ? vac_data : nothing) + return result end """ diff --git a/src/DCON/WriteImas.jl b/src/DCON/WriteImas.jl new file mode 100644 index 00000000..a523161c --- /dev/null +++ b/src/DCON/WriteImas.jl @@ -0,0 +1,77 @@ +""" + write_imas(dd::IMASdd.dd, result::NamedTuple) -> IMASdd.dd + +Write JPEC/DCON linear stability results into the `mhd_linear` IDS of an IMAS data +dictionary `dd`. + +## What gets written + +One `time_slice` is created (time taken from `dd.equilibrium.time_slice[1]`). +Inside that slice, one `toroidal_mode` entry is written **per computed eigenmode** +(there are `numpert_total = mpert × npert` eigenmodes in total). + +For each eigenmode `i`: +- `n_tor` — toroidal mode number, assigned by block index `nlow + (i-1) ÷ mpert`. + Exact for single-n runs (`npert == 1`); approximate for multi-n + runs because the global eigenvalue sort breaks strict block ordering. +- `energy_perturbed` — real part of the `i`-th total (plasma+vacuum) energy eigenvalue + `real(et[i])` in DCON's internal normalisation units. + Negative value → mode is MHD-unstable (δW < 0 criterion). + Only written when `vac_flag = true` and vacuum data is available. + +## Top-level metadata written to `dd.mhd_linear` +- `ideal_flag = 1` (ideal MHD) +- `ids_properties.homogeneous_time = 1` +- `ids_properties.comment` (code name + run range) +- `code.name = "JPEC"` +""" +function write_imas(dd::IMASdd.dd, result::NamedTuple) + + ctrl = result.ctrl + intr = result.intr + vac_data = result.vac_data + + # ── Time: take from the equilibrium already stored in dd ───────────────── + t = dd.equilibrium.time_slice[1].time + + # ── Top-level mhd_linear metadata ──────────────────────────────────────── + mhd = dd.mhd_linear + + mhd.ids_properties.comment = "JPEC DCON ideal MHD linear stability, n = $(intr.nlow):$(intr.nhigh)" + mhd.ids_properties.homogeneous_time = 1 + mhd.ideal_flag = 1 + mhd.code.name = "JPEC" + mhd.time = [t] + + # ── Create one time slice ───────────────────────────────────────────────── + resize!(mhd.time_slice, 1) + ts = mhd.time_slice[1] + ts.time = t + + # ── One toroidal_mode entry per eigenmode ───────────────────────────────── + # numpert_total = mpert × npert eigenmodes were computed. For each one we + # record the toroidal mode number (from its block position) and, when the + # vacuum calculation was run, the total perturbed energy eigenvalue. + resize!(ts.toroidal_mode, intr.numpert_total) + + for i in 1:intr.numpert_total + + tm = ts.toroidal_mode[i] + + # Block-position estimate of the toroidal mode number. + # For npert == 1 this is exact; for npert > 1 it is approximate because + # the global eigenvalue sort may interleave modes from different n-blocks. + ipert_n = (i - 1) ÷ intr.mpert + 1 + tm.n_tor = intr.nlow + ipert_n - 1 + + # Perturbed energy: real part of total (plasma + vacuum) energy eigenvalue. + # Eigenvalues are sorted ascending — et[1] is the most unstable (most negative). + # Units are DCON-internal (normalised by kinetic energy / dV/dψ at the edge). + if ctrl.vac_flag && vac_data !== nothing + tm.energy_perturbed = real(vac_data.et[i]) + end + + end + + return dd +end diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 4376be25..abb470b4 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -6,6 +6,7 @@ import ..Spl using Printf, OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra, HDF5 using TOML import FastInterpolations +import IMASdd using FastInterpolations: cubic_interp, deriv1, deriv2, deriv3, LinearBinary, CubicFit import StaticArrays: @MMatrix, SVector @@ -66,6 +67,20 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin end eq_input = sol_run(eq_config, additional_input) + elseif eq_type == "imas" + + # For IMAS input, additional_input must be a populated dd (data dictionary). + # There is no file to read — the equilibrium data comes directly from dd. + # Usage: + # config = EquilibriumConfig(eq_type="imas", eq_filename="N/A") + # dd = IMAS.json2imas("myshot.json") + # plasma_eq = setup_equilibrium(config, dd) + if additional_input === nothing + error("eq_type=\"imas\" requires a dd object passed as additional_input.\n" * + "Usage: setup_equilibrium(config, dd)") + end + + eq_input = read_imas(eq_config, additional_input) else error("Equilibrium type $(equil_in.eq_type) is not implemented") end diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index a0309d6d..f0e8e3f1 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -379,3 +379,177 @@ function read_chease_ascii(config::EquilibriumConfig) println(" Magnetic axis at (ro=$ro, zo=$zo), psio=$psio") return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio) end + +""" + read_imas(config, dd) + +Reads equilibrium data from an IMAS data dictionary (`dd`), builds the same +1D and 2D splines that the direct solver expects, and returns a `DirectRunInput`. + +This function is the IMAS equivalent of `read_efit`: instead of parsing a text +file (g-file), it pulls data directly from the standardized IMAS data structure. +The output format is identical so the rest of the solver chain is unchanged. + +## Why DirectRunInput and not InverseRunInput? + +IMAS equilibrium data always contains a 2D ψ(R,Z) map (just like an EFIT g-file). +The direct solver starts from that map and integrates along field lines to build +flux coordinates. The inverse solver (CHEASE path) starts from R(ψ,θ), Z(ψ,θ) +which IMAS does not provide in a convenient form. + +## Arguments + + - `config`: The `EquilibriumConfig` object containing grid and solver settings. + - `dd`: An IMAS data dictionary. Must have `dd.equilibrium` populated with at + least one time slice containing `profiles_1d` and `profiles_2d` data. + +## Returns + + - A `DirectRunInput` object ready to be passed to `equilibrium_solver`. + +## IMAS fields used + + - `dd.equilibrium.time_slice[].global_quantities.psi_axis` — ψ at magnetic axis [Wb/rad] + - `dd.equilibrium.time_slice[].global_quantities.psi_boundary` — ψ at LCFS [Wb/rad] + - `dd.equilibrium.time_slice[].profiles_1d.psi` — 1D ψ grid [Wb/rad] + - `dd.equilibrium.time_slice[].profiles_1d.f` — F = R·Bt [T·m] + - `dd.equilibrium.time_slice[].profiles_1d.pressure` — plasma pressure [Pa] + - `dd.equilibrium.time_slice[].profiles_1d.q` — safety factor + - `dd.equilibrium.time_slice[].profiles_2d[1].grid.dim1` — R grid [m] + - `dd.equilibrium.time_slice[].profiles_2d[1].grid.dim2` — Z grid [m] + - `dd.equilibrium.time_slice[].profiles_2d[1].psi` — 2D ψ(R,Z) [Wb/rad] +""" +function read_imas(config::EquilibriumConfig, dd) + println("--> Processing IMAS equilibrium at global_time = $(dd.global_time) s") + + # ------------------------------------------------------------------ + # BLOCK 1: Get the equilibrium time slice + # ------------------------------------------------------------------ + # The [] syntax means "at the current dd.global_time" — the standard + # IMAS way to access time-dependent arrays of structures. + # All equilibrium data for this moment in time lives inside eqt. + eqt = dd.equilibrium.time_slice[] + + # ------------------------------------------------------------------ + # BLOCK 2: Global flux values — the fundamental scale of the equilibrium + # ------------------------------------------------------------------ + # IMAS (COCOS 11) stores ψ as 2π × the value that JPEC's internal solver + # expects (COCOS 2 convention, matching the raw gEQDSK / read_efit units). + # We therefore divide all ψ values by 2π before handing them to the solver, + # and multiply q by 2π for the same reason (q_COCOS11 = q_COCOS2 / 2π). + cocos11_to_internal = 1.0 / (2π) + + # psi_axis: ψ at the magnetic axis (O-point) [Wb, JPEC internal] + # psi_boundary: ψ at the last closed flux surface (LCFS) [Wb, JPEC internal] + # psio: total flux swing = |ψ_axis - ψ_boundary| [Wb] + # + # psio is the normalization scale used throughout the solver. + # Every normalized flux coordinate psi_norm ∈ [0,1] is built from psio. + psi_axis = eqt.global_quantities.psi_axis * cocos11_to_internal + psi_boundary = eqt.global_quantities.psi_boundary * cocos11_to_internal + psio = abs(psi_boundary - psi_axis) + + if psio < 1e-10 + error("read_imas: |psi_axis - psi_boundary| = $psio is too small. " * + "Check that dd.equilibrium is properly populated.") + end + + # ------------------------------------------------------------------ + # BLOCK 3: Extract and normalize the 1D profiles + # ------------------------------------------------------------------ + # IMAS profiles_1d stores quantities as 1D arrays indexed by ψ, + # running from the axis (ψ_axis) to the boundary (ψ_boundary). + # + # The solver's sq_in spline uses psi_norm ∈ [0,1]: + # psi_norm = 0 at the magnetic axis + # psi_norm = 1 at the plasma boundary + # + # So we compute: psi_norm = (ψ - ψ_axis) / (ψ_boundary - ψ_axis) + psi_1d = eqt.profiles_1d.psi .* cocos11_to_internal # ψ [Wb, JPEC internal] + f_1d = eqt.profiles_1d.f # F(ψ) = R·Bt [T·m], COCOS-independent + p_1d = eqt.profiles_1d.pressure # plasma pressure P(ψ) [Pa], COCOS-independent + q_1d = eqt.profiles_1d.q ./ cocos11_to_internal # q [JPEC internal] = q_IMAS × 2π + + nw = length(psi_1d) + psi_norm_grid = (psi_1d .- psi_axis) ./ (psi_boundary - psi_axis) + # clamp protects against tiny floating-point overshoot past [0, 1] + psi_norm_grid = clamp.(psi_norm_grid, 0.0, 1.0) + + # Build the 4-column table that sq_in expects — same format as read_efit: + # column 1: |F| toroidal flux function [T·m] + # abs() enforces the sign convention used by the solver + # column 2: μ₀·P pressure in magnetic units [T²] + # IMAS gives P in [Pa]; μ₀ = 4π×10⁻⁷ [T²/Pa] converts to [T²] + # max(..., 0) prevents unphysical negative pressure from spline overshoot + # column 3: q safety factor [dimensionless] + # column 4: √(psi_norm) square root of normalized flux + sq_fs_nodes = hcat( + abs.(f_1d), + max.(p_1d .* mu0, 0.0), + q_1d, + sqrt.(psi_norm_grid) + ) + sq_in = cubic_interp(psi_norm_grid, sq_fs_nodes; bc=CubicFit(), extrap=:extension) + + # ------------------------------------------------------------------ + # BLOCK 4: Extract the 2D ψ(R,Z) map and apply the solver's convention + # ------------------------------------------------------------------ + # IMAS stores the full 2D poloidal flux map on a rectangular (R,Z) grid + # in profiles_2d[1] (grid_type.index = 1 in the IMAS standard). + # + # Grid layout: + # dim1 → R values [m] (first array dimension) + # dim2 → Z values [m] (second array dimension) + # psi → ψ(R,Z) [Wb/rad], array of shape [nR × nZ] + if isempty(eqt.profiles_2d) + error("read_imas: no profiles_2d found in equilibrium time slice. " * + "Ensure the 2D ψ(R,Z) map is stored in dd.equilibrium.") + end + prof2d = eqt.profiles_2d[1] + r_grid = prof2d.grid.dim1 # R coordinates [m], 1D array, length nR + z_grid = prof2d.grid.dim2 # Z coordinates [m], 1D array, length nZ + psi_rz = prof2d.psi .* cocos11_to_internal # ψ(R,Z) [Wb, JPEC internal] + + # The direct solver's psi_in convention (see DirectRunInput docstring): + # + # psi_in(R,Z) = ψ_boundary - ψ(R,Z) + # + # This makes: + # psi_in = 0 at the plasma boundary (ψ = ψ_boundary) + # psi_in = psio at the magnetic axis (ψ = ψ_axis) + # + # The solver then internally computes psi_norm = 1 - psi_in/psio, + # giving psi_norm = 0 at the axis and 1 at the boundary. + # + # If ψ_boundary < ψ_axis (ψ decreases outward, common in many machines), + # then ψ_boundary - ψ(R,Z) is negative at the axis, so we flip the sign + # to ensure the axis value is positive — matching the solver's expectation. + psi_proc = psi_boundary .- psi_rz + if psi_boundary - psi_axis < 0.0 + psi_proc .*= -1.0 + end + + rmin, rmax = extrema(r_grid) + zmin, zmax = extrema(z_grid) + + psi_in_xs = collect(r_grid) + psi_in_ys = collect(z_grid) + + # CubicFit boundary conditions, no periodic BC (open rectangular grid). + # Omitting `search` uses the same default as read_efit, keeping the + # DirectRunInput type consistent across all input methods. + psi_in = cubic_interp( + (psi_in_xs, psi_in_ys), psi_proc; + bc = (CubicFit(), CubicFit()), + extrap = (:extension, :extension) + ) + + println("--> IMAS equilibrium loaded:") + println(" psio = $(round(psio; sigdigits=5)) Wb") + println(" 1D profile points: nw = $nw") + println(" 2D grid: nR = $(length(r_grid)), nZ = $(length(z_grid))") + println(" R ∈ [$(round(rmin; sigdigits=4)), $(round(rmax; sigdigits=4))] m") + println(" Z ∈ [$(round(zmin; sigdigits=4)), $(round(zmax; sigdigits=4))] m") + + return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio) +end diff --git a/test/runtests.jl b/test/runtests.jl index 5393faeb..285ccfb5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,6 +16,7 @@ else include("./runtests_fastinterp.jl") include("./runtests_vacuum_julia.jl") include("./runtests_equil.jl") + include("./runtests_imas.jl") include("./runtests_solovev.jl") include("./runtests_eulerlagrange.jl") include("./runtests_sing.jl") diff --git a/test/runtests_imas.jl b/test/runtests_imas.jl new file mode 100644 index 00000000..25902e9c --- /dev/null +++ b/test/runtests_imas.jl @@ -0,0 +1,156 @@ +# Tests that the IMAS equilibrium loading path (read_imas, Task 1) produces +# the same result as reading the same equilibrium directly from a gEQDSK file. +# +# Strategy (as recommended in the task description): +# 1. Read the existing test gEQDSK with the standard EFIT path → reference equil +# 2. Convert that same gEQDSK to an IMAS dd using EFIT.jl's geqdsk2imas! +# 3. Load via the new IMAS path → test equil +# 4. Compare key quantities — they should match to within interpolation tolerance + +import EFIT +import EFIT.IMASdd + +@testset "IMAS equilibrium: read_imas matches read_efit" begin + + data_dir = joinpath(@__DIR__, "test_data", "regression_equilibrium_example") + geqdsk_path = joinpath(data_dir, "EQDSK_COCOS_02") + + # ── 1. Reference: load via the standard EFIT g-file path ───────────────── + efit_control = JPEC.Equilibrium.EquilibriumControl(; + eq_type = "efit", + eq_filename = geqdsk_path, + jac_type = "boozer", + psilow = 0.01, + psihigh = 0.994, + ) + efit_config = JPEC.Equilibrium.EquilibriumConfig(efit_control, JPEC.Equilibrium.EquilibriumOutput()) + equil_efit = JPEC.Equilibrium.setup_equilibrium(efit_config) + + # ── 2. Convert the same gEQDSK → IMAS dd using EFIT.jl ─────────────────── + # set_time=0.0 is required because EQDSK_COCOS_02 has a non-standard header + # that EFIT.jl cannot automatically parse a time value from. + g = EFIT.readg(geqdsk_path; set_time=0.0) + dd = IMASdd.dd() + + # Manually set up the time structure so time_slice[] resolves correctly + # in read_imas (which uses dd.global_time to select the active slice). + dd.global_time = g.time + dd.equilibrium.time = [g.time] + resize!(dd.equilibrium.time_slice, 1) + dd.equilibrium.time_slice[1].time = g.time + + # Fill all equilibrium profile data from the gEQDSK. + # wall=nothing: wall data is not needed for the equilibrium comparison. + # COCOS is auto-detected (geqdsk_cocos=0); dd is written in COCOS 11. + EFIT.geqdsk2imas!(g, dd.equilibrium.time_slice[1]; wall=nothing) + + # ── 3. Load via the new IMAS path with the same grid settings ───────────── + # eq_filename is unused for IMAS data loading; its directory is used for + # any optional diagnostic output files (disabled by default). + imas_control = JPEC.Equilibrium.EquilibriumControl(; + eq_type = "imas", + eq_filename = joinpath(data_dir, "imas_placeholder"), + jac_type = "boozer", + psilow = 0.01, + psihigh = 0.994, + ) + imas_config = JPEC.Equilibrium.EquilibriumConfig(imas_control, JPEC.Equilibrium.EquilibriumOutput()) + equil_imas = JPEC.Equilibrium.setup_equilibrium(imas_config, dd) + + @test equil_imas isa JPEC.Equilibrium.PlasmaEquilibrium + + # ── 4. One-to-one comparison ────────────────────────────────────────────── + # Both results come from the same physical data, so differences are only + # from COCOS conversion (EFIT.jl) and our IMAS reader interpolation. + # Tolerances are set at ~0.1% for scalar quantities and ~1% for profiles. + + # Magnetic axis R position + @test isapprox(equil_efit.ro, equil_imas.ro, rtol=1e-3) + + # Magnetic axis Z position (may be near zero; use absolute tolerance) + @test isapprox(equil_efit.zo, equil_imas.zo, atol=1e-3) + + # Normalising flux |psi_boundary - psi_axis| + @test isapprox(equil_efit.psio, equil_imas.psio, rtol=1e-3) + + # Safety-factor profile q(psi): compare grid-point values + @test isapprox(equil_efit.profiles.q_spline.y, + equil_imas.profiles.q_spline.y, rtol=1e-2) + + # Pressure profile p(psi) + @test isapprox(equil_efit.profiles.P_spline.y, + equil_imas.profiles.P_spline.y, rtol=1e-2) + +end + +@testset "IMAS write: write_imas populates dd.mhd_linear" begin + + # Build a minimal synthetic DCON result to test write_imas without running + # the full solver. We just verify that the IMAS data dictionary is filled + # correctly from the result NamedTuple. + + # ── 1. Create a dd with an equilibrium time slice already set ───────────── + dd_out = IMASdd.dd() + dd_out.equilibrium.time = [0.0] + resize!(dd_out.equilibrium.time_slice, 1) + dd_out.equilibrium.time_slice[1].time = 0.0 + + # ── 2. Build minimal control / internal structs ─────────────────────────── + # n = 1 run, 7 poloidal modes (mlow = -3 … mhigh = 3) + n_modes = 7 + ctrl = JPEC.DCON.DconControl(; + nn_low = 1, + nn_high = 1, + vac_flag = true, + verbose = false, + ) + intr = JPEC.DCON.DconInternal(; + nlow = 1, + nhigh = 1, + npert = 1, + mlow = -3, + mhigh = 3, + mpert = n_modes, + numpert_total = n_modes, + ) + + # ── 3. Synthetic VacuumData — eigenvalues in ascending real order ───────── + # (negative real part → unstable, positive → stable) + vac_data = JPEC.DCON.VacuumData(10, n_modes, n_modes) + vac_data.et .= ComplexF64[-0.3, 0.1, 0.2, 0.4, 0.6, 0.8, 1.2] # et[1] < 0: unstable + + fake_result = (ctrl=ctrl, equil=nothing, intr=intr, + ffit=nothing, odet=nothing, vac_data=vac_data) + + # ── 4. Run write_imas and check the dd ──────────────────────────────────── + JPEC.DCON.write_imas(dd_out, fake_result) + + mhd = dd_out.mhd_linear + + # Top-level metadata + @test mhd.ideal_flag == 1 + @test mhd.ids_properties.homogeneous_time == 1 + @test mhd.code.name == "JPEC" + + # One time slice with the right time + @test length(mhd.time_slice) == 1 + @test mhd.time_slice[1].time == 0.0 + + # One toroidal_mode entry per eigenmode + @test length(mhd.time_slice[1].toroidal_mode) == n_modes + + # All entries should have n_tor == 1 (single n = 1 run) + for tm in mhd.time_slice[1].toroidal_mode + @test tm.n_tor == 1 + end + + # First eigenmode: most unstable, energy_perturbed < 0 + @test mhd.time_slice[1].toroidal_mode[1].energy_perturbed < 0 + + # Energy values should match the real parts of et + for i in 1:n_modes + @test isapprox(mhd.time_slice[1].toroidal_mode[i].energy_perturbed, + real(vac_data.et[i]); atol=1e-10) + end + +end