Calin/imas

Project.toml

@@ -1,39 +1,69 @@
 name = "JPEC"
 uuid = "462872dd-e066-4d2e-b993-6468b5239634"
 license = "MIT"
-authors = ["Matthew Pharr <matthewcpharr@gmail.com>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Jacob Halpern <jmh2363@columbia.edu>"]
 version = "0.1.0"
+authors = ["Matthew Pharr <matthewcpharr@gmail.com>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Jacob Halpern <jmh2363@columbia.edu>"]
 
 [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"]

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

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,14 +210,25 @@ 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")
     println("Normal termination.")
 
     # 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
 
 """

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

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