diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md
index c9b9860ab9..279e946bda 100644
--- a/docs/src/systems/fluid.md
+++ b/docs/src/systems/fluid.md
@@ -204,6 +204,24 @@ Pages = [joinpath("schemes", "fluid", "viscosity.jl")]
 
 ## [Corrections](@id corrections)
 
+Gradient corrections generally make the two kernel gradients of a particle pair asymmetric.
+The corresponding eight-argument pressure formulations are therefore selected from the
+configured correction method, even when a particular pair happens to satisfy
+``\nabla W_b = -\nabla W_a``. In that symmetric case, the asymmetric formulation reduces to
+the standard symmetric formulation.
+
+The antisymmetric combination of corrected gradients preserves pairwise linear momentum.
+Since corrected gradients are generally not parallel to the particle separation, the resulting
+force is not necessarily central and does not in general preserve angular momentum. The usual
+linear- and angular-momentum guarantee applies to symmetric radial kernel gradients.
+
+When EDAC average-pressure reduction is enabled, two interacting EDAC particles use the
+arithmetic mean of their local pressure offsets. The shared pair offset ensures that both
+directed evaluations use identical reduced pair pressures, which is required for the
+antisymmetric corrected-gradient formulation to preserve linear momentum. Interactions between
+schemes using different pressure formulations do not gain a conservation guarantee from this
+construction.
+
 ```@autodocs
 Modules = [TrixiParticles]
 Pages = [joinpath("general", "corrections.jl")]
diff --git a/src/general/abstract_system.jl b/src/general/abstract_system.jl
index 79b52a81df..954f836107 100644
--- a/src/general/abstract_system.jl
+++ b/src/general/abstract_system.jl
@@ -163,6 +163,24 @@ end
                                         system_correction(system), system, particle)
 end
 
+# Hydrodynamic corrections of structure systems are stored in their boundary model and are
+# independent of corrections used by the structural scheme itself.
+@inline hydrodynamic_correction(system) = system_correction(system)
+
+@inline function hydrodynamic_smoothing_kernel_grad(system, pos_diff, distance, particle)
+    h = smoothing_length(system, particle)
+    correction = hydrodynamic_correction(system)
+    compact_support_ = compact_support(system_smoothing_kernel(system), h)
+
+    if distance >= compact_support_ ||
+       (skip_zero_distance(correction) && distance^2 < eps(h^2))
+        return zero(pos_diff)
+    end
+
+    return corrected_kernel_grad_unsafe(system_smoothing_kernel(system), pos_diff,
+                                        distance, h, correction, system, particle)
+end
+
 # System updates do nothing by default, but can be dispatched if needed
 function update_positions!(system, v, u, v_ode, u_ode, semi, t)
     return system
@@ -172,6 +190,10 @@ function update_quantities!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
+function update_density_correction!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
 function update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
@@ -180,6 +202,14 @@ function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
+function update_gradient_correction!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
+function update_surface_quantities!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
 function update_final!(system, v, u, v_ode, u_ode, semi, t; kwargs...)
     return system
 end
diff --git a/src/general/corrections.jl b/src/general/corrections.jl
index d97eeaf341..19c3a11853 100644
--- a/src/general/corrections.jl
+++ b/src/general/corrections.jl
@@ -52,7 +52,8 @@ end
     ShepardKernelCorrection()
 
 Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation
-to obtain a 0-th order accurate result, which was first proposed by [Li et al. (1996)](@cite Li1996).
+to obtain a zeroth-order consistent result (exact reproduction of constants), which was first
+proposed by [Li et al. (1996)](@cite Li1996).
 
 The kernel correction coefficient is determined by
 ```math
@@ -61,7 +62,10 @@ c(x) = \sum_{b=1} V_b W_b(x),
 where ``V_b = m_b / \rho_b`` is the volume of particle ``b``.
 
 This correction is applied with [`SummationDensity`](@ref) to correct the density and leads
-to an improvement, especially at free surfaces.
+to an improvement, especially at free surfaces. With summation density, the current one-pass
+implementation uses the provisional density in ``V_b`` and therefore reduces the free-surface
+error without guaranteeing convergence. [`DensityReinitializationCallback`](@ref) instead uses
+the independently evolved continuity density and realizes the consistent Shepard operator.
 
 !!! note
     - It is also referred to as "0th order correction".
@@ -73,7 +77,8 @@ struct ShepardKernelCorrection end
     KernelCorrection()
 
 Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation
-to obtain a 0-th order accurate result, which was first proposed by Li et al.
+to obtain a zeroth-order consistent kernel gradient (an exact zero gradient for constants),
+which was first proposed by Li et al.
 This can be further extended to obtain a kernel corrected gradient as shown by [Basa et al. (2008)](@cite Basa2008).
 
 The kernel correction coefficient is determined by
@@ -100,7 +105,8 @@ struct KernelCorrection end
     MixedKernelGradientCorrection()
 
 Combines [`GradientCorrection`](@ref) and [`KernelCorrection`](@ref),
-which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1999)).
+which results in a first-order consistent kernel gradient reproducing both constant and affine
+fields exactly (see [Bonet, 1999](@cite Bonet1999)).
 
 # Notes:
 - Stability issues, especially when particles separate into small clusters.
@@ -108,6 +114,16 @@ which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1
 """
 struct MixedKernelGradientCorrection end
 
+correction_density(::Any) = nothing
+correction_density(correction::ShepardKernelCorrection) = correction
+
+correction_gradient(::Nothing) = nothing
+correction_gradient(::ShepardKernelCorrection) = nothing
+correction_gradient(::AkinciFreeSurfaceCorrection) = nothing
+correction_gradient(correction) = correction
+
+correction_force(correction) = correction
+
 function kernel_correction_coefficient(system::AbstractFluidSystem, particle)
     return system.cache.kernel_correction_coefficient[particle]
 end
@@ -166,9 +182,24 @@ function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi,
         end
     end
 
+    sanitize_kernel_correction_coefficient!(kernel_correction_coefficient, system, semi)
+
     return kernel_correction_coefficient
 end
 
+function sanitize_kernel_correction_coefficient!(coefficient, system, semi)
+    minimum_coefficient = sqrt(eps(eltype(coefficient)))
+
+    @threaded semi for particle in eachparticle(system)
+        value = coefficient[particle]
+        if !isfinite(value) || value <= minimum_coefficient
+            coefficient[particle] = one(value)
+        end
+    end
+
+    return coefficient
+end
+
 function dw_gamma(system::AbstractFluidSystem, particle)
     return extract_svector(system.cache.dw_gamma, system, particle)
 end
@@ -255,9 +286,22 @@ function compute_correction_values!(system,
         end
     end
 
-    for particle in eachparticle(system), i in axes(dw_gamma, 1)
-        dw_gamma[i, particle] /= kernel_correction_coefficient[particle]
+    minimum_coefficient = sqrt(eps(eltype(kernel_correction_coefficient)))
+    @threaded semi for particle in eachparticle(system)
+        coefficient = kernel_correction_coefficient[particle]
+        if !isfinite(coefficient) || coefficient <= minimum_coefficient
+            kernel_correction_coefficient[particle] = one(coefficient)
+            for i in axes(dw_gamma, 1)
+                dw_gamma[i, particle] = zero(eltype(dw_gamma))
+            end
+        else
+            for i in axes(dw_gamma, 1)
+                dw_gamma[i, particle] /= coefficient
+            end
+        end
     end
+
+    return kernel_correction_coefficient
 end
 
 @doc raw"""
@@ -284,6 +328,9 @@ The gradient correction, as commonly proposed, involves multiplying this gradien
 
 The correction matrix  $\bm{L}_a$ is computed based on the provided particle configuration,
 aiming to make the corrected gradient more accurate, especially near domain boundaries.
+It gives a first-order consistent gradient by differentiating every affine field exactly.
+For smooth fields, the local truncation error is generally ``O(h)`` on asymmetric supports and
+``O(h^2)`` on symmetric interior supports.
 
 To satisfy
 ```math
@@ -313,6 +360,8 @@ This calculates the following,
 \tilde\nabla A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i
 ```
 with ``0 \leq \lambda \leq 1`` being the blending factor.
+For a fixed ``\lambda < 1``, the uncorrected first-moment error remains and no asymptotic order
+improvement is guaranteed.
 
 # Arguments
 - `blending_factor`: Blending factor between corrected and regular SPH gradient.
@@ -321,6 +370,10 @@ struct BlendedGradientCorrection{ELTYPE <: Real}
     blending_factor::ELTYPE
 
     function BlendedGradientCorrection(blending_factor)
+        if !(zero(blending_factor) <= blending_factor <= one(blending_factor))
+            throw(ArgumentError("`blending_factor` must be between 0 and 1"))
+        end
+
         return new{eltype(blending_factor)}(blending_factor)
     end
 end
@@ -376,8 +429,10 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
                                    semi) do particle, neighbor, pos_diff, distance
                 function kernel_grad_local(correction, smoothing_kernel, pos_diff, distance,
                                            smoothing_length_, system, particle)
-                    return smoothing_kernel_grad_unsafe(system, pos_diff, distance,
-                                                        particle)
+                    # Do not dispatch through `system`: the correction matrix being used
+                    # by that path is the matrix currently being assembled here.
+                    return kernel_grad_unsafe(smoothing_kernel, pos_diff, distance,
+                                              smoothing_length_)
                 end
 
                 # Compute gradient of corrected kernel
@@ -426,8 +481,9 @@ function correction_matrix_inversion_step!(corr_matrix, system, semi)
     @threaded semi for particle in eachparticle(system)
         L = extract_smatrix(corr_matrix, system, particle)
 
-        # The matrix `L` only becomes singular when the particle and all neighbors
-        # are collinear (in 2D) or lie all in the same plane (in 3D).
+        # The matrix `L` becomes singular when the particle and all neighbors are collinear
+        # (in 2D) or lie all in the same plane (in 3D). Nearly singular matrices are also
+        # rejected below to avoid amplifying particle disorder.
         # This happens only when two (in 2D) or three (in 3D) particles are isolated,
         # or in cases where there is only one layer of fluid particles on a wall.
         # In these edge cases, we just disable the correction and set the corrected
@@ -441,7 +497,12 @@ function correction_matrix_inversion_step!(corr_matrix, system, semi)
         # so `L` is singular if and only if the position vectors X_ab don't span the
         # full space, i.e., particle a and all neighbors lie on the same line (in 2D)
         # or plane (in 3D).
-        if abs(det(L)) < 1.0f-9
+        scale = maximum(abs, L)
+        relative_determinant = abs(det(L)) / scale^ndims(system)
+        minimum_relative_determinant = sqrt(eps(eltype(L)))
+
+        if !isfinite(relative_determinant) ||
+           relative_determinant < minimum_relative_determinant
             L_inv = I
         else
             L_inv = inv(L)
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index a9ed769dc6..bb5fc5363d 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -633,7 +633,13 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
 
     update_implicit_sph!(semi, v_ode, u_ode, t)
 
-    # Perform correction and pressure calculation
+    # Correction moments can use densities from every interacting system, so density
+    # correction has to be a global phase.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_density_correction!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
+    # Fluid pressure must be available before boundary pressure interpolation.
     foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
         update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     end
@@ -644,6 +650,18 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t)
         update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     end
 
+    # Boundary interpolation can update boundary density. Assemble all gradient corrections
+    # only after every interacting system exposes its final density.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_gradient_correction!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
+    # Surface quantities can depend on corrected gradients and must be complete for every
+    # system before curvature and stress are computed in `update_final!`.
+    foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
+        update_surface_quantities!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
     # Final update step for all remaining systems
     foreach_system_wrapped(semi, v_ode, u_ode) do system, v, u
         update_final!(system, v, u, v_ode, u_ode, semi, t)
diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl
index fcde467f2e..88a3a2f3ed 100644
--- a/src/schemes/boundary/wall_boundary/dummy_particles.jl
+++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl
@@ -237,19 +237,19 @@ function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particl
 end
 
 function create_cache_model(::KernelCorrection, density, NDIMS, n_particles)
-    dw_gamma = Array{Float64}(undef, NDIMS, n_particles)
+    dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles)
     return (; kernel_correction_coefficient=similar(density), dw_gamma)
 end
 
 function create_cache_model(::Union{GradientCorrection, BlendedGradientCorrection}, density,
                             NDIMS, n_particles)
