Tearing - FEATURE - Add critical resonant field calculations to the tearing workflow

.gitignore

@@ -38,3 +38,4 @@ scratch/
 
 # Local profiling scratch (one-off study scripts, not part of the package)
 profiling/
+src/Tearing/CriticalResonantField/CRF Dev/

Project.toml

@@ -1,7 +1,7 @@
 name = "GeneralizedPerturbedEquilibrium"
 uuid = "462872dd-e066-4d2e-b993-6468b5239634"
 license = "MIT"
-authors = ["Nikolas Logan <ncl2128@columbia.edu>", "Jong-Kyu Park <jkpark@snu.ac.kr>", "Matthew Pharr <m.pharr@protonmail.com>", "Jacob Halpern <jmh2363@columbia.edu>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Daniel Burgess <dab2245@columbia.edu>", "Min-Gu Yoo <yoom@fusion.gat.com>"]
+authors = ["Nikolas Logan <ncl2128@columbia.edu>", "Jong-Kyu Park <jkpark@snu.ac.kr>", "Matthew Pharr <m.pharr@protonmail.com>", "Jacob Halpern <jmh2363@columbia.edu>", "Rithik Banerjee <rb3736@columbia.edu>", "Jaebeom Cho <aspire1019@snu.ac.kr>", "Daniel Burgess <dab2245@columbia.edu>", "Min-Gu Yoo <yoom@fusion.gat.com>","Evan Bursch <emb2333@columbia.edu>"]
 version = "0.1.0"
 
 [deps]
@@ -69,5 +69,3 @@ Statistics = "1"
 TOML = "1"
 Test = "1"
 julia = "1.11"
-
-

benchmarks/benchmark_against_fortran_run.jl

@@ -109,7 +109,7 @@ end
 
 function _find_indexed2_float(text::String, name::String, i::Int, j::Int; default=nothing)
     m = match(Regex(
-            "\\b$(name)\\s*\\(\\s*$(i)\\s*,\\s*$(j)\\s*\\)\\s*=\\s*([+-]?[\\d.]+(?:[eE][+-]?\\d+)?)", "i"),
+        "\\b$(name)\\s*\\(\\s*$(i)\\s*,\\s*$(j)\\s*\\)\\s*=\\s*([+-]?[\\d.]+(?:[eE][+-]?\\d+)?)", "i"),
         text)
     isnothing(m) ? default : parse(Float64, m.captures[1])
 end
@@ -548,17 +548,17 @@ function compare_phix(
 
     grid = ForcingTerms.sample_boundary_grid(equil, p.mtheta_coil, p.nzeta_coil; psi=psilim)
     nobs = p.mtheta_coil * p.nzeta_coil
-    obs_R = zeros(nobs);
-    obs_phi = zeros(nobs);
+    obs_R = zeros(nobs)
+    obs_phi = zeros(nobs)
     obs_Z = zeros(nobs)
     for j in 1:p.nzeta_coil, i in 1:p.mtheta_coil
         idx = i + (j - 1) * p.mtheta_coil
         obs_R[idx] = grid.R[i]
         obs_phi[idx] = grid.phi_grid[j] + grid.phi_offset[i]
         obs_Z[idx] = grid.Z[i]
     end
-    B_R = zeros(nobs);
-    B_phi = zeros(nobs);
+    B_R = zeros(nobs)
+    B_phi = zeros(nobs)
     B_Z = zeros(nobs)
     ForcingTerms.compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, coil_sets)
     bn = zeros(p.mtheta_coil, p.nzeta_coil)
@@ -639,9 +639,9 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
     push!(lines, "")
 
     push!(lines, "--- Mercier Criterion ---")
-    fdi = fort["di"];
+    fdi = fort["di"]
     fdr = fort["dr"]
-    jdi = julia["di"];
+    jdi = julia["di"]
     jdr = julia["dr"]
     if !isempty(jdi) && !isempty(fdi)
         jdi_on_fg = [_interp1(jq_psi, jdi, p) for p in fq_psi]
@@ -857,7 +857,7 @@ function build_comparison_table(fort, julia, fortran_dir, bench_dir, nn)
 end
 
 function _interp1(x::AbstractVector, y::AbstractVector, xi::Real)
