Loop boundary conditions from Eq. (37) through Riccati.jl#321
Draft
viaweber98 wants to merge 183 commits into
Draft
Loop boundary conditions from Eq. (37) through Riccati.jl#321viaweber98 wants to merge 183 commits into
viaweber98 wants to merge 183 commits into
Conversation
) Store the full merged TOML, raw equilibrium arrays, and (when enabled) forcing-mode data inside gpec.h5 so a run can be replayed from any directory without the original g-file / CHEASE / auxiliary TOML / forcing inputs. `main("run.h5")` dispatches to `main_from_h5`, which rebuilds the `DirectRunInput` / `InverseRunInput` / `SolovevConfig` / `LargeAspectRatioConfig` directly from the snapshot and hands it to `setup_equilibrium` via the existing `additional_input` channel. The rerun output is written as `<basename>_rerun.h5` to guarantee the source is never overwritten; `--output-dir`, `--output-name`, `--override-file`, and repeatable `--override key=value` CLI flags let users retarget or tweak TOML fields without editing the snapshot. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…te galerkin solver to support coil-response columns
…ure/rdcon-galerkin
Add a Julia-specific "Sign Conventions Reference" page to the docs site, mirroring the Fortran GPEC sign-conventions reference but written against the Julia code with traceable src/ references. Registered in the nav after Workflow. Covers the coordinate system and Fourier kernel, poloidal flux, poloidal/toroidal angles, working-coordinate Jacobian, helicity, F=R*Bphi, safety factor q, mode numbers, output sign flips (SURFMN/VACUUM), KineticForces rotation conventions, and COCOS compatibility. Two corrections relative to the Fortran reference, verified against the Julia code and Sauter & Medvedev (2013): - q is computed by field-line integration for direct equilibria and is positive in standard operation, not read directly from the EQDSK. - The internal "COCOS 2" label is not a certified index; document the practical IMAS imas_cocos input requirement and flag full COCOS certification as future work. The KineticForces rotation section was verified against PR #112, correcting the effective bounce harmonic (l_eff = l + nq passing / l trapped, not l - sigma*nq) and the bounce/precession frequency description. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Delete the empty tracked .codex scratch file (and gitignore it) and the Bal_salpha_delta_di_summary.ipynb notebook; its plotting is ported to analyze_example.jl. Notebooks are no longer maintained in the repo. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… annotate IMAS toml Finish the mer_flag/bal_flag -> local_stability_flag rename in the example and regression TOMLs that still used the old flags after the develop merge. Revert psilow 1e-2 -> 1e-4 in the DIIID and Solovev examples (the new ballooning ODE converges fine at 1e-4, verified by running the example). Add per-parameter annotations to the IMAS example to match the other examples. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…iguate ballooning Delta' Restore the resistive interchange D_R (locstab/dr), which became empty when the standalone Mercier.jl was folded into the rewritten Bal.jl. resistive_interchange() recomputes D_R = D_I + (H - 1/2)^2 from flux-surface averages (the former mercier_scan! formulation) and fills locstab column 2; write_outputs_to_HDF5 emits it under local_stability_flag like locstab/di. Rename the ballooning output locstab/delta_prime -> locstab/ballooning_Delta_prime to avoid colliding with the Riccati tearing Delta' (singular_coupling/delta_prime), and drop the always-empty locstab/ca1. NOTE for reviewer/author: locstab/di reports det(d0bar) from the new approach, while D_R is built on the surface-average D_I (old Mercier route). The numerical equivalence of those two D_I values should be confirmed against the Bal derivation. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…alyze_example.jl Add a local-stability section to the DIIID analyze script that reproduces the deleted notebook: s and alpha profiles, Mercier D_I and ballooning Delta' profiles, and 2D s-alpha (p', q') maps of Delta' and D_I with their zero contours, via salpha_reference, compute_ballooning_stability! and scan_delta_prime_map. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ooning/s-alpha) section Document the local_stability_flag scan and the locstab/ HDF5 keys (di, dr, ballooning_Delta_prime), the s-alpha diagram, and the distinction between the local ballooning Delta' and the resistive tearing Delta'. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… Conventions Reference Folds the PerturbedEquilibrium unit-convention documentation (previously the abandoned PR #261) into the conventions reference page and expands its scope: - Rename sign_conventions.md -> conventions.md; retitle "Conventions Reference" and update the nav label, since the page now covers sign, coordinate, and field-amplitude conventions. - Add a "Field Amplitudes and Units" section clarifying the two amplitudes users most often confuse: the resonant field (resonant flux / scalar surface area) and the power-normalized field (sqrt-area-weighted normal field / scalar sqrt-area), both in tesla by scalar-area normalization, with a compact units table for the reported quantities. - Add Pharr (2026), "Coordinate-invariant flux-surface Fourier analysis in tokamaks," to docs/resources and the Citations page as the reference for the coordinate-invariant area weightings. Documentation only; no source or calculations changed. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…ure/rerun-from-h5
#229) Adds modern coil-geometry inputs alongside the legacy ASCII .dat format: - Analytic generators make_pf_hoop, make_window_pane (axis-aligned, two [R,Z] corners), and make_helical (circular-torus), each producing a CoilSet that flows through the existing currents + apply_transforms (shift/tilt) path. - HDF5 coil I/O: group-based, multi-set, content-tolerant (save_coils_to_h5 / load_coils_from_h5_group! / read_coil_h5 / write_coil_h5) plus dat<->h5 converters. One subgroup per set; unrelated content ignored. - CoilSetConfig gains a `source` discriminator (file/pf_hoop/window_pane/helical) and analytic params; load_coil_sets dispatches via _build_raw_coil_set with .h5/.dat resolved by extension. Legacy .dat behavior unchanged. - Rerun snapshot: a coil run stores the geometry used under input/raw_inputs/coils; main_from_h5 gains --coil-source {forcing-modes|coils} so a run can be replayed from gpec.h5 with no original .dat (recomputing the field from stored geometry, optionally against a new equilibrium). - Fixes a pre-existing stale-key bug in runtests_rerun_from_h5.jl (vacuum/et -> FreeBoundaryStability/eigenmode_energies). - Tests for generators, h5 round-trip/multi-set/junk-tolerance, dispatch, and an end-to-end coil snapshot+replay that deletes the .dat between runs. - Docs: forcing_terms.md coil-source + HDF5 schema sections; coil_geometries README HDF5 schema and corrected toroidal-rotation note. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- CoilSetConfig docstring: `[rad] (...)` parsed as a markdown link; reword to plain prose so makedocs cross-reference check passes. - forcing_terms.md / coil_geometries README: label the HDF5 schema fence as `text` and drop bracket/paren unit annotations. Verified: docs build clean (no missing_docs / invalid-link); regression harness diiid_n1 base vs HEAD shows 30/30 quantities unchanged (.dat coil path identical). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…rices to coordinate-invariant field space Implements the first slice of issue #233: form the PerturbedEquilibrium control-surface response matrices and resonant coupling matrices in the coordinate-invariant power-normalized field (b-tilde) space of Pharr (2026), extending the convention PR #225 introduced for the DCON energy eigenvalues. - Extract the shared power-norm transform (compute_sqrtamat, compute_sqrt_jac_delpsi, new control_surface_ptof = sqrtamat*sqrt(jarea)) from ForceFreeStates/PowerNorm.jl into Equilibrium/CoordinateInvariant.jl as one source of truth for FFS and PE. - Conform the control-surface matrices in place: P~=ptof^-1*P*ptof (similarity), L~/Lambda~=ptof^-1*(.)*ptof^-dagger and rho~=ptof^dagger*rho*ptof (congruence). Internal physics (singular coupling, field reconstruction, energies) stays in flux space, so all physical quantities are bit-identical. - Expose ptof in HDF5 (response_matrices/powernorm_field_to_flux_operator) with documented recovery formulas so users can recover the flux-space matrices. - Conform the resonant coupling matrices on input (C~=C*ptof); evaluate the applied resonant scalars from flux-space C first so they carry no round-trip noise. Add LHS coordinate-invariance audit annotations (resonant field, penetrated field, Delta', shielding current). - Rename tesla quantities flux->field per Pharr 2026 to remove units confusion: resonant_flux->resonant_field, C_resonant_flux->C_resonant_field, plot_resonant_flux->plot_resonant_field_amplitude. Genuine flux quantities (psi, Phi_x, chi1) are untouched. Regression case and Fortran benchmark updated in lockstep. Verification: physics review PASS (all transform rules exact and consistent); new unit test runtests_coordinate_invariant.jl 9/9 (recovery round-trip, consistency, energy invariance); regression diiid_n1 develop vs local 28/28 quantities 0.0e+00 (the 2 N/A rows are the intended resonant_flux->resonant_field key rename). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Unify main() so each input source builds a ready (inputs, eq_config, additional_input) and hands it to the renamed main_from_inputs: build_inputs_from_toml and build_inputs_from_h5. This also fixes a latent IMAS bug where main(; dd) never threaded dd into the pipeline body. - equil.toml snapshots are now rerunnable: serialize the built eq_config into inputs[Equilibrium] so gpec_toml_raw is complete (equilibrium_config_to_dict). - parse_override_flag warns instead of silently storing an unparseable value as a String. - Rerun.jl file-level docstring -> comments. - Slim rerun tests to one source + one replay: analytic Solovev replay now diffs all output datasets; direct (EFIT) recovery is checked through the real h5 reader without a pipeline run; override flags and parse logic are covered pipeline-free. Drop the heavy regression-harness rerun cases. - Add an IMAS dd-dispatch regression test. - CLAUDE.md: forbid PR/issue references in source comments. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…control_3d analog) Adds two reusable Analysis.CoilForcing entry points that replicate GPEC's plot_control_3d: - plot_surface_3d(equil; color_by, coil_sets, ...) — the 3D plasma boundary surface, optionally coloured by ANY control-surface output and/or with coils overlaid. One function covers surface-only / +field / +coils / +field+coils; coils-only stays plot_coil_geometry_3d. - control_surface_scalar(equil, modes, n; ...) — inverse-Fourier-transforms spectral modes (b_n, b_n_x, xi_n, …) to a real-space [θ,φ] field, reusing FourierTransforms.inverse and the ν/helicity correction from PerturbedEquilibriumModes (no field recomputation, no Biot-Savart). Conventions: - Shared :RdBu symmetric scale: a blue (+current) coil reads the same colour as the blue (+) field it drives. Sign/handedness chosen and VERIFIED against coil currents so a positive coil produces positive b_n_x directly beneath it (the legacy Python mis-tracks toroidally). - No new Project.toml deps. Coloured 3D surfaces need a PlotlyJS-style backend (already a dep); on GR the tool warns and degrades to an uncoloured surface rather than mis-colouring by Z. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This contains good structure: calculate in Bal.jl, output in h5, and helper function plots in Analysis + a use case example in examples. But details need to be improved: second stability boundary inclusion, plot cleanup to show the window properly, etc.
…uation for big and small solutions, reducing memory allocation and improving performance by 3x
…ile to improve performance and reduce memory allocation by 5x
…E integration by reducing memory allocations and reusing parameters
…ing ODE integration by using SVector and minimizing interpolated coefficients
…previous transcription bug in GW2020 Eq. 55
…ree variants for improved performance in resonant quadrature
…flatten benchmark dir
…llooning stability computation to enhance performance
…pha boundary calculations
…erimental function solve_inner_converged
…-ode-allocation
Galerkin Integrator
Rename the custom Jacobian selector jac_type = "other" to "custom" and add user-facing jac_custom_power_bp/b/r/rc inputs to EquilibriumConfig, read only when jac_type = "custom". Named types (hamada/pest/boozer/equal_arc/park) keep forcing the internal power_* exponents; downstream Jacobian formulas are unchanged and still read the internal power_* fields as the single source of truth. The four power_* constructor slots are now derived-only (marked `_`). Deprecate power_bp/power_b/power_r/power_rc as [Equilibrium] TOML input keys by reusing the centralized deprecation mechanism: generalize _drop_deprecated_ffs_keys! to _drop_deprecated_keys!(table, keys, section) and add _DEPRECATED_EQUIL_KEYS, applied on both the gpec.toml and h5-replay input paths. A stale power_* key now warns loudly and is dropped instead of silently overriding the coordinate system. Also fix the _DEPRECATED_FFS_KEYS one-tuple bug (was a string, so the FFS deprecation never fired). Update the EquilibriumConfig docstring and conventions.md, and clean the now-inactive power_* lines out of the example/test gpec.toml files while standardizing the jac_type option-list comment. Regression: diiid_n1, solovev_n1, solovev_multi_n all 0.0e+00 diff vs develop. Closes #288 Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…er-knobs EQUIL - IMPROVEMENT - Scope custom-Jacobian knobs to jac_custom_power_*
…tic-profiles # Conflicts: # examples/DIIID-like_ideal_example/gpec.toml
…profiles Add D3D-like kinetic profiles and update example equilibrium
…rams to reduce multithreading overhead and add a pre-allocated quadgk buffer
…memory allocation in QuadGK inner loop and update Energy Integration to support QuadGK buffer pre-allocation
…-surface scratch buffers to reduce memory allocation during bounce-averaging computations
… into gpec.toml and remove deprecated sol.toml
…tialized entries in bounce integration
…buffer, BounceScratch field docs fvals stays ComplexF64: the kinetic-matrix (wt) path builds fbnce from complex op_wmats data, so a real buffer throws InexactError in the pitch kernel. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TTu6sczod1Qtjy9e9pB6jP
…eallocation Pentrc statictyping preallocation
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
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.
The intent of this branch is to add a loop that solves the Δ' boundary-value problem from Eq. (37) in the paper "A robust..." (Glasser & Koleman 2018). The idea from the paper states that you build the big BVP matrix once, then solve it once for each boundary condition while changing the RHS each time. For N singular surfaces there are 2N boundary conditions, and each solve gives back the small-solution coefficients (which fill in the D' matrix).
The plan is to assemble the Eq. (37) matrix once, factorize it, loop all over the 2N boundary conditions (setting the big-solution coefficient to 1, and everything else to 0), solve the system for each boundary condition, and collect the small-solution coefficients into the D' matrix.