PerturbedEquilibrium - BUGFIX - change singular coupling calculations to match Fortran, using du_store instead of splined xi' #348
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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 | ||
|
Collaborator
Author
There was a problem hiding this comment. Minor: this comment states the reason slightly wrong. Generated by Claude Code |
||
|
|
||
| # 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`; 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) | ||
|
Collaborator
Author
There was a problem hiding this comment. Observation (non-blocking): slot 1 of This array is Also: Generated by Claude Code |
||
| 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) | ||
|
|
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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 | ||
|
Comment on lines
+73
to
+78
Collaborator
Author
There was a problem hiding this comment. Two minor robustness notes on the stencil selection:
Neither is blocking given the current save density; flagging for awareness. Generated by Claude Code |
||
|
|
||
| # 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, | ||
| 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, | ||
| 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 | ||
| continue | ||
| end | ||
|
|
||
|
|
@@ -201,36 +307,46 @@ function compute_singular_coupling_metrics!( | |
| lpsi = sing_surf.psifac - spot_psi | ||
| rpsi = sing_surf.psifac + spot_psi | ||
|
|
||
| # 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] | ||
| # 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 | ||
|
|
||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The old block kept a comment explaining why this call exists; the
elseifbranch now looks like dead code (itsdu_bufferresult is discarded). Suggest restoring a one-liner, e.g.# Gaussian reduction modified u: recompute for the odet.ud side effect (Fortran ode_output.f sing_der before euler.bin write)— otherwise a future cleanup pass will delete it.Generated by Claude Code