From 47cc30bebaa4e89bc606c50b0d7bc01eb496e139 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Mon, 9 Mar 2026 00:53:28 -0700 Subject: [PATCH 1/3] Add MXH to IMAS.jl --- src/physics/fluxsurfaces.jl | 61 ++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/src/physics/fluxsurfaces.jl b/src/physics/fluxsurfaces.jl index 74cf4e04..115a4e9f 100644 --- a/src/physics/fluxsurfaces.jl +++ b/src/physics/fluxsurfaces.jl @@ -1650,7 +1650,7 @@ function flux_surfaces(eqt::equilibrium__time_slice{T1}, wall_r::AbstractVector{ # trace flux surfaces Br, Bz = Br_Bz(eqt2d) - surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z) + surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z;refine_extrema=false) # calculate flux surface averaged and geometric quantities N = length(eqt1d.psi) @@ -2342,3 +2342,62 @@ end @compat public areal_elongation push!(document[Symbol("Physics flux-surfaces")], :areal_elongation) + + + +""" + rec_to_mxh_coeffs(dd::IMAS.dd) + +fits rectangular flux surface mesh in dd to MXH +""" +function rec_to_mxh_coeffs(dd; MXH_modes=4, λ=10.0) + return rec_to_mxh_coeffs(dd.equilibrium.time_slice[], dd.wall; MXH_modes, λ) +end + +function rec_to_mxh_coeffs(eqt::IMAS.equilibrium__time_slice, wall::IMAS.wall; MXH_modes=4, λ=10.0) + + function extrap(x::Vector, u::Vector) + m = (u[3] - u[2]) / (x[3] - x[2]) + return u[3] - m * x[3] + end + fw = first_wall(wall) + + pnorm = eqt.profiles_1d.psi_norm + surfaces = trace_surfaces(eqt, fw.r, fw.z) + ns = length(surfaces) + + nimas = 5 + 2*MXH_modes + coeffs = zeros(nimas, ns) + mxh_fits = Vector{Any}(nothing, ns) + + for isurface in 1:ns + r,z = surfaces[isurface].r, surfaces[isurface].z + # Inside-out with extrapolated prior + mxh_prev = if isurface <3 + nothing + else + # Extrapolate from two already-fitted inner surfaces + x0, x1, x2 = pnorm[isurface], pnorm[isurface-1], pnorm[isurface-2] + mxh1, mxh2 = mxh_fits[isurface-1], mxh_fits[isurface-2] + c0_ext = mxh1.c0 + (mxh1.c0 - mxh2.c0) / (x1 - x2) * (x0 - x1) + c_ext = mxh1.c .+ (mxh1.c .- mxh2.c) ./ (x1 - x2) .* (x0 - x1) + s_ext = mxh1.s .+ (mxh1.s .- mxh2.s) ./ (x1 - x2) .* (x0 - x1) + MXH(mxh1.R0, mxh1.Z0, mxh1.ϵ, mxh1.κ, c0_ext, c_ext, s_ext) + end + mxh_fits[isurface] = MXH(r, z, MXH_modes; spline=false, optimize_fit=false, mxh_prev, λ) + coeffs[:, isurface] = flat_coeffs(mxh_fits[isurface]) + end + + rnorm = sqrt.(pnorm) + coeffs[1, 1] = extrap(pnorm, coeffs[1, :]) # R0: linear in pnorm + coeffs[2, 1] = extrap(rnorm, coeffs[2, :]) # Z0: linear in rnorm + coeffs[3, 1] = 0.0 # ϵ = 0 at axis + coeffs[4, 1] = extrap(rnorm, coeffs[4, :]) # κ: linear in rnorm + coeffs[5, 1] = extrap(pnorm, coeffs[5, :]) # c0: linear in pnorm + coeffs[6:end, 1] .= 0.0 # Fourier terms = 0 + + return coeffs +end + +@compat public rec_to_mxh_coeffs +push!(document[Symbol("Physics flux-surfaces")], :rec_to_mxh_coeffs) From 2e7a52451580f963dd8cf00d1e87b5c9d4771263 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Mon, 16 Mar 2026 16:07:10 -0700 Subject: [PATCH 2/3] Update fluxsurfaces.jl Revert refine_extrema back to defaulting to true --- src/physics/fluxsurfaces.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/fluxsurfaces.jl b/src/physics/fluxsurfaces.jl index 115a4e9f..53533264 100644 --- a/src/physics/fluxsurfaces.jl +++ b/src/physics/fluxsurfaces.jl @@ -1650,7 +1650,7 @@ function flux_surfaces(eqt::equilibrium__time_slice{T1}, wall_r::AbstractVector{ # trace flux surfaces Br, Bz = Br_Bz(eqt2d) - surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z;refine_extrema=false) + surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z;refine_extrema=true) # calculate flux surface averaged and geometric quantities N = length(eqt1d.psi) From 9e11fc70096a4fdff9f346e8bee675714fb305b6 Mon Sep 17 00:00:00 2001 From: jmcclena Date: Tue, 17 Mar 2026 22:50:40 -0700 Subject: [PATCH 3/3] Update interferometer.jl --- src/physics/interferometer.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/interferometer.jl b/src/physics/interferometer.jl index b37d52cc..751576c0 100644 --- a/src/physics/interferometer.jl +++ b/src/physics/interferometer.jl @@ -61,7 +61,7 @@ function line_average( @assert ϕ1 == ϕ2 "line_average cannot handle different ϕ yet" # Create 1D interpolant for the quantity q - q_interp = cubic_interp1d(rho_tor_norm, q) + q_interp = linear_interp1d(rho_tor_norm, q) # Find intersections between line and LCFS intersections = intersection([r1, r2], [z1, z2], lcfs_r, lcfs_z) @@ -83,7 +83,7 @@ function line_average( rho_seg = rho_interp.RHO_interpolant.(r_seg, z_seg) # Interpolate quantity q along the segment - q_seg = q_interp.(rho_seg) + q_seg = max.(q_interp.(rho_seg), zero(T)) # Calculate path length inside LCFS path_length_inside_lcfs = sqrt((r_exit - r_entry)^2 + (z_exit - z_entry)^2)