From e33bdde98320ac1312ef15570524f6b5122a3eb8 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 2 Nov 2026 00:14:16 -0500 Subject: [PATCH 01/10] Docs - DOCS - Rebuild the workflow map, condense setup and order the module nav by pipeline Rebuild workflow.md as the physics-first map of every stage main runs, with a which-stage-do-I-need table, headline-output definitions, the complete deck section index, all example directories and the output-file groups. Give index.md a short stage table that points into it. Collapse set_up.md to one install procedure for Linux, macOS and WSL with the correct repository URL and a first run that reads three datasets; move its developer tips to developer_notes.md. Rename the "API Reference" nav group to "Modules" in pipeline order. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- docs/make.jl | 10 +- docs/src/developer_notes.md | 100 ++++++++ docs/src/index.md | 30 ++- docs/src/set_up.md | 468 ++++-------------------------------- docs/src/workflow.md | 309 +++++++++++++----------- 5 files changed, 342 insertions(+), 575 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 8b5fe99dbc..1bd8818aa2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,15 +28,15 @@ makedocs(; "Workflow" => "workflow.md", "Scripting API" => "api.md", "Conventions Reference" => "conventions.md", - "API Reference" => [ - "Vacuum" => "vacuum.md", + "Modules" => [ "Equilibrium" => "equilibrium.md", - "Stability Analysis" => "stability.md", + "Vacuum" => "vacuum.md", + "Stability" => "stability.md", "Galerkin Solver" => "galerkin.md", - "Ballooning Local Stability" => "ballooning.md", - "KineticForces" => "kinetic_forces.md", + "Local Stability" => "ballooning.md", "Forcing Terms" => "forcing_terms.md", "Perturbed Equilibrium" => "perturbed_equilibrium.md", + "Kinetic Forces" => "kinetic_forces.md", "Error Fields" => "error_fields.md", "Tearing" => "inner_layer.md", "Analysis" => "analysis.md", diff --git a/docs/src/developer_notes.md b/docs/src/developer_notes.md index f25c828e1a..b5c55c963b 100644 --- a/docs/src/developer_notes.md +++ b/docs/src/developer_notes.md @@ -239,3 +239,103 @@ regress --cases solovev_n1 --ref-range develop~10..develop The full flag list, the environment-pinning behaviour that keeps source code the only variable in a comparison, and the exit-status contract are documented in [`docs/development/regression-harness.md`](https://github.com/OpenFUSIONToolkit/GPEC/blob/develop/docs/development/regression-harness.md). + +## Developer tips + +### Precompile workload + +Building the package runs a small Solovev case through the solver so that a fresh Julia process +skips most first-call compilation, at the cost of a longer one-time build. Code compiled into the +package image can round differently in the last bit from code compiled on first use, and +near-marginal quantities such as Δ′ amplify that, so results are reproducible only within one +build mode. The regression harness runs with the workload disabled. To disable it yourself (for a +faster build, or to match harness numbers), set the preference by UUID from the project, without +loading GPEC (which would first build it with the workload on). `Preferences` must be installed, +e.g. in your global environment: + +```julia +using Preferences, UUIDs +set_preferences!(UUID("462872dd-e066-4d2e-b993-6468b5239634"), "precompile_workload" => false; force=true) +``` + +This writes `LocalPreferences.toml` next to the active `Project.toml` (the file is gitignored); +delete the entry to re-enable. The regression harness overrides this file for its own runs. + +The workload runs only when Julia rebuilds the package image, not at the start of each run. That +happens after any change under `src/` (a pull, a branch switch or a local edit) and after a change +of dependency versions, Julia version or this preference. Repeated runs on an unchanged checkout +reuse the cached image. Julia keeps a separate image per set of compile flags, so the first +`Pkg.test()` after such a change (which checks array bounds) runs the workload once more. If you edit `src/` and restart Julia often, disabling the workload keeps +each rebuild short; within one Revise session nothing is rebuilt. + +### Revise.jl + +When iterating on code, use [Revise.jl](https://timholy.github.io/Revise.jl/stable/) to avoid full recompilation on every change. It tracks source file modifications and recompiles only the affected code. + +Install Revise in your **global** Julia environment (not in this project's `Project.toml`): + +```julia +# In the Julia REPL (no --project flag) +using Pkg +Pkg.add("Revise") +``` + +Then load it before `GeneralizedPerturbedEquilibrium` in any session: + +```julia +using Revise +using GeneralizedPerturbedEquilibrium +``` + +For it to load automatically, add `using Revise` to your [Julia startup file](https://timholy.github.io/Revise.jl/stable/config/#Using-Revise-by-default): + +```julia +# ~/.julia/config/startup.jl +try + using Revise +catch e + @warn "Could not load Revise" exception=e +end +``` + +### Pre-commit hooks + +The repository uses pre-commit hooks to maintain code quality. They run automatically on `git commit`, and +the same hooks run in CI on every pull request (the `pre-commit` check), which fails if they would change a file. +They: +- Format Julia code according to `.JuliaFormatter.toml` settings +- Remove trailing whitespace and fix line endings +- Validate YAML/TOML syntax +- Prevent accidentally committing large files (>5MB) +- Check config-style TOML annotation conventions (header block, inline `# description` + comments, no deprecated keys — see `docs/development/toml-conventions.md`) + +**Installation:** + +```bash +# Install pre-commit (requires Python/pip) +pip install pre-commit + +# Install JuliaFormatter 1.0.62 in the default environment: the hook runs whatever version is there, +# and CI pins this one, so another version formats differently and fails the check +julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="1.0.62"))' + +# Install the git hooks in your local repository +cd /path/to/GPEC +pre-commit install +``` + +You can also run the hooks manually without committing: + +```bash +pre-commit run --all-files # run on all files +pre-commit run # run on staged files only +``` + +In rare cases where you need to bypass the hooks: + +```bash +git commit --no-verify +``` + +This is discouraged as it may introduce formatting inconsistencies. diff --git a/docs/src/index.md b/docs/src/index.md index 287c9e46e8..f2e2ae295a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -13,18 +13,32 @@ readme = replace(readme, r"^Full documentation:.*$"m => "") Markdown.parse(readme) ``` -## Overview -- [**Introduction Slides**](https://docs.google.com/presentation/d/1Iey7bsvpDG1WDRow24OD0NOWyKYOXMs0U7QLc7Y-wj4/edit?slide=id.g2798d4bd85b_0_15#slide=id.g2798d4bd85b_0_15) - GPEC Introduction slides - +## What GPEC computes + +One deck, `gpec.toml`, runs the stages below in order, each when its section is present, into one output file, `gpec.h5`. +The [Workflow](workflow.md) page maps every stage to its physics, deck sections and datasets. + +| Stage | Deck section | Main output | +|---|---|---| +| `Equilibrium` — the axisymmetric equilibrium, its flux-surface coordinates and rational surfaces | `[Equilibrium]` | `Equilibrium/` | +| `ForceFreeStates` — ideal (or drift-kinetic) MHD stability by the Newcomb criterion, with the `Vacuum` response of the `[Wall]` | `[ForceFreeStates]` | `ForceFreeStates/FreeBoundaryStability/`, `SingularSurfaces/Delta_prime_matrix` | +| `LocalStability` — Mercier, resistive-interchange and ballooning criteria | `[LocalStability]` | `LocalStability/` | +| `Match` — the resistive inner layer matched to the outer-region Δ′ | `[Match]` | `SingularSurfaces/Match/` | +| `PerturbedEquilibrium` — the plasma response to the field of `[ForcingTerms]`: resonant fields, islands, overlap | `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/` | +| `KineticForces` — the neoclassical toroidal viscosity torque of that response | `[KineticForces]` | `KineticForces/` | +| `ErrorFields` — coil tolerances, locking risk and correction limits | `[ErrorFields]` | `ErrorFields/Risk/` | +| `Tearing` — tearing growth rates and critical resonant fields from an inner-layer model | `[Tearing]` | `Tearing/Roots/` | ## Documentation -- [**Setup**](set_up.md) — Installation instructions for macOS and Windows (WSL) -- [**Workflow**](workflow.md) — Full pipeline description: inputs, outputs, and key physics per module -- [**API Reference**](vacuum.md) — Detailed function and type documentation for each module +- [**Setup**](set_up.md) — Installation on Linux, macOS and Windows (WSL), and a first run +- [**Workflow**](workflow.md) — The pipeline: physics of each stage, deck sections, examples and the output file +- [**Scripting API**](api.md) — The same pipeline from Julia scripts +- [**Conventions**](conventions.md) — Coordinates, signs, normalizations and field amplitudes +- [**Modules**](equilibrium.md) — Governing equations, deck keys, outputs and function reference per module - [**Citations**](citations.md) — Papers underpinning GPEC's algorithms -- [**Developer Notes**](developer_notes.md) — Git workflow, commit conventions, and the regression-testing harness - +- [**Developer Notes**](developer_notes.md) — Git workflow, commit conventions and the regression-testing harness +- [**Introduction slides**](https://docs.google.com/presentation/d/1Iey7bsvpDG1WDRow24OD0NOWyKYOXMs0U7QLc7Y-wj4/edit?slide=id.g2798d4bd85b_0_15#slide=id.g2798d4bd85b_0_15) — An overview of GPEC and its physics ## Contact diff --git a/docs/src/set_up.md b/docs/src/set_up.md index c63df74f4b..4ca77d190b 100644 --- a/docs/src/set_up.md +++ b/docs/src/set_up.md @@ -1,460 +1,88 @@ -# Setting up GPEC +# Setup -- [Setting up GPEC](#setting-up-gpec) - - [On Windows via WSL (Ubuntu)](#on-windows-via-wsl-ubuntu) - - [On macOS](#on-macos) - - [Developer Tips](#developer-tips) +GPEC is a Julia package. +Installing it means installing Julia, cloning the repository and letting Julia's package manager fetch the dependencies; the same steps apply on Linux, macOS and Windows (WSL). -## On Windows via WSL (Ubuntu) -1. Install WSL and Ubuntu - If you don't already have WSL and Ubuntu installed, set this up. [This page](https://learn.microsoft.com/en-us/windows/wsl/install) gives detailed instructions on how to complete the installation. In the Windows Powershell, +## Install - 1. Make sure WSL is installed: - ```PowerShell - wsl --install - ``` - 2. Set Ubuntu as your default WSL distro: - ```PowerShell - wsl --set-default Ubuntu - ``` - 3. Launch Ubuntu and update: - ```PowerShell - sudo apt update && sudo apt upgrade -y - ``` -2. Install build tools in WSL - ```shell - sudo apt install build-essential cmake -y - ``` +1. **Julia 1.11 or later.** The simplest route on every platform is [juliaup](https://github.com/JuliaLang/juliaup), which also keeps Julia updated: + - *Linux*: `curl -fsSL https://install.julialang.org | sh`, or unpack a tarball from [julialang.org/downloads](https://julialang.org/downloads/) and add its `bin/` to `PATH`. On a cluster, check `module avail julia` first. + - *macOS*: `brew install juliaup` (or `brew install julia`), or the `.dmg` from the downloads page. Install the Xcode command-line tools (`xcode-select --install`) for `git`. + - *Windows*: install WSL with an Ubuntu distribution (`wsl --install` in PowerShell), then follow the Linux steps inside it. Keep the repository in the WSL file system, not under `/mnt/c`. - `build-essential` → GCC, make + Check with `julia --version`. - `cmake` → sometimes needed by dependencies - -3. Install Julia in WSL - - 1. Download the latest Linux tarball from the official site [Julia downloads](https://julialang.org/downloads/). It will look like - ```shell - wget https://julialang-s3.julialang.org/bin/linux/x64/1.11/julia-1.11.3-linux-x86_64.tar.gz - ``` - - ☆ Replace the URL with the latest stable version. - - 2. Extract and move it to /opt (or any path): - - ```shell - tar -xvzf julia-1.11.3-linux-x86_64.tar.gz - sudo mv julia-1.11.3 /opt/ - ``` - - ☆ Ensure these commands match the tarball you installed. These commands match the above tarball and might need to be modified for your installation. - - 3. Add Julia to PATH: - - ```shell - echo 'export PATH=/opt/julia-1.11.3/bin:$PATH' >> ~/.bashrc - source ~/.bashrc - ``` - - 4. Test it is properly installed - - ```shell - julia --version - ``` - -4. Clone GPEC into your WSL home folder: - ```shell - git clone https://github.com/OpenFUSIONToolkit/GeneralizedPerturbedEquilibrium.git - cd GeneralizedPerturbedEquilibrium - ``` - -5. Install the Julia packages for GPEC - 1. Launch Julia: - ```shell - julia - ``` - 2. In Julia REPL: - ```julia - using Pkg - Pkg.instantiate() # install recorded dependencies - Pkg.add("Preferences") # install missing dependency if needed - Pkg.precompile() # precompile all packages - probably unnecessary - ``` - -6. At this point, you should be able to run the code or connect VS Code to your WSL session. - 1. (Optionally) Integrate WSL with VS Code - 1. Install **Remote - WSL** extension in VS Code. - 2. Open VS Code → Connect To → Connect to WSL. - 3. Click Open Folder and then navigate to the GPEC folder on your VM. Open your GPEC folder from WSL: ~/GeneralizedPerturbedEquilibrium. - - If this is not working, you can launch vscode from the WSL shell you have using the command `code .` - - 4. Open a terminal inside VS Code — it will automatically use WSL/Ubuntu. - 5. You can now run Julia scripts directly from the terminal. - 2. Run GPEC - 1. Launch Julia and run your scripts as usual: - ```shell - include("path/to/script.jl") - ``` - -## On macOS - -### Prerequisites - -Before starting, you'll need a Terminal app to enter commands. You can find it by: -- Press `Cmd + Space` to open Spotlight -- Type "Terminal" and press Enter - -Keep this Terminal window open throughout the installation process. - -### 1. Install Xcode Command Line Tools - -These tools provide the compilers and build utilities needed by GPEC's dependencies. - -1. Open Terminal and run: +2. **Clone the repository.** ```bash - xcode-select --install + git clone https://github.com/OpenFUSIONToolkit/GPEC.git + cd GPEC ``` -2. A dialog will appear asking you to install the tools. Click "Install" and wait for it to complete (this may take several minutes). - -3. Verify installation: +3. **Install the dependencies** into the project environment and check that the package loads: ```bash - gcc --version + julia --project=. -e 'using Pkg; Pkg.instantiate()' + julia --project=. -e 'using GeneralizedPerturbedEquilibrium' ``` - You should see output showing the GCC version. - -### 2. Install Homebrew (Package Manager) - -Homebrew makes it easy to install software on macOS. If you already have Homebrew installed, skip to step 3. + The first call downloads and precompiles every dependency, which takes several minutes; the second builds GPEC itself, which also runs a small precompile workload. -1. Install Homebrew by running this command in Terminal: - ```bash - /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)" - ``` +Every `julia` command on GPEC needs `--project=.` (or `--project=/path/to/GPEC`): without it Julia uses the default environment and cannot find the package. +Never run Julia with `sudo`. -2. Follow the on-screen instructions. You may need to enter your Mac password. +## First run -3. After installation completes, the installer will show you two commands to run to add Homebrew to your PATH. They will look something like: - ```bash - echo 'eval "$(/opt/homebrew/bin/brew shellenv)"' >> ~/.zprofile - eval "$(/opt/homebrew/bin/brew shellenv)" - ``` - **Important:** Copy and run these exact commands from your Terminal output. +Run one of the shipped examples; the deck, the equilibrium file and the coil geometry are all in the directory: -4. Verify Homebrew is installed: - ```bash - brew --version - ``` - -### 3. Install Julia - -Julia is the programming language GPEC is written in. - -**Option A: Install via Homebrew (Recommended for beginners)** - -1. Install Julia: - ```bash - brew install julia - ``` - -2. Verify Julia is installed: - ```bash - julia --version - ``` - You should see something like `julia version 1.11.x`. - -**Option B: Install via Official Installer** - -1. Go to [https://julialang.org/downloads/](https://julialang.org/downloads/) - -2. Download the macOS installer (`.dmg` file) for the latest stable version (1.11 or higher) - -3. Open the downloaded `.dmg` file and drag Julia to your Applications folder - -4. Add Julia to your PATH by running in Terminal: - ```bash - sudo mkdir -p /usr/local/bin - sudo ln -s /Applications/Julia-1.11.app/Contents/Resources/julia/bin/julia /usr/local/bin/julia - ``` - Replace `1.11` with your actual version if different. - -5. Verify Julia is installed: - ```bash - julia --version - ``` - -### 4. Clone the GPEC Repository - -Now we'll download the GPEC code from GitHub. - -1. Choose where you want to put GPEC. For example, your home directory: - ```bash - cd ~ - ``` - Or create a Code folder: - ```bash - mkdir -p ~/Code - cd ~/Code - ``` - -2. Clone GPEC from GitHub: - ```bash - git clone https://github.com/OpenFUSIONToolkit/GeneralizedPerturbedEquilibrium.git - ``` - If you don't have `git` installed, macOS will prompt you to install it. - -3. Enter the GPEC directory: - ```bash - cd GeneralizedPerturbedEquilibrium - ``` - -### 5. Install Julia Packages - -Now we'll install all the Julia packages that GPEC depends on. - -1. Launch Julia from the GPEC directory: - ```bash - julia --project=. - ``` - The `--project=.` flag tells Julia to use the GPEC project environment. - -2. You should now see the Julia prompt: `julia>` - -3. Install all dependencies by typing these commands in the Julia prompt: - ```julia - using Pkg - Pkg.instantiate() - ``` - This will download and install all required packages. It may take several minutes the first time. - -4. Precompile all packages (optional, but speeds up first use): - ```julia - Pkg.precompile() - ``` - -5. Test that GPEC loads correctly: - ```julia - using GeneralizedPerturbedEquilibrium - ``` - If you see no errors, everything is working! - -6. Exit Julia: - ```julia - exit() - ``` - -### 6. Run an Example - -Now you're ready to run GPEC! See the [Running GPEC](#running-gpec) section below for full details. As a quick test: - -1. Make sure you're in the GPEC directory: - ```bash - cd ~/Code/GeneralizedPerturbedEquilibrium - ``` - -2. Run GPEC on one of the included examples: - ```bash - ./gpec examples/DIIID-like_ideal_example - ``` - - The first time you run GPEC, it will take a few minutes to compile. Subsequent runs will be faster. - -### 7. Troubleshooting - -**If Julia can't find packages:** - -Make sure you're running Julia with the project environment: -```bash -cd ~/Code/GeneralizedPerturbedEquilibrium -julia --project=. -``` - -**If you get permission errors:** - -Make sure you have write permissions in the GPEC directory. You may need to use `sudo` for some Homebrew commands, but avoid using `sudo` with Julia commands. - -### Alternative: Using VS Code (Optional) - -VS Code provides a convenient IDE for developing with GPEC: - -1. Install VS Code from [https://code.visualstudio.com/](https://code.visualstudio.com/) - -2. Install the Julia extension: - - Open VS Code - - Click the Extensions icon (or press `Cmd + Shift + X`) - - Search for "Julia" and install the official Julia extension - -3. Open the GPEC folder in VS Code: - - Click File → Open Folder - - Navigate to and select your GPEC directory - -4. Open a terminal inside VS Code (Terminal → New Terminal) and run scripts directly: - ```bash - julia --project=. path/to/script.jl - ``` - -## Running GPEC - -Once GPEC is installed and built, you can run it in two ways: - -### As a Command-Line Script - -GPEC includes an executable script (`gpec`) in the project root directory. To run GPEC on a directory containing a `gpec.toml` configuration file: - -```bash -./gpec path/to/directory -``` - -**Example:** ```bash -# Run GPEC on one of the included examples ./gpec examples/DIIID-like_ideal_example - -# Run in the current directory (must contain gpec.toml) -./gpec ``` -The script will: -1. Read configuration from `gpec.toml` in the specified directory -2. Load or generate the equilibrium based on the `[Equilibrium]` section -3. Compute force-free states (stability analysis) based on the `[ForceFreeStates]` section -4. If a `[PerturbedEquilibrium]` section exists, compute the plasma response to external perturbations -5. Save every stage into the one output file named by `[Run] output_file` (default `gpec.h5`) - -**Early Termination:** a stage runs only when its section is present. A deck without -`[ForceFreeStates]` stops after the equilibrium; one without `[PerturbedEquilibrium]` stops after -the stability analysis. - -### As a Julia Library - -You can also use GPEC programmatically in your own Julia scripts: +The first run compiles the code it needs and takes a few minutes; later runs are faster. +The run logs each stage and writes `examples/DIIID-like_ideal_example/gpec.h5`. +A quick look at three results, from the shell or from Julia: -```julia -using GeneralizedPerturbedEquilibrium - -# Run the full GPEC analysis pipeline -GeneralizedPerturbedEquilibrium.main(["path/to/directory"]) - -# Or access individual modules -using GeneralizedPerturbedEquilibrium.Equilibrium -using GeneralizedPerturbedEquilibrium.Vacuum -using GeneralizedPerturbedEquilibrium.ForceFreeStates - -# Set up equilibrium only -equil = Equilibrium.PlasmaEquilibrium("path/to/g-file"; eq_type="efit") - -# Access equilibrium data -println("q at axis: ", equil.params.q0) -println("Beta-N: ", equil.params.betan) +```bash +h5ls -r examples/DIIID-like_ideal_example/gpec.h5 | less ``` -### Configuration Files - -GPEC uses TOML configuration files (`gpec.toml`) with the following main sections: - -- **`[Equilibrium]`**: Equilibrium solver settings (input file, grid resolution, coordinate system, etc.) -- **`[Wall]`**: Optional wall geometry for vacuum calculations (shape and its parameters); omit for no wall -- **`[ForceFreeStates]`**: Stability analysis settings (mode numbers, tolerances, flags) -- **`[PerturbedEquilibrium]`**: Plasma response settings (forcing data, output options) - -See the example directories for complete configuration file templates. - -## Developer Tips - -### Precompile workload - -Building the package runs a small Solovev case through the solver so that a fresh Julia process -skips most first-call compilation, at the cost of a longer one-time build. Code compiled into the -package image can round differently in the last bit from code compiled on first use, and -near-marginal quantities such as Δ′ amplify that, so results are reproducible only within one -build mode. The regression harness runs with the workload disabled. To disable it yourself (for a -faster build, or to match harness numbers), set the preference by UUID from the project, without -loading GPEC (which would first build it with the workload on). `Preferences` must be installed, -e.g. in your global environment: - ```julia -using Preferences, UUIDs -set_preferences!(UUID("462872dd-e066-4d2e-b993-6468b5239634"), "precompile_workload" => false; force=true) +using HDF5 +h5open("examples/DIIID-like_ideal_example/gpec.h5") do f + read(f, "Equilibrium/q_95") + read(f, "ForceFreeStates/FreeBoundaryStability/eigenmode_energies") # positive first entry: ideal-kink stable + read(f, "PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field") # resonant field per rational surface +end ``` -This writes `LocalPreferences.toml` next to the active `Project.toml` (the file is gitignored); -delete the entry to re-enable. The regression harness overrides this file for its own runs. - -The workload runs only when Julia rebuilds the package image, not at the start of each run. That -happens after any change under `src/` (a pull, a branch switch or a local edit) and after a change -of dependency versions, Julia version or this preference. Repeated runs on an unchanged checkout -reuse the cached image. Julia keeps a separate image per set of compile flags, so the first -`Pkg.test()` after such a change (which checks array bounds) runs the workload once more. If you edit `src/` and restart Julia often, disabling the workload keeps -each rebuild short; within one Revise session nothing is rebuilt. +The [Workflow](workflow.md) page explains what each stage computed and which datasets to read next, and `examples/DIIID-like_ideal_example/analyze_example.jl` plots the output with the [Analysis](analysis.md) module. -### Revise.jl - -When iterating on code, use [Revise.jl](https://timholy.github.io/Revise.jl/stable/) to avoid full recompilation on every change. It tracks source file modifications and recompiles only the affected code. +## Running GPEC -Install Revise in your **global** Julia environment (not in this project's `Project.toml`): +The `gpec` script in the repository root runs the pipeline on a directory holding a `gpec.toml` deck (the current directory when no path is given) and writes the `[Run] output_file` next to it: -```julia -# In the Julia REPL (no --project flag) -using Pkg -Pkg.add("Revise") +```bash +./gpec path/to/directory ``` -Then load it before `GeneralizedPerturbedEquilibrium` in any session: +The same pipeline is available from Julia, either whole or stage by stage: ```julia -using Revise using GeneralizedPerturbedEquilibrium -``` - -For it to load automatically, add `using Revise` to your [Julia startup file](https://timholy.github.io/Revise.jl/stable/config/#Using-Revise-by-default): - -```julia -# ~/.julia/config/startup.jl -try - using Revise -catch e - @warn "Could not load Revise" exception=e -end -``` - -### Pre-commit Hooks -The repository uses pre-commit hooks to maintain code quality. They run automatically on `git commit`, and -the same hooks run in CI on every pull request (the `pre-commit` check), which fails if they would change a file. -They: -- Format Julia code according to `.JuliaFormatter.toml` settings -- Remove trailing whitespace and fix line endings -- Validate YAML/TOML syntax -- Prevent accidentally committing large files (>5MB) -- Check config-style TOML annotation conventions (header block, inline `# description` - comments, no deprecated keys — see `docs/development/toml-conventions.md`) +# The full pipeline of the deck, returning every stage result +results = GeneralizedPerturbedEquilibrium.main(["path/to/directory"]) -**Installation:** - -```bash -# Install pre-commit (requires Python/pip) -pip install pre-commit - -# Install JuliaFormatter 1.0.62 in the default environment: the hook runs whatever version is there, -# and CI pins this one, so another version formats differently and fails the check -julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="1.0.62"))' - -# Install the git hooks in your local repository -cd /path/to/GPEC -pre-commit install +# One stage: the equilibrium of a g-file +using GeneralizedPerturbedEquilibrium.Equilibrium +eq = PlasmaEquilibrium("path/to/g-file"; eq_type="efit") ``` -You can also run the hooks manually without committing: +A stage runs only when its deck section is present; the section index on the [Workflow](workflow.md) page lists them, and the [Scripting API](api.md) page shows the `solve` calls behind each stage. -```bash -pre-commit run --all-files # run on all files -pre-commit run # run on staged files only -``` +## Editor -In rare cases where you need to bypass the hooks: +[VS Code](https://code.visualstudio.com/) with the Julia extension is a convenient environment: open the `GPEC` folder, and the integrated terminal runs `julia --project=.` or `./gpec` as above. +On Windows, use the *Remote - WSL* extension and open the folder from inside WSL. -```bash -git commit --no-verify -``` +## Developing -This is discouraged as it may introduce formatting inconsistencies. +Pre-commit hooks, the precompile workload and Revise are described under *Developer tips* on the [Developer Notes](developer_notes.md) page. diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 02e227150e..d42ec4e062 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -1,195 +1,220 @@ # Workflow -GPEC follows a five-stage analysis pipeline driven by a single `gpec.toml` configuration file. Each stage produces structured output consumed by the next. +GPEC reads one deck, `gpec.toml`, runs every stage whose section the deck contains, and writes every result into one output file, `gpec.h5`. +This page is the map: what each stage computes physically, which deck sections run it, and which datasets to read first. +The module pages give the governing equations and the full key and output tables; the [Scripting API](api.md) page shows the same pipeline called from Julia. ``` - gpec.toml - (user options for every module) - │ - ┌────────────┬───────────┴───┬─────────────────────┐ - ▼ ▼ ▼ ▼ - Equilibrium ──► Vacuum ──► ForceFreeStates ──► PerturbedEquilibrium ──► gpec.h5 - ▲ │ ▲ ▲ - │ └────────────────────────────────────┘ │ - │ │ - equilibrium file ForcingTerms - kinetic profiles (planned) (external perturbation file) + equilibrium file ──► Equilibrium ──► ForceFreeStates ──► PerturbedEquilibrium ──► KineticForces + kinetic profiles │ (Vacuum) ▲ │ ErrorFields + ├──► LocalStability │ │ + ├──► Match │ └────────────► gpec.h5 + └──► Tearing ForcingTerms (coils or file) ``` -The single `gpec.toml` file supplies user-selected options to every module. The two primary input files — the equilibrium file and (planned) the kinetic profiles file — are read by the `Equilibrium` module, while the external perturbation specification enters at `PerturbedEquilibrium` via `ForcingTerms`. `Vacuum` feeds both `ForceFreeStates` and `PerturbedEquilibrium` directly: the response matrices it builds are consumed by both downstream stages. +`Vacuum` runs inside `ForceFreeStates` whenever the solve is free-boundary, and `ForcingTerms` is read when `PerturbedEquilibrium` needs a drive, so neither is a stage of its own. +Every stage after `Equilibrium` reads the force-free states, and `KineticForces` and `ErrorFields` read the plasma response, so a section runs only when the sections it depends on are present. ---- +## Which stage do I need -## Stage 1: Equilibrium +| Question | Deck sections | Read first | +|---|---|---| +| Is the equilibrium read and mapped correctly; what are q and β_N? | `[Equilibrium]` (+ `[ForceFreeStates]`, which triggers the write) | `Equilibrium/q_95`, `Equilibrium/beta_N`, `Equilibrium/Profiles/` | +| Is the plasma ideal-kink stable, with or without a wall? | `[ForceFreeStates]` (+ `[Wall]`) | `ForceFreeStates/FreeBoundaryStability/W_freeboundary`, `…/eigenmode_energies` | +| Is it Mercier, resistive-interchange or ballooning unstable? | `[LocalStability]` (+ `[Ballooning]`) | `LocalStability/D_I`, `LocalStability/D_R`, `LocalStability/ballooning_Delta_prime` | +| What is the tearing stability index Δ′ at each rational surface? | `[ForceFreeStates] integrator = "Riccati"` or `"Galerkin"` (+ `[Riccati]` or `[Galerkin]`) | `SingularSurfaces/Delta_prime_matrix` | +| How does the plasma respond to an RMP or error field; do islands overlap? | `[ForcingTerms]` + `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `…/island_half_width`, `…/chirikov_parameter` | +| What NTV torque does that field drive? | `[KineticForces]` (+ `[Equilibrium] kinetic_file`) | `KineticForces//total_torque` | +| How tight must the coil tolerances be to avoid locking? | `[ErrorFields]` | `ErrorFields/Risk/locking_probability_percent` | +| What are the tearing growth rate and critical resonant field? | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `Tearing/Roots/gamma`, `Tearing/Roots/omega`, `Tearing/CriticalResonantField/` | +| What does a resistive inner layer do to the response? | `[Match]` + `[GGJ]` | `SingularSurfaces/Match/` | -**Module**: `Equilibrium` +### Headline outputs -**Purpose**: Reconstruct the axisymmetric MHD equilibrium. Solves the Grad-Shafranov equation on a poloidal flux grid, computes field-aligned coordinates, and builds bicubic splines for all field quantities used by downstream modules. +- **Free-boundary energy** `ForceFreeStates/FreeBoundaryStability/W_freeboundary` is the power-normalized free-boundary energy matrix, the sum of `W_plasma` and `W_vacuum` for the wall of the deck (no wall when `[Wall]` is absent). `eigenmode_energies` lists its generalized eigenvalues, the ideal δW per unit displacement of each free-boundary eigenmode, sorted most unstable first, so a negative first entry means the plasma is ideal-kink unstable and a positive one is the stability margin. +- **Tearing index** `SingularSurfaces/Delta_prime_matrix` is the Δ′ matching matrix in the PEST-3 convention, from which the outer-region tearing drive of each rational surface follows; a positive Δ′ means the outer region drives the tearing mode and the inner layer (`Tearing`) sets its growth rate. +- **Resonant field** `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field` is the area-weighted resonant normal field b̄ʳ `[T]` at each rational surface after the plasma response, the quantity that drives islands and sets locking thresholds. `island_half_width` is the vacuum-superposition island each resonant flux would open, w_s ∝ √|δψ_s|, and `chirikov_parameter` is the ratio of adjacent half-widths to the surface separation, so σ > 1 means the islands overlap. +- **Torque** `KineticForces//total_torque` is the neoclassical toroidal viscosity (NTV) torque T_φ `[N·m]` that the perturbed equilibrium exerts on the plasma, with the sign convention of the [Conventions](conventions.md) page. -**Configuration**: `[Equilibrium]` section of `gpec.toml` +## What each stage computes -**Inputs**: -- Equilibrium data file in one of the supported formats: - - `efit` — EFIT g-file from experimental reconstruction - - `chease` / `chease2` — CHEASE equilibrium code output - - `ifile` — inverse i-file (e.g. TokaMaker `save_ifile`) - - `lar` — Large aspect ratio analytical model - - `sol` — Solov'ev analytical equilibrium -- Kinetic profiles file (planned) — temperature and density profiles for the kinetic analysis path -- Grid resolution and solver settings +### Equilibrium -**Outputs** (`PlasmaEquilibrium`): -- Bicubic splines for R(ψ,θ), Z(ψ,θ), and all field quantities -- Safety factor profile q(ψ) and singular surface locations -- Global parameters: magnetic axis (R₀, Z₀), β, q₀, q_a, plasma volume -- Flux surface geometry and separatrix location +Reads or constructs the axisymmetric equilibrium, maps it to straight-field-line coordinates (`[Equilibrium] jac_type`), builds the flux-surface splines every later stage evaluates, and locates the rational surfaces q = m/n. +`[Equilibrium] eq_type` selects a file format: `efit`, `efit_arclength` and `efit_by_inversion` (EFIT g-files, three flux-surface tracing strategies), `chease` and `chease2` (CHEASE binary and ASCII), `ifile` (inverse i-file, as written by TokaMaker) and `imas` (an `IMASdd.dd` passed in from a script). +An analytic equilibrium is selected by its own section instead: `[Solovev]`, `[LargeAspectRatio]` or `[TJAnalytic]`. +The optional `[Equilibrium] kinetic_file` adds the density, temperature and rotation profiles that the kinetic solve, `KineticForces` and the `SLAYER` layer model need; `[[Equilibrium.ion_species]]` declares the ion composition. +Outputs: `Equilibrium/` scalars (`beta_N`, `q_axis`, `q_95`, …), `Equilibrium/Profiles/` and `Equilibrium/Geometry/`. +Page: [Equilibrium](equilibrium.md). ---- +### ForceFreeStates -## Stage 2: Vacuum Response +Determines ideal MHD stability by the Newcomb criterion of the Fortran DCON: the Euler–Lagrange equations of δW are integrated in ψ for all poloidal harmonics at once, crossing each rational surface with the prescribed small-solution jump, to give the force-free states ξ(ψ, θ) and the plasma energy matrix. +With `vac_flag = true` the `Vacuum` module adds the vacuum energy for the wall of the `[Wall]` section (shapes `conformal`, `elliptical`, `dee`, `mod_dee` and `file`; no wall when the section is absent), and the free-boundary eigenvalue problem gives the kink stability. +`[ForceFreeStates] integrator` selects the radial method and its section: `"Riccati"` (the default; `[Riccati]`) also produces the Δ′ matrix, `"Forward"` (`[Forward]`) is the plain ODE integration, and `"Galerkin"` (`[Galerkin]`) is the resistive outer-region solve on a finite-element basis. +`kinetic = "KineticForces"` replaces the ideal energy matrices by their drift-kinetic form, which needs a kinetic file. +Outputs: `ForceFreeStates/FreeBoundaryStability/`, `ForceFreeStates/Solutions/`, `ForceFreeStates/EdgeScan/`, `SingularSurfaces/` and `SurfaceGeometries/`. +Key references: [Glasser (2016)](citations.md#ForceFreeStates-Module), [Glasser (2018)](citations.md#ForceFreeStates-Module), [Chance et al. (1997)](citations.md#Vacuum-Module). +Pages: [Stability](stability.md), [Galerkin solver](galerkin.md), [Vacuum](vacuum.md). -**Module**: `Vacuum` +### LocalStability -**Purpose**: Compute the magnetostatic vacuum response matrices that couple the plasma boundary to external fields. These matrices are used by ForceFreeStates and PerturbedEquilibrium to apply boundary conditions and decompose fields in poloidal mode space. +Evaluates the Mercier criterion D_I, the resistive-interchange criterion D_R and the high-n ballooning Δ′ on every flux surface of the solve, and with `[LocalStability] alpha_boundary_scan` the critical normalized pressure gradient α_crit. +Runs only with `[ForceFreeStates]`, on the same equilibrium grid; `[Ballooning]` sets the ballooning integration. +Outputs: `LocalStability/`. +Page: [Local stability](ballooning.md). -**Configuration**: optional `[Wall]` section of `gpec.toml`; omit it for no wall +### Match -**Inputs**: -- `PlasmaEquilibrium` from Stage 1 (plasma boundary coordinates, toroidal mode number n) -- Wall geometry: shape (conformal, elliptical, dee-shaped, or custom) and dimensions +Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, section `[GGJ]`) and rebuilds the force-free states with the resulting layer response, so that the later stages see a resistive rather than perfectly shielded plasma; `ideal = true` keeps the shielded limit for comparison. +Outputs: `SingularSurfaces/Match/`. +Page: [Tearing](inner_layer.md). -**Outputs**: -- `wv` — Vacuum response matrix (scaled by the singular factor (m - nq)(m' - nq), see Chance 1997) -- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_Iv=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L` +### PerturbedEquilibrium -**Key references**: [Chance et al. (1997)](citations.md#Vacuum-Module), [Chance et al. (2007)](citations.md#Vacuum-Module) +Computes the plasma response to an external non-axisymmetric field: `plasma_response` builds the permeability, inductance and reluctance operators on the control surface from the force-free states and the vacuum matrices, and `perturbed_equilibrium` applies them to the drive that `ForcingTerms` supplies, either a normal-field spectrum from a file or the Biot–Savart field of `[[ForcingTerms.coil_set]]` coils. +The resonant fields, island half-widths and Chirikov parameter at each rational surface follow from the singular coupling (`compute_singular_coupling`), and `regularization_width` controls the regularization of the ideal response there. +Outputs: `PerturbedEquilibrium/ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`. +Key references: [Park et al. (2007a)](citations.md#PerturbedEquilibrium-Module), [Park et al. (2009)](citations.md#PerturbedEquilibrium-Module). +Pages: [Forcing terms](forcing_terms.md), [Perturbed equilibrium](perturbed_equilibrium.md). ---- +### KineticForces -## Stage 3: Ideal MHD Stability +Integrates the NTV torque and kinetic energy of the perturbed equilibrium over ψ for each torque `method`, by species when the ion composition has more than one. +Runs after `PerturbedEquilibrium` and needs n_e, T_i, T_e and ω_E in the kinetic file; without `[PerturbedEquilibrium]` the section only configures the kinetic solve of `ForceFreeStates`. +Outputs: `KineticForces//` and `KineticForces/PerSpecies/`. +Key references: [Park and Logan (2017)](citations.md#Kinetic-Forces), [Logan et al. (2013)](citations.md#Kinetic-Forces). +Page: [Kinetic forces](kinetic_forces.md). -**Module**: `ForceFreeStates` +### ErrorFields -**Purpose**: Determine ideal MHD stability via Newcomb's direct criterion. Integrates the Euler-Lagrange equations across the plasma volume, crossing singular surfaces (where q = m/n) with prescribed jump conditions. Computes the potential and kinetic energy matrices that define the eigenvalue problem for the free-boundary stability. +Assesses error-field tolerances for a coil-driven response: the sensitivity of the resonant field to the shift and tilt of each coil set, a Monte Carlo sample of the tolerances of `[ErrorFields.MonteCarlo]`, the locking probability of `[ErrorFields.Risk]` against a threshold scenario, and the NTV coupling of the correction coils of `[ErrorFields.NTV]`. +Needs `[PerturbedEquilibrium]` driven by `[[ForcingTerms.coil_set]]` coils. +Outputs: `ErrorFields/CoilSensitivities/`, `MonteCarlo/`, `Risk/`, `NTV/`. +Page: [Error fields](error_fields.md). -**Configuration**: `[ForceFreeStates]` section of `gpec.toml` +### Tearing -**Inputs**: -- `PlasmaEquilibrium` from Stage 1 -- Vacuum response matrices from Stage 2 -- Poloidal mode range (m_low, m_pert) and toroidal mode number n +Solves the tearing dispersion relation of each rational surface by matching the outer-region Δ′ to an inner-layer model: the two-fluid slab `SLAYER` (`[SLAYER]`) or the toroidal resistive-MHD `GGJ` (`[GGJ]`), selected by `[Tearing] inner_model`. +`[Tearing.search]` sets the root search in the complex frequency plane, and `[Tearing.CriticalResonantField]` adds the resonant field at which the layer's shielding is overcome. +Needs a Δ′-producing integrator (`"Riccati"` or `"Galerkin"`) and, for `SLAYER`, a kinetic file. +Outputs: `Tearing/PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`. +Key references: [Glasser, Wang and Park (2016)](citations.md#InnerLayer-Module), [Park (2022)](citations.md#InnerLayer-Module). +Page: [Tearing](inner_layer.md). -**Outputs**: -- Force-free eigenfunctions ξ(ψ, θ) — the normal displacement field -- Fixed boundary energy W_fixed (negative → unstable without wall) -- Free boundary energy eigenvalues `et` (sign determines ideal stability with wall) -- Mercier stability criterion as a function of ψ -- Singular surface data at each rational surface q = m/n: - - Flux surface location ψ_s - - Tearing stability parameter Δ' - - Small solution coefficients (used by PerturbedEquilibrium) +## The deck -**Key references**: [Glasser (2016) Newcomb](citations.md#ForceFreeStates-Module), [Glasser (2018) Riccati](citations.md#ForceFreeStates-Module) - ---- - -## Stage 4: External Forcing - -**Module**: `ForcingTerms` - -**Purpose**: Load the external magnetic perturbation specification — the amplitude and phase of each (m, n) Fourier component of the applied field at the plasma boundary. This is typically computed externally (e.g. from a coil geometry code or measured from sensors). - -**Configuration**: `[ForcingTerms]` section of `gpec.toml` - -**Inputs**: -- External perturbation data file (ASCII or HDF5 format) - -**Outputs** (array of `ForcingMode`): -- For each (m, n) mode: poloidal mode number, toroidal mode number, complex amplitude (magnitude + phase) - ---- - -## Stage 5: Perturbed Equilibrium - -**Module**: `PerturbedEquilibrium` - -**Purpose**: Compute the self-consistent plasma response to the external perturbation specified in Stage 4, in two steps: -- `plasma_response(ffs)` builds the source-independent operators once from the force-free eigenfunctions and vacuum matrices: the permeability, inductance and reluctance matrices on the control surface, and the singular coupling matrices at each rational surface (`PlasmaResponse`). -- `perturbed_equilibrium(resp, rmp)` applies them to one forcing field: it projects the response into mode space and evaluates the singular coupling diagnostics for that field. One response serves any number of fields. - -**Configuration**: `[PerturbedEquilibrium]` section of `gpec.toml` - -**Inputs**: -- ForceFreeStates results from Stage 3 (eigenfunctions, singular surface data), for `plasma_response` -- Vacuum response matrices from Stage 2, for `plasma_response` -- The forcing field from Stage 4, for `perturbed_equilibrium` - -**Outputs** (`PerturbedEquilibriumState`): -- Mode-space normal displacement ξ(m, ψ) across the plasma -- Mode-space magnetic perturbation b(m, ψ) -- At each rational surface q = m_s/n: - - Resonant flux δψ_s - - Resonant current sheet amplitude - - Island half-width w_s (proportional to √|δψ_s|) - - Chirikov overlap parameter σ (ratio of adjacent island widths to their separation) - -**Key references**: [Park et al. (2007a)](citations.md#PerturbedEquilibrium-Module), [Park et al. (2009)](citations.md#PerturbedEquilibrium-Module), [Park et al. (2017)](citations.md#Kinetic-Forces) - ---- - -## Configuration File: `gpec.toml` - -All stages are controlled by a single TOML file. A minimal example structure: +A stage runs when its section is present, so a deck with only `[Equilibrium]` and `[ForceFreeStates]` stops after the stability solve, and one without `[ForceFreeStates]` writes nothing. +The abridged `examples/DIIID-like_ideal_example/gpec.toml` below runs the equilibrium, the ideal free-boundary solve with local stability, and the coil-driven perturbed equilibrium: ```toml [Run] -output_file = "gpec.h5" verbose = true [Equilibrium] eq_type = "efit" -eq_filename = "g-file.txt" +eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" +kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5" +jac_type = "hamada" +psihigh = 0.995 -[Wall] -shape = "conformal" -standoff_fraction = 0.3 +[LocalStability] +alpha_boundary_scan = true [ForceFreeStates] nn = 1 -integrator = "Forward" vac_flag = true +integrator = "Forward" # the default is "Riccati" + +[ForcingTerms] +forcing_data_format = "coil" +machine = "d3d" + +[[ForcingTerms.coil_set]] +name = "c" +currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] [PerturbedEquilibrium] compute_singular_coupling = true - -[ForcingTerms] -forcing_data_file = "forcing.dat" -forcing_data_format = "ascii" +regularization_width = 0.05 ``` -A stage runs when its section is present, so omitting the later sections stops the pipeline early -(an equilibrium-only deck has no `[ForceFreeStates]`). Every stage saves into the `[Run] output_file`. - -Example configuration files are provided in: -- `examples/Solov'ev_ideal_example/gpec.toml` -- `examples/DIIID-like_ideal_example/gpec.toml` -- `examples/DIIID-like_error_field_example/gpec.toml` (coil-forced run with the `[ErrorFields]` sensitivity stage) +Every key is documented on the module page of its section, and every example deck carries an inline comment per key. + +### Deck section index + +| Section | Runs or configures | Page | +|---|---|---| +| `[Run]` | `output_file` (default `gpec.h5`, written next to the deck) and `verbose` | this page | +| `[Equilibrium]` | the equilibrium stage: file, format, coordinates, grid, kinetic file | [Equilibrium](equilibrium.md) | +| `[[Equilibrium.ion_species]]` | the ion composition of the kinetic profiles | [Equilibrium](equilibrium.md) | +| `[Solovev]`, `[LargeAspectRatio]`, `[TJAnalytic]` | an analytic equilibrium in place of a file | [Equilibrium](equilibrium.md) | +| `[Wall]` | the conducting wall of the vacuum response; omit for no wall | [Vacuum](vacuum.md) | +| `[ForceFreeStates]` | the stability stage: `nn`, the mode window, `integrator`, `vac_flag`, `kinetic` | [Stability](stability.md) | +| `[Forward]`, `[Riccati]` | the radial ODE integrators | [Stability](stability.md) | +| `[Galerkin]` | the finite-element outer-region solver | [Galerkin solver](galerkin.md) | +| `[LocalStability]`, `[Ballooning]` | the local stability stage and its ballooning integration | [Local stability](ballooning.md) | +| `[Match]`, `[GGJ]` | the inner-layer matching stage and the GGJ layer model | [Tearing](inner_layer.md) | +| `[ForcingTerms]`, `[[ForcingTerms.coil_set]]` | the external field: a file or coil sets | [Forcing terms](forcing_terms.md) | +| `[PerturbedEquilibrium]` | the plasma response stage and its singular-coupling options | [Perturbed equilibrium](perturbed_equilibrium.md) | +| `[KineticForces]` | the NTV torque stage and the kinetic solve of `ForceFreeStates` | [Kinetic forces](kinetic_forces.md) | +| `[Tearing]`, `[Tearing.search]`, `[Tearing.CriticalResonantField]`, `[SLAYER]` | the tearing stage, its root search, the critical field and the SLAYER layer model | [Tearing](inner_layer.md) | +| `[ErrorFields]`, `[ErrorFields.MonteCarlo]`, `[ErrorFields.scenario]`, `[ErrorFields.Risk]`, `[ErrorFields.NTV]` | the error-field tolerance stage | [Error fields](error_fields.md) | + +## Running it + +From a shell, on a directory holding a `gpec.toml`: + +```bash +./gpec examples/DIIID-like_ideal_example +``` ---- +From Julia, the same pipeline is `GeneralizedPerturbedEquilibrium.main(["examples/DIIID-like_ideal_example"])`, which returns the stage results as a named tuple; each stage is also a `solve` call on a problem type, so a script can run any subset. See the [Scripting API](api.md). +`examples/*/run_example.jl` and `analyze_example.jl` show the scripted run and the [Analysis](analysis.md) plots of its output. -## Output File: `gpec.h5` +### Examples -All results are written to a single HDF5 file (default: `gpec.h5`). The top-level groups are organized by physics topic (see the schema conventions in `docs/development/hdf5-conventions.md`): +| Directory | Exercises | +|---|---| +| `DIIID-like_ideal_example` | DIII-D-like H-mode, n = 1: free-boundary ideal stability, local stability, C-coil-driven perturbed equilibrium and NTV torque; the default benchmark case | +| `DIIID-like_ideal_multi_n_example` | the same plasma with n = 1 and 2 solved together and driven by mixed-n coil currents | +| `DIIID-like_ideal_example_IMAS` | the same equilibrium read from an `IMASdd.dd`, run from `run_imas_example.jl` | +| `DIIID-like_riccati_deltaprime_example` | the Δ′ matrix on the `"Riccati"` integrator, no forcing | +| `DIIID-like_gal_resistive_example` | the `"Galerkin"` outer-region Δ′ solve, no forcing | +| `DIIID-like_gal_resistive_pe_example` | Galerkin Δ′ with `[Match]` inner-layer matching and coil forcing | +| `DIIID-like_SLAYER_example` | tearing growth rates and critical resonant fields with the `SLAYER` layer model | +| `DIIID-like_truncation_example` | the resistive-layer overlap cap on the edge of the integration domain | +| `DIIID-like_error_field_example` | error-field sensitivities, Monte Carlo tolerances and locking risk for the DIII-D poloidal-field coils | +| `Solovev_ideal_example` | analytic Solov'ev equilibrium from `[Solovev]`: a near-marginal n = 1 kink held by a conformal wall | +| `Solovev_ideal_example_3D` | the same plasma with the vacuum response on a toroidal grid | +| `Solovev_ideal_example_multi_n` | the same plasma swept over n = 1 and 2 in one run, with single-n cross-checks | +| `LAR_reversed_shear_example` | a `[LargeAspectRatio]` equilibrium with reversed shear, q = 2 crossed twice | +| `LAR_ideal_match_test` | a `[TJAnalytic]` equilibrium with `[Match] ideal = true`, the shielded reference for the resistive case | +| `LAR_resistive_match_test` | the same equilibrium with resistive GGJ matching | +| `LAR_beta_scan`, `LAR_epsilon_scan` | `run_scan.jl` sweeps of the `[TJAnalytic]` pressure factor and inverse aspect ratio | +| `TJ_ideal_example` | a `[TJAnalytic]` equilibrium: n = 1 stability and a file-driven perturbed equilibrium | +| `TJ_kinetic_NTV_example` | adds `[KineticForces]`: the NTV torque of the perturbed equilibrium | +| `TJ_kinetic_calculated_example` | the drift-kinetic stability solve (`kinetic = "KineticForces"`) with a D–T main ion | +| `TJ_kinetic_multiion_example` | NTV summed over D, T, the impurity and electrons | +| `a10_kinetic_example` | drift-kinetic stability of a low-aspect-ratio circular tokamak | +| `efit_fixedbdy_separatrix_example` | a fixed-boundary EFIT file whose box hugs the separatrix, for the separatrix search | + +## The output file + +All results are written to one HDF5 file (`[Run] output_file`, default `gpec.h5`, next to the deck), with top-level groups by physics topic; the schema conventions are in `docs/development/hdf5-conventions.md` of the repository. | Group | Contents | |---|---| -| `Info/` | Run metadata: git version, mode-number ranges, ψ limit, `Runtimes/` (per-stage wall-clock seconds) | -| `Input/` | Self-contained rerun snapshot: the recorded deck, raw equilibrium/forcing/coil inputs | -| `Equilibrium/` | Equilibrium scalars (`beta_N`, `q_axis`, `q_95`, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) | -| `ForceFreeStates/` | Stability solve: `Solutions/{ForwardIntegration,GalerkinIntegration}`, `EulerLagrangeMatrices/` (opt-in: `write_el_matrices = true`), `FreeBoundaryStability/`, `EdgeScan/` | -| `LocalStability/` | Mercier D_I, resistive interchange D_R, ballooning Δ' profiles | -| `SingularSurfaces/` | Per-rational-surface data: ψ_s, q, m/n, GGJ coefficients, Δ'/PEST-3 matching matrices, kinetic surfaces (`Kinetic/`) | -| `PerturbedEquilibrium/` | Plasma response: `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/` | -| `KineticForces/` | NTV torque per method: energy integrals, kinetic matrices | -| `Tearing/` | SLAYER/GGJ inner-layer growth rates: `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `Scan/` | +| `Info/` | Run metadata: git version, mode-number ranges, `psilim`, `qlim`, `Runtimes/` (per-stage wall-clock seconds) | +| `Input/` | Rerun snapshot: the recorded deck and the raw equilibrium, forcing and coil inputs | +| `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) | +| `ForceFreeStates/` | `Solutions/`, `FreeBoundaryStability/`, `EdgeScan/`, and `EulerLagrangeMatrices/` when `write_el_matrices = true` | +| `LocalStability/` | `D_I`, `D_R`, `ballooning_Delta_prime` on `psi`; the α boundary on `ballooning_psi` | +| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`, `rational_q`, GGJ coefficients, `Delta_prime_matrix` and the PEST-3 matrices, `Match/`, `Kinetic/` | +| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/` | +| `KineticForces/` | `/` torque and energy profiles; `PerSpecies///` for multi-ion runs | +| `ErrorFields/` | `CoilSensitivities/`, `MonteCarlo/`, `Risk/`, `NTV/` | +| `Tearing/` | `PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/` | | `SurfaceGeometries/` | Plasma and wall surface point clouds for visualization | + +Rerunning `main(["gpec.h5"])` on an output file reproduces the run from its `Input/` snapshot; see the [Scripting API](api.md). From c3923246904adf96d87da1e5edf2db620534605f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 2 Nov 2026 00:27:19 -0500 Subject: [PATCH 02/10] Docs - DOCS - Rewrite the stability, local stability, Galerkin and perturbed equilibrium pages around their physics Give stability.md a what-it-computes opener, structural governing equations with the explicit (m - nq) convention, an outputs table and a short domain and layer-overlap section in place of the integrator and diffusivity essays; the per-surface Delta-prime stub and the Notes that contradicted the example go. Make ballooning.md the Local Stability page: criteria, condensed ballooning theory, the run and outputs block moved from stability.md, and numerical notes without the equivalence proof or the implementation history. Correct the Glasser 2016 article numbers on the Galerkin page, state the Q convention once and add its deck snippet and outputs. Prepend the response physics (P, L, reluctance, resonant field, islands, Chirikov, energies) and a units table to perturbed_equilibrium.md. Add the LocalStability citations, move the SLAYER diffusivity inputs to the Tearing page, and fix the Match and Tearing prerequisites on the workflow page. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- docs/src/ballooning.md | 935 ++---------------------------- docs/src/citations.md | 16 + docs/src/galerkin.md | 88 ++- docs/src/inner_layer.md | 20 +- docs/src/perturbed_equilibrium.md | 104 +++- docs/src/stability.md | 543 ++++------------- docs/src/workflow.md | 14 +- 7 files changed, 336 insertions(+), 1384 deletions(-) diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md index 8ef6cd10bc..b21c217891 100644 --- a/docs/src/ballooning.md +++ b/docs/src/ballooning.md @@ -1,918 +1,99 @@ -# Ballooning and Mercier Local Stability +# Local Stability -This note describes the physical and mathematical foundations of `Bal.jl`, which -computes local Mercier and high-``n`` ballooning stability. It covers the -ballooning-mode equation, the asymptotic boundary treatment, and the local -perturbation framework used to generate the ``s``-``\alpha`` diagram. +`LocalStability` evaluates the local, high-``n`` ideal and resistive criteria on every flux surface of the equilibrium: the Mercier interchange index ``D_I``, the Glasser–Greene–Johnson resistive-interchange index ``D_R`` and the ballooning stability index, with an optional scan for the critical pressure gradient. +It runs as its own stage, after `ForceFreeStates` and on the same equilibrium, and is the local counterpart of the global analysis on the [Stability](stability.md) page. -## Ballooning-mode equation +## What it computes -### Ballooning representation +- **Mercier criterion** ``D_I``: the ideal interchange index of each surface. A surface is Mercier unstable where ``D_I > 0``; the same number sets the exponents ``-1/2 \pm \sqrt{-D_I}`` of the large and small solutions that the global integration matches at every rational surface. +- **Resistive interchange** ``D_R = D_I + (H - 1/2)^2``, with ``H`` from the flux-surface averages of [Glasser, Greene and Johnson (1975)](citations.md#LocalStability-Module). ``D_R > 0`` indicates resistive interchange instability. +- **Ballooning stability** ``\Delta'``: the high-``n`` ballooning index, obtained by integrating the marginal ballooning equation along the field line and matching the solution to its asymptotic forms at both ends ([Connor, Hastie and Taylor (1978)](citations.md#LocalStability-Module)). It is distinct from the tearing Δ′ of `SingularSurfaces/`; ``\Delta' = 0`` is the marginal ballooning boundary. +- **Critical pressure gradient**: with `alpha_boundary_scan = true` the pressure gradient of each surface is scaled until the ballooning boundary is reached, giving the normalized gradient ``\alpha`` and its critical value ``\alpha_\mathrm{crit}``, the s–α diagram collapsed to a profile. -In a tokamak equilibrium, physical quantities are periodic in the poloidal angle -``\theta`` and toroidal angle ``\zeta``. The perturbation function ``\varphi`` -satisfies +## Governing equations -```math -\varphi(\psi, \theta, \beta) -= \varphi(\psi, \theta, \beta + 1) -= \varphi(\psi, \theta + 1, \beta - q), -``` - -where ``q`` is the safety factor. - -Ballooning modes have long wavelengths along magnetic field lines and short -wavelengths perpendicular to the field. They are represented using a -covering-space function ``\hat{\varphi}``, defined on -``-\infty < \theta < \infty``: - -```math -\varphi(\psi, \theta, \beta) -= \sum_{l=-\infty}^{\infty} -\hat{\varphi}(\psi, \theta - l) -\exp[-2\pi i n(\beta + lq)]. -``` - -Here - -```math -\beta = \zeta - q\theta + \int q'\theta_0(\psi)\,d\psi -``` - -is the Clebsch coordinate. The wave-number vector contains a secular term from -the magnetic shear ``q'``: - -```math -|\nabla\beta|^2 -= |\nabla\zeta - q\nabla\theta|^2 -+ q'^2(\theta-\theta_0)^2|\nabla\psi|^2 -- 2q'(\theta-\theta_0)(\nabla\zeta-q\nabla\theta)\cdot\nabla\psi. -``` - -### Ideal marginal ballooning equation - -At the stability boundary, ``\gamma = 0``. Variation of the ideal-MHD energy -principle gives the ballooning ODE - -```math -\mathbf{B}\cdot\nabla -\left( - \frac{|\nabla\beta|^2}{B^2} - \mathbf{B}\cdot\nabla\hat{\varphi} -\right) -+ 2\frac{P'}{\chi'}\kappa_w\hat{\varphi} -= 0, -``` - -where ``\chi`` is the poloidal flux and ``P'`` is the pressure-gradient term used -by the local-stability equations. Written using the Jacobian ``\mathcal{J}``, - -```math -\frac{\chi'}{\mathcal{J}}\frac{\partial}{\partial\theta} -\left( - \frac{|\nabla\beta|^2}{B^2} - \frac{\chi'}{\mathcal{J}} - \frac{\partial\hat{\varphi}}{\partial\theta} -\right) -+ 2\frac{P'}{\chi'}\kappa_w\hat{\varphi} -= 0. -``` - -`Bal.jl` converts this to a first-order matrix equation, - -```math -\frac{\partial\mathbf{u}}{\partial\theta} = \mathbf{L}\mathbf{u}, -\qquad -\mathbf{u} = -\begin{pmatrix} -\hat{\varphi} \\ -\dfrac{|\nabla\beta|^2}{B^2} -\dfrac{\chi'}{\mathcal{J}} -\dfrac{\partial\hat{\varphi}}{\partial\theta} -\end{pmatrix}, -``` - -with - -```math -\mathbf{L} -= \frac{\mathcal{J}}{\chi'} -\begin{pmatrix} -0 & \dfrac{B^2}{|\nabla\beta|^2} \\ --2\dfrac{P'}{\chi'}\kappa_w & 0 -\end{pmatrix}. -``` - -### Curvature decomposition - -The instability drive is the projected curvature - -```math -\kappa_w -\equiv -\frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\beta. -``` - -Using -``\nabla\beta = \nabla\zeta - q\nabla\theta - q'\theta\nabla\psi``, -the curvature separates into non-secular and secular pieces: - -```math -\kappa_w = \kappa_n - q'\theta\kappa_s, -``` - -where - -```math -\kappa_s -\equiv -\frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\psi -= -\frac{\pi f}{\mathcal{J}} -\frac{\partial}{\partial\theta} -\left(\frac{1}{B^2}\right). -``` - -This decomposition is the main structure used to assemble the ballooning -coefficient matrix. - -## Asymptotic basis and Mercier exponent - -In the ``\theta\to\infty`` limit, the coefficient matrix diverges because -``|\nabla\beta|^2 \simeq q'^2\theta^2|\nabla\psi|^2``. `Bal.jl` applies the -canonical transformation - -```math -\mathbf{R} = \operatorname{diag}(\theta^{-1/2}, \theta^{1/2}) -``` - -and solves - -```math -\frac{\partial\mathbf{v}}{\partial\theta} = \mathbf{M}\mathbf{v}, -\qquad -\mathbf{M} -= -\begin{pmatrix} -\dfrac{1}{2\theta} -& -\dfrac{\mathcal{J}}{\chi'} -\dfrac{B^2\theta}{|\nabla\beta|^2} -\\ --\dfrac{\mathcal{J}}{\chi'} -\dfrac{2P'\kappa_w}{\chi'\theta} -& --\dfrac{1}{2\theta} -\end{pmatrix}. -``` - -To keep the lowest-order matrix finite, the code applies - -```math -\mathbf{S} -= -\begin{pmatrix} -1 & 0 \\ -\sigma & 1 -\end{pmatrix}, -\qquad -\sigma = --\frac{2\pi f P'q'}{\chi'^2 B^2}. -``` - -The transformed matrices ``\mathbf{N}_k = \mathbf{S}^{-1}\mathbf{M}_k\mathbf{S}`` -enter the recurrence relation - -```math -\frac{\partial\mathbf{w}^{(k+1)}}{\partial\theta} -= -[\mathbf{N}_0-(\alpha+k)\mathbf{I}]\mathbf{w}^{(k)} -+ \sum_{l=1}^{k}\mathbf{N}_l\mathbf{w}^{(k-l)}. -``` - -The Mercier exponent ``\alpha`` follows from the zeroth-order solvability -condition: - -```math -\det[\langle\mathbf{N}_0\rangle-\alpha\mathbf{I}] = 0, -\qquad -\alpha = \pm\sqrt{-D_I}. -``` - -### Equivalence with the surface-average Mercier formula - -The previous `Mercier.jl` evaluated the Mercier -index from flux-surface averages. The two routes are mathematically equivalent: -they are different reductions of the same large-``\theta`` ballooning -solvability condition. - -Let the flux-surface average used by the old implementation be - -```math -\langle X\rangle_{\mathrm{fs}} -\equiv -\int_0^1 \frac{\mathcal{J}}{V'} X(\theta)\,d\theta, -``` - -and define - -```math -A_1 = \left\langle\frac{B^2}{|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_2 = \left\langle\frac{1}{|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -``` - -```math -A_3 = \left\langle\frac{1}{B^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_4 = \left\langle\frac{1}{B^2|\nabla\psi|^2}\right\rangle_{\mathrm{fs}}, -\qquad -A_5 = \left\langle B^2\right\rangle_{\mathrm{fs}}. -``` - -Below, an overbar denotes the period integral used to form the `Bal.jl` -``\bar{d}_0`` matrix: - -```math -\overline{X} \equiv \int_0^1 X(\theta)\,d\theta. -``` - -With ``F = 2\pi f``, the old code evaluated - -```math -T = -\frac{F P' V'}{q'\chi'^3} A_2, -``` - -and then - -```math -D_I^{\mathrm{avg}} -= --\frac14 -+ T(1-T) -+ P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -This is the surface-average Mercier expression that appeared in the old -`Mercier.jl` and, temporarily, in the earlier `resistive_interchange` helper. -The associated GGJ resistive-interchange quantity used - -```math -H = -\frac{F P' V'}{q'\chi'^3} -\left( - A_2 - \frac{A_1}{A_5} -\right), -\qquad -D_R = D_I + \left(H-\frac12\right)^2. -``` - -The determinant route starts from the same ballooning equation, takes the -large-``\theta`` limit, applies the ``\mathbf{R}`` and ``\mathbf{S}`` -transformations, and obtains the periodic lowest-order system - -```math -\frac{\partial\mathbf{w}^{(1)}}{\partial\theta} -= -(\mathbf{N}_0-\alpha\mathbf{I})\mathbf{w}^{(0)}. -``` - -Solvability requires the mean operator to have a null vector: - -```math -[\overline{\mathbf{N}_0}-\alpha\mathbf{I}]\mathbf{w}^{(0)}=0. -``` - -The transformed matrix is traceless by construction, - -```math -\operatorname{tr}\overline{\mathbf{N}_0} = 0, -``` - -so the characteristic equation reduces to - -```math -\alpha^2 + \det\overline{\mathbf{N}_0} = 0. -``` - -Comparing this with the Mercier exponent definition -``\alpha^2 = -D_I`` gives - -```math -D_I = \det\overline{\mathbf{N}_0}. -``` - -The reduction can be written explicitly. For the unperturbed Mercier scan, -``\mathbf{N}_0`` has the code form - -```math -\mathbf{N}_0(\theta) -= -\begin{pmatrix} -\frac12 + M_{12}\sigma & M_{12} \\ -M_{21}-\sigma-M_{12}\sigma^2 & -\frac12-M_{12}\sigma -\end{pmatrix}, -``` - -where - -```math -M_{12} -= -\frac{\mathcal{J}}{\chi'} -\frac{B^2}{q'^2|\nabla\psi|^2}, -\qquad -M_{21} -= --2\frac{\mathcal{J}}{\chi'^2}P'\kappa_n, -\qquad -\sigma -= --\frac{F P' q'}{\chi'^2 B^2}. -``` - -The averaged matrix can therefore be written as - -```math -\bar{d}_0 -\equiv -\overline{\mathbf{N}_0} -= -\begin{pmatrix} -\frac12+B & A \\ -K & -\frac12-B -\end{pmatrix}, -``` - -with - -```math -A = \overline{M_{12}}, -\qquad -B = \overline{M_{12}\sigma}, -\qquad -K = \overline{M_{21}-\sigma-M_{12}\sigma^2}. -``` - -The first two averages are direct: - -```math -A -= -\frac{V'}{q'^2\chi'}A_1, -``` - -```math -B -= -\overline{ - \frac{\mathcal{J}}{\chi'} - \frac{B^2}{q'^2|\nabla\psi|^2} - \left( - -\frac{F P' q'}{\chi'^2B^2} - \right) -} -= --\frac{F P' V'}{q'\chi'^3}A_2 -= --T. -``` - -For ``K``, separate the three pieces: - -```math -K -= -\overline{M_{21}} --\overline{\sigma} --\overline{M_{12}\sigma^2}. -``` - -The last term is also direct: - -```math -\overline{M_{12}\sigma^2} -= -\overline{ - \frac{\mathcal{J}}{\chi'} - \frac{B^2}{q'^2|\nabla\psi|^2} - \frac{F^2P'^2q'^2}{\chi'^4B^4} -} -= -\frac{F^2P'^2V'}{\chi'^5}A_4. -``` - -The remaining curvature combination reduces by the axisymmetric equilibrium -identity used in the Mercier derivation: - -```math -\int_0^1 -\left( - 2\mathcal{J}\kappa_n - - \frac{Fq'}{B^2} -\right)d\theta -= -\frac{V'}{\chi'} -\left( - P'A_3 - \frac{V''}{V'} -\right). -``` - -Using ``M_{21}=-2\mathcal{J}P'\kappa_n/\chi'^2`` and -``\sigma=-FP'q'/(\chi'^2B^2)``, this gives - -```math -\overline{M_{21}}-\overline{\sigma} -= --\frac{P'}{\chi'^2} -\int_0^1 -\left( - 2\mathcal{J}\kappa_n - - \frac{Fq'}{B^2} -\right)d\theta -= --\frac{P'V'}{\chi'^3} -\left( - P'A_3 - \frac{V''}{V'} -\right). -``` - -Therefore - -```math -K -= --\frac{P'V'}{\chi'^3} -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -Now take the determinant: - -```math -\det\bar{d}_0 -= -\left(\frac12+B\right)\left(-\frac12-B\right)-AK. -``` - -Substituting ``B=-T`` gives the first part of the old Mercier expression: - -```math -\left(\frac12-T\right)\left(-\frac12+T\right) -= --\frac14+T(1-T). -``` - -Substituting ``A`` and ``K`` gives the second part: - -```math --AK -= -\frac{V'}{q'^2\chi'}A_1 -\frac{P'V'}{\chi'^3} -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right] -``` - -```math -= -P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right]. -``` - -Thus the determinant route gives exactly - -```math -\det\bar{d}_0 -= --\frac14 -+ T(1-T) -+ P'\left(\frac{V'}{q'\chi'^2}\right)^2 A_1 -\left[ - P'\left(A_3+\frac{F^2}{\chi'^2}A_4\right) - - \frac{V''}{V'} -\right] -= -D_I^{\mathrm{avg}}. -``` - -So, under the same normalization and exact equilibrium identities, - -```math -D_I^{\mathrm{det}} -\equiv -\det(\bar{d}_0) -= -D_I^{\mathrm{avg}}. -``` - -The implementation difference is numerical rather than physical. The old -surface-average route evaluates the closed Mercier formula directly, while the -new `Bal.jl` route constructs the same Mercier index from the averaged -large-``\theta`` local-mode matrix. In finite-precision calculations the two can -separate near the magnetic axis or other ill-conditioned surfaces because both -forms contain factors such as ``1/q'`` and ``1/q'^2``. Away from those regions, -they should agree to the accuracy of the equilibrium derivatives and surface -quadrature. - -The divergent big solution and decaying small solution have leading forms - -```math -\mathbf{u}_{\mathrm{large}}^{(0)} -\simeq -\begin{pmatrix} -\theta^{\alpha-1/2} \\ -\theta^{\alpha+1/2}(\sigma+w_{21}) -\end{pmatrix}, -\qquad -\mathbf{u}_{\mathrm{small}}^{(0)} -\simeq -\begin{pmatrix} -\theta^{-\alpha-1/2} \\ -\theta^{-\alpha+1/2}(\sigma+w_{22}) -\end{pmatrix}. -``` - -## First-order asymptotic matching - -The zeroth-order large-``\theta`` solution alone is not accurate enough at a -finite integration boundary ``\theta_{\max}``. `Bal.jl` therefore includes the -first-order correction ``\mathbf{w}^{(1)}``: - -```math -\frac{\partial\mathbf{w}^{(1)}}{\partial\theta} -= -(\mathbf{N}_0-\alpha\mathbf{I})\mathbf{w}^{(0)} -= -\tilde{\mathbf{N}}_0(\theta)\mathbf{w}^{(0)}. -``` - -The first-order solution is written as - -```math -\mathbf{w}^{(1)}(\theta) -= -\left(\int\tilde{\mathbf{N}}_0(\theta)\,d\theta\right)\mathbf{w}^{(0)} -+ \mathbf{c}^{(1)} -\equiv -\tilde{\mathbf{w}}^{(1)}(\theta) + \mathbf{c}^{(1)}. -``` - -The constant vector is fixed by the first-order solvability condition, - -```math -[\langle\mathbf{N}_0\rangle-(\alpha+1)\mathbf{I}]\mathbf{c}^{(1)} -= -- -\left\langle - [\mathbf{N}_0-(\alpha+1)\mathbf{I}] - \tilde{\mathbf{w}}^{(1)}(\theta) - + \mathbf{N}_1\mathbf{w}^{(0)} -\right\rangle. -``` - -This gives the asymptotic boundary form - -```math -\mathbf{u}_{\mathrm{boundary}} -\simeq -\mathbf{R}\mathbf{S} -\left( - \mathbf{w}^{(0)} - + \frac{1}{\theta}\mathbf{w}^{(1)}(\theta) -\right). -``` - -## Local thin-layer ordering - -The ``s``-``\alpha`` diagram scans local instability near a flux surface -``\psi_b`` without recomputing the full equilibrium. A thin-layer parameter -``0 < \mu \ll 1`` is introduced through - -```math -s = \frac{\psi-\psi_b}{\mu}, -\qquad -\frac{\partial}{\partial\psi} -= -\frac{1}{\mu}\frac{\partial}{\partial s}. -``` - -Profiles are expanded as background plus local perturbation, - -```math -P(\psi) = P^{(0)}(\psi) + \mu P^{(1)}(s) + \cdots, -\qquad -q(\psi) = q^{(0)}(\psi) + \mu q^{(1)}(s) + \cdots. -``` - -After differentiation, the ``\mu`` amplitude cancels the ``1/\mu`` derivative -scale. Thus the derivative-level perturbations are order unity: +Ballooning modes have long wavelength along the field and short wavelength across it. +A perturbation ``\varphi`` periodic in the poloidal angle ``\theta`` and the Clebsch angle ``\beta = \zeta - q\theta + \int q' \theta_0\, d\psi`` is represented on the covering space ``-\infty < \theta < \infty``, ```math -P_\psi = P_\psi^{(0)} + P_s^{(1)}, -\qquad -q_\psi = q_\psi^{(0)} + q_s^{(1)}. +\varphi(\psi, \theta, \beta) = \sum_{l=-\infty}^{\infty} \hat{\varphi}(\psi, \theta - l)\, \exp[-2\pi i n(\beta + lq)], ``` -The subscript ``s`` therefore represents the promoted local derivative -perturbation used as the practical scan variable in the ``s``-``\alpha`` diagram. -Retained first-order contributions are denoted by a superscript ``[1]``; for -example, ``I_{\mathrm{per}}^{[1]}`` is the retained local periodic variation of -the integrated shear. - -## Local ``P'`` and ``q'`` perturbations - -Using the integrated-shear function +so that the secular growth of the perpendicular wave number with the shear ``q'`` enters through ``|\nabla\beta|^2``. +At marginal stability the ideal energy principle gives the ballooning equation along the field line, ```math -I_{\mathrm{tot}} = q'(\theta-\theta_k) + I_{\mathrm{per}}, +\mathbf{B}\cdot\nabla\left(\frac{|\nabla\beta|^2}{B^2}\, \mathbf{B}\cdot\nabla\hat{\varphi}\right) + 2\frac{P'}{\chi'}\kappa_w \hat{\varphi} = 0 , ``` -the wave-number space and curvature are written - -```math -\left(\frac{|\nabla\beta|}{B}\right)^2 -= -\frac{1}{\chi'^2|\nabla\psi|^2} -+ \frac{|\nabla\psi|^2}{B^2}I_{\mathrm{tot}}^2, -``` +with ``\chi`` the poloidal flux, ``P'`` the pressure gradient and the projected curvature ```math -\kappa_w -= -\frac{\boldsymbol{\kappa}\cdot\nabla\psi} -{\chi'|\nabla\psi|^2} -- \kappa_s I_{\mathrm{tot}}. -``` - -With - -```math -I_{\mathrm{tot}} -= I_{\mathrm{tot}}^{(0)} + I_{\mathrm{tot}}^{[1]}, +\kappa_w \equiv \frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\beta = \kappa_n - q'\theta\,\kappa_s , \qquad -I_{\mathrm{tot}}^{[1]} -= q_s^{(1)}(\theta-\theta_k) + I_{\mathrm{per}}^{[1]}, -``` - -the retained order-unity perturbation terms are - -```math -\left(\frac{|\nabla\beta|^2}{B^2}\right)_{O(1),\mathrm{ret}} -= -\frac{|\nabla\psi|_0^2}{B_0^2} -\left[ - 2I_{\mathrm{tot}}^{(0)}I_{\mathrm{tot}}^{[1]} - + \left(I_{\mathrm{tot}}^{[1]}\right)^2 -\right], -``` - -```math -\kappa_w^{[1]} -= --\kappa_s^{(0)}I_{\mathrm{tot}}^{[1]}. +\kappa_s = \frac{\boldsymbol{\kappa}\times\mathbf{B}}{B^2}\cdot\nabla\psi = -\frac{\pi f}{\mathcal{J}}\frac{\partial}{\partial\theta}\left(\frac{1}{B^2}\right), ``` -The cross term and nonlinear ``(I_{\mathrm{tot}}^{[1]})^2`` term are retained in -the matrix assembly rather than discarded by linearization. - -## ``\mathbf{N}_0`` matrix construction +split into its non-secular part ``\kappa_n`` and the secular part driven by the shear. +The code integrates it as the first-order system ``\partial_\theta \mathbf{u} = \mathbf{L}\mathbf{u}`` in ``\mathbf{u} = (\hat{\varphi},\ |\nabla\beta|^2 \chi' \partial_\theta\hat{\varphi}/(B^2\mathcal{J}))``. -The lower-triangular element of the canonical transformation combines the local -``P'`` and ``q'`` perturbations: +As ``\theta \to \infty`` the coefficients diverge like ``q'^2\theta^2|\nabla\psi|^2``. +The canonical transformation ``\mathbf{R} = \operatorname{diag}(\theta^{-1/2}, \theta^{1/2})`` followed by a shear ``\mathbf{S}`` with ``\sigma = -2\pi f P' q'/(\chi'^2 B^2)`` makes the leading-order matrix ``\mathbf{N}_0`` finite, and the period average of that matrix fixes the asymptotic exponent: ```math -\sigma -= --\frac{2\pi f}{\chi_0'^2B_0^2} -\left(P_\psi^{(0)} + P_s^{(1)}\right) -\left(q_\psi^{(0)} + q_s^{(1)}\right). +\det[\langle\mathbf{N}_0\rangle - \alpha\mathbf{I}] = 0, \qquad \alpha = \pm\sqrt{-D_I} . ``` -The ``(1,2)`` component of the asymptotic matrix responds to ``q_s^{(1)}``: +This is the Mercier index: ``D_I = -\alpha^2`` evaluated from ``\langle\mathbf{N}_0\rangle`` equals the classical surface-average Mercier expression, and ``D_I > 0`` (complex exponents) is interchange instability. +The large and small solutions behave as ``\theta^{\alpha - 1/2}`` and ``\theta^{-\alpha - 1/2}``; the boundary condition at the finite integration limit ``\theta_\max`` includes the first-order correction of the asymptotic series, and the ballooning index is the jump in the logarithmic derivative of the solution between the two ends. -```math -(\mathbf{M}_0)_{12} -= -\frac{\mathcal{J}_0}{\chi_0'} -\frac{B_0^2} -{\left(q_\psi^{(0)}+q_s^{(1)}\right)^2|\nabla\psi|_0^2}. -``` +The s–α diagram scans the local stability of one surface without recomputing the equilibrium: a thin-layer ordering promotes local perturbations of the pressure gradient ``P'_s`` and shear ``q'_s`` to order unity in the derivatives, and the integrated shear ``I_\mathrm{tot} = q'(\theta - \theta_k) + I_\mathrm{per}`` carries them into ``|\nabla\beta|^2`` and ``\kappa_w``. +The cross term and the quadratic term in the perturbed shear are kept in the matrix assembly, and the ballooning angle ``\theta_k`` is the gauge where the total accumulated shear vanishes, ``I_\mathrm{tot}(\theta_k) = 0``. -The ``(2,1)`` component uses the expanded non-secular curvature, +## Running it -```math -\kappa_n(\theta) -= -\kappa_n^{(0)}(\theta) -+ \kappa_s^{(0)}(\theta)q_s^{(1)}\theta_k -- \kappa_s^{(0)}(\theta)I_{\mathrm{per}}^{[1]}(\theta), -``` +A `[LocalStability]` section runs the stage; it needs `[ForceFreeStates]` in the same deck. +Its keys are the problem settings, and the optional `[Ballooning]` section holds the [`Ballooning`](@ref) numerical settings, as in `examples/DIIID-like_ideal_example/gpec.toml`: -which gives +```toml +[LocalStability] +theta_k = 0.0 # ballooning angle +alpha_boundary_scan = true # also find the critical alpha of each surface -```math -(\mathbf{M}_0)_{21} -= --\frac{2\mathcal{J}_0}{\chi_0'^2} -\left(P_\psi^{(0)} + P_s^{(1)}\right) -\left[ - \kappa_n^{(0)}(\theta) - + \kappa_s^{(0)}(\theta)q_s^{(1)}\theta_k - - \kappa_s^{(0)}(\theta)I_{\mathrm{per}}^{[1]}(\theta) -\right]. +[Ballooning] +ode_tol = 1e-8 +n_scan = 24 # alpha scan points +max_alpha_scale = 8.0 +alpha_tol = 1e-3 +psi_window = [0.1, 0.99] # surfaces scanned for the alpha boundary ``` -The final matrix -``\mathbf{N}_0(\theta)=\mathbf{S}^{-1}\mathbf{M}_0(\theta)\mathbf{S}`` is -period-averaged to form ``\langle\mathbf{N}_0\rangle`` in the Mercier solvability -equation. +A script runs the same stage with `solve(LocalStabilityProblem(eq; theta_k=0.0, alpha_boundary_scan=true), Ballooning(); verbose=false)`; see the [Scripting API](api.md). +`examples/DIIID-like_ideal_example/analyze_example.jl` plots the ``s`` and ``\alpha`` profiles, the ``D_I`` and ballooning ``\Delta'`` profiles, and the two-dimensional s–α maps with their zero contours, using `salpha_reference`, `compute_ballooning_stability!` and `scan_delta_prime_map`. -## ``\mathbf{N}_1`` correction +## Outputs -The first-order asymptotic boundary condition requires the -``O(\theta^{-2})`` matrix ``\mathbf{N}_1``. In the large-``\theta`` limit, only the -upper-right component of ``\mathbf{M}_1`` remains: - -```math -\mathbf{M}_1(\theta) -= -\begin{pmatrix} -0 & a_1(\theta) \\ -0 & 0 -\end{pmatrix}, -``` - -where - -```math -a_1(\theta) -= --2\frac{\mathcal{J}_0}{\chi_0'} -\frac{B_0^2} -{|\nabla\psi|_0^2\left(q_\psi^{(0)}+q_s^{(1)}\right)^3} -\left[ - I_{\mathrm{per}}^{(0)}(\theta) - + I_{\mathrm{per}}^{[1]}(\theta) - - \left(q_\psi^{(0)}+q_s^{(1)}\right)\theta_k -\right]. -``` - -Applying the ``\mathbf{S}^{-1}`` transformation gives - -```math -\mathbf{N}_1(\theta) -= -a_1(\theta) -\begin{pmatrix} -\sigma & 1 \\ --\sigma^2 & -\sigma -\end{pmatrix}. -``` - -Thus ``\mathbf{N}_1`` is closed in terms of -``I_{\mathrm{per}}^{[1]}``, ``P_s^{(1)}``, and ``q_s^{(1)}``. The -``I_{\mathrm{per}}^{[1]}`` and ``q_s^{(1)}`` terms control ``a_1(\theta)``, while -``P_s^{(1)}`` enters through ``\sigma``. - -## Practical implementation in `Bal.jl` - -### Geometric curvature formula - -The non-secular curvature calculation uses the geometric expression - -```math -\kappa_n^{(0)} -= --\frac{1}{\mathcal{J}_0B_0} -\frac{\partial}{\partial\psi} -\left(\frac{\mathcal{J}_0B_0}{\chi'_0}\right) -+ -\frac{1}{\mathcal{J}_0B_0} -\frac{\partial}{\partial\theta} -\left( - \mathcal{J}_0 - \frac{\mathbf{B}_0\cdot\nabla\theta\times\nabla\zeta}{B_0} -\right) -+ -\frac{2\pi f q'_0}{B_0^2\mathcal{J}_0}. -``` - -This avoids relying entirely on Grad-Shafranov identities, which can amplify -numerical equilibrium errors near the edge or in other regions where the input -equilibrium does not satisfy the Grad-Shafranov equation exactly. - -### Two-loop structure - -The code separates the calculation into two loops to avoid repeatedly -differentiating noisy geometric quantities. - -- **Loop 1:** Use only unperturbed background quantities, such as ``P'`` and - ``q'``, to precompute the base structural curvature ``\kappa_n^{(0)}``. -- **Loop 2:** Insert the requested local perturbations, such as - ``P_s^{(1)}``, ``q_s^{(1)}``, and ``I_{\mathrm{per}}^{[1]}``. The code adds only - the analytically derived algebraic perturbation - ``\kappa_s^{(0)}q_s^{(1)}\theta_k - - \kappa_s^{(0)}I_{\mathrm{per}}^{[1]}`` to the precomputed - ``\kappa_n^{(0)}``. - -This preserves the perturbation structure while avoiding unstable repeated -geometric differentiation. - -### ``\theta_{\max}`` and `Baloo.f`-style matching - -For large ``\theta``, the wave-number factor behaves like a quadratic polynomial, - -```math -a_2\theta^2 + a_1\theta + a_0 -= -a_2\left(\theta+\frac{a_1}{2a_2}\right)^2 -+ \left(a_0-\frac{a_1^2}{4a_2}\right). -``` - -The residual scale is ``c = a_0-a_1^2/(4a_2)``. To make the asymptotic -approximation valid, the code chooses an integration boundary where the -``\theta^2`` term dominates, using the scale ``\sqrt{|c/a_2|}``. -Because the shooting boundary is a single flux-surface quantity, `Bal.jl` -evaluates this scale on the ballooning theta grid and uses the largest finite -value on that surface. The base integration boundary is then capped as -``\theta_{\max}=\min(16.5,10\max_\theta\sqrt{|c/a_2|})``. - -Directly extracting ``c_{a1}`` or ``\Delta'`` from the small asymptotic solution is -numerically fragile. The comparison below shows that the `Baloo.f` style, which -sets the displacement to vanish at ``\pm\theta_{\max}`` and computes ``\Delta'`` -from the resulting solution, gives the most stable continuous profile. - -![Comparison of numerical stability among methods for extracting ``\Delta'`` and ``c_{a1}``.](figures/ballooning/delta_prime_matching.png) - -`Bal.jl` therefore adopts the `Baloo.f`-style boundary replacement for the -ballooning ``\Delta'`` calculation. - -### Ballooning angle and integral-gauge correction - -The ballooning angle ``\theta_k`` is the integration constant specifying where the -total accumulated local shear is zero: - -```math -I_{\mathrm{tot}}(\theta_k) = 0. -``` - -In `Bal.jl`, - -```math -I_{\mathrm{tot}} -= -\left(q'^{(0)}+q_s^{(1)}\right)(\theta-\theta_k) -+ \left(I_{\mathrm{per}}^{(0)}(\theta) - + I_{\mathrm{per}}^{[1]}(\theta)\right). -``` - -To satisfy the gauge exactly, the periodic shear functions are shifted so that -they vanish at ``\theta_k``: - -```math -I_{\mathrm{per}}^{(0)}(\theta) -\leftarrow -I_{\mathrm{per}}^{(0)}(\theta) -- I_{\mathrm{per}}^{(0)}(\theta_k), -``` - -```math -I_{\mathrm{per}}^{[1]}(\theta) -\leftarrow -I_{\mathrm{per}}^{[1]}(\theta) -- I_{\mathrm{per}}^{[1]}(\theta_k). -``` +| Dataset | Meaning | +|---|---| +| `LocalStability/psi` | the flux surfaces, ψ_N | +| `LocalStability/D_I` | Mercier index; ``D_I > 0`` is interchange unstable | +| `LocalStability/D_R` | resistive interchange index; ``D_R > 0`` is unstable | +| `LocalStability/ballooning_Delta_prime` | ballooning stability index along `psi`; ``\Delta' = 0`` is marginal | +| `LocalStability/ballooning_psi`, `alpha`, `alpha_critical` | with `alpha_boundary_scan = true`: the scanned surfaces, their normalized pressure gradient and its critical value | -This aligns the numerical implementation with the theoretical gauge used in the -ballooning formulation. +`SingularSurfaces/D_I` repeats the Mercier index at the rational surfaces of the global solve. -### `locstab_fs` outputs +## Numerical notes -`Bal.jl` now computes both the Mercier criterion and the ballooning ``\Delta'``. -The previous standalone `Mercier.jl` path can therefore be removed from the local -stability calculation. The comparison below shows agreement between the previous -Mercier calculation and the new `Bal.jl` calculation. +- The non-secular curvature ``\kappa_n`` is evaluated from the geometric expression in ``\mathcal{J}B`` and its ``\psi`` and ``\theta`` derivatives rather than from Grad–Shafranov identities, which amplify equilibrium errors near the edge. +- Background quantities are differentiated once, in a first loop; the s–α perturbations are then added algebraically, so noisy geometric derivatives are never repeated. +- The integration limit is set where the quadratic term of the wave-number factor dominates, ``\theta_\max = \min(16.5,\ 10\max_\theta\sqrt{|c/a_2|})``, and the ballooning index is taken from a solution forced to vanish at ``\pm\theta_\max`` rather than from the coefficient of the small solution, which is numerically fragile: -![Comparison of Mercier criterion calculations from `Mercier.jl` and `Bal.jl`.](figures/ballooning/mercier_comparison.png) +![Comparison of the numerical stability of the methods for extracting the ballooning index.](figures/ballooning/delta_prime_matching.png) -The local-stability output now stores ballooning ``\Delta'`` in the fourth -`locstab_fs` entry. In the HDF5 output this is written as -`LocalStability/ballooning_Delta_prime`, distinct from the tearing ``\Delta'`` outputs -under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/`. +- The periodic shear functions are shifted to vanish at ``\theta_k`` so the gauge ``I_\mathrm{tot}(\theta_k) = 0`` holds exactly. -## API +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.LocalStability] diff --git a/docs/src/citations.md b/docs/src/citations.md index 9bf7f5a03e..7f8ac07b46 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -110,6 +110,22 @@ Establishes the area weightings that make flux-surface Fourier spectra of pertur --- +## LocalStability Module + +> A. H. Glasser, J. M. Greene, and J. L. Johnson, "Resistive instabilities in general toroidal plasma configurations," +> *Physics of Fluids* **18**, 875 (1975). + +The resistive interchange criterion ``D_R`` and the layer coefficients E, F, G, H, K, M of each rational surface, evaluated by `LocalStability` and stored under `SingularSurfaces/`. + +--- + +> J. W. Connor, R. J. Hastie, and J. B. Taylor, "Shear, periodicity, and plasma ballooning modes," +> *Physical Review Letters* **40**, 396 (1978). + +The ballooning representation on the covering space and the marginal ballooning equation that `LocalStability` integrates along the field line. + +--- + ## InnerLayer Module > A. H. Glasser, Z. R. Wang, and J.-K. Park, "Computation of resistive instabilities by matched asymptotic expansions," diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md index 05ecabf956..9fc7d97af0 100644 --- a/docs/src/galerkin.md +++ b/docs/src/galerkin.md @@ -1,57 +1,43 @@ -# Galerkin Δ′ Solver (RDCON outer region) - -The Galerkin solver is the Julia port of RDCON's `gal.f` — Dewar's singular Galerkin method for -the outer-region Δ′ computation. Instead of integrating the Euler-Lagrange ODE radially, it -discretizes the displacement on a packed radial grid of Hermite-cubic finite elements and solves -a single global banded system. Cells adjacent to each rational surface ("resonant" and -"extension" cells) carry the large/small power-law asymptotic solutions analytically, so the -singular behavior is built into the basis rather than resolved numerically. - -The ideal coefficient is reconstructed as ``F = (QL)(QL)^\dagger`` from the interpolated -Cholesky factor ``L``, with ``Q = \operatorname{diag}(m-nq)``. This uses the same coefficient -representation as the Euler–Lagrange integration and singular-surface expansions, including -between spline knots, and preserves Hermitian positive semidefiniteness of ``F``. - -The solve produces the inter-surface Δ′ matrix and the PEST-3 matching blocks -(``A'``, ``B'``, ``\Gamma'``, ``\Delta'``), written to the HDF5 output under the same -`SingularSurfaces/` paths the Riccati Δ′ uses (`Delta_prime_matrix`, `Delta_prime_raw`, -`Delta_coil`, `pest3_A`, `pest3_B`, `pest3_Gamma`) — one set of names per quantity, -whichever formalism produced it. These are the outer-region inputs to resistive matched-asymptotic stability analysis -[Glasser 2016, Phys. Plasmas **23**, 072505]. - -Select it with `integrator = "Galerkin"` in `[ForceFreeStates]`; its settings go in `[Galerkin]`. It replaces the radial ODE -integration rather than supplementing it: the run computes its own vacuum response at the -control surface (when `vac_flag`) and produces no free-boundary energies or ODE trace. -A `[Match]` section additionally matches the inner layer, giving a driven ξ solution. The -composite inner-region profiles need the cut solution, `[Galerkin] cut_solution = true` (or -`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(...))`). - -The banded solve automatically tries a Hermitian Cholesky -factorization of the lower band (LAPACK `zpbtrf`/`zpbtrs`, with the Dirichlet rows' columns moved -to the right-hand side). The free block's off-diagonal entries must be Hermitian within -`1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj))`; the diagonal's imaginary part is deliberately -discarded. Unequal bandwidths, off-diagonal asymmetry, or failure of positive definiteness -trigger a reported fallback to Jacobi-scaled LU. The factorization actually used is -recorded as `ForceFreeStates/Solutions/GalerkinIntegration/solver` in `gpec.h5`. - -The implementation lives in `src/ForceFreeStates/Galerkin/`: - -| File | Content | -|------|---------| -| `GalerkinStructs.jl` | Cell, interval, and workspace types mirroring the Fortran derived types | -| `GalerkinGrid.jl` | Packed grid construction and local→global DOF mapping | -| `GalerkinAssembly.jl` | Element-level assembly: Hermite basis, Gauss-Lobatto stiffness, resonant and extension cells, boundary conditions | -| `GalerkinSolution.jl` | Reconstruct ξ(ψ) and analytic ξ′(ψ) on the gal-native grid | -| `GalerkinSolve.jl` | Top-level driver `galerkin_solve`, banded solve, Δ′ extraction, PEST-3 blocks, HDF5 output | - -## API Reference +# Galerkin Solver + +The Galerkin solver is the third `[ForceFreeStates] integrator`: the singular Galerkin method of the Fortran RDCON for the outer-region Δ′ of resistive stability ([Glasser, Wang and Park (2016)](citations.md#InnerLayer-Module)). +Instead of integrating the Euler–Lagrange system radially, it discretizes the displacement on a packed radial grid of Hermite-cubic finite elements and solves one global banded system. +The cells adjacent to each rational surface carry the large and small power-law asymptotic solutions analytically, so the singular behavior is built into the basis rather than resolved numerically. + +The ideal coefficient is reconstructed as ``F = (QL)(QL)^\dagger`` from the interpolated Cholesky factor ``L`` with the resonance factor ``Q = \operatorname{diag}(m - nq)``, the same representation the Euler–Lagrange integration and the singular-surface expansions use, so ``F`` stays Hermitian positive semidefinite between spline knots. + +## Running it + +Select it in `[ForceFreeStates]` and set its options in `[Galerkin]`, as in `examples/DIIID-like_gal_resistive_example`: + +```toml +[ForceFreeStates] +integrator = "Galerkin" + +[Galerkin] +rpec_flag = true # edge coil-response columns (Delta_coil), needed by [Match] and [Tearing] +cut_solution = false # true: the cut solution that composite inner-region profiles need +``` + +It replaces the radial integration rather than supplementing it: the run computes its own vacuum response at the control surface (with `vac_flag`) but produces no free-boundary energies, no ODE trace and no `crit` scan, and it does not support the kinetic solve. +A `[Match]` section matches the inner layer to its Δ′ and gives a driven ``\xi`` solution; the composite inner-region profiles need `cut_solution = true` (`Galerkin(; cut_solution=true)` followed by `solve(MatchProblem(ffs))` in a script). + +## Outputs + +| Dataset | Meaning | +|---|---| +| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, on the same path as the Riccati result, so downstream stages never branch on which formalism ran | +| `SingularSurfaces/Delta_prime_raw`, `Delta_coil` | the side-major matching matrix and the coil-response columns (`rpec_flag = true`) | +| `SingularSurfaces/pest3_A`, `pest3_B`, `pest3_Gamma` | the PEST-3 interchange and mixed blocks ``A'``, ``B'``, ``\Gamma'`` | +| `ForceFreeStates/Solutions/GalerkinIntegration/` | the closed ``\xi`` profiles on the Galerkin grid, the surface list, the factorization used (`solver`: Cholesky, or an LU fallback when the band is not Hermitian positive definite) and, with `basis_output = true`, the `Basis/` group | + +Nothing is written under `ForceFreeStates/FreeBoundaryStability/`. + +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinSolve.jl"] ``` -## See also - -- `docs/src/stability.md` — the ForceFreeStates module this solver belongs to (EL system, metric matrices, singular surfaces) -- `docs/src/inner_layer.md` — inner-layer solutions to be matched against the outer-region Δ′ +See [Stability](stability.md) for the module this solver belongs to and [Tearing](inner_layer.md) for the inner-layer models matched to its Δ′. diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index dc8d93fa90..1b91bf9810 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -1,4 +1,4 @@ -# Tearing Module +# Tearing The `Tearing` module groups the resistive tearing-mode analysis stack: `InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and @@ -8,6 +8,24 @@ error-field penetration threshold), and the tearing solve itself, `solve(TearingProblem(ffs; ...), model; config=TearingConfig(...))`, with its `TearingResult` and HDF5 output. +## SLAYER inputs: perpendicular diffusivities + +The layer widths and the critical island width of `SLAYER` take perpendicular diffusivities from the kinetic file (`chi_e`, `chi_i`, `chi_phi` profiles in m²/s). +Two formulas want different quantities: the layer width goes through `P_perp`, the diffusion of the total pressure, while the critical island width ``W_d \propto (\chi_\perp/\chi_\parallel)^{1/4}`` is set by the electron heat transport that flattens the electron temperature across the island. +`[SLAYER] P_perp_model` selects how `P_perp` is built: from the perpendicular particle diffusivity `D_⊥` that weights the electron and ion channels by their share of the pressure gradient (`"D_perp"`, the default), from `chi_e` alone (`"chi_perp_e"`, the single-χ convention of Fitzpatrick's TJ code), from `chi_i`, from the toroidal Prandtl number, from its ``c_\beta^2`` floor, or from a whole-plasma confinement time (`"tau_E"`, for a plasma without χ profiles). +`[SLAYER] chi_perp_Dc` sets the electron diffusivity of the critical island width alone, one value for every surface, for a run in which one χ⊥,e cannot serve both formulas; where an interpretive power balance cannot separate the species, treat it as a modeling choice and scan it. +Since the widths scale as ``\chi^{1/4}``, a factor-two change in χ moves them by about 19%, which can be enough to move the resistive-layer overlap cut of [Stability](stability.md) by one surface. + +The DIII-D-like kinetic file of the examples carries analytic profiles, + +``` +chi_e = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²) +chi_i = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²) +chi_phi = 3/2 − (5/4) (1 − ψ)⁴ +``` + +which are of the order of power-balance values for DIII-D H-modes but are not fits to a measured discharge and do not resolve the pedestal transport barrier; `chi_e` is raised by a factor of two so that the default deck exercises the layer-overlap cut. + ## Tearing solve and layer inputs ```@autodocs diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index 016c70d471..a06ef21858 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/docs/src/perturbed_equilibrium.md @@ -1,9 +1,103 @@ # Perturbed Equilibrium -The `PerturbedEquilibrium` module computes the plasma response to external magnetic perturbations, -in two stages: `plasma_response(ffs)` builds the source-independent response and singular-coupling -operators of a force-free-states result once, and `perturbed_equilibrium(resp, rmp)` drives them with -one external field. A scan over fields reuses one response. +`PerturbedEquilibrium` computes the response of the plasma to an external non-axisymmetric field: the three-dimensional perturbed equilibrium that the force-free states and the applied field together form. +It is the stage that turns a stability calculation into the quantities of error-field and RMP physics: the field that actually reaches each rational surface, the islands it would open, and the energy and torque the perturbation carries. + +## What it computes + +The force-free states of [Stability](stability.md) give the displacement inside the plasma for any boundary value, and the [Vacuum](vacuum.md) response gives the field outside it. +On the control surface at `psilim` the two are joined by requiring continuity of the normal field, which fixes the plasma response to an applied spectrum uniquely ([Park et al. (2007a)](citations.md#PerturbedEquilibrium-Module), [Park et al. (2009)](citations.md#PerturbedEquilibrium-Module)). +`plasma_response(ffs)` builds the operators of that response once, with no reference to a drive; `perturbed_equilibrium(resp, rmp)` applies them to one external field from [Forcing Terms](forcing_terms.md), so a scan over fields reuses one response. +`compute_perturbed_equilibrium` is the lower-level form of the same call. + +The resonant part of the response is reported at every rational surface as the singular coupling: the pitch-resonant field and current, the island half-width and the Chirikov overlap parameter. +In an ideal solve the plasma shields every rational surface perfectly, so the stored resonant field is the effective one, the field the ideal response would deliver to the surface before shielding; a kinetic solve (`kinetic = "KineticForces"`) gives a finite penetrated field. See the [Conventions](conventions.md) page. + +## Governing equations + +The applied field is specified as its normal-flux spectrum ``\Phi_x`` on the control surface, one entry per (m, n). +The plasma inductance ``\Lambda`` relates the total flux to the surface current the plasma response induces, and the surface inductance ``L`` is the vacuum inductance of the same surface; the plasma response is the permeability + +```math +\Phi = P\,\Phi_x, \qquad P = \Lambda L^{-1}, +``` + + +oindent and the reluctance `` arrho = L^{-1}(\Lambda^\dagger - L)L^{-1}`` gives the surface current per unit applied flux. +The perturbed plasma energy and the toroidal torque of the response are quadratic forms in the total flux, + +```math +\delta W_\mathrm{plasma} = frac14\,\Phi^\dagger \Lambda^{-1}\Phi, \qquad +T_\phi = - frac{2n}{4}\,\operatorname{Im}\!\left(\Phi^\dagger \Lambda^{-1}\Phi +ight), +``` + + +oindent with the vacuum and surface energies written alongside them. +The stored matrices act on the coordinate-invariant root-area-weighted field `` ilde b``; the area-weighted field ``ar b`` and the flux follow from ``ar b = S ilde b`` and ``\Phi = Aar b`` with the operator `S` and the control-surface area `A` stored next to them (see *Field amplitudes* on the [Conventions](conventions.md) page). + +At each rational surface q = m/n the singular coupling evaluates the resonant flux ``\Phi^r`` of the response and reports the area-weighted resonant field ``ar b^r = \Phi^r / A^r`` in tesla, the pitch-resonant current, and the island the resonant flux would open in vacuum superposition, with half-width + +```math +w_s \propto \sqrt{|\delta\psi_s|} +``` + + +oindent in normalized flux and the Chirikov parameter ``\sigma_s = w_s / d_s``, the half-width over the half-distance to the neighbouring rational surface, so ``\sigma > 1`` means the islands overlap. +The coupling matrices `C_*` map the applied `` ilde b`` to each of these quantities, and their singular-value decomposition ranks the applied spectra by resonant drive (see below). + +## Running it + +From `examples/DIIID-like_ideal_example/gpec.toml`, after the `[ForceFreeStates]` solve and a `[ForcingTerms]` drive: + +```toml +[PerturbedEquilibrium] +compute_singular_coupling = true # resonant fields, islands and Chirikov parameter +regularization_width = 0.05 # width of the regularization of the ideal response at each rational surface +``` + +`compute_torque_response = true` adds the torque response matrices of a kinetic forward solve (below). +In a script the same two steps are + +```julia +using GeneralizedPerturbedEquilibrium + +resp = plasma_response(ffs; compute_singular_coupling=true) +pe = perturbed_equilibrium(resp, rmp; regularization_width=0.05) # rmp: an RMPField, see the Scripting API +save_to_hdf5("gpec.h5", pe) +``` + +with the `RMPField` built from coils or a file as on the [Scripting API](api.md) page; computing the response does not write a file. + +## Outputs + +| Dataset | Meaning | Units | +|---|---|---| +| `PerturbedEquilibrium/mode_m`, `mode_n` | the (m, n) of each entry of the mode axis, m fastest, one block per n | | +| `PerturbedEquilibrium/ForcingModes/m`, `n`, `amplitude` | the applied spectrum | T | +| `PerturbedEquilibrium/forcing_b`, `forcing_b_root_area`, `forcing_b_area` | the applied control-surface spectrum as bare, root-area-weighted and area-weighted normal field | T | +| `PerturbedEquilibrium/response_b`, `response_b_root_area`, `response_b_area` | the total (applied plus plasma) spectrum in the same three forms | T | +| `PerturbedEquilibrium/ResponseMatrices/plasma_inductance`, `surface_inductance`, `permeability`, `reluctance` | `` ilde\Lambda``, `` ilde L``, `` ilde P``, `` ilde arrho`` in root-area-weighted field space | | +| `PerturbedEquilibrium/ResponseMatrices/rootarea_to_area_weight_operator`, `surface_area` | ``S`` and ``A`` of the control surface | m² (`surface_area`) | +| `PerturbedEquilibrium/Response/psi` | the ψ_N grid of the profiles | | +| `PerturbedEquilibrium/Response/xi_n`, `xi_R`, `xi_Z`, `xi_phi` | the displacement in mode space: normal and cylindrical components | m | +| `PerturbedEquilibrium/Response/xi_psi`, `Jxi_*`, `xi_cov_*`, `xi_clebsch_*` | its contravariant, Jacobian-weighted, covariant and Clebsch components | | +| `PerturbedEquilibrium/Response/b_n`, `b_R`, `b_Z`, `b_phi`, `b_psi_area_weighted` | the perturbed field in mode space: normal, cylindrical and area-weighted normal components | T | +| `PerturbedEquilibrium/Response/Jb_*`, `b_cov_*` | its Jacobian-weighted and covariant components | | +| `PerturbedEquilibrium/SingularCoupling/rational_psi`, `rational_q`, `rational_m`, `rational_n`, `rational_area` | the resonant (surface, n) rows and their surface areas | m² (`rational_area`) | +| `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field` | ``ar b^r`` at each rational surface | T | +| `PerturbedEquilibrium/SingularCoupling/penetrated_area_weighted_field` | the penetrated area-weighted field | T | +| `PerturbedEquilibrium/SingularCoupling/resonant_current` | the pitch-resonant current | A | +| `PerturbedEquilibrium/SingularCoupling/island_width_sq`, `island_half_width` | ``w_s^2`` and ``w_s`` in ψ_N | | +| `PerturbedEquilibrium/SingularCoupling/chirikov_parameter` | ``\sigma_s`` | | +| `PerturbedEquilibrium/SingularCoupling/Delta_prime` | the forcing-driven tearing Δ′ per surface (Riccati solves) | | +| `PerturbedEquilibrium/SingularCoupling/C_*` | the coupling matrices from the applied `` ilde b`` to each quantity above | | +| `PerturbedEquilibrium/SingularCoupling/DominantMode/` | the core-window SVD of the resonant-field coupling (below) | | +| `PerturbedEquilibrium/Energies/plasma_energy`, `vacuum_energy`, `surface_energy` | the perturbed energies | J | +| `PerturbedEquilibrium/Energies/toroidal_torque` | the boundary-response torque ``T_\phi`` | N·m | +| `PerturbedEquilibrium/TorqueResponse/` | the torque response matrices of a kinetic forward solve (below) | N·m/T² | + +The `Analysis` module plots these; see [Analysis](analysis.md). ## Types @@ -19,8 +113,6 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.plasma_response GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_perturbed_equilibrium ``` -Save a computed state with [`save_to_hdf5`](@ref); computing the response does not write a file. - ## Several toroidal modes With several toroidal modes (`nn = [1, 2, ...]`) every `mode` axis in `PerturbedEquilibrium/` runs over the full diff --git a/docs/src/stability.md b/docs/src/stability.md index ca8ccc1911..475a26f8df 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -1,494 +1,153 @@ -# Ideal MHD Stability (ForceFreeStates) +# Stability -The `ForceFreeStates` module implements ideal MHD stability analysis for axisymmetric toroidal -plasmas following the direct Newcomb criterion described in [Glasser 2016]. It solves the -Euler-Lagrange (EL) system derived from the potential energy functional, identifies singular -(rational) surfaces where resonant coupling occurs, and returns eigenmode energies, the -tearing stability parameters Δ', and the full inter-surface Δ' matrix. +`ForceFreeStates` determines the ideal MHD stability of the equilibrium and produces the force-free states, the displacement fields that every later stage builds on. +It implements the direct Newcomb criterion of the Fortran DCON ([Glasser (2016)](citations.md#ForceFreeStates-Module)): rather than solving an eigenvalue problem for the growth rate, it integrates the Euler–Lagrange equations of the perturbed potential energy across the plasma and reads stability from the solutions themselves. -## Physical background +## What it computes -Ideal MHD stability is determined by the sign of the perturbed potential energy +For a toroidal mode number n, the perturbation is expanded in poloidal harmonics m over a window around the resonant range, and the displacement ``\xi(\psi)`` becomes a vector of harmonic amplitudes. +The ideal MHD potential energy ``\delta W[\xi]`` is a quadratic functional of ``\xi``, and its extrema satisfy a linear system of ODEs in ``\psi`` whose coefficients are singular at every rational surface q = m/n, where the resonant harmonic decouples from the field-line bending. +Integrating that system from the axis outward with the prescribed crossing of each rational surface gives the force-free states: the set of displacements that satisfy force balance everywhere inside the plasma for arbitrary boundary values. -```math -\delta W[\xi] = \int_0^{\psi_\mathrm{lim}} \mathcal{F}(\xi, \xi') \, d\psi, -``` - -where ``\xi(\psi)`` is the poloidal displacement vector. The extremum of ``\delta W`` over all -admissible ``\xi`` satisfies the Euler-Lagrange system [Glasser 2016, Eq. 24]: - -```math -\frac{d}{d\psi} -\begin{pmatrix} U_1 \\ U_2 \end{pmatrix} -= -\begin{pmatrix} A & B \\ C & D \end{pmatrix} -\begin{pmatrix} U_1 \\ U_2 \end{pmatrix}, -\quad -A = -Q\bar{F}^{-1}\bar{K}, \; -B = Q\bar{F}^{-1}Q, \; -C = \bar{G} - \bar{K}^\dagger\bar{F}^{-1}\bar{K}, \; -D = \bar{K}^\dagger\bar{F}^{-1}Q, -``` - -where ``\bar{F}``, ``\bar{K}``, ``\bar{G}`` are the MHD metric matrices in Fourier-mode space -and ``Q = \mathrm{diag}(1/(m - nq))`` is the singular factor. The Newcomb criterion states -that the plasma is stable if and only if this system admits a regular solution that remains -finite across every rational surface. - -**Key references** - -| Paper | Content | -|-------|---------| -| [Glasser 2016] Phys. Plasmas **23**, 112506 | Newcomb criterion, EL system, standard DCON integration | -| [Glasser 2018a] Phys. Plasmas **25**, 032507 | Riccati reformulation, reduced stiffness near singular surfaces | -| [Glasser 2018b] Phys. Plasmas **25**, 032501 | STRIDE code: parallel FM integration, inter-surface Δ' matrix | - -## Integration methods - -Three formalisms are available, selected by `integrator`. Forward and Riccati solve the same -EL system and differ in numerical strategy and in what they leave behind for the rest of the -pipeline: the forward driver returns dense displacement profiles, the Riccati driver returns -the inter-surface ``\Delta'`` matrix. Galerkin solves the outer region variationally instead, -and is documented in `docs/src/galerkin.md`. +Three results follow from them: -### Forward integration +- **Fixed-boundary stability.** A zero of the determinant of the solution matrix between the axis and the edge means a displacement vanishing at the boundary with negative energy, so the plasma is unstable even with a perfectly conducting wall on its surface (the `crit` scan). +- **Free-boundary stability.** With `vac_flag = true`, the [Vacuum](vacuum.md) module supplies the vacuum energy of the `[Wall]` section (no wall when the section is absent), and the total energy ``W = W_\mathrm{plasma} + W_\mathrm{vacuum}`` on the control surface is a Hermitian matrix in the boundary harmonics. Its generalized eigenvalues with respect to the power norm are the free-boundary eigenmode energies; a negative least eigenvalue means an ideal external kink. +- **The tearing index Δ′.** Across each rational surface the solution has a large and a small asymptotic part. The Δ′ matrix of the outer region, assembled from the force-free states, gives the jump of the small solution driven by a unit tearing perturbation at each surface; it is the outer-region input of every resistive calculation ([Tearing](inner_layer.md), `[Match]`). -`forward_eulerlagrange_integration` is the baseline driver — our implementation of the -standard DCON radial integration [Glasser 2016]. It integrates the EL ODE directly -in ``(U_1, U_2)`` using the adaptive 9th-order `Vern9` solver. Near each rational surface the -columns of ``U_2`` that correspond to resonant modes are zeroed via Gaussian reduction (GR), -keeping the solution bounded; `transform_u!` undoes the reduction at the end, so `u_store` -comes back dense in the axis (Euler-Lagrange) basis. This is the reference path for -correctness comparisons, the only path whose solution the perturbed-equilibrium stage can -consume, and the only path that supports a kinetic source. +## Governing equations -Enable with: -```toml -[ForceFreeStates] -integrator = "Forward" -``` - -### Riccati integration - -`riccati_eulerlagrange_integration` (the default) is our implementation of the STRIDE -approach [Glasser 2018b], built on the dual Riccati reformulation [Glasser 2018a]. It -decomposes the radial domain into independent chunks, integrates each chunk's -fundamental-matrix (FM) propagator in parallel -using `Threads.@threads`, then multiplies the propagators in order and applies each -singular-surface crossing serially. It is the only driver that produces the inter-surface -``\Delta'`` matrix. Because chunk endpoints are all it stores, `u_store` is sparse and stays -in the Riccati basis — dense ``\xi`` profiles and the -`ForceFreeStates/Solutions/ForwardIntegration/xi_*` datasets require the forward driver. - -Crossings and the outer-plasma re-integration use the dual Riccati matrix -``S = U_1 \cdot U_2^{-1}`` [Glasser 2018a, Eq. 19]: +Ideal stability is decided by the sign of the perturbed potential energy ```math -\frac{dS}{d\psi} = w^\dagger \bar{F}^{-1} w - S\bar{G}S, \qquad -w = Q - \bar{K}S. +\delta W[\xi] = \int_0^{\psi_\mathrm{lim}} \mathcal{F}(\xi, \xi')\, d\psi , ``` -``S`` remains bounded near rational surfaces (where ``U_1, U_2`` grow exponentially). Rather -than integrating the quadratic Riccati ODE directly (which blows up when ``|S|`` is large), -the code integrates the linear EL system with `sing_der!` as the RHS and recovers -``S = U_1 U_2^{-1}`` via periodic renormalization — an approach that is mathematically -equivalent to O(Δψ) but uses `Vern9`'s full 9th-order accuracy. Renormalization is -triggered whenever ``\max(|U_1|)`` or ``\max(|U_2|)`` exceeds the threshold `renorm_threshold`, and is -forced at the end of each chunk. At singular surface crossings, -`riccati_cross_ideal_singular_surf!` applies the small-asymptotic matching directly in column -`ipert_res` — without Gaussian reduction — and renormalizes to ``(S, I)``. - -Enable with: -```toml -[ForceFreeStates] -integrator = "Riccati" +whose extrema over all admissible ``\xi`` satisfy a first-order system in the harmonic amplitudes ``U_1 = \xi`` and their conjugate momenta ``U_2`` (Glasser 2016, Eq. 24): -[Riccati] -nchunks = 0 # 0 = auto: derived from the singular-surface count alone +```math +\frac{d}{d\psi}\begin{pmatrix} U_1 \\ U_2 \end{pmatrix} += \mathbf{L}(\psi)\begin{pmatrix} U_1 \\ U_2 \end{pmatrix}, ``` -#### Chunking and thread independence - -`balance_integration_chunks` splits the base chunks until the count reaches a target derived -from the number of singular surfaces, `max(2 m_s + 3, 8(m_s + 1) + m_s)`, where ``m_s`` is -`msing`. Setting `nchunks` overrides that target; a value below the ``2 m_s + 3`` floor is -clamped up with a warning. The target never consults `Threads.nthreads()`, and every chunk -integrates independently from identity initial conditions, so the results do not depend on how -many threads `julia -t` provides — threads change wall-clock only. +where ``\mathbf{L}`` is built from the Fourier-space metric matrices ``\bar{F}``, ``\bar{K}`` and ``\bar{G}`` of the equilibrium and the resonance factor ``Q = \mathrm{diag}(m - nq)``. +The resonance factor vanishes at each rational surface, so the solutions there behave as ``|\psi - \psi_s|^{-1/2 \pm \sqrt{-D_I}}``, a large and a small solution governed by the Mercier index ``D_I``; GPEC assembles the right-hand side from reduced matrices that carry explicit ``(m - nq)`` factors rather than their inverses, the convention used on every page. +The Newcomb criterion states that the plasma is fixed-boundary stable if and only if the solution that is regular at the axis remains nonsingular across every rational surface and up to the edge. -#### Bidirectional integration for large N +The Riccati integrator works with the matrix ``S = U_1 U_2^{-1}`` ([Glasser, Kolemen and Glasser (2018)](citations.md#ForceFreeStates-Module), Eq. 19), which obeys a quadratic matrix ODE and stays bounded where ``U_1`` and ``U_2`` grow exponentially, so that the rational surfaces can be crossed without the Gaussian reduction of the forward integration. +Its global boundary-value problem over all surfaces gives the Δ′ matrix of [Glasser and Kolemen (2018)](citations.md#InnerLayer-Module), stored in the PEST-3 convention: the diagonal entry of each surface is its outer-region Δ′, positive when the outer region drives the tearing mode, and the off-diagonal entries couple the surfaces. -For large mode counts the FM propagator for a chunk ending near a rational surface is -ill-conditioned: the EL solutions grow exponentially toward the rational surface, so the -forward FM amplifies numerical errors. GPEC follows the STRIDE approach [Glasser 2018b, -Sec. III.A]: the crossing chunk (the last sub-chunk before each rational surface) is -integrated *backward* — from the rational surface toward the interior — producing a -well-conditioned backward FM ``\Phi_L``. The forward propagation is recovered as -``\Phi_L^{-1}`` via an LU solve in serial assembly, which is accurate precisely because -``\Phi_L`` is well-conditioned. +## Integrators -The implementation uses a `direction` field on `IntegrationChunk`: +`[ForceFreeStates] integrator` selects the radial method; its numerical settings sit in a section of the same name. -- `direction = +1`: standard forward integration, `tspan = (ψ_start, ψ_end)`. -- `direction = -1`: backward integration, `tspan = (ψ_end, ψ_start)` (reversed). +- `"Riccati"` (the default; `[Riccati]`) splits the domain into independent chunks integrated in parallel, crosses each rational surface in the Riccati basis and assembles the Δ′ matrix. It stores the chunk endpoints only, so it yields no dense ``\xi`` profiles. +- `"Forward"` (`[Forward]`) integrates the system directly in ``(U_1, U_2)`` with Gaussian reduction at each rational surface. It returns dense displacement profiles, is the only integrator the perturbed equilibrium stage consumes for ``\xi``-resolved output, and the only one that supports the kinetic solve. +- `"Galerkin"` (`[Galerkin]`) solves the outer region variationally on a finite-element basis, giving the resistive Δ′ matrix and the PEST-3 blocks but no free-boundary energies; see [Galerkin Solver](galerkin.md). -`chunk_el_integration_bounds(...; bidirectional=true)` assigns `direction = -1` to every -crossing chunk. `balance_integration_chunks` preserves this: the sub-chunk closest to the -rational surface inherits `direction`, while the earlier sub-chunk always gets `direction=+1`. +`kinetic = "KineticForces"` replaces the ideal energy matrices by their drift-kinetic form, which needs `[Equilibrium] kinetic_file`; see [Kinetic Forces](kinetic_forces.md). -**Accuracy** (N=26, DIIID-like example): energy eigenvalue within 2% of the forward path. -The residual ~2% gap comes from the different crossing convention (Riccati-style direct -zeroing vs GR), not from ODE tolerance; it is present at every thread count. +## Running it -### Galerkin +From `examples/DIIID-like_ideal_example/gpec.toml` (abridged), with the `[Riccati]` keys of `examples/DIIID-like_riccati_deltaprime_example`: -`integrator = "Galerkin"` solves the same EL system variationally instead of integrating it: -the outer region is discretized on packed Hermite-cubic elements and solved as one global banded -system, giving the RDCON resistive ``\Delta'`` matrix and the PEST-3 matching blocks. It computes its own vacuum response and -returns no free-boundary energies, no ODE trace, and no fixed-boundary `crit` scan. With -a `[Match]` section it also matches the inner layer, producing a driven ``\xi`` solution. -Kinetic runs are not supported. See -`docs/src/galerkin.md` for the solver and its `[Galerkin]` settings. - -Enable with: ```toml [ForceFreeStates] -integrator = "Galerkin" -``` - -## Local stability: Mercier and ballooning (s–α) - -A `[LocalStability]` section in the deck runs a local high-``n`` stability scan over every -flux surface, as its own stage before the global ideal analysis above. The section's keys are -the problem settings (`theta_k`, `alpha_boundary_scan`); an optional `[Ballooning]` section holds the -[`Ballooning`](@ref) numerical settings, and `[Run] verbose` the solve verbosity. A script runs the same stage with `solve(LocalStabilityProblem(eq), Ballooning(); verbose=false)` -(see the [scripting API](api.md)). -For the derivation and implementation details behind these diagnostics, see -[Ballooning and Mercier Local Stability](ballooning.md). -Three diagnostics are produced and stored under the `LocalStability/` HDF5 group, each a profile -in normalized poloidal flux ``\psi``: - -- **Mercier criterion ``D_I``** (`LocalStability/D_I`) — the ideal interchange criterion. A surface - is Mercier-unstable where ``D_I > 0``. It is evaluated from the ``\det(\bar{d}_0)`` of the - integrated local-mode matrix. -- **Resistive interchange ``D_R``** (`LocalStability/D_R`) — the Glasser–Greene–Johnson resistive - interchange criterion ``D_R = D_I + (H - 1/2)^2``. The ``D_I`` term is the same - ``\det(\bar{d}_0)`` value reported in `LocalStability/D_I`; ``H`` is computed from the legacy - Mercier/GGJ flux-surface averages of the field and metric quantities. ``D_R > 0`` - indicates resistive interchange instability. -- **Ballooning ``\Delta'``** (`LocalStability/ballooning_Delta_prime`) — the high-``n`` ballooning - stability index, obtained by integrating the ballooning equation along the field line and - taking the jump in the logarithmic derivative of the solution between the two asymptotic - ends. - -With `alpha_boundary_scan = true` the stage also writes the critical ballooning -``\alpha`` boundary: `LocalStability/ballooning_psi`, `LocalStability/alpha` and -`LocalStability/alpha_critical` (see the [scripting API](api.md)). - -!!! note "Two different Δ' quantities" - `LocalStability/ballooning_Delta_prime` is the **local high-``n`` ballooning** index and is - distinct from the **resistive tearing** ``\Delta'`` described in the next section, which - is written under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/Delta_prime`. - They measure different instabilities; do not confuse them. - -### s–α diagram - -The local criteria can be mapped over a two-parameter ``(p', q')`` scan at a single flux -surface — the classic s–α (shear–pressure) stability diagram. For a chosen surface the -pressure gradient ``p'`` and shear ``q'`` are scaled away from their equilibrium values, and -``\Delta'`` and ``D_I`` are recomputed on the resulting grid. The ``\Delta' = 0`` and -``D_I = 0`` contours bound the ballooning- and Mercier-stable regions, with the equilibrium -operating point marked inside. - -The example script `examples/DIIID-like_ideal_example/analyze_example.jl` demonstrates the -full workflow: it plots the ``s`` and ``\alpha`` profiles, the ``D_I`` and ballooning -``\Delta'`` profiles, and the 2-D s–α maps with their zero contours, using -`salpha_reference`, `compute_ballooning_stability!`, and `scan_delta_prime_map`. - -## Δ' tearing stability parameter - -### Per-surface Δ' (`delta_prime`) - -At each rational surface the asymptotic matching condition gives the tearing stability -parameter [Glasser 2016]: - -```math -\Delta'_s = \frac{c_{a,r}[i_s,i_s,2] - c_{a,l}[i_s,i_s,2]}{4\pi^2 \psi_0}, -``` - -where ``c_{a,l}`` and ``c_{a,r}`` are the left and right asymptotic coefficients at surface -``s``, and ``i_s`` is the column index of the resonant mode. Positive ``\Delta' > 0`` -indicates a tearing-unstable surface. - -The Riccati and parallel FM paths populate `intr.sing[s].delta_prime` (a length-``n_\mathrm{res}`` -vector) inline during each crossing. A companion vector `delta_prime_col` (length N) stores -the coupling of all poloidal modes to the resonant mode at surface ``s``: - -```math -(\Delta'_\mathrm{col})_{j,i} = \frac{c_{a,r}[j,i_s,2] - c_{a,l}[j,i_s,2]}{4\pi^2 \psi_0}. -``` +nn = 1 # toroidal mode number, or consecutive numbers as a list ([1, 2, 3]) +vac_flag = true # free-boundary energies with the [Wall] vacuum response +qlow = 1.02 # inner edge of the integration domain +qhigh = 1e3 # outer bound of qlim (see Integration domain) +delta_mlow = 8 # extra poloidal modes below the resonant range +delta_mhigh = 8 # extra poloidal modes above it +integrator = "Forward" # "Riccati" is the default +truncate_at_rational_offset = true +rational_offset_fraction = 0.2 + +[Forward] +eulerlagrange_tolerance = 1e-10 -The diagonal element ``(\Delta'_\mathrm{col})_{i_s,i}`` equals `delta_prime[i]` exactly by -construction. - -### Inter-surface Δ' matrix (`delta_prime_matrix`) - -`compute_delta_prime_matrix!` assembles an ``m_\mathrm{sing} \times m_\mathrm{sing}`` -inter-surface tearing matrix following the STRIDE global BVP [Glasser 2018b, Sec. III.B]. -Internally, the solver builds a raw ``2 m_\mathrm{sing} \times 2 m_\mathrm{sing}`` matrix -whose rows/columns index the *left* and *right* inner-layer boundaries of every rational -surface; the stored PEST3-convention ``\Delta'`` is the four-term combination -``\text{dp\_raw}[2i, 2j] - \text{dp\_raw}[2i, 2j{-}1] - \text{dp\_raw}[2i{-}1, 2j] + \text{dp\_raw}[2i{-}1, 2j{-}1]`` -that folds the raw block into a per-surface response. The BVP unknowns are the plasma -state at the left and right inner-layer boundaries of every rational surface; the driving -terms are unit-amplitude asymptotic solutions at each boundary. The resulting matrix -encodes the full plasma response between all pairs of surfaces and is required for -resistive stability analysis of multi-surface configurations. - -The BVP is well-conditioned because it is formulated using the split ``(\Phi_R, \Phi_L)`` -propagator blocks from bidirectional integration rather than the monolithic forward product -``\Phi_L^{-1} \Phi_R`` (which is ill-conditioned for large N): - -```math -\Phi_R[j] \cdot x_R[j-1] - \Phi_L[j] \cdot x_L[j] = 0 -\quad \text{(junction at } \psi_m[j]\text{)}, +[Riccati] +eulerlagrange_tolerance = 1e-10 +save_interval = 1 +renorm_threshold = 1e3 # renormalize the Riccati basis above this solution norm ``` -where ``\Phi_R[j]`` is the forward FM product from ``\psi_{R,j-1}`` to the junction, and -``\Phi_L[j]`` is the backward crossing FM from ``\psi_{L,j}`` to the junction. - -The matrix is written to the HDF5 output under `SingularSurfaces/Delta_prime_matrix`. -The Galerkin integrator computes the same quantity in the same PEST-3 convention and -publishes it on the same path, so downstream consumers (SLAYER among them) never branch -on which formalism ran. - -## Integration domain truncation - -The Euler–Lagrange integration runs from the axis out to `psilim`, which is set from the -target safety factor `qlim`: - -- `qlim` starts at `min(qmax, qhigh)`, where `qmax` is q at the equilibrium's `psihigh`. -- With `truncate_at_rational_offset = true` (the default, for diverted equilibria where q → ∞ at the - separatrix), `qlim` moves to `(last rational + rational_offset_fraction) / n`, so the integration ends a fixed - fraction of a rational spacing past the last surface instead of close to one. The offset is - then stepped back by `1/n` until it no longer exceeds the bound it started from, so `rational_offset_fraction` - never extends the domain past `qmax`, `qhigh`, or the layer-overlap cap below. A bound that - falls on a rational, or less than `rational_offset_fraction / n` above one, therefore excludes it: with - `rational_offset_fraction = 0.2`, both `qhigh = 4.0` and `qhigh = 4.1` give `qlim = 3.2`. -- With `truncate_at_rational_offset = false` (limited or analytic equilibria with finite edge q), - `qhigh` and `psihigh` set the domain directly. +The same solve from a script takes the deck keys as keywords of the problem and the integrator: -### Resistive-layer overlap - -Near the separatrix the rational surfaces crowd together faster than their resistive layers -narrow, so eventually adjacent layers overlap. Past that point no surface keeps a -well-separated inner region and the matched-asymptotic treatment is not defined; Sect. 5.9 of -Fitzpatrick (2025), listed under [InnerLayer Module](@ref) in the citations, restricts the -response calculation to the region inside it. - -GPEC locates that point with a scan of the rational surfaces of each toroidal mode number in -the run; for multi-`n` runs the innermost overlap point, over all `n`, is the one used. Each surface gets a layer solve from the kinetic profiles, and its -diffusive-resistive width (Fitzpatrick 2025, Eq. 100) is converted to a normalized-flux width -so it can be compared with the spacing between neighbouring surfaces. The cut is placed at the -last surface before two adjacent layers first touch. - -Surfaces beyond the equilibrium grid are located on the logarithmic divergence of q at the -separatrix, taken from the equilibrium's own poloidal flux rather than fitted. Each X-point on -the last closed flux surface contributes one term, - -```math -q(\psi) = q_K + \sum_i A_i \ln\frac{\psi_{X_i} - \psi_K}{\psi_{X_i} - \psi}, -\qquad A_i = \frac{F}{2\pi R_{X_i}\sqrt{-\det H_i}}, -``` +```julia +using GeneralizedPerturbedEquilibrium -where ``\psi_{X_i}`` is the flux at the saddle, ``H_i`` the Hessian of the poloidal flux there, -``F = R B_\phi``, and ``q_K`` the traced q on the outermost grid surface ``\psi_K``. One term -describes a single null and two a double null, connected or not. The X-points are found by -following each minimum of the poloidal field along the outermost surface to the saddle it leads -to; on a single-null shape this usually also reaches the far secondary X-point, whose term is -small. The law is used only if it reproduces the traced shear over the outer knots, and an -equilibrium with no X-point nearby is not extrapolated. A limited plasma that has an X-point -nearby is not told apart from a diverted one: one rational is placed past its boundary and the run -warns that no overlap could be established there, leaving the domain unchanged. Divertors with a -second null beside the primary one or on the leg (snowflake, x-point target) are not modelled. - -The X-point flux need not coincide with the equilibrium's boundary: a reconstruction that places -its boundary slightly inside the separatrix has finite q at ``\psi = 1``, and rational surfaces -between that boundary and the X-point lie in the real plasma but outside the equilibrium GPEC -holds. The scan places only the first of them, to judge whether the outermost surface inside is -clean, and scores it only when the kinetic profiles extend that far. `psihigh` remains the upper -limit throughout: the overlap point can only narrow the domain, and when no overlap can be -established the domain keeps its requested limit. - -TokaMaker's `save_eqdsk` writes such files when `lcfs_pad` (``p``) is set with the default -`truncate_eq=True`: the equilibrium is cut at ``1 - p`` and renormalised, so the X-point sits at -``\psi = 1/(1 - p)``. Writing with `truncate_eq=False` keeps the whole plasma; EFIT-written files -place their boundary on the X-point flux. - -The scan runs only when the cap is switched on, and it needs kinetic profiles on the equilibrium. -It reads the plasma from the equilibrium's kinetic data and the layer closures from a -`SLAYERModel`, the same inputs the SLAYER tearing solve uses, so the layer widths cannot disagree -between them. In a script, `psilim_from_layer_overlap = true` scans with the default closures and a -model sets others: +eq = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada") +ffs = solve(eq, Riccati(; eulerlagrange_tolerance=1e-10); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true, qlow=1.02) -```julia -ffs = solve(eq, Riccati(); nn=1, psilim_from_layer_overlap=SLAYERModel(; P_perp_model=:chi_perp_e)) -ffs.layer_overlap.psihigh # the overlap point +real(ffs.free_boundary.et[1]) # least free-boundary eigenmode energy: negative is unstable +for (s, sing) in enumerate(ffs.surfaces) + println("q = $(sing.m[1])/$(sing.n[1]) at ψ_N = $(sing.psifac): Δ′ = $(real(ffs.delta_prime.matrix[s, s]))") +end ``` -In a deck, `psilim_from_layer_overlap = true` takes its closures from the `[SLAYER]` section. The -scan is carried on the result as `layer_overlap` and written to `ForceFreeStates/LayerOverlap/`. It -caps `qlim` before the `rational_offset_fraction` step, so `rational_offset_fraction` still selects -the final surface from inside the cap. A cap beyond `psihigh` has no effect. - -The flag is off by default, so every shipped deck is unchanged by it. With the default -`P_perp_model = "D_perp"` the layers on the shipped DIII-D-like profiles do not overlap inside -ψ = 1. With the electron-channel `"chi_perp_e"` they are about 40% wider and the overlap point is -ψ ≈ 0.9985: outside the ideal and Riccati decks' `psihigh = 0.995`, but inside the SLAYER deck's -0.9995, where the cap would exclude the q = 7 surface. `examples/DIIID-like_truncation_example` -pins `"chi_perp_e"` and raises `psihigh` to 0.9999 with the flag on: the cap excludes q = 7 and -q = 8, and the resulting domain (`psilim ≈ 0.99498`, last surface q = 6) reproduces the hand-chosen -`psihigh = 0.995` of the ideal deck. - -**Which perpendicular diffusivity.** The two formulas that take a perpendicular diffusivity want -different quantities. The layer width goes through `P_perp`, which multiplies the diffusion of the -total pressure; under the closure `T_s = T_s(n)` it is set by the perpendicular particle diffusivity -`D_⊥` of Fitzpatrick, Phys. Plasmas **29**, 032507 (2022), Eq. 16, which weights the electron and ion energy diffusivities by -`n T_s′/p′`. The critical-Δ island width `W_d ∝ (χ_⊥/χ_∥)^{1/4}` is the -width above which parallel transport, through an electron parallel conductivity, makes the electron -temperature a flux function (Fitzpatrick, Phys. Plasmas **32**, 062509 (2025), Eqs. 230–234). That -paper and the `"lar"` and `"toroidal"` forms of Connor et al., Plasma Phys. Control. Fusion **57**, -065001 (2015), Eqs. 59 and 61, are single-fluid: they have one χ_⊥, which in Connor et al. diffuses -the total pressure. GPEC reads it as the electron channel and uses `chi_e` for every `D_c_type`, -since the width is set by electron heat transport; that reading is a modeling choice, not a result of -those papers. `SLAYERModel`'s `chi_perp_Dc` (`[SLAYER] chi_perp_Dc`, optional) sets that electron diffusivity for the critical Δ -alone, leaving `P_perp` on its own χ inputs. The deck key is one value for every surface, over -any `chi_e` profile; `build_slayer_inputs` also takes a function of ψ. It moves the critical -resonant field through the critical Δ. It is for a run where one χ⊥,e cannot serve -both formulas. `D_⊥` needs only the total conducted heat, so a run can use a one-fluid χ_eff in both -channels when an interpretive power balance cannot separate the species; where T_i ≈ T_e the -electron-ion exchange term makes that split uncertain, and it can return χ⊥,i < 0. On the electron reading, -the critical Δ still needs an electron diffusivity, which χ_eff is not, and `chi_perp_Dc` makes the choice -explicit. Where the split is uncertain the electron χ is uncertain too, so treat the value as a -modeling choice and scan it. The same split can often be reached without this key, by choosing -`chi_e` and `chi_i` to keep the `D_⊥` total, but not when the wanted `chi_e` would leave `chi_i ≤ 0`. -`SLAYERModel`'s `P_perp_model` (`[SLAYER] P_perp_model` in a deck) selects how `P_perp` is built: from `D_⊥` (`"D_perp"`, -the default), from `chi_e` alone (the single-χ convention of Fitzpatrick's TJ code), from -`chi_i`, from the toroidal Prandtl number, from its `c_β²` floor, or from a whole-plasma -confinement time (`"tau_E"`); see `slayer_parameters`. `"tau_E"` is for a plasma with no χ -profiles: its uniform χ_E sets both `P_perp` and the toroidal `P_tor`, so the file's `chi_phi` -and the kinetic data's `chi_phi_default` are not used in that mode. On the DIII-D-like -profiles `D_⊥` is 0.22–0.32 of `chi_e` at the outer rational surfaces, because the weights sum to -less than one where density gradients carry part of `p′`. Since the widths scale as `χ^{1/4}`, a -factor-2 change in χ moves them by about 19%, which can be enough to move the overlap cut by a -surface. - -**The shipped diffusivities.** The DIII-D-like kinetic file carries analytic χ profiles in the -normalized poloidal flux ψ, in m²/s: +`main(["examples/DIIID-like_riccati_deltaprime_example"]).ffs` returns the same result from a deck. +The number of Julia threads is set at startup (`-t N` or `JULIA_NUM_THREADS`); the chunking never depends on it, so threads change the wall-clock time only. -``` -chi_e = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²) -chi_i = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²) -chi_phi = 3/2 − (5/4) (1 − ψ)⁴ -``` +## Outputs -They are not fits to a measured discharge. Their magnitudes are of the same order as the TRANSP -power-balance χ_e and χ_i reported for DIII-D H-modes by Hong *et al.*, Phys. Plasmas **30**, -072512 (2023), Fig. 5; `chi_e` includes a factor-2 raise so that the default deck exercises the -resistive-layer overlap cut. `chi_phi` is set to track `chi_i` (Prandtl number of order one, rising -outward), as reported for DIII-D by Yoshida *et al.*, Nucl. Fusion **52**, 123005 (2012). The forms -do not resolve the H-mode pedestal transport barrier, where power-balance χ falls to about -0.1–0.3 m²/s (Callen *et al.*, Nucl. Fusion **50**, 064004 (2010)). +| Dataset | Meaning | +|---|---| +| `ForceFreeStates/FreeBoundaryStability/W_plasma`, `W_vacuum`, `W_freeboundary` | plasma, vacuum and total energy matrices on the control surface, power-normalized (per unit ⟨\|ξ\|²⟩) | +| `ForceFreeStates/FreeBoundaryStability/eigenmode_energies` | generalized eigenvalues of ``W``, most unstable first; `eigenmode_plasma_energies` and `eigenmode_vacuum_energies` split each one | +| `ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes` | the eigenvectors, unit power norm | +| `ForceFreeStates/FreeBoundaryStability/vacuum_eigenvalue` | least eigenvalue of the vacuum energy matrix | +| `ForceFreeStates/EdgeScan/` | `total_energy`, `plasma_energy`, `vacuum_energy` and `vacuum_eigenvalue` against `psi` and `q` as the control surface moves outward from `dW_edge_scan_start` | +| `ForceFreeStates/Solutions/ForwardIntegration/` | the forward integrator's `psi`, `q`, `xi_psi`, `xi_s`, `dxi_psidpsi` profiles and the `crit` scan whose zero crossings mark fixed-boundary instability | +| `ForceFreeStates/EulerLagrangeMatrices/` | the ``\bar{F}``, ``\bar{K}``, ``\bar{G}`` matrices on `psi`, written when `write_el_matrices = true` | +| `ForceFreeStates/LayerOverlap/` | the resistive-layer overlap scan, when `psilim_from_layer_overlap` is set (table below) | +| `SingularSurfaces/rational_psi`, `rational_q`, `rational_m`, `rational_n`, `dqdpsi` | the rational surfaces of the domain, core to edge | +| `SingularSurfaces/D_I`, `E`, `F`, `G`, `H`, `K`, `M` | the Mercier index and the Glasser–Greene–Johnson layer coefficients at each surface | +| `SingularSurfaces/ca_left`, `ca_right` | the asymptotic (large, small) coefficients on either side of each surface, from ideal crossings | +| `SingularSurfaces/Delta_prime_matrix` | the Δ′ matrix in the PEST-3 convention, from `"Riccati"` or `"Galerkin"`; `Delta_prime_raw` is the side-major (left, right) matrix it is folded from, `Delta_coil` the edge coil-response columns | +| `SingularSurfaces/Kinetic/` | the surfaces of a kinetic solve | + +`Info/psilim` and `Info/qlim` record the outer edge of the integration domain. + +## Integration domain + +The integration runs from the axis (`qlow` sets the inner edge) to `psilim`, the flux of the target safety factor `qlim`. +`qlim` starts at `min(qmax, qhigh)`, with `qmax` the value of q at the equilibrium's `psihigh`. +With `truncate_at_rational_offset = true` (the default, for diverted equilibria where q diverges at the separatrix), `qlim` moves to `(m_last + rational_offset_fraction)/n`, a fixed fraction of a rational spacing past the last surface, stepped back by `1/n` until it no longer exceeds the bound it started from; a bound that falls on a rational, or less than `rational_offset_fraction/n` above one, therefore excludes that surface. +With `truncate_at_rational_offset = false` (limited or analytic equilibria with finite edge q), `qhigh` and `psihigh` set the domain directly. + +### Resistive-layer overlap cap + +Near the separatrix the rational surfaces crowd faster than their resistive layers narrow, and past the point where adjacent layers touch the matched-asymptotic treatment is not defined (Fitzpatrick (2025), listed under [InnerLayer Module](@ref)). +`psilim_from_layer_overlap = true` locates that point: each rational surface of each n in the run, including surfaces beyond the equilibrium grid placed on the logarithmic divergence of q at the X-points, gets a layer width from the kinetic profiles with the closures of `[SLAYER]` (a `SLAYERModel` in a script), and the domain is cut at the last surface before two adjacent layers overlap. +The cap needs kinetic profiles, acts before the `rational_offset_fraction` step, and can only narrow the domain; `examples/DIIID-like_truncation_example` exercises it. +The scan is written to `ForceFreeStates/LayerOverlap/`: | Dataset | Contents | |---|---| | `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius √(2ψ_t/B₀) | | `delta_s_abs`, `width_delta_s` | Riccati layer thickness, in metres and as a Δψ width | -| `width_visco`, `width_dr` | visco-resistive and diffusive-resistive (Eq. 100) Δψ widths; `width_dr` sets the cut | +| `width_visco`, `width_dr` | visco-resistive and diffusive-resistive Δψ widths; `width_dr` sets the cut | | `extrapolated` | 1 where the surface lies beyond the equilibrium grid | | `psilim_overlap`, `first_overlap_index` | the cut, and the first surface whose layer overlaps its neighbour | -| `applied` | 1 when the cap was passed to the truncation (`rational_offset_fraction` / `qhigh` may still cut deeper) | +| `applied` | 1 when the cap was passed to the truncation | -## Configuration reference +The perpendicular diffusivities that set the layer widths are described with the SLAYER inputs on the [Tearing](inner_layer.md) page. -The `[ForceFreeStates]` section holds the [`ForceBalanceProblem`](@ref) keywords and names the -integrator; the integrator's numerical settings sit in its own section, `[Forward]`, `[Riccati]` or -`[Galerkin]`, with the fields of that type. +## Numerical notes -```toml -[ForceFreeStates] -integrator = "Riccati" # "Forward" for dense xi profiles and kinetic runs; "Galerkin" for the RDCON outer-region solve - -# Mode space -nn = 1 # toroidal mode number, or consecutive numbers as a list ([1, 2, 3]) -delta_mlow = 0 # extra low poloidal modes (m < mlow) -delta_mhigh = 0 # extra high poloidal modes (m > mhigh) - -# Output -write_el_matrices = false # also write ForceFreeStates/EulerLagrangeMatrices (mpert² × npsi per matrix) +- The Riccati integrator splits the domain into chunks whose count is derived from the number of singular surfaces (`[Riccati] nchunks` overrides the target, `0` is automatic); every chunk integrates from identity initial conditions, so the result is independent of the thread count. +- The chunk ending at each rational surface is integrated backward, from the surface inward, because the forward propagator is ill-conditioned there; the Δ′ boundary-value problem is assembled from the split forward and backward propagators for the same reason. +- Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati basis is recovered by renormalizing whenever the solution exceeds `renorm_threshold`. +- The two crossing conventions differ slightly, so Riccati and Forward free-boundary energies agree to within a few percent, not to machine precision. -[Riccati] -nchunks = 0 # chunk-count target (0 = auto, from msing alone) -eulerlagrange_tolerance = 1e-8 # ODE relative tolerance -``` +Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md). -The number of Julia threads is controlled at startup via `-t N` or the `JULIA_NUM_THREADS` -environment variable; it is not a runtime parameter. +## API reference -## API Reference - -The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separately in -`docs/src/galerkin.md`. +The Galerkin solver's functions are documented on the [Galerkin Solver](galerkin.md) page. ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"] ``` - -## Example usage - -### Run stability analysis from a TOML configuration - -```julia -using GeneralizedPerturbedEquilibrium - -ffs = main(["examples/DIIID-like_riccati_deltaprime_example"]).ffs -println("Energy eigenvalue et[1] = ", real(ffs.free_boundary.et[1])) -``` - -The same solve from a script, with the deck's keys as keywords: - -```julia -eq = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada") -ffs = solve(eq, Riccati(; eulerlagrange_tolerance=1e-10); nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true, qlow=1.02) -``` - -### Inspect Δ' at singular surfaces - -```julia -for (s, sing) in enumerate(ffs.surfaces) - println("Surface $s: ψ = $(sing.psifac), m/n = $(sing.m[1])/$(sing.n[1])") - println(" Δ' = $(real(ffs.delta_prime.matrix[s, s]))") -end -``` - -## Notes - -- The standard path does not populate `delta_prime`; the canonical Δ' is the STRIDE BVP - `SingularSurfaces/Delta_prime_matrix` from the parallel FM path. `ca_l`/`ca_r` are filled - only by ideal surface crossings (kinetic and galerkin-matched runs emit zero-extent - `ca_left`/`ca_right` sentinels). -- The Riccati and parallel FM paths compute Δ' inline at each crossing, using the - direct diagonal formula (no GR permutation). The result in `delta_prime_col[ipert_res, i]` - equals `delta_prime[i]` to machine precision. -- `delta_prime_matrix` contains raw BVP coefficients, not asymptotic-normalized values; - its diagonal elements do **not** in general equal `delta_prime`. -- ODE step counts depend on the equilibrium profile and mode count; the `numsteps_init` - parameter sets the initial allocation but the solver adapts automatically. - -## See also - -- `docs/src/galerkin.md` — RDCON outer-region Galerkin Δ′ solver (part of this module) -- `docs/src/equilibrium.md` — build the `PlasmaEquilibrium` object required by this module -- `docs/src/vacuum.md` — vacuum response computed from the EL solution in `free_run` -- `docs/src/perturbed_equilibrium.md` — downstream singular coupling analysis using Δ' diff --git a/docs/src/workflow.md b/docs/src/workflow.md index d42ec4e062..bc69c32f7d 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -26,14 +26,14 @@ Every stage after `Equilibrium` reads the force-free states, and `KineticForces` | How does the plasma respond to an RMP or error field; do islands overlap? | `[ForcingTerms]` + `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `…/island_half_width`, `…/chirikov_parameter` | | What NTV torque does that field drive? | `[KineticForces]` (+ `[Equilibrium] kinetic_file`) | `KineticForces//total_torque` | | How tight must the coil tolerances be to avoid locking? | `[ErrorFields]` | `ErrorFields/Risk/locking_probability_percent` | -| What are the tearing growth rate and critical resonant field? | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `Tearing/Roots/gamma`, `Tearing/Roots/omega`, `Tearing/CriticalResonantField/` | -| What does a resistive inner layer do to the response? | `[Match]` + `[GGJ]` | `SingularSurfaces/Match/` | +| What are the tearing growth rate and critical resonant field? | `[Tearing]` (+ `[SLAYER]` or `[GGJ]`) | `Tearing/Roots/gamma`, `Tearing/Roots/omega`, `Tearing/CriticalResonantField/` | +| What does a resistive inner layer do to the response? | `[Match]` (+ `[GGJ]`) | `SingularSurfaces/Match/` | ### Headline outputs - **Free-boundary energy** `ForceFreeStates/FreeBoundaryStability/W_freeboundary` is the power-normalized free-boundary energy matrix, the sum of `W_plasma` and `W_vacuum` for the wall of the deck (no wall when `[Wall]` is absent). `eigenmode_energies` lists its generalized eigenvalues, the ideal δW per unit displacement of each free-boundary eigenmode, sorted most unstable first, so a negative first entry means the plasma is ideal-kink unstable and a positive one is the stability margin. -- **Tearing index** `SingularSurfaces/Delta_prime_matrix` is the Δ′ matching matrix in the PEST-3 convention, from which the outer-region tearing drive of each rational surface follows; a positive Δ′ means the outer region drives the tearing mode and the inner layer (`Tearing`) sets its growth rate. -- **Resonant field** `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field` is the area-weighted resonant normal field b̄ʳ `[T]` at each rational surface after the plasma response, the quantity that drives islands and sets locking thresholds. `island_half_width` is the vacuum-superposition island each resonant flux would open, w_s ∝ √|δψ_s|, and `chirikov_parameter` is the ratio of adjacent half-widths to the surface separation, so σ > 1 means the islands overlap. +- **Tearing index** `SingularSurfaces/Delta_prime_matrix` is the Δ′ matching matrix in the PEST-3 convention, whose diagonal is the outer-region Δ′ of each rational surface; a positive Δ′ means the outer region drives the tearing mode and the inner layer (`Tearing`) sets its growth rate. +- **Resonant field** `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field` is the area-weighted resonant normal field b̄ʳ `[T]` at each rational surface after the plasma response, the quantity that drives islands and sets locking thresholds. `island_half_width` is the vacuum-superposition island each resonant flux would open, w_s ∝ √|δψ_s|, and `chirikov_parameter` is that half-width over the half-distance to the neighbouring rational surface, so σ > 1 means the islands overlap. - **Torque** `KineticForces//total_torque` is the neoclassical toroidal viscosity (NTV) torque T_φ `[N·m]` that the perturbed equilibrium exerts on the plasma, with the sign convention of the [Conventions](conventions.md) page. ## What each stage computes @@ -66,7 +66,7 @@ Page: [Local stability](ballooning.md). ### Match -Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, section `[GGJ]`) and rebuilds the force-free states with the resulting layer response, so that the later stages see a resistive rather than perfectly shielded plasma; `ideal = true` keeps the shielded limit for comparison. +Matches the outer-region Δ′ to a resistive inner layer (`inner_model = "GGJ"`, the only model, with its options in an optional `[GGJ]` section) and rebuilds the force-free states with the resulting layer response, so that the later stages see a resistive rather than perfectly shielded plasma; `ideal = true` keeps the shielded limit for comparison. Outputs: `SingularSurfaces/Match/`. Page: [Tearing](inner_layer.md). @@ -97,7 +97,7 @@ Page: [Error fields](error_fields.md). Solves the tearing dispersion relation of each rational surface by matching the outer-region Δ′ to an inner-layer model: the two-fluid slab `SLAYER` (`[SLAYER]`) or the toroidal resistive-MHD `GGJ` (`[GGJ]`), selected by `[Tearing] inner_model`. `[Tearing.search]` sets the root search in the complex frequency plane, and `[Tearing.CriticalResonantField]` adds the resonant field at which the layer's shielding is overcome. -Needs a Δ′-producing integrator (`"Riccati"` or `"Galerkin"`) and, for `SLAYER`, a kinetic file. +Needs the Δ′ matrix of `"Riccati"` with `vac_flag = true` or `"Galerkin"` with `[Galerkin] rpec_flag = true`, and the kinetic profiles of `[Equilibrium] kinetic_file` (`SLAYER` reads its densities, temperatures, rotation and diffusivities there). Outputs: `Tearing/PerSurface/`, `Roots/`, `LayerWidths/`, `Diagnostics/`, `CriticalResonantField/`. Key references: [Glasser, Wang and Park (2016)](citations.md#InnerLayer-Module), [Park (2022)](citations.md#InnerLayer-Module). Page: [Tearing](inner_layer.md). @@ -154,7 +154,7 @@ Every key is documented on the module page of its section, and every example dec | `[Forward]`, `[Riccati]` | the radial ODE integrators | [Stability](stability.md) | | `[Galerkin]` | the finite-element outer-region solver | [Galerkin solver](galerkin.md) | | `[LocalStability]`, `[Ballooning]` | the local stability stage and its ballooning integration | [Local stability](ballooning.md) | -| `[Match]`, `[GGJ]` | the inner-layer matching stage and the GGJ layer model | [Tearing](inner_layer.md) | +| `[Match]`, `[GGJ]` | the inner-layer matching stage and the options of its GGJ layer model | [Tearing](inner_layer.md) | | `[ForcingTerms]`, `[[ForcingTerms.coil_set]]` | the external field: a file or coil sets | [Forcing terms](forcing_terms.md) | | `[PerturbedEquilibrium]` | the plasma response stage and its singular-coupling options | [Perturbed equilibrium](perturbed_equilibrium.md) | | `[KineticForces]` | the NTV torque stage and the kinetic solve of `ForceFreeStates` | [Kinetic forces](kinetic_forces.md) | From c04fd25aee3624882c29da1f080293b966204be4 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 2 Nov 2026 00:38:08 -0500 Subject: [PATCH 03/10] Docs - DOCS - Rewrite the remaining module pages around their physics, inputs and outputs Give kinetic_forces.md the NTV physics, the torque methods, the collision and distribution options, a deck snippet and an outputs table, and drop the Fortran line-number warning; move the kinetic file schema and ion composition to equilibrium.md, which gains a format table, a condensed radial-grid section and an outputs table. Resolve the coil normalization statement on forcing_terms.md (the Biot-Savart pipeline produces the unit-norm flux directly) and condense the page. Rebuild inner_layer.md as the Tearing page: outer-inner matching, SLAYER and GGJ, the dispersion search, the critical resonant field, a deck snippet, outputs, and the solver internals and benchmarks as numerical notes. Give error_fields.md a workflow list, drop the stale "later releases" sentence and condense the sampling, convergence, correction-design and naming sections. Give vacuum.md the vacuum energy, the wall shapes and the outputs. Correct the DIIID-like_ideal_example deck header, which named itself the Delta-prime fixture. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF --- docs/src/citations.md | 2 +- docs/src/equilibrium.md | 260 +++++--------- docs/src/error_fields.md | 112 ++---- docs/src/forcing_terms.md | 187 +++------- docs/src/inner_layer.md | 372 ++++++-------------- docs/src/kinetic_forces.md | 275 +++++---------- docs/src/vacuum.md | 121 ++++--- docs/src/workflow.md | 2 +- examples/DIIID-like_ideal_example/gpec.toml | 5 +- 9 files changed, 410 insertions(+), 926 deletions(-) diff --git a/docs/src/citations.md b/docs/src/citations.md index 7f8ac07b46..60e8c46c98 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -188,7 +188,7 @@ Combines the toroidal outer-region calculation with the two-fluid slab layer mod > *Nuclear Fusion* (2025). > DOI: [10.1088/1741-4326/ae4fdd](https://doi.org/10.1088/1741-4326/ae4fdd) (open access: arXiv:2511.07666) -Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain truncation](@ref). +Gives the diffusive-resistive layer width (Eq. 100) and, in Sect. 5.9, restricts the response calculation to the region inside the point where adjacent resonant layers first overlap. Both are implemented by the layer-overlap scan in `src/Tearing/LayerOverlap.jl`, which feeds the `psilim_from_layer_overlap` cap described in [Integration domain](@ref). --- diff --git a/docs/src/equilibrium.md b/docs/src/equilibrium.md index a7215f1de1..b1ed932a14 100644 --- a/docs/src/equilibrium.md +++ b/docs/src/equilibrium.md @@ -1,192 +1,102 @@ -# Equilibrium Module - -The Equilibrium module provides tools to read, construct, and analyze -magnetohydrodynamic (MHD) plasma equilibria. It wraps a variety of -equilibrium file formats (EFIT, CHEASE, SOL, LAR, and others), runs the -appropriate direct or inverse solver, and returns a processed -`PlasmaEquilibrium` object ready for downstream calculations. - -## Overview - -Key responsibilities of the module: - -- Read equilibrium input files and TOML configuration (see `EquilibriumConfig`). -- Provide convenient constructors for analytic / model equilibria (Large - Aspect Ratio, Solovev). -- Build spline representations used throughout the code (1D cubic and - 2D bicubic splines). -- Run the direct or inverse equilibrium solver and post-process results - (global parameters, q-profile finding, separatrix finding, GSE checks). - -The module exposes a small public API that covers setup, configuration, -and common analyses used by other GPEC components (e.g. `ForceFreeStates`, vacuum -interfaces). - -## EFIT solver strategies - -When `eq_type` is one of the EFIT-based options, three solver strategies are available: - -| `eq_type` | Method | Best for | -|-----------|--------|----------| -| `"efit"` | Geometric-angle field-line ODE | Default; fast and robust for standard cases | -| `"efit_arclength"` | Arc-length field-line ODE | Near-separatrix surfaces where the geometric-angle ODE becomes singular | -| `"efit_by_inversion"` | Contour.jl marching-squares → inverse solver | Highest geometric accuracy; avoids ODE singularities entirely | - -**`efit`** integrates field lines using the geometric angle η (polar angle from the -magnetic axis, 0→2π) as the independent variable. Position at each step is computed as -`R = R₀ + rfac·cos(η)`, so the RHS denominator is `Bz·cos(η) − Br·sin(η)` — the -projection of B_pol onto the radial direction. Near an X-point where Bp → 0 this -denominator passes through zero, causing a coordinate singularity: the ODE steps -become extremely small or the solver fails to converge. - -**`efit_arclength`** uses arc length along the flux surface as the independent variable -instead. The tangent direction `(dR/ds, dZ/ds) = ∇ψ⊥/|∇ψ|` has unit magnitude by -construction, so the position equations have no denominator and remain well-behaved -all the way to the separatrix. The 1/Bp factors in the accumulated integrals still -diverge near X-points, but their solver tolerances are set large so they do not -restrict the step size — only the position tracking drives adaptivity. The two methods -produce identical results when both succeed; `efit_arclength` extends the reliable -range of psihigh closer to 1. - -**`efit_by_inversion`** traces flux surface level sets directly from ψ(R,Z) using marching -squares, resamples each closed curve to a uniform geometric-angle grid, and feeds the -result into the inverse equilibrium solver (the same path used by CHEASE input). -The Cartesian evaluation grid is clipped to the separatrix bounding box and its -resolution is set adaptively from a bilinear interpolation error bound, so no manual -tuning is needed. - -## Radial grid packing - -With `grid_type = "auto"` and `mpsi = 0` (the defaults), the radial grid is -built by a **two-pass measured-curvature refinement** driven by the single accuracy knob -`psi_accuracy` (τ): - -1. **Pass 1** forms the equilibrium on a coarse three-region layout (geometric in - log(ψ) at the core, uniform in the middle, geometric in log(1−ψ) at the edge). -2. A knot density is derived from **one sizing rule** — the cubic-spline derivative - error model err(f′) ≈ h³|f''''|/24 ≤ τ·|f| — applied to three sources of - curvature |f''''|: - - **Measured** (mid-radius and edge): divided differences on the pass-1 nodes of the - 1D profiles (F, P, dV/dψ, q), the 2D geometry channels (r², η-offset, ν, Jacobian) - along sampled θ-lines (the bicubic ψ-axis shares the same knots), and the kinetic - profiles (n, T, ω_E) when loaded, so steep pedestal gradients attract knots. Nodal - values come from independent field-line integrations, so the estimate is immune to - inter-knot spline ringing. Curvature is measured against ρ = √ψ, in which the - equilibrium is regular at the magnetic axis — in ψ itself the geometry channels - behave as ψ^(k/2) (R−R₀ ~ √ψ), so their ψ-curvature diverges under refinement while - their ρ-curvature converges; the ρ-spacing maps back through dψ/dρ = 2√ψ, which by - itself packs √ψ-tight toward the axis. - - **Separatrix model** (edge floor, ψ ≥ 0.9): the same rule on q ≈ −A·ln(1−ψ) gives - geometric-in-log(1−ψ) packing with uniform relative q′ error, independent of A. - Normally inactive — the pass-1 layout is already log-packed at the edge, so the - measured density dominates — it only guards a pass-1 that under-sampled the - divergence. - - **Axis model** (core, ψ ≤ 0.03): the same rule on the power-law axis form gives - geometric-in-log(ψ) packing (constant-ratio spacing ∝ ψ down to `psilow`). Here the - model *replaces* measurement: nodal data from the smallest flux surfaces is - dominated by integration and axis-extrapolation noise, which grid refinement - amplifies rather than resolves. - - **Rational surfaces**: local packing around every q = m/n in the requested n range - (pinned as mandatory knots by the main driver), at any ψ including inside the - modeled core region. This local resolution of the ψ-splined Euler-Lagrange - coefficient matrices is what converges the Δ′ boundary-value problem. -3. The density integral is equidistributed into the knot vector, mandatory - rational-surface knots are inserted with a minimum-spacing snap guard, and - **pass 2** re-forms the equilibrium on the refined grid from the in-memory input - (no file re-read). - -Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines -the core, pedestal, and edge proportionally. The fixed `grid_type = "core_edge_packed"` -(sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed -layout) are still supported. A `PlasmaEquilibrium` built directly with `mpsi = 0` keeps -the coarse pass-1 grid until `solve` re-forms it; to refine it yourself, pass the nodes from -`refined_psi_grid` to `Equilibrium.reform`. - -The packing on the DIII-D-like example (n=1) compared to fixed `core_edge_packed` grids — note the -coarse spacing across the smooth mid-radius, the spacing dips at each rational surface, -and the core/pedestal/edge packing (`benchmarks/plot_grid_knot_placement.jl` regenerates -this figure): - -![Radial knot packing: auto two-pass vs core_edge_packed](assets/grid_knot_placement.png) - -Decomposing the density by source on the same example shows the pedestal band -(ψ_N ≈ 0.85–0.98) is driven by *measured* curvature, not the edge floor: the pressure, -q, and dV/dψ profiles contribute comparably, and the rzphi geometry channels — the -Grad-Shafranov response to the same pedestal p′ (Shafranov-shift / angle-offset -steepening) — contribute the most at every node in the band. Because the density is -measured from the formed solution, the packing follows the pedestal wherever it sits: - -![Knot-density decomposition by source](assets/density_decomposition.png) - -## API Reference +# Equilibrium -```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Equilibrium] -``` +`Equilibrium` reads or constructs the axisymmetric MHD equilibrium, maps it to the straight-field-line coordinates every later stage works in, and returns a `PlasmaEquilibrium` holding the flux-surface splines, the profiles and the global parameters of the plasma. +It is the first stage of every run and the only one that reads physics input files: the equilibrium itself and, optionally, the kinetic profiles. -## Important types - -- `EquilibriumConfig` — configuration container, built from the `[Equilibrium]` - section of `gpec.toml` or from keywords. -- `PlasmaEquilibrium` — the runtime structure containing spline fields, - geometry, profiles, and computed diagnostics (q-profile, separatrix, - etc.). -- `SolovevConfig`, `LargeAspectRatioConfig`, `TJAnalyticConfig` — the parameters of the - analytic equilibria. A deck states one as a `[Solovev]`, `[LargeAspectRatio]` or - `[TJAnalytic]` section, which selects the equilibrium, so it sets no `eq_type`. - `[TJAnalytic] solver = "inverse" | "direct"` picks the TJ pipeline. - -## Key functions - -- `PlasmaEquilibrium(file; eq_type="efit", kwargs...)`, `PlasmaEquilibrium(analytic_config; kwargs...)` - and `PlasmaEquilibrium(dd; kwargs...)` — read or solve the equilibrium and return a - `PlasmaEquilibrium`; `kwargs` are `EquilibriumConfig` fields. -- `PlasmaEquilibrium(config::EquilibriumConfig, source=nothing)` — the form every other - constructor ends in, with `source` an analytic config, an IMAS `dd`, or `nothing` to read `eq_filename`. -- `equilibrium_separatrix_find!(pe::PlasmaEquilibrium)` — locate - separatrix and related boundary geometry in-place. -- `equilibrium_global_parameters!(pe::PlasmaEquilibrium)` — populate - common global parameters (major radius, magnetic axis, volumes, etc.). -- `equilibrium_qfind!(pe::PlasmaEquilibrium)` — compute safety factor - (q) information across the grid. -- `equilibrium_gse!(pe::PlasmaEquilibrium)` — diagnostics on the - Grad–Shafranov solution. - -## Example usage - -Basic example: read a g-file and build an equilibrium +## What it computes -```julia -using GeneralizedPerturbedEquilibrium +From the poloidal flux ``\psi(R, Z)`` and the profiles ``F = R B_\phi`` and ``p`` of a Grad–Shafranov solution, the module traces the closed flux surfaces, finds the magnetic axis and the separatrix, integrates the safety factor ``q(\psi)`` and the surface quantities ``dV/d\psi``, and builds bicubic splines of the surface geometry ``R(\psi, \theta)``, ``Z(\psi, \theta)`` and of the toroidal-angle offset ``\nu(\psi, \theta)`` in the coordinates selected by `jac_type` (`"hamada"`, `"pest"`, `"boozer"`, `"equal_arc"`, …; see the [Conventions](conventions.md) page). +The normalized flux ``\psi_N \in [0, 1]`` runs from the axis to the boundary, and the rational surfaces q = m/n of the run are located on it. + +## Inputs + +`[Equilibrium] eq_type` selects the file format of `eq_filename`: + +| `eq_type` | Format | Notes | +|---|---|---| +| `"efit"` | EFIT g-file | flux surfaces traced by a field-line ODE in the geometric angle; the default | +| `"efit_arclength"` | EFIT g-file | the same traced in arc length, which stays regular near the separatrix and reaches higher `psihigh` | +| `"efit_by_inversion"` | EFIT g-file | level sets traced by marching squares and fed to the inverse solver; the most accurate geometry | +| `"chease"` (`"chease_binary"`) | CHEASE binary output | inverse solver | +| `"chease2"` (`"chease_ascii"`) | CHEASE ASCII output | inverse solver | +| `"ifile"` | inverse i-file, as TokaMaker writes | inverse solver | +| `"imas"` | an `IMASdd.dd` passed from a script | `PlasmaEquilibrium(dd; kwargs...)`; see `examples/DIIID-like_ideal_example_IMAS` | + +An analytic equilibrium is selected by its own section instead of `eq_type`: `[Solovev]` (the Solov'ev solution), `[LargeAspectRatio]` (a large-aspect-ratio model with prescribed profiles) or `[TJAnalytic]` (the profile family of Fitzpatrick's TJ code, `solver = "inverse" | "direct"`); their keys are the fields of `SolovevConfig`, `LargeAspectRatioConfig` and `TJAnalyticConfig`. + +### Kinetic profiles + +`[Equilibrium] kinetic_file` attaches the density, temperature, rotation and transport profiles that the kinetic solve, [Kinetic Forces](kinetic_forces.md) and the `SLAYER` layer of [Tearing](inner_layer.md) read, so every kinetic consumer sees one plasma. +Two formats are read, by extension: + +- **HDF5** (`.h5`, `.hdf5`), the GPEC kinetic schema: datasets at the file root, each with a `units` attribute, on the normalized flux `psi`: required `n_e`, `T_i`, `T_e`, `omega_E` (and `n_i`, defaulting to `n_e`); optional `omega_tor`, `chi_e`, `chi_i`, `chi_phi` (perpendicular heat and momentum diffusivities) and `zeff`; further `n_*` datasets are named per-species densities. Densities in m⁻³, temperatures in eV, frequencies in rad/s, diffusivities in m²/s. `write_kinetic_h5` writes the format. +- **ASCII** (`.gpeckf`, `.kin`, `.dat`): the six-column table `psi_n n_i n_e T_i T_e omega_E`, header rows skipped. -pe = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada") +The ion composition is part of the same data: `[[Equilibrium.ion_species]]` lists the main ions, each with `z`, `m` and either `fraction` (its share of the file's total `n_i`) or `density` (a named `n_*` dataset), and `impurity_charge` and `impurity_mass` describe the impurity that closes quasineutrality; the default is one deuterium ion carrying all of `n_i` with a carbon impurity. +``Z_\mathrm{eff}(\psi)`` is taken from the `zeff` dataset when present, else from `n_i/n_e` with the impurity charge. +In a script the same data is attached with `attach_kinetic_profiles!(equil, "kinetic.h5"; ion_species=[IonSpecies(; z=1, m=2, fraction=0.5), IonSpecies(; z=1, m=3, fraction=0.5)], impurity_charge=6, impurity_mass=12)`. -println("Magnetic axis: ", pe.ro, ", ", pe.zo) -println("q(0) = ", pe.params.q0) +### Radial grid -# Find separatrix (in-place) and inspect results -GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_separatrix_find!(pe) -println("rsep = ", pe.params.rsep) +With `grid_type = "auto"` and `mpsi = 0` (the defaults) the radial knots are placed by a two-pass refinement driven by one accuracy knob, `psi_accuracy`: the equilibrium is formed on a coarse grid, the curvature of the profiles, of the surface geometry and of the kinetic profiles is measured on it, and the knot density that keeps the cubic-spline derivative error below the tolerance is equidistributed, with geometric packing toward the axis and the separatrix and a mandatory knot at every rational surface of the run. +Every region's knot count scales as `psi_accuracy^(-1/3)`, so tightening it refines the core, the pedestal and the edge together; the rational-surface knots are what converge the Δ′ boundary-value problem. +The fixed layouts `grid_type = "core_edge_packed"`, `"pow1"` and `"uniform"` with an explicit `mpsi` remain available. + +![Radial knot packing on the DIII-D-like example: the automatic two-pass grid against a fixed core_edge_packed grid.](assets/grid_knot_placement.png) + +## Running it + +From `examples/DIIID-like_ideal_example/gpec.toml`: + +```toml +[Equilibrium] +eq_type = "efit" +eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" +kinetic_file = "TkMkr_D3Dlike_Hmode_kinetic.h5" +jac_type = "hamada" +grid_type = "auto" +psilow = 1e-4 # innermost flux surface +psihigh = 0.995 # outermost flux surface +mpsi = 0 # 0: automatic radial grid +psi_accuracy = 0.001 # spline accuracy that sizes the automatic grid +mtheta = 256 # poloidal resolution +impurity_charge = 6 +impurity_mass = 12 + +[[Equilibrium.ion_species]] +z = 1 +m = 2 +fraction = 1.0 ``` -Analytic / testing example: construct a large-aspect-ratio model +The same from a script, with the deck keys as keywords: ```julia using GeneralizedPerturbedEquilibrium -larcfg = GeneralizedPerturbedEquilibrium.Equilibrium.LargeAspectRatioConfig(r0=10.0, a=1.0, beta0=1e-3) -pe = PlasmaEquilibrium(larcfg) # the config type implies the equilibrium kind - -println("Built LAR equilibrium with a = ", larcfg.a) +eq = PlasmaEquilibrium("examples/DIIID-like_ideal_example/TkMkr_D3Dlike_Hmode.geqdsk"; jac_type="hamada", psihigh=0.995) +eq = PlasmaEquilibrium(LargeAspectRatioConfig(; r0=10.0, a=1.0, beta0=1e-3)) # an analytic equilibrium ``` -## Notes and Caveats +Relative paths in the deck are resolved against the deck directory. +A `PlasmaEquilibrium` built with `mpsi = 0` keeps the coarse first-pass grid until a `solve` re-forms it; `Equilibrium.reform` with the nodes of `refined_psi_grid` does the same by hand. -- `EquilibriumConfig` is constructed from the `[Equilibrium]` section of `gpec.toml`. - Paths that are not absolute are resolved relative to the TOML file location. -- The Equilibrium module contains readers for EFIT and CHEASE formats. Ensure the - required data files are present and paths are set correctly in `gpec.toml`. +## Outputs -## See also +The group is written once `[ForceFreeStates]` has run, on the grid of the solve. -- `docs/src/stability.md` — ideal MHD stability analysis built on top of the equilibrium +| Dataset | Meaning | Units | +|---|---|---| +| `Equilibrium/R_axis`, `Z_axis`, `B_T_axis`, `R_mean`, `a_mean`, `aspect_ratio`, `kappa`, `delta_upper`, `delta_lower` | axis position and field, mean radii and shape parameters | m, T | +| `Equilibrium/I_p`, `volume` | plasma current and volume | A, m³ | +| `Equilibrium/beta_N`, `beta_t`, `beta_p_1`, `l_i_1`, … | normalized, toroidal and poloidal beta and internal inductance | | +| `Equilibrium/q_axis`, `q_95`, `q_edge`, `q_min`, `q_max`, `q_extrema_psi` | the safety factor at the axis, the 95% surface and the edge, and its extrema | | +| `Equilibrium/psi_boundary`, `psi_total`, `I_p_sign`, `B_T_sign` | the flux normalization and the sign conventions of the file | | +| `Equilibrium/Profiles/psi`, `q`, `2piF`, `mu0p`, `dVdpsi` | the 1-D profiles on ψ_N: q, 2πF, μ₀p and dV/dψ_N | T·m (`2piF`), T² (`mu0p`), m³ (`dVdpsi`) | +| `Equilibrium/Geometry/psi`, `theta`, `rcoords`, `offset`, `nu`, `jac` | the 2-D surface geometry on (ψ_N, θ): the radial coordinate, the angle offset, the toroidal offset ν and the Jacobian | | + +## API reference + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Equilibrium] +``` diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index 993b4a7471..86adfc0c93 100644 --- a/docs/src/error_fields.md +++ b/docs/src/error_fields.md @@ -1,4 +1,4 @@ -# ErrorFields Module +# Error Fields The `ErrorFields` module quantifies how sensitive a perturbed equilibrium's resonant drive is to the placement of each coil set. It is the foundation of an error-field tolerance @@ -7,6 +7,14 @@ misaligned coil produce?"* is linear in the coil's motion, so the module lineari and everything downstream (tolerance Monte Carlo, locking risk, allowable-tolerance scans) becomes linear algebra on a small table. +## The workflow + +1. **Coil sensitivities** (`[ErrorFields]`, `ErrorFields/CoilSensitivities/`): the resonant spectrum of each coil set as built and its derivatives with respect to rigid shifts and tilts. +2. **Tolerances** (`tolerance_file`): the manufacturing tolerances of each coil set, in their own TOML file. +3. **Monte Carlo** (`[ErrorFields.MonteCarlo]`, `ErrorFields/MonteCarlo/`): the distribution of the dominant-mode overlap ``|\delta|`` over sampled misalignments, intrinsic and with error-field correction. +4. **Locking risk** (`[ErrorFields.scenario]`, `[ErrorFields.Risk]`, `ErrorFields/Risk/`): the overlap distribution against the empirical penetration-threshold scaling, giving the locking probability and its scan against tolerance scale, from which the allowable tolerance follows. +5. **Correction limits** (`[ErrorFields.NTV]`, `ErrorFields/NTV/`): how much correction current each array may carry before its non-resonant field's NTV torque spends the rotation budget, and how much it needs. + ## What is computed For every named coil set of a coil-forced run, the module evaluates the root-area-weighted @@ -17,6 +25,7 @@ about the machine axes in degrees. The spectra are placed on the run's ``(m, n)` conformed with the control-surface operator, so they are exactly the vectors the resonant coupling matrix and its singular vectors act on (see [Dominant resonant-coupling mode](perturbed_equilibrium.md#Dominant-resonant-coupling-mode)). +In the snippets below `EF` is `GeneralizedPerturbedEquilibrium.ErrorFields` and `AEF` is `GeneralizedPerturbedEquilibrium.Analysis.ErrorFields`. The overlap of a coil set with singular mode ``k`` of the coupling matrix, normalized by the axis toroidal field, is the dimensionless @@ -77,7 +86,7 @@ coils by error field per millimetre of shift and per tenth of a degree of tilt, rational surfaces and over the edge only. The tilt pivot matters for multi-filament winding packs: `"conductor"` rotates each filament -about its own arc-length centre, the Fortran `coil_read` convention inherited by the OMFIT +about its own arc-length centre, the convention of the Fortran GPEC coil loader inherited by the OMFIT tolerance tool, while `"set"` rotates the pack rigidly about its common centre, which is what an engineering axis-line tolerance constrains. Single-conductor sets give identical results either way. @@ -88,8 +97,8 @@ Manufacturing tolerances are engineering data, reviewed and versioned independen run, so they live in their own TOML file named by `tolerance_file` in `[ErrorFields]` (a path relative to the deck directory). The run validates the file, echoes its text into `Input/RawInputs/ErrorFields/tolerance_toml_raw` for replay, and leaves the device numbers -where they belong: outside the repository. The sampling and Monte Carlo stages that consume -these tolerances follow in later releases; this release fixes the format. +where they belong: outside the repository. The sampling and Monte Carlo stages below consume +these tolerances. ```toml # Fallbacks for keys a coil block leaves out (every key optional) @@ -160,23 +169,9 @@ tight = EF.update(ts; other_field=EF.update(ts.other_field; magnitude=1e-5)) ## Sampling a tolerance -A tolerance is one number, the radius of the disk the coil centre may lie in; the direction is -random and the radial density is a `RadialDistribution`: `Flat` (uniform in radius, the OMFIT -tool's `flat`), `UniformArea` (uniform over the disk), `Hollow` (peaked toward the edge, the -OMFIT default), `Ring` (always on the edge), or `PowerLaw(p)`. `sample_disk` draws a point as -`Δx + iΔy`, `sample_uncertainty` adds the Gaussian uncertainty on where the coil actually sits, -and the two tolerance models combine them: `sample_additive` draws a shift (metres) and a tilt -(degrees) independently, while `sample_cylinder` confines the coil's axis line to a cylinder -of radius `R` and half-height `z_top`, drawing its two endpoints and deriving the correlated -midplane shift and lean. Tilts are the rotation angles `θx + iθy` about the machine axes in -the sense `apply_transforms` uses. Every sampler takes the random generator and optional fixed -directions, so coherent groups sharing a direction pass one phase and draw their own radii. - -```julia -rng = Random.Xoshiro(1) -Δ = EF.sample_disk(rng, 0.5e-3, EF.Hollow()) # a point in a 0.5 mm disk, m -Δ, θ = EF.sample_cylinder(rng, 0.5e-3, 1.5, EF.randpow(EF.Flat())) # correlated shift [m] and tilt [deg] -``` +A tolerance is the radius of the disk the coil centre may lie in; the direction is random and the radial density is a `RadialDistribution` (`Flat`, `UniformArea`, `Hollow`, `Ring` or `PowerLaw(p)`), with a Gaussian placement uncertainty added on top. +The `"additive"` model draws a shift and a tilt independently; the `"cylinder"` model confines the coil's axis line to a cylinder of the tolerance radius and the given half-height and derives the correlated midplane shift and lean from its two endpoints. +Coherent groups share one random direction and draw their own radii. ## Tolerance Monte Carlo @@ -295,18 +290,7 @@ A threshold scaling is fitted for one toroidal mode number, so the in-run risk a file-path entry point refuse a run that spans several (`single_toroidal_mode`); project onto one n's coupling first. -A number quoted at the target risk level should be shown free of sampling and binning bias. -`risk_convergence` repeats the Monte Carlo over a list of sample counts and a list of bin -counts, holding the other at its control value, and returns the locking probability with its -batch spread for each: the spread must fall as `1/√nsample`, and the value must not move with -`nbins` by more than the spread. The spread alone cannot see a binning bias, which is why the bin -sweep is there. - -```julia -conv = EF.risk_convergence(table, ts, coil_sets, sc, scen; nsamples=[10^5, 10^6, 10^7], nbins_list=[100, 300, 3000]) -conv.locking_probability_percent_by_nsample, conv.locking_probability_spread_percent_by_nsample -conv.locking_probability_percent_by_nbins -``` +A number quoted at the target risk level should be free of sampling and binning bias: `risk_convergence` repeats the Monte Carlo over lists of sample and bin counts, and the batch spread must fall as ``1/\sqrt{n_\mathrm{sample}}`` while the value must not move with the bin count by more than the spread. ## Plots and coil-array phasing @@ -400,25 +384,11 @@ ARC paper 5–20 N·m for ARC; pick a value for the machine being analysed. ### Designing the correction -The NTV limits say how much correction current an array *may* carry; `correction_requirement` -says how much it *needs*, and on which surfaces the correction falls short. For one array the -dominant-mode answer is the ratio `−δ_source/δ_array`; for several arrays the dominant-mode -condition alone does not fix the currents (one complex equation in several unknowns leaves a family -of solutions in which any two arrays can cancel each other, of which the result names the -member with the least total ampere-turns squared); the question worth a function is the joint one over every rational surface: the complex factors on all arrays that minimize the resonant -field `C·(b̃_source + Σ_k f_k b̃_k)` in the least-squares sense, the field left on each surface -after that and after each array's dominant-mode correction alone, and how much of the source's -spectrum each array can see. Each current comes twice: as a complex factor on the array's -current as given, and in kilo-ampere-turns through the array's `ampere_turns_kat` -(`|winding multiplier| × max|I| / 1000`, the same normalization as the EFC couplings), which is the -form to compare arrays whose patterns were given at different currents and to set beside the -NTV-limited current. Any `CoilOverlap` can be a source or an array, so an as-built assembly -(`combine_overlaps`) can be corrected with arrays that were never part of the run. Over the tolerance samples the needed -current is a distribution, and `uncorrectable_probability` is the fraction of the Monte Carlo's -overlap beyond what the NTV-limited array can correct at all: the explicit-correction counterpart -of the corrected locking probability's `efc_factor` model. The Monte Carlo leaves the -correction arrays out of its overlap, so its as-designed current and a `correction_requirement` -whose source combines the error-field coils describe the same field. +The NTV limits say how much correction current an array *may* carry; `correction_requirement` says how much it *needs*, and on which surfaces the correction falls short. +For one array the dominant-mode answer is the ratio ``-\delta_\mathrm{source}/\delta_\mathrm{array}``; for several arrays the dominant-mode condition alone leaves a family of solutions (the result names the member with the least total ampere-turns squared), so the joint question is solved over every rational surface: the complex factors on all arrays that minimize the resonant field ``C\,(\tilde b_\mathrm{source} + \sum_k f_k \tilde b_k)`` in the least-squares sense, the field left on each surface, and how much of the source's spectrum each array can see. +Each current is given twice, as a complex factor on the array's current as given and in kilo-ampere-turns (`ampere_turns_kat`, the normalization of the EFC couplings), which is the form to compare arrays given at different currents and to set beside the NTV-limited current. +Any `CoilOverlap` can be a source or an array, so an as-built assembly (`combine_overlaps`) can be corrected with arrays that were never part of the run. +Over the tolerance samples the needed current is a distribution, and `uncorrectable_probability` is the fraction of the Monte Carlo overlap beyond what the NTV-limited array can correct, the explicit-correction counterpart of the `efc_factor` model. ```julia req = EF.correction_requirement("gpec.h5", "F_coils_as_built", ["d3d_c", "iu", "il"]) @@ -471,11 +441,7 @@ AEF.plot_linearity_check(lin) # stems per coil, log axis AEF.plot_linearity_residuals("gpec.h5"; at=:step, yscale=:log10) # the stored step residuals, readable below 1 % ``` -A coil set whose field at the run's toroidal mode is round-off (an axisymmetric hoop at n = 1) has -no resonant fraction: `coil_overlaps`, `combine_overlaps` and `plot_coil_sensitivities(quantity=:fraction)` -report `NaN` for it rather than the ratio of two noise vectors, which would read as a fraction -between 1 % and 30 %; the floor is `1e-8` of the largest field norm among the coil sets judged -together (`resonant_fraction_percent`). +A coil set whose field at the run's toroidal mode is round-off (an axisymmetric hoop at n = 1) has no resonant fraction, and `coil_overlaps`, `combine_overlaps` and `plot_coil_sensitivities(quantity=:fraction)` report `NaN` for it rather than the ratio of two noise vectors. ## Comparing coil revisions @@ -518,36 +484,8 @@ when only the overlaps are wanted. ## Result names -Every result struct in `ErrorFields` (`CoilSensitivities`, `SensitivityTable`, `CoilOverlap`, -`MonteCarloResult`, `RiskResult`, `ToleranceScan`, `PhasingMap`, `EFCCoupling`) names its fields by one -grammar, and every field name is also its HDF5 dataset name under `ErrorFields/`: - -``` -[_][_efc][_][_] -``` - -- **quantity** is the physical thing, spelled out: `delta` is the dimensionless dominant-mode overlap - `v·b̃ / B_T0`, a complex phasor per coil set, and `abs_delta` its magnitude; `threshold` is the - penetration threshold in the same units; `locking_probability` is the probability of locking (never - abbreviated, and never called a risk in a name); `resonant_fraction` is `100·|v·b̃| / ‖b̃‖`; - `field`, `spectrum`, `torque`, `current`, `shift_sensitivity` and `tilt_sensitivity` are what they say. -- **instance** says which geometry or value the quantity is evaluated for: `_as_designed`, the - unperturbed design; `_sampled`, over the Monte Carlo misalignment samples; `_worst_case`, every term - at its tolerance edge and in phase; `_fit`, a scaling-law value at the fitted exponents without the - fit's scatter; `_total`, the coherent sum over coil sets (a per-coil vector carries no qualifier and - documents its shape as `[ncoil_set]`). A result evaluated on whatever geometry the caller passed - (`CoilOverlap`, `EFCCoupling`) carries no instance qualifier. The word *nominal* is not used in names. -- **`_efc`** marks a value with error-field correction applied, as in the TOML keys `efc_factor` and - `efc_coils`. -- **statistic**: `_mean`, `_pdf`, `_batches`, `_spread` (range over the Monte Carlo batches), - `_samples`, `_bin_edges`. -- **unit** is present whenever the quantity has one: `_percent`, `_t` (tesla), `_m`, `_mm`, `_deg`, - `_per_m`, `_per_deg`, `_per_mm`, `_per_kat`, `_per_kat2`. `delta` and probabilities given as fractions - carry none; a probability stored as a percentage always carries `_percent`. - -The order is fixed, so `locking_probability_efc_batches_percent` reads as "the locking probability, -corrected, per batch, in percent", and `abs_delta_total_as_designed` as "the magnitude of the total -overlap of the design as built". +Every result field of `ErrorFields` is also its dataset name under `ErrorFields/`, built as `[_][_efc][_][_]`: `delta` is the dimensionless dominant-mode overlap and `abs_delta` its magnitude, `_as_designed`, `_sampled` and `_worst_case` say which geometry it is evaluated for, `_efc` marks error-field correction applied, and a probability stored as a percentage always carries `_percent`. +The grammar is specified in `docs/development/hdf5-conventions.md` of the repository. ## API Reference diff --git a/docs/src/forcing_terms.md b/docs/src/forcing_terms.md index 256fbfebc4..a33c4a0d3c 100644 --- a/docs/src/forcing_terms.md +++ b/docs/src/forcing_terms.md @@ -1,141 +1,68 @@ -# ForcingTerms Module +# Forcing Terms -The ForcingTerms module specifies external magnetic field perturbations for the perturbed -equilibrium calculation. It supports reading from ASCII or HDF5 files, as well as computing -perturbations directly from coil Biot-Savart calculations. +`ForcingTerms` specifies the external magnetic perturbation that drives the perturbed equilibrium: the spectrum of the normal field on the control surface at `psilim`, given either as a file of (m, n) amplitudes or computed by Biot–Savart from coil geometry. +The `[ForcingTerms]` section is read when `[PerturbedEquilibrium]` runs; its result is the `RMPField` that `perturbed_equilibrium` applies to the plasma response. -## Normalization conventions +## Normalization -Two normalization tags are supported, specified inside the file (not in `gpec.toml`): +Internally every forcing is the unit-norm flux spectrum ``\Phi_x``, the Fourier amplitude of ``2\pi R\,(B_R\,\partial Z/\partial\theta - B_Z\,\partial R/\partial\theta)`` over the unit-norm angles ``\theta, \zeta \in [0, 1]``, the quantity the Fortran GPEC calls `Phi_x`. +A forcing file declares its own normalization by a tag; a coil set needs none, because the Biot–Savart pipeline produces ``\Phi_x`` directly. -| Tag | Quantity | Units | Notes | -|-----|----------|-------|-------| -| `"normal_field_T"` | Fourier modes of B·n̂ | Tesla | **Default** — most intuitive for users | -| `"sfl_flux_Wb"` | Fourier modes of R×(B_R ∂Z/∂θ − B_Z ∂R/∂θ) | T·m² | 2π-angle SFL flux (user input) | +| Tag | Quantity | Units | On load | +|---|---|---|---| +| `"normal_field_T"` (default) | Fourier modes of ``B\cdot\hat n`` over ``\theta \in [0, 2\pi]`` | T | weighted by the area element ``J\lvert\nabla\psi\rvert`` and projected onto every harmonic of the run; a poloidally varying weight mixes modes, so the file must cover its sidebands (a warning asks for a wider m range when an edge harmonic carries more than 10% of the flux) | +| `"sfl_flux_Wb"` | the straight-field-line flux ``R\,(B_R\,\partial Z/\partial\theta - B_Z\,\partial R/\partial\theta)`` over ``\theta \in [0, 2\pi]`` | T·m² | multiplied by ``(2\pi)^2`` | +| `"unit_flux"` | ``\Phi_x`` itself, as a run records its forcing for replay | T·m² | loaded as is | -### Angle convention +## File formats -Julia's ForcingTerms use the **unit-norm angle convention** (θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1]) -internally, matching Fortran GPEC and the Equilibrium/ForceFreeStates modules. - -The coil Biot-Savart pipeline produces mode amplitudes directly equal to Fortran's `Phi_x`: - -``` -Julia output ≡ Fortran Phi_x -``` - -The formula is `Phi_x(θ_norm, ζ_norm) = 2π × R × (B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm)`, where -the `2π` comes from the toroidal Jacobian ∂r/∂ζ_norm = 2π·R·ê_φ. - -Users who provide forcing files may use the more intuitive **2π-angle Fourier amplitudes** -(e.g. the amplitude of cos(mθ − nζ) where θ ∈ [0, 2π]). These are tagged `"sfl_flux_Wb"` -and multiplied by (2π)² internally to reach unit-norm convention. - -### Automatic conversion - -Files tagged `normal_field_T` are automatically converted to unit-norm Phi_x convention -during loading by weighting with the control-surface area element J|∇ψ|. This weight varies -poloidally, so the conversion is a mode-mixing operation and the flux is projected onto -every harmonic of the run (`mlow:mhigh`), sidebands included. - -The file's spectrum must also extend far enough that the harmonics it omits are negligible, -since those would mix into the run's range too. GPEC estimates this from the file's edge -harmonics: if either one alone carries more than 10% of the converted flux, a warning asks for -a wider m range. Files with fewer than three harmonics are not checked. - -Files tagged `sfl_flux_Wb` (2π-angle SFL flux) are scaled by (2π)² to reach unit-norm. - -The coil Biot-Savart code (`forcing_data_format = "coil"`) always produces unit-norm -amplitudes directly — no conversion is applied. - -## ASCII file format +ASCII (`forcing_data_format = "ascii"`): an optional `# normalization: ` line, comment lines starting with `#`, and one `n m real(amplitude) imag(amplitude)` row per mode (the imaginary part defaults to 0): ``` # normalization: normal_field_T -# Format: n m real(amplitude) imag(amplitude) 1 6 1.0e-4 0.0 1 7 5.0e-5 2.0e-5 ``` -- The `# normalization: ` line is optional; default is `normal_field_T`. -- Lines starting with `#` are comments and are ignored by the data parser. -- Columns: toroidal mode n, poloidal mode m, real part of amplitude, imaginary part (optional, default 0). -- Amplitudes are in Tesla (for `normal_field_T`) or T·m² (for `sfl_flux_Wb`). - -## HDF5 file format - -Required datasets: -- `n`: integer array of toroidal mode numbers -- `m`: integer array of poloidal mode numbers -- `amplitude_real`: float array of real parts -- `amplitude_imag`: float array of imaginary parts (optional) - -Optional root attribute: -- `normalization`: string, either `"normal_field_T"` (default) or `"sfl_flux_Wb"` +HDF5 (`forcing_data_format = "hdf5"`): datasets `n`, `m`, `amplitude_real` and optional `amplitude_imag`, with the optional root attribute `normalization`: -Example (Julia): ```julia using HDF5 h5open("forcing.h5", "w") do f f["n"] = [1, 1, 1] f["m"] = [5, 6, 7] f["amplitude_real"] = [1e-4, 2e-4, 1e-4] - f["amplitude_imag"] = [0.0, 0.0, 0.0] HDF5.attrs(f)["normalization"] = "normal_field_T" end ``` -## Coil Biot-Savart format - -Set `forcing_data_format = "coil"` in `[ForcingTerms]` TOML to compute the field -directly from coil geometry. Output is always in `sfl_flux_Wb` convention. - -The coil pipeline: -1. `sample_boundary_grid`: evaluate plasma boundary geometry at (mtheta × nzeta) grid -2. `compute_biot_savart_boundary!`: compute B at all grid points via Biot-Savart law -3. `project_normal_flux!`: compute Phi_x = 2π×R×(B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm) -4. `fourier_decompose_bn`: Fourier decompose to get bmn in unit-norm (= Phi_x) convention - -Steps 2–4 are linear in the coils, so a `CoilForcingGrid` (step 1 plus the observation-point -layout, built once per equilibrium and `n`) can be reused to evaluate coil sets one at a time -with `coil_forcing_modes`; the per-set spectra sum to the assembly's. Use this when many coil -variants must be compared against one plasma — moved or re-wound coils, or one coil at a time — -without re-sampling the boundary: - -```julia -forcing_grid = CoilForcingGrid(equil, cfg, n; psi=ffs.psilim) -spectra = [coil_forcing_modes(cs, forcing_grid, n, m_low, m_high) for cs in coil_sets] -``` - -The toroidal angle for each observation point at SFL grid coordinate (θ_i, ζ_j) is: -``` -φ_phys(i, j) = -helicity × (2π × ζ_j + ν(ψ, θ_i)) -``` -where `helicity = sign(Bt) × sign(Ip)` and `ν(ψ, θ)` is the SFL toroidal coordinate correction -stored in `equil.rzphi_nu`. This matches Fortran's `phi = -helicity*(twopi*czeta + crzphi_f(3))`. -For DIII-D (Bt < 0, Ip > 0 → helicity = -1), ν ranges from about -3.4 to +3.4 radians at the -boundary — omitting it completely scrambles the Fourier mode decomposition. - -### Coil geometry sources +## Coils -Each `[[ForcingTerms.coil_set]]` block supplies one coil set from one of several sources, -selected by `source`. Shifts and tilts (`shiftx/y/z`, `tiltx/y/z`, `n_tilt`, ...) apply to -every source, so analytic coils can be perturbed exactly like file-based ones. +`forcing_data_format = "coil"` computes the field from coil geometry: the plasma boundary is sampled on an `mtheta_coil × nzeta_coil` grid, the Biot–Savart field of every coil set is evaluated there, projected onto the normal flux and Fourier decomposed. +The toroidal angle of each observation point includes the straight-field-line offset ``\nu(\psi, \theta)`` and the helicity ``\mathrm{sign}(B_T)\,\mathrm{sign}(I_p)`` of the equilibrium, without which the decomposition would be wrong. +Each `[[ForcingTerms.coil_set]]` block supplies one coil set, selected by `source`; shifts and tilts (`shiftx/y/z`, `tiltx/y/z`, `n_tilt`, …) apply to every source, which is what the [Error Fields](error_fields.md) stage sweeps. +From `examples/DIIID-like_ideal_example/gpec.toml` (the first block) and the other sources: ```toml -# Legacy ASCII .dat file (unchanged; resolves to {dat_dir}/{machine}_{name}.dat) +[ForcingTerms] +forcing_data_format = "coil" +machine = "d3d" +mtheta_coil = 480 +nzeta_coil = 32 + +# ASCII .dat file, resolved to {dat_dir}/{machine}_{name}.dat [[ForcingTerms.coil_set]] -name = "il" -currents = [1e4, 1e4, 1e4, 1e4, 1e4, 1e4] +name = "c" +currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0] -# Modern HDF5 file; `name` selects one set from a multi-set file (omit if it holds one) +# HDF5 file; `name` selects one set of a multi-set file [[ForcingTerms.coil_set]] source = "file" h5_file = "/path/to/my_coils.h5" name = "il" currents = [1e4, 1e4, 1e4, 1e4, 1e4, 1e4] -# Analytic PF hoop — a horizontal circle; tilt breaks axisymmetry to give n ≠ 0 forcing +# Analytic PF hoop; a tilt breaks axisymmetry to give n ≠ 0 forcing [[ForcingTerms.coil_set]] source = "pf_hoop" radius = 1.8 @@ -143,7 +70,7 @@ height = 0.9 currents = [5e3] tiltx = [2.0] -# Analytic window-pane (picture-frame) array — corner mode: two opposite [R, Z] corners +# Analytic window-pane array, by two opposite [R, Z] corners ... [[ForcingTerms.coil_set]] source = "window_pane" n_coils = 8 @@ -151,10 +78,7 @@ rz_corners = [[2.2, -0.6], [2.2, 0.6]] gap_fraction = 0.15 currents = [1e3, -1e3, 1e3, -1e3, 1e3, -1e3, 1e3, -1e3] -# Analytic window-pane — standoff mode: placed by physical distance off the plasma surface. -# The frame center sits `standoff` m along the outward surface normal at cylindrical -# `poloidal_angle` (degrees); legs span `poloidal_length` m, tilted by `poloidal_tilt` -# (degrees, 0 = tangential). Requires an equilibrium (the main forcing pipeline supplies it). +# ... or by standoff off the plasma surface at a poloidal angle (degrees) [[ForcingTerms.coil_set]] source = "window_pane" n_coils = 2 @@ -165,8 +89,7 @@ poloidal_tilt = 0.0 gap_fraction = 0.15 currents = [1e4, -1e4] -# Analytic helical array on a circular cross-section torus (pitch = m/n). -# theta_lo/theta_hi (poloidal extent, degrees) default to a full 0–360 turn. +# Analytic helical array on a circular torus, pitch m_hel/n_hel [[ForcingTerms.coil_set]] source = "helical" R0 = 1.7 @@ -177,47 +100,21 @@ n_coils = 4 currents = [1e3, 1e3, 1e3, 1e3] ``` -### HDF5 coil geometry format +A coil HDF5 file stores one subgroup per coil set under a parent group (default `coils`): `//x`, `y`, `z` as `ncoil × s × nsec` Cartesian metres, an optional `currents` dataset, and the winding multiplier as the attribute `nw`; `convert_coil_dat_to_h5` and `convert_coil_h5_to_dat` translate the ASCII format, `read_coil_h5` and `write_coil_h5` read and write it. -A coil HDF5 file stores one **subgroup per coil set** under a parent group (default `coils`), -so a single file can hold many sets and the reader ignores any unrelated content elsewhere: +The field is linear in the coils, so a `CoilForcingGrid` built once per equilibrium and n evaluates coil sets one at a time, for comparing coil variants against one plasma without re-sampling the boundary: -```text -//x Float64 ncoil×s×nsec Cartesian metres -//y Float64 ncoil×s×nsec -//z Float64 ncoil×s×nsec -//currents Float64 ncoil optional -/ attr nw Float64 winding multiplier, default 1.0 +```julia +forcing_grid = CoilForcingGrid(equil, cfg, n; psi=ffs.psilim) +spectra = [coil_forcing_modes(cs, forcing_grid, n, m_low, m_high) for cs in coil_sets] ``` -Use `convert_coil_dat_to_h5` / `convert_coil_h5_to_dat` to migrate legacy files, and -`read_coil_h5` / `write_coil_h5` for direct I/O. - -### Rerun snapshot - -When a coil run writes `gpec.h5`, the coil geometry actually used (after shifts/tilts and with -currents) is captured under `Input/RawInputs/Coils/`. The run can then be replayed with -`main_from_h5` even if the original `.dat`/`.h5` files are gone: - -- `--coil-source forcing-modes` (default) recomputes the field from the stored TOML coil config - (needs the original geometry files to be present). -- `--coil-source coils` recomputes the field from the geometry stored in `gpec.h5`, so no - original files are required — and the equilibrium may be overridden on replay. - -An ASCII or HDF5 forcing file is captured raw, with its `normalization` tag, under -`Input/RawInputs/ForcingTerms/`. The replay writes it back out as `_forcing.h5` beside its -output and reads it like any forcing file, so the normalization is redone against the replay's -own equilibrium. - -## Comparison with Fortran GPEC +## Outputs -| Quantity | Julia convention | Fortran GPEC convention | -|----------|-----------------|------------------------| -| Angle range | θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1] | θ_norm ∈ [0, 1], ζ_norm ∈ [0, 1] | -| Coil forcing output | unit-norm = 2π×R×(B_R ∂Z/∂θ_norm − B_Z ∂R/∂θ_norm) | `Phi_x` [Wb] | -| Ratio | Julia output ≡ Fortran Phi_x (directly comparable) | — | +The applied spectrum is written with the response: `PerturbedEquilibrium/ForcingModes/` (`m`, `n`, `amplitude`) and the control-surface spectra `PerturbedEquilibrium/forcing_b`, `forcing_b_root_area` and `forcing_b_area` (see [Perturbed Equilibrium](perturbed_equilibrium.md)). +A forcing file is captured raw with its tag under `Input/RawInputs/ForcingTerms/`, and the coil geometry as used, after shifts and tilts, under `Input/RawInputs/Coils/`, so that `main(["gpec.h5"])` can replay the run without the original files; see the [Scripting API](api.md). -## API Reference +## API reference ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForcingTerms] diff --git a/docs/src/inner_layer.md b/docs/src/inner_layer.md index 1b91bf9810..42da0fdc77 100644 --- a/docs/src/inner_layer.md +++ b/docs/src/inner_layer.md @@ -1,301 +1,149 @@ # Tearing -The `Tearing` module groups the resistive tearing-mode analysis stack: -`InnerLayer` (per-surface inner-layer matching data Δ(Q) for the GGJ and -SLAYER models), `Dispersion` (physics-agnostic complex-plane scan and -contour-intersection root extraction), `CriticalResonantField` (torque-balance -error-field penetration threshold), and the tearing solve itself, -`solve(TearingProblem(ffs; ...), model; config=TearingConfig(...))`, with its -`TearingResult` and HDF5 output. +`Tearing` solves the resistive tearing stability of each rational surface by matching the outer-region Δ′ of the force-free states to the response of a resistive inner layer. +The outer region, ideal MHD everywhere except at the rational surfaces, gives the tearing drive Δ′; the inner layer, where resistivity, two-fluid effects and rotation act, gives the response Δ(Q) as a function of the complex mode frequency; the roots of Δ′ = Δ(Q) are the tearing modes, with their growth rates and real frequencies. -## SLAYER inputs: perpendicular diffusivities +## What it computes -The layer widths and the critical island width of `SLAYER` take perpendicular diffusivities from the kinetic file (`chi_e`, `chi_i`, `chi_phi` profiles in m²/s). -Two formulas want different quantities: the layer width goes through `P_perp`, the diffusion of the total pressure, while the critical island width ``W_d \propto (\chi_\perp/\chi_\parallel)^{1/4}`` is set by the electron heat transport that flattens the electron temperature across the island. -`[SLAYER] P_perp_model` selects how `P_perp` is built: from the perpendicular particle diffusivity `D_⊥` that weights the electron and ion channels by their share of the pressure gradient (`"D_perp"`, the default), from `chi_e` alone (`"chi_perp_e"`, the single-χ convention of Fitzpatrick's TJ code), from `chi_i`, from the toroidal Prandtl number, from its ``c_\beta^2`` floor, or from a whole-plasma confinement time (`"tau_E"`, for a plasma without χ profiles). -`[SLAYER] chi_perp_Dc` sets the electron diffusivity of the critical island width alone, one value for every surface, for a run in which one χ⊥,e cannot serve both formulas; where an interpretive power balance cannot separate the species, treat it as a modeling choice and scan it. -Since the widths scale as ``\chi^{1/4}``, a factor-two change in χ moves them by about 19%, which can be enough to move the resistive-layer overlap cut of [Stability](stability.md) by one surface. +Two layer models are available, selected by `[Tearing] inner_model`: -The DIII-D-like kinetic file of the examples carries analytic profiles, +- `"SLAYER"`, the two-fluid drift-MHD slab layer of Fitzpatrick ([Park (2022)](citations.md#InnerLayer-Module), [Burgess et al. (2026)](citations.md#InnerLayer-Module)), with resistivity, perpendicular and toroidal transport, diamagnetic drifts and E×B rotation taken from the kinetic profiles of the equilibrium; its options are in `[SLAYER]`. +- `"GGJ"`, the toroidal resistive-MHD layer of Glasser, Greene and Johnson ([Glasser, Wang and Park (2016)](citations.md#InnerLayer-Module), [Glasser and Wang (2020)](citations.md#InnerLayer-Module)), whose coefficients E, F, G, H, K are the flux-surface averages stored under `SingularSurfaces/`; its solver options are in `[GGJ]`. -``` -chi_e = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²) -chi_i = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²) -chi_phi = 3/2 − (5/4) (1 − ψ)⁴ -``` +For each surface the `Dispersion` search scans the complex frequency plane, extracts the roots of the dispersion relation and classifies them, and `[Tearing.search]` sets the scan. +With `[Tearing.CriticalResonantField]` the stage also computes, for the `SLAYER` layer, the resonant field at which a static error field overcomes the layer's shielding: the error-field penetration threshold. -which are of the order of power-balance values for DIII-D H-modes but are not fits to a measured discharge and do not resolve the pedestal transport barrier; `chi_e` is raised by a factor of two so that the default deck exercises the layer-overlap cut. +The same GGJ layer serves the `[Match]` stage of the [Workflow](workflow.md), which does the complementary job: instead of finding the free modes it imposes the layer response on the driven force-free states, so that the plasma response of [Perturbed Equilibrium](perturbed_equilibrium.md) is resistive rather than perfectly shielded. -## Tearing solve and layer inputs +## Governing equations -```@autodocs -Modules = [GeneralizedPerturbedEquilibrium.Tearing] -``` +The GGJ layer equations couple the perturbed flux ``\Psi`` and the layer variables ``\Xi``, ``\Upsilon`` in the scaled distance ``x`` from the rational surface (Glasser, Wang and Park (2016), Eq. 11). +A physical growth rate ``\gamma`` maps to the scaled frequency ``Q`` by [`inner_Q`](@ref InnerLayer.inner_Q), the model requires Mercier stability ``D_I < 0`` ([`mercier_di`](@ref InnerLayer.mercier_di)), and at large ``x`` the solutions contain the two power laws of the outer region. +The matching datum ``\Delta`` is the ratio of the large to the small power-law coefficient; imposing the two parities at ``x = 0`` gives the pair ``(\Delta_+, \Delta_-)``: ``\Delta_+`` even in the displacement ``\Xi``, the **interchange** channel, ``\Delta_-`` odd, the **tearing** channel. +Positional pairs are always ordered ``(\Delta_+, \Delta_-)``, and `InnerLayerResponse` names them `interchange` and `tearing`. +The dispersion relation of a surface is then ``\Delta'_s = \Delta_-(Q)`` with ``\Delta'_s`` the surface's diagonal of `SingularSurfaces/Delta_prime_matrix` converted to the layer's normalization (`PerSurface/Delta_prime_norm`); `coupling_mode = "coupled"` solves the determinant of the full matrix over up to `msing_max` surfaces instead. -The module provides an abstract [`InnerLayerModel`](@ref InnerLayer.InnerLayerModel) interface and two -concrete models: the **GGJ (Glasser–Greene–Johnson)** layer documented below, with three -interchangeable solver backends, and the pressureless **SLAYER** Fitzpatrick Riccati layer. -Further inner-layer models (kinetic layers) plug in through the same interface. - -Two papers define the equations and the asymptotic construction; the code and this page cite -them by equation number rather than restating them: - -- **GWP2016** — A. H. Glasser, Z. R. Wang and J.-K. Park, *Computation of - resistive instabilities by matched asymptotic expansions*, Phys. Plasmas - **23**, 112506 (2016). -- **GW2020** — A. H. Glasser and Z. R. Wang, *Asymptotic solutions and - convergence studies of the resistive inner region equations*, Phys. Plasmas - **27**, 012506 (2020). - -## Equations and matching data - -The GGJ layer equations (GWP2016 Eq. 11 ≡ GW2020 Eq. 1) couple the perturbed flux ``\Psi`` and -the layer variables ``\Xi``, ``\Upsilon`` in the scaled distance ``x`` from the rational surface. -Their equilibrium coefficients ``E, F, G, H, K`` (GWP2016 Eq. A8) are held in -[`GGJParameters`](@ref InnerLayer.GGJParameters). A physical growth rate ``\gamma`` maps to the -scaled ``Q`` by [`inner_Q`](@ref InnerLayer.inner_Q) (GWP2016 Eqs. A11, A14–A16), and -[`rescale_delta`](@ref InnerLayer.rescale_delta) returns the matching data to physical units. The -model requires Mercier stability, ``D_I < 0`` ([`mercier_di`](@ref InnerLayer.mercier_di), -GWP2016 Eq. A9), so that at large ``x`` the solutions include two power laws (GW2020 Eq. 49). - -The matching datum ``\Delta`` is the ratio of the large to the small power-law coefficient. -Imposing the two parities at ``x = 0`` gives the pair ``(\Delta_+, \Delta_-)`` of GWP2016 -Eqs. (34)–(35). Parity is that of the displacement ``\Xi`` (GWP2016 Eq. 34; Pletzer, Bondeson & -Dewar, J. Comput. Phys. **115**, 530 (1994), Eqs. 32–33): ``\Delta_+`` is even, -``\Psi(0)=\Xi'(0)=\Upsilon'(0)=0``, the **interchange** channel; ``\Delta_-`` is odd, -``\Psi'(0)=\Xi(0)=\Upsilon(0)=0``, the **tearing** channel. Positional pairs, such as `deltar` in -the matching and `RaySolveResult.Δ`, are always ordered ``(\Delta_+, \Delta_-)``, like the -unknowns ``d_\pm`` of GWP2016 Eq. (38) and Pletzer, Bondeson & Dewar Eqs. 83, 85; -`InnerLayerResponse` names them `interchange` and `tearing`. - -## Numerical method - -!!! note "Benchmark surface used in the figures" - The figures use the ``q = 4`` DIII-D-like surface of - [`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark) - (``S \approx 4.6\times10^{6}``, ``D_I \approx -0.31``). - -### Solver backends - -The `GGJ` model exposes three solvers through the `solver` type parameter of -[`GGJModel`](@ref InnerLayer.GGJModel); all use the same large-``x`` asymptotic basis and return -``\Delta`` in the same convention: +The SLAYER layer is a slab with the Lundquist number ``S``, the temperature ratio ``\tau = T_i/T_e``, the compressibility ``c_\beta``, the perpendicular and toroidal Prandtl numbers ``P_\perp``, ``P_\mathrm{tor}`` and the diamagnetic frequencies ``Q_e``, ``Q_i`` of each surface (all stored under `Tearing/PerSurface/`), solved by a Riccati integration across the layer; its layer thickness ``\delta_s`` is written to `Tearing/LayerWidths/` and feeds the layer-overlap cap of [Stability](stability.md). -| backend | method | robust regime | -|:--------|:-------|:--------------| -| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` | -| `:galerkin` | Hermite-cubic finite-element (real axis) | surface- and ``\arg Q``-dependent; see [Choosing a backend](@ref) | -| `:ray` (default) | rotated contour: inward-marched far-field condition + Hermite-cubic BVP | ``\lvert Q\rvert \sim 500`` on/near the imaginary axis | +### Critical resonant field -Rotation and resistivity scans live on the imaginary-``Q`` axis, where the pseudo-resonance at -``x^2 = -Q^2(G+KF)`` sits on the real-``x`` path. On the ``q = 4`` surface there `:shooting` and -`:galerkin` lose accuracy beyond ``|Q| \sim 1`` and ``|Q| \sim 4`` respectively (the poles of ``\Delta(Q)`` -themselves lie near the real-``Q`` axis): +With `[Tearing.CriticalResonantField]` each `SLAYER` surface gets the critical resonant field for error-field penetration, from the steady-state torque balance of Cole and Fitzpatrick (2006): the electromagnetic torque on the layer balances the viscous restoring torque when -![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q=4 benchmark surface. Each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png) - -### The rotated ray +```math +\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2} + = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2}, +\qquad +\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r}, +``` -`:ray` writes the layer equations as ``dv/dx = M(x)\,v`` in the plain state -``v = (\Psi, \Xi, \Upsilon, \Psi', \Xi', \Upsilon')``. ``M`` is polynomial in ``x``, so -``x = 0`` is an ordinary point and the system continues analytically onto the ray +and balance is lost above ```math -x = e^{i\theta}\, s, \qquad s \in [0, S], \qquad \theta = \tfrac14\arg Q . +\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit} + = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} . ``` -This angle makes the parabolic-cylinder exponent of the outer solutions real and lifts the -pseudo-resonance off the contour; ``\theta = 0`` is a real-axis solve. +``P`` is the layer's `P_tor`, ``S`` its Lundquist number, ``Q_0 = \tau_k n \omega_E`` the natural E×B rotation of the mode (from the kinetic file's `omega_E` or the `omega_E_kHz` override; a surface without one is skipped), and the viscosity integral is taken as ``\hat\kappa = 2/s^2``. +Each surface uses its own stability parameter ``\alpha = S^{-1/3}(\Delta_\mathrm{crit} - r_s\Delta')``, so the torque balance shares its poles with the free modes; the threshold needs a tearing-stable surface, and one with ``\alpha \le 0`` is skipped with a warning. +Penetration is the first resolved saddle-node bifurcation on the path from ``Q_0`` toward the electron fluid at ``Q_e``, above which the shielded branch no longer exists; small numerical extrema and the unshielded zone next to ``Q_e`` are excluded, and a surface without a resolved bifurcation reports NaN. +`br_crit` is normalized by the toroidal field, and `br_crit_tesla = br_crit·|B0|` is the quantity to compare with `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, the total resonant field of the ideal response; the comparison is approximate, since the theory's ``b_r`` is a cylindrical field at the layer. + +## Running it + +From `examples/DIIID-like_SLAYER_example/gpec.toml`, after a `[ForceFreeStates]` solve with `integrator = "Riccati"` and `vac_flag = true` (or `"Galerkin"` with `rpec_flag = true`), which provides the Δ′ matrix, and with a kinetic file on the equilibrium: + +```toml +[Tearing] +inner_model = "SLAYER" # or "GGJ" +coupling_mode = "uncoupled" # one dispersion relation per surface; "coupled" uses the full Δ′ matrix + +[Tearing.search] +scan_mode = "amr" # adaptively refined scan of the complex Q plane +Q_re_range = [-2.0, 2.0] +Q_im_range = [-0.5, 3.0] +nre = 41 +nim = 31 +pole_threshold_adaptive = true +filter_above_poles = true +filter_outside_re = true +polish_roots = true +store_scan = false # true: write Tearing/Scan/Surface_/ + +[Tearing.CriticalResonantField] +n = 2000 # Q samples of the torque balance +store_scan = false + +[SLAYER] +D_c_type = "none" # critical-Δ offset from parallel transport ("none": zero) +P_perp_model = "D_perp" # how the perpendicular Prandtl number is built (below) +``` -![The rotated integration ray in the complex layer-coordinate plane at Q = 500i. The real-axis contour runs through the pseudo-resonance x² = −Q²(G+KF); the ray at θ = arg(Q)/4 = 22.5° clears it.](figures/inner_layer/rotated_ray_contour.png) +In a script, `solve(TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl()), SLAYERModel(); config=TearingConfig(...))` runs the same stage, and `solve_inner(GGJModel(), p, γ)` evaluates one GGJ layer at a growth rate ``\gamma``. -### Hermite-Galerkin BVP +### SLAYER inputs: perpendicular diffusivities -On ``[0, s_m]`` the system is solved with the `:galerkin` backend's Hermite-cubic assembly -(GWP2016 Eq. 32) written on the ray, ``A u_{ss} + e^{i\theta}B u_s + e^{2i\theta}C u = 0``, -with the parity conditions at ``s = 0``. The far end carries the point condition +The layer widths and the critical island width of `SLAYER` take perpendicular diffusivities from the kinetic file (`chi_e`, `chi_i`, `chi_phi` profiles in m²/s). +Two formulas want different quantities: the layer width goes through `P_perp`, the diffusion of the total pressure, while the critical island width ``W_d \propto (\chi_\perp/\chi_\parallel)^{1/4}`` is set by the electron heat transport that flattens the electron temperature across the island. +`[SLAYER] P_perp_model` selects how `P_perp` is built: from the perpendicular particle diffusivity ``D_\perp`` that weights the electron and ion channels by their share of the pressure gradient (`"D_perp"`, the default), from `chi_e` alone (`"chi_perp_e"`, the single-χ convention of Fitzpatrick's TJ code), from `chi_i`, from the toroidal Prandtl number, from its ``c_\beta^2`` floor, or from a whole-plasma confinement time (`"tau_E"`, for a plasma without χ profiles). +`[SLAYER] chi_perp_Dc` sets the electron diffusivity of the critical island width alone, one value for every surface, for a run in which one χ⊥,e cannot serve both formulas; where an interpretive power balance cannot separate the species, treat it as a modeling choice and scan it. +Since the widths scale as ``\chi^{1/4}``, a factor-two change in χ moves them by about 19%, which can be enough to move the resistive-layer overlap cut of [Stability](stability.md) by one surface. +The DIII-D-like kinetic file of the examples carries analytic χ profiles of the order of power-balance values for DIII-D H-modes; they are not fits to a measured discharge, do not resolve the pedestal transport barrier, and `chi_e` is raised by a factor of two so that the default deck exercises the layer-overlap cut. -```math -v(s_m) = U_s + \Delta\, U_b + c_1 E_1 + c_2 E_2 , -``` +## Outputs -where ``U_s, U_b`` are the small and large power-like solutions and ``E_{1,2}`` the -forward-decaying exponential pair; ``\Delta, c_1, c_2`` replace the last node's six unknowns. -Each cell resolves the local WKB rate, with a minimum cell count for small ``|Q|``. - -### Far-field boundary and the inward march - -The boundary data come from the same `inps` Wasow basis as the other backends (GW2020 -Eqs. 3–53), evaluated at the radius ``S`` where its residual is below tolerance -([`pick_smax`](@ref InnerLayer.pick_smax)). At large ``|Q|``, ``S`` lies far beyond ``s_m``, and -the power pair is marched inward by 2-stage Radau IIA modulo the decaying pair, the subspace in -which ``\Delta`` is defined. The decaying pair grows under backward integration at a rate -``\sim\rho s``. Radau's stability function decays as ``2/z``, so steps with every growing mode at -``z > 8`` damp those directions while the slowly varying power pair stays resolved; the step count -grows only like ``\ln(S/s_m)``. Nearer ``s_m`` the march resolves the growing directions, carries -the decaying frame, and projects it out periodically. The damped zone runs in `Complex{Double64}` -because its ill-conditioned implicit solves mix the power pair at a level ``\sim\epsilon\rho S^2``. - -On the real axis with imaginary ``Q`` the switch to resolved steps lands on the -pseudo-resonance, a turning point where the decaying directions cannot be identified; the -default ``\theta = \arg(Q)/4`` keeps it off the contour. - -For matching, `solve_inner_profile` returns ``\Delta`` and the layer-center fields -``\Psi(0), \Xi(0)`` from the rotated solve; ``x = 0`` lies on every contour, so these are the -real-axis values. Fields across the layer on the real ``\psi`` axis are an opt-in diagnostic -(`GGJModel(; real_profiles=true)`, not exposed in the deck): they come from an extra ``\theta = 0`` -solve, while the ``x = 0`` values, and so the penetrated field, stay those of the rotated solve. -The ``\theta = 0`` solve's ``\Delta_\pm`` are checked against the rotated ones (`certΔ`, which -warns on disagreement, as can happen at large ``|Q|`` near the imaginary axis). - -![Inner-layer fields Ψ, Ξ, Υ on the rotated ray for the q=4 surface at Q = 2i. The BVP solution (solid) joins the asymptotic representation (dashed) seamlessly at the match point S = s_m.](figures/inner_layer/solution_profiles.png) - -## Validation and benchmarks - -**Cross-check against the Fortran `rmatch` solver.** At the Glasser & Wang -(2020) Eq. (55) operating point ``Q = 0.1234`` (real) the `:galerkin` backend -reproduces the Fortran `rmatch deltac` solver to ``\sim 10^{-8}``: - -| quantity | `:galerkin` (= Fortran to ``10^{-8}``) | `:ray` | -|:---------|:---------------------------------------|:-------| -| ``\Delta_+`` (interchange) | ``3.698368\times 10^{4}`` | ``3.69789\times 10^{4}`` | -| ``\Delta_-`` (tearing) | ``14.759721`` | ``14.759715`` | - -This point sits near a pole of ``\Delta_+(Q)``, which amplifies every solver's error; hence the -``1.3\times10^{-4}`` `:galerkin`↔`:ray` difference in ``\Delta_+`` against ``4\times10^{-7}`` in ``\Delta_-``. - -**Physical benchmark on the imaginary axis.** On the ``q=4`` rational-surface -benchmark ([`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark), ``S \approx 4.58\times10^6``, -``D_I \approx -0.312``) the `:ray` backend is pinned at ``Q = 500i``, a regime -entirely beyond `:galerkin`: - -| quantity | value at ``Q = 500i`` | -|:---------|:----------------------| -| ``\Delta_+`` (interchange) | ``2.4720 + 13.3540\,i`` | -| ``\Delta_-`` (tearing) | ``0.13750 + 0.74275\,i`` | - -**Convergence and contour invariance.** ``\Delta`` does not depend on the contour angle, so -re-solving with each numerical knob perturbed in turn -([`delta_convergence`](@ref InnerLayer.delta_convergence)) gives an error bar: at ``Q = 500i`` the -worst case is ``\sim 5\times10^{-6}`` for ``\Delta_+`` and ``\sim 6\times10^{-7}`` for ``\Delta_-``. -Across machines and BLAS builds the values reproduce to ``\sim 10^{-5}`` relative, hence the five -significant figures above. - -![Relative change of Δ at Q = 500i under independent perturbations of each numerical knob (contour angle, mesh resolution, series order/radius, march tolerance, BVP edge s_m, decaying-pair separation). The worst-case spread is the reported error bar.](figures/inner_layer/convergence_Sinvariance.png) - -## Choosing a backend - -`:ray` is the **default** — `GGJModel()` constructs `GGJModel{:ray}()` — and the recommended -backend for every ``Q``. It uses the same Hermite-cubic weak form as `:galerkin`; they differ in -where and how the far-field condition is imposed. `:galerkin` imposes it on the real axis through -resonant and extension cells at the radius where the series residual reaches ``10^{-7}``, on a -fixed `nx`-cell grid. `:ray` takes the series data where the residual reaches `smax_tol` = -``10^{-9}`` ([`pick_smax`](@ref InnerLayer.pick_smax)), marches them inward to ``s_m``, where every -exponential family has decayed 45 e-folds from the origin, and resolves the local WKB rate in -every cell on ``[0, s_m]``. Both radii are computed per solve; `S`, `smax_tol` and `sm_fac` -override them, and [`delta_convergence`](@ref InnerLayer.delta_convergence) perturbs them to give -an error bar. - -`:galerkin` takes a few milliseconds per solve, 3–60 times less than `:ray`, but its accuracy -depends on the surface and on ``\arg Q`` even at small ``|Q|``. Against a tightened `:ray` -reference, on the GW2020 Eq. (55) surface it errs by ``7\times10^{-5}`` at ``Q = 0.1`` and -``3\times10^{-3}`` at ``Q = 0.3``; near ``\arg Q = 157.5^\circ`` it errs by order unity at -``|Q| = 0.1`` on that surface and by 0.4 at ``|Q| = 0.3`` on the ``q = 4`` surface. `:ray` held -``2\times10^{-6}`` on all three surfaces tested, for ``|Q|`` from 0.01 to 1. - -The backends take disjoint numerical-knob keywords: pass `GGJModel(solver=:galerkin)` explicitly -when supplying the Galerkin `nx`/`xfac` knobs; passing them to `:ray` throws a `MethodError`. +| Dataset | Meaning | +|---|---| +| `Tearing/PerSurface/rational_psi`, `rational_q`, `m`, `n` | the analyzed surfaces | +| `Tearing/PerSurface/S`, `tau`, `c_beta`, `P_perp`, `P_tor`, `Q_e`, `Q_i`, `tau_R`, `tau_A`, `eta`, `rs`, … | the layer parameters of each surface (SLAYER), and `E`, `F`, `G`, `H`, `K`, `M` (GGJ) | +| `Tearing/PerSurface/Delta_prime_norm`, `delta_prime_conversion`, `Delta_prime_raw` | the Δ′ normalization of the layer and the outer matrix it was matched to | +| `Tearing/Roots/Q_root`, `omega`, `gamma`, `no_root` | the dispersion root of each surface in the lab frame (`frequency_frame` attribute), its real frequency and growth rate, and the flag of a surface with no usable root | +| `Tearing/LayerWidths/delta_s`, `delta_s_abs`, `d_beta` | the complex resistive layer thickness, its magnitude in metres, and the ion drift scale | +| `Tearing/Diagnostics/ValidRoots`, `Poles`, `FilteredRoots` | the root classification of the scan | +| `Tearing/Scan/Surface_/` | the scanned Δ(Q) with `store_scan = true` | +| `Tearing/CriticalResonantField/br_crit`, `br_crit_tesla`, `Q_peak`, `Q0` | the penetration threshold as ``b_r/B_\phi`` and in tesla, the Q of the balance maximum and the natural rotation | -```julia -using GeneralizedPerturbedEquilibrium.InnerLayer +## Numerical notes -p = GGJ.q4_surface_benchmark() # GGJParameters for the q=4 benchmark surface -γ = 500im * GGJ.q0(p) # physical rate placing Q on the imaginary axis at 500i +The `GGJ` model exposes three solvers through the `solver` type parameter of [`GGJModel`](@ref InnerLayer.GGJModel), all returning ``\Delta`` in the same convention: -Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend +| backend | method | robust regime | +|:--------|:-------|:--------------| +| `:shooting` | backward stable shoot from ``X_\mathrm{max}\to 0`` | ``\lvert Q\rvert \ll 1`` | +| `:galerkin` | Hermite-cubic finite-element solve on the real axis | surface- and ``\arg Q``-dependent; see [Choosing a backend](@ref) | +| `:ray` (default) | rotated contour: inward-marched far-field condition and a Hermite-cubic BVP | ``\lvert Q\rvert \sim 500`` on or near the imaginary axis | -# Full result (raw Δ, contour, mesh, nodal solution, boundary conditioning) via solve_ray: -res = solve_ray(p, GGJ.inner_Q(p, γ)) -res.Δ, res.breaks, res.bc_cond # res.Δ = (Δ₊, Δ₋) = (interchange, tearing) -``` +`:ray` continues the layer equations analytically onto the ray ``x = e^{i\theta}s`` with ``\theta = \arg(Q)/4``, which makes the exponent of the outer solutions real and lifts the pseudo-resonance at ``x^2 = -Q^2(G + KF)`` off the contour; it takes the far-field data from the asymptotic series where its residual is below tolerance ([`pick_smax`](@ref InnerLayer.pick_smax)), marches the power pair inward by an implicit Radau scheme that damps the growing exponential directions, and solves the Hermite-cubic boundary-value problem with the parity conditions at ``s = 0``. +``\Delta`` does not depend on the contour angle, so [`delta_convergence`](@ref InnerLayer.delta_convergence) re-solves with each numerical knob perturbed and reports the spread as an error bar. -## Critical resonant field +![Relative error of the :shooting and :galerkin backends against the :ray reference along the imaginary-Q axis, on the q = 4 DIII-D-like benchmark surface; each curve's crossing of the 1% line marks that method's practical reach.](figures/inner_layer/backend_regime_map.png) -With a `[Tearing.CriticalResonantField]` section (`TearingProblem(ffs; critical_resonant_field=CriticalResonantFieldControl())`), each SLAYER surface also gets the -critical resonant field for error-field penetration, from the steady-state torque balance of -Cole and Fitzpatrick, Phys. Plasmas **13**, 032503 (2006). The electromagnetic torque on the -layer balances the viscous restoring torque when (Cole Eq. 61) +Validation: at the real operating point ``Q = 0.1234`` of Glasser and Wang (2020) the `:galerkin` backend reproduces the Fortran `rmatch` solver to ``\sim 10^{-8}`` and `:ray` agrees to ``4\times10^{-7}`` in ``\Delta_-`` (``1.3\times10^{-4}`` in ``\Delta_+``, which sits near a pole); on the ``q = 4`` benchmark surface ([`q4_surface_benchmark`](@ref InnerLayer.GGJ.q4_surface_benchmark), ``S \approx 4.6\times10^6``, ``D_I \approx -0.31``) `:ray` is pinned at ``Q = 500i`` to ``\Delta_+ = 2.4720 + 13.3540\,i``, ``\Delta_- = 0.13750 + 0.74275\,i``, with a knob-perturbation error bar of ``\sim 5\times10^{-6}`` and ``6\times10^{-7}``. +The critical-field torque balance reproduces the Fortran SLAYER's balance to all printed digits when fed the same Δ(Q) table; the thresholds differ through the layer model (Park's drift-MHD layer against Fitzpatrick's) and the selection rule above. -```math -\frac{\operatorname{Im}\hat\Delta(Q)}{|\alpha + \hat\Delta(Q)|^2} - = \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,(b_r/B_\phi)^2}, -\qquad -\hat\kappa = \left[\frac{2}{s}\right]^2 \int_{r_s}^{a} \frac{\mu(r_s)}{\mu(r)}\,\frac{dr}{r}, -``` +![Fortran and Julia inner-layer Δ(Q) and torque balance on the DIII-D-like example.](assets/crf_fortran_benchmark.png) -and balance is lost above (Cole Eq. 62) +### Choosing a backend -```math -\left(\frac{b_r}{B_\phi}\right)^2_\mathrm{crit} - = \max_Q \frac{2P\,(Q_0 - Q)}{S\hat\kappa\,\operatorname{Im}[-1/(\alpha + \hat\Delta(Q))]} . -``` +`:ray` is the default and the recommended backend for every ``Q``: it resolves the local WKB rate in every cell and holds ``\sim 2\times10^{-6}`` for ``|Q|`` from 0.01 to 1 on every surface tested. +`:galerkin` takes a few milliseconds per solve, 3–60 times less than `:ray`, but its accuracy depends on the surface and on ``\arg Q`` even at small ``|Q|``, erring by order unity near ``\arg Q = 157.5^\circ``. +The backends take disjoint numerical keywords: pass `GGJModel(solver=:galerkin)` explicitly when supplying its `nx`/`xfac` knobs. -- **Inputs.** P is each layer's `P_tor` (so χ_φ from the kinetic file, or the scalar - `chi_tor`), S its Lundquist number, and Q0 = τ_k·n·ω_E the natural E×B rotation of the - m/n mode (Cole's ω0; the diamagnetic drifts enter through Q_e and Q_i in Δ̂). Ω_E is the - kinetic file's `omega_E` or the `omega_E_kHz` override; the analysis needs it, and a surface - without one is skipped. -- **Approximations.** As in Fortran `gslayer.f`, the viscosity integral is taken as 1/2 - (κ̂ = 2/s²). -- **Stability α.** Each surface uses its own α = S^(−1/3)(Δ_crit − r_s Δ'), so α + Δ̂ is the - dispersion relation's residual and the torque balance shares its poles with SLAYER's free - modes (Fortran fixes α = 10⁻²). Cole and Fitzpatrick have no Δ_crit; including it treats the - curvature offset as independent of frequency. It is zero unless `D_c_type` is set. Cole's - threshold needs a tearing-stable surface, so a surface with α ≤ 0 is skipped with a warning. -- **Q axis.** Q is on Cole's axis, with the diamagnetic poles at Q_e and Q_i; it runs opposite to the dispersion - scans' Re Q, so the scan uses the layer response at −Q. -- **Threshold.** Penetration is associated with a saddle-node bifurcation (Fitzpatrick, Nucl. Fusion 33, 1049 (1993), Sec. 7.4): - as the field grows, Q moves from Q0 toward the electron fluid at Q_e, where the torque - vanishes, and the threshold is the first resolved saddle-node bifurcation on that path, above - which the shielded branch followed from Q0 no longer exists. Both the rise to and fall from - an interior maximum must exceed 0.1 % of its value within a contiguous retained interval. - This excludes monotonic clipped windows and small numerical extrema; excluded samples do - not connect separate intervals. A later, higher maximum does not replace the first resolved - saddle-node bifurcation. Points where |Δ̂| ≤ α, the - unshielded zone next to Q_e, are not candidates, so the balance spike there is never taken. - With no resolved saddle-node bifurcation no discrete threshold is reported (NaN). Fortran takes the global maximum on a - window instead, which picks that spike. The search assumes Q_i ≤ 0 ≤ Q_e; a surface where - that fails is skipped with a warning. - -Outputs go to `Tearing/CriticalResonantField/` (`br_crit` = b_r/B_φ, `br_crit_tesla`, `Q_peak`, `Q0`, with P in -`PerSurface/P_tor`, and with `store_scan` the per-surface `Scan/Surface_/{Q, balance, Delta}`). - -`br_crit` is normalized by the toroidal field, as in Cole and Fitzpatrick (2006, Eqs. 61–68) and -Fitzpatrick, Phys. Plasmas **30**, 092512 (2023), Eq. 120. `br_crit_tesla` = `br_crit`·|B0|, with -the layer's B0, is for comparison with `SingularCoupling/resonant_area_weighted_field`, the total -resonant field with the ideal plasma response, which Park et al., Phys. Rev. Lett. **99**, 195003 -(2007) identify as the field that opens islands. The comparison is approximate: the theory's b_r is -a cylindrical field at the layer, and no reference maps it onto the area-weighted toroidal field. -`penetrated_area_weighted_field` already includes the layer response and is not the matching quantity. - -### Benchmark against Fortran - -On the DIII-D-like example (2/1–6/1, n = 1, P = 1, κ̂ = 2/s², Q ∈ [−5, 5] with 20001 -points), both codes were given Fortran's normalized inputs and Fortran's fixed α = 10⁻². - -- **Torque balance.** Fed Fortran's own Δ(Q) table, the Julia routine reproduces Fortran's - balance to all printed digits on every surface, so the thresholds differ only through the - selection rule above. -- **Layer model.** Fortran solves Park's drift-MHD layer and Julia Fitzpatrick's, so Δ(Q) - differs by a median 24–78 % between Q_i and Q_e on 2/1, 3/1, 4/1 and 6/1 (7 % on 5/1). - -The dots in the figure mark each code's balance maximum between Q_e and Q0 at that fixed α, -not the fold threshold above. - -![Fortran vs Julia inner-layer Δ(Q) and torque balance](assets/crf_fortran_benchmark.png) +```julia +using GeneralizedPerturbedEquilibrium.InnerLayer + +p = GGJ.q4_surface_benchmark() # GGJParameters for the q = 4 benchmark surface +γ = 500im * GGJ.q0(p) # physical rate placing Q on the imaginary axis at 500i +Δ = solve_inner(GGJModel(), p, γ) # Δ.interchange, Δ.tearing with the default :ray backend +``` ## API Reference +### Tearing + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Tearing] +``` + ### InnerLayer ```@autodocs diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 091492b543..8ba011db09 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -1,217 +1,110 @@ -# KineticForces Module - -Kinetic torque and energy calculations for perturbed equilibria. -Implements neoclassical toroidal viscosity (NTV) from the PENTRC formulation -([Logan & Park, 2013](citations.md); [Logan, 2015](citations.md)). - -## Kinetic profile file formats - -Kinetic profiles are read by `read_kinetic_file`, which dispatches on the file -extension: - - - **HDF5 (`.h5`/`.hdf5`)** — the GPEC kinetic schema (recommended). Datasets - at the file root: required `psi` (normalized poloidal flux), `n_e`, `T_i`, - `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional - `omega_tor`, `chi_e` and `chi_i` (electron and ion perpendicular heat diffusivities - ``\chi_{\perp,e}``, ``\chi_{\perp,i}``), `chi_phi` (toroidal momentum diffusivity - ``\chi_\phi``), and `zeff` - (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset - carries a `units` attribute; the root carries `schema_version` and - `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s, - diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`) - are named per-species densities for multi-ion runs; names outside `n_*` - are reserved for future schema fields and ignored. Write files with - `write_kinetic_h5` (which round-trips the per-species densities). - - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table - `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward - compatibility. Header rows are skipped. - -The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`, -when present, are carried for the resistive-layer (SLAYER) analysis and ignored -here. That analysis takes ``Z_\mathrm{eff}(\psi)`` from the `zeff` dataset when the file -has one, otherwise from `n_i/n_e` with the kinetic data's `impurity_charge`, otherwise -from its `zeff_default`. Supply `zeff` whenever the impurity content is not a -single species fixed by `n_i/n_e`. - -The plasma's species composition (main ions, impurity charge and mass), the defaults that stand -in for absent χ and Z_eff profiles, the confinement time and the resistivity closure belong to the -equilibrium's kinetic data, not to any one physics stage, so KineticForces and the tearing layers -see one plasma. In a script, set them where the profiles are attached: +# Kinetic Forces -```julia -attach_kinetic_profiles!(equil, "kinetic.h5"; ion_species=[IonSpecies(; z=1, m=2, fraction=0.5), - IonSpecies(; z=1, m=3, fraction=0.5)], impurity_charge=6, impurity_mass=12) -``` +`KineticForces` computes the neoclassical toroidal viscosity (NTV) of a perturbed equilibrium: the toroidal torque and the kinetic energy that a non-axisymmetric field exerts on the plasma through the drift-kinetic response of its particles. +It implements the PENTRC formulation of [Park et al. (2009)](citations.md#Kinetic-Forces), [Logan et al. (2013)](citations.md#Kinetic-Forces) and [Park and Logan (2017)](citations.md#Kinetic-Forces). -or copy loaded data with new settings, `KineticProfileData(equil.kinetic; impurity_charge=6)`. -The default is one deuterium main ion carrying all of `n_i`, with a carbon impurity. In a deck -they are `[Equilibrium]` keys: `[[Equilibrium.ion_species]]`, `impurity_charge` and `impurity_mass`. +## What it computes -## Multi-ion runs +A non-axisymmetric perturbation ``\delta B`` breaks the toroidal symmetry that conserves the canonical angular momentum of each particle. +The resulting non-ambipolar radial transport, dominated by trapped particles whose bounce and precession motion resonates with the E×B rotation, carries toroidal momentum out of the plasma: the NTV torque ``T_\phi``, which brakes or spins up the rotation, together with a kinetic contribution ``\delta W_k`` to the perturbed energy that stabilizes or destabilizes the response. +Both are halves of one complex quantity, ``T = T_\phi + 2in\,\delta W_k``, and the module integrates it over the plasma for one or more torque `method`s. -A plasma may declare an arbitrary list of main-ion species; the NTV is computed -per species under one shared full-composition ``Z_\mathrm{eff}`` and summed -(``\tau = \sum_s \tau_s``) over the main ions, the quasineutrality-closing -impurity, and (with `electron = true`) the electrons. This applies to both NTV -paths: the post-PE ψ-quadrature diagnostic and the self-consistent -kinetic matrices of a `kinetic = "KineticForces"` force-balance solve. +The torque enters the pipeline in two ways: -```toml -[Equilibrium] -kinetic_file = "kinetic.h5" # n_i column/dataset = TOTAL main-ion density -impurity_charge = 6 # impurity charge (closes quasineutrality) -impurity_mass = 12 # impurity mass - -[[Equilibrium.ion_species]] -z = 1 -m = 2 -fraction = 0.5 # this species' share of the total n_i -[[Equilibrium.ion_species]] -z = 1 -m = 3 -density = "n_T" # or: explicit n_* dataset from the HDF5 kinetic file +| Workflow | Deck | What it gives | +|---|---|---| +| Torque of a perturbed equilibrium | `[KineticForces]` after `[PerturbedEquilibrium]` | ``T(\psi)`` integrated over the stored response; the stage of the [Workflow](workflow.md) | +| Self-consistent kinetic solve | `[ForceFreeStates] kinetic = "KineticForces"` | the drift-kinetic energy matrices replace the ideal ones inside the stability solve, so the response itself carries the torque (see [Stability](stability.md) and `TorqueResponse/` on the [Perturbed Equilibrium](perturbed_equilibrium.md) page) | -[KineticForces] -electron = true # add electron NTV in addition to the ion species -``` +Both need the kinetic profiles of `[Equilibrium] kinetic_file` (n_e, T_i, T_e, ω_E, and n_i when it differs from n_e) and the ion composition of `[[Equilibrium.ion_species]]`; the file schema is on the [Equilibrium](equilibrium.md) page. -The composition describes the plasma, so it sits on the equilibrium's kinetic data beside -`kinetic_file`, where every kinetic consumer (NTV, the tearing layers) reads it. - -Each species sets exactly one of `fraction` (share of the file's total `n_i`) -or `density` (a named `n_*` dataset). An all-fraction list must sum to 1; in a -mixed list the impurity content is set by the file's `n_i`/`n_e` deficit, not -by a fraction shortfall, and fractions may sum below (never above) 1. Every -main-ion charge must satisfy `z < impurity_charge`. A single main ion carrying all of `n_i` -runs the single-species NTV (with `electron = true` still adding the electron species — the -electron flag, a `KineticForcesControl` setting, always means *in addition to* the ions). Without -an `ion_species` list the plasma is one deuterium main ion. - -The summed total is written to `KineticForces//` exactly as in a single-ion -run, and each species' own contribution to `KineticForces/PerSpecies/