-    correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles)
+    correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles)
     return (; correction_matrix)
 end
 
 function create_cache_model(::MixedKernelGradientCorrection, density, NDIMS, n_particles)
-    dw_gamma = Array{Float64}(undef, NDIMS, n_particles)
-    correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles)
+    dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles)
+    correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles)
     return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix)
 end
 
@@ -393,21 +393,62 @@ end
 
 @inline function update_pressure!(boundary_model::BoundaryModelDummyParticles,
                                   system, v, u, v_ode, u_ode, semi)
-    (; correction, density_calculator) = boundary_model
+    (; density_calculator) = boundary_model
 
     compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
 
-    # These are only computed when using corrections
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
-    compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode,
-                                        semi)
-    # `kernel_correct_density!` only performed for `SummationDensity`
-    kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction,
+    return boundary_model
+end
+
+@inline function update_density_correction!(boundary_model::BoundaryModelDummyParticles,
+                                            system, v, u, v_ode, u_ode, semi)
+    (; correction, density_calculator) = boundary_model
+    density_correction = correction_density(correction)
+
+    compute_boundary_correction_values!(boundary_model, system, density_correction, u,
+                                        v_ode, u_ode, semi)
+    kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi,
+                            density_correction,
                             density_calculator)
 
     return boundary_model
 end
 
+@inline function update_gradient_correction!(boundary_model::BoundaryModelDummyParticles,
+                                             system, v, u, v_ode, u_ode, semi)
+    gradient_correction = correction_gradient(boundary_model.correction)
+
+    compute_boundary_correction_values!(boundary_model, system, gradient_correction, u,
+                                        v_ode, u_ode, semi)
+    compute_gradient_correction_matrix!(gradient_correction, boundary_model, system, u,
+                                        v_ode, u_ode, semi)
+
+    return boundary_model
+end
+
+@inline function compute_boundary_correction_values!(boundary_model, system, correction, u,
+                                                     v_ode, u_ode, semi)
+    return boundary_model
+end
+
+function compute_boundary_correction_values!(boundary_model, system,
+                                             ::ShepardKernelCorrection, u,
+                                             v_ode, u_ode, semi)
+    return compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode,
+                                  semi,
+                                  boundary_model.cache.kernel_correction_coefficient)
+end
+
+function compute_boundary_correction_values!(boundary_model, system,
+                                             correction::Union{KernelCorrection,
+                                                               MixedKernelGradientCorrection},
+                                             u, v_ode, u_ode, semi)
+    return compute_correction_values!(system, correction, current_coordinates(u, system),
+                                      v_ode, u_ode, semi,
+                                      boundary_model.cache.kernel_correction_coefficient,
+                                      boundary_model.cache.dw_gamma)
+end
+
 function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi,
                                  correction, density_calculator)
     return boundary_model
@@ -428,13 +469,13 @@ function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
                                                          MixedKernelGradientCorrection},
                                              boundary_model,
                                              system, u, v_ode, u_ode, semi)
-    (; cache, correction, smoothing_kernel) = boundary_model
+    (; cache, smoothing_kernel) = boundary_model
     (; correction_matrix) = cache
 
     system_coords = current_coordinates(u, system)
 
     compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
-                                        v_ode, u_ode, semi, correction, smoothing_kernel)
+                                        v_ode, u_ode, semi, corr, smoothing_kernel)
 end
 
 function compute_density!(boundary_model, ::SummationDensity, system, v, u, v_ode, u_ode,
diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl
index f9864ecc05..1e8c7a4e03 100644
--- a/src/schemes/boundary/wall_boundary/system.jl
+++ b/src/schemes/boundary/wall_boundary/system.jl
@@ -218,6 +218,13 @@ function update_quantities!(system::WallBoundarySystem, v, u, v_ode, u_ode, semi
     return system
 end
 
+function update_density_correction!(system::WallBoundarySystem{<:BoundaryModelDummyParticles},
+                                    v, u, v_ode, u_ode, semi, t)
+    update_density_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 # This update depends on the computed quantities of the fluid system and therefore
 # has to be in `update_boundary_interpolation!` after `update_quantities!`.
 function update_boundary_interpolation!(system::WallBoundarySystem, v, u, v_ode, u_ode,
@@ -231,6 +238,13 @@ function update_boundary_interpolation!(system::WallBoundarySystem, v, u, v_ode,
     return system
 end
 
+function update_gradient_correction!(system::WallBoundarySystem{<:BoundaryModelDummyParticles},
+                                     v, u, v_ode, u_ode, semi, t)
+    update_gradient_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 function write_u0!(u0, ::WallBoundarySystem)
     return u0
 end
@@ -329,7 +343,7 @@ function system_smoothing_kernel(system::WallBoundarySystem{<:BoundaryModelDummy
 end
 
 function system_correction(system::WallBoundarySystem{<:BoundaryModelDummyParticles})
-    return system.boundary_model.correction
+    return correction_gradient(system.boundary_model.correction)
 end
 
 @inline function density_calculator(system::WallBoundarySystem)
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index 2ea1e4dc4d..5be0122510 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -4,6 +4,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                    particle_system::EntropicallyDampedSPHSystem,
                    neighbor_system, semi)
     (; sound_speed, density_calculator, correction, nu_edac) = particle_system
+    gradient_correction = correction_gradient(correction)
 
     system_coords = current_coordinates(u_particle_system, particle_system)
     neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
@@ -54,7 +55,8 @@ function interact!(dv, v_particle_system, u_particle_system,
         # It results in significant improvement for EDAC, especially with TVF,
         # but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2.
         # Note that the return value is zero when not using average pressure reduction.
-        p_avg = @inbounds average_pressure(particle_system, particle)
+        p_avg = @inbounds pair_pressure_offset(particle_system, neighbor_system, particle,
+                                               neighbor)
 
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
@@ -63,7 +65,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             particle, neighbor,
                                             m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a,
                                             rho_b, pos_diff, distance, grad_kernel,
-                                            correction)
+                                            gradient_correction)
 
         dv_particle = Ref(dv_pressure)
         @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
@@ -77,7 +79,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                particle_system, neighbor_system,
                                v_particle_system, v_neighbor_system,
                                particle, neighbor, m_a, m_b, rho_a, rho_b, v_a, v_b,
-                               pos_diff, distance, grad_kernel, correction)
+                               pos_diff, distance, grad_kernel, gradient_correction)
 
         @inbounds surface_tension_force!(dv_particle, surface_tension_a,
                                          surface_tension_b,
diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl
index 6e03cf7837..64ff0c741f 100644
--- a/src/schemes/fluid/entropically_damped_sph/system.jl
+++ b/src/schemes/fluid/entropically_damped_sph/system.jl
@@ -8,7 +8,7 @@
                                 acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing,
                                 surface_normal_method=nothing, buffer_size=nothing,
                                 reference_particle_spacing=0.0, color_value=1,
-                                source_terms=nothing)
+                                correction=nothing, source_terms=nothing)
 
 System for particles of a fluid.
 As opposed to the [weakly compressible SPH scheme](@ref wcsph), which uses an equation of state,
@@ -35,7 +35,10 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                                 formulation](@ref transport_velocity_formulation) to use
                                 with this system. Default is no shifting.
 - `average_pressure_reduction`: Whether to subtract the average pressure of neighboring particles
-                                from the local pressure (default: `true` when using shifting, `false` otherwise).
+                                 from the local pressure (default: `true` when using shifting, `false` otherwise).
+                                 Interacting EDAC particles use the arithmetic mean of their local
+                                 pressure offsets so both directed pair evaluations use the same
+                                 reduced pressures and preserve linear momentum.
 - `buffer_size`:                Number of buffer particles.
                                 This is needed when simulating with [`OpenBoundarySystem`](@ref).
 - `correction`:                 Correction method used for this system. (default: no correction, see [Corrections](@ref corrections))
@@ -110,6 +113,9 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth
     mass = copy(initial_condition.mass)
     n_particles = length(initial_condition.mass)
 
+    density_correction_ = correction_density(correction)
+    gradient_correction_ = correction_gradient(correction)
+
     if ndims(smoothing_kernel) != NDIMS
         throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
     end
@@ -127,7 +133,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth
         throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model"))
     end
 
-    if correction isa ShepardKernelCorrection &&
+    if density_correction_ isa ShepardKernelCorrection &&
        density_calculator isa ContinuityDensity
         throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`"))
     end
@@ -135,7 +141,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth
     pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
                                                                      density_calculator,
                                                                      NDIMS, ELTYPE,
-                                                                     correction)
+                                                                     gradient_correction_)
 
     avg_pressure_reduction = Val(average_pressure_reduction)
 
@@ -251,7 +257,9 @@ end
 
 @inline buffer(system::EntropicallyDampedSPHSystem) = system.buffer
 
-system_correction(system::EntropicallyDampedSPHSystem) = system.correction
+function system_correction(system::EntropicallyDampedSPHSystem)
+    correction_gradient(system.correction)
+end
 
 @inline function current_velocity(v, system::EntropicallyDampedSPHSystem)
     return view(v, 1:ndims(system), :)
@@ -273,6 +281,22 @@ end
 
 @inline average_pressure(system, ::Val{false}, particle) = zero(eltype(system))
 
+@propagate_inbounds function interaction_pressure_offset(system::EntropicallyDampedSPHSystem,
+                                                         particle)
+    return average_pressure(system, particle)
+end
+
+# Both directed evaluations of an EDAC pair must use the same pressure offset to preserve
+# pairwise linear momentum. This also handles pairs where only one system enables reduction.
+@propagate_inbounds function pair_pressure_offset(system::EntropicallyDampedSPHSystem,
+                                                  neighbor_system::EntropicallyDampedSPHSystem,
+                                                  particle, neighbor)
+    pressure_offset_a = average_pressure(system, particle)
+    pressure_offset_b = average_pressure(neighbor_system, neighbor)
+
+    return (pressure_offset_a + pressure_offset_b) / 2
+end
+
 @inline function current_density(v, system::EntropicallyDampedSPHSystem)
     return current_density(v, system.density_calculator, system)
 end
@@ -298,12 +322,69 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u,
     compute_density!(system, u, u_ode, semi, system.density_calculator)
 end
 
-function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t)
+function update_density_correction!(system::EntropicallyDampedSPHSystem, v, u, v_ode,
+                                    u_ode, semi, t)
+    (; correction, density_calculator) = system
+    density_correction = correction_density(correction)
+
+    compute_correction_values!(system, density_correction, u, v_ode, u_ode, semi)
+    kernel_correct_density!(system, v, u, v_ode, u_ode, semi, density_correction,
+                            density_calculator)
+
+    return system
+end
+
+function update_gradient_correction!(system::EntropicallyDampedSPHSystem, v, u, v_ode,
+                                     u_ode, semi, t)
+    gradient_correction = correction_gradient(system.correction)
+
+    compute_correction_values!(system, gradient_correction, u, v_ode, u_ode, semi)
+    compute_gradient_correction_matrix!(gradient_correction, system, u, v_ode, u_ode, semi)
+
+    return system
+end
+
+function update_surface_quantities!(system::EntropicallyDampedSPHSystem, v, u, v_ode,
+                                    u_ode, semi, t)
     compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi,
                             t)
     compute_surface_delta_function!(system, system.surface_tension, semi)
 end
 
+function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode,
+                                 semi, correction, density_calculator)
+    return system
+end
+
+function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode,
+                                 semi, ::ShepardKernelCorrection, ::SummationDensity)
+    system.cache.density ./= system.cache.kernel_correction_coefficient
+end
+
+function compute_gradient_correction_matrix!(correction,
+                                             system::EntropicallyDampedSPHSystem, u,
+                                             v_ode, u_ode, semi)
+    return system
+end
+
+function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
+                                                         BlendedGradientCorrection,
+                                                         MixedKernelGradientCorrection},
+                                             system::EntropicallyDampedSPHSystem, u,
+                                             v_ode, u_ode, semi)
+    (; cache, smoothing_kernel) = system
+    (; correction_matrix) = cache
+
+    system_coords = current_coordinates(u, system)
+
+    compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
+                                        v_ode, u_ode, semi, corr, smoothing_kernel)
+end
+
+@inline function correction_matrix(system::EntropicallyDampedSPHSystem, particle)
+    extract_smatrix(system.cache.correction_matrix, system, particle)
+end
+
 function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t;
                        kwargs...)
     (; surface_tension) = system
@@ -387,6 +468,8 @@ function restart_with!(system::EntropicallyDampedSPHSystem, v, u)
     for particle in each_integrated_particle(system)
         system.initial_condition.coordinates[:, particle] .= u[:, particle]
         system.initial_condition.velocity[:, particle] .= v[1:ndims(system), particle]
-        system.initial_condition.pressure[particle] = v[end, particle]
+        system.initial_condition.pressure[particle] = v[ndims(system) + 1, particle]
     end
+
+    return restart_with!(system, system.density_calculator, v, u)
 end
diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl
index b6114c0bc1..e20448e844 100644
--- a/src/schemes/fluid/pressure_acceleration.jl
+++ b/src/schemes/fluid/pressure_acceleration.jl
@@ -76,8 +76,10 @@ end
 # different inter-particle averages or to assume different inter-particle distributions.
 # Ramachandran (2019) and Adami (2012) use this formulation for the pressure acceleration.
 #
-# However, the tests show that the formulation is only linear and angular momentum conserving
-# but not energy conserving.
+# With a symmetric radial kernel gradient, the formulation conserves linear and angular
+# momentum, but not energy. With asymmetric corrected gradients, the formulation below still
+# conserves linear momentum, while the generally non-central pair force does not conserve
+# angular momentum.
 #
 # Note that the authors also used this formulation for an ISPH method in (https://doi.org/10.1016/j.jcp.2007.07.013)
 #
@@ -93,6 +95,21 @@ end
     return -volume_term * pressure_tilde * W_a
 end
 
+# Linear-momentum-conserving extension for correction methods with asymmetric kernel gradients.
+# The asymmetric overload is selected based on the configured correction, even when a particular
+# particle pair happens to have symmetric gradients. It reduces to the symmetric formulation
+# above when `W_b == -W_a`.
+@inline function inter_particle_averaged_pressure(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a,
+                                                  W_b)
+    volume_a = m_a / rho_a
+    volume_b = m_b / rho_b
+    volume_term = (volume_a^2 + volume_b^2) / m_a
+    pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b)
+
+    half = oftype(volume_term, 0.5)
+    return -half * volume_term * pressure_tilde * (W_a - W_b)
+end
+
 function choose_pressure_acceleration_formulation(pressure_acceleration,
                                                   density_calculator, NDIMS, ELTYPE,
                                                   correction)
@@ -143,6 +160,15 @@ end
 
 @inline pressure_acceleration_formulation(system) = system.pressure_acceleration_formulation
 
+# Use the local pressure offset by default. Specialized pair methods can combine offsets from
+# both systems when both directed interaction evaluations support the same reduction.
+@inline interaction_pressure_offset(system, particle) = zero(eltype(system))
+
+@propagate_inbounds function pair_pressure_offset(system, neighbor_system, particle,
+                                                  neighbor)
+    return interaction_pressure_offset(system, particle)
+end
+
 # Formulation using symmetric gradient formulation for corrections not depending on local neighborhood.
 @inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor,
                                        m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
@@ -161,7 +187,7 @@ end
                                                          GradientCorrection,
                                                          BlendedGradientCorrection,
                                                          MixedKernelGradientCorrection})
-    W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
+    W_b = hydrodynamic_smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor)
 
     # With correction, the kernel gradient is not necessarily symmetric, so call the
     # asymmetric version of the pressure acceleration formulation.
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index 8ceb29d27f..e63fabbace 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -8,6 +8,8 @@ function interact!(dv, v_particle_system, u_particle_system,
                    eachparticle=each_integrated_particle(particle_system),
                    kwargs...)
     (; density_calculator, correction) = particle_system
+    gradient_correction = correction_gradient(correction)
+    force_correction = correction_force(correction)
 
     sound_speed = system_sound_speed(particle_system)
 
@@ -72,7 +74,7 @@ function interact!(dv, v_particle_system, u_particle_system,
             # Determine correction factors.
             # This can usually be ignored, as these are all 1 when no correction is used.
             (viscosity_correction, pressure_correction,
-             surface_tension_correction) = free_surface_correction(correction,
+             surface_tension_correction) = free_surface_correction(force_correction,
                                                                    particle_system,
                                                                    rho_a, rho_b)
 
@@ -81,7 +83,7 @@ function interact!(dv, v_particle_system, u_particle_system,
             dv_pressure = pressure_acceleration(particle_system, neighbor_system,
                                                 particle, neighbor,
                                                 m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
-                                                distance, grad_kernel, correction)
+                                                distance, grad_kernel, gradient_correction)
             dv_particle[] += dv_pressure * pressure_correction
 
             # Propagate `@inbounds` to the viscosity function, which accesses particle data
@@ -96,7 +98,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                    particle_system, neighbor_system,
                                    v_particle_system, v_neighbor_system,
                                    particle, neighbor, m_a, m_b, rho_a, rho_b, v_a, v_b,
-                                   pos_diff, distance, grad_kernel, correction)
+                                   pos_diff, distance, grad_kernel, gradient_correction)
 
             @inbounds surface_tension_force!(dv_particle,
                                              surface_tension_a, surface_tension_b,
diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl
index eea0607d7d..7d5fa2c2e8 100644
--- a/src/schemes/fluid/weakly_compressible_sph/system.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/system.jl
@@ -112,6 +112,9 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
     ELTYPE = eltype(initial_condition)
     n_particles = nparticles(initial_condition)
 
+    density_correction_ = correction_density(correction)
+    gradient_correction_ = correction_gradient(correction)
+
     mass = copy(initial_condition.mass)
     pressure = similar(initial_condition.pressure)
 
@@ -125,7 +128,7 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
         throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
     end
 
-    if correction isa ShepardKernelCorrection &&
+    if density_correction_ isa ShepardKernelCorrection &&
        density_calculator isa ContinuityDensity
         throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`"))
     end
@@ -141,7 +144,7 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
     pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
                                                                      density_calculator,
                                                                      NDIMS, ELTYPE,
-                                                                     correction)
+                                                                     gradient_correction_)
 
     cache = (; create_cache_density(initial_condition, density_calculator)...,
              create_cache_correction(correction, initial_condition.density, NDIMS,
@@ -243,7 +246,9 @@ end
 
 @inline buffer(system::WeaklyCompressibleSPHSystem) = system.buffer
 
-system_correction(system::WeaklyCompressibleSPHSystem) = system.correction
+function system_correction(system::WeaklyCompressibleSPHSystem)
+    correction_gradient(system.correction)
+end
 
 @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem)
     return current_velocity(v, system.density_calculator, system)
@@ -320,17 +325,37 @@ end
     return system
 end
 
-function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
-    (; density_calculator, correction, surface_normal_method, surface_tension) = system
+function update_density_correction!(system::WeaklyCompressibleSPHSystem, v, u, v_ode,
+                                    u_ode, semi, t)
+    (; density_calculator, correction) = system
+    density_correction = correction_density(correction)
+
+    compute_correction_values!(system, density_correction, u, v_ode, u_ode, semi)
+    kernel_correct_density!(system, v, u, v_ode, u_ode, semi, density_correction,
+                            density_calculator)
+
+    return system
+end
 
+function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t)
     compute_pressure!(system, v, semi)
 
-    # These are only computed when using corrections
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
-    compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi)
-    # `kernel_correct_density!` only performed for `SummationDensity`
-    kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction,
-                            density_calculator)
+    return system
+end
+
+function update_gradient_correction!(system::WeaklyCompressibleSPHSystem, v, u, v_ode,
+                                     u_ode, semi, t)
+    gradient_correction = correction_gradient(system.correction)
+
+    compute_correction_values!(system, gradient_correction, u, v_ode, u_ode, semi)
+    compute_gradient_correction_matrix!(gradient_correction, system, u, v_ode, u_ode, semi)
+
+    return system
+end
+
+function update_surface_quantities!(system::WeaklyCompressibleSPHSystem, v, u, v_ode,
+                                    u_ode, semi, t)
+    (; surface_normal_method, surface_tension) = system
 
     # These are only computed when using surface tension
     compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t)
@@ -369,13 +394,13 @@ function compute_gradient_correction_matrix!(corr::Union{GradientCorrection,
                                                          MixedKernelGradientCorrection},
                                              system::WeaklyCompressibleSPHSystem, u,
                                              v_ode, u_ode, semi)
-    (; cache, correction, smoothing_kernel) = system
+    (; cache, smoothing_kernel) = system
     (; correction_matrix) = cache
 
     system_coords = current_coordinates(u, system)
 
     compute_gradient_correction_matrix!(correction_matrix, system, system_coords,
-                                        v_ode, u_ode, semi, correction, smoothing_kernel)
+                                        v_ode, u_ode, semi, corr, smoothing_kernel)
 end
 
 function reinit_density!(vu_ode, semi)
@@ -397,14 +422,16 @@ end
 
 function reinit_density!(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity, v, u,
                          v_ode, u_ode, semi)
+    # Use the independently evolved density to determine particle volumes before replacing it
+    # with the reinitialized summation density.
+    kernel_correction_coefficient = similar(v, size(v, 2))
+    compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi,
+                           kernel_correction_coefficient)
+
     # Compute density with `SummationDensity` and store the result in `v`,
     # overwriting the previous integrated density.
     summation_density!(system, semi, u, u_ode, v[end, :])
 
-    # Apply `ShepardKernelCorrection`
-    kernel_correction_coefficient = zeros(size(v[end, :]))
-    compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi,
-                           kernel_correction_coefficient)
     @threaded semi for particle in eachparticle(system)
         v[end, particle] /= kernel_correction_coefficient[particle]
     end
diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl
index dccacc2668..9dea18c51d 100644
--- a/src/schemes/structure/rigid_body/system.jl
+++ b/src/schemes/structure/rigid_body/system.jl
@@ -267,7 +267,25 @@ end
 end
 
 @inline function system_correction(system::RigidBodySystem{<:BoundaryModelDummyParticles})
-    return system.boundary_model.correction
+    return correction_gradient(system.boundary_model.correction)
+end
+
+@inline function hydrodynamic_correction(system::RigidBodySystem{<:BoundaryModelDummyParticles})
+    return correction_gradient(system.boundary_model.correction)
+end
+
+@inline function kernel_correction_coefficient(system::RigidBodySystem{<:BoundaryModelDummyParticles},
+                                               particle)
+    return system.boundary_model.cache.kernel_correction_coefficient[particle]
+end
+
+@inline function dw_gamma(system::RigidBodySystem{<:BoundaryModelDummyParticles}, particle)
+    return extract_svector(system.boundary_model.cache.dw_gamma, system, particle)
+end
+
+@inline function correction_matrix(system::RigidBodySystem{<:BoundaryModelDummyParticles},
+                                   particle)
+    return extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle)
 end
 
 function initialize!(system::RigidBodySystem, semi)
@@ -360,6 +378,13 @@ function restart_with!(system::RigidBodySystem, v, u)
     return system
 end
 