-    n = length(x);
+    n = length(x)
     n < 2 && return y[1]
     xi <= x[1] && return y[1]
     xi >= x[end] && return y[end]
@@ -1451,15 +1451,15 @@ function main(argv=ARGS)
     println()
     table_lines = build_comparison_table(fort, julia, fortran_dir, bench_dir, p.nn)
     for line in table_lines
-        ;
-        println(line);
+
+        println(line)
     end
 
     table_path = joinpath(bench_dir, "comparison_table.txt")
     open(table_path, "w") do io
         for line in table_lines
-            ;
-            println(io, line);
+
+            println(io, line)
         end
     end
     println("\nComparison table saved to: ", abspath(table_path))

benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl

@@ -25,7 +25,7 @@ using NCDatasets
 using Plots
 using Printf
 
-const OUTPUT_DIR  = joinpath(@__DIR__, "coil_pipeline")
+const OUTPUT_DIR = joinpath(@__DIR__, "coil_pipeline")
 const COIL_DAT_DIR = joinpath(@__DIR__, "..", "src", "ForcingTerms", "coil_geometries")
 
 # ---------------------------------------------------------------------------
@@ -89,38 +89,38 @@ end
 
 function parse_fortran_run(dir::String)::FortranRunParams
     equil_file = joinpath(dir, "equil.in")
-    dcon_file  = joinpath(dir, "dcon.in")
-    coil_file  = joinpath(dir, "coil.in")
+    dcon_file = joinpath(dir, "dcon.in")
+    coil_file = joinpath(dir, "coil.in")
 
     isfile(equil_file) || error("equil.in not found in $dir")
-    isfile(dcon_file)  || error("dcon.in not found in $dir")
-    isfile(coil_file)  || error("coil.in not found in $dir")
+    isfile(dcon_file) || error("dcon.in not found in $dir")
+    isfile(coil_file) || error("coil.in not found in $dir")
 
-    eq_text   = _strip_fortran_comments(read(equil_file, String))
-    dcon_text = _strip_fortran_comments(read(dcon_file,  String))
-    coil_text = _strip_fortran_comments(read(coil_file,  String))
+    eq_text = _strip_fortran_comments(read(equil_file, String))
+    dcon_text = _strip_fortran_comments(read(dcon_file, String))
+    coil_text = _strip_fortran_comments(read(coil_file, String))
 
     # Equilibrium parameters
-    eq_type   = _find_string(eq_text, "eq_type";   default="efit")
-    eq_file   = _find_string(eq_text, "eq_filename"; default="")
-    jac_type  = _find_string(eq_text, "jac_type";  default="hamada")
-    psilow    = _find_scalar(eq_text, "psilow";    default=1e-4)
-    psihigh   = _find_scalar(eq_text, "psihigh";   default=0.99)
-    mtheta_eq = _find_int(eq_text,    "mtheta";    default=256)
+    eq_type = _find_string(eq_text, "eq_type"; default="efit")
+    eq_file = _find_string(eq_text, "eq_filename"; default="")
+    jac_type = _find_string(eq_text, "jac_type"; default="hamada")
+    psilow = _find_scalar(eq_text, "psilow"; default=1e-4)
+    psihigh = _find_scalar(eq_text, "psihigh"; default=0.99)
+    mtheta_eq = _find_int(eq_text, "mtheta"; default=256)
     grid_type = _find_string(eq_text, "grid_type"; default="ldp")
 
     # Toroidal mode number and m-range expansion
-    nn          = _find_int(dcon_text, "nn";          default=1)
-    delta_mlow  = _find_int(dcon_text, "delta_mlow";  default=8)
+    nn = _find_int(dcon_text, "nn"; default=1)
+    delta_mlow = _find_int(dcon_text, "delta_mlow"; default=8)
     delta_mhigh = _find_int(dcon_text, "delta_mhigh"; default=8)
 
     # Coil grid parameters
     mtheta_coil = _find_int(coil_text, "cmtheta"; default=480)
-    nzeta_coil  = _find_int(coil_text, "cmzeta";  default=40)
-    machine     = _find_string(coil_text, "machine"; default="d3d")
-    coil_num    = _find_int(coil_text, "coil_num"; default=1)
+    nzeta_coil = _find_int(coil_text, "cmzeta"; default=40)
+    machine = _find_string(coil_text, "machine"; default="d3d")
+    coil_num = _find_int(coil_text, "coil_num"; default=1)
 
