Top level organization of unified Perturbed Equilibrium code

[logan-nc]

This branch has already merged in the changes from #118, and supersedes that PR.

I think it should be merged quickly, and not necessarily wait for full PerturbedEquilibrium functionality / benchmarking. This is because the organizational changes will quickly cause many conflicts if others are working on feature branches from the old develop. Let's get the organization reviewed by @matt-pharr and merged, then keep going on the details of the perturbed equilibrium

[logan-nc]

@claude please summarize the changes of this PR. Be succinct. Do it only once and do not repeat it when new commits are pushed unless explicitly asked to do so.

[claude]

Claude finished @logan-nc's task in 53s —— View job

---

Summary of PR Changes

This PR reorganizes JPEC to support unified perturbed equilibrium calculations with a cleaner architecture:

Major Structural Changes:

• Module Rename: DCON → ForceFreeStates (more descriptive name for stability analysis)
• New Modules: Added PerturbedEquilibrium/ (945 lines) and ForcingTerms/ for response calculations
• Config Consolidation: Replaced separate equil.toml, dcon.toml, and vac.in files with unified jpec.toml containing [Equilibrium], [WALL], [ForceFreeStates], and [PerturbedEquilibrium] sections
• Main Refactor: JPEC.jl now orchestrates full workflow: equilibrium → stability → perturbed equilibrium response

Code Organization:

• Moved all DCON functionality to src/ForceFreeStates/ with consistent naming
• Added perturbed equilibrium skeleton with data structures (PerturbedEquilibriumControl, PerturbedEquilibriumState) and core functions for plasma response and singular coupling calculations
• Updated Vacuum.jl to support new data flow

Cleanup:

• Removed 14,741-line test GEQDSK file from benchmarks (moved to external data)
• Deleted legacy config files across all examples and tests
• Removed obsolete DCON/Main.jl after migration

Documentation:

• Added PERTURBED_EQUILIBRIUM_PLAN.md and GPEC_RESPONSE_IMPLEMENTATION.md design docs
• Updated CLAUDE.md to reflect new module structure

The changes are primarily organizational to establish clean foundations for perturbed equilibrium work, with minimal functional changes to existing equilibrium/stability code.

---

Tasks

• Read key changed files to understand reorganization
• Provide succinct summary of changes
• Post summary to comment

---

[matt-pharr]

@logan-nc Would you like my review of this now?

P.S. thanks claude :^)

[matt-pharr]

Also, for what it is worth, this is what I came up with after discussing with copilot of a good layout for pull requests to implement GPEC functionality.

GPEC Julia Conversion Plan

Executive Summary

This document outlines an incremental conversion plan for GPEC (General Perturbed Equilibrium Code) from Fortran to Julia, organized into bite-sized pull requests. The conversion eliminates file I/O between DCON and GPEC by passing data directly through memory via Julia structs, with parallelization-first design.

Key Goal: Get a working GPEC Julia implementation that can:

1. Receive force-free states from DCON in memory
2. Apply external forcing terms (from simple file or eigenmode)
3. Calculate plasma response properties (Lambda, permeability, reluctance)
4. Compute resonant metrics at singular surfaces
5. Output perturbed equilibrium profiles

Approach: Implement via small, testable PRs. Each PR is self-contained and independently mergeable.

Pull Request Roadmap

PR1: Directory Restructuring & Module Stubs

Goal: Set up new directory structure with minimal working modules

Tasks:

1. Rename src/DCON/ → src/ForceFreeStates/
2. Update module name: DCON → ForceFreeStates
3. Create src/ForcingTerms/ directory with stub module
4. Create src/PerturbedEquilibrium/ directory with stub module
5. Move src/ForceFreeStates/Main.jl → src/Main.jl
6. Update all imports and module references
7. Verify all existing tests pass

Deliverable: New directory structure compiles, existing functionality unchanged

Files to create:

• src/ForcingTerms/ForcingTerms.jl (empty module stub)
• src/PerturbedEquilibrium/PerturbedEquilibrium.jl (empty module stub)

Files to modify:

• Rename src/DCON/ → src/ForceFreeStates/
• Move and update src/Main.jl
• Update src/JPEC.jl exports

---

PR2: ForcingTerms Module - Data Structures

Goal: Define data structures for external field specification

Tasks:

1. Create ExternalField struct
2. Create ForcingConfig struct for TOML parsing
3. Add docstrings and examples
4. Add unit tests for data structure creation

Deliverable: ForcingTerms module with data types, no functionality yet

Files to create:

• src/ForcingTerms/ExternalField.jl
• test/test_forcing_terms_structs.jl

Example struct:

@kwdef struct ExternalField
    finmn::Vector{ComplexF64}   # External field spectrum [mpert]
    source::String              # "file", "eigenmode", or "coil"
    mpert::Int
    mlow::Int
    mhigh::Int
end

---

PR3: ForcingTerms Module - Field Specification

Goal: Implement reading external fields from files or eigenmodes

Tasks:

1. Implement read_field_from_ascii(filename, mpert, mlow, mhigh)
2. Implement extract_eigenmode(odet, mode_idx, intr)
3. Implement set_forcing_terms(config, intr, odet)
4. Add comprehensive tests with example input files
5. Document ASCII file format

Deliverable: Can read/specify external forcing terms

Files to create:

• src/ForcingTerms/FieldSpecification.jl
• test/test_forcing_terms_io.jl
• test/fixtures/external_field_example.dat

ASCII format: Simple m real(b_x) imag(b_x) per line

---

PR4: PerturbedEquilibrium Module - Data Structures

Goal: Define data structures for plasma response and resonant metrics

Tasks:

1. Create PlasmaResponseData struct
2. Create ResonantMetrics struct
3. Create PerturbedEquilibriumConfig struct for TOML parsing
4. Add docstrings
5. Add unit tests

Deliverable: PerturbedEquilibrium module with data types

Files to create:

• src/PerturbedEquilibrium/Structs.jl
• test/test_perturbed_equilibrium_structs.jl

Example structs:

@kwdef struct PlasmaResponseData
    mpert::Int
    Lambda::Matrix{ComplexF64}   # Plasma inductance
    P::Matrix{ComplexF64}        # Permeability
    reluct::Matrix{ComplexF64}   # Reluctance
    P_evals::Vector{ComplexF64}
    P_evecs::Matrix{ComplexF64}
end

@kwdef struct ResonantMetrics
    msing::Int
    mpert::Int
    coupling::Array{ComplexF64, 3}  # [5 models, msing, mpert]
    sing_field::Array{ComplexF64, 2} # [msing, mpert]
end

---

PR5: PlasmaResponse - Inductance Calculations

Goal: Implement Lambda (plasma inductance) calculation

Tasks:

1. Convert gpresp_pinduct → compute_plasma_inductance(odet, intr, equil)
2. Implement helper functions for eigenfunction reconstruction
3. Add comprehensive unit tests against Fortran benchmarks
4. Document the physics/math

Deliverable: Can compute plasma inductance Lambda from force-free eigenfunctions

Files to create:

• src/PerturbedEquilibrium/PlasmaResponse.jl (partial)
• test/test_plasma_inductance.jl
• test/benchmarks/fortran_lambda_benchmark.h5

Key function:

function compute_plasma_inductance(
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium
) -> Matrix{ComplexF64}
    # Convert from gpresp_pinduct
end

---

PR6: PlasmaResponse - Permeability & Reluctance

Goal: Complete plasma response calculations

Tasks:

1. Implement compute_permeability(Lambda, L) (P = Lambda * inv(L))
2. Implement compute_reluctance(P) (reluct = inv(P - I))
3. Implement compute_fixed_boundary_inductance(equil, intr) for fixed boundary case
4. Implement calculate_plasma_response_properties() main entry point
5. Add eigendecomposition of P
6. Test against Fortran benchmarks

Deliverable: Complete plasma response calculation (Lambda, P, reluctance)

Files to modify:

• src/PerturbedEquilibrium/PlasmaResponse.jl (complete)
• test/test_plasma_response.jl

Key function:

function calculate_plasma_response_properties(
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    vac::Union{VacuumData, Nothing},
    external_field::ExternalField,
    config::Dict
) -> PlasmaResponseData
end

---

PR7: ResonantMetrics - Singular Coupling

Goal: Implement singular coupling calculations

Tasks:

1. Convert gpout_singcoup → compute_coupling_single()
2. Implement 5 coupling models (flux, current, island width, penetrated field, delta)
3. Implement calculate_resonant_metrics() with parallel loop
4. Add tests against Fortran benchmarks
5. Document physics of each coupling model

Deliverable: Can compute singular coupling at rational surfaces

Files to create:

• src/PerturbedEquilibrium/ResonantMetrics.jl (partial)
• test/test_singular_coupling.jl

Key function:

function calculate_resonant_metrics(
    resp::PlasmaResponseData,
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    external_field::ExternalField,
    config::Dict
) -> ResonantMetrics
    # Parallel over (ising, imode) pairs
end

---

PR8: FieldProfiles - Profile Calculations

Goal: Implement at least one perturbed field profile type

Tasks:

1. Implement evaluate_perturbed_field(psi, ...) for eigenfunction reconstruction
2. Implement compute_normal_field_profile() (xbnormal) as primary profile
3. Optionally add compute_clebsch_profile() (xclebsch) or compute_mod_b_perturbation() (pmodb)
4. Add tests
5. Document profile physics

Deliverable: Can compute at least xbnormal profile

Files to create:

• src/PerturbedEquilibrium/FieldProfiles.jl
• test/test_field_profiles.jl

Key function:

function calculate_perturbed_equil_profiles(
    resp::PlasmaResponseData,
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    external_field::ExternalField,
    profile_type::String
) -> Dict
end

---

PR9: Output Module & HDF5 Writing

Goal: Implement HDF5 output for all results

Tasks:

1. Implement write_perturbed_equil_to_file() with comprehensive HDF5 structure
2. Include all response matrices, resonant metrics, profiles
3. Add metadata (parameters, timestamps, etc.)
4. Add tests for output file structure
5. Document HDF5 schema

Deliverable: Single comprehensive HDF5 output file

Files to create:

• src/PerturbedEquilibrium/Output.jl
• test/test_output.jl

HDF5 structure:

perturbed_equilibrium.h5
├── plasma_response/
│   ├── Lambda
│   ├── permeability
│   ├── reluctance
│   └── P_eigenvalues
├── resonant/
│   ├── coupling [5, msing, mpert]
│   ├── singular_field [msing, mpert]
│   ├── psi [msing]
│   └── q [msing]
├── profiles/
│   └── xbnormal [...]
└── metadata/
    ├── timestamp
    ├── forcing_source
    └── config

---

PR10: Main.jl Integration & End-to-End Pipeline

Goal: Connect all modules in working end-to-end pipeline

Tasks:

1. Update src/Main.jl to call all new modules in sequence
2. Add TOML configuration parsing for perturbed_equilibrium.toml
3. Add verbose output and timing
4. Add error handling and validation
5. Create end-to-end integration test
6. Document usage with example

Deliverable: Complete working MWE pipeline from Main.jl

Files to modify:

• src/Main.jl
• test/test_integration_mwe.jl
• Create examples/perturbed_equilibrium.toml
• Create examples/README.md with usage

Main.jl pseudo-code:

function main(path::String="./")
    println("JPEC: Julia Perturbed Equilibrium Code")
    
    # 1. Setup equilibrium
    equil = Equilibrium.setup_equilibrium(...)
    
    # 2. Compute force-free states
    odet, intr, vac = ForceFreeStates.compute_force_free_states(...)
    
    # 3. Set forcing terms
    config = load_perturbed_equilibrium_config(...)
    external_field = ForcingTerms.set_forcing_terms(config, intr, odet)
    
    # 4. Compute perturbed equilibrium
    resp = PerturbedEquilibrium.calculate_plasma_response_properties(...)
    metrics = PerturbedEquilibrium.calculate_resonant_metrics(...)
    profiles = PerturbedEquilibrium.calculate_perturbed_equil_profiles(...)
    
    # 5. Write output
    PerturbedEquilibrium.write_perturbed_equil_to_file(...)
    
    println("Normal termination.")
end

---

Post-MWE Pull Requests (Optional Enhancement)

PR11: CoilGeometry Module (Post-MWE)

Goal: Add coil geometry calculations for coil_flag

Tasks:

1. Create src/ForcingTerms/CoilGeometry.jl
2. Implement coil-to-plasma field calculation
3. Add tests and examples

---

PR12: Additional Profile Types (Post-MWE)

Goal: Add xclebsch, pmodb, and other profile types

Tasks:

1. Extend FieldProfiles.jl with additional profile calculations
2. Add tests for each type

---

PR13: (R,Z) Grid Outputs (Post-MWE)

Goal: Add spatial grid outputs for visualization

Tasks:

1. Create src/PerturbedEquilibrium/GridOutputs.jl
2. Implement (R,Z) grid evaluation
3. Add to HDF5 output

---

PR14: Coordinate Transformations (Post-MWE)

Goal: Support jac_in/jac_out different from equilibrium

Tasks:

1. Create src/PerturbedEquilibrium/CoordinateTransforms.jl
2. Implement Hamada, PEST, Boozer, etc. transforms
3. Add tests

---

PR15: Parallelization Optimization (Post-MWE)

Goal: Optimize parallel performance

Tasks:

1. Profile code to identify bottlenecks
2. Optimize memory access patterns
3. Tune thread counts
4. Benchmark speedup
5. Document parallel configuration

---

Configuration Simplification

Minimal TOML for MWE (perturbed_equilibrium.toml):

[FORCING_TERMS]
# External field specification
data_flag = true                # Read from file
infile = "external_field.dat"   # Simple ASCII: m, real(b_x), imag(b_x)
# OR
mode_flag = true                # Use force-free eigenmode
mode = 0                        # Which eigenmode (0 = most unstable)

[PLASMA_RESPONSE]
fixed_boundary_flag = false     # true = fixed boundary, false = free boundary
sing_spot = 5e-4                # Spot size around singularities for integration
sing_npsi = 100                 # Number of psi points for singular field
reg_flag = false                # Regularization (usually false)
reg_spot = 5e-2                 # Regularization parameter

[OUTPUT]
singcoup_flag = true            # Compute singular coupling
singfld_flag = true             # Compute singular field  
profile_type = "xbnormal"       # "xbnormal", "xclebsch", or "pmodb"
verbose = true
timeit = true

[PARALLEL]
num_threads = 0                 # 0 = use all available

What to Skip for MWE (add later):

• ❌ Coil geometry calculations (coil_flag)
• ❌ Complex field specifications (displacement_flag, surfmn format)
• ❌ Galerkin eigenmodes (gal_flag)
• ❌ Coordinate transformations (jac_in/jac_out != equilibrium type)
• ❌ Filtering (filter_types, filter_modes)
• ❌ Chebyshev modes (chebyshev_flag)
• ❌ Complex output formats (NetCDF, ASCII tables, binary)
• ❌ (R,Z) grid outputs (add after MWE works)
• ❌ Multiple singular field analysis variants

Configuration Simplification

Minimal TOML for MWE (perturbed_equilibrium.toml):

[FORCING_TERMS]
# External field specification
data_flag = true                # Read from file
infile = "external_field.dat"   # Simple ASCII: m, b_x (one per line)
# OR
mode_flag = true                # Use force-free eigenmode
mode = 0                        # Which eigenmode (0 = most unstable)

[PLASMA_RESPONSE]
fixed_boundary_flag = false     # true = fixed boundary, false = free boundary
sing_spot = 5e-4                # Spot size around singularities for integration
sing_npsi = 100                 # Number of psi points for singular field
reg_flag = false                # Regularization (usually false)
reg_spot = 5e-2                 # Regularization parameter

[OUTPUT]
singcoup_flag = true            # Compute singular coupling
singfld_flag = true             # Compute singular field  
profile_type = "xbnormal"       # "xbnormal", "xclebsch", or "pmodb"
verbose = true
timeit = true

[PARALLEL]
num_threads = 0                 # 0 = use all available

Deprecated/Hardcoded for MWE:

• jac_in, jac_out → Always use equilibrium coordinate type, error if different requested
• resp_index → Always 0
• jsurf_in, jsurf_out, tmag_in, tmag_out, mlim_out → Placeholders, not implemented
• mthsurf → Use sensible default from equilibrium
• netcdf_flag, ascii_flag, bin_flag, bin_2d_flag → Always false, HDF5 only
• filter_flag, filter_types, filter_modes → Not implemented for MWE
• Output flags for different field types → Pick one (xbnormal) for MWE

Implementation Phases

Phase 1: Directory Restructuring & Data Structures

Objective: Set up new directory structure and minimal data types

Tasks:

1. Rename DCON → ForceFreeStates

   • Move src/DCON/ to src/ForceFreeStates/
   • Update module names and exports
   • Update Main.jl imports

2. Create ForcingTerms directory

   src/ForcingTerms/
   ├── ForcingTerms.jl          # Module exports
   ├── ExternalField.jl         # Data structure for finmn
   └── FieldSpecification.jl    # Read from file or eigenmode
   

3. Create PerturbedEquilibrium directory

   src/PerturbedEquilibrium/
   ├── PerturbedEquilibrium.jl  # Module exports
   ├── PlasmaResponse.jl        # Stub for Phase 2
   ├── ResonantMetrics.jl       # Stub for Phase 2
   ├── FieldProfiles.jl         # Stub for Phase 3
   └── Output.jl                # Stub for Phase 3
   

