[aelligp]

What's the purpose of this PR?

• Bug fix
• New feature
• Documentation update
• Other, please explain: test coverage for previously untested kernels, and removal of dead code

Describe it in more detail below:

This PR adds kernel-level tests for code paths that had no coverage. Writing them turned up a number of bugs, which are fixed here. It also removes code that is never loaded or never called.

Bug fixes

Edge-to-edge averaging in 3D stress kernels (src/stokes/StressKernels.jl, used by src/stokes, src/variational_stokes and src/DYREL)

The six av_clamped_{yz,xz,xy}_{x,y,z} helpers move a shear component from one edge family to another. For example, εxy lives at (x vertex, y vertex, z center) and is needed on a yz edge at (x center, y vertex, z vertex). Along each axis the source-to-target step is one of three kinds, and only one of them can meet a missing neighbor:

step along an axis	what the stencil needs	boundary treatment
center → vertex	the two cells meeting at the vertex, (i0, ic)	one cell is missing on a boundary vertex: clamp (zero-gradient) or wrap (periodic)
vertex → center	the two faces of the cell, (i, i+1)	none: every cell, including the last, has both faces
same location	the node itself, i	none: the boundary node ni + 1 is a real node of both arrays

All indices came from clamped_indices, which clamps to the number of cells, ni. That is right for the center → vertex step, but the two other kinds index arrays with ni + 1 entries along that axis:

• vertex → center: the neighbor i + 1 was clamped to nx. At i = nx the stencil averaged A[nx] with itself, so the value came from x_v(nx) instead of x_c(nx), half a cell off.
• same location: jc was clamped to ny. On the last vertex plane (j = ny + 1) the stencil read the plane one step inside the domain.

The helpers now take the unclamped target index I as well as Ic. They use the raw index for the vertex → center and same-location axes, and the clamped/wrapped cell pair only for the center → vertex axis. Callers already guard I by the size of the target edge array, so every raw access is in range. In a periodic direction, vertex planes 1 and ni + 1 are the same physical point, so these axes never have to cross the seam; only the center pair wraps. clamped_indices no longer returns the unused i1, j1, k1.

Principal stresses (src/stokes/PrincipalStresses.jl)

• 2D: the eigenvalue radius used (τxx - τyy)^2 / 2 instead of / 4. The angle used atan(y / x), which loses the quadrant, and now uses atan(y, x).
• 3D: the unshifted Hessenberg-QR solver is replaced by a cyclic Jacobi solver for symmetric 3×3 matrices (eigen_symmetric_3x3). It converges quadratically, including for degenerate spectra, and throws if it hasn't converged after max_sweeps sweeps instead of returning an unconverged result.

Viscosity (src/rheology/Viscosity.jl)

• The 2D single-rheology kernel read the vertex-located Axy at the cell index. It now averages the four surrounding vertices, as the 3D kernel does.
• compute_viscosity_τII without phase ratios now uses the stress tensor (select_tensor) instead of always using the strain rate.
• The 3D single-rheology kernel was commented out; it is restored.
• update_viscosity! passed an extra positional argument to the dispatched function.

3D solver for a single MaterialParams (src/stokes/Stokes3D.jl)
err was never initialized (a bare 1.0 line), and the loop called phase-ratio methods even though no phase ratios exist in this call form. The pressure update now receives G, η_vep is used in the stress update, and the vertex shear stresses come from center2vertex!. Vorticity uses the velocity spacings (grid._di.velocity).

Buoyancy (src/rheology/BuoyancyForces.jl)
compute_ρg! called first(rheology), which fails for a single MaterialParams. gravity_of now handles both a single MaterialParams and a tuple of them.

