matt-pharr · Oct 26, 2026 · Oct 20, 2026 · Oct 20, 2026 · Oct 20, 2026 · Oct 20, 2026
diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
@@ -514,7 +514,7 @@ function setup_equil(p::FortranRunParams, fortran_dir::String)
         "flux_surface_rtol" => 1e-7
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, fortran_dir)
-    return Equilibrium.setup_equilibrium(eq_config)
+    return Equilibrium.PlasmaEquilibrium(eq_config)
 end
 
 # ─── Phi_x comparison (coil mode) ────────────────────────────────────────────

diff --git a/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl b/benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl
@@ -207,7 +207,7 @@ t_equil = @elapsed begin
         "flux_surface_rtol" => 1e-7,
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, run_dir)
-    equil = Equilibrium.setup_equilibrium(eq_config)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config)
 end
 @printf "    Done in %.1f s  (bt_sign=%d)\n" t_equil equil.params.bt_sign
 

diff --git a/benchmarks/benchmark_delta_prime_methods.jl b/benchmarks/benchmark_delta_prime_methods.jl
@@ -28,7 +28,7 @@ function setup_and_run_solovev()
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)

diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl
@@ -23,7 +23,7 @@ function reference_q()
     inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
     equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "core_edge_packed", "mpsi" => 1024))
     eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR)
-    equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    equil = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
     # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable.
     return equil, equil_dict
 end

diff --git a/benchmarks/benchmark_riccati_der.jl b/benchmarks/benchmark_riccati_der.jl
@@ -29,7 +29,7 @@ function setup_solovev()
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)

diff --git a/benchmarks/benchmark_threads.jl b/benchmarks/benchmark_threads.jl
@@ -17,7 +17,7 @@ function run_ffs(ex; integrator)
     ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
     eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(eq_config)
+    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)

diff --git a/benchmarks/equil_method_comparison.jl b/benchmarks/equil_method_comparison.jl
@@ -50,9 +50,9 @@ for method in methods
     local t2 = 0.0
     local success = true
     try
-        setup_equilibrium(cfg)  # JIT warmup
-        t1 = @elapsed pe = setup_equilibrium(cfg)
-        t2 = @elapsed setup_equilibrium(cfg)
+        PlasmaEquilibrium(cfg)  # JIT warmup
+        t1 = @elapsed pe = PlasmaEquilibrium(cfg)
+        t2 = @elapsed PlasmaEquilibrium(cfg)
     catch e
         @warn "  FAILED: $e"
         success = false

diff --git a/benchmarks/equil_numerical_params.jl b/benchmarks/equil_numerical_params.jl
@@ -56,7 +56,7 @@ end
 # JIT warmup
 println("JIT warmup...")
 warmup_cfg = make_config(config_path, "efit", 64, 128)
-setup_equilibrium(warmup_cfg)
+PlasmaEquilibrium(warmup_cfg)
 
 rows = []
 
@@ -72,7 +72,7 @@ for method in ["efit", "efit_arclength", "efit_by_inversion"]
         local runtime = NaN
         local rt_err = NaN
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             raw = Equilibrium.read_efit(cfg)
             rt_err = roundtrip_error(pe, raw)
@@ -98,7 +98,7 @@ for method in ["efit", "efit_arclength", "efit_by_inversion"]
         local runtime = NaN
         local rt_err = NaN
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             raw = Equilibrium.read_efit(cfg)
             rt_err = roundtrip_error(pe, raw)
@@ -123,7 +123,7 @@ println("=" ^ 65)
 
 # ── Contour grid refinement scan (efit_by_inversion only) ─────────────────────
 # The refine factor in equilibrium_solver_by_inversion is a keyword arg.
-# We call it directly here, bypassing the setup_equilibrium wrapper.
+# We call it directly here, bypassing the PlasmaEquilibrium constructor.
 println("\n" * "=" ^ 65)
 println("Contour grid refinement scan (efit_by_inversion only, mpsi=128, mtheta=256)")
 println("=" ^ 65)

diff --git a/benchmarks/equil_psihigh_scan.jl b/benchmarks/equil_psihigh_scan.jl
@@ -97,7 +97,7 @@ for method in methods
         local success = false
         local err_msg = ""
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             success = true
         catch e

diff --git a/benchmarks/equil_psilow_scan.jl b/benchmarks/equil_psilow_scan.jl
@@ -50,7 +50,7 @@ end
 # JIT warmup
 println("JIT warmup...")
 warmup_cfg = make_config(config_path, "efit", 1e-2)
