From 3e2e4f49202610ecb602991f774054a9b9be9589 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Mon, 2 Feb 2026 20:01:55 -0500 Subject: [PATCH 1/2] DCON - NEW - added complex smoothing at singularities --- src/DCON/DconStructs.jl | 2 +- src/DCON/Ode.jl | 4 ++-- src/DCON/Sing.jl | 14 +++++++++++++- src/Equilibrium/DirectEquilibrium.jl | 2 +- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/DCON/DconStructs.jl b/src/DCON/DconStructs.jl index 95378c1f..460d85e9 100644 --- a/src/DCON/DconStructs.jl +++ b/src/DCON/DconStructs.jl @@ -397,7 +397,7 @@ and a small set of temporary matrices and factors used to compute singular-layer gmat::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total^2) tmp::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, numpert_total, numpert_total) Afact::Union{Cholesky{ComplexF64,Matrix{ComplexF64}},Nothing} = nothing - singfac_vec::Vector{Float64} = Vector{Float64}(undef, numpert_total) + singfac_vec::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) end # Initialize function for OdeState with relevant parameters for array initialization diff --git a/src/DCON/Ode.jl b/src/DCON/Ode.jl index 32ccc100..81a02d49 100644 --- a/src/DCON/Ode.jl +++ b/src/DCON/Ode.jl @@ -152,7 +152,7 @@ function ode_axis_init!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.Pl end end # Determine psimax and classify next integration limit type - if odet.ising > intr.msing || intr.psilim < intr.sing[odet.ising].psifac || ctrl.singfac_min == 0 + if odet.ising > intr.msing || intr.psilim < intr.sing[odet.ising].psifac || ctrl.singfac_min == 0 || true odet.psimax = intr.psilim * (1 - eps) odet.next = "finish" else @@ -294,7 +294,7 @@ function ode_step!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaE # Callback to be run at every step, handles fixups, tolerances, and data storage cb = DiscreteCallback((u, t, integrator) -> true, integrator_callback!) - + @info "step!: Integrating from ψ = $((@sprintf "%.3f" odet.psifac)) to ψ = $((@sprintf "%.3f" odet.psimax))" # Advance differential equation to next singular surface or edge rtol = compute_tols(ctrl, intr, odet) # initial tolerances prob = ODEProblem(sing_der!, odet.u, (odet.psifac, odet.psimax), (ctrl, equil, ffit, intr, odet)) diff --git a/src/DCON/Sing.jl b/src/DCON/Sing.jl index f55ecb42..65679e7f 100644 --- a/src/DCON/Sing.jl +++ b/src/DCON/Sing.jl @@ -705,7 +705,19 @@ function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, # Compute singfac = 1 / (m - nq) odet.q = Spl.spline_eval!(equil.sq, psieval)[4] - odet.singfac_vec .= vec(1.0 ./ ((intr.mlow:intr.mhigh) .- odet.q .* (intr.nlow:intr.nhigh)')) + singfac_den = ComplexF64.((intr.mlow:intr.mhigh) .- odet.q .* (intr.nlow:intr.nhigh)') + # ϵ_res = 1e-10 + ϵ_res = 1e-2 + # Add 1j * eps * norm to denominators near resonance to avoid singularity + for (j, n) in enumerate(intr.nlow:intr.nhigh) + mres = clamp(round(Int, odet.q * n), intr.mlow, intr.mhigh) + i = mres - intr.mlow + 1 + x = abs(singfac_den[i, j]) + gauss_norm = exp(- (x / (5 * ϵ_res))^2) + singfac_den[i, j] += 1im * ϵ_res * gauss_norm + end + + odet.singfac_vec .= vec(1.0 ./ singfac_den) # kinetic stuff - skip for now if false #(TODO: kin_flag) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 20c7e794..bbb280a6 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -296,7 +296,7 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro:: prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params) sol = solve(prob, BS5(); callback=callback, reltol=1e-6, abstol=1e-8, dt=2π / 200, adaptive=true) - if sol.retcode != :Success && sol.retcode != :Terminated + if sol.retcode != SciMLBase.ReturnCode.Success && sol.retcode != SciMLBase.ReturnCode.Terminated error("ODE integration failed for psi = $psifac with code: $(sol.retcode)") end From e201a1646566f11991612fe0d69ff082a3761d80 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Tue, 3 Feb 2026 22:12:06 -0500 Subject: [PATCH 2/2] DCON - MINOR - renamed smoothing factor to eps_res and added to namelist --- examples/Solovev_ideal_example/dcon.toml | 2 ++ src/DCON/DconStructs.jl | 2 ++ src/DCON/Sing.jl | 8 ++++---- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/examples/Solovev_ideal_example/dcon.toml b/examples/Solovev_ideal_example/dcon.toml index dde9cd83..6ad80cf8 100644 --- a/examples/Solovev_ideal_example/dcon.toml +++ b/examples/Solovev_ideal_example/dcon.toml @@ -38,6 +38,8 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Maximum fraction of solutions allowed before re-normalized force_wv_symmetry = true # Forces vacuum energy matrix symmetry +eps_res = 0.001 + [WALL] shape = "conformal" # String selecting wall shape ["nowall", "conformal", "elliptical", "dee", "mod_dee", "from_file"] a = 0.2415 # The distance of the wall from the plasma in units of major radius (conformal), or minor radius parameter (others). diff --git a/src/DCON/DconStructs.jl b/src/DCON/DconStructs.jl index 460d85e9..06ad156a 100644 --- a/src/DCON/DconStructs.jl +++ b/src/DCON/DconStructs.jl @@ -170,6 +170,7 @@ A mutable struct containing control parameters for stability analysis, set by th - `write_outputs_to_HDF5::Bool` - Write results to HDF5 format - `HDF5_filename::String` - Name of HDF5 output file - `force_wv_symmetry::Bool` - Boolean flag to enforce symmetry in the vacuum response matrix + - `eps_res::Float64` - Small parameter for numerical stability in resonance calculations """ @kwdef mutable struct DconControl verbose::Bool = true @@ -223,6 +224,7 @@ A mutable struct containing control parameters for stability analysis, set by th write_outputs_to_HDF5::Bool = true HDF5_filename::String = "euler.h5" force_wv_symmetry::Bool = true + eps_res::Float64 = 1e-2 end @kwdef mutable struct FourFitVars diff --git a/src/DCON/Sing.jl b/src/DCON/Sing.jl index 65679e7f..021ea57f 100644 --- a/src/DCON/Sing.jl +++ b/src/DCON/Sing.jl @@ -706,15 +706,15 @@ function sing_der!(du::Array{ComplexF64,3}, u::Array{ComplexF64,3}, # Compute singfac = 1 / (m - nq) odet.q = Spl.spline_eval!(equil.sq, psieval)[4] singfac_den = ComplexF64.((intr.mlow:intr.mhigh) .- odet.q .* (intr.nlow:intr.nhigh)') - # ϵ_res = 1e-10 - ϵ_res = 1e-2 + # eps_res = 1e-10 + # eps_res = 1e-1 # Add 1j * eps * norm to denominators near resonance to avoid singularity for (j, n) in enumerate(intr.nlow:intr.nhigh) mres = clamp(round(Int, odet.q * n), intr.mlow, intr.mhigh) i = mres - intr.mlow + 1 x = abs(singfac_den[i, j]) - gauss_norm = exp(- (x / (5 * ϵ_res))^2) - singfac_den[i, j] += 1im * ϵ_res * gauss_norm + gauss_norm = exp(- (x / (5 * ctrl.eps_res))^2) + singfac_den[i, j] += 1im * ctrl.eps_res * gauss_norm end odet.singfac_vec .= vec(1.0 ./ singfac_den)