diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl
index 00eff1581b..58b171d0f7 100644
--- a/src/general/interpolation.jl
+++ b/src/general/interpolation.jl
@@ -593,7 +593,9 @@ end
                 other_density[point] += m_b * W_ab
 
                 if include_wall_velocity
-                    velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor)
+                    velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor)
+                    velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor,
+                                                         velocity_neighbor_)
                     for i in axes(velocity_neighbor, 1)
                         cache.velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab
                     end
diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl
index 971f3c8a83..768ccc2212 100644
--- a/src/schemes/boundary/open_boundary/rhs.jl
+++ b/src/schemes/boundary/open_boundary/rhs.jl
@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
 
@@ -58,15 +61,17 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_pressure_boundary = 2 * p_boundary * (m_b / (rho_a * rho_b)) * grad_kernel
 
         # Propagate `@inbounds` to the viscosity function, which accesses particle data
-        dv_viscosity_ = @inbounds dv_viscosity(viscosity_model(fluid_system,
-                                                               neighbor_system),
-                                               particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
-
-        dv_particle = dv_pressure + dv_viscosity_ + dv_pressure_boundary
+        dv_viscosity_ = Ref(zero(pos_diff))
+        @inbounds dv_viscosity!(dv_viscosity_,
+                                viscosity_model(fluid_system,
+                                                neighbor_system),
+                                particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
+
+        dv_particle = dv_pressure + dv_viscosity_[] + dv_pressure_boundary
 
         for i in 1:ndims(particle_system)
             @inbounds dv[i, particle] += dv_particle[i]
diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
index 2b205e3cfb..75d3ba6074 100644
--- a/src/schemes/boundary/open_boundary/system.jl
+++ b/src/schemes/boundary/open_boundary/system.jl
@@ -654,8 +654,10 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
             W_ab = kernel(smoothing_kernel, distance, smoothing_length)
             @inbounds shepard_coefficient[point] += volume_b * W_ab
 
+            velocity_neighbor_ = @inbounds current_velocity(v_neighbor, neighbor_system,
+                                                            neighbor)
             velocity_neighbor = @inbounds viscous_velocity(v_neighbor, neighbor_system,
-                                                           neighbor)
+                                                           neighbor, velocity_neighbor_)
             for i in axes(velocity_neighbor, 1)
                 @inbounds sample_velocity[i,
                                           point] += velocity_neighbor[i] * volume_b * W_ab
diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl
index dbfa8613cd..a176be0327 100644
--- a/src/schemes/boundary/wall_boundary/system.jl
+++ b/src/schemes/boundary/wall_boundary/system.jl
@@ -145,15 +145,20 @@ end
     return zero(SVector{ndims(system), eltype(system)})
 end
 
-@inline function viscous_velocity(v, system::WallBoundarySystem, particle)
-    return viscous_velocity(v, system.boundary_model.viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, system::WallBoundarySystem,
+                                              particle, v_particle)
+    return viscous_velocity(v, system.boundary_model.viscosity, system,
+                            particle, v_particle)
 end
 
-@inline function viscous_velocity(v, ::Nothing, system, particle)
-    return current_velocity(v, system, particle)
+@inline function viscous_velocity(v, ::Nothing, system, particle, v_particle)
+    # Regular particle velocity is used for the viscosity calculation by default
+    return v_particle
 end
 
-@inline function viscous_velocity(v, viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, viscosity, system, particle, v_particle)
+    # Wall velocity in the viscosity calculation contains the physical wall velocity
+    # and an interpolated velocity when a wall viscosity (no-slip BC) is used.
     return extract_svector(system.boundary_model.cache.wall_velocity, system, particle)
 end
 
diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl
index a6c98a5d64..5f04efda51 100644
--- a/src/schemes/fluid/entropically_damped_sph/rhs.jl
+++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl
@@ -42,6 +42,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         p_a = @inbounds current_pressure(v_particle_system, particle_system, particle)
         p_b = @inbounds current_pressure(v_neighbor_system, neighbor_system, neighbor)
 
@@ -62,13 +65,12 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             rho_b, pos_diff, distance, grad_kernel,
                                             correction)
 
-        dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
-
-        dv_particle = Ref(dv_pressure + dv_viscosity_)
+        dv_particle = Ref(dv_pressure)
+        @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
 
         # Extra terms in the momentum equation when using a shifting technique
         @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
diff --git a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl
index 8cebf8d20c..eedf81bb15 100644
--- a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl
+++ b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl
@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
         rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
 
@@ -62,14 +65,15 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             distance, grad_kernel, nothing)
 
         # Propagate `@inbounds` to the viscosity function, which accesses particle data
-        dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
+        dv_viscosity_ = Ref(zero(pos_diff))
+        @inbounds dv_viscosity!(dv_viscosity_, particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel)
 
         for i in 1:ndims(particle_system)
-            @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i]
+            @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[][i]
         end
     end
     return dv
diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl
index 8f02018bbf..218dd96f87 100644
--- a/src/schemes/fluid/implicit_incompressible_sph/system.jl
+++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl
@@ -281,16 +281,21 @@ function calculate_predicted_velocity_and_d_ii_values!(system::ImplicitIncompres
             rho_a = @inbounds current_density(v_particle_system, system, particle)
             rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
 
+            v_a = @inbounds current_velocity(v_particle_system, system, particle)
+            v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
             grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
 
-            dv_viscosity_ = @inbounds dv_viscosity(system, neighbor_system,
-                                                   v_particle_system, v_neighbor_system,
-                                                   particle, neighbor, pos_diff, distance,
-                                                   sound_speed, m_a, m_b, rho_a, rho_b,
-                                                   grad_kernel)
+            dv_viscosity_ = Ref(zero(pos_diff))
+            @inbounds dv_viscosity!(dv_viscosity_, system, neighbor_system,
+                                    v_particle_system, v_neighbor_system,
+                                    particle, neighbor, pos_diff, distance,
+                                    sound_speed, m_a, m_b, rho_a, rho_b,
+                                    v_a, v_b, grad_kernel)
             # Add all other non-pressure forces
             for i in 1:ndims(system)
-                @inbounds advection_velocity[i, particle] += time_step * dv_viscosity_[i]
+                @inbounds advection_velocity[i,
+                                             particle] += time_step * dv_viscosity_[][i]
             end
             # Calculate d_ii with eq. 9 in Ihmsen et al. (2013)
             for i in 1:ndims(system)
diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl
index 770e3194be..b39ccb0ef4 100644
--- a/src/schemes/fluid/viscosity.jl
+++ b/src/schemes/fluid/viscosity.jl
@@ -1,33 +1,52 @@
+@propagate_inbounds function viscous_velocity(v, system, particle, v_particle)
+    # Regular particle velocity is used for the viscosity calculation by default
+    return v_particle
+end
+
 # Unpack the neighboring systems viscosity to dispatch on the viscosity type.
 # This function is only necessary to allow `nothing` as viscosity.
 # Otherwise, we could just apply the viscosity as a function directly.
-@propagate_inbounds function dv_viscosity(particle_system, neighbor_system,
-                                          v_particle_system, v_neighbor_system,
-                                          particle, neighbor, pos_diff, distance,
-                                          sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+@propagate_inbounds function dv_viscosity!(dv_particle,
+                                           particle_system::AbstractSystem, neighbor_system,
+                                           v_particle_system, v_neighbor_system,
+                                           particle, neighbor, pos_diff, distance,
+                                           sound_speed, m_a, m_b, rho_a, rho_b,
+                                           v_a, v_b, grad_kernel,
+                                           viscosity_correction=1)
     viscosity = viscosity_model(particle_system, neighbor_system)
 
-    return dv_viscosity(viscosity, particle_system, neighbor_system,
-                        v_particle_system, v_neighbor_system,
-                        particle, neighbor, pos_diff, distance,
-                        sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+    return dv_viscosity!(dv_particle, viscosity,
+                         particle_system, neighbor_system,
+                         v_particle_system, v_neighbor_system,
+                         particle, neighbor, pos_diff, distance,
+                         sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel,
+                         viscosity_correction)
 end
 
-@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system,
-                                          v_particle_system, v_neighbor_system,
-                                          particle, neighbor, pos_diff, distance,
-                                          sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
-    return viscosity(particle_system, neighbor_system,
-                     v_particle_system, v_neighbor_system,
-                     particle, neighbor, pos_diff, distance,
-                     sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+@propagate_inbounds function dv_viscosity!(dv_particle,
+                                           viscosity, particle_system, neighbor_system,
+                                           v_particle_system, v_neighbor_system,
+                                           particle, neighbor, pos_diff, distance,
+                                           sound_speed, m_a, m_b, rho_a, rho_b,
+                                           v_a, v_b, grad_kernel,
+                                           viscosity_correction=1)
+    viscosity(dv_particle, particle_system, neighbor_system,
+              v_particle_system, v_neighbor_system,
+              particle, neighbor, pos_diff, distance,
+              sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel,
+              viscosity_correction)
+
+    return dv_particle
 end
 
-@inline function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system,
-                              v_particle_system, v_neighbor_system,
-                              particle, neighbor, pos_diff, distance,
-                              sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
-    return zero(pos_diff)
+@inline function dv_viscosity!(dv_particle,
+                               viscosity::Nothing,
+                               particle_system, neighbor_system,
+                               v_particle_system, v_neighbor_system,
+                               particle, neighbor, pos_diff, distance,
+                               sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b,
+                               grad_kernel, viscosity_correction=1)
+    return dv_particle
 end
 
 @doc raw"""
@@ -56,6 +75,62 @@ struct ArtificialViscosityMonaghan{ELTYPE}
     end
 end
 
+# See, e.g.,
+# Joseph J. Monaghan. "Smoothed Particle Hydrodynamics".
+# In: Reports on Progress in Physics (2005), pages 1703-1759.
+# [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01)
+@inline function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan,
+                                     smoothing_length, sound_speed)
+    (; alpha) = viscosity
+
+    return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4)
+end
+
+@propagate_inbounds function (viscosity::ArtificialViscosityMonaghan)(dv_particle,
+                                                                      particle_system,
+                                                                      neighbor_system,
+                                                                      v_particle_system,
+                                                                      v_neighbor_system,
+                                                                      particle, neighbor,
+                                                                      pos_diff, distance,
+                                                                      sound_speed,
+                                                                      m_a, m_b,
+                                                                      rho_a, rho_b,
+                                                                      v_a, v_b, grad_kernel,
+                                                                      viscosity_correction=1)
+    v_visc_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_visc_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
+    v_diff = v_visc_a - v_visc_b
+
+    # v_ab ⋅ r_ab
+    vr = dot(v_diff, pos_diff)
+
+    # Monaghan 2005 p. 1741 (doi: 10.1088/0034-4885/68/8/r01):
+    # "In the case of shock tube problems, it is usual to turn the viscosity on for
+    # approaching particles and turn it off for receding particles. In this way, the
+    # viscosity is used for shocks and not rarefactions."
+    if vr < 0
+        h_a = smoothing_length(particle_system, particle)
+        h_b = smoothing_length(neighbor_system, neighbor)
+        h = (h_a + h_b) / 2
+
+        rho_mean = (rho_a + rho_b) / 2
+        (; alpha, beta, epsilon) = viscosity
+
+        # Since this is one of the most performance critical functions, using fast divisions
+        # here gives a significant speedup on GPUs.
+        # See the docs page "Development" for more details on `div_fast`.
+        mu = div_fast(h * vr, distance^2 + epsilon * h^2)
+        c = sound_speed
+        # TODO why is m_b inside the `div_fast` faster on H100 than `m_b * div_fast(...)`?
+        dv_viscosity = div_fast(m_b * alpha * c * mu + m_b * beta * mu^2, rho_mean) *
+                       grad_kernel
+        dv_particle[] += viscosity_correction * dv_viscosity
+    end
+
+    return dv_particle
+end
+
 @doc raw"""
     ViscosityMorris(; nu, epsilon=0.01)
 
@@ -76,30 +151,29 @@ struct ViscosityMorris{ELTYPE}
     end
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityMorris, smoothing_length,
-                             sound_speed)
+@inline function kinematic_viscosity(system, viscosity::ViscosityMorris, smoothing_length,
+                                     sound_speed)
     return viscosity.nu
 end
 
-@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan,
-                                               ViscosityMorris})(particle_system,
-                                                                 neighbor_system,
-                                                                 v_particle_system,
-                                                                 v_neighbor_system,
-                                                                 particle, neighbor,
-                                                                 pos_diff, distance,
-                                                                 sound_speed,
-                                                                 m_a, m_b, rho_a, rho_b,
-                                                                 grad_kernel)
-    rho_mean = (rho_a + rho_b) / 2
-
-    v_a = viscous_velocity(v_particle_system, particle_system, particle)
-    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
-    v_diff = v_a - v_b
+@propagate_inbounds function (viscosity::ViscosityMorris)(dv_particle,
+                                                          particle_system,
+                                                          neighbor_system,
+                                                          v_particle_system,
+                                                          v_neighbor_system,
+                                                          particle, neighbor,
+                                                          pos_diff, distance,
+                                                          sound_speed,
+                                                          m_a, m_b, rho_a, rho_b,
+                                                          v_a, v_b, grad_kernel,
+                                                          viscosity_correction=1)
+    v_visc_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_visc_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
+    v_diff = v_visc_a - v_visc_b
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
-    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
-    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
+    smoothing_length_neighbor = smoothing_length(neighbor_system, neighbor)
+    h = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
                                viscosity_model(neighbor_system, particle_system),
@@ -108,38 +182,6 @@ end
                                viscosity_model(particle_system, neighbor_system),
                                smoothing_length_neighbor, sound_speed)
 
