Complex regularization of singularities instead of enforcing ideal crossing

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).

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
@@ -397,7 +399,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

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))

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)')
+    # 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 * ctrl.eps_res))^2)
+        singfac_den[i, j] += 1im * ctrl.eps_res * gauss_norm
+    end
+
+    odet.singfac_vec .= vec(1.0 ./ singfac_den)
 
     # kinetic stuff - skip for now
     if false #(TODO: kin_flag)

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