Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions docs/development/references.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ The ForceFreeStates module (ideal MHD stability analysis) implements methods fro

- **Glasser (2016)**: "The direct criterion of Newcomb for the ideal MHD stability of an axisymmetric toroidal plasma"
- Location: `docs/resources/2016-Glasser-The_direct_criterion_of_Newcomb_for_the_ideal_MHD_stability_of_an_axisymmetric_toroidal_plasma.pdf`
- Published: Physics of Plasmas **23**, 112506 (2016)
- Published: Physics of Plasmas **23**, 072505 (2016)
- **Describes: FUNDAMENTAL PAPER - Newcomb's criterion for ideal MHD stability (current implementation)**

- **Glasser (2018)**: "A Riccati solution for the ideal MHD plasma response with applications to real-time stability control"
Expand Down Expand Up @@ -97,7 +97,7 @@ GPEC will eventually implement resistive MHD stability analysis based on:

- **Glasser (2016)**: "Computation of resistive instabilities by matched asymptotic expansions"
- Location: `docs/resources/2016-Glasser-Computation_of_resistive_instabilities_by_matched_asymptotic_expansions-compressed.pdf`
- Published: Physics of Plasmas **23**, 072505 (2016)
- Published: Physics of Plasmas **23**, 112506 (2016)
- Describes: Resistive stability analysis and Δ' calculation via matched asymptotic expansions

- **Glasser (2018)**: "A robust solution for the resistive MHD toroidal Δ′ matrix in near real-time"
Expand Down
16 changes: 9 additions & 7 deletions src/PerturbedEquilibrium/FieldReconstruction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,9 @@ field perturbations in mode space, following the GPEC gpeq module approach.
Displacement components from ODE integration (u_store/du_store/xi_s_store):
- ξ_ψ: radial displacement (u_store[:,:,1,:])
- dξ_ψ/dψ: radial derivative (du_store)
- ξ_s: toroidal displacement (xi_s_store, Glasser 2016 eq. 18)
- ξ_s: surface displacement ξ_s = J(ξ×B)·(∇θ×∇ζ) = χ₁(q ξ^θ - ξ^ζ), the in-surface
displacement across B (binormal Clebsch component); it is not a toroidal displacement
(xi_s_store; Glasser 2016 Phys. Plasmas 23, 072505, eqs. 9 and 18)

Contravariant perturbed field from ideal MHD (matches Fortran gpeq_sol):
b^ψ = χ₁·(m - n·q)·2πi·ξ_ψ
Expand Down Expand Up @@ -227,13 +229,13 @@ alpha = flux_matrix \\ response_vector
Then sum eigenmode contributions at each radial point (matches Fortran gpeq_sol):
xi_psi[ipsi, :] = u_store[:, :, 1, ipsi] * alpha # Ξ_ψ
xi_psi1[ipsi, :] = du_store[:, :, ipsi] * alpha # dΞ_ψ/dψ
xi_s[ipsi, :] = xi_s_store[:, :, ipsi] * alpha # Ξ_s (toroidal, Glasser 2016 eq. 18)
xi_s[ipsi, :] = xi_s_store[:, :, ipsi] * alpha # Ξ_s (surface displacement, Glasser 2016 eqs. 9, 18)

# Returns

- `xi_psi_modes`: Radial displacement ξ_ψ(ψ, m) [npsi, numpert_total]
- `xi_psi1_modes`: Radial derivative dξ_ψ/dψ(ψ, m) [npsi, numpert_total]
- `xi_s_modes`: Toroidal displacement ξ_s(ψ, m) = -A⁻¹(B·dξ_ψ/dψ + C·ξ_ψ) [npsi, numpert_total]
- `xi_s_modes`: Surface displacement ξ_s(ψ, m) = χ₁(q ξ^θ - ξ^ζ) = -A⁻¹(B·dξ_ψ/dψ + C·ξ_ψ) [npsi, numpert_total]
"""
function sum_eigenmode_contributions(
response_vector::Vector{ComplexF64},
Expand Down Expand Up @@ -262,7 +264,7 @@ function sum_eigenmode_contributions(
mul!(view(xi_psi1_modes, ipsi, :),
@view(solution.du_store[:, :, ipsi]),
alpha)
# xi_s_store = Ξ_s = -A⁻¹(B·Ξ'_ψ + C·Ξ_ψ) (toroidal displacement, Glasser 2016 eq. 18)
# xi_s_store = Ξ_s = -A⁻¹(B·Ξ'_ψ + C·Ξ_ψ) (surface displacement, Glasser 2016 eqs. 9, 18)
mul!(view(xi_s_modes, ipsi, :),
@view(solution.xi_s_store[:, :, ipsi]),
alpha)
Expand Down Expand Up @@ -422,7 +424,7 @@ function compute_clebsch_displacements(
# xms = -(A\B)*xmp1 - (A\C)*xsp
xsp_vec = view(xi_psi_modes, ipsi, :)
mul!(xms_vec, bmat, xmp1_vec) # xms = B*xmp1
mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im) # xms += C*xsp
mul!(xms_vec, cmat_buf, xsp_vec, 1.0 + 0.0im, 1.0 + 0.0im) # xms += C*xsp
# ideal A only; the kinetic A is non-Hermitian and is treated separately above
amat_fact = cholesky!(Hermitian(amat, :L))
ldiv!(amat_fact, xms_vec) # xms = A\(B*xmp1 + C*xsp)
Expand Down Expand Up @@ -831,7 +833,7 @@ function compute_b_n_xi_n_modes(
jd = jacs .* delpsis

for in in 1:npert
blk = ((in - 1) * mpert + 1):(in * mpert)
blk = ((in-1)*mpert+1):(in*mpert)
# IDFT: mode space → theta space
bwp_fun = phase_fwd * b_psi_modes[ipsi, blk] # [mthsurf] — b^ψ (not J-weighted)
xwp_fun = phase_fwd * xwp_modes[ipsi, blk] # [mthsurf] — J·ξ^ψ (J-weighted from gpeq_contra)
Expand Down Expand Up @@ -1020,7 +1022,7 @@ function _apply_rzphi_transform(
zeta_fn = buf.ze

for in in 1:npert
blk = ((in - 1) * mpert + 1):(in * mpert)
blk = ((in-1)*mpert+1):(in*mpert)
# Inverse DFT: modes → theta-space (in place)
Utilities.FourierTransforms.inverse_transform!(psi_fun, ft, view(psi_input, ipsi, blk))
Utilities.FourierTransforms.inverse_transform!(theta_fn, ft, view(theta_input, ipsi, blk))
Expand Down
Loading