+function update_density_correction!(system::RigidBodySystem{<:BoundaryModelDummyParticles},
+                                    v, u, v_ode, u_ode, semi, t)
+    update_density_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 function update_boundary_interpolation!(system::RigidBodySystem, v, u, v_ode, u_ode,
                                         semi, t)
     return update_boundary_interpolation!(system.boundary_model, system, v, u, v_ode,
@@ -377,6 +402,13 @@ function update_boundary_interpolation!(boundary_model, system::RigidBodySystem,
     return system
 end
 
+function update_gradient_correction!(system::RigidBodySystem{<:BoundaryModelDummyParticles},
+                                     v, u, v_ode, u_ode, semi, t)
+    update_gradient_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 function reset_interaction_caches!(system::RigidBodySystem)
     set_zero!(system.force_per_particle)
     system.cache.contact_count[] = 0
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index dd7e8d38ea..f13e5e520c 100644
--- a/src/schemes/structure/structure.jl
+++ b/src/schemes/structure/structure.jl
@@ -43,10 +43,7 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
         # Note that `return` only exits the closure, i.e., skips the current neighbor.
         skip_zero_distance(neighbor_system) && distance < almostzero && return
 
-        # Now that we know that `distance` is not zero, we can safely call the unsafe
-        # version of the kernel gradient to avoid redundant zero checks.
-        # Note that we use the `neighbor_system` to compute the kernel gradient
-        # to obtain the same force as in the fluid-structure interaction.
+        # The structure-oriented gradient is used by viscosity and adhesion below.
         grad_kernel = smoothing_kernel_grad_unsafe(neighbor_system, pos_diff,
                                                    distance, neighbor)
 
@@ -66,19 +63,29 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
 
         # In fluid-structure interaction, use the "hydrodynamic pressure" of the structure
         # particles corresponding to the chosen boundary model.
-        p_a = current_pressure(v_particle_system, particle_system, particle)
-        p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)
-
-        # Particle and neighbor (and the corresponding systems and particle quantities) are
-        # switched in the following two calls. This yields the exact same pair force as in the
-        # fluid-structure interaction, but with flipped sign because `pos_diff` is reversed.
-        dv_boundary = pressure_acceleration(neighbor_system, particle_system,
-                                            neighbor, particle,
-                                            m_b, m_a, p_b, p_a, rho_b, rho_a,
-                                            pos_diff, distance, grad_kernel,
-                                            system_correction(neighbor_system))
-
-        dv_particle = Ref(dv_boundary)
+        p_fluid = current_pressure(v_neighbor_system, neighbor_system, neighbor)
+        p_boundary = neighbor_pressure(v_particle_system, particle_system, particle,
+                                       p_fluid)
+        p_avg = pair_pressure_offset(neighbor_system, particle_system, neighbor, particle)
+
+        # Reconstruct the fluid-oriented pair exactly as in the fluid-structure interaction.
+        # Corrected gradients are generally not odd, so evaluating the fluid gradient at the
+        # reversed displacement would not yield the reaction force. Instead, compute the fluid
+        # acceleration with the same orientation and apply its exact negative to the structure.
+        fluid_pos_diff = -pos_diff
+        fluid_grad_kernel = smoothing_kernel_grad_unsafe(neighbor_system, fluid_pos_diff,
+                                                         distance, neighbor)
+        dv_fluid_pressure = pressure_acceleration(neighbor_system, particle_system,
+                                                  neighbor, particle,
+                                                  m_b, m_a, p_fluid - p_avg,
+                                                  p_boundary - p_avg, rho_b, rho_a,
+                                                  fluid_pos_diff, distance,
+                                                  fluid_grad_kernel,
+                                                  system_correction(neighbor_system))
+        pressure_correction = interaction_pressure_correction(neighbor_system, rho_b,
+                                                              rho_a)
+
+        dv_particle = Ref(-dv_fluid_pressure * pressure_correction)
         dv_viscosity!(dv_particle, neighbor_system, particle_system,
                       v_neighbor_system, v_particle_system,
                       neighbor, particle, pos_diff, distance,
@@ -101,6 +108,14 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
     return dv
 end
 
+@inline interaction_pressure_correction(system, rho_a, rho_b) = one(rho_a)
+
+@inline function interaction_pressure_correction(system::WeaklyCompressibleSPHSystem,
+                                                 rho_a, rho_b)
+    return free_surface_correction(correction_force(system.correction),
+                                   system, rho_a, rho_b)[2]
+end
+
 @inline function continuity_equation!(drho_particle,
                                       particle_system::AbstractStructureSystem,
                                       neighbor_system::AbstractFluidSystem,
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index a1cd48fb54..7255eb85b2 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -347,10 +347,29 @@ end
     return system.boundary_model.hydrodynamic_mass[particle]
 end
 
-@propagate_inbounds function correction_matrix(system, particle)
+@propagate_inbounds function tlsph_correction_matrix(system, particle)
     extract_smatrix(system.correction_matrix, system, particle)
 end
 
+@inline function hydrodynamic_correction(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles})
+    return correction_gradient(system.boundary_model.correction)
+end
+
+@inline function kernel_correction_coefficient(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
+                                               particle)
+    return system.boundary_model.cache.kernel_correction_coefficient[particle]
+end
+
+@inline function dw_gamma(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
+                          particle)
+    return extract_svector(system.boundary_model.cache.dw_gamma, system, particle)
+end
+
+@inline function correction_matrix(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
+                                   particle)
+    return extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle)
+end
+
 @propagate_inbounds function deformation_gradient(system, particle)
     extract_smatrix(system.deformation_grad, system, particle)
 end
@@ -457,6 +476,13 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode
     return system
 end
 
+function update_density_correction!(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
+                                    v, u, v_ode, u_ode, semi, t)
+    update_density_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 function update_boundary_interpolation!(system::TotalLagrangianSPHSystem, v, u,
                                         v_ode, u_ode, semi, t)
     (; boundary_model) = system
@@ -465,6 +491,13 @@ function update_boundary_interpolation!(system::TotalLagrangianSPHSystem, v, u,
     update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi)
 end
 
+function update_gradient_correction!(system::TotalLagrangianSPHSystem{<:BoundaryModelDummyParticles},
+                                     v, u, v_ode, u_ode, semi, t)
+    update_gradient_correction!(system.boundary_model, system, v, u, v_ode, u_ode, semi)
+
+    return system
+end
+
 @inline function compute_pk1_corrected!(system, semi)
     (; deformation_grad, pk1_rho2, material_density) = system
 
@@ -473,7 +506,7 @@ end
     @threaded semi for particle in eachparticle(system)
         pk1_particle = @inbounds pk1_stress_tensor(system, particle)
         pk1_particle_corrected = pk1_particle *
-                                 @inbounds correction_matrix(system, particle)
+                                 @inbounds tlsph_correction_matrix(system, particle)
         rho2_inv = 1 / @inbounds material_density[particle]^2
 
         for j in 1:ndims(system), i in 1:ndims(system)
@@ -505,7 +538,7 @@ end
         # We are looping over the particles of `system`, so it is guaranteed
         # that `particle` is in bounds of `system`.
         current_coords_a = @inbounds current_coords(system, particle)
-        L_a = @inbounds correction_matrix(system, particle)
+        L_a = @inbounds tlsph_correction_matrix(system, particle)
 
         # Accumulate the contributions over all neighbors before writing
         # to `deformation_grad` to reduce the number of memory writes.
diff --git a/test/callbacks/density_reinit.jl b/test/callbacks/density_reinit.jl
index bb65c1d133..b972e07c21 100644
--- a/test/callbacks/density_reinit.jl
+++ b/test/callbacks/density_reinit.jl
@@ -216,4 +216,73 @@
                                                                         (;
                                                                          systems=(MockNoDensityReinitSystem(:boundary),)))
     end
+
+    @testset "nonuniform multi-system reinitialization with buffer" begin
+        spacing = 0.1
+        kernel = WendlandC6Kernel{2}()
+        smoothing_length = 2spacing
+        state_equation = StateEquationCole(; sound_speed=10.0,
+                                           reference_density=1000.0, exponent=1)
+        initial1 = RectangularShape(spacing, (3, 3), (0.0, 0.0); density=1000.0)
+        initial2 = RectangularShape(spacing, (3, 3), (0.35, 0.0); density=1000.0)
+        system1 = WeaklyCompressibleSPHSystem(initial1; smoothing_kernel=kernel,
+                                              smoothing_length,
+                                              density_calculator=ContinuityDensity(),
+                                              state_equation, buffer_size=2)
+        system2 = WeaklyCompressibleSPHSystem(initial2; smoothing_kernel=kernel,
+                                              smoothing_length,
+                                              density_calculator=ContinuityDensity(),
+                                              state_equation)
+        semi = Semidiscretization(system1, system2; neighborhood_search=nothing,
+                                  parallelization_backend=SerialBackend())
+        callback = DensityReinitializationCallback(system1, semi; interval=1,
+                                                   reinit_initial_solution=false)
+        ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+        semi = ode.p.semi
+        system1, system2 = semi.systems
+        vu_ode = deepcopy(ode.u0)
+        v_ode, u_ode = vu_ode.x
+        v1 = TrixiParticles.wrap_v(v_ode, system1, semi)
+        v2 = TrixiParticles.wrap_v(v_ode, system2, semi)
+        u1 = TrixiParticles.wrap_u(u_ode, system1, semi)
+        u2 = TrixiParticles.wrap_u(u_ode, system2, semi)
+        active1 = collect(TrixiParticles.eachparticle(system1))
+        active2 = collect(TrixiParticles.eachparticle(system2))
+        v1[end, active1] .= range(800.0, 1200.0; length=length(active1))
+        v2[end, active2] .= range(900.0, 1100.0; length=length(active2))
+        density2_before = copy(v2[end, :])
+
+        TrixiParticles.update_nhs!(semi, u_ode)
+        coefficient = zeros(size(v1, 2))
+        TrixiParticles.compute_shepard_coeff!(system1,
+                                              TrixiParticles.current_coordinates(u1,
+                                                                                 system1),
+                                              v_ode, u_ode, semi, coefficient)
+        summation = zeros(size(v1, 2))
+        TrixiParticles.summation_density!(system1, semi, u1, u_ode, summation)
+        expected = summation[active1] ./ coefficient[active1]
+
+        cross_contribution = zeros(size(v1, 2))
+        coords1 = TrixiParticles.current_coordinates(u1, system1)
+        coords2 = TrixiParticles.current_coordinates(u2, system2)
+        TrixiParticles.foreach_point_neighbor(system1, system2, coords1, coords2,
+                                              semi) do particle, neighbor, pos_diff,
+                                                       distance
+            volume = TrixiParticles.hydrodynamic_mass(system2, neighbor) /
+                     TrixiParticles.current_density(v2, system2, neighbor)
+            cross_contribution[particle] += volume *
+                                            TrixiParticles.smoothing_kernel(system1,
+                                                                            distance,
+                                                                            particle)
+        end
+
+        integrator = MockDensityReinitIntegrator((; semi), vu_ode, 0.05)
+        callback.affect!(integrator)
+        inactive1 = setdiff(axes(v1, 2), active1)
+
+        @test v1[end, active1]≈expected rtol=2e-14 atol=2e-14
+        @test all(iszero, v1[end, inactive1])
+        @test v2[end, :] == density2_before
+        @test any(>(0), cross_contribution[active1])
+    end
 end
diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl
index f63bca3b45..03f2b11fe6 100644
--- a/test/examples/gpu.jl
+++ b/test/examples/gpu.jl
@@ -77,6 +77,121 @@ end
     end
 end
 