-    pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
-                      smoothing_length_average, grad_kernel, nu_a, nu_b)
-
-    return m_b * pi_ab
-end
-
-@inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance,
-                                                          rho_mean, rho_a, rho_b, h,
-                                                          grad_kernel, nu_a, nu_b)
-    (; alpha, beta, epsilon) = viscosity
-
-    # v_ab ⋅ r_ab
-    vr = dot(v_diff, pos_diff)
-
-    # Monaghan 2005 p. 1741 (doi: 10.1088/0034-4885/68/8/r01):
-    # "In the case of shock tube problems, it is usual to turn the viscosity on for
-    # approaching particles and turn it off for receding particles. In this way, the
-    # viscosity is used for shocks and not rarefactions."
-    if vr < 0
-        # Since this is one of the most performance critical functions, using fast divisions
-        # here gives a significant speedup on GPUs.
-        # See the docs page "Development" for more details on `div_fast`.
-        mu = div_fast(h * vr, distance^2 + epsilon * h^2)
-        return div_fast(alpha * c * mu + beta * mu^2, rho_mean) * grad_kernel
-    end
-
-    return zero(v_diff)
-end
-
-@inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean,
-                                              rho_a, rho_b, h, grad_kernel, nu_a,
-                                              nu_b)
     epsilon = viscosity.epsilon
 
     mu_a = nu_a * rho_a
@@ -148,19 +190,11 @@ end
     # Since this is one of the most performance critical functions, using fast divisions
     # here gives a significant speedup on GPUs.
     # See the docs page "Development" for more details on `div_fast`.
-    return div_fast((mu_a + mu_b) * dot(pos_diff, grad_kernel),
-                    rho_a * rho_b * (distance^2 + epsilon * h^2)) * v_diff
-end
-
-# See, e.g.,
-# Joseph J. Monaghan. "Smoothed Particle Hydrodynamics".
-# In: Reports on Progress in Physics (2005), pages 1703-1759.
-# [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01)
-function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan,
-                             smoothing_length, sound_speed)
-    (; alpha) = viscosity
+    dv_viscosity = div_fast(m_b * (mu_a + mu_b) * dot(pos_diff, grad_kernel),
+                            rho_a * rho_b * (distance^2 + epsilon * h^2)) * v_diff
+    dv_particle[] += viscosity_correction * dv_viscosity
 
-    return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4)
+    return dv_particle
 end
 
 @doc raw"""
@@ -183,8 +217,10 @@ struct ViscosityAdami{ELTYPE}
     end
 end
 
-function adami_viscosity_force(h, pos_diff, distance, grad_kernel,
-                               m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+@inline function adami_viscosity_force!(dv_particle, h, pos_diff, distance,
+                                        grad_kernel, m_a, m_b, rho_a, rho_b,
+                                        v_diff, nu_a, nu_b, epsilon,
+                                        viscosity_correction=1)
     eta_a = nu_a * rho_a
     eta_b = nu_b * rho_b
 
@@ -209,18 +245,22 @@ function adami_viscosity_force(h, pos_diff, distance, grad_kernel,
     # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394
     visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp
 
-    return visc * v_diff
+    dv_particle[] += viscosity_correction * visc * v_diff
+
+    return dv_particle
 end
 
-@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system,
+@inline function (viscosity::ViscosityAdami)(dv_particle, particle_system,
+                                             neighbor_system,
                                              v_particle_system, v_neighbor_system,
                                              particle, neighbor, pos_diff,
                                              distance, sound_speed, m_a, m_b,
-                                             rho_a, rho_b, grad_kernel)
+                                             rho_a, rho_b, v_a, v_b, grad_kernel,
+                                             viscosity_correction=1)
     epsilon = viscosity.epsilon
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
-    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_neighbor = smoothing_length(neighbor_system, neighbor)
     smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
@@ -230,23 +270,20 @@ end
                                viscosity_model(particle_system, neighbor_system),
                                smoothing_length_neighbor, sound_speed)
 
-    v_a = viscous_velocity(v_particle_system, particle_system, particle)
-    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
     v_diff = v_a - v_b
 
-    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
-                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+    return adami_viscosity_force!(dv_particle, smoothing_length_average, pos_diff,
+                                  distance, grad_kernel, m_a, m_b, rho_a, rho_b,
+                                  v_diff, nu_a, nu_b, epsilon, viscosity_correction)
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length,
-                             sound_speed)
+@inline function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length,
+                                     sound_speed)
     return viscosity.nu
 end
 
-@propagate_inbounds function viscous_velocity(v, system, particle)
-    return current_velocity(v, system, particle)
-end
-
 @doc raw"""
     ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01)
 
@@ -291,24 +328,30 @@ This model is appropriate for turbulent flows where unresolved scales contribute
 - `epsilon=0.01`: Parameter to prevent singularities
 """
 struct ViscosityAdamiSGS{ELTYPE}
-    nu      :: ELTYPE      # kinematic viscosity [e.g., 1e-6 m²/s]
-    C_S     :: ELTYPE     # Smagorinsky constant [e.g., 0.1-0.2]
+    nu      :: ELTYPE # Kinematic viscosity [e.g., 1e-6 m²/s]
+    C_S     :: ELTYPE # Smagorinsky constant [e.g., 0.1-0.2]
     epsilon :: ELTYPE # Epsilon for singularity prevention [e.g., 0.001]
 end
 
-ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon)
+function ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001)
+    # The eltype is determined by the type of `nu`.
+    return ViscosityAdamiSGS{typeof(nu)}(nu, C_S, epsilon)
+end
 
-@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system,
+@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(dv_particle,
+                                                            particle_system,
                                                             neighbor_system,
                                                             v_particle_system,
                                                             v_neighbor_system,
                                                             particle, neighbor, pos_diff,
                                                             distance, sound_speed, m_a, m_b,
-                                                            rho_a, rho_b, grad_kernel)
+                                                            rho_a, rho_b,
+                                                            v_a, v_b, grad_kernel,
+                                                            viscosity_correction=1)
     epsilon = viscosity.epsilon
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
-    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_neighbor = smoothing_length(neighbor_system, neighbor)
     smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
@@ -318,8 +361,8 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps
                                viscosity_model(particle_system, neighbor_system),
                                smoothing_length_neighbor, sound_speed)
 
-    v_a = viscous_velocity(v_particle_system, particle_system, particle)
-    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
     v_diff = v_a - v_b
 
     # ------------------------------------------------------------------------------
@@ -353,12 +396,13 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps
     nu_a = nu_a + nu_SGS
     nu_b = nu_b + nu_SGS
 
-    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
-                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+    return adami_viscosity_force!(dv_particle, smoothing_length_average, pos_diff,
+                                  distance, grad_kernel, m_a, m_b, rho_a, rho_b,
+                                  v_diff, nu_a, nu_b, epsilon, viscosity_correction)
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length,
-                             sound_speed)
+@inline function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS,
+                                     smoothing_length, sound_speed)
     return viscosity.nu
 end
 
@@ -406,24 +450,29 @@ This model is appropriate for turbulent flows where unresolved scales contribute
 - `epsilon=0.01`: Parameter to prevent singularities
 """
 struct ViscosityMorrisSGS{ELTYPE}
-    nu::ELTYPE      # kinematic viscosity [e.g., 1e-6 m²/s]
-    C_S::ELTYPE     # Smagorinsky constant [e.g., 0.1-0.2]
-    epsilon::ELTYPE # Epsilon for singularity prevention [e.g., 0.001]
+    nu      :: ELTYPE # Kinematic viscosity [e.g., 1e-6 m²/s]
+    C_S     :: ELTYPE # Smagorinsky constant [e.g., 0.1-0.2]
+    epsilon :: ELTYPE # Epsilon for singularity prevention [e.g., 0.001]
 end
 
-ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon)
+function ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001)
+    return ViscosityMorrisSGS{typeof(nu)}(nu, C_S, epsilon)
+end
 
-@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system,
+@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(dv_particle,
+                                                             particle_system,
                                                              neighbor_system,
                                                              v_particle_system,
                                                              v_neighbor_system,
                                                              particle, neighbor, pos_diff,
                                                              distance, sound_speed, m_a,
-                                                             m_b, rho_a, rho_b, grad_kernel)
+                                                             m_b, rho_a, rho_b,
+                                                             v_a, v_b, grad_kernel,
+                                                             viscosity_correction=1)
     epsilon = viscosity.epsilon
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
-    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_neighbor = smoothing_length(neighbor_system, neighbor)
     h = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
     nu_a = kinematic_viscosity(particle_system,
@@ -433,8 +482,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e
                                viscosity_model(particle_system, neighbor_system),
                                smoothing_length_neighbor, sound_speed)
 
-    v_a = viscous_velocity(v_particle_system, particle_system, particle)
-    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
     v_diff = v_a - v_b
 
     # SGS part: Compute the subgrid-scale eddy viscosity.
@@ -456,12 +505,15 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e
     # Since this is one of the most performance critical functions, using fast divisions
     # here gives a significant speedup on GPUs.
     # See the docs page "Development" for more details on `div_fast`.
-    return div_fast(m_b * (mu_a + mu_b) * dot(pos_diff, grad_kernel),
-                    rho_a * rho_b * (distance^2 + epsilon * h^2)) * v_diff
+    dv_viscosity = div_fast(m_b * (mu_a + mu_b) * dot(pos_diff, grad_kernel),
+                            rho_a * rho_b * (distance^2 + epsilon * h^2)) * v_diff
+    dv_particle[] += viscosity_correction * dv_viscosity
+
+    return dv_particle
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length,
-                             sound_speed)
+@inline function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS,
+                                     smoothing_length, sound_speed)
     return viscosity.nu
 end
 
@@ -490,27 +542,28 @@ struct ViscosityCarreauYasuda{ELTYPE}
 end
 
 function ViscosityCarreauYasuda(; nu0, nu_inf, lambda, a, n, epsilon=0.01)
-    ViscosityCarreauYasuda(nu0, nu_inf, lambda, a, n, epsilon)
+    ViscosityCarreauYasuda{typeof(nu0)}(nu0, nu_inf, lambda, a, n, epsilon)
 end
 
-@propagate_inbounds function (viscosity::ViscosityCarreauYasuda)(particle_system,
+@propagate_inbounds function (viscosity::ViscosityCarreauYasuda)(dv_particle,
+                                                                 particle_system,
                                                                  neighbor_system,
                                                                  v_particle_system,
                                                                  v_neighbor_system,
                                                                  particle, neighbor,
                                                                  pos_diff, distance,
                                                                  sound_speed,
-                                                                 m_a, m_b,
-                                                                 rho_a, rho_b,
-                                                                 grad_kernel)
+                                                                 m_a, m_b, rho_a, rho_b,
+                                                                 v_a, v_b, grad_kernel,
+                                                                 viscosity_correction=1)
     epsilon = viscosity.epsilon
 
     smoothing_length_particle = smoothing_length(particle_system, particle)
-    smoothing_length_neighbor = smoothing_length(particle_system, neighbor)
+    smoothing_length_neighbor = smoothing_length(neighbor_system, neighbor)
     smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
 
-    v_a = viscous_velocity(v_particle_system, particle_system, particle)
-    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor)
+    v_a = viscous_velocity(v_particle_system, particle_system, particle, v_a)
+    v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor, v_b)
     v_diff = v_a - v_b
 
     # Since this is one of the most performance critical functions, using fast divisions
@@ -524,11 +577,12 @@ end
     nu_a = nu_eff
     nu_b = nu_eff
 
-    return adami_viscosity_force(smoothing_length_average, pos_diff, distance, grad_kernel,
-                                 m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon)
+    return adami_viscosity_force!(dv_particle, smoothing_length_average, pos_diff,
+                                  distance, grad_kernel, m_a, m_b, rho_a, rho_b,
+                                  v_diff, nu_a, nu_b, epsilon, viscosity_correction)
 end
 
-function kinematic_viscosity(system, viscosity::ViscosityCarreauYasuda,
-                             smoothing_length, sound_speed)
+@inline function kinematic_viscosity(system, viscosity::ViscosityCarreauYasuda,
+                                     smoothing_length, sound_speed)
     return viscosity.nu0
 end
diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
index 96915d8b14..4a0350fa90 100644
--- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl
+++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl
@@ -60,6 +60,9 @@ function interact!(dv, v_particle_system, u_particle_system,
         m_a = @inbounds hydrodynamic_mass(particle_system, particle)
         m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
 
+        v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
+        v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         # The following call is equivalent to
         #     `p_a = current_pressure(v_particle_system, particle_system, particle)`
         #     `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
@@ -78,15 +81,14 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
                                             distance, grad_kernel, correction)
 
+        dv_particle = Ref(dv_pressure)
+
         # Propagate `@inbounds` to the viscosity function, which accesses particle data
-        dv_viscosity_ = viscosity_correction *
-                        @inbounds dv_viscosity(particle_system, neighbor_system,
-                                               v_particle_system, v_neighbor_system,
-                                               particle, neighbor, pos_diff, distance,
-                                               sound_speed, m_a, m_b, rho_a, rho_b,
-                                               grad_kernel)
-
-        dv_particle = Ref(dv_pressure + dv_viscosity_)
+        @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
+                                v_particle_system, v_neighbor_system,
+                                particle, neighbor, pos_diff, distance,
+                                sound_speed, m_a, m_b, rho_a, rho_b,
+                                v_a, v_b, grad_kernel, viscosity_correction)
 
         # Extra terms in the momentum equation when using a shifting technique
         @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl
index 20a726b7b1..0ba34b1fce 100644
--- a/src/schemes/structure/rigid_body/system.jl
+++ b/src/schemes/structure/rigid_body/system.jl
@@ -509,9 +509,12 @@ end
     return neighbor_system.boundary_model.viscosity
 end
 
-@inline function viscous_velocity(v, system::RigidBodySystem, particle)
-    # This function is only used in fluid-structure interaction, so it is never called when `boundary_model` is `nothing`
-    return viscous_velocity(v, system.boundary_model.viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, system::RigidBodySystem,
+                                              particle, v_particle)
+    # This function is only used in fluid-structure interaction,
+    # so it is never called when `boundary_model` is `nothing`
+    return viscous_velocity(v, system.boundary_model.viscosity, system,
+                            particle, v_particle)
 end
 
 @inline acceleration_source(system::RigidBodySystem) = system.acceleration
diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl
index ff872c04fc..0c508c2e71 100644
--- a/src/schemes/structure/structure.jl
+++ b/src/schemes/structure/structure.jl
@@ -58,6 +58,9 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
         rho_a = current_density(v_particle_system, particle_system, particle)
         rho_b = current_density(v_neighbor_system, neighbor_system, neighbor)
 
+        v_a = current_velocity(v_particle_system, particle_system, particle)
+        v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor)
+
         surface_tension = surface_tension_model(neighbor_system)
 
         # In fluid-structure interaction, use the "hydrodynamic mass" of the structure particles
@@ -78,12 +81,13 @@ function interact_structure_fluid!(dv, v_particle_system, u_particle_system,
                                             pos_diff, distance, grad_kernel,
                                             system_correction(neighbor_system))
 
-        dv_viscosity_ = dv_viscosity(neighbor_system, particle_system,
-                                     v_neighbor_system, v_particle_system,
-                                     neighbor, particle, pos_diff, distance,
-                                     sound_speed, m_b, m_a, rho_b, rho_a, grad_kernel)
+        dv_particle = Ref(dv_boundary)
+        dv_viscosity!(dv_particle, neighbor_system, particle_system,
+                      v_neighbor_system, v_particle_system,
+                      neighbor, particle, pos_diff, distance,
+                      sound_speed, m_b, m_a, rho_b, rho_a,
+                      v_b, v_a, grad_kernel)
 
-        dv_particle = Ref(dv_boundary + dv_viscosity_)
         adhesion_force!(dv_particle, surface_tension, neighbor_system, particle_system,
                         neighbor, particle, pos_diff, distance)
 
diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl
index 3383929c28..874f1c92f7 100644
--- a/src/schemes/structure/total_lagrangian_sph/rhs.jl
+++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl
@@ -67,14 +67,13 @@ end
                                                        current_pos_diff, current_distance,
                                                        system, m_a, m_b, rho_a, rho_b)
 
-        dv_viscosity = @inbounds dv_viscosity_tlsph(system, v_system, particle, neighbor,
-                                                    current_pos_diff, current_distance,
-                                                    m_a, m_b, rho_a, rho_b, grad_kernel)
-
-        dv_particle = dv_stress + dv_penalty_force_ + dv_viscosity
+        dv_particle = Ref(dv_stress + dv_penalty_force_)
+        @inbounds dv_viscosity_tlsph!(dv_particle, system, v_system, particle, neighbor,
+                                      current_pos_diff, current_distance,
+                                      m_a, m_b, rho_a, rho_b, grad_kernel)
 
         for i in 1:ndims(system)
-            @inbounds dv[i, particle] += dv_particle[i]
+            @inbounds dv[i, particle] += dv_particle[][i]
         end
 
         # TODO continuity equation for boundary model with `ContinuityDensity`?
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index b4fb02f31a..81237caf14 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -333,9 +333,12 @@ end
     return current_density(v, system.boundary_model, system)
 end
 
-@inline function viscous_velocity(v, system::TotalLagrangianSPHSystem, particle)
-    # This function is only used in fluid-structure interaction, so it is never called when `boundary_model` is `nothing`
-    return viscous_velocity(v, system.boundary_model.viscosity, system, particle)
+@propagate_inbounds function viscous_velocity(v, system::TotalLagrangianSPHSystem,
+                                              particle, v_particle)
+    # This function is only used in fluid-structure interaction,
+    # so it is never called when `boundary_model` is `nothing`
+    return viscous_velocity(v, system.boundary_model.viscosity, system,
+                            particle, v_particle)
 end
 
 # In fluid-structure interaction, use the "hydrodynamic pressure" of the structure particles
diff --git a/src/schemes/structure/total_lagrangian_sph/viscosity.jl b/src/schemes/structure/total_lagrangian_sph/viscosity.jl
index 50629e8966..b19aec0117 100644
--- a/src/schemes/structure/total_lagrangian_sph/viscosity.jl
+++ b/src/schemes/structure/total_lagrangian_sph/viscosity.jl
@@ -1,73 +1,83 @@
 # Unpack the neighboring systems viscosity to dispatch on the viscosity type.
 # This function is only necessary to allow `nothing` as viscosity.
 # Otherwise, we could just apply the viscosity as a function directly.
