Add fusion particle production diagnostic - #6954
Merged
dpgrote merged 110 commits intoAug 28, 2026
Merged
Conversation
for more information, see https://pre-commit.ci
Contributor
|
Thanks Dave for this PR. It looks good to me, but because it touches so many files it's a good idea to also have it reviewed by someone else. One question. What happens on a restart? |
Member
Author
|
The PR #6966 was added in association. This allows new Multifabs to be added from Python and be included in the diagnostic output, taking advantage of the new capability added here of allowing any MultiFab to be written out. |
ax3l
reviewed
Jun 23, 2026
ax3l
added a commit
that referenced
this pull request
Jun 23, 2026
This adds a new call back to the Python interface. This call back, `allocdata`, is called during initialization, both from scratch and from restart. New MultiFabs should be allocated at this point. This is associated with the PR #6954 which allows writing of any MultiFabs to the diagnostics and to checkpoint files. The new MultiFabs would need to be added in the `allocdata` call back to be included in the diagnostics. --------- Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
This PR does some clean up in `AddPlasmaFlux`: - It removes use of `ppos` which was doing somewhat arbitrary conversions between `Real` and `ParticleReal` (this fixes a bug in `RSPHERE` where `pos` was being set instead of `ppos`) - It changes `pu` to type `XDim3` for the same reason, avoiding arbitrary conversions between `Real` and `ParticleReal` - Since it is no longer used, the `PDim3` struct is removed - Rename `u` to `gamma_beta` for clarity (maybe do the same in other routines?) - Move the geometry `ifdef`s into `insideBoundsInclusive` to reduce code duplication A general note is that this is setting up particle quantities, but most of the variables are type `Real`, with the conversion to `ParticleReal` only at the end of the loop when the particle arrays are being set. Is this an oversight or intended? --------- Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
Add new paper that uses WarpX to the highlights section. https://iopscience.iop.org/article/10.1088/1741-4326/ae96c0 Signed-off-by: Roelof Groenewald <rgroenewald@realtafusion.com>
## Stack of PRs (2/3) 1. BLAST-WarpX#7142 — `Python: Fix WarpX Singleton Lifetime and Leaking Statics` — **merged** 2. **This PR** — `Python: Clear Input State in warpx.finalize()` 3. BLAST-WarpX#7144 — `Tests: pytest Unit Tests for Charge and Current Deposition`, which is what actually exercises this ## Summary Even with the WarpX singleton properly finalizable (BLAST-WarpX#7142), running more than one simulation in a Python process still does not work: the input deck lives in module-level `Bucket` objects (`pywarpx.warpx`, `pywarpx.geometry`, `pywarpx.particles`, ...) that accumulate across simulations, so a second simulation inherits everything the first one set. * User-facing: **`warpx.finalize()`** (`Python/pywarpx/WarpX.py`) now also restores all module-level buckets and lists to their construction-time defaults, after tearing down WarpX and AMReX. No new API is added. The compiled `warpx_pybind_*` module deliberately stays loaded: multiple AMReX/WarpX geometries still cannot coexist in one process, so the dimensionality remains fixed for the lifetime of the process. To run another simulation afterwards, construct new PICMI objects. * Internal: **`Bucket.set_default_attrs()`** factors out the loop that applies the construction-time defaults, shared by `__init__` and the new **`Bucket.clear()`**. The defaults are held as a deep-copied snapshot, so that mutating a mutable default — for example `pywarpx.particles.species_names.append(...)`, which `picmi.Species` does — cannot leak into the next simulation. ## Details Per the review discussion below, the clearing lives in `finalize()` instead of a separate `reset()`, and it is not optional: what is cleared is the input deck, not results. Once the C++ side is gone, a half-populated deck is not something a script can act on, it is only a way for the settings of one simulation to leak into the next one. Two notes on the implementation: * it is in `WarpX.finalize()` and not in `LibWarpX.finalize()`, which is the function registered with `atexit`. There is nothing to gain from resetting Python state while the interpreter is shutting down, and doing so would run imports during teardown. * `picmi.Simulation.finalize()` needs no change. A `Simulation` builds its deck from PICMI objects, and the grid writes `geometry.dims`, `geometry.prob_lo` and `geometry.prob_hi` when it is constructed rather than in `grid_initialize_inputs()`, so a further simulation means new PICMI objects in any case. ## Testing `Bucket.clear()` was checked in isolation for the mutable-default aliasing case, and `warpx.finalize()` over three consecutive PICMI simulation cycles: species names, geometry and the dynamic sub-buckets all come back empty each time. End to end, the pytest suite in BLAST-WarpX#7144 calls `warpx.finalize()` between every test and passes 7/7, parametrizing over three particle shape orders and two current deposition algorithms, each of which needs a genuinely fresh simulation. Verified there that all seven tests really do build and tear down a full WarpX instance. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
This PR adds a semi-implicit Darwin field solver.
## Overview of the model
In the Darwin approximation, the transverse displacement current is
dropped from the Maxwell-Ampere equation, which removes light waves from
the system while retaining the inductive (low-frequency magnetic)
physics. In other words, the Maxwell-Ampere equation is replaced by:
$$\nabla\times\mathbf{B} = \mu_0\left(\mathbf{J} +
\epsilon_0\frac{\partial\mathbf{E}_{irr}}{\partial t}\right)$$
where the electric field has been decomposed into its irrotational and
solenoidal part:
$$\mathbf{E} = \mathbf{E}_{irr} + \mathbf{E}_{sol} \qquad
\mathbf{E}_{irr} = -\nabla\phi \qquad \mathbf{E}_{sol} = -
\frac{\partial\mathbf{A}}{\partial t} \qquad \mathbf{B} =
\nabla\times\mathbf{A}$$
## User interface
The solver is selected as an evolve scheme, on top of an electrostatic
solver:
```
algo.evolve_scheme = semi_implicit_darwin
warpx.do_electrostatic = labframe # required: the Darwin scheme adds the inductive
# field on top of an electrostatic solve
algo.maxwell_solver = yee # required (default)
amrex_gmres.relative_tolerance = 1.e-4 # magnetostatic (GMRES) solve controls
amrex_gmres.max_iterations = 1000
```
From PICMI:
```python
simulation.evolve_scheme = picmi.SemiImplicitDarwinEvolveScheme(
linear_solver=picmi.GMRESLinearSolver(relative_tolerance=1e-4)
)
simulation.solver = picmi.ElectrostaticSolver(...)
```
Unlike a pure electrostatic run, setting this scheme does *not* disable
the electromagnetic solver, since the magnetic field is still evolved.
The new parameters are documented in `Docs/source/usage/parameters.rst`.
## The algorithm, and how it maps onto `OneStep`
### Time staggering and updates during one timestep
| Quantity | Time level |
|---|---|
| position $\mathbf{x}$ | Updated $n \rightarrow n+1$ (integer) |
| momentum $\mathbf{u}$ | Updated $n-1/2 \rightarrow n+1/2$
(half-integer) |
| $\phi$, $\mathbf{E}_{irr} = -\nabla\phi$ | Computed from scratch at
$n$ |
| $\mathbf{A}$, $\mathbf{B}$ | Updated $n-1/2 \rightarrow n+1/2$ |
| $\mathbf{E}_{sol}$ | Computed at $n$ from $\mathbf{A}$ |
| current $\mathbf{J}$ | $n$ (integer — **not** $n\pm1/2$) |
The time-centring of $\mathbf{J}$ is the main departure from the
standard explicit scheme: $\mathbf{J}$ is needed at the *same* time as
the electric field being solved for, so it is deposited at the integer
time $n$ from the time-centred velocity $\mathbf{u}^n =
(\mathbf{u}^{n+1/2} + \mathbf{u}^{n-1/2})/2$. This is what forces the
predictor/corrector structure, since $\mathbf{u}^{n+1/2}$ is not yet
known when $\mathbf{J}^n$ must be deposited.
### The discretized field equation
Writing the vector potential increment as $\Delta\mathbf{A}^n =
\mathbf{A}^{n+1/2} - \mathbf{A}^{n-1/2}$, so that $\mathbf{E}_{sol}^n =
-\Delta\mathbf{A}^n/\Delta t$, and using $\mathbf{A}^n =
\mathbf{A}^{n-1/2} + \Delta\mathbf{A}^n/2$ in the elliptic equation
above:
$$-\nabla^2\Delta\mathbf{A}^n = 2\mu_0\mathbf{J}_{sol}^n +
2\nabla^2\mathbf{A}^{n-1/2}$$
Taking the curl eliminates the irrotational part of the current (since
$\nabla\times\mathbf{J}_{sol} = \nabla\times\mathbf{J}$), which avoids
ever having to perform the Helmholtz decomposition explicitly, and
introduces $\mathbf{B}^{n-1/2} = \nabla\times\mathbf{A}^{n-1/2}$. An
auxiliary field $\mathbf{Z}^n$ is then introduced through
$\Delta\mathbf{A}^n \equiv \nabla\times\mathbf{Z}^n$, which guarantees
by construction that the increment is divergence-free, i.e. that the
Coulomb gauge is preserved. This yields the single fourth-order equation
that is actually solved:
$$\nabla^4\mathbf{Z}^n +
\nabla\times\left(\boldsymbol{\chi}\,\nabla\times\mathbf{Z}^n\right) =
2\mu_0\nabla\times\hat{\mathbf{J}}^n + 2\nabla^2\mathbf{B}^{n-1/2}$$
Here $\hat{\mathbf{J}}^n$ is the *predicted* current and
$\boldsymbol{\chi}$ is the mass-matrix (susceptibility) tensor; both are
explained next.
### Predictor/corrector, and where the mass matrix comes in
The Boris push is affine in the electric field, so writing
$\boldsymbol{\Theta} = \frac{q\Delta t}{2m}\mathbf{B}^{n-1/2}$ for the
usual Boris rotation vector, the velocity update splits *exactly* into a
part driven by the (known) electrostatic field and a part driven by the
(unknown) inductive field:
$$\mathbf{u}^{n+1/2} + \mathbf{u}^{n-1/2} = \mathbf{u}^{\ast} +
\mathbf{u}^{\dagger}$$
$$\mathbf{u}^{\ast} = 2\,\frac{\mathbb{I} - \boldsymbol{\Theta}\times +
\boldsymbol{\Theta}\boldsymbol{\Theta}^T}{1+\Theta^2}\left(\mathbf{u}^{n-1/2}
+ \frac{q\Delta t}{2m}\mathbf{E}_{irr}^n\right), \qquad
\mathbf{u}^{\dagger} = 2\,\frac{\mathbb{I} - \boldsymbol{\Theta}\times +
\boldsymbol{\Theta}\boldsymbol{\Theta}^T}{1+\Theta^2}\left(\frac{q\Delta
t}{2m}\mathbf{E}_{sol}^n\right)$$
The predictor push computes $\mathbf{u}^{\ast}$, which is all that is
knowable before the field solve. Taking the charge-weighted moments of
the two terms and using $\mathbf{u}^n =
(\mathbf{u}^{n+1/2}+\mathbf{u}^{n-1/2})/2$ gives the current at the
integer time,
$$2\mathbf{J}^n = \mathbf{J}^{\ast} + \mathbf{J}^{\dagger}, \qquad
\mathbf{J}^{\ast} = 2\hat{\mathbf{J}}^n = \sum_p q_p w_p
\mathbf{u}_p^{\ast}\,S(\mathbf{x}_p^n)$$
**This is where the mass matrix enters.** Because $\mathbf{u}^{\dagger}$
is *linear* in $\mathbf{E}_{sol}^n$, the moment $\mathbf{J}^{\dagger}$
can be written in closed form as a linear operator acting on the
still-unknown field, rather than requiring an iterative solve:
$$\mathbf{J}^{\dagger} = \frac{\Delta
t}{\mu_0}\boldsymbol{\chi}\,\mathbf{E}_{sol}^n, \qquad
\boldsymbol{\chi}_{ii'} = \sum_s\sum_{p\in s}\frac{\mu_0 q_p^2
w_p}{m_p}\,\frac{\mathbb{I} +
\boldsymbol{\Theta}_p\boldsymbol{\Theta}_p^T -
\boldsymbol{\Theta}_p\times}{1+\Theta_p^2}\,S_i(\mathbf{x}_p)S_{i'}(\mathbf{x}_p)$$
$\boldsymbol{\chi}$ is the mass matrix (susceptibility): a
grid-point-to-grid-point tensor accumulated over particles in the same
sweep as the current deposition. It is exactly the linear response of
the deposited current to the field the solve is about to produce.
Finally, substituting $\mathbf{E}_{sol}^n = -\Delta\mathbf{A}^n/\Delta t
= -\nabla\times\mathbf{Z}^n/\Delta t$ turns the response term into an
operator on the unknown $\mathbf{Z}^n$,
$$\mathbf{J}^{\dagger} =
-\frac{1}{\mu_0}\boldsymbol{\chi}\,\nabla\times\mathbf{Z}^n$$
so that $2\mu_0\nabla\times\mathbf{J}^n =
2\mu_0\nabla\times\hat{\mathbf{J}}^n -
\nabla\times(\boldsymbol{\chi}\,\nabla\times\mathbf{Z}^n)$. Moving the
response term to the left-hand side produces the
$\nabla\times(\boldsymbol{\chi}\,\nabla\times\mathbf{Z}^n)$ term of the
field equation above, leaving the system **linear** — hence
"semi-implicit": one linear solve per step, with no Newton iteration. In
the code, that term is evaluated by `ApplySusceptibility`, which applies
the deposited mass matrices to $\nabla\times\mathbf{Z}$ with the
appropriate $\mu_0/\Delta t$ scaling.
### `SemiImplicitDarwin::OneStep`
At the start of a step the fields hold $\mathbf{E} = -\nabla\phi^n$
(left over from the electrostatic solve at the end of the previous step)
and $\mathbf{B}^{n-1/2}$; the particles hold $\mathbf{x}^n$ and
$\mathbf{u}^{n-1/2}$. The step then proceeds as:
1. **Predictor (electrostatic) push.** `PushP` applies a full Boris push
with the electrostatic field only, giving $\hat{\mathbf{u}}^{n+1/2}$.
`PrepareVelocitiesForCurrentDeposition` then replaces `u` by the
time-centred average $(\hat{\mathbf{u}}^{n+1/2} + \mathbf{u}^{n-1/2})/2$
and stashes $\hat{\mathbf{u}}^{n+1/2}$ in `u_n` for later.
2. **Current and mass-matrix deposition**
(`AccumulateCurrentAndSusceptibility`). $\hat{\mathbf{J}}^n$ is
deposited from those time-centred velocities, and $\boldsymbol{\chi}$
from the same `DepositMassMatrices` machinery already used by the
implicit EM solvers.
3. **Magnetostatic solve.** `CalculateSourceVector` builds the
right-hand side $2\mu_0\nabla\times\hat{\mathbf{J}}^n +
2\nabla^2\mathbf{B}^{n-1/2}$, and `amrex::GMRES` solves the fourth-order
equation for $\mathbf{Z}^n$. The operator $\nabla^4 +
\nabla\times(\boldsymbol{\chi}\,\nabla\times\,\cdot\,)$ is applied
matrix-free by `ComputeRHS`.
4. **Inductive field.** `ComputeInductiveEfromdA` sets
$\mathbf{E}_{sol}^n = -\Delta\mathbf{A}^n/\Delta t =
-\nabla\times\mathbf{Z}^n/\Delta t$, overwriting `Efield_fp`, which from
here on holds the inductive component only.
5. **Corrector (solenoidal) push.** The velocities are zeroed and
`PushP` is called again with the inductive field alone, which isolates
$\delta\mathbf{u}(\mathbf{E}_{sol}^n)$; `FinishVelocityUpdate` then adds
back $\hat{\mathbf{u}}^{n+1/2}$ from `u_n` to form the complete
$\mathbf{u}^{n+1/2}$. Positions are advanced to $\mathbf{x}^{n+1}$.
6. **Magnetic field update.** `EvolveB` advances $\mathbf{B}^{n-1/2}$ to
$\mathbf{B}^{n+1/2}$ with $\partial\mathbf{B}/\partial t =
-\nabla\times\mathbf{E}$. Since `Efield_fp` holds
$-\nabla\times\mathbf{Z}^n/\Delta t$, this is equivalent to
$\mathbf{A}^{n+1/2} = \mathbf{A}^{n-1/2} + \nabla\times\mathbf{Z}^n$
followed by $\mathbf{B} = \nabla\times\mathbf{A}$, without ever storing
$\mathbf{A}$.
The electrostatic solve for $\phi^{n+1}$ is then performed by the
existing electrostatic branch of `WarpX::Evolve`, which leaves
`Efield_fp` holding $-\nabla\phi^{n+1}$ ready for the next step.
## Code structure
- **`Source/FieldSolver/ImplicitSolvers/SemiImplicitDarwin.{H,cpp}`** —
the new solver, an `ImplicitSolver` subclass. `OneStep` drives the
sequence above; `ComputeRHS` evaluates the linear operator for GMRES.
- **`Source/NonlinearSolvers/LinearFunctionMF.H`** — matrix-free linear
operator that lets `amrex::GMRES` call back into `ComputeRHS`.
- **Finite-difference operators** — `ComputeCurlB.cpp` (curl of a
B-staggered field, output on E staggering), `ComputeLaplacian.cpp`
(scalar and vector Laplacian, plus a single-pass vector bi-Laplacian),
and the fourth-derivative stencils
`Dxxxx`/`Dyyyy`/`Dzzzz`/`Dxxyy`/`Dyyzz`/`Dxxzz` in
`CartesianYeeAlgorithm.H`.
- **`Source/Fields.H`** — new `dA_fp` vector field holding the
vector-potential increment over a step.
- **`Source/Parallelization/WarpXComm.cpp`** —
`WarpX::SyncMassMatrices()`, the boundary summation of the deposited
mass matrices.
-
**`Source/Diagnostics/ComputeDiagFunctors/DarwinEfieldFunctor.{H,cpp}`**
— since `Efield_fp` holds only the electrostatic component at diagnostic
time, this functor reconstructs the full field as `E_es + (-dA/dt)` for
output; wired up in `FullDiagnostics.cpp`.
- **`Source/WarpX.cpp` / `WarpXEvolve.cpp`** — the new
`Semi_Implicit_Darwin` evolve scheme, its input validation, and the
separation of the hybrid-PIC branch from the electrostatic branch of the
PIC loop (Darwin needs the electrostatic solve without the accompanying
B-field reset).
- **`Python/pywarpx/picmi.py`** — `SemiImplicitDarwinEvolveScheme`.
## Tests
The existing `Examples/Tests/ohm_solver_em_modes` directory has been
generalised into `Examples/Tests/magnetized_plasma_modes`, with a single
input script that runs the same magnetized-plasma EM-mode setup with
either solver (`--darwin` / `--ohm`). Two new tests were added there:
- `test_1d_darwin_solver_em_modes_picmi`
- `test_2d_darwin_solver_em_modes_es_picmi` (Darwin combined with the
effective-potential electrostatic solver)
The implementation was verified by reproducing the dispersion of left-
and right-hand circularly polarized Alfvén waves propagating parallel to
an applied magnetic field (see comment further below).
---------
Signed-off-by: roelof-groenewald <regroenewald@gmail.com>
Signed-off-by: Roelof Groenewald <rgroenewald@realtafusion.com>
Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Cursor <cursoragent@cursor.com>
…T-WarpX#7156) `FieldDiagnostic._get_diagnostic_data()` maps each entry of `data_list` onto a `fields_to_plot` name through a chain of `elif` branches. A name matching none of them was silently dropped. Since BLAST-WarpX#7025 ("Allow any MultiFab to be written to the diagnostics") the C++ diagnostics resolve any name present in the `MultiFabRegister`. That covers fields registered from Python and internal fields that have no short alias in the Python branches -- e.g. "hybrid_current_fp" or "vector_potential_fp". (Fields that do have an alias, such as "Pe" and "Te" added in BLAST-WarpX#7081, are matched earlier and never reach the new branch.) Without a fallthrough the name never survives the Python layer, so the diagnostic silently produces no such field and gives no indication why. Add a final else that forwards the name unchanged. Validation stays in C++, where `FullDiagnostics` already raises a descriptive error for a name that is neither a known field type nor in the register, so a genuine typo is still reported -- and now reported with a message, rather than silently ignored. Verified via the 3D ohm-solver cylinder-compression test with an unaliased registered field name in `data_list` runs to completion and the field appears in the plotfile Header.
The `flag_info_face` iMultiFab used by the ECT solver stored bare integer codes (0, 1, 2 and -1) whose meaning had to be looked up in a comment at every use site. Introduce `FaceInfo::Flag` in `WarpXFaceInfoBox.H`, with one enumerator per category (`extended`, `available`, `intruded`, `bck_stabilized`), and use it in the initialization, in the face extensions and in the ECT B push. The enumeration is unscoped and has a fixed `int` underlying type, because the values are read from and written to an `amrex::iMultiFab` and a scoped enumeration would require a cast at every comparison. This is a pure readability change: the stored values are unchanged. Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
…re_output.py (BLAST-WarpX#7185) This is a small fix to `update_benchmarks_from_azure_output.py` to explicitly add an end-of-line character to the last line of the JSON files written out. Whenever I use this script, it was always showing a change to the last line due to the lack of the end-of-line character there. Having the end-of-line character is much cleaner (and is the POSIX standard).
Automated via .github/workflows/weekly_update.yml. --------- Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com> Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>
This should be merged after - BLAST-WarpX#6293 In this PR, the energy conserving PIC scheme is recovered for the Darwin implementation with `direct` current deposition. This is done by forcing the electrostatic particle push to use the "Galerkin" interpolation method unless the user specifically sets the run to use "momentum-conserving" mode. --------- Signed-off-by: roelof-groenewald <regroenewald@gmail.com> Signed-off-by: Roelof Groenewald <rgroenewald@realtafusion.com> Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Bumps [github/codeql-action](https://github.com/github/codeql-action) from 4.37.7 to 4.37.8. <details> <summary>Release notes</summary> <p><em>Sourced from <a href="https://github.com/github/codeql-action/releases">github/codeql-action's releases</a>.</em></p> <blockquote> <h2>v4.37.8</h2> <p>No user facing changes.</p> </blockquote> </details> <details> <summary>Changelog</summary> <p><em>Sourced from <a href="https://github.com/github/codeql-action/blob/main/CHANGELOG.md">github/codeql-action's changelog</a>.</em></p> <blockquote> <h2>4.37.8 - 21 Aug 2026</h2> <p>No user facing changes.</p> </blockquote> </details> <details> <summary>Commits</summary> <ul> <li><a href="https://github.com/github/codeql-action/commit/db488ddef3bf6cb639b32c2e9a7c0a7ea8271d28"><code>db488dd</code></a> Merge pull request <a href="https://redirect.github.com/github/codeql-action/issues/4102">#4102</a> from github/update-v4.37.8-9ee088e13</li> <li><a href="https://github.com/github/codeql-action/commit/1845f5ba8b4057590f49ee8e246c95ef2ba4b53f"><code>1845f5b</code></a> Update changelog for v4.37.8</li> <li><a href="https://github.com/github/codeql-action/commit/9ee088e13615f8d1eaef4766f9dde95d3356a8f6"><code>9ee088e</code></a> Merge pull request <a href="https://redirect.github.com/github/codeql-action/issues/4080">#4080</a> from github/henrymercer/studious-giggle</li> <li><a href="https://github.com/github/codeql-action/commit/1aef003397c876c0ab5bd118e1b1f34c175622e9"><code>1aef003</code></a> Address review feedback on overlay disk flags</li> <li><a href="https://github.com/github/codeql-action/commit/508b83bc415e8df76ce8ea08c0cf42c2529ebc63"><code>508b83b</code></a> Merge main into overlay minimum disk feature branch</li> <li><a href="https://github.com/github/codeql-action/commit/d97b3428e8eebbb1810cf454d6397886d136b4ba"><code>d97b342</code></a> Merge pull request <a href="https://redirect.github.com/github/codeql-action/issues/4098">#4098</a> from github/mbg/permission-error-as-configuration-error</li> <li><a href="https://github.com/github/codeql-action/commit/47fa6222231b12097f83215dd7a6b4a0915841fd"><code>47fa622</code></a> Make <code>EACCES</code> a <code>ConfigurationError</code></li> <li><a href="https://github.com/github/codeql-action/commit/45693cc6882bb175b58a06818c91876e201037c7"><code>45693cc</code></a> Refactor <code>ENOSPC</code> check into <code>isDiskConfigurationError</code> function</li> <li><a href="https://github.com/github/codeql-action/commit/c2fd8f54d19fa46c94ed79cb92e6dd6606d61762"><code>c2fd8f5</code></a> Merge pull request <a href="https://redirect.github.com/github/codeql-action/issues/4081">#4081</a> from github/mario-campos/version-cache-to-disk</li> <li><a href="https://github.com/github/codeql-action/commit/c56f48e9bd458a387eb68a68534459e503e56b17"><code>c56f48e</code></a> Log unexpected conditions during caching CLI output</li> <li>Additional commits viewable in <a href="https://github.com/github/codeql-action/compare/v4.37.7...v4.37.8">compare view</a></li> </ul> </details> <br /> [](https://docs.github.com/en/github/managing-security-vulnerabilities/about-dependabot-security-updates#about-compatibility-scores) Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) --- <details> <summary>Dependabot commands and options</summary> <br /> You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot show <dependency name> ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself) </details> Signed-off-by: dependabot[bot] <support@github.com> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
## What this fixes
`setMassMatricesKernels` built the mass-matrix (susceptibility)
prefactor as `alpha*rhop`, where both factors carried a power of the
Lorentz factor:
$$
\mathtt{alpha} = \frac{q \Delta t}{2m \bar\gamma},
\qquad
\mathtt{rhop} = \frac{q w}{V} \frac{1}{\bar\gamma},
\qquad
\bar\gamma = \tfrac{1}{2}\left(\gamma^{n}+\gamma^{n+1}\right),
$$
so the kernel applied $1/\bar\gamma$ **twice**. The correct response
carries a single power (see
[relativisticPICMCC.pdf](https://github.com/user-attachments/files/31387241/relativisticPICMCC.pdf))
## Impact
- **Theta-implicit EM schemes**: the mass matrices enter only the JFNK
Jacobian approximation and the preconditioner, so the converged
nonlinear solution is unchanged in exact arithmetic. With finite solver
tolerances the iterates differ, which produces small checksum shifts
(relative changes of order $10^{-8}$ in the CI tests, which run at
$\gamma-1 \sim 3\times 10^{-6}$).
- **Semi-implicit Darwin solver**: the mass matrix enters the linear
operator of the magnetoinductive solve directly, so this changes the
converged solution.
---------
Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
EZoni
reviewed
Aug 27, 2026
EZoni
reviewed
Aug 27, 2026
dpgrote
commented
Aug 28, 2026
| @@ -1,24 +1,12 @@ | |||
| { | |||
Member
Author
There was a problem hiding this comment.
Note, the only changes in this file are adding the DTF1_particle_production and DTF2_particle_production benchmark values. Otherwise the changes below are there since the file has been alphabetized. The values are unchanged.
dpgrote
enabled auto-merge (squash)
August 28, 2026 16:10
dpgrote
added a commit
that referenced
this pull request
Aug 28, 2026
The calculation of the global Debye length during binary collisions was added in PR #5763. It was recently realized that the access was being done incorrectly, not taking into account that the MFIter loop was using tiling. With tiling, the `i_cell` index is relative to the tile and cannot be directly used to access MultiFab data. The fix is to use the `global_index` that was add in PR #6954. Note that this error was not detected since the CI test does not use tiling.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
16 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This adds a diagnostic for the nuclear fusion collisions that save the integrated particle production on the grid, turned on by setting the
<collision_name>.save_particle_productionflag. A corresponding flag was added,<collision_name>.create_productsthat allows the creation of particles to be turned off, allowing fusion to be done effectively as a diagnostic and not affect the simulation.Along the way, several extensive changes were needed. A method was needed to create the MultiFab holding the data. This could not be done when the collisions are first setup since the grid is not defined at that point. It could not be done during the collision operation since this done inside an MFIter loop and creating a new MultiFab at that point is problematic. The solution was to add
AllocDatamethods to the collision classes and to the collisionFuncclasses that is called fromMultiParticleContainer::AllocData.When adding the
AllocDatamethod to the collisionFuncclasses, an easy way was to create a new base classCollisionFuncBaseand have the classes all inherit from it. This could be extended, moving a few more things to the base class.Then in
NuclearFusionFunc, when the executor is fetched during the collision, it grabs a pointer to the data from that MultiFab needed in collision. Note that this same method could be used to do a fair amount of clean up in the binary collisions, for example having the density and temperature precalculated and only referenced when needed and not cluttering up the code when they are not needed.Some of the capability implemented here has been pulled out into separate PRs, #7023 and #7025. This PR should only be merged after those two.
Also see PR #6889 which implements similar capabilities for more generic reading and writing MultiFabs.