Clarify Euler-Lagrange integration bounds and refactor for explicitness

[Copilot]

Plan for Clarifying Euler-Lagrange Integration Bounds

Based on the discussion in the issue, this PR implements the following changes:

• Rename Ode.jl to EulerLagrange.jl
• Update module include in DCON.jl to reference new filename
• Rename functions with clear, descriptive names
• Create set_el_integration_bounds! function to extract boundary setting logic
• Refactor initialize_el_at_axis to only initialize state
• Refactor cross_ideal_singular_surf to extract bound setting logic
• Replace odet.next with clearer boolean odet.needs_crossing
• Make integration bounds explicit arguments to integrate_el_region
• Update all call sites in Main.jl
• Update test file references
• Fix bug in set_el_integration_bounds
• Replace while loop with for loop over pre-computed integration chunks
• Create chunk_el_integration_bounds() function to pre-compute all chunks
• Add IntegrationChunk struct to hold chunk metadata
• Remove set_el_integration_bounds! (replaced by chunking function)
• NEW: Separate singular surface struct into basic parameters (SingType) and asymptotics (SingAsymptotics)
• NEW: Compute asymptotics on-demand in cross_ideal_singular_surf! instead of pre-computing
• NEW: Remove sing_scan! function (no longer needed)
• NEW: Refactor sing_vmat! → compute_sing_asymptotics to return struct
• NEW: Update sing_get_ua and sing_get_ca to use SingAsymptotics
• Run full test suite to ensure no regressions (all tests passing)

Summary of Changes

This PR makes the Euler-Lagrange integration bounds more clear and explicit, and clarifies that asymptotic calculations are specific to ideal DCON:

1. Integration Bounds Clarification (previous commits)

• Renamed functions for clarity (ode_run → eulerlagrange_integration, etc.)
• Pre-computed integration chunks with chunk_el_integration_bounds()
• Made integration bounds explicit arguments
• Replaced while loop with for loop over chunks

2. On-Demand Asymptotic Computation (this commit)

• Separated SingType into two structs:

  • SingType - Basic singular surface parameters (psifac, m, n, q, q1, di)
  • SingAsymptotics - Asymptotic expansion data (alpha, vmat, mmat, r1, r2, power, etc.)

• On-demand computation: Asymptotics are now computed during each singular surface crossing in cross_ideal_singular_surf! instead of being pre-computed for all surfaces upfront

• Benefits:

  • Makes it crystal clear that asymptotics are only computed and used in ideal DCON, not kinetic
  • Follows the pattern: compute when needed, use immediately, discard after crossing
  • Avoids storing large asymptotic arrays that may never be used
  • Easier to understand: asymptotics computation is localized to where it's used

Implementation Details

• compute_sing_asymptotics() replaces sing_vmat!() and returns SingAsymptotics struct
• sing_get_ua() and sing_get_ca() now take SingAsymptotics and SingType separately
• sing_scan!() removed (was only used to pre-compute asymptotics)
• compute_sing_mmat!() and sing_solve!() updated to work with individual arrays instead of mutating singp

This architecture supports future parallelization and makes the code much clearer.

Original prompt

---

This section details on the original issue you should resolve

<issue_title>Make ode integration bounds more clear</issue_title>
<issue_description>Right now the legacy code has us doing a lot of things hidden deep inside of places one might not expect. For example, the ode_ideal_cross doesn't just cross, it also finds the next singular surface and sets an internal struct parameter that then determines the bounds of the next ode_step integration chunk. Lots of badness here:

• ode_step doesn't do one step, it actually integrates across a big chunk of psin
• The bounds of the integration are a mystery. One has to click in from ode_run to ode_step to see ode.psimax is used and then search the repo for psimax to find that is gets set in cross, which then one can search to see that it gets called in ode_run after each ode_step and then, maybe, if you are very familiar with the code, you get it.
• ode_ideal_cross does more than the name implies

I propose we:

• Rename all ode* things to be eulerlagrange* things
• Rename ode_step to eulerlagrange_integration and make the bounds explicit arguments at that level
  • Note, this could come in handy when implementing STRIDE, which integrates backwards and forwards in many (parallelized) chunks
  • An alternative to explicit bounds arguments would be to pull out a separate set_el_integration_bounds function from the crossing/init functions and call it right before each integration chunk - so hopefully it is clear where that is done
• Rename ode_run to be something like form_eulerlagrange_basis or some such

This is NOT pressing.
Thoughts @jhalpern30?</issue_description>

<agent_instructions>Please create a PR that addresses this issue based on the discussion below. Make sure to consider which points were agreed upon in the discussion and those that were abandoned in favor of other changes (i.e. function names, implementation strategy)</agent_instructions>

Comments on the Issue (you are @copilot in this section)

@jhalpern30 - I agree `ode` can be renamed to `eulerlagrange` things. Renaming the file `EulerLagrange.jl` would also be fine. Unlike DCON, I don't think every function needs to be `eulerlagrange_...` since it just makes the names long. I'd rather just do all functions in the file are descriptive, and if they're in this file, you know they are related to the EL equation solve. - I would rather rename `ode_run` to `eulerlagrange_integration`, since this is really the function that integrates the entire equation and forms the basis. I'd like to avoid using basis in function names since it took me over a month of working in detail with this code to actually understand what the basis we were forming was. I think a more appropriate name for `ode_step` is something like `integrate_el_region` or `integrate_between_singular` since it doesn't integrate the entire EL equation in one go unless there are zero rationals - I support adding in a `set_el_integration_bounds` function. There is a ton of reused logic that occurs in `ode_axis_init`, `ode_sing_init`, `ode_kin_cross`, and `ode_ideal_cross`, that I very much believe can be turned into one function with control logic for `kin_flag` and `init=true/false`. (i.e. the compute conditions at next singular surface logic). This function will find the next singular surface `ising` and set `psimax`, the next integration upper bound, and `next` which should be renamed but denotes if we've finished crossing or not. Perhaps a bool `need_to_cross_rational`? Then with these changes, the logic becomes the below segment. Note that this also assumes the logic was taken out of `ode_axis/sing_init`, so that these functions just initialize the starting `psifac`, `ising`, and `u`, which I think is appropriate for an init function. Depending on how much code it ends up being, I could see combining `set_el_integration_bounds` with `ode_step` (or its new name), since #82 made `ode_step` extremely short now, so I don't think this logic would clutter it and make it "choose integration bounds and integration" all in one function ```julia # Always integrate once, even if no rational surfaces are crossed set_el_integration_bounds() ode_step!(odet, ctrl, equil, ffit, intr)

If at a rational surface, do the appropriate crossing routine, then integrate again

while odet.ising != ctrl.ksing && odet.need_to_cross_rational == "true"
if ctrl.kin_flag
error("kin_flag = true not implemented yet!")
else
ode_ideal_cross!(odet, ctrl, equil, ffit, intr)
end

set_el_integration_bounds()
ode_step!(odet, ctrl, equil, ffit, intr)

end

  - Note that I've simplified/clarified a lot of this next surface logic in my WIP PR OpenFUSIONToolkit/JPEC#75, probably most importantly deprecating `m1` and `singfac` which further simplifies the integration bounds logic. This can be included in this PR if I have not finished multi-n...

</details>



<!-- START COPILOT CODING AGENT SUFFIX -->

- Fixes OpenFUSIONToolkit/JPEC#87

<!-- START COPILOT CODING AGENT TIPS -->
---

✨ Let Copilot coding agent [set things up for you](https://github.com/OpenFUSIONToolkit/JPEC/issues/new?title=✨+Set+up+Copilot+instructions&body=Configure%20instructions%20for%20this%20repository%20as%20documented%20in%20%5BBest%20practices%20for%20Copilot%20coding%20agent%20in%20your%20repository%5D%28https://gh.io/copilot-coding-agent-tips%29%2E%0A%0A%3COnboard%20this%20repo%3E&assignees=copilot) — coding agent works faster and does higher quality work when set up for your repo.

[jhalpern30]

@copilot Please address the comments in the github issue relating to precomputing all chunks before integrating and getting rid of the while loop. Specifically:
"Adding in a set_el_integration_bounds function → sounds like we are in agreement here! Good! I vote it stay separate from integrate_el_region since, as mentioned in 2, the logic for chunking the integration could get more complex and i'd like to be able to not find it internally within the actual integration of said chunks. For me personally, while loops are aways tough to read/understand and bug prone. It seems legacy from the one-step-at-a-time logic. With our modern approach, I think it's cleanest if we just chunk up the integration region and then do a for loop through the chunks. What do you think @jhalpern30 ?
Sudo code replacement of the code block you had:

chunks = chunk_el_integration_bounds()

for chunk in chunks:
if chunk.init_type == 'ideal_crossing'
ode_ideal_cross!(odet, ctrl, equil, ffit, intr). # no need for lugging around all the sing structs here? Just do msing = int(n*q(odet.psifac)) inside this function and we know we are jumping from odet.psifac to chunk.psimin
end
integrate_el_region!(odet, ctrl, equil, ffit, intr, chunk.psistart, chunk.psiend)
end
"

[Copilot]

@copilot Please address the comments in the github issue relating to precomputing all chunks before integrating and getting rid of the while loop. Specifically:
"Adding in a set_el_integration_bounds...

Addressed in commit 278e21b.

The while loop has been replaced with a for loop that iterates over pre-computed integration chunks. I've added:

• IntegrationChunk struct to hold chunk metadata (psi_start, psi_end, needs_crossing, ising)
• chunk_el_integration_bounds() function that pre-computes all chunks upfront
• For loop in eulerlagrange_integration() that iterates through chunks

The logic now follows the pattern you suggested:

chunks = chunk_el_integration_bounds(odet, ctrl, intr)

for chunk in chunks
    odet.ising = chunk.ising  # Track which surface we're near for tolerance calculations
    integrate_el_region!(odet, ctrl, equil, ffit, intr, chunk.psi_start, chunk.psi_end)
    
    if chunk.needs_crossing
        cross_ideal_singular_surf!(odet, ctrl, equil, ffit, intr)
    end
end

All tests passing (ODE: 73/73, Full runs: 2/2).

[jhalpern30]

@copilot Now, please address the last part of the discussion relating to the singular asymptotics:
Feel free to ignore the comment for now, since things are already working.
My instinct is that it's clunky to be doing "where are all the singularities, lets get all the asymptotics there, store them... ok now I am at a crossing! which singularity am I at? I'll look up the pre-calculated asymptotics there". We could easily just chunk to known psi values near singularities, then in our init we can just say "which singularity am I near (m=n*q). Let me calculate the asymptotics here."

In sudo code, my point is that I prefer

for i in range(9)
x = i**2
y = 2*x
end
over

x = zeros(9)
for i in range(9)
x[i] = i**2
end

various other code here that makes me forget how x was formed

for j in range(9)
y = 2*x[j]
end

Agreed that this would be helpful and would come later. For my own reference or whoeever ends up adding this, I think the code would look something like

function ode_ideal_cross!(odet::OdeState, ctrl::DconControl, equil::Equilibrium.PlasmaEquilibrium, ffit::FourFitVars, intr::DconInternal)

Fixup solution at singular surface

ode_unorm!(odet.u, odet, ctrl, intr, true)

Compute asymptotic power series terms alpha/v^k in eq. 45 of Glasser 2016

sing_vmat!() # rename this function to something better

OR sing_asymptotics = sing_vmat(), return a new struct with the relevant asympotic data so we don't have to hide it within the sing stuff, I like this better

Compute full asymptotic series expansion ua at current psi (V in eq. 50) and invert to obtain coefficients c_a just left of rational surface

odet.ca_l[:, :, :, odet.ising] .= sing_get_ca(ctrl, odet, sing_asymptotics)

some other computations that use the asymptotics we just computed

end # asymptotic struct falls out of scope at this function, recompute at next singular surface
I think this also makes it very clear that the asymptotics are only computed and used in ideal DCON and not kinetic DCON (which I think is true?). This would also deprecate sing_scan, which is the function responsible for doing the "compute all x and use it later" logic Nik brings up above.

Please do this by separating the singular surface struct into the existing paramters and the asymptotics that are computed at each singular surface crossing

[Copilot]

@copilot Now, please address the last part of the discussion relating to the singular asymptotics:
Feel free to ignore the comment for now, since things are already working.
My instinct is that it's ...

Addressed in commit 98a6d5d.

I've separated the singular surface data into two structs:

• SingType - Basic parameters (psifac, m, n, q, q1, di)
• SingAsymptotics - Asymptotic data (alpha, vmat, mmat, power, r1, r2, n1, n2, m0mat)

Asymptotics are now computed on-demand in cross_ideal_singular_surf! via compute_sing_asymptotics(), which returns a SingAsymptotics struct. This struct is used immediately and falls out of scope after the crossing function completes.

Key changes:

• Removed sing_scan!() (no longer needed)
• compute_sing_asymptotics() replaces sing_vmat!() and returns struct instead of mutating
• sing_get_ua() and sing_get_ca() now take SingAsymptotics + SingType separately

This makes it crystal clear that asymptotics are only computed/used in ideal DCON. The pattern is now: compute → use → discard, matching your preferred style.

All tests passing (ODE: 73/73, Full runs: 2/2).

[jhalpern30]

@claude please review this PR. Please look for and provide suggestions for any possible bugs that are introduced, any sections where code can be made more concise, improvements to documentation or unit testing. Please provide a concise summary of the entire PR including the modifications and your recommendations

[claude]

Claude encountered an error —— View job

---

I'll analyze this and get back to you.

[jhalpern30]

@copilot please review this PR. Please look for and provide suggestions for any possible bugs that are introduced, any sections where code can be made more concise, improvements to documentation or unit testing. Please provide a concise summary of the entire PR including the modifications and your recommendations

[Copilot]

@copilot please review this PR. Please look for and provide suggestions for any possible bugs that are introduced, any sections where code can be made more concise, improvements to documentation or un...

Code Review Summary

Overview

This PR successfully refactors the Euler-Lagrange integration to make bounds explicit, replace while-loop logic with for-loops over pre-computed chunks, and separate singular surface data from asymptotic calculations. All tests pass (ODE: 73/73, Full runs: 2/2).

Key Modifications

1. File rename: Ode.jl → EulerLagrange.jl
2. Function renames: ode_run → eulerlagrange_integration, ode_step → integrate_el_region, etc.
3. Integration chunking: New chunk_el_integration_bounds() pre-computes all regions upfront
4. Data separation: Split SingType (basic params) from SingAsymptotics (computed on-demand)
5. Simplified sing_get_ua/ca: Now take z::Float64 directly instead of full state objects

Potential Issues & Recommendations

1. Edge Case in chunk_el_integration_bounds (line 52)

while !(ising_current > intr.msing || intr.psilim < intr.sing[ising_current].psifac || ctrl.singfac_min == 0)

• Issue: If ising_current > intr.msing, accessing intr.sing[ising_current] will fail
• Fix: Check bounds first: while ising_current <= intr.msing && intr.psilim >= intr.sing[ising_current].psifac && ctrl.singfac_min != 0

2. Duplicated Logic (lines 36-47 and 69-81)
The surface-skipping logic appears twice. Consider extracting to helper function:

function find_next_resonant_surface!(ising::Int, intr::DconInternal)
    while ising <= intr.msing
        if intr.psilim < intr.sing[ising].psifac || 
           any(m -> intr.mlow <= m <= intr.mhigh, intr.sing[ising].m)
            break
        end
        ising += 1
    end
    return ising
end

3. SingAsymptotics Default Values
Default-initialized arrays (line 47-56 in DconStructs.jl) waste memory. Consider using constructor instead of @kwdef.

4. Documentation

• Line 190 in EulerLagrange.jl references deleted set_el_integration_bounds!
• sing_get_ca docstring (line 621) still references removed parameters

5. Unit Testing
Consider adding tests for:

• chunk_el_integration_bounds with zero/one/multiple singular surfaces
• Edge case: ctrl.singfac_min == 0
• sing_get_ua with positive/negative z values

6. Code Conciseness
Line 283 in EulerLagrange.jl:

ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert

This computation appears in multiple places. Consider adding to SingAsymptotics struct.

Minor Issues

• Line 143: Progress output should be conditional on ctrl.verbose
• Line 272: dpsi computed but could be moved closer to first use (line 275)
• Several TODOs for kin_flag support - ensure consistent across all locations

Strengths

• Clear separation of concerns
• Explicit integration bounds
• On-demand asymptotic computation
• Better parallelization potential
• Good test coverage

[jhalpern30]

@copilot Please implement the unit test suggestions you have above

[jhalpern30]

@logan-nc I think this is ready to review, AI didn't do as well as I'd like so it took longer than I thought but I think I knocked off all the TODOs and this will close #87. A little larger than expected but the main points are the function renaming, moving integration bounds to be precomputing in the chunk_el_bounds function, and restructuring the asymptotics logic to do everything at the surface crossing instead of precomputing everything

[logan-nc]

skimming it, this looks great!
I'll need some more time to read through in detail. In the meantime, can you post some screenshots documenting that the results are unchanged?

[jhalpern30]

I'll need some more time to read through in detail. In the meantime, can you post some screenshots documenting that the results are unchanged?

@logan-nc Comparison of outputs from DIIID-like_ideal_example using BenchmarkTools - identical outputs and no significant change in runtime or allocations, as expected since this is mainly a reorg for clarity.
-- Develop branch

-- This PR

Since I also changed the asymptotics, here is a demonstration that the coefficient matrices just inside and outside the rational surfaces (i.e. ca_right and ca_left in the HDF5 output) are identical between develop and this PR too

[claude]

PR Review: Clarify Euler-Lagrange integration bounds and refactor for explicitness

This PR successfully addresses issue #87 by clarifying integration bounds and separating concerns in the DCON stability analysis code. The changes improve code readability and maintainability.

Strengths

Integration bounds clarity: The pre-computation of integration chunks via chunk_el_integration_bounds() and the IntegrationChunk struct makes the integration flow explicit and easier to understand. The for-loop over chunks replaces the previous while-loop logic, which is more predictable.

Separation of concerns: Splitting SingType into basic parameters and SingAsymptotics for computation-specific data is a good architectural decision. Computing asymptotics on-demand in cross_ideal_singular_surf! clarifies that these are specific to ideal DCON, not inherent to singular surfaces.

Naming improvements: Renaming from ode_* to eulerlagrange_* and Ode.jl to EulerLagrange.jl improves clarity about what these functions do. The file-level rename is consistent with the project's physics-oriented naming.

Test coverage: The addition of integration chunk tests in runtests_eulerlagrange.jl provides good coverage of the new chunking logic.

Issues and Recommendations

1. Documentation inconsistencies (Sing.jl:660-697)

The docstring for sing_der! at line 660 has outdated parameter documentation:

• States params is a 5-tuple but it's actually a 6-tuple (includes IntegrationChunk)
• Line 704 comment says "note the two items not used here" but should clarify which two items

Fix: Update docstring to match actual signature and clarify which params are used where.

2. Function signature comment mismatch (Sing.jl:704)

# note the two items not used here are needed in the integrator params for use in the integrator_callbackcallback

Has a typo: "integrator_callbackcallback" should be "integrator callback"

3. Potential performance consideration (Sing.jl:258)

compute_sing_asymptotics() allocates large arrays (vmat, mmat) on every singular surface crossing. For runs with many singular surfaces, this could impact performance. Consider:

• Profiling memory allocation impact
• If significant, pre-allocate and reuse buffers
• Document if this is acceptable trade-off for clarity

4. Missing validation in chunk_el_integration_bounds (EulerLagrange.jl:179-234)

The function doesn't validate that chunks are monotonically increasing in psi or that they don't overlap. While the algorithm should ensure this, defensive checks would catch bugs earlier.

Suggestion: Add assertions like:

@assert chunk.psi_start < chunk.psi_end "Invalid chunk bounds"
@assert isempty(chunks) || chunk.psi_start >= chunks[end].psi_end "Overlapping chunks"

5. Incomplete TODO removal (EulerLagrange.jl:131-142)

The TODO comment about moving ising_start logic to chunk_el_integration_bounds should be addressed or clarified. The current placement in initialize_el_at_axis! seems reasonable since it depends on the starting psifac, but the comment suggests uncertainty.

6. Code comment violates project standards (multiple locations)

The CLAUDE.md guidelines explicitly state: "No step numbering in code comments - Avoid annotations like 'Step 1: do this' followed by 'Step 2: do that'."

However, I don't see violations of this in the changed code - this is good adherence to the guidelines.

7. sing_order field duplication (DconStructs.jl:43)

SingAsymptotics stores sing_order which is already in ctrl. While this makes the struct self-contained, it creates potential for inconsistency. Consider whether this field is necessary or if functions using SingAsymptotics should also receive ctrl.

Minor Issues

• EulerLagrange.jl:76: Comment says "step was incremented one extra time in ode_step!" but should reference integrate_el_region! (old function name)
• Sing.jl:625-637: Docstring for sing_get_ca has parameter list that doesn't match the actual signature (lists 5 params, function takes 3)

Testing

The PR description indicates all tests pass. Test additions for chunking logic are appropriate. Consider adding a test that verifies asymptotics computed on-demand match pre-computed values (for regression testing).

Performance

No obvious performance regressions expected. The on-demand asymptotic computation replaces pre-computation, which should be roughly equivalent in total cost but more memory-efficient.

Architecture

The refactoring supports future work:

• Pre-computed chunks enable STRIDE-style parallelization
• Explicit bounds make debugging easier
• Separated asymptotics facilitate kinetic DCON implementation

Recommendation

Approve with minor changes. The core refactoring is sound and achieves the stated goals. Address the documentation inconsistencies and consider the performance/validation suggestions before merging.