Stress rotation (src/stress_rotation/stress_rotation_particles.jl, src/ext/*/3D.jl)

• GeoParams.rotate_elastic_stress3D expects the curl of the velocity, which is twice the stored vorticity components.
• The Jaumann update had the wrong sign on τxx, and its kernel was declared (I) instead of (I...).
• The CUDA/AMDGPU 3D rotate_stress_particles! dispatched to the 2D-only rotation-matrix and Jaumann kernels. It now always uses the GeoParams rotation and errors on any other method.

Thermal diffusion (src/thermal_diffusion/)

• Interior fluxes use the distance between adjacent cell centers (_di.center). At a physical boundary, thermal_bcs! reflects the ghost temperature about the boundary, so the ghost-to-center distance is the boundary cell width (_di.vertex). The flux kernels now receive the complete grid._di object, select the appropriate staggered spacing at each face, and use @dx/@dy/@dz for uniform and non-uniform grids.
• update_T! / check_res! divide the flux divergence by the cell width _di.vertex, not _di.center. _di.center has one entry fewer, and on a non-uniform grid this threw a BoundsError.
• Face-located phase ratios (PhaseRatios.Vx/Vy/Vz) are now read at the face, and cell-centered ones at the two neighboring cells (phase_at_face). Any other size throws a DimensionMismatch.
• compute_phase and compute_radioactive_heating for an integer phase now index rheology[phase].

Variational Stokes 3D (src/variational_stokes/)

• In the momentum equations, the shear-stress divergence now uses face differences at each row's own staggered location. The masked 3D _d_xi/_d_yi/_d_zi stencils are corrected.
• Strain rates of inactive cells and edges are set to zero instead of left stale.
• Plastic-strain resets no longer overwrite the edge components owned by the edge blocks.
• The plastic stress correction now includes the missing factor dt.
• The launch range of compute_V! is corrected.
• A bounds guard is added on the edge branches.
• compute_viscosity! is called with air_phase as a keyword.

Removed code

Nothing removed is exported, documented, or used by a miniapp or test.

removed	why
src/thermal_diffusion/DiffusionExplicit.jl (ThermalDiffusion1D/2D/3D)	never included, so never loaded
src/thermal_diffusion/Rheology.jl	never included; duplicates of DiffusionPT_GeoParams.jl
src/stress_rotation/stress_rotation_grid.jl	never included
src/types/constructors/phases.jl (PhaseRatio)	never included; superseded by JustPIC's PhaseRatios
2D _solve! taking K without G (and its di wrapper)	calls compute_P! with 10 arguments and compute_V! with 13; no such methods exist, so it could not run
2D update_τ_o!	refers to the undefined variable τxy
update_stress! (dispatcher and both trait methods)	no callers; the nonlinear method passed its arguments to compute_τ_nonlinear! in an order matching neither of its methods
compute_viscosity! / _compute_viscosity! with air_phase but no phase_ratios	no callers, including the GPU extensions
3D compute_vorticity!(…, _di)	every caller uses the per-component velocity spacings
compute_buoyancy(::MaterialParams, args, phase_ratios)	covered by the generic method

Tests

New files:

• test_stokes_kernels2D.jl, test_stokes_kernels3D.jl: divergence, strain rate, stress (visco-elastic, nonlinear, center/vertex), momentum residual, principal stresses, and all solve! call forms. They include affine-field checks of the 3D edge-to-edge averages at every edge except the one-sided center → boundary-vertex planes, plus boundary-plane checks of the stress update.
• test_variational_stokes3D.jl: masked stencils, strain rate, stress and velocity kernels, and the 3D variational solver (CPU only).
• test_thermal_kernels2D.jl, test_thermal_kernels3D.jl: flux against a hand-written stencil on non-uniform T, K and θ, including physical boundary faces; the energy residual; the pseudo-transient update; Dirichlet pinning; rheology-based flux; and phase ratios of mismatched size.
• test_rheology_kernels2D.jl, test_rheology_kernels3D.jl: viscosity from the strain-rate and stress invariants (relaxation, cutoff), harmonic phase averaging without the air phase, vertex arguments, and buoyancy.
• test_stress_rotation2D.jl, test_stress_rotation3D.jl: vorticity of a rigid rotation, and the per-particle rotation kernels.
• test_bc_kernels2D.jl, test_bc_kernels3D.jl: no-slip velocity boundaries on every face, and pure/simple-shear background fields.
• test_weno5_kernels.jl: WENO5 weights and convergence order.

Extended: test_diffusion2D/3D.jl (non-uniform grid), test_types.jl, test_rheology.jl, test_rockratio.jl, test_mask.jl, test_phase_ratios3D.jl.

The full CPU Pkg.test() suite passes, including the MPI tests. CUDA and AMDGPU were not run locally.

Checklist

• The PR title is descriptive and starts with the appropriate tag: [BUGFIX], [ADDITION], [DOC], etc.
• New tests (either assessing the correct behavior of new internal functions or the correctness of a miniapps or benchmarks) were added, or old tests were updated
• Affected miniapps have also been updated
• The new feature was added in a way that does not break public API
• New documentation related to the new feature was added
• The new code follows the contributor guidelines, in particular the Runic Style

🤖 Generated with Claude Code

[github-actions]

Your PR no longer requires formatting changes. Thank you for your contribution!

[aelligp]

@albert-de-montserrat gpu still needs fixing

[codecov]

Codecov Report

❌ Patch coverage is 94.63087% with 8 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/rheology/Viscosity.jl	82.60%	4 Missing ⚠️
src/thermal_diffusion/DiffusionPT_GeoParams.jl	81.81%	2 Missing ⚠️
src/variational_stokes/StressKernels.jl	50.00%	2 Missing ⚠️

📢 Thoughts on this report? Let us know!

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.