-setup_equilibrium(warmup_cfg)
+PlasmaEquilibrium(warmup_cfg)
 
 println("=" ^ 70)
 println("psilow near-axis accuracy scan")
@@ -67,7 +67,7 @@ for method in methods
         local success = false
         local err_msg = ""
         try
-            t = @elapsed pe = setup_equilibrium(cfg)
+            t = @elapsed pe = PlasmaEquilibrium(cfg)
             runtime = t
             success = true
         catch e

diff --git a/benchmarks/equil_spline_comparison.jl b/benchmarks/equil_spline_comparison.jl
@@ -82,7 +82,7 @@ for m in all_methods
     println("--- Running: $m ---")
     cfg = make_config(config_path, m, psihigh_arg)
     try
-        pes[m] = setup_equilibrium(cfg)
+        pes[m] = PlasmaEquilibrium(cfg)
         println("  Done.\n")
     catch e
         println("  FAILED: $e\n")

diff --git a/benchmarks/plot_grid_knot_placement.jl b/benchmarks/plot_grid_knot_placement.jl
@@ -19,7 +19,7 @@ function build_grids()
     tau = eq_config.psi_accuracy
 
     # Pass 1 (coarse layout) and the two-pass refined grid, as the driver builds them
-    equil1 = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    equil1 = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
     pass1 = copy(equil1.profiles.xs)
     rationals = GPE.ForceFreeStates.rational_psi_nodes(equil1; nlow=1, nhigh=1)
     auto = GPE.Equilibrium.refined_psi_grid(equil1; tau, mandatory=rationals)

diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -26,7 +26,7 @@ Splines are provided by the external `FastInterpolations` package rather than by
 ### Core Physics Modules
 
 2. **Equilibrium** (`src/Equilibrium/`) - MHD equilibrium solvers
-   - Main entry point: `setup_equilibrium(path)` or `setup_equilibrium(config)`
+   - Main entry point: `PlasmaEquilibrium(file; eq_type)`, `PlasmaEquilibrium(analytic_config)`, `PlasmaEquilibrium(dd)`, or `PlasmaEquilibrium(config, source)`
    - Supports multiple equilibrium types:
      - `efit` - EFIT g-file format
      - `chease`, `chease2` - CHEASE equilibrium code formats
@@ -149,14 +149,14 @@ Example configuration files are provided in:
 - `examples/Solovev_ideal_example/gpec.toml`
 - `examples/DIIID-like_ideal_example/gpec.toml`
 
-**Note**: Legacy configuration files (`equil.toml`, `vac.in`) are deprecated.
+**Note**: Legacy configuration files (`equil.toml`, `vac.in`) are not read; `gpec.toml` is the only input deck.
 
 ## Data Flow
 
 The complete GPEC analysis pipeline:
 
 1. **Equilibrium Setup**:
-   - `setup_equilibrium(config)` reads configuration from `gpec.toml`
+   - `PlasmaEquilibrium(config, source)` forms the equilibrium from the `[Equilibrium]` section of `gpec.toml`
    - Parses equilibrium data (EFIT, CHEASE, or analytical)
    - Runs Grad-Shafranov solver (direct or inverse)
    - Computes global parameters: q-profile, pressure, current density, β

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -28,9 +28,16 @@ ffs  = solve(prob, Riccati())
 The four objects map one-to-one onto the pipeline stages:
 
   - [`PlasmaEquilibrium`](@ref) reads and processes the equilibrium — the `[Equilibrium]`
-    section. Analytic equilibria (`sol`, `lar`, `tj_analytic`) take their parameters from a
-    separate config object and are built with
-    `setup_equilibrium(eq_config, analytic_config)` instead.
+    section. Analytic equilibria take their parameters from their config object, with the
+    `eq_type` implied, and IMAS equilibria from a data dictionary:
+
+    ```julia
+    eq = PlasmaEquilibrium(Equilibrium.SolovevConfig(; elongation=1.5); mpsi=0)
+    eq = PlasmaEquilibrium(dd; imas_cocos=11)
+    ```
+
+    Every equilibrium remembers its source, so with `mpsi=0` `solve` re-forms it on the same
+    refined two-pass grid the deck uses; `Equilibrium.reform(eq, psi_nodes)` re-forms it on any grid.
   - `solve` runs the force-free-states stage — the `[ForceFreeStates]` section — and returns
     a `ForceFreeStatesResult`, exactly the object the TOML driver publishes.
   - [`RMPField`](@ref) describes the external field — the `[ForcingTerms]` section — without

diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md
@@ -98,9 +98,9 @@ built by a **two-pass measured-curvature refinement** driven by the single accur
 Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines
 the core, pedestal, and edge proportionally. The fixed `grid_type = "core_edge_packed"`
 (sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed
-layout) are still supported. Library users calling `setup_equilibrium` directly with
-`mpsi = 0` receive the coarse pass-1 grid; use `refined_psi_grid` and the
-`override_psi_nodes` keyword to apply the refinement manually.
+layout) are still supported. A `PlasmaEquilibrium` built directly with `mpsi = 0` keeps
+the coarse pass-1 grid until `solve` re-forms it; to refine it yourself, pass the nodes from
+`refined_psi_grid` to `Equilibrium.reform`.
 
 The packing on the DIII-D-like example (n=1) compared to fixed `core_edge_packed` grids — note the
 coarse spacing across the smooth mid-radius, the spacing dips at each rational surface,
@@ -126,11 +126,8 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Important types
 
-- `EquilibriumConfig` — top-level configuration container parsed from a
-	TOML file (outer constructor `EquilibriumConfig(path::String)` is
-	provided).
-- `EquilibriumControl` — low-level control parameters (grid, jacobian
-	type, tolerances, etc.).
+- `EquilibriumConfig` — configuration container, built from the `[Equilibrium]`
+	section of `gpec.toml` or from keywords.
 - `PlasmaEquilibrium` — the runtime structure containing spline fields,
 	geometry, profiles, and computed diagnostics (q-profile, separatrix,
 	etc.).
@@ -140,9 +137,11 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Key functions
 
-- `setup_equilibrium(path::String = "equil.toml")`
-	— main entry point that reads configuration, builds the equilibrium,
-	runs the solver, and returns a `PlasmaEquilibrium` instance.
+- `PlasmaEquilibrium(file; eq_type="efit", kwargs...)`, `PlasmaEquilibrium(analytic_config; kwargs...)`
+	and `PlasmaEquilibrium(dd; kwargs...)` — read or solve the equilibrium and return a
+	`PlasmaEquilibrium`; `kwargs` are `EquilibriumConfig` fields.
+- `PlasmaEquilibrium(config::EquilibriumConfig, source=nothing)` — the form every other
+	constructor ends in, with `source` an analytic config, an IMAS `dd`, or `nothing` to read `eq_filename`.
 - `equilibrium_separatrix_find!(pe::PlasmaEquilibrium)` — locate
 	separatrix and related boundary geometry in-place.
 - `equilibrium_global_parameters!(pe::PlasmaEquilibrium)` — populate
@@ -154,15 +153,14 @@ Modules = [GeneralizedPerturbedEquilibrium.Equilibrium]
 
 ## Example usage
 
-Basic example: read a TOML config and build an equilibrium
+Basic example: read a g-file and build an equilibrium
 
 ```julia
 using GeneralizedPerturbedEquilibrium
 
-# Build from a TOML file (searches relative paths if needed)
-pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium("docs/examples/ForceFreeStates.toml")
+pe = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada")
 
-println("Magnetic axis: ", pe.params.r0, ", ", pe.params.z0)
+println("Magnetic axis: ", pe.ro, ", ", pe.zo)
 println("q(0) = ", pe.params.q0)
 
 # Find separatrix (in-place) and inspect results
@@ -175,9 +173,8 @@ Analytic / testing example: construct a large-aspect-ratio model
 ```julia
 using GeneralizedPerturbedEquilibrium
 
-# Create a LAR config from a small TOML fragment or file
 larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(r0=10.0, a=1.0, beta0=1e-3)