-@propagate_inbounds function dv_viscosity_tlsph(system, v_system, particle, neighbor,
-                                                current_pos_diff, current_distance,
-                                                m_a, m_b, rho_a, rho_b, grad_kernel)
+@propagate_inbounds function dv_viscosity_tlsph!(dv_particle, system, v_system,
+                                                 particle, neighbor,
+                                                 current_pos_diff, current_distance,
+                                                 m_a, m_b, rho_a, rho_b, grad_kernel)
     viscosity = system.viscosity
 
-    return dv_viscosity_tlsph(viscosity, system, v_system, particle, neighbor,
-                              current_pos_diff, current_distance,
-                              m_a, m_b, rho_a, rho_b, grad_kernel)
+    return dv_viscosity_tlsph!(dv_particle, viscosity, system, v_system,
+                               particle, neighbor, current_pos_diff, current_distance,
+                               m_a, m_b, rho_a, rho_b, grad_kernel)
 end
 
-@propagate_inbounds function dv_viscosity_tlsph(viscosity, system,
-                                                v_system, particle, neighbor,
-                                                current_pos_diff, current_distance,
-                                                m_a, m_b, rho_a, rho_b, grad_kernel)
-    return viscosity(system, v_system, particle, neighbor,
+@propagate_inbounds function dv_viscosity_tlsph!(dv_particle, viscosity, system,
+                                                 v_system, particle, neighbor,
+                                                 current_pos_diff, current_distance,
+                                                 m_a, m_b, rho_a, rho_b, grad_kernel)
+    return viscosity(dv_particle, system, v_system, particle, neighbor,
                      current_pos_diff, current_distance,
                      m_a, m_b, rho_a, rho_b, grad_kernel)
 end
 
-@inline function dv_viscosity_tlsph(viscosity::Nothing, system,
-                                    v_system, particle, neighbor,
-                                    current_pos_diff, current_distance,
-                                    m_a, m_b, rho_a, rho_b, grad_kernel)
+@inline function dv_viscosity_tlsph!(dv_particle, viscosity::Nothing, system,
+                                     v_system, particle, neighbor,
+                                     current_pos_diff, current_distance,
+                                     m_a, m_b, rho_a, rho_b, grad_kernel)
     return zero(current_pos_diff)
 end
 
 # Applying the viscosity according to Lin et al. (2015):
 # "Geometrically nonlinear analysis of two-dimensional structures using an improved
 # smoothed particle hydrodynamics method"
-@propagate_inbounds function (viscosity::ArtificialViscosityMonaghan)(system::TotalLagrangianSPHSystem,
+@propagate_inbounds function (viscosity::ArtificialViscosityMonaghan)(dv_particle,
+                                                                      system::TotalLagrangianSPHSystem,
                                                                       v_system,
                                                                       particle, neighbor,
                                                                       current_pos_diff,
                                                                       current_distance,
                                                                       m_a, m_b, rho_a,
                                                                       rho_b, grad_kernel)
-    rho_mean = (rho_a + rho_b) / 2
-
     v_a = current_velocity(v_system, system, particle)
     v_b = current_velocity(v_system, system, neighbor)
     v_diff = v_a - v_b
 
-    smoothing_length_particle = smoothing_length(system, particle)
-    smoothing_length_neighbor = smoothing_length(system, neighbor)
-    smoothing_length_average = (smoothing_length_particle + smoothing_length_neighbor) / 2
+    # v_ab ⋅ r_ab
+    vr = dot(v_diff, current_pos_diff)
 
-    # Compute bulk modulus from Young's modulus and Poisson's ratio.
-    # See the table at the end of https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
-    E = young_modulus(system, particle)
-    K = E / (ndims(system) * (1 - 2 * poisson_ratio(system, particle)))
+    # Monaghan 2005 p. 1741 (doi: 10.1088/0034-4885/68/8/r01):
+    # "In the case of shock tube problems, it is usual to turn the viscosity on for
+    # approaching particles and turn it off for receding particles. In this way, the
+    # viscosity is used for shocks and not rarefactions."
+    if vr < 0
+        # Compute bulk modulus from Young's modulus and Poisson's ratio.
+        # See the table at the end of https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
+        E = young_modulus(system, particle)
+        K = E / (ndims(system) * (1 - 2 * poisson_ratio(system, particle)))
 
-    # Newton–Laplace equation
-    sound_speed = sqrt(K / rho_a)
+        # Newton–Laplace equation
+        sound_speed = sqrt(K / rho_a)
 
-    # This is not needed for `ArtificialViscosityMonaghan`
-    nu_a = nu_b = 0
+        h_a = smoothing_length(system, particle)
+        h_b = smoothing_length(system, neighbor)
+        h = (h_a + h_b) / 2
 
-    pi_ab = viscosity(sound_speed, v_diff, current_pos_diff, current_distance,
-                      rho_mean, rho_a, rho_b, smoothing_length_average,
-                      grad_kernel, nu_a, nu_b)
+        rho_mean = (rho_a + rho_b) / 2
 
-    # See eq. 26 of Lin et al. (2015)
-    F = deformation_gradient(system, particle)
+        (; alpha, beta, epsilon) = viscosity
+        mu = h * vr / (current_distance^2 + epsilon * h^2)
+        c = sound_speed
+        pi_ab = (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel
 
-    if abs(det(F)) < 1.0f-9
-        return zero(grad_kernel)
+        F = deformation_gradient(system, particle)
+        det_F = det(F)
+        if abs(det_F) < 1.0f-9
+            return dv_particle
+        end
+        # See eq. 26 of Lin et al. (2015)
+        dv_particle[] += m_b * det_F * inv(F)' * pi_ab
     end
 
-    return m_b * det(F) * inv(F)' * pi_ab
+    return dv_particle
 end
diff --git a/test/general/density_calculator.jl b/test/general/density_calculator.jl
index 2ded394766..6bf895f7be 100644
--- a/test/general/density_calculator.jl
+++ b/test/general/density_calculator.jl
@@ -1,5 +1,3 @@
-using OrdinaryDiffEq
-
 # Setup a single particle and calculate its density
 @testset verbose=true "DensityCalculators" begin
     @testset verbose=true "SummationDensity" begin
diff --git a/test/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl
index 0e359ef5f1..d4ead2cbcc 100644
--- a/test/schemes/fluid/viscosity.jl
+++ b/test/schemes/fluid/viscosity.jl
@@ -25,12 +25,19 @@
         grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
                                                            distance, 1)
 
-        dv = viscosity(sound_speed, v_diff, pos_diff, distance,
-                       rho_mean, rho_a, rho_b, smoothing_length,
-                       grad_kernel, 0.0, 0.0)
-
-        @test isapprox(dv[1], -0.02049217623299368, atol=6e-15)
-        @test isapprox(dv[2], 0.03073826434949052, atol=6e-15)
+        v = nothing
+        v_a = v_diff
+        v_b = zeros(2)
+        m_a = 1.0
+        m_b = 1.0
+
+        dv = Ref(zero(v_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+
+        @test isapprox(dv[][1], -0.02049217623299368, atol=6e-15)
+        @test isapprox(dv[][2], 0.03073826434949052, atol=6e-15)
     end
     @testset verbose=true "`ViscosityMorris`" begin
         nu = 7e-3
@@ -51,12 +58,16 @@
         v[1, 2] = 0.0
         v[2, 2] = 0.0
 
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
+
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
-        @test isapprox(dv[1], -1.0895602048035404e-5, atol=6e-15)
-        @test isapprox(dv[2], 3.631867349345135e-5, atol=6e-15)
+        @test isapprox(dv[][1], -1.0895602048035404e-5, atol=6e-15)
+        @test isapprox(dv[][2], 3.631867349345135e-5, atol=6e-15)
     end
     @testset verbose=true "`ViscosityAdami`" begin
         viscosity = ViscosityAdami(nu=7e-3)
@@ -76,12 +87,16 @@
         v[1, 2] = 0.0
         v[2, 2] = 0.0
 
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
 
-        @test isapprox(dv[1], -1.089560204803541e-5, atol=6e-15)
-        @test isapprox(dv[2], 3.6318673493451364e-5, atol=6e-15)
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+
+        @test isapprox(dv[][1], -1.089560204803541e-5, atol=6e-15)
+        @test isapprox(dv[][2], 3.6318673493451364e-5, atol=6e-15)
     end
     @testset verbose=true "`ViscosityMorrisSGS`" begin
         nu = 7e-3
@@ -103,12 +118,16 @@
         v[1, 2] = 0.0
         v[2, 2] = 0.0
 
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
+
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
-        @test isapprox(dv[1], -2.032835697804103e-5, atol=6e-15)
-        @test isapprox(dv[2], 6.776118992680343e-5, atol=6e-15)
+        @test isapprox(dv[][1], -2.032835697804103e-5, atol=6e-15)
+        @test isapprox(dv[][2], 6.776118992680343e-5, atol=6e-15)
     end
     @testset verbose=true "`ViscosityAdamiSGS`" begin
         viscosity = ViscosityAdamiSGS(nu=7e-3)
@@ -128,12 +147,16 @@
         v[1, 2] = 0.0
         v[2, 2] = 0.0
 
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
 
-        @test isapprox(dv[1], -2.0328356978041036e-5, atol=6e-15)
-        @test isapprox(dv[2], 6.776118992680346e-5, atol=6e-15)
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
+
+        @test isapprox(dv[][1], -2.0328356978041036e-5, atol=6e-15)
+        @test isapprox(dv[][2], 6.776118992680346e-5, atol=6e-15)
     end
 
     @testset verbose=true "`ViscosityCarreauYasuda`" begin
@@ -161,12 +184,16 @@
 
         grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
                                                            distance, 1)
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
+
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
-        @test isapprox(dv[1], -1.0895602048035410e-5, atol=6e-15)
-        @test isapprox(dv[2], 3.6318673493451364e-5, atol=6e-15)
+        @test isapprox(dv[][1], -1.0895602048035410e-5, atol=6e-15)
+        @test isapprox(dv[][2], 3.6318673493451364e-5, atol=6e-15)
 
         # ------------------------------------------------------------------
         # 2) Shear-thinning case: fixed (precomputed) values
@@ -183,11 +210,15 @@
 
         grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff,
                                                            distance, 1)
-        dv = viscosity(system_wcsph, system_wcsph,
-                       v, v, 1, 2, pos_diff, distance,
-                       sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel)
+        v_a = v[:, 1]
+        v_b = v[:, 2]
+
+        dv = Ref(zero(pos_diff))
+        viscosity(dv, system_wcsph, system_wcsph,
+                  v, v, 1, 2, pos_diff, distance,
+                  sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
 
-        @test isapprox(dv[1], -5.33743497379846e-9, atol=6e-15)
-        @test isapprox(dv[2], 1.7791449912661534e-8, atol=6e-15)
+        @test isapprox(dv[][1], -5.33743497379846e-9, atol=6e-15)
+        @test isapprox(dv[][2], 1.7791449912661534e-8, atol=6e-15)
     end
 end