+@testset verbose=true "Correction lifecycle $TRIXIPARTICLES_TEST_" begin
+    function correction_fluid(kind, initial_condition, smoothing_kernel,
+                              smoothing_length, density_calculator, correction)
+        if kind == :wcsph
+            state_equation = StateEquationCole(; sound_speed=10.0f0,
+                                               reference_density=1000.0f0, exponent=1)
+            return WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel,
+                                               smoothing_length, density_calculator,
+                                               state_equation, correction)
+        end
+
+        return EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel,
+                                           smoothing_length, sound_speed=10.0f0,
+                                           pressure_acceleration=nothing,
+                                           density_calculator, correction)
+    end
+
+    function correction_cache_is_valid(system, backend)
+        return all((:kernel_correction_coefficient, :dw_gamma,
+                    :correction_matrix)) do name
+            hasproperty(system.cache, name) || return true
+            values = getproperty(system.cache, name)
+            return eltype(values) == Float32 && all(isfinite, Array(values)) &&
+                   TrixiParticles.KernelAbstractions.get_backend(values) == backend
+        end
+    end
+
+    function correction_rhs_is_valid(kind, correction, density_calculator, backend)
+        spacing = 0.1f0
+        initial_condition = RectangularShape(spacing, (4, 4), (0.0f0, 0.0f0);
+                                             density=1000.0f0,
+                                             velocity=pos -> SVector(pos[1], -pos[2]),
+                                             coordinates_eltype=Float32)
+        system = correction_fluid(kind, initial_condition, WendlandC6Kernel{2}(),
+                                  2spacing, density_calculator, correction)
+        semi = Semidiscretization(system; neighborhood_search=nothing,
+                                  parallelization_backend=backend)
+        ode = semidiscretize(semi, (0.0f0, 0.1f0); reset_threads=false)
+        v_ode, u_ode = ode.u0.x
+        dv_ode = similar(v_ode)
+        fill!(dv_ode, 0.0f0)
+        TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0f0)
+
+        return eltype(dv_ode) == Float32 && all(isfinite, Array(dv_ode)) &&
+               correction_cache_is_valid(first(ode.p.semi.systems), backend)
+    end
+
+    backend = Main.parallelization_backend
+    for kind in (:wcsph, :edac)
+        @test correction_rhs_is_valid(kind, ShepardKernelCorrection(), SummationDensity(),
+                                      backend)
+        @test correction_rhs_is_valid(kind, KernelCorrection(), ContinuityDensity(),
+                                      backend)
+        @test correction_rhs_is_valid(kind, GradientCorrection(), ContinuityDensity(),
+                                      backend)
+        @test correction_rhs_is_valid(kind, BlendedGradientCorrection(0.4f0),
+                                      ContinuityDensity(), backend)
+        @test correction_rhs_is_valid(kind, MixedKernelGradientCorrection(),
+                                      ContinuityDensity(), backend)
+    end
+
+    spacing = 0.1f0
+    initial_condition = RectangularShape(spacing, (2, 2), (0.0f0, 0.0f0);
+                                         density=1000.0f0,
+                                         coordinates_eltype=Float32)
+    system = correction_fluid(:wcsph, initial_condition, WendlandC6Kernel{2}(),
+                              2spacing, SummationDensity(), ShepardKernelCorrection())
+    semi = Semidiscretization(system; neighborhood_search=nothing,
+                              parallelization_backend=backend)
+    ode = semidiscretize(semi, (0.0f0, 0.1f0); reset_threads=false)
+    coefficient = Adapt.adapt(backend, Float32[0.0, NaN, 1.0, 2.0])
+    TrixiParticles.sanitize_kernel_correction_coefficient!(coefficient,
+                                                           first(ode.p.semi.systems),
+                                                           ode.p.semi)
+    @test Array(coefficient) == Float32[1.0, 1.0, 1.0, 2.0]
+
+    coordinates = Float32[0.0 0.1 0.2; 0.0 0.0 0.0]
+    collinear = InitialCondition(; coordinates, velocity=zeros(Float32, 2, 3),
+                                 density=fill(1000.0f0, 3), particle_spacing=spacing)
+    system = correction_fluid(:wcsph, collinear, WendlandC6Kernel{2}(), 2spacing,
+                              ContinuityDensity(), GradientCorrection())
+    semi = Semidiscretization(system; neighborhood_search=nothing,
+                              parallelization_backend=backend)
+    ode = semidiscretize(semi, (0.0f0, 0.1f0); reset_threads=false)
+    dv_ode = similar(ode.u0.x[1])
+    fill!(dv_ode, 0.0f0)
+    TrixiParticles.kick!(dv_ode, ode.u0.x[1], ode.u0.x[2], ode.p, 0.0f0)
+    matrix = Array(first(ode.p.semi.systems).cache.correction_matrix)
+    identity = Matrix{Float32}(I, 2, 2)
+    @test all(particle -> matrix[:, :, particle] == identity, axes(matrix, 3))
+
+    initial_condition = RectangularShape(spacing, (4, 4), (0.0f0, 0.0f0);
+                                         density=1000.0f0,
+                                         coordinates_eltype=Float32)
+    state_equation = StateEquationCole(; sound_speed=10.0f0,
+                                       reference_density=1000.0f0, exponent=1)
+    system = WeaklyCompressibleSPHSystem(initial_condition;
+                                         smoothing_kernel=WendlandC6Kernel{2}(),
+                                         smoothing_length=2spacing,
+                                         density_calculator=ContinuityDensity(),
+                                         state_equation)
+    semi = Semidiscretization(system; neighborhood_search=nothing,
+                              parallelization_backend=backend)
+    ode = semidiscretize(semi, (0.0f0, 0.1f0); reset_threads=false)
+    v_ode, u_ode = ode.u0.x
+    system = first(ode.p.semi.systems)
+    v = TrixiParticles.wrap_v(v_ode, system, ode.p.semi)
+    u = TrixiParticles.wrap_u(u_ode, system, ode.p.semi)
+    density = collect(range(800.0f0, 1200.0f0; length=size(v, 2)))
+    v[end, :] .= Adapt.adapt(backend, density)
+    TrixiParticles.update_nhs!(ode.p.semi, u_ode)
+    TrixiParticles.reinit_density!(system, v, u, v_ode, u_ode, ode.p.semi)
+    @test all(isfinite, Array(v[end, :]))
+end
+
 @testset verbose=true "Examples $TRIXIPARTICLES_TEST_" begin
     @testset verbose=true "Fluid" begin
         @trixi_testset "fluid/dam_break_2d_gpu.jl Float64" begin
diff --git a/test/general/corrections.jl b/test/general/corrections.jl
new file mode 100644
index 0000000000..f4fee9b772
--- /dev/null
+++ b/test/general/corrections.jl
@@ -0,0 +1,9 @@
+@trixi_testset "Correction Consistency" begin
+    include("corrections/common.jl")
+    include("corrections/lifecycle.jl")
+    include("corrections/shepard.jl")
+    include("corrections/kernel.jl")
+    include("corrections/gradient.jl")
+    include("corrections/mixed.jl")
+    include("corrections/coupling.jl")
+end
diff --git a/test/general/corrections/common.jl b/test/general/corrections/common.jl
new file mode 100644
index 0000000000..703f3e02df
--- /dev/null
+++ b/test/general/corrections/common.jl
@@ -0,0 +1,159 @@
+function correction_setup(correction=nothing; n=9, perturbation=false,
+                          density_calculator=ContinuityDensity(), edac=false,
+                          pressure_acceleration=:default,
+                          velocity=(pos -> SVector(pos[1], pos[2])))
+    particle_spacing = 1.0 / n
+    smoothing_length = 2.0 * particle_spacing
+    smoothing_kernel = WendlandC6Kernel{2}()
+    fluid = RectangularShape(particle_spacing, (n, n), (0.0, 0.0);
+                             density=1000.0, velocity,
+                             coordinates_perturbation=perturbation ? 0.1 : nothing)
+
+    if edac
+        if pressure_acceleration === :default
+            system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length, sound_speed=10.0,
+                                                 density_calculator, correction)
+        else
+            system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length, sound_speed=10.0,
+                                                 density_calculator, correction,
+                                                 pressure_acceleration)
+        end
+    else
+        state_equation = StateEquationCole(; sound_speed=10.0,
+                                           reference_density=1000.0, exponent=1)
+        if pressure_acceleration === :default
+            system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length, density_calculator,
+                                                 state_equation, correction)
+        else
+            system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length, density_calculator,
+                                                 state_equation, correction,
+                                                 pressure_acceleration)
+        end
+    end
+
+    semi = Semidiscretization(system; parallelization_backend=SerialBackend())
+    ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+    v_ode = Array(ode.u0.x[1])
+    u_ode = Array(ode.u0.x[2])
+    semi = ode.p.semi
+    system = first(semi.systems)
+
+    return (; system, semi, v_ode, u_ode, particle_spacing)
+end
+
+function fill_correction_cache!(system, value)
+    for name in (:kernel_correction_coefficient, :dw_gamma, :correction_matrix)
+        hasproperty(system.cache, name) || continue
+        fill!(getproperty(system.cache, name), value)
+    end
+    return system
+end
+
+function update_correction!(setup)
+    (; system, semi, v_ode, u_ode) = setup
+    fill_correction_cache!(system, NaN)
+    TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0)
+    return setup
+end
+
+function correction_moments(setup; field=(pos -> 1.0))
+    (; system, semi, v_ode, u_ode) = setup
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+    coordinates = Array(TrixiParticles.current_coordinates(u, system))
+    values = [field(SVector{2}(view(coordinates, :, particle)))
+              for particle in TrixiParticles.eachparticle(system)]
+    n_particles = TrixiParticles.nparticles(system)
+
+    zeroth_gradient_moment = zeros(2, n_particles)
+    first_gradient_moment = zeros(2, 2, n_particles)
+    direct_gradient = zeros(2, n_particles)
+    difference_gradient = zeros(2, n_particles)
+
+    GC.@preserve v_ode u_ode begin
+        TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates,
+                                              semi) do particle, neighbor, pos_diff,
+                                                       distance
+            pos_diff_ = SVector(pos_diff)
+            volume = TrixiParticles.hydrodynamic_mass(system, neighbor) /
+                     TrixiParticles.current_density(v, system, neighbor)
+            gradient = TrixiParticles.smoothing_kernel_grad(system, pos_diff_, distance,
+                                                            particle)
+            neighbor_offset = -pos_diff_
+
+            for i in 1:2
+                zeroth_gradient_moment[i, particle] += volume * gradient[i]
+                direct_gradient[i, particle] += volume * values[neighbor] * gradient[i]
+                difference_gradient[i,
+                                    particle] += volume *
+                                                 (values[neighbor] - values[particle]) *
+                                                 gradient[i]
+                for j in 1:2
+                    first_gradient_moment[i, j,
+                                          particle] += volume * gradient[i] *
+                                                       neighbor_offset[j]
+                end
+            end
+        end
+    end
+
+    return (; zeroth_gradient_moment, first_gradient_moment, direct_gradient,
+            difference_gradient)
+end
+
+function corner_particle(system)
+    coordinates = TrixiParticles.initial_coordinates(system)
+    return argmin(eachindex(axes(coordinates, 2))) do particle
+        coordinates[1, particle] + coordinates[2, particle]
+    end
+end
+
+function correction_restart_result(correction; edac, density_calculator)
+    direct = correction_setup(correction; edac, density_calculator,
+                              pressure_acceleration=nothing)
+    (; system, semi, v_ode, u_ode) = direct
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+
+    for particle in TrixiParticles.eachparticle(system)
+        v[1, particle] = 0.01particle
+        v[2, particle] = -0.02particle
+        u[1, particle] += 1.0e-3 * sin(particle)
+        u[2, particle] += 1.0e-3 * cos(particle)
+    end
+    if edac
+        v[3, :] .= range(1.0, 2.0; length=size(v, 2))
+    end
+    if density_calculator isa ContinuityDensity
+        v[end, :] .= range(900.0, 1100.0; length=size(v, 2))
+    end
+
+    dv_direct = zero(v_ode)
+    TrixiParticles.kick!(dv_direct, v_ode, u_ode,
+                         (; semi, split_integration_data=nothing), 0.0)
+
+    restarted = correction_setup(correction; edac, density_calculator,
+                                 pressure_acceleration=nothing)
+    mock_solution = (; u=[(; x=(copy(v_ode), copy(u_ode)))])
+    restart_with!(restarted.semi, mock_solution; reset_threads=false)
+    ode_restart = semidiscretize(restarted.semi, (0.0, 1.0); reset_threads=false)
+    v_restart = Array(ode_restart.u0.x[1])
+    u_restart = Array(ode_restart.u0.x[2])
+    dv_restart = zero(v_restart)
+    TrixiParticles.kick!(dv_restart, v_restart, u_restart,
+                         (; semi=ode_restart.p.semi, split_integration_data=nothing), 0.0)
+
+    cache = first(ode_restart.p.semi.systems).cache
+    cache_finite = all((:kernel_correction_coefficient, :dw_gamma,
+                        :correction_matrix)) do name
+        return !hasproperty(cache, name) || all(isfinite, getproperty(cache, name))
+    end
+
+    return (; state_equal=v_restart == v_ode && u_restart == u_ode,
+            rhs_equal=isapprox(dv_restart, dv_direct; rtol=2e-13, atol=2e-13),
+            cache_finite)
+end
diff --git a/test/general/corrections/coupling.jl b/test/general/corrections/coupling.jl
new file mode 100644
index 0000000000..956c6e7280
--- /dev/null
+++ b/test/general/corrections/coupling.jl
@@ -0,0 +1,234 @@
+struct CustomForceCorrection end
+
+function TrixiParticles.free_surface_correction(::CustomForceCorrection,
+                                                particle_system, rho_a, rho_b)
+    return 2, 3, 4
+end
+
+@testset "Correction role routing" begin
+    setup = correction_setup()
+    correction = CustomForceCorrection()
+    selected = TrixiParticles.correction_force(correction)
+    @test selected === correction
+    @test TrixiParticles.free_surface_correction(selected, setup.system,
+                                                 1000.0, 1000.0) == (2, 3, 4)
+end
+
+@testset "Supported pressure variation matrix" begin
+    function set_pressure_field!(setup, edac)
+        pressure = range(1.0, 2.0; length=TrixiParticles.nparticles(setup.system))
+        if edac
+            v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi)
+            v[3, :] .= pressure
+        elseif setup.system.density_calculator isa ContinuityDensity
+            v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi)
+            v[end, :] .= 1000.0 .+ pressure
+        else
+            setup.system.pressure .= pressure
+        end
+        return setup
+    end
+
+    summation_corrections = (nothing, ShepardKernelCorrection(), KernelCorrection(),
+                             GradientCorrection(), BlendedGradientCorrection(0.5),
+                             MixedKernelGradientCorrection())
+    continuity_corrections = (nothing, KernelCorrection(), GradientCorrection(),
+                              BlendedGradientCorrection(0.5),
+                              MixedKernelGradientCorrection())
+    summation_pressure = (nothing,
+                          TrixiParticles.pressure_acceleration_summation_density,
+                          TrixiParticles.inter_particle_averaged_pressure)
+    continuity_pressure = (nothing,
+                           TrixiParticles.pressure_acceleration_continuity_density,
+                           TrixiParticles.inter_particle_averaged_pressure)
+
+    for edac in (false, true), correction in summation_corrections,
+        pressure_acceleration in summation_pressure
+        setup = correction_setup(correction; n=4, edac,
+                                 density_calculator=SummationDensity(),
+                                 pressure_acceleration)
+        set_pressure_field!(setup, edac)
+        dv_ode = zero(setup.v_ode)
+        TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode,
+                             (; semi=setup.semi, split_integration_data=nothing), 0.0)
+        @test all(isfinite, dv_ode)
+        @test any(!iszero, view(dv_ode, 1:2, :))
+    end
+
+    for edac in (false, true), correction in continuity_corrections,
+        pressure_acceleration in continuity_pressure
+        setup = correction_setup(correction; n=4, edac,
+                                 density_calculator=ContinuityDensity(),
+                                 pressure_acceleration)
+        set_pressure_field!(setup, edac)
+        dv_ode = zero(setup.v_ode)
+        TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode,
+                             (; semi=setup.semi, split_integration_data=nothing), 0.0)
+        @test all(isfinite, dv_ode)
+        @test any(!iszero, view(dv_ode, 1:2, :))
+    end
+
+    for edac in (false, true)
+        setup = correction_setup(; n=4, edac,
+                                 density_calculator=ContinuityDensity(),
+                                 pressure_acceleration=tensile_instability_control)
+        set_pressure_field!(setup, edac)
+        dv_ode = zero(setup.v_ode)
+        TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode,
+                             (; semi=setup.semi, split_integration_data=nothing), 0.0)
+        @test all(isfinite, dv_ode)
+        @test any(!iszero, view(dv_ode, 1:2, :))
+
+        for correction in continuity_corrections[2:end]
+            @test_throws ArgumentError correction_setup(correction; n=4, edac,
+                                                        density_calculator=ContinuityDensity(),
+                                                        pressure_acceleration=tensile_instability_control)
+        end
+    end
+end
+
+@testset "Corrected structure coupling" begin
+    function cache_is_finite(system)
+        cache = system isa Union{RigidBodySystem, TotalLagrangianSPHSystem} ?
+                system.boundary_model.cache : system.cache
+        return all((:kernel_correction_coefficient, :dw_gamma,
+                    :correction_matrix)) do name
+            !hasproperty(cache, name) || all(isfinite, getproperty(cache, name))
+        end
+    end
+
+    function coupled_result(kind, structure_kind, correction;
+                            boundary_correction=correction, reverse_order=false,
+                            average_pressure_reduction=false)
+        spacing = 0.1
+        density = 1000.0
+        kernel = WendlandC6Kernel{2}()
+        smoothing_length = 2spacing
+        state_equation = kind == :wcsph ?
+                         StateEquationCole(; sound_speed=10.0,
+                                           reference_density=density,
+                                           exponent=1) : nothing
+        fluid_initial = RectangularShape(spacing, (4, 3), (0.0, 0.0); density)
+        fluid_initial.coordinates[1, 2] += 0.013
+        fluid_initial.coordinates[2, 5] -= 0.009
+        if kind == :wcsph
+            fluid = WeaklyCompressibleSPHSystem(fluid_initial; smoothing_kernel=kernel,
+                                                smoothing_length,
+                                                density_calculator=SummationDensity(),
+                                                state_equation, correction)
+        else
+            fluid = EntropicallyDampedSPHSystem(fluid_initial;
+                                                smoothing_kernel=kernel,
+                                                smoothing_length,
+                                                sound_speed=10.0,
+                                                density_calculator=SummationDensity(),
+                                                correction,
+                                                average_pressure_reduction)
+        end
+
+        structure_initial = RectangularShape(spacing, (4, 2), (0.0, -0.2);
+                                             density=1200.0)
+        hydrodynamic_density = fill(density, TrixiParticles.nparticles(structure_initial))
+        hydrodynamic_mass = hydrodynamic_density .* spacing^2
+        boundary_model = BoundaryModelDummyParticles(hydrodynamic_density,
+                                                     hydrodynamic_mass,
+                                                     AdamiPressureExtrapolation(),
+                                                     kernel, smoothing_length;
+                                                     state_equation,
+                                                     correction=boundary_correction)
+        if structure_kind == :rigid
+            structure = RigidBodySystem(structure_initial; boundary_model,
+                                        particle_spacing=spacing)
+        else
+            structure = TotalLagrangianSPHSystem(structure_initial;
+                                                 smoothing_kernel=kernel,
+                                                 smoothing_length,
+                                                 young_modulus=0.0,
+                                                 poisson_ratio=0.0,
+                                                 boundary_model)
+        end
+
+        systems = reverse_order ? (structure, fluid) : (fluid, structure)
+        semi = Semidiscretization(systems...; neighborhood_search=nothing,
+                                  parallelization_backend=SerialBackend())
+        ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+        v_ode = Array(ode.u0.x[1])
+        u_ode = Array(ode.u0.x[2])
+        fluid = only(system
+                     for system in ode.p.semi.systems
+                     if system isa Union{WeaklyCompressibleSPHSystem,
+                              EntropicallyDampedSPHSystem})
+        structure = only(system
+                         for system in ode.p.semi.systems
+                         if system isa Union{RigidBodySystem,
+                                             TotalLagrangianSPHSystem})
+        v_fluid = TrixiParticles.wrap_v(v_ode, fluid, ode.p.semi)
+        if kind == :wcsph
+            v_fluid .= 0.0
+        else
+            v_fluid[3, :] .= range(1.0, 2.0; length=size(v_fluid, 2))
+        end
+
+        dv_ode = zero(v_ode)
+        TrixiParticles.kick!(dv_ode, v_ode, u_ode,
+                             (; semi=ode.p.semi, split_integration_data=nothing), 0.0)
+        dv_fluid = TrixiParticles.wrap_v(dv_ode, fluid, ode.p.semi)
+        fluid_force = vec(sum(fluid.mass' .* view(dv_fluid, 1:2, :); dims=2))
+        if structure_kind == :rigid
+            structure_force = vec(sum(structure.force_per_particle; dims=2))
+        else
+            dv_structure = TrixiParticles.wrap_v(dv_ode, structure, ode.p.semi)
+            structure_force = vec(sum(structure.mass' .* view(dv_structure, 1:2, :);
+                                      dims=2))
+        end
+        force_scale = norm(fluid_force) + norm(structure_force)
+        relative_residual = norm(fluid_force + structure_force) /
+                            max(force_scale, eps())
+
+        fluid_correction = hasproperty(fluid.cache, :correction_matrix) ?
+                           copy(fluid.cache.correction_matrix) : nothing
+        structure_cache = structure.boundary_model.cache
+        structure_correction = hasproperty(structure_cache, :correction_matrix) ?
+                               copy(structure_cache.correction_matrix) : nothing
+
+        return (; relative_residual, force_scale, fluid_force, structure_force,
+                fluid_rhs=copy(dv_fluid), fluid_correction, structure_correction,
+                finite=all(isfinite, dv_ode) && cache_is_finite(fluid) &&
+                       cache_is_finite(structure))
+    end
+
+    corrections = (KernelCorrection(), GradientCorrection(),
+                   BlendedGradientCorrection(0.4), MixedKernelGradientCorrection())
+    for kind in (:wcsph, :edac), structure_kind in (:rigid, :tlsph),
+        correction in corrections
+        result = coupled_result(kind, structure_kind, correction)
+        @test result.finite
+        @test result.force_scale > eps()
+        @test result.relative_residual < 2e-13
+    end
+
+    for structure_kind in (:rigid, :tlsph)
+        forward = coupled_result(:edac, structure_kind, GradientCorrection())
+        reverse = coupled_result(:edac, structure_kind, GradientCorrection();
+                                 reverse_order=true)
+        @test forward.fluid_correction≈reverse.fluid_correction rtol=5e-13 atol=5e-13
+        @test forward.structure_correction≈reverse.structure_correction rtol=5e-13 atol=5e-13
+        @test forward.fluid_rhs≈reverse.fluid_rhs rtol=1e-11 atol=1e-10
+        @test forward.fluid_force≈reverse.fluid_force rtol=1e-11 atol=1e-10
+        @test forward.structure_force≈reverse.structure_force rtol=1e-11 atol=1e-10
+    end
+
+    result = coupled_result(:edac, :rigid, GradientCorrection();
+                            average_pressure_reduction=true)
+    @test result.force_scale > eps()
+    @test result.relative_residual < 2e-13
+
+    for structure_kind in (:rigid, :tlsph),
+        correction in (KernelCorrection(), MixedKernelGradientCorrection())
+        result = coupled_result(:wcsph, structure_kind, correction;
+                                boundary_correction=nothing)
+        @test result.finite
+        @test result.force_scale > eps()
+        @test result.relative_residual < 2e-13
+    end
+end
diff --git a/test/general/corrections/gradient.jl b/test/general/corrections/gradient.jl
new file mode 100644
index 0000000000..7ef9dba05c
--- /dev/null
+++ b/test/general/corrections/gradient.jl
@@ -0,0 +1,138 @@
+@testset "Gradient and blended corrections" begin
+    identity_matrix = Matrix{Float64}(I, 2, 2)
+    linear_field(pos) = 2.0 + 3.0 * pos[1] - 2.0 * pos[2]
+    exact_gradient = [3.0, -2.0]
+
+    @test_throws ArgumentError BlendedGradientCorrection(-0.1)
+    @test_throws ArgumentError BlendedGradientCorrection(1.1)
+
+    for correction in (GradientCorrection(), BlendedGradientCorrection(0.4)),
+        edac in (false, true)
+        setup = correction_setup(correction; edac, pressure_acceleration=nothing)
+        update_correction!(setup)
+        @test all(isfinite, setup.system.cache.correction_matrix)
+    end
+
+    for perturbation in (false, true)
+        raw_setup = update_correction!(correction_setup(nothing; perturbation))
+        raw_moments = correction_moments(raw_setup; field=linear_field)
+
+        gradient_setup = update_correction!(correction_setup(GradientCorrection();
+                                                             perturbation))
+        gradient_moments = correction_moments(gradient_setup; field=linear_field)
+        @test maximum(particle -> norm(gradient_moments.first_gradient_moment[:, :,
+                                                                              particle] -
+                                       identity_matrix),
+                      TrixiParticles.eachparticle(gradient_setup.system)) < 2e-12
+        @test maximum(particle -> norm(gradient_moments.difference_gradient[:, particle] -
+                                       exact_gradient),
+                      TrixiParticles.eachparticle(gradient_setup.system)) < 5e-12
+
+        blending_factor = 0.4
+        blended_setup = update_correction!(correction_setup(BlendedGradientCorrection(blending_factor);
+                                                            perturbation))
+        blended_moments = correction_moments(blended_setup; field=linear_field)
+        expected = (1 - blending_factor) * raw_moments.first_gradient_moment
+        for particle in TrixiParticles.eachparticle(blended_setup.system)
+            expected[:, :, particle] .+= blending_factor * identity_matrix
+        end
+        @test maximum(abs, blended_moments.first_gradient_moment - expected) < 2e-12
+
+        corner = corner_particle(raw_setup.system)
+        @test norm(raw_moments.first_gradient_moment[:, :, corner] - identity_matrix) > 1e-2
+    end
+
+    density32 = fill(1000.0f0, 4)
+    mass32 = fill(10.0f0, 4)
+    state_equation32 = StateEquationCole(; sound_speed=10.0f0,
+                                         reference_density=1000.0f0, exponent=1)
+    boundary = BoundaryModelDummyParticles(density32, mass32, SummationDensity(),
+                                           WendlandC6Kernel{2}(), 0.2f0;
+                                           state_equation=state_equation32,
+                                           correction=GradientCorrection())
+    @test eltype(boundary.cache.correction_matrix) == Float32
+
+    particle_spacing = 0.25
+    particles = RectangularShape(particle_spacing, (4, 4, 4), (0.0, 0.0, 0.0);
+                                 density=1000.0)
+    state_equation = StateEquationCole(; sound_speed=10.0,
+                                       reference_density=1000.0, exponent=1)
+    system = WeaklyCompressibleSPHSystem(particles;
+                                         smoothing_kernel=WendlandC6Kernel{3}(),
+                                         smoothing_length=2particle_spacing,
+                                         density_calculator=ContinuityDensity(),
+                                         state_equation,
+                                         correction=GradientCorrection())
+    semi = Semidiscretization(system; parallelization_backend=SerialBackend())
+    ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+    v_ode = Array(ode.u0.x[1])
+    u_ode = Array(ode.u0.x[2])
+    system = first(ode.p.semi.systems)
+    TrixiParticles.update_systems_and_nhs(v_ode, u_ode, ode.p.semi, 0.0)
+    v = TrixiParticles.wrap_v(v_ode, system, ode.p.semi)
+    u = TrixiParticles.wrap_u(u_ode, system, ode.p.semi)
+    coordinates = Array(TrixiParticles.current_coordinates(u, system))
+    first_moment = zeros(3, 3)
+    GC.@preserve v_ode u_ode begin
+        TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates,
+                                              ode.p.semi;
+                                              points=1:1) do particle,
+                                                             neighbor,
+                                                             pos_diff,
+                                                             distance
+            volume = TrixiParticles.hydrodynamic_mass(system, neighbor) /
+                     TrixiParticles.current_density(v, system, neighbor)
+            gradient = TrixiParticles.smoothing_kernel_grad(system, SVector(pos_diff),
+                                                            distance, particle)
+            for j in 1:3, i in 1:3
+                first_moment[i, j] -= volume * gradient[i] * pos_diff[j]
+            end
+        end
+    end
+    @test first_moment ≈ Matrix{Float64}(I, 3, 3) atol = 3e-12
+
+    for y_offset in (0.0, 1.0e-12)
+        coordinates = [0.0 0.1 0.2; 0.0 y_offset 0.0]
+        initial = InitialCondition(; coordinates, velocity=zeros(2, 3),
+                                   density=fill(1000.0, 3), particle_spacing=0.1)
+        system = WeaklyCompressibleSPHSystem(initial;
+                                             smoothing_kernel=WendlandC6Kernel{2}(),
+                                             smoothing_length=0.2,
+                                             density_calculator=ContinuityDensity(),
+                                             state_equation,
+                                             correction=GradientCorrection())
+        semi = Semidiscretization(system; parallelization_backend=SerialBackend())
+        ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+        v_ode = Array(ode.u0.x[1])
+        u_ode = Array(ode.u0.x[2])
+        system = first(ode.p.semi.systems)
+        TrixiParticles.update_systems_and_nhs(v_ode, u_ode, ode.p.semi, 0.0)
+        for particle in TrixiParticles.eachparticle(system)
+            @test TrixiParticles.correction_matrix(system, particle) == I
+        end
+    end
+
+    analytic_density_rate = -2000.0
+    errors = Dict{Any, Float64}()
+    for correction in (nothing, GradientCorrection(), BlendedGradientCorrection(0.4))
+        setup = correction_setup(correction)
+        dv_ode = zero(setup.v_ode)
+        TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode,
+                             (; semi=setup.semi, split_integration_data=nothing), 0.0)
+        dv = TrixiParticles.wrap_v(dv_ode, setup.system, setup.semi)
+        error = dv[end, :] .- analytic_density_rate
+        errors[correction] = sqrt(sum(abs2, error) / length(error))
+    end
+    @test errors[GradientCorrection()] < 2e-10
+    @test errors[BlendedGradientCorrection(0.4)] < errors[nothing]
+    @test errors[nothing] > 1.0
+
+    for correction in (GradientCorrection(), BlendedGradientCorrection(0.4)),
+        edac in (false, true),
+        density_calculator in (SummationDensity(), ContinuityDensity())
+        result = correction_restart_result(correction; edac, density_calculator)
+        @test result.state_equal
+        @test result.rhs_equal
+        @test result.cache_finite
+    end
+end
diff --git a/test/general/corrections/kernel.jl b/test/general/corrections/kernel.jl
new file mode 100644
index 0000000000..ec0acac585
--- /dev/null
+++ b/test/general/corrections/kernel.jl
@@ -0,0 +1,35 @@
+@testset "Kernel correction" begin
+    for edac in (false, true)
+        setup = correction_setup(KernelCorrection(); edac,
+                                 pressure_acceleration=nothing)
+        update_correction!(setup)
+
+        @test all(isfinite, setup.system.cache.kernel_correction_coefficient)
+        @test all(isfinite, setup.system.cache.dw_gamma)
+    end
+
+    for perturbation in (false, true)
+        setup = update_correction!(correction_setup(KernelCorrection(); perturbation))
+        moments = correction_moments(setup)
+        @test maximum(abs, moments.zeroth_gradient_moment) < 2e-12
+    end
+
+    density32 = fill(1000.0f0, 4)
+    mass32 = fill(10.0f0, 4)
+    state_equation = StateEquationCole(; sound_speed=10.0f0,
+                                       reference_density=1000.0f0, exponent=1)
+    boundary = BoundaryModelDummyParticles(density32, mass32, SummationDensity(),
+                                           WendlandC6Kernel{2}(), 0.2f0;
+                                           state_equation,
+                                           correction=KernelCorrection())
+    @test eltype(boundary.cache.dw_gamma) == Float32
+
+    for edac in (false, true),
+        density_calculator in (SummationDensity(),
+                               ContinuityDensity())
+        result = correction_restart_result(KernelCorrection(); edac, density_calculator)
+        @test result.state_equal
+        @test result.rhs_equal
+        @test result.cache_finite
+    end
+end
diff --git a/test/general/corrections/lifecycle.jl b/test/general/corrections/lifecycle.jl
new file mode 100644
index 0000000000..1c871d92e7
--- /dev/null
+++ b/test/general/corrections/lifecycle.jl
@@ -0,0 +1,118 @@
+@testset "Cross-system update ordering" begin
+    function ordered_correction_result(reverse_order; edac)
+        spacing = 0.1
+        smoothing_length = 2spacing
+        smoothing_kernel = WendlandC6Kernel{2}()
+        density = 1000.0
+        velocity(pos) = SVector(0.1 + pos[1], -0.2 - pos[2])
+        pressure(pos) = 1.0 + 2pos[1] - pos[2]
+
+        gradient_initial = RectangularShape(spacing, (3, 3), (0.0, 0.0);
+                                            density, velocity, pressure)
+        shepard_initial = RectangularShape(spacing, (3, 3), (0.05, 0.025);
+                                           density, velocity, pressure)
+
+        if edac
+            gradient_system = EntropicallyDampedSPHSystem(gradient_initial;
+                                                          smoothing_kernel,
+                                                          smoothing_length,
+                                                          sound_speed=10.0,
+                                                          pressure_acceleration=nothing,
+                                                          density_calculator=ContinuityDensity(),
+                                                          correction=GradientCorrection())
+            shepard_system = EntropicallyDampedSPHSystem(shepard_initial;
+                                                         smoothing_kernel,
+                                                         smoothing_length,
+                                                         sound_speed=10.0,
+                                                         density_calculator=SummationDensity(),
+                                                         correction=ShepardKernelCorrection())
+        else
+            state_equation = StateEquationCole(; sound_speed=10.0,
+                                               reference_density=density, exponent=1)
+            gradient_system = WeaklyCompressibleSPHSystem(gradient_initial;
+                                                          smoothing_kernel,
+                                                          smoothing_length,
+                                                          state_equation,
+                                                          density_calculator=ContinuityDensity(),
+                                                          correction=GradientCorrection())
+            shepard_system = WeaklyCompressibleSPHSystem(shepard_initial;
+                                                         smoothing_kernel,
+                                                         smoothing_length,
+                                                         state_equation,
+                                                         density_calculator=SummationDensity(),
+                                                         correction=ShepardKernelCorrection())
+        end
+
+        systems = reverse_order ? (shepard_system, gradient_system) :
+                  (gradient_system, shepard_system)
+        semi = Semidiscretization(systems...; neighborhood_search=nothing,
+                                  parallelization_backend=SerialBackend())
+        ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+        v_ode = Array(ode.u0.x[1])
+        u_ode = Array(ode.u0.x[2])
+        dv_ode = zero(v_ode)
+        TrixiParticles.kick!(dv_ode, v_ode, u_ode,
+                             (; semi=ode.p.semi, split_integration_data=nothing), 0.0)
+
+        gradient_system = only(system
+                               for system in ode.p.semi.systems
+                               if system.correction isa GradientCorrection)
+        shepard_system = only(system
+                              for system in ode.p.semi.systems
+                              if system.correction isa ShepardKernelCorrection)
+        v_gradient = TrixiParticles.wrap_v(v_ode, gradient_system, ode.p.semi)
+        v_shepard = TrixiParticles.wrap_v(v_ode, shepard_system, ode.p.semi)
+        dv_gradient = TrixiParticles.wrap_v(dv_ode, gradient_system, ode.p.semi)
+        dv_shepard = TrixiParticles.wrap_v(dv_ode, shepard_system, ode.p.semi)
+
+        return (;
+                gradient_density=copy(TrixiParticles.current_density(v_gradient,
+                                                                     gradient_system)),
+                shepard_density=copy(TrixiParticles.current_density(v_shepard,
+                                                                    shepard_system)),
+                gradient_pressure=copy(TrixiParticles.current_pressure(v_gradient,
+                                                                       gradient_system)),
+                shepard_pressure=copy(TrixiParticles.current_pressure(v_shepard,
+                                                                      shepard_system)),
+                correction_matrix=copy(gradient_system.cache.correction_matrix),
+                shepard_coefficient=copy(shepard_system.cache.kernel_correction_coefficient),
+                gradient_rhs=copy(dv_gradient), shepard_rhs=copy(dv_shepard))
+    end
+
+    for edac in (false, true)
+        forward = ordered_correction_result(false; edac)
+        reverse = ordered_correction_result(true; edac)
+
+        @test forward.gradient_density≈reverse.gradient_density rtol=5e-13 atol=5e-13
+        @test forward.shepard_density≈reverse.shepard_density rtol=5e-13 atol=5e-13
+        @test forward.gradient_pressure≈reverse.gradient_pressure rtol=5e-13 atol=5e-13
+        @test forward.shepard_pressure≈reverse.shepard_pressure rtol=5e-13 atol=5e-13
+        @test forward.correction_matrix≈reverse.correction_matrix rtol=5e-13 atol=5e-13
+        @test forward.shepard_coefficient≈reverse.shepard_coefficient rtol=5e-13 atol=5e-13
+        @test forward.gradient_rhs≈reverse.gradient_rhs rtol=1e-11 atol=1e-10
+        @test forward.shepard_rhs≈reverse.shepard_rhs rtol=1e-11 atol=1e-10
+    end
+end
+
+@testset "Boundary density before pressure" begin
+    n = 5
+    particle_spacing = 1.0 / n
+    smoothing_kernel = WendlandC6Kernel{2}()
+    particles = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); density=1000.0)
+    state_equation = StateEquationCole(; sound_speed=10.0,
+                                       reference_density=1000.0, exponent=1)
+    boundary_model = BoundaryModelDummyParticles(particles.density, particles.mass,
+                                                 SummationDensity(), smoothing_kernel,
+                                                 2particle_spacing; state_equation,
+                                                 correction=ShepardKernelCorrection())
+    boundary = WallBoundarySystem(particles, boundary_model)
+    semi = Semidiscretization(boundary; parallelization_backend=SerialBackend())
+    ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+    v_ode = Array(ode.u0.x[1])
+    u_ode = Array(ode.u0.x[2])
+    boundary = first(ode.p.semi.systems)
+    TrixiParticles.update_systems_and_nhs(v_ode, u_ode, ode.p.semi, 0.0)
+
+    @test boundary.boundary_model.pressure ≈
+          state_equation.(boundary.boundary_model.cache.density)
+end
diff --git a/test/general/corrections/mixed.jl b/test/general/corrections/mixed.jl
new file mode 100644
index 0000000000..dc84c2b83f
--- /dev/null
+++ b/test/general/corrections/mixed.jl
@@ -0,0 +1,56 @@
+@testset "Mixed kernel-gradient correction" begin
+    identity_matrix = Matrix{Float64}(I, 2, 2)
+    linear_field(pos) = 2.0 + 3.0 * pos[1] - 2.0 * pos[2]
+    exact_gradient = [3.0, -2.0]
+
+    for edac in (false, true)
+        setup = correction_setup(MixedKernelGradientCorrection(); edac,
+                                 pressure_acceleration=nothing)
+        update_correction!(setup)
+        @test all(isfinite, setup.system.cache.kernel_correction_coefficient)
+        @test all(isfinite, setup.system.cache.dw_gamma)
+        @test all(isfinite, setup.system.cache.correction_matrix)
+    end
+
+    for perturbation in (false, true)
+        setup = update_correction!(correction_setup(MixedKernelGradientCorrection();
+                                                    perturbation))
+        moments = correction_moments(setup; field=linear_field)
+        @test maximum(abs, moments.zeroth_gradient_moment) < 3e-12
+        @test maximum(particle -> norm(moments.first_gradient_moment[:, :, particle] -
+                                       identity_matrix),
+                      TrixiParticles.eachparticle(setup.system)) < 3e-12
+        @test maximum(particle -> norm(moments.direct_gradient[:, particle] -
+                                       exact_gradient),
+                      TrixiParticles.eachparticle(setup.system)) < 1e-11
+    end
+
+    setup = correction_setup(MixedKernelGradientCorrection())
+    dv_ode = zero(setup.v_ode)
+    TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode,
+                         (; semi=setup.semi, split_integration_data=nothing), 0.0)
+    dv = TrixiParticles.wrap_v(dv_ode, setup.system, setup.semi)
+    density_error = dv[end, :] .+ 2000.0
+    @test sqrt(sum(abs2, density_error) / length(density_error)) < 2e-10
+
+    for edac in (false, true),
+        density_calculator in (SummationDensity(),
+                               ContinuityDensity())
+        result = correction_restart_result(MixedKernelGradientCorrection();
+                                           edac, density_calculator)
+        @test result.state_equal
+        @test result.rhs_equal
+        @test result.cache_finite
+    end
+
+    density32 = fill(1000.0f0, 4)
+    mass32 = fill(10.0f0, 4)
+    state_equation = StateEquationCole(; sound_speed=10.0f0,
+                                       reference_density=1000.0f0, exponent=1)
+    boundary = BoundaryModelDummyParticles(density32, mass32, SummationDensity(),
+                                           WendlandC6Kernel{2}(), 0.2f0;
+                                           state_equation,
+                                           correction=MixedKernelGradientCorrection())
+    @test eltype(boundary.cache.dw_gamma) == Float32
+    @test eltype(boundary.cache.correction_matrix) == Float32
+end
diff --git a/test/general/corrections/shepard.jl b/test/general/corrections/shepard.jl
new file mode 100644
index 0000000000..ece7ea0bf3
--- /dev/null
+++ b/test/general/corrections/shepard.jl
@@ -0,0 +1,55 @@
+@testset "Shepard correction" begin
+    setup = update_correction!(correction_setup(ShepardKernelCorrection();
+                                                density_calculator=SummationDensity()))
+    v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi)
+    density = TrixiParticles.current_density(v, setup.system)
+
+    @test all(isfinite, setup.system.cache.kernel_correction_coefficient)
+    @test setup.system.pressure ≈ setup.system.state_equation.(density)
+
+    setup_edac = update_correction!(correction_setup(ShepardKernelCorrection();
+                                                     density_calculator=SummationDensity(),
+                                                     edac=true))
+    @test all(isfinite, setup_edac.system.cache.kernel_correction_coefficient)
+    @test all(isfinite, setup_edac.system.cache.density)
+
+    coefficients = ones(TrixiParticles.nparticles(setup.system))
+    coefficients[1] = 0.0
+    coefficients[2] = NaN
+    TrixiParticles.sanitize_kernel_correction_coefficient!(coefficients, setup.system,
+                                                           setup.semi)
+    @test coefficients[1:2] == ones(2)
+end
+
+@testset "Shepard partition of unity" begin
+    setup = correction_setup(nothing)
+    (; system, semi, v_ode, u_ode) = setup
+    v = TrixiParticles.wrap_v(v_ode, system, semi)
+    u = TrixiParticles.wrap_u(u_ode, system, semi)
+    coefficient = zeros(TrixiParticles.nparticles(system))
+    numerator = zero(coefficient)
+
+    TrixiParticles.compute_shepard_coeff!(system,
+                                          TrixiParticles.current_coordinates(u, system),
+                                          v_ode, u_ode, semi, coefficient)
+    coordinates = TrixiParticles.current_coordinates(u, system)
+    TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates,
+                                          semi) do particle, neighbor, pos_diff, distance
+        numerator[particle] += TrixiParticles.hydrodynamic_mass(system, neighbor) *
+                               TrixiParticles.smoothing_kernel(system, distance, particle)
+    end
+
+    @test numerator ./ coefficient ≈ fill(1000.0, length(numerator)) atol = 2e-12
+    @test TrixiParticles.current_density(v, system) == fill(1000.0, length(numerator))
+end
+
+@testset "Continuity density reinitialization" begin
+    setup = correction_setup()
+    v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi)
+    u = TrixiParticles.wrap_u(setup.u_ode, setup.system, setup.semi)
+    TrixiParticles.reinit_density!(setup.system, v, u, setup.v_ode, setup.u_ode,
+                                   setup.semi)
+
+    @test TrixiParticles.current_density(v, setup.system) ≈ fill(1000.0, 81) atol = 2e-12
+    @test maximum(abs, setup.system.pressure) < 2e-10
+end
diff --git a/test/general/general.jl b/test/general/general.jl
index acb07de8b9..2ab768000e 100644
--- a/test/general/general.jl
+++ b/test/general/general.jl
@@ -1,6 +1,7 @@
 include("initial_condition.jl")
 include("smoothing_kernels.jl")
 include("density_calculator.jl")
+include("corrections.jl")
 include("semidiscretization.jl")
 include("interpolation.jl")
 include("buffer.jl")
diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl
index 5b92d133ba..9b7ddf3cb9 100644
--- a/test/general/semidiscretization.jl
+++ b/test/general/semidiscretization.jl
@@ -342,7 +342,7 @@
             u = TrixiParticles.wrap_u(u_ode, system, semi)
 
             TrixiParticles.compute_correction_values!(system,
-                                                      TrixiParticles.system_correction(system),
+                                                      TrixiParticles.correction_density(system.correction),
                                                       u, v_ode, u_ode, semi)
 
             return copy(system.cache.kernel_correction_coefficient), semi
diff --git a/test/schemes/fluid/pressure_acceleration.jl b/test/schemes/fluid/pressure_acceleration.jl
index 80de3d9a66..c3d18bcbe3 100644
--- a/test/schemes/fluid/pressure_acceleration.jl
+++ b/test/schemes/fluid/pressure_acceleration.jl
@@ -12,6 +12,61 @@
         @test f_2 == TrixiParticles.pressure_acceleration_continuity_density
     end
 
+    @testset "Algebraic formulations and asymmetric conservation" begin
+        m_a, m_b = 1.2, 0.8
+        rho_a, rho_b = 1000.0, 980.0
+        p_a, p_b = 2.0, 3.0
+        W_a = SVector(0.2, -0.1)
+        W_b = -W_a
+        W_b_asymmetric = SVector(-0.13, 0.17)
+
+        summation = TrixiParticles.pressure_acceleration_summation_density
+        continuity = TrixiParticles.pressure_acceleration_continuity_density
+        interparticle = TrixiParticles.inter_particle_averaged_pressure
+
+        @test summation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈
+              -m_b * (p_a / rho_a^2 + p_b / rho_b^2) * W_a
+        @test continuity(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈
+              -m_b * (p_a + p_b) / (rho_a * rho_b) * W_a
+
+        volume_term = ((m_a / rho_a)^2 + (m_b / rho_b)^2) / m_a
+        pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b)
+        @test interparticle(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈
+              -volume_term * pressure_tilde * W_a
+        @test tensile_instability_control(m_a, m_b, rho_a, rho_b, -p_a, p_b, W_a) ≈
+              -m_b * (p_a + p_b) / (rho_a * rho_b) * W_a
+
+        for pressure_formulation in (summation, continuity, interparticle)
+            # Asymmetric formulations are selected based on the configured correction and
+            # must reduce to the symmetric formulation when a pair has `W_b == -W_a`.
+            symmetric = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a)
+            asymmetric = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b,
+                                              W_a, W_b)
+            @test asymmetric ≈ symmetric
+            @test pressure_formulation(m_a, m_b, rho_a, rho_b, 0.0, 0.0,
+                                       W_a, W_b) == zero(W_a)
+
+            acceleration_a = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b,
+                                                  W_a, W_b)
+            acceleration_b = pressure_formulation(m_b, m_a, rho_b, rho_a, p_b, p_a,
+                                                  W_b, W_a)
+            @test m_a * acceleration_a + m_b * acceleration_b ≈ zero(W_a) atol = eps()
+
+            acceleration_a = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b,
+                                                  W_a, W_b_asymmetric)
+            acceleration_b = pressure_formulation(m_b, m_a, rho_b, rho_a, p_b, p_a,
+                                                  W_b_asymmetric, W_a)
+            @test m_a * acceleration_a + m_b * acceleration_b ≈ zero(W_a) atol = eps()
+        end
+
+        result32 = @inferred interparticle(1.2f0, 0.8f0, 1000.0f0, 980.0f0,
+                                           2.0f0, 3.0f0, SVector(0.2f0, -0.1f0),
+                                           SVector(-0.13f0, 0.17f0))
+        @test result32 isa SVector{2, Float32}
+        @test tensile_instability_control(m_a, m_b, rho_a, rho_b, 0.0, 0.0,
+                                          W_a) == zero(W_a)
+    end
+
     @testset verbose=true "Illegal Inputs" begin
         correction_dict_1 = Dict(
             "KernelCorrection" => KernelCorrection(),
diff --git a/test/schemes/fluid/rhs.jl b/test/schemes/fluid/rhs.jl
index 20e7e6330a..c0f622d750 100644
--- a/test/schemes/fluid/rhs.jl
+++ b/test/schemes/fluid/rhs.jl
@@ -107,6 +107,39 @@
         end
     end
 
+    @testset "EDAC average-pressure momentum conservation" begin
+        particle_spacing = 0.1
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+        smoothing_length = 1.6particle_spacing
+
+        for correction in (nothing, GradientCorrection())
+            fluid = rectangular_patch(particle_spacing, (4, 3); pressure=1000.0, seed=7)
+            system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel,
+                                                 smoothing_length, sound_speed=10.0,
+                                                 correction,
+                                                 average_pressure_reduction=true)
+            semi = Semidiscretization(system; parallelization_backend=SerialBackend())
+            ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false)
+            v_ode = Array(ode.u0.x[1])
+            u_ode = Array(ode.u0.x[2])
+            dv_ode = zero(v_ode)
+            TrixiParticles.kick!(dv_ode, v_ode, u_ode,
+                                 (; semi=ode.p.semi, split_integration_data=nothing), 0.0)
+
+            system = first(ode.p.semi.systems)
+            dv = TrixiParticles.wrap_v(dv_ode, system, ode.p.semi)
+            acceleration = view(dv, 1:2, :)
+            net_force = vec(sum(system.mass' .* acceleration; dims=2))
+            force_scale = sum(TrixiParticles.eachparticle(system)) do particle
+                norm(system.mass[particle] * acceleration[:, particle])
+            end
+
+            @test force_scale > eps()
+            @test norm(net_force) / force_scale < 5e-13
+            @test all(isfinite, dv_ode)
+        end
+    end
+
     # The following tests for linear and angular momentum and total energy conservation
     # are based on Sections 3.3.4 and 3.4.2 of
     # Daniel J. Price. "Smoothed Particle Hydrodynamics and Magnetohydrodynamics."