4. Move Main.jl to top level

   • Move src/ForceFreeStates/Main.jl to src/Main.jl
   • Make it orchestrate: Equilibrium, ForceFreeStates, ForcingTerms, PerturbedEquilibrium

5. Define minimal data structures

   # In ForcingTerms/ExternalField.jl
   @kwdef struct ExternalField
       finmn::Vector{ComplexF64}   # External field spectrum
       source::String              # "file", "eigenmode", or "coil"
       mpert::Int
   end
   
   # In PerturbedEquilibrium/PlasmaResponse.jl
   @kwdef struct PlasmaResponseData
       mpert::Int
       Lambda::Matrix{ComplexF64}   # Plasma inductance
       P::Matrix{ComplexF64}        # Permeability
       reluct::Matrix{ComplexF64}   # Reluctance
       P_evals::Vector{ComplexF64}  # Permeability eigenvalues
       P_evecs::Matrix{ComplexF64}  # Permeability eigenvectors
   end
   
   # In PerturbedEquilibrium/ResonantMetrics.jl
   @kwdef struct ResonantMetrics
       msing::Int
       mpert::Int
       coupling::Array{ComplexF64, 3}  # [5 models, msing, mpert]
       sing_field::Array{ComplexF64, 2} # [msing, mpert]
   end

Deliverable: New directory structure, minimal types, stubs compiling

Phase 2: Core Plasma Response

Objective: Implement plasma response calculations (Lambda, P, reluctance)

Key Functions (in PlasmaResponse.jl):

"""
    calculate_plasma_response_properties(
        odet::OdeState,
        intr::DconInternal, 
        equil::PlasmaEquilibrium,
        vac::Union{VacuumData, Nothing},
        external_field::ExternalField
    ) -> PlasmaResponseData

Main entry point for plasma response calculations.
Corresponds to gpout_response in Fortran GPEC.
"""
function calculate_plasma_response_properties(
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    vac::Union{VacuumData, Nothing},
    external_field::ExternalField
)
    # 1. Compute plasma inductance Lambda from force-free eigenfunctions
    Lambda = compute_plasma_inductance(odet, intr, equil)
    
    # 2. Get surface inductance L from vacuum data
    L = vac !== nothing ? vac.wt : compute_fixed_boundary_inductance(equil, intr)
    
    # 3. Compute permeability P = Lambda * inv(L)
    P = compute_permeability(Lambda, L)
    
    # 4. Compute reluctance = inv(P - I) with power normalization
    reluct = compute_reluctance(P)
    
    # 5. Eigendecomposition of P
    P_evals, P_evecs = eigen(Hermitian(P))
    
    return PlasmaResponseData(
        mpert = intr.mpert,
        Lambda = Lambda,
        P = P,
        reluct = reluct,
        P_evals = P_evals,
        P_evecs = P_evecs
    )
end

# Helper functions (convert from gpresp.f):
compute_plasma_inductance(odet, intr, equil)
compute_fixed_boundary_inductance(equil, intr)  # For fixed_boundary_flag=true
compute_permeability(Lambda, L)
compute_reluctance(P)

Conversion Notes:

• gpresp_pinduct → compute_plasma_inductance: Build Lambda from eigenfunctions in OdeState
• gpresp_sinduct → Use VacuumData.wt directly (already computed)
• gpresp_permeab → compute_permeability: Matrix operations
• gpresp_reluct → compute_reluctance: Matrix inversion and normalization
• gpresp_eigen → Use Julia's eigen() from LinearAlgebra

Deliverable: Working plasma response calculation, tested against Fortran

Phase 3: Resonant Metrics & Profiles

Objective: Compute singular coupling/field and at least one profile type

Part A: Resonant Metrics (in ResonantMetrics.jl):

"""
    calculate_resonant_metrics(
        resp::PlasmaResponseData,
        odet::OdeState,
        intr::DconInternal,
        equil::PlasmaEquilibrium,
        external_field::ExternalField,
        config::Dict
    ) -> ResonantMetrics

Combines singular coupling and singular field analysis.
Simpler than separate singcoup_flag and singfld_flag.
"""
function calculate_resonant_metrics(
    resp::PlasmaResponseData,
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    external_field::ExternalField,
    config::Dict
)
    sing_spot = config[:sing_spot]
    sing_npsi = config[:sing_npsi]
    
    # Pre-allocate
    coupling = zeros(ComplexF64, 5, intr.msing, intr.mpert)
    sing_field = zeros(ComplexF64, intr.msing, intr.mpert)
    
    # Parallel loop over singular surfaces and modes
    Threads.@threads for idx in CartesianIndices((1:intr.msing, 1:intr.mpert))
        ising, imode = idx[1], idx[2]
        
        # Compute for unit field at this mode
        coupling[:, ising, imode] = compute_coupling_single(
            ising, imode, resp, odet, intr, equil, sing_spot
        )
        
        # Singular field at rational surface
        sing_field[ising, imode] = extract_singular_field(
            ising, imode, resp, odet, intr, equil, sing_spot
        )
    end
    
    return ResonantMetrics(
        msing = intr.msing,
        mpert = intr.mpert,
        coupling = coupling,
        sing_field = sing_field
    )
end

# Helper functions (convert from gpout.f::gpout_singcoup and gpout_singfld):
compute_coupling_single(ising, imode, resp, odet, intr, equil, sing_spot)
extract_singular_field(ising, imode, resp, odet, intr, equil, sing_spot)
compute_surface_properties(psi, equil)  # Area, shear, etc.

Part B: Field Profiles (in FieldProfiles.jl):

"""
    calculate_perturbed_equil_profiles(
        resp::PlasmaResponseData,
        odet::OdeState,
        intr::DconInternal,
        equil::PlasmaEquilibrium,
        external_field::ExternalField,
        profile_type::String
    ) -> Dict

Compute perturbed equilibrium profiles.
Start with xbnormal for MWE, add others later.
"""
function calculate_perturbed_equil_profiles(
    resp::PlasmaResponseData,
    odet::OdeState,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    external_field::ExternalField,
    profile_type::String
)
    if profile_type == "xbnormal"
        return compute_normal_field_profile(resp, odet, intr, equil, external_field)
    elseif profile_type == "xclebsch"
        return compute_clebsch_profile(resp, odet, intr, equil, external_field)
    elseif profile_type == "pmodb"
        return compute_mod_b_perturbation(resp, odet, intr, equil, external_field)
    else
        error("Unknown profile type: $profile_type")
    end
end

# For MWE, implement one:
compute_normal_field_profile(resp, odet, intr, equil, external_field)

Deliverable: Resonant metrics and at least one profile type working

Phase 4: I/O & Integration

Objective: Connect everything in Main.jl and output to HDF5

Part A: ForcingTerms implementation:

# In ForcingTerms/FieldSpecification.jl

"""
    set_forcing_terms(
        config::Dict,
        intr::DconInternal,
        odet::Union{OdeState, Nothing} = nothing
    ) -> ExternalField

Read external field from file or select eigenmode.
"""
function set_forcing_terms(
    config::Dict,
    intr::DconInternal,
    odet::Union{OdeState, Nothing} = nothing
)
    if config[:data_flag]
        # Read from simple ASCII file: m, b_x (one per line)
        finmn = read_field_from_file(config[:infile], intr.mpert, intr.mlow, intr.mhigh)
        return ExternalField(finmn=finmn, source="file", mpert=intr.mpert)
        
    elseif config[:mode_flag]
        # Use eigenmode from DCON
        mode_idx = config[:mode]
        finmn = extract_eigenmode(odet, mode_idx, intr)
        return ExternalField(finmn=finmn, source="eigenmode", mpert=intr.mpert)
        
    else
        error("Must set either data_flag or mode_flag")
    end
end

function read_field_from_file(filename::String, mpert::Int, mlow::Int, mhigh::Int)
    # Simple format: one line per mode
    # m  real(b_x)  imag(b_x)
    finmn = zeros(ComplexF64, mpert)
    open(filename) do f
        for line in eachline(f)
            parts = split(line)
            m = parse(Int, parts[1])
            idx = m - mlow + 1
            if 1 <= idx <= mpert
                finmn[idx] = parse(Float64, parts[2]) + im * parse(Float64, parts[3])
            end
        end
    end
    return finmn
end

Part B: Output implementation:

# In PerturbedEquilibrium/Output.jl

"""
    write_perturbed_equil_to_file(
        filename::String,
        resp::PlasmaResponseData,
        metrics::ResonantMetrics,
        profiles::Dict,
        intr::DconInternal,
        equil::PlasmaEquilibrium,
        external_field::ExternalField
    )

Write all perturbed equilibrium data to single HDF5 file.
"""
function write_perturbed_equil_to_file(
    filename::String,
    resp::PlasmaResponseData,
    metrics::ResonantMetrics,
    profiles::Dict,
    intr::DconInternal,
    equil::PlasmaEquilibrium,
    external_field::ExternalField
)
    h5open(filename, "w") do h5
        # Plasma response
        h5["plasma_response/Lambda"] = resp.Lambda
        h5["plasma_response/permeability"] = resp.P
        h5["plasma_response/reluctance"] = resp.reluct
        h5["plasma_response/P_eigenvalues"] = resp.P_evals
        
        # Resonant metrics
        h5["resonant/coupling"] = metrics.coupling
        h5["resonant/singular_field"] = metrics.sing_field
        h5["resonant/psi"] = [s.psifac for s in intr.sing]
        h5["resonant/q"] = [s.q for s in intr.sing]
        
        # Profiles
        for (key, val) in profiles
            h5["profiles/$key"] = val
        end
        
        # External field
        h5["forcing/finmn"] = external_field.finmn
        h5["forcing/source"] = external_field.source
    end
end

Part C: Main.jl integration:

# In src/Main.jl (top level)

function main(path::String="./")
    println("=" ^ 60)
    println("JPEC: Julia Perturbed Equilibrium Code")
    println("=" ^ 60)
    
    # 1. Setup equilibrium
    println("\n[1/5] Setting up equilibrium...")
    equil = Equilibrium.setup_equilibrium(joinpath(path, "equil.toml"))
    
    # 2. Compute force-free states (DCON)
    println("\n[2/5] Computing force-free states...")
    odet, intr, vac = ForceFreeStates.compute_force_free_states(path)
    
    # 3. Set forcing terms
    println("\n[3/5] Setting forcing terms...")
    config = load_config(joinpath(path, "perturbed_equilibrium.toml"))
    external_field = ForcingTerms.set_forcing_terms(config, intr, odet)
    
    # 4. Compute perturbed equilibrium
    println("\n[4/5] Computing perturbed equilibrium...")
    
    # 4a. Plasma response properties
    t_start = time()
    resp = PerturbedEquilibrium.calculate_plasma_response_properties(
        odet, intr, equil, vac, external_field
    )
    if config[:verbose]
        println("  → Plasma response: $(round(time()-t_start, digits=3))s")
    end
    
    # 4b. Resonant metrics
    t_start = time()
    metrics = PerturbedEquilibrium.calculate_resonant_metrics(
        resp, odet, intr, equil, external_field, config
    )
    if config[:verbose]
        println("  → Resonant metrics: $(round(time()-t_start, digits=3))s")
    end
    
    # 4c. Field profiles
    t_start = time()
    profiles = PerturbedEquilibrium.calculate_perturbed_equil_profiles(
        resp, odet, intr, equil, external_field, config[:profile_type]
    )
    if config[:verbose]
        println("  → Field profiles: $(round(time()-t_start, digits=3))s")
    end
    
    # 5. Write output
    println("\n[5/5] Writing output...")
    PerturbedEquilibrium.write_perturbed_equil_to_file(
        joinpath(path, "perturbed_equilibrium.h5"),
        resp, metrics, profiles, intr, equil, external_field
    )
    
    println("\n" * "=" ^ 60)
    println("Normal termination.")
    println("=" ^ 60)
end

Deliverable: Full MWE pipeline working end-to-end

• Transformation matrices

3. Perturbed Field Components:

   • xsp_mn, xsp1_mn - Radial displacement and derivative
   • xss_mn, xms_mn - Perpendicular displacements
   • bwp_mn, bwt_mn, bwz_mn - Contravariant B-field components
   • Similar arrays for covariant, normal, tangential components

4. Singular Coupling Data:

   • singcoup - Coupling matrices (5 types: flux, current, island width, penetrated field, delta)
   • singcoup_out_vecs - Right singular vectors in output coordinates
   • fldflxmat - Field to flux transformation matrix

5. Coil/External Field Data:

   • coil_indmat - Coil inductance matrix
   • finmn, foutmn - External/total flux
   • xspmn - External field spectrum

Step 1.2: Design GpecStructs.jl

Create analogous structures to DconStructs.jl:

# File: src/GPEC/GpecStructs.jl

"""
    GpecControl
    
Configuration parameters for GPEC, read from gpec.toml
"""
@kwdef mutable struct GpecControl
    # Input/Output settings
    verbose::Bool = true
    dir_path::String = ""
    
    # Parallelization settings
    num_threads::Int = Threads.nthreads()  # Number of threads to use

### Phase 5: Post-MWE Enhancements

**Objective**: Add features beyond MWE

**Priority enhancements**:
1. **Coil geometry** (`coil_flag`) - Generate fields from coil positions
2. **Torque coupling matrix** - For rotation studies
3. **(R,Z) grid outputs** - For visualization
4. **Additional profile types** - xclebsch, pmodb, etc.
5. **Coordinate transformations** - jac_in/jac_out support
6. **Filtering** - Mode filtering for numerical noise

**Lower priority** (implement if needed):
- Displacement-based input (`displacement_flag`)
- Complex field formats (`surfmn`)
- Galerkin eigenmode support (`gal_flag`)
- Chebyshev modes
- Multiple output formats (NetCDF, ASCII)

### Phase 6: Parallelization Optimization

**Objective**: Maximize parallel performance

**Key optimizations**:
1. **Profile serial code** - Identify bottlenecks
2. **Tune thread counts** - Test scaling 1, 2, 4, 8, 16 threads
3. **Optimize memory access** - Use views, reduce allocations
4. **Parallel BLAS** - Ensure eigendecompositions use all cores
5. **Benchmark vs Fortran** - Compare single-thread and multi-thread

**Expected performance**:
- Resonant metrics: 8-16x speedup (embarrassingly parallel)
- Overall: 3-6x speedup

## Detailed Conversion Notes

### Fortran → Julia Module Mapping

| Fortran File | → | Julia Module/File | Priority |
|--------------|---|-------------------|----------|
| `idcon.f` | → | Removed (in-memory handoff) | - |
| `gpresp.f` | → | `PerturbedEquilibrium/PlasmaResponse.jl` | **MWE** |
| `gpout.f::gpout_response` | → | `PlasmaResponse.jl::calculate_plasma_response_properties` | **MWE** |
| `gpout.f::gpout_singcoup` | → | `ResonantMetrics.jl::calculate_resonant_metrics` | **MWE** |
| `gpout.f::gpout_singfld` | → | `ResonantMetrics.jl` (merged with singcoup) | **MWE** |
| `gpeq.f::gpeq_sol` | → | `FieldProfiles.jl::evaluate_perturbed_field` | **MWE** |
| `gpeq.f::gpeq_normal` | → | `FieldProfiles.jl::compute_normal_field_profile` | **MWE** |
| `gpcoil.f` | → | `ForcingTerms/CoilGeometry.jl` | Post-MWE |
| `gpeq.f::gpeq_fcoords` | → | `CoordinateTransforms.jl` | Post-MWE |
| `gpout.f::gpout_*` (others) | → | `Output.jl` or future files | Post-MWE |

### Key Simplifications for MWE

1. **No coordinate transformations**: Always use equilibrium coordinate type
2. **No filtering**: Implement later if needed for noise reduction
3. **Single output file**: HDF5 only, comprehensive data structure
4. **Merged singular analysis**: Combine singcoup and singfld for cleaner code
5. **Simple forcing**: ASCII table or eigenmode, no complex formats
6. **Hardcoded defaults**: Fewer TOML parameters, sensible built-in values

### Data That DCON Must Retain

After Phase 2-3 implementation, if we discover GPEC needs additional data from DCON:

**Checklist to verify DCON retains** (most likely already available):
- ✅ Eigenfunction solutions at all radial points (`u_store`)
- ✅ Singular surface locations and properties (`sing`)
- ✅ Asymptotic coefficients (`ca_l`, `ca_r`)
- ✅ Vacuum response matrices if computed (`wt`, `wv`, `wt0`)
- ✅ Mode numbers (`mlow`, `mhigh`, `mpert`)
- ✅ Equilibrium splines (`sq`, `rzphi`)

**If something is missing**: Add field to appropriate struct (OdeState, DconInternal, VacuumData) and populate during DCON calculation.

## Testing Strategy