-pe = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(control=Dict("eq_filename"=>"unused","eq_type"=>"lar")), larcfg)
+pe = PlasmaEquilibrium(larcfg)  # eq_type = "lar" is implied by the config type
 
 println("Built LAR equilibrium with a = ", larcfg.a)
 ```

diff --git a/docs/src/set_up.md b/docs/src/set_up.md
@@ -341,7 +341,7 @@ using GeneralizedPerturbedEquilibrium.Vacuum
 using GeneralizedPerturbedEquilibrium.ForceFreeStates
 
 # Set up equilibrium only
-equil = Equilibrium.setup_equilibrium("path/to/gpec.toml")
+equil = Equilibrium.PlasmaEquilibrium("path/to/g-file"; eq_type="efit")
 
 # Access equilibrium data
 println("q at axis: ", equil.params.q0)

diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -414,7 +414,7 @@ inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
 
 ctrl  = FFS.ForceFreeStatesControl(;
             (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-equil = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(
+equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(
             GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex))
 
 intr  = FFS.ForceFreeStatesInternal(; dir_path=ex)

diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl
@@ -58,7 +58,7 @@ for geqdsk_file in geqdsk_cases
     eq_dict = deepcopy(base_toml["Equilibrium"])
     eq_dict["eq_filename"] = geqdsk_file
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, @__DIR__)
-    equil = Equilibrium.setup_equilibrium(eq_config)
+    equil = Equilibrium.PlasmaEquilibrium(eq_config)
     psi_norm = Vector(equil.profiles.xs)
 
     # s and alpha profiles vs normalized flux

diff --git a/examples/LAR_beta_scan/run_scan.jl b/examples/LAR_beta_scan/run_scan.jl
@@ -16,7 +16,7 @@ using Pkg
 Pkg.activate(joinpath(@__DIR__, "../.."))
 
 using GeneralizedPerturbedEquilibrium
-using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig, setup_equilibrium
+using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig
 using HDF5
 using TOML
 using Printf

diff --git a/examples/LAR_epsilon_scan/run_scan.jl b/examples/LAR_epsilon_scan/run_scan.jl
@@ -16,7 +16,7 @@ using Pkg
 Pkg.activate(joinpath(@__DIR__, "../.."))
 
 using GeneralizedPerturbedEquilibrium
-using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig, setup_equilibrium
+using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig
 using HDF5
 using TOML
 using Printf

diff --git a/examples/efit_fixedbdy_separatrix_example/README.md b/examples/efit_fixedbdy_separatrix_example/README.md
@@ -13,7 +13,7 @@ loads it regardless of box size.
 
 This g-file is the fixture for the `efit_fixedbdy_separatrix` regression case
 (`regression-harness/cases/efit_fixedbdy_separatrix.toml`, a `kind = "computed"` case that calls
-`setup_equilibrium` directly). It is **not** a runnable example deck — there is no `gpec.toml`.
+`PlasmaEquilibrium` directly). It is **not** a runnable example deck — there is no `gpec.toml`.
 
 Free-boundary equilibria (real coil currents) and standard diverted EFITs are unaffected: their
 outboard midplane flux keeps dropping past the LCFS, so the bracket was never ambiguous.
diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl
@@ -151,7 +151,7 @@ end
 # computational box hugs the LCFS (eps=0.05 TokaMaker aspect-scan g-file): outside the prescribed
 # LCFS the coil-vacuum flux turns back above the boundary value before the grid edge, so the old
 # bracketed-Brent separatrix finder could not bracket psi=sibry and equilibrium setup threw. Calls
-# setup_equilibrium directly and writes the leading equilibrium scalars: errors on the buggy code,
+# PlasmaEquilibrium directly and writes the leading equilibrium scalars: errors on the buggy code,
 # passes with the Newton finder. Runtime and environment are written to ARGS[2] by the run-info epilogue.
 const COMPUTED_SEPARATRIX_SCRIPT_TEMPLATE = """
 using Pkg
@@ -161,9 +161,8 @@ using GeneralizedPerturbedEquilibrium.Equilibrium
 using HDF5
 root = dirname(Base.active_project())
 gfile = joinpath(root, "examples", "efit_fixedbdy_separatrix_example", "eq_eps0.0500000_k1.000_d0.000.geqdsk")
-cfg = Equilibrium.EquilibriumConfig(; eq_type="efit", eq_filename=gfile)
 t_start = time()
-pe = Equilibrium.setup_equilibrium(cfg)
+pe = Equilibrium.PlasmaEquilibrium(gfile; eq_type="efit")
 elapsed = time() - t_start
 h5open(ARGS[1], "w") do fid
     fid["Equilibrium/psi_total"]  = pe.psio

diff --git a/src/Equilibrium/DirectEquilibriumByInversion.jl b/src/Equilibrium/DirectEquilibriumByInversion.jl
@@ -703,7 +703,7 @@ function equilibrium_solver_by_inversion(
     rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_table; itp_opts2d...)
 
     # Intermediate inverse input; the captured DirectIngest rides on the original eq_input,
-    # which setup_equilibrium forwards onto the equilibrium, so this one carries ingest=nothing.
+    # which PlasmaEquilibrium forwards onto the equilibrium, so this one carries ingest=nothing.
     inv_input = InverseRunInput(raw_profile.config, raw_profile.sq_in,
         rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing, raw_profile.psihigh_resolved)