ForceFreeStates: physics-driven auto-computation of mlow/mhigh and mthvac - #193
ForceFreeStates: physics-driven auto-computation of mlow/mhigh and mthvac#193logan-nc wants to merge 6 commits into
Conversation
…ow/mhigh and mthvac Replace fixed delta_mlow/delta_mhigh padding and manual mthvac with physics-driven auto-computation: - _estimate_m_bandwidth: DFT of Jacobian at psilim gives metric Fourier bandwidth; used as poloidal mode padding instead of hardcoded delta_m values - _estimate_arc_nonuniformity: mean/min arc-length ratio at psilim detects SFL coordinate clustering (e.g. near X-point in Hamada coords); scales mthvac to compensate for reduced vacuum integration accuracy - Auto-mthvac: max(4*mpert, max(96, 6*mhigh)*arc_factor) rounded to next multiple of 32 - mlow/mhigh: explicit user overrides (nonzero bypasses auto) - m_accuracy: single tuning knob for mode range (replaces delta_mlow/delta_mhigh) - Remove ctrl.delta_mhigh *= 2 legacy hack - Deprecate delta_mlow/delta_mhigh with @warn; fields retained for backward compat Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
PlotlyJS was added to Project.toml in the develop merge but is not imported anywhere in src/ or test/. Its presence breaks Pkg.instantiate() due to a JSON version conflict (PlotlyJS 0.18.17 requires JSON < 0.22, incompatible with JSON 1.4.0 required elsewhere). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…uto-computed) equil.config.mpsi stays 0 when auto-mpsi is used (make_optimal_mpsi returns the value but never writes it back to config). sing_lim! used equil.config.mpsi - 1 as a clamp, giving jpsi = -1 and a BoundsError. Use length(profiles.xs) - 1 instead, which is always correct regardless of how mpsi was determined. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…dwidth estimation
In Hamada coordinates the Jacobian is constant by construction (equal-volume cells),
so _estimate_m_bandwidth returned 0, collapsing mlow to 0 and giving completely wrong
et[1] values for DIIID-like H-mode equilibria.
r² (rzphi_rsquared) is the correct proxy: it is proportional to g^{θθ} and carries the
angular variation of the metric that drives inter-mode coupling in the EL equations,
regardless of which SFL coordinate system is used.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…gh/mlow floors Generalize _estimate_m_bandwidth to accept any 2D field callable, normalizing by peak amplitude (rather than DC) so zero-mean fields (e.g. metric cross-terms) are handled correctly. Add _estimate_metric_bandwidth: scans all 5 equilibrium metric fields (rzphi_rsquared, rzphi_jac, eqfun_B, eqfun_metric1, eqfun_metric2), takes the maximum single-field DFT bandwidth, then multiplies by 4 to account for convolution of metric products in the EL equations. Add physics-motivated floors: - mhigh ≥ 2*nhigh*qmax: mode coupling in perturbed equilibria extends to ~2× the highest resonant surface (Park et al. rule of thumb) - mlow ≤ -0.5*mhigh: maintain ~2:1 positive/negative m asymmetry Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
|
@logan-nc is this PR going anywhere? |
|
@matt-pharr I believe this was put on the back burner since until things "settled down" and we became more ok with Julia diverging from fortran. If implemented then it will become harder to directly compare to fortran (still possible - you would just run fortran, read the m range, then set julia mlow and mhigh based on that). If you and @jhalpern30 think it's time - I can revisit this to bring it up to date and clean it up a bit. Rereading the description, one thing that would be a good add is to do the ForcingTerms calculation earlier in the workflow and include that knowledge when setting m ranges - a classic issue is that the plasma physics (DCON) sets some optimized m-range but the simple boxy 3D coils have some broad -m spread and the m-range has to be expanded to fully capture the forcing terms (otherwise weird aliasing can happen). That's probably the number one failure mode that makes people have to re-run cases with more delta_m*, so we should automate the logic. It would be good to test on some more cases too. I feel like the AIs always over optimize local conditions... perhaps the shaping scan? That is a script in benchmarks now, right? A q95 scan would be good too. |
|
@logan-nc the scan tests will be much easier after I finish the tokamaker integration but that is still a ways out, will likely not want to finish that until we make our big refactor push |
|
I think this is useful at some point but probably a distraction on the path to GPEC 2.0 and can wait |
@matt-pharr Don't we have a aspect ratio scan in a benchmark script from @d-burg? Is that on develop now? I guess it uses TJ, so can do q95 scans and aspect ratio but not shaping? Solovev can do some shaping I suppose. |
3 participants
Summary
Test results
The DIIID `et[1]` values differ from `develop` even when mlow/mhigh are set identically to develop's values. Isolation testing confirmed this is not caused by the auto-m changes: setting `mlow=-12, mhigh=21` (identical to develop) on this branch still produces a different `et[1]`.
The source is pre-existing equilibrium differences inherited from `direct_equilibrium_improvements`: the auto-mpsi feature produces only 3 edge grid points at `psilim=0.989` (vs develop's 128 uniform points), causing under-resolved equilibrium near the separatrix. This is a separate issue to be addressed in `direct_equilibrium_improvements` before merging to develop.
The auto-m and auto-mthvac features are self-consistent: DIIID gives bandwidth=4, arc_factor=4.01 (correctly reflecting Hamada coordinate clustering near the X-point), and Solovev validates cleanly against its own baseline.
Test plan
🤖 Generated with Claude Code