### Unit Tests (per module)
```julia
# test/test_plasma_response.jl
@testset "Plasma Response" begin
    # Load test case
    odet, intr, equil, vac = load_test_case("solovev_n1")
    
    # Compute response
    resp = calculate_plasma_response_properties(odet, intr, equil, vac, dummy_field)
    
    # Check against Fortran benchmark
    fortran_resp = load_fortran_benchmark("solovev_n1_response.h5")
    
    @test resp.Lambda ≈ fortran_resp.Lambda rtol=1e-12
    @test resp.P ≈ fortran_resp.P rtol=1e-12
    @test resp.reluct ≈ fortran_resp.reluct rtol=1e-12
    @test resp.P_evals ≈ fortran_resp.P_evals rtol=1e-10
end

Integration Tests (full workflow)

# test/test_full_workflow.jl
@testset "Full MWE Workflow" begin
    path = "test/test_case_solovev_n1/"
    
    # Run full pipeline
    main(path)
    
    # Check output file exists and has expected structure
    @test isfile(joinpath(path, "perturbed_equilibrium.h5"))
    
    h5open(joinpath(path, "perturbed_equilibrium.h5")) do h5
        @test haskey(h5, "plasma_response/Lambda")
        @test haskey(h5, "resonant/coupling")
        @test haskey(h5, "profiles/xbnormal")
    end
    
    # Compare with Fortran reference
    julia_results = load_julia_results(joinpath(path, "perturbed_equilibrium.h5"))
    fortran_results = load_fortran_results(joinpath(path, "reference.h5"))
    
    compare_results(julia_results, fortran_results, rtol=1e-10)
end

Physics Validation Tests

@testset "Physics Validation" begin
    # Test 1: Energy conservation
    @test check_energy_conservation(resp, odet)
    
    # Test 2: Hermiticity of matrices
    @test ishermitian(resp.Lambda)
    @test ishermitian(resp.P)
    
    # Test 3: Singular coupling reciprocity
    @test check_reciprocity(metrics.coupling)
    
    # Test 4: Field continuity
    @test check_field_continuity(profiles)
end

Success Criteria

For MWE (Phases 1-4):

1. ✅ Restructured directories (ForceFreeStates, ForcingTerms, PerturbedEquilibrium)
2. ✅ Can read external field from simple ASCII file or eigenmode
3. ✅ Computes Lambda, P, reluctance matching Fortran to machine precision
4. ✅ Computes singular coupling/field at rational surfaces
5. ✅ Outputs at least one profile type (xbnormal)
6. ✅ Single HDF5 output file with all results
7. ✅ Runs end-to-end from Main.jl
8. ✅ Comprehensive unit and integration tests pass

For Full Implementation:

1. ✅ All MWE criteria above
2. ✅ Coil geometry calculations
3. ✅ (R,Z) grid outputs for visualization
4. ✅ Multiple profile types
5. ✅ Parallel speedup 3-6x overall, 8-16x for resonant metrics
6. ✅ Documentation and examples

Next Steps (Immediate)

1. Review and merge this plan
2. Phase 1 starts: Create branches, restructure directories

   git checkout -b feature/perturbed-equilibrium-phase1
   mv src/DCON src/ForceFreeStates
   mkdir src/ForcingTerms
   mkdir src/PerturbedEquilibrium
   mv src/ForceFreeStates/Main.jl src/Main.jl

3. Create stub files with minimal types and function signatures
4. Test compilation of new structure
5. Begin Phase 2: Start converting gpresp.f → PlasmaResponse.jl

Key Architectural Decisions Summary

1. In-memory data flow: Zero file I/O between modules (DCON → GPEC via structs)
2. MWE first: Get working code fast, add features incrementally
3. Single HDF5 output: Deprecate all other formats initially
4. Simplified config: Minimal TOML, hardcode sensible defaults
5. Merged analysis: Combine singcoup+singfld for cleaner implementation
6. Parallelization-ready: Thread-safe from start, optimize in Phase 6
7. Directory structure: Mirrors physics organization (ForceFreeStates, ForcingTerms, PerturbedEquilibrium)
8. Top-level Main.jl: Orchestrates all modules

References

• Fortran GPEC source: src/GPEC/gpec_fortran/ (This is important for using claude to code so that it can have a reference of the original fortran code)
• DCON Julia (to be renamed ForceFreeStates): src/DCON/
• Equilibrium module: src/Equilibrium/
• Splines module: src/Splines/
• Vacuum module: src/Vacuum/

---

Last Updated: After incorporating detailed prioritization feedback (MWE approach, directory restructuring, simplified configuration)

Deprecated/Hardcoded for MWE:

• jac_in, jac_out → Always use equilibrium coordinate type, error if different requested
• resp_index → Always 0
• jsurf_in, jsurf_out, tmag_in, tmag_out, mlim_out → Not implemented for MWE
• mthsurf → Use sensible default from equilibrium
• netcdf_flag, ascii_flag, bin_flag → Always false, HDF5 only
• filter_flag, filter_types, filter_modes → Not implemented for MWE
• Output flags for different field types → Always compute, configurable via profile_type

---

Testing Strategy for Each PR

Each PR should include appropriate tests:

Unit Tests: Test individual functions in isolation

@testset "PlasmaInductance" begin
    odet, intr, equil = load_test_case("solovev_n1")
    Lambda = compute_plasma_inductance(odet, intr, equil)
    
    # Compare with Fortran benchmark
    fortran_Lambda = load_benchmark("lambda_solovev_n1.h5")
    @test Lambda ≈ fortran_Lambda rtol=1e-12
end

Integration Tests: Test module interactions

@testset "PlasmaResponse" begin
    resp = calculate_plasma_response_properties(...)
    
    # Check internal consistency
    @test ishermitian(resp.Lambda)
    @test ishermitian(resp.P)
    
    # Check against Fortran
    @test resp.P_evals ≈ fortran_evals rtol=1e-10
end

End-to-End Tests (PR10 only): Test full pipeline

@testset "FullPipeline" begin
    main("test/test_case_solovev/")
    
    # Check output file
    @test isfile("test/test_case_solovev/perturbed_equilibrium.h5")
    
    # Validate results
    results = load_results("test/test_case_solovev/perturbed_equilibrium.h5")
    compare_with_fortran(results, "test/benchmarks/fortran_reference.h5")
end

---

Fortran → Julia Conversion Reference

Fortran File/Subroutine	→	Julia Module/Function	PR
idcon.f	→	Removed (in-memory handoff)	PR1
gpresp.f::gpresp_pinduct	→	PlasmaResponse::compute_plasma_inductance	PR5
gpresp.f::gpresp_permeab	→	PlasmaResponse::compute_permeability	PR6
gpresp.f::gpresp_reluct	→	PlasmaResponse::compute_reluctance	PR6
gpout.f::gpout_response	→	PlasmaResponse::calculate_plasma_response_properties	PR6
gpout.f::gpout_singcoup	→	ResonantMetrics::calculate_resonant_metrics	PR7
gpout.f::gpout_singfld	→	Merged with singcoup in ResonantMetrics	PR7
gpeq.f::gpeq_sol	→	FieldProfiles::evaluate_perturbed_field	PR8
gpeq.f::gpeq_normal	→	FieldProfiles::compute_normal_field_profile	PR8
gpout.f::output_routines	→	Output::write_perturbed_equil_to_file	PR9

---

Key Data That DCON Must Retain

Already available (verify during implementation):

• ✅ Eigenfunction solutions at all radial points (u_store)
• ✅ Singular surface locations and properties (sing)
• ✅ Asymptotic coefficients (ca_l, ca_r)
• ✅ Vacuum response matrices if computed (wt, wv, wt0)
• ✅ Mode numbers (mlow, mhigh, mpert)
• ✅ Equilibrium splines (sq, rzphi)

If something is missing during implementation: Add field to appropriate struct (OdeState, DconInternal, VacuumData) and populate during ForceFreeStates calculation. No file format changes needed with in-memory flow.

---

Success Criteria

For MWE (PR1-10):

1. ✅ All 10 PRs merged and tests passing
2. ✅ Restructured directories (ForceFreeStates, ForcingTerms, PerturbedEquilibrium)
3. ✅ Can read external field from ASCII file or eigenmode
4. ✅ Computes Lambda, P, reluctance matching Fortran to machine precision
5. ✅ Computes singular coupling/field at rational surfaces
6. ✅ Outputs at least xbnormal profile
7. ✅ Single HDF5 output with all results
8. ✅ Runs end-to-end from Main.jl with example config
9. ✅ Comprehensive test coverage (unit + integration tests)
10. ✅ Documentation and usage examples

For Full Implementation (with optional PR11-15):

1. ✅ All MWE criteria above
2. ✅ Coil geometry calculations (PR11)
3. ✅ Multiple profile types (PR12)
4. ✅ (R,Z) grid outputs (PR13)
5. ✅ Coordinate transformations (PR14)
6. ✅ Parallel optimization with 3-6x speedup (PR15)

---

Key Architectural Decisions

1. In-memory data flow: Zero file I/O between modules (ForceFreeStates → ForcingTerms → PerturbedEquilibrium)
2. PR-based development: 10 bite-sized PRs for MWE, each independently mergeable
3. Test-driven: Each PR includes comprehensive tests before merging
4. Single HDF5 output: Deprecate all other formats initially
5. Simplified config: Minimal TOML parameters (~10 vs 40+ in Fortran)
6. Merged analysis: Combine singcoup+singfld for cleaner implementation
7. Parallelization-ready: Thread-safe from start, optimize later
8. Directory structure: Mirrors physics organization (ForceFreeStates, ForcingTerms, PerturbedEquilibrium)
9. Top-level Main.jl: Orchestrates all modules

---

PR Dependency Graph

PR1 (Directory Restructuring)
  ├─> PR2 (ForcingTerms Structs)
  │     └─> PR3 (ForcingTerms I/O)
  └─> PR4 (PerturbedEquilibrium Structs)
        ├─> PR5 (PlasmaResponse: Lambda)
        │     └─> PR6 (PlasmaResponse: P, reluct)
        ├─> PR7 (ResonantMetrics)
        └─> PR8 (FieldProfiles)
              └─> PR9 (Output)
                    └─> PR10 (Main.jl Integration)

Optional Post-MWE:
PR10 ─> PR11, PR12, PR13, PR14, PR15 (any order)

---

Quick Start for AI Agent

To implement PR1:

1. git checkout -b feature/perturbed-equilibrium-pr1
2. mv src/DCON src/ForceFreeStates
3. Update module names in all files
4. Create stub modules for ForcingTerms and PerturbedEquilibrium
5. Move Main.jl to top level
6. Update imports
7. Run tests to verify nothing broken
8. Commit and create PR

To implement PR2:

1. Create src/ForcingTerms/ExternalField.jl
2. Define ExternalField struct with docstrings
3. Create test/test_forcing_terms_structs.jl
4. Add tests for struct creation
5. Commit and create PR

... and so on for each PR.

---

References

• Fortran GPEC source: src/GPEC/gpec_fortran/
• ForceFreeStates (formerly DCON): src/ForceFreeStates/ (after PR1)
• Equilibrium module: src/Equilibrium/
• Splines module: src/Splines/
• Vacuum module: src/Vacuum/

[logan-nc]

wow! extensive! I think there are some good ideas from copilot that claude didn't implement. so let's flag the ones we want.

yes @matt-pharr please review now. I know it's frustrating to review during active dev, but I think it's necessary with this top level organization. you can gloss over the perturbed equilibrium calculation details for now - I'll have claude work on getting to an actual singular coupling calc today (I'm iterating with it slowly in between family activities).

Note, I like the idea of a top level math utilities directory... but the fourier transforms there are currently 2-10x SLOWER than the FFTW package. This contradicts @priyanshlunia findings that the vacuum transforms were faster than FFTW. idk why. So maybe we just want to use FFTW in the end for gpec stuff? the custom ones are still needed for vacuum because they include an nqd phase shift. whatever we do, other things that can go in the utilities might be the grid packing options like ldp grid making, etc.