diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 1f8fe3bc..48e1e704 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -645,6 +645,7 @@ function cross_ideal_singular_surf!( odet.q_store[odet.step] = odet.q odet.u_store[:, :, :, odet.step] = odet.u odet.ud_store[:, :, :, odet.step] = odet.ud + odet.du_store[:, :, :, odet.step] = du1 odet.step += 1 end @@ -696,6 +697,7 @@ function cross_kinetic_singular_surf!( odet.q_store[odet.step] = odet.q odet.u_store[:, :, :, odet.step] = odet.u odet.ud_store[:, :, :, odet.step] = odet.ud + odet.du_store[:, :, :, odet.step] = du1 odet.step += 1 end @@ -757,12 +759,6 @@ function integrate_el_region!( compute_solution_norms!(integrator.u, odet, ctrl, intr, false) - # If Gaussian reduction modified u, recompute ud to keep it consistent. - # Matches Fortran ode_output.f which calls sing_der before writing euler.bin. - if odet.new - sing_der!(du_buffer, integrator.u, integrator.p, integrator.t) - end - # Save near segment boundaries (symmetric, in q not psi) and every Nth step. # The step-count fallback (== 1) guarantees the first step is always saved # even for near-degenerate segments where q_range ≈ 0. @@ -775,11 +771,15 @@ function integrate_el_region!( if odet.step >= size(odet.u_store, 4) resize_storage!(odet) end + sing_der!(du_buffer, integrator.u, integrator.p, integrator.t) odet.psi_store[odet.step] = integrator.t @views odet.u_store[:, :, :, odet.step] .= integrator.u odet.q_store[odet.step] = odet.q @views odet.ud_store[:, :, :, odet.step] .= odet.ud + @views odet.du_store[:, :, :, odet.step] .= du_buffer odet.step += 1 + elseif odet.new + sing_der!(du_buffer, integrator.u, integrator.p, integrator.t) end end @@ -794,10 +794,12 @@ function integrate_el_region!( if odet.step >= size(odet.u_store, 4) resize_storage!(odet) end + sing_der!(du_buffer, sol.u[end], (ctrl, equil, ffit, intr, odet, chunk), sol.t[end]) odet.psi_store[odet.step] = sol.t[end] @views odet.u_store[:, :, :, odet.step] .= sol.u[end] odet.q_store[odet.step] = odet.q @views odet.ud_store[:, :, :, odet.step] .= odet.ud + @views odet.du_store[:, :, :, odet.step] .= du_buffer odet.step += 1 end @@ -1035,6 +1037,10 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal) odet.ud_store[:, :, 1, istep] .= gauss_buffer mul!(gauss_buffer, odet.ud_store[:, :, 2, istep], transforms[:, :, ifix]) odet.ud_store[:, :, 2, istep] .= gauss_buffer + mul!(gauss_buffer, odet.du_store[:, :, 1, istep], transforms[:, :, ifix]) + odet.du_store[:, :, 1, istep] .= gauss_buffer + mul!(gauss_buffer, odet.du_store[:, :, 2, istep], transforms[:, :, ifix]) + odet.du_store[:, :, 2, istep] .= gauss_buffer end jfix = kfix + 1 end diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 7f727ea5..7f52cdf5 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -346,6 +346,8 @@ end # Kinetic torque matrix splines: 6 components ktmats::Vector{S} = [_empty_series_interp_complex(numpert_total^2, itp_opts) for _ in 1:6] + kinetic_populated::Bool = false # set by make_kinetic_matrix; fmats_lower/kmats are then not the EL operators + # Pre-computed FKG kinetic matrices (populated by make_kinetic_matrix) f0mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) pmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts) @@ -492,7 +494,9 @@ and a small set of temporary matrices and factors used to compute singular-layer - `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape `(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver). - - `ud_store::Array{ComplexF64,4}` - Stored derivatives of the solution at each saved step with same shape as `u_store`. + - `ud_store::Array{ComplexF64,4}` - Stored derivatives of the solution at each saved step with same shape as `u_store`; slot 2 holds the Clebsch Ξ_s, not du₂/dψ. + + - `du_store::Array{ComplexF64,4}` - du₁/dψ and du₂/dψ at the accepted point of each saved step, populated only by the serial EL path. - `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`). @@ -560,6 +564,7 @@ and a small set of temporary matrices and factors used to compute singular-layer q_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init) u_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init) ud_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init) + du_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init) crit_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init) ca_r::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing) ca_l::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 52acb9ba..cb2d3555 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -130,6 +130,10 @@ and data dumping. odet.ud_store[:, :, 1, istep] .= tmp_mat mul!(tmp_mat, odet.ud_store[:, :, 2, istep], coeffs) odet.ud_store[:, :, 2, istep] .= tmp_mat + mul!(tmp_mat, odet.du_store[:, :, 1, istep], coeffs) + odet.du_store[:, :, 1, istep] .= tmp_mat + mul!(tmp_mat, odet.du_store[:, :, 2, istep], coeffs) + odet.du_store[:, :, 2, istep] .= tmp_mat end # Write energies to screen diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index a1434b85..da1202ba 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -257,6 +257,7 @@ function _compute_fkg_matrices!( ffit.amats = cubic_interp(xs, Series(ak_flat); ffit.itp_opts...) ffit.bmats = cubic_interp(xs, Series(bk_flat); ffit.itp_opts...) ffit.cmats = cubic_interp(xs, Series(ck_flat); ffit.itp_opts...) + ffit.kinetic_populated = true return nothing end diff --git a/src/ForceFreeStates/Utils.jl b/src/ForceFreeStates/Utils.jl index 99410935..95f420d1 100644 --- a/src/ForceFreeStates/Utils.jl +++ b/src/ForceFreeStates/Utils.jl @@ -12,18 +12,21 @@ function resize_storage!(odet::OdeState) # Allocate new arrays u_new = Array{ComplexF64,4}(undef, odet.numpert_total, odet.numpert_total, 2, newlen) ud_new = Array{ComplexF64,4}(undef, odet.numpert_total, odet.numpert_total, 2, newlen) + du_new = zeros(ComplexF64, odet.numpert_total, odet.numpert_total, 2, newlen) psi_new = Vector{Float64}(undef, newlen) q_new = Vector{Float64}(undef, newlen) # Copy old data u_new[:, :, :, 1:odet.step] = odet.u_store[:, :, :, 1:odet.step] ud_new[:, :, :, 1:odet.step] = odet.ud_store[:, :, :, 1:odet.step] + du_new[:, :, :, 1:odet.step] = odet.du_store[:, :, :, 1:odet.step] psi_new[1:odet.step] = odet.psi_store[1:odet.step] q_new[1:odet.step] = odet.q_store[1:odet.step] # Replace old arrays odet.u_store = u_new odet.ud_store = ud_new + odet.du_store = du_new odet.psi_store = psi_new odet.q_store = q_new end @@ -40,4 +43,5 @@ function trim_storage!(odet::OdeState) resize!(odet.q_store, odet.step) odet.u_store = odet.u_store[:, :, :, 1:odet.step] odet.ud_store = odet.ud_store[:, :, :, 1:odet.step] + odet.du_store = odet.du_store[:, :, :, 1:odet.step] end diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index 5f1fe83b..f94eabbc 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -127,7 +127,7 @@ function compute_perturbed_equilibrium( if vac_data === nothing @warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section." else - compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl) + compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, ffit) end end diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index d9cd1ebc..a7290d4f 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -43,6 +43,106 @@ end # Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf). _reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) +""" + _solution_at(psi, psi_surf, odet, equil, ffs_intr, nstep) -> (u1, du1) + +Evaluate Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution: cubic Hermite for the value, +Lagrange cubic through the singfac-weighted `du_store` nodes for the derivative. Works for +any formulation; ideal runs use the more accurate `_el_solution_at`. The Ξ′ stencil stays +on `psi`'s side of the singular surface `psi_surf`, where the stored solution is +discontinuous. Requires a populated `odet.du_store`. +""" +function _solution_at( + psi::Float64, + psi_surf::Float64, + odet::OdeState, + equil::Equilibrium.PlasmaEquilibrium, + ffs_intr::ForceFreeStatesInternal, + nstep::Int +) + il, ir, _ = _psi_bracket(odet.psi_store, psi, nstep) + psi_a, psi_b = odet.psi_store[il], odet.psi_store[ir] + + u1_a = @view odet.u_store[:, :, 1, il] + u1_b = @view odet.u_store[:, :, 1, ir] + du1_a = @view odet.du_store[:, :, 1, il] + du1_b = @view odet.du_store[:, :, 1, ir] + + u1_e = _hermite_cubic_val(u1_a, u1_b, du1_a, du1_b, psi_a, psi_b, psi) + + # Same-side candidate nodes around the bracket, trimmed to the 4 nearest psi. + side = sign(psi - psi_surf) + idxs = [j for j in max(1, il - 3):min(nstep, ir + 3) if sign(odet.psi_store[j] - psi_surf) == side] + while length(idxs) > 4 + abs(odet.psi_store[idxs[1]] - psi) > abs(odet.psi_store[idxs[end]] - psi) ? popfirst!(idxs) : pop!(idxs) + end + + # Interpolate singfac·Ξ′ and divide the ~1/singfac pole back out at psi. + singfac_at(p) = vec([m - equil.profiles.q_spline(p) * n for m in ffs_intr.mlow:ffs_intr.mhigh, n in ffs_intr.nlow:ffs_intr.nhigh]) + du1_e = zeros(ComplexF64, size(odet.du_store, 1), size(odet.du_store, 2)) + for k in eachindex(idxs) + w = 1.0 + for j in eachindex(idxs) + j == k && continue + w *= (psi - odet.psi_store[idxs[j]]) / (odet.psi_store[idxs[k]] - odet.psi_store[idxs[j]]) + end + du1_e .+= (w .* singfac_at(odet.psi_store[idxs[k]])) .* @view(odet.du_store[:, :, 1, idxs[k]]) + end + du1_e ./= singfac_at(psi) + + return u1_e, du1_e +end + +""" + _el_solution_at(psi, odet, ffit, equil, ffs_intr, nstep) -> (u1, du1) + +Evaluate Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution via the ideal Euler-Lagrange +relation Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) [Glasser 2016 eqs. 22-24], with u₁, u₂ Hermite- +interpolated to `psi`. Only valid for ideal runs where `ffit.fmats_lower` and `kmats` +generated the solution. Requires a populated `odet.du_store`. +""" +function _el_solution_at( + psi::Float64, + odet::OdeState, + ffit::FourFitVars, + equil::Equilibrium.PlasmaEquilibrium, + ffs_intr::ForceFreeStatesInternal, + nstep::Int +) + npert = ffs_intr.numpert_total + il, ir, _ = _psi_bracket(odet.psi_store, psi, nstep) + psi_a, psi_b = odet.psi_store[il], odet.psi_store[ir] + + u1_a = @view odet.u_store[:, :, 1, il] + u1_b = @view odet.u_store[:, :, 1, ir] + u2_a = @view odet.u_store[:, :, 2, il] + u2_b = @view odet.u_store[:, :, 2, ir] + du1_a = @view odet.du_store[:, :, 1, il] + du1_b = @view odet.du_store[:, :, 1, ir] + du2_a = @view odet.du_store[:, :, 2, il] + du2_b = @view odet.du_store[:, :, 2, ir] + + hint = Ref(1) + kmat = Matrix{ComplexF64}(undef, npert, npert) + + u1_e = _hermite_cubic_val(u1_a, u1_b, du1_a, du1_b, psi_a, psi_b, psi) + u2_e = _hermite_cubic_val(u2_a, u2_b, du2_a, du2_b, psi_a, psi_b, psi) + + # Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) with Q⁻¹ = diag(1/(m − n·q)) + q_e = equil.profiles.q_spline(psi) + singfac_inv = vec([1.0 / (m - q_e * n) for m in ffs_intr.mlow:ffs_intr.mhigh, n in ffs_intr.nlow:ffs_intr.nhigh]) + fmat_lower = Matrix{ComplexF64}(undef, npert, npert) + ffit.fmats_lower(vec(fmat_lower), psi; hint=hint) + ffit.kmats(vec(kmat), psi; hint=hint) + du1_e = u2_e .* singfac_inv + du1_e .-= kmat * u1_e + ldiv!(LowerTriangular(fmat_lower), du1_e) + ldiv!(UpperTriangular(fmat_lower'), du1_e) + du1_e .*= singfac_inv + + return u1_e, du1_e +end + """ compute_singular_coupling_metrics!( state::PerturbedEquilibriumState, @@ -51,7 +151,8 @@ _reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) vac_data::VacuumData, ffs_intr::ForceFreeStatesInternal, intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl + ctrl::PerturbedEquilibriumControl, + ffit::FourFitVars ) Compute singular layer coupling matrices and applied resonant vectors. @@ -82,7 +183,8 @@ function compute_singular_coupling_metrics!( vac_data::VacuumData, ffs_intr::ForceFreeStatesInternal, intr::PerturbedEquilibriumInternal, - ctrl::PerturbedEquilibriumControl + ctrl::PerturbedEquilibriumControl, + ffit::FourFitVars ) ctrl.verbose && @info "Computing singular coupling metrics (GPEC method)" @@ -159,6 +261,10 @@ function compute_singular_coupling_metrics!( # @threads region. psi_store_all = ForceFreeStates_results.psi_store nstep = ForceFreeStates_results.step + # ξ′ evaluation preference: ideal EL relation, then interpolated stored RHS for kinetic, + # then chord slope for paths that never populate du_store like the gal-matched odet. + use_du_store = any(!iszero, ForceFreeStates_results.du_store) + use_el = use_du_store && !ffit.kinetic_populated _blas_nthreads = BLAS.get_num_threads() BLAS.set_num_threads(1) try @@ -169,7 +275,7 @@ function compute_singular_coupling_metrics!( resnum = findfirst(j -> intr.m_modes[j] == m_res && intr.n_modes[j] == nn, 1:numpert_total) if resnum === nothing - @warn "Could not find index for resonant mode (m=$m_res, n=$nn)" maxlog=1 + @warn "Could not find index for resonant mode (m=$m_res, n=$nn)" maxlog = 1 continue end @@ -201,36 +307,46 @@ function compute_singular_coupling_metrics!( lpsi = sing_surf.psifac - spot_psi rpsi = sing_surf.psifac + spot_psi - # Interpolate u and du/dψ at lpsi and rpsi using stored ODE solution. - # Value (u): Hermite cubic using ud_store for smoother interpolation than linear. - # Derivative (ud): chord slope from u_store only — see comment below. - il_l, ir_l, _ = _psi_bracket(psi_store_all, lpsi, nstep) - il_r, ir_r, _ = _psi_bracket(psi_store_all, rpsi, nstep) - - u_node = ForceFreeStates_results.u_store - ud_node = ForceFreeStates_results.ud_store - ua_l = u_node[resnum, :, 1, il_l] - ub_l = u_node[resnum, :, 1, ir_l] - ua_r = u_node[resnum, :, 1, il_r] - ub_r = u_node[resnum, :, 1, ir_r] - dua_l = ud_node[resnum, :, 1, il_l] - dub_l = ud_node[resnum, :, 1, ir_l] - dua_r = ud_node[resnum, :, 1, il_r] - dub_r = ud_node[resnum, :, 1, ir_r] - - psi_il_l = psi_store_all[il_l] - psi_ir_l = psi_store_all[ir_l] - psi_il_r = psi_store_all[il_r] - psi_ir_r = psi_store_all[ir_r] - - u_l = _hermite_cubic_val(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi) - u_r = _hermite_cubic_val(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi) - # Derivative (ud): chord slope from u_store only — ud_store can be systematically off near - # outer surfaces (q=4/5) where the ODE solution varies rapidly; near-cancellation in bwp1 - # then amplifies even a ~10% ud_store error into a large Delta' error. Chord slope avoids - # this by using only u values, which are accurately stored by the ODE integrator. - ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l) - ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r) + # Evaluate u and dξ/dψ at lpsi and rpsi from the stored ODE solution + if !use_du_store + il_l, ir_l, _ = _psi_bracket(psi_store_all, lpsi, nstep) + il_r, ir_r, _ = _psi_bracket(psi_store_all, rpsi, nstep) + + u_node = ForceFreeStates_results.u_store + ud_node = ForceFreeStates_results.ud_store + ua_l = u_node[resnum, :, 1, il_l] + ub_l = u_node[resnum, :, 1, ir_l] + ua_r = u_node[resnum, :, 1, il_r] + ub_r = u_node[resnum, :, 1, ir_r] + dua_l = ud_node[resnum, :, 1, il_l] + dub_l = ud_node[resnum, :, 1, ir_l] + dua_r = ud_node[resnum, :, 1, il_r] + dub_r = ud_node[resnum, :, 1, ir_r] + + psi_il_l = psi_store_all[il_l] + psi_ir_l = psi_store_all[ir_l] + psi_il_r = psi_store_all[il_r] + psi_ir_r = psi_store_all[ir_r] + + u_l = _hermite_cubic_val(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi) + u_r = _hermite_cubic_val(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi) + ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l) + ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r) + elseif use_el + u1m_l, du1m_l = _el_solution_at(lpsi, ForceFreeStates_results, ffit, equil, ffs_intr, nstep) + u1m_r, du1m_r = _el_solution_at(rpsi, ForceFreeStates_results, ffit, equil, ffs_intr, nstep) + u_l = u1m_l[resnum, :] + u_r = u1m_r[resnum, :] + ud_l = du1m_l[resnum, :] + ud_r = du1m_r[resnum, :] + else + u1m_l, du1m_l = _solution_at(lpsi, sing_surf.psifac, ForceFreeStates_results, equil, ffs_intr, nstep) + u1m_r, du1m_r = _solution_at(rpsi, sing_surf.psifac, ForceFreeStates_results, equil, ffs_intr, nstep) + u_l = u1m_l[resnum, :] + u_r = u1m_r[resnum, :] + ud_l = du1m_l[resnum, :] + ud_r = du1m_r[resnum, :] + end q_l = equil.profiles.q_spline(lpsi) q1_l = equil.profiles.q_deriv(lpsi) @@ -242,10 +358,11 @@ function compute_singular_coupling_metrics!( jump_vec = Vector{ComplexF64}(undef, numpert_total) for k in 1:numpert_total ck = @view C_coeffs[:, k] - xsp_l = dot(u_l, ck) - xsp1_l = dot(ud_l, ck) - xsp_r = dot(u_r, ck) - xsp1_r = dot(ud_r, ck) + # unconjugated contraction, dot would conjugate u + xsp_l = transpose(u_l) * ck + xsp1_l = transpose(ud_l) * ck + xsp_r = transpose(u_r) * ck + xsp1_r = transpose(ud_r) * ck bwp1_l = 2π * im * chi1 * (singfac_l * xsp1_l - nn * q1_l * xsp_l) bwp1_r = 2π * im * chi1 * (singfac_r * xsp1_r - nn * q1_r * xsp_r) jump_vec[k] = bwp1_r - bwp1_l