-    coil_names    = String[]
+    coil_names = String[]
     coil_currents = Vector{Float64}[]
     for ci in 1:coil_num
         name = _find_indexed_string(coil_text, "coil_name", ci; default="coil$ci")
@@ -135,9 +135,9 @@ function parse_fortran_run(dir::String)::FortranRunParams
     end
 
     return FortranRunParams(eq_type, eq_file, jac_type, psilow, psihigh, mtheta_eq,
-                            grid_type, nn, delta_mlow, delta_mhigh,
-                            mtheta_coil, nzeta_coil, machine,
-                            coil_names, coil_currents)
+        grid_type, nn, delta_mlow, delta_mhigh,
+        mtheta_coil, nzeta_coil, machine,
+        coil_names, coil_currents)
 end
 
 function find_output_nc(dir::String, nn::Int)
@@ -177,16 +177,16 @@ end
 nc_path = find_output_nc(run_dir, p.nn)
 
 # Read Fortran output (psilim and Phi_x spectrum)
-fortran_m      = Int[]
-fortran_Phix   = Float64[]
+fortran_m = Int[]
+fortran_Phix = Float64[]
 fortran_psilim = NaN
 NCDatasets.Dataset(nc_path, "r") do ds
-    global fortran_m      = Int.(ds["m"][:])
+    global fortran_m = Int.(ds["m"][:])
     global fortran_psilim = Float64(ds.attrib["psilim"])
     phi_x_raw = ds["Phi_x"][:, :]   # (mpert, 2); col 1=cos, col 2=sin
-    global fortran_Phix = sqrt.(phi_x_raw[:, 1].^2 .+ phi_x_raw[:, 2].^2)
+    global fortran_Phix = sqrt.(phi_x_raw[:, 1] .^ 2 .+ phi_x_raw[:, 2] .^ 2)
 end
-mlow  = minimum(fortran_m)
+mlow = minimum(fortran_m)
 mhigh = maximum(fortran_m)
 @printf "    Fortran psilim=%.6f,  m=%d:%d\n" fortran_psilim mlow mhigh
 
@@ -196,15 +196,15 @@ mhigh = maximum(fortran_m)
 println("\n[2/4] Loading equilibrium ...")
 t_equil = @elapsed begin
     eq_dict = Dict{String,Any}(
-        "eq_type"      => p.eq_type,
-        "eq_filename"  => p.eq_filename,
-        "jac_type"     => p.jac_type,
-        "psilow"       => p.psilow,
-        "psihigh"      => p.psihigh,
-        "mtheta"       => p.mtheta_equil,
-        "grid_type"    => p.grid_type,
+        "eq_type" => p.eq_type,
+        "eq_filename" => p.eq_filename,
+        "jac_type" => p.jac_type,
+        "psilow" => p.psilow,
+        "psihigh" => p.psihigh,
+        "mtheta" => p.mtheta_equil,
+        "grid_type" => p.grid_type,
         "psi_accuracy" => 0.001,
-        "etol"         => 1e-7,
+        "etol" => 1e-7
     )
     eq_config = Equilibrium.EquilibriumConfig(eq_dict, run_dir)
     equil = Equilibrium.setup_equilibrium(eq_config)
@@ -223,34 +223,38 @@ for (name, currents) in zip(p.coil_names, p.coil_currents)
     isfile(dat_path) || error("Coil geometry file not found: $dat_path")
     raw = read_coil_dat(dat_path)
     push!(coil_sets, CoilSet(raw.name, raw.ncoil, raw.s, raw.nw, raw.nsec,
-                             raw.x, raw.y, raw.z, Float64.(currents[1:raw.ncoil])))
+        raw.x, raw.y, raw.z, Float64.(currents[1:raw.ncoil])))
 end
 for cs in coil_sets
     @printf "    %-10s %d conductors  currents: %s A\n" cs.name cs.ncoil string(round.(cs.currents; digits=1))
 end
 
 julia_modes = ForcingMode[]
 bnd_grid = nothing
-bn       = nothing
-t_biot   = @elapsed begin
+bn = nothing
+t_biot = @elapsed begin
     bnd_grid = sample_boundary_grid(equil, p.mtheta_coil, p.nzeta_coil; psi=fortran_psilim)
     nobs = p.mtheta_coil * p.nzeta_coil
-    obs_R   = zeros(nobs); obs_phi = zeros(nobs); obs_Z = zeros(nobs)
+    obs_R = zeros(nobs)
+    obs_phi = zeros(nobs)
+    obs_Z = zeros(nobs)
     for j in 1:p.nzeta_coil, i in 1:p.mtheta_coil
         idx = i + (j - 1) * p.mtheta_coil
-        obs_R[idx]   = bnd_grid.R[i]
+        obs_R[idx] = bnd_grid.R[i]
         obs_phi[idx] = bnd_grid.phi_grid[j] + bnd_grid.phi_offset[i]
-        obs_Z[idx]   = bnd_grid.Z[i]
+        obs_Z[idx] = bnd_grid.Z[i]
     end
-    B_R = zeros(nobs); B_phi = zeros(nobs); B_Z = zeros(nobs)
+    B_R = zeros(nobs)
+    B_phi = zeros(nobs)
+    B_Z = zeros(nobs)
     compute_biot_savart_boundary!(B_R, B_phi, B_Z, obs_R, obs_phi, obs_Z, coil_sets)
     bn = zeros(p.mtheta_coil, p.nzeta_coil)
     project_normal_flux!(bn, B_R, B_Z, bnd_grid)
     append!(julia_modes, fourier_decompose_bn(bn, bnd_grid, p.nn, mlow, mhigh))
 end
 @printf "    Done in %.1f s,  max|Phi_x| = %.3e T·m²\n" t_biot maximum(abs, bn)
 
-julia_m    = [md.m for md in sort(julia_modes; by=md -> md.m)]
+julia_m = [md.m for md in sort(julia_modes; by=md -> md.m)]
 julia_amps = abs.([md.amplitude for md in sort(julia_modes; by=md -> md.m)])
 
 # Print comparison table
@@ -259,7 +263,7 @@ println("    " * "-"^55)
 for (jm, ja) in zip(julia_m, julia_amps)
     k = findfirst(==(jm), fortran_m)
     if !isnothing(k)
-        fp    = fortran_Phix[k]
+        fp = fortran_Phix[k]
         ratio = ja > 1e-40 ? fp / ja : NaN
         @printf "    %3d    %12.4e         %12.4e      %.4f\n" jm ja fp ratio
     end
@@ -272,7 +276,7 @@ println("    Expected ratio ≈ 1.000 for all modes.")
 println("\n[4/4] Generating figure ...")
 
 function step_series(m_vals, amps)
-    m_ext   = [m_vals[1] - 1; m_vals; m_vals[end] + 1]
+    m_ext = [m_vals[1] - 1; m_vals; m_vals[end] + 1]
     amp_ext = [0.0; amps; 0.0]
     return m_ext, amp_ext
 end
@@ -286,19 +290,20 @@ p2 = plot_bn_contour(bn, p.mtheta_coil, p.nzeta_coil; n=p.nn)
 title!(p2, "Normal flux Φₓ(θ, ζ)  at ψ=$(round(fortran_psilim; digits=4))")
 
 # Panel 3: Julia vs Fortran Phi_x spectrum
-jm_ext, ja_ext   = step_series(julia_m, julia_amps)
-fm_ext, fa_ext   = step_series(fortran_m, fortran_Phix)
+jm_ext, ja_ext = step_series(julia_m, julia_amps)
+fm_ext, fa_ext = step_series(fortran_m, fortran_Phix)
 p3 = plot(; xlabel="Poloidal mode m", ylabel="|Phi_x| [T·m²]",
-            title="Spectrum: Julia vs Fortran  (n=$(p.nn))",
-            legend=:topright)
+    title="Spectrum: Julia vs Fortran  (n=$(p.nn))",
+    legend=:topright)
 plot!(p3, jm_ext, ja_ext; seriestype=:steppre, lw=2, color=:blue,
-      label="Julia ψ=$(round(fortran_psilim; digits=4))")
+    label="Julia ψ=$(round(fortran_psilim; digits=4))")
 plot!(p3, fm_ext, fa_ext; seriestype=:steppre, lw=2, color=:orange,
-      linestyle=:dash, label="Fortran Phi_x")
-ylims!(p3, (0, Inf)); xlims!(p3, mlow - 2, mhigh + 2)
+    linestyle=:dash, label="Fortran Phi_x")
+ylims!(p3, (0, Inf));
+xlims!(p3, mlow - 2, mhigh + 2)
 
 # Panel 4: Phi_x / Julia ratio
-ratio_m    = Int[]
+ratio_m = Int[]
 ratio_vals = Float64[]
 for (jm, ja) in zip(julia_m, julia_amps)
     k = findfirst(==(jm), fortran_m)
@@ -308,8 +313,8 @@ for (jm, ja) in zip(julia_m, julia_amps)
     end
 end
 p4 = plot(; xlabel="Poloidal mode m", ylabel="Phi_x / Julia",
-            title="Ratio (flat at 1.0 → correct)",
-            legend=:topright)
+    title="Ratio (flat at 1.0 → correct)",
+    legend=:topright)
 if !isempty(ratio_m)
     scatter!(p4, ratio_m, ratio_vals; color=:blue, label="ratio", markersize=5)
     hline!(p4, [1.0]; color=:red, linestyle=:dash, label="expected = 1.0")
@@ -318,8 +323,8 @@ xlims!(p4, mlow - 2, mhigh + 2)
 
 run_label = basename(run_dir)
 fig = plot(p1, p2, p3, p4;
-           layout=(2, 2), size=(1400, 900),
-           plot_title="$run_label  |  $(p.machine)  n=$(p.nn)  psilim=$(round(fortran_psilim; digits=4))")
+    layout=(2, 2), size=(1400, 900),
+    plot_title="$run_label  |  $(p.machine)  n=$(p.nn)  psilim=$(round(fortran_psilim; digits=4))")
 
 out_png = joinpath(OUTPUT_DIR, "coil_pipeline_$(run_label)_n$(p.nn).png")
 savefig(fig, out_png)

benchmarks/benchmark_delta_prime_methods.jl

@@ -32,8 +32,8 @@ function setup_and_run_solovev()
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)
-    intr.nlow = ctrl.nn_low;
-    intr.nhigh = ctrl.nn_high;
+    intr.nlow = ctrl.nn_low
+    intr.nhigh = ctrl.nn_high
     intr.npert = 1
     FFS.sing_find!(intr, equil)
     intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow

benchmarks/benchmark_diiid_ideal_ntv_torque.jl

@@ -35,7 +35,9 @@ const KF = GPE.KineticForces
 const Eq = GPE.Equilibrium
 const PE = GPE.PerturbedEquilibrium
 
-"Fortran GPEC DIII-D ideal example run directory, from the environment (no on-disk default)."
+"""
+Fortran GPEC DIII-D ideal example run directory, from the environment (no on-disk default).
+"""
 default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID, or pass the Fortran run directory as the first argument"),
     ENV, "GPEC_FORTRAN_DIIID")
 

benchmarks/benchmark_diiid_kinetic_stability.jl

@@ -36,7 +36,9 @@ using GeneralizedPerturbedEquilibrium
 const GPE = GeneralizedPerturbedEquilibrium
 const AnalysisFFS = GPE.Analysis.ForceFreeStates
 
-"Fortran GPEC DIII-D kinetic example run directory, from the environment (no on-disk default)."
+"""
+Fortran GPEC DIII-D kinetic example run directory, from the environment (no on-disk default).
+"""
 default_fortran_dir() = get(() -> error("Set GPEC_FORTRAN_DIIID_DCON, or pass the Fortran run directory as the first argument"),
     ENV, "GPEC_FORTRAN_DIIID_DCON")
 

benchmarks/benchmark_riccati_der.jl

@@ -33,8 +33,8 @@ function setup_solovev()
     intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
         (Symbol(k) => v for (k, v) in inputs["Wall"])...)
     FFS.sing_lim!(intr, ctrl, equil)
-    intr.nlow = ctrl.nn_low;
-    intr.nhigh = ctrl.nn_high;
+    intr.nlow = ctrl.nn_low
+    intr.nhigh = ctrl.nn_high
     intr.npert = 1
     FFS.sing_find!(intr, equil)
     intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow