From b454780d4628640eeb4a4f4be126715092e54b33 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 23 Mar 2026 13:28:55 +0100 Subject: [PATCH 01/12] Avoid duplicate velocity loads for viscosity calculation --- src/general/interpolation.jl | 4 +- src/schemes/boundary/open_boundary/rhs.jl | 5 +- src/schemes/boundary/open_boundary/system.jl | 4 +- src/schemes/boundary/wall_boundary/system.jl | 14 ++-- .../fluid/entropically_damped_sph/rhs.jl | 5 +- .../fluid/implicit_incompressible_sph/rhs.jl | 5 +- .../implicit_incompressible_sph/system.jl | 5 +- src/schemes/fluid/viscosity.jl | 79 ++++++++++--------- .../fluid/weakly_compressible_sph/rhs.jl | 5 +- src/schemes/structure/rigid_body/system.jl | 8 +- src/schemes/structure/structure.jl | 6 +- .../structure/total_lagrangian_sph/system.jl | 8 +- 12 files changed, 92 insertions(+), 56 deletions(-) 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 186a9e7b7d..ac104925ed 100644 --- a/src/schemes/boundary/open_boundary/rhs.jl +++ b/src/schemes/boundary/open_boundary/rhs.jl @@ -23,6 +23,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) + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) m_a = @inbounds hydrodynamic_mass(particle_system, particle) @@ -48,7 +51,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) + v_a, v_b, grad_kernel) dv_particle = dv_pressure + dv_viscosity_ + dv_pressure_boundary diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 94aee1b048..dc5b58f7b7 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -644,8 +644,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..6fc3f48d2c 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -145,15 +145,19 @@ 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) +@inline 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) +@inline 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 cea18484d4..3a93816e40 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -24,6 +24,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) @@ -50,7 +53,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) + v_a, v_b, grad_kernel) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(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 ca7c4b2c35..f42ca7f18c 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl @@ -23,6 +23,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) + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) m_a = @inbounds hydrodynamic_mass(particle_system, particle) @@ -50,7 +53,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) + v_a, v_b, grad_kernel) for i in 1:ndims(particle_system) @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index 8f02018bbf..b685911ff6 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -281,13 +281,16 @@ 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) + 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] diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 2a6c51ba71..a8ed4d2014 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -1,32 +1,39 @@ +@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) + sound_speed, m_a, m_b, rho_a, rho_b, + v_a, v_b, grad_kernel) 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) + sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel) 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) + sound_speed, m_a, m_b, rho_a, rho_b, + v_a, v_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) + sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel) 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) + sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel) return zero(pos_diff) end @@ -76,8 +83,8 @@ 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 @@ -90,12 +97,12 @@ end pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) + v_a, v_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 + 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) @@ -150,8 +157,8 @@ end # 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) +@inline function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan, + smoothing_length, sound_speed) (; alpha) = viscosity return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4) @@ -177,8 +184,9 @@ struct ViscosityAdami{ELTYPE} end end -function 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) +@inline function 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) eta_a = nu_a * rho_a eta_b = nu_b * rho_b @@ -207,7 +215,7 @@ end 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) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -221,23 +229,19 @@ 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) 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) @@ -295,7 +299,8 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps 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) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -309,8 +314,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 # ------------------------------------------------------------------------------ @@ -407,7 +412,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e 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) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -421,8 +427,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. @@ -483,17 +489,16 @@ end 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) epsilon = viscosity.epsilon 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 - 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 gamma_dot = norm(v_diff) / (distance + epsilon) @@ -508,7 +513,7 @@ end m_a, m_b, rho_a, rho_b, v_diff, nu_a, nu_b, epsilon) 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 acb0f06894..55b9eb93e9 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -44,6 +44,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)` @@ -68,7 +71,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - grad_kernel) + v_a, v_b, grad_kernel) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index c9281aa5f9..82279cf949 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -506,9 +506,11 @@ 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) +@inline 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 1ac304d79f..649ef5f30a 100644 --- a/src/schemes/structure/structure.jl +++ b/src/schemes/structure/structure.jl @@ -46,6 +46,9 @@ function interact_structure_fluid!(dv, v_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 @@ -69,7 +72,8 @@ function interact_structure_fluid!(dv, v_particle_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) + sound_speed, m_b, m_a, rho_b, rho_a, + v_a, v_b, grad_kernel) dv_adhesion = adhesion_force(surface_tension, neighbor_system, particle_system, neighbor, particle, pos_diff, distance) diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index cac4ed1f12..37b86ed85c 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -333,9 +333,11 @@ 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) +@inline 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 From f79a933e47bac932882ebdfb0e6b2470b153a690 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 1 Apr 2026 11:30:27 +0200 Subject: [PATCH 02/12] Make viscosity inplace --- src/schemes/boundary/open_boundary/rhs.jl | 20 ++-- .../fluid/entropically_damped_sph/rhs.jl | 13 +-- .../fluid/implicit_incompressible_sph/rhs.jl | 13 +-- .../implicit_incompressible_sph/system.jl | 14 +-- src/schemes/fluid/viscosity.jl | 89 ++++++++++------- .../fluid/weakly_compressible_sph/rhs.jl | 15 +-- src/schemes/structure/structure.jl | 13 +-- test/general/density_calculator.jl | 2 - test/schemes/fluid/viscosity.jl | 95 ++++++++++++------- 9 files changed, 161 insertions(+), 113 deletions(-) diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl index ac104925ed..4fd09f5f87 100644 --- a/src/schemes/boundary/open_boundary/rhs.jl +++ b/src/schemes/boundary/open_boundary/rhs.jl @@ -45,15 +45,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, - v_a, v_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/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3a93816e40..9022c23218 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -49,11 +49,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, - v_a, v_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) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), @@ -70,7 +71,7 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension + + dv_particle = dv_pressure + dv_viscosity_[] + dv_tvf + dv_surface_tension + dv_adhesion for i in 1:ndims(particle_system) diff --git a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl index f42ca7f18c..25b3d09097 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl @@ -49,14 +49,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, - v_a, v_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 b685911ff6..a07b52786b 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -286,14 +286,16 @@ function calculate_predicted_velocity_and_d_ii_values!(system::ImplicitIncompres 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, - v_a, v_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 a8ed4d2014..5046814293 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -6,35 +6,38 @@ 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, +@propagate_inbounds function 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 = viscosity_model(particle_system, neighbor_system) - return dv_viscosity(viscosity, particle_system, neighbor_system, + 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) end -@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system, +@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) - return viscosity(particle_system, neighbor_system, + return 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) end -@inline function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system, +@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) - return zero(pos_diff) + return dv_particle end @doc raw""" @@ -89,7 +92,8 @@ end end @propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, - ViscosityMorris})(particle_system, + ViscosityMorris})(dv_particle, + particle_system, neighbor_system, v_particle_system, v_neighbor_system, @@ -115,15 +119,16 @@ 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) + viscosity(dv_particle, sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, + smoothing_length_average, grad_kernel, nu_a, nu_b, m_b) - return m_b * pi_ab + return dv_particle end -@inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance, +@inline function (viscosity::ArtificialViscosityMonaghan)(dv_particle, c, v_diff, + pos_diff, distance, rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b) + grad_kernel, nu_a, nu_b, m_b) (; alpha, beta, epsilon) = viscosity # v_ab ⋅ r_ab @@ -135,22 +140,24 @@ end # viscosity is used for shocks and not rarefactions." if vr < 0 mu = h * vr / (distance^2 + epsilon * h^2) - return (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel + dv_particle[] += m_b * (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel end - return zero(v_diff) + return dv_particle end -@inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean, - rho_a, rho_b, h, grad_kernel, nu_a, - nu_b) +@inline function (viscosity::ViscosityMorris)(dv_particle, c, v_diff, pos_diff, + distance, rho_mean, rho_a, rho_b, h, + grad_kernel, nu_a, nu_b, m_b) epsilon = viscosity.epsilon mu_a = nu_a * rho_a mu_b = nu_b * rho_b - return (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) / - (distance^2 + epsilon * h^2) * v_diff + dv_particle[] += m_b * (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) / + (distance^2 + epsilon * h^2) * v_diff + + return dv_particle end # See, e.g., @@ -184,9 +191,9 @@ struct ViscosityAdami{ELTYPE} end end -@inline function 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) +@inline function 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) eta_a = nu_a * rho_a eta_b = nu_b * rho_b @@ -208,10 +215,13 @@ end # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a - return visc .* v_diff + dv_particle[] += 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, @@ -233,8 +243,9 @@ end 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) end @inline function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, @@ -293,7 +304,8 @@ end ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, epsilon) -@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(particle_system, +@propagate_inbounds function (viscosity::ViscosityAdamiSGS)(dv_particle, + particle_system, neighbor_system, v_particle_system, v_neighbor_system, @@ -346,8 +358,9 @@ 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) end function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, @@ -406,7 +419,8 @@ end ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, epsilon) -@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(particle_system, +@propagate_inbounds function (viscosity::ViscosityMorrisSGS)(dv_particle, + particle_system, neighbor_system, v_particle_system, v_neighbor_system, @@ -444,9 +458,10 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e mu_a = nu_a_eff * rho_a mu_b = nu_b_eff * rho_b - force_Morris = (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / - (distance^2 + epsilon * smoothing_length_average^2) * v_diff - return m_b * force_Morris + dv_particle[] += m_b * (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / + (distance^2 + epsilon * smoothing_length_average^2) * v_diff + + return dv_particle end function kinematic_viscosity(system, viscosity::ViscosityMorrisSGS, smoothing_length, @@ -482,7 +497,8 @@ function ViscosityCarreauYasuda(; nu0, nu_inf, lambda, a, n, epsilon=0.01) ViscosityCarreauYasuda(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, @@ -509,8 +525,9 @@ 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) end @inline function kinematic_viscosity(system, viscosity::ViscosityCarreauYasuda, diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 55b9eb93e9..000d08ff57 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -66,12 +66,12 @@ function interact!(dv, v_particle_system, u_particle_system, distance, grad_kernel, correction) # 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, - v_a, v_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) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), @@ -89,7 +89,8 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - dv_particle = dv_pressure + dv_viscosity_ + dv_tvf + dv_surface_tension + + dv_particle = dv_pressure + viscosity_correction * dv_viscosity_[] + dv_tvf + + dv_surface_tension + dv_adhesion for i in 1:ndims(particle_system) diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl index 649ef5f30a..a3b216657d 100644 --- a/src/schemes/structure/structure.jl +++ b/src/schemes/structure/structure.jl @@ -69,16 +69,17 @@ function interact_structure_fluid!(dv, v_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, - v_a, v_b, grad_kernel) + dv_viscosity_ = Ref(zero(pos_diff)) + 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, + v_a, v_b, grad_kernel) dv_adhesion = adhesion_force(surface_tension, neighbor_system, particle_system, neighbor, particle, pos_diff, distance) - dv_fs = dv_boundary + dv_viscosity_ + dv_adhesion + dv_fs = dv_boundary + dv_viscosity_[] + dv_adhesion accumulate_structure_fluid_pair!(dv, dv_fs, particle_system, particle, m_b) 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..ca50e6d5c1 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -25,12 +25,13 @@ 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) + dv = Ref(zero(v_diff)) + viscosity(dv, sound_speed, v_diff, pos_diff, distance, + rho_mean, rho_a, rho_b, smoothing_length, + grad_kernel, 0.0, 0.0, 1.0) - @test isapprox(dv[1], -0.02049217623299368, atol=6e-15) - @test isapprox(dv[2], 0.03073826434949052, atol=6e-15) + @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 +52,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.0895602048035404e-5, atol=6e-15) - @test isapprox(dv[2], 3.631867349345135e-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.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 +81,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.089560204803541e-5, atol=6e-15) - @test isapprox(dv[2], 3.6318673493451364e-5, atol=6e-15) + @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 +112,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 +141,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 +178,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 +204,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 From 327b47c0616ea9e62883786317ff446b9c073ada Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 2 Apr 2026 17:41:10 +0200 Subject: [PATCH 03/12] Rename `dv_viscosity` and add free surface correction --- src/schemes/boundary/open_boundary/rhs.jl | 16 +-- .../fluid/entropically_damped_sph/rhs.jl | 10 +- .../fluid/implicit_incompressible_sph/rhs.jl | 10 +- .../implicit_incompressible_sph/system.jl | 10 +- src/schemes/fluid/viscosity.jl | 97 +++++++++++-------- .../fluid/weakly_compressible_sph/rhs.jl | 15 ++- src/schemes/structure/structure.jl | 10 +- 7 files changed, 94 insertions(+), 74 deletions(-) diff --git a/src/schemes/boundary/open_boundary/rhs.jl b/src/schemes/boundary/open_boundary/rhs.jl index 4fd09f5f87..c74f779a9a 100644 --- a/src/schemes/boundary/open_boundary/rhs.jl +++ b/src/schemes/boundary/open_boundary/rhs.jl @@ -46,14 +46,14 @@ function interact!(dv, v_particle_system, u_particle_system, # Propagate `@inbounds` to the viscosity function, which accesses particle data 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) + @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 diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 9022c23218..838e12e2ac 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -50,11 +50,11 @@ function interact!(dv, v_particle_system, u_particle_system, correction) 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) + @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) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(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 25b3d09097..6da8b716a4 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl @@ -50,11 +50,11 @@ function interact!(dv, v_particle_system, u_particle_system, # Propagate `@inbounds` to the viscosity function, which accesses particle data 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) + @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] diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index a07b52786b..218dd96f87 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -287,11 +287,11 @@ function calculate_predicted_velocity_and_d_ii_values!(system::ImplicitIncompres grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) 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) + @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, diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 5046814293..87da250f9e 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -6,37 +6,46 @@ 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(dv_particle, - particle_system, neighbor_system, +@propagate_inbounds function 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) + v_a, v_b, grad_kernel, + viscosity_correction=1) viscosity = viscosity_model(particle_system, neighbor_system) - 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) + 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(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) - return 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) +@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(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) +@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 @@ -101,7 +110,8 @@ end pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - v_a, v_b, grad_kernel) + v_a, v_b, grad_kernel, + viscosity_correction=1) rho_mean = (rho_a + rho_b) / 2 v_visc_a = viscous_velocity(v_particle_system, particle_system, particle, v_a) @@ -120,7 +130,8 @@ end smoothing_length_neighbor, sound_speed) viscosity(dv_particle, sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length_average, grad_kernel, nu_a, nu_b, m_b) + smoothing_length_average, grad_kernel, nu_a, nu_b, m_b, + viscosity_correction) return dv_particle end @@ -128,7 +139,8 @@ end @inline function (viscosity::ArtificialViscosityMonaghan)(dv_particle, c, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b, m_b) + grad_kernel, nu_a, nu_b, m_b, + viscosity_correction=1) (; alpha, beta, epsilon) = viscosity # v_ab ⋅ r_ab @@ -140,7 +152,8 @@ end # viscosity is used for shocks and not rarefactions." if vr < 0 mu = h * vr / (distance^2 + epsilon * h^2) - dv_particle[] += m_b * (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel + dv_particle[] += viscosity_correction * m_b * + (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel end return dv_particle @@ -148,13 +161,15 @@ end @inline function (viscosity::ViscosityMorris)(dv_particle, c, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b, m_b) + grad_kernel, nu_a, nu_b, m_b, + viscosity_correction=1) epsilon = viscosity.epsilon mu_a = nu_a * rho_a mu_b = nu_b * rho_b - dv_particle[] += m_b * (mu_a + mu_b) / (rho_a * rho_b) * dot(pos_diff, grad_kernel) / + dv_particle[] += viscosity_correction * m_b * (mu_a + mu_b) / + (rho_a * rho_b) * dot(pos_diff, grad_kernel) / (distance^2 + epsilon * h^2) * v_diff return dv_particle @@ -193,7 +208,8 @@ end @inline function 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) + v_diff, nu_a, nu_b, epsilon, + viscosity_correction=1) eta_a = nu_a * rho_a eta_b = nu_b * rho_b @@ -215,7 +231,7 @@ end # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 visc = (volume_a^2 + volume_b^2) * dot(grad_kernel, pos_diff) * tmp / m_a - dv_particle[] += visc .* v_diff + dv_particle[] += viscosity_correction * visc .* v_diff return dv_particle end @@ -225,7 +241,8 @@ end 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) + rho_a, rho_b, v_a, v_b, grad_kernel, + viscosity_correction=1) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -245,7 +262,7 @@ end 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) + v_diff, nu_a, nu_b, epsilon, viscosity_correction) end @inline function kinematic_viscosity(system, viscosity::ViscosityAdami, smoothing_length, @@ -312,7 +329,8 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - v_a, v_b, grad_kernel) + v_a, v_b, grad_kernel, + viscosity_correction=1) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -360,7 +378,7 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps 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) + v_diff, nu_a, nu_b, epsilon, viscosity_correction) end function kinematic_viscosity(system, viscosity::ViscosityAdamiSGS, smoothing_length, @@ -427,7 +445,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - v_a, v_b, grad_kernel) + v_a, v_b, grad_kernel, + viscosity_correction=1) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -458,7 +477,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e mu_a = nu_a_eff * rho_a mu_b = nu_b_eff * rho_b - dv_particle[] += m_b * (mu_a + mu_b) / (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / + dv_particle[] += viscosity_correction * m_b * (mu_a + mu_b) / + (rho_a * rho_b) * (dot(pos_diff, grad_kernel)) / (distance^2 + epsilon * smoothing_length_average^2) * v_diff return dv_particle @@ -506,7 +526,8 @@ end pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, - v_a, v_b, grad_kernel) + v_a, v_b, grad_kernel, + viscosity_correction=1) epsilon = viscosity.epsilon smoothing_length_particle = smoothing_length(particle_system, particle) @@ -527,7 +548,7 @@ end 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) + v_diff, nu_a, nu_b, epsilon, viscosity_correction) end @inline function kinematic_viscosity(system, viscosity::ViscosityCarreauYasuda, diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 000d08ff57..075850b9d2 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -67,11 +67,11 @@ function interact!(dv, v_particle_system, u_particle_system, # Propagate `@inbounds` to the viscosity function, which accesses particle data 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) + @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, viscosity_correction) # Extra terms in the momentum equation when using a shifting technique dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), @@ -89,9 +89,8 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - dv_particle = dv_pressure + viscosity_correction * dv_viscosity_[] + dv_tvf + - dv_surface_tension + - dv_adhesion + dv_particle = dv_pressure + dv_viscosity_[] + dv_tvf + + dv_surface_tension + dv_adhesion for i in 1:ndims(particle_system) @inbounds dv[i, particle] += dv_particle[i] diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl index a3b216657d..4a6091877f 100644 --- a/src/schemes/structure/structure.jl +++ b/src/schemes/structure/structure.jl @@ -70,11 +70,11 @@ function interact_structure_fluid!(dv, v_particle_system, system_correction(neighbor_system)) dv_viscosity_ = Ref(zero(pos_diff)) - 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, - v_a, v_b, grad_kernel) + 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, + v_a, v_b, grad_kernel) dv_adhesion = adhesion_force(surface_tension, neighbor_system, particle_system, neighbor, particle, pos_diff, distance) From ba1551fbe639ee741546d81e9e39ce0ddaa5f80c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 17:34:04 +0200 Subject: [PATCH 04/12] Optimize viscosity --- src/schemes/fluid/viscosity.jl | 141 +++++++++++++++++++-------------- 1 file changed, 80 insertions(+), 61 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 87da250f9e..cf7c8cf449 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -75,6 +75,57 @@ 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 + (; alpha, beta, epsilon) = viscosity + + h_a = smoothing_length(particle_system, particle) + h_b = smoothing_length(particle_system, neighbor) + h = (h_a + h_b) / 2 + + rho_mean = (rho_a + rho_b) / 2 + + mu = h * vr / (distance^2 + epsilon * h^2) + c = sound_speed + dv_viscosity = m_b * (alpha * c * mu + 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) @@ -100,27 +151,24 @@ end return viscosity.nu end -@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, - 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) - rho_mean = (rho_a + rho_b) / 2 - +@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 + h = (smoothing_length_particle + smoothing_length_neighbor) / 2 nu_a = kinematic_viscosity(particle_system, viscosity_model(neighbor_system, particle_system), @@ -129,40 +177,6 @@ end viscosity_model(particle_system, neighbor_system), smoothing_length_neighbor, sound_speed) - viscosity(dv_particle, sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length_average, grad_kernel, nu_a, nu_b, m_b, - viscosity_correction) - - return dv_particle -end - -@inline function (viscosity::ArtificialViscosityMonaghan)(dv_particle, c, v_diff, - pos_diff, distance, - rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b, m_b, - viscosity_correction=1) - (; 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 - mu = h * vr / (distance^2 + epsilon * h^2) - dv_particle[] += viscosity_correction * m_b * - (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel - end - - return dv_particle -end - -@inline function (viscosity::ViscosityMorris)(dv_particle, c, v_diff, pos_diff, - distance, rho_mean, rho_a, rho_b, h, - grad_kernel, nu_a, nu_b, m_b, - viscosity_correction=1) epsilon = viscosity.epsilon mu_a = nu_a * rho_a @@ -314,12 +328,15 @@ 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)(dv_particle, particle_system, @@ -381,8 +398,8 @@ ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityAdamiSGS(nu, C_S, eps 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 @@ -430,12 +447,14 @@ 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)(dv_particle, particle_system, @@ -484,8 +503,8 @@ ViscosityMorrisSGS(; nu, C_S=0.1, epsilon=0.001) = ViscosityMorrisSGS(nu, C_S, e 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 @@ -514,7 +533,7 @@ 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)(dv_particle, From 0a537b0b5119968fb95c41e5fa0e4ebc773b2053 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 17:54:26 +0200 Subject: [PATCH 05/12] Fix dispatch --- src/schemes/fluid/viscosity.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index cf7c8cf449..dca03e02f8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -7,12 +7,12 @@ end # 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!(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) + 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!(dv_particle, viscosity, From 09092b669ea914e121abed9f74168cb792e9f50d Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 18:47:41 +0200 Subject: [PATCH 06/12] Remove duplicate function --- src/schemes/fluid/viscosity.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index dca03e02f8..f308299b6a 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -189,17 +189,6 @@ end return dv_particle 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 - @doc raw""" ViscosityAdami(; nu, epsilon=0.01) From bd13b6b38ebb959fe002532545a0833bb64c5107 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 18:56:04 +0200 Subject: [PATCH 07/12] Fix tests --- test/schemes/fluid/viscosity.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/test/schemes/fluid/viscosity.jl b/test/schemes/fluid/viscosity.jl index ca50e6d5c1..d4ead2cbcc 100644 --- a/test/schemes/fluid/viscosity.jl +++ b/test/schemes/fluid/viscosity.jl @@ -25,10 +25,16 @@ grad_kernel = TrixiParticles.smoothing_kernel_grad(system_wcsph, pos_diff, distance, 1) + 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, sound_speed, v_diff, pos_diff, distance, - rho_mean, rho_a, rho_b, smoothing_length, - grad_kernel, 0.0, 0.0, 1.0) + 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) From 7e3d5c49e83fa54fb13d71fbaab2c9c2cdbacb95 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 18:58:41 +0200 Subject: [PATCH 08/12] Propagate inbounds into `viscous_velocity` functions --- src/schemes/boundary/wall_boundary/system.jl | 5 +++-- src/schemes/structure/rigid_body/system.jl | 3 ++- src/schemes/structure/total_lagrangian_sph/system.jl | 3 ++- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index 6fc3f48d2c..a176be0327 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -145,7 +145,8 @@ end return zero(SVector{ndims(system), eltype(system)}) end -@inline function viscous_velocity(v, system::WallBoundarySystem, particle, v_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 @@ -155,7 +156,7 @@ end return v_particle end -@inline function viscous_velocity(v, viscosity, system, particle, v_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) diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index 82279cf949..0fda0f1cea 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -506,7 +506,8 @@ end return neighbor_system.boundary_model.viscosity end -@inline function viscous_velocity(v, system::RigidBodySystem, particle, v_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, diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl index 37b86ed85c..aea2cade98 100644 --- a/src/schemes/structure/total_lagrangian_sph/system.jl +++ b/src/schemes/structure/total_lagrangian_sph/system.jl @@ -333,7 +333,8 @@ end return current_density(v, system.boundary_model, system) end -@inline function viscous_velocity(v, system::TotalLagrangianSPHSystem, particle, v_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, From fada74277078beadea3ec2d3273242a64df97b5a Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 19:01:16 +0200 Subject: [PATCH 09/12] Reformat --- src/schemes/fluid/viscosity.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index f308299b6a..795e151b57 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -94,7 +94,8 @@ end particle, neighbor, pos_diff, distance, sound_speed, - m_a, m_b, rho_a, rho_b, + 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) From 798a97e2ec8d727f8e881190e3bd5f389bc5c301 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 7 Apr 2026 19:16:13 +0200 Subject: [PATCH 10/12] Implement copilot suggestions --- src/schemes/fluid/viscosity.jl | 12 ++++++------ src/schemes/structure/structure.jl | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 795e151b57..a8a2abd581 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -113,7 +113,7 @@ end (; alpha, beta, epsilon) = viscosity h_a = smoothing_length(particle_system, particle) - h_b = smoothing_length(particle_system, neighbor) + h_b = smoothing_length(neighbor_system, neighbor) h = (h_a + h_b) / 2 rho_mean = (rho_a + rho_b) / 2 @@ -168,7 +168,7 @@ end 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_neighbor = smoothing_length(neighbor_system, neighbor) h = (smoothing_length_particle + smoothing_length_neighbor) / 2 nu_a = kinematic_viscosity(particle_system, @@ -250,7 +250,7 @@ end 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, @@ -341,7 +341,7 @@ end 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, @@ -459,7 +459,7 @@ end 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, @@ -540,7 +540,7 @@ end 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_a) diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl index 4a6091877f..3d904acfcf 100644 --- a/src/schemes/structure/structure.jl +++ b/src/schemes/structure/structure.jl @@ -74,7 +74,7 @@ function interact_structure_fluid!(dv, v_particle_system, v_neighbor_system, v_particle_system, neighbor, particle, pos_diff, distance, sound_speed, m_b, m_a, rho_b, rho_a, - v_a, v_b, grad_kernel) + v_b, v_a, grad_kernel) dv_adhesion = adhesion_force(surface_tension, neighbor_system, particle_system, neighbor, particle, pos_diff, distance) From 8e81dc040ea3a05dbad8b2ac113f35eaa1d80eb3 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 8 Apr 2026 10:23:24 +0200 Subject: [PATCH 11/12] Fix TLSPH viscosity --- src/schemes/fluid/viscosity.jl | 3 +- .../structure/total_lagrangian_sph/rhs.jl | 11 ++- .../total_lagrangian_sph/viscosity.jl | 84 +++++++++++-------- 3 files changed, 53 insertions(+), 45 deletions(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index a8a2abd581..e3f4832ae6 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -110,14 +110,13 @@ end # approaching particles and turn it off for receding particles. In this way, the # viscosity is used for shocks and not rarefactions." if vr < 0 - (; alpha, beta, epsilon) = viscosity - 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 mu = h * vr / (distance^2 + epsilon * h^2) c = sound_speed dv_viscosity = m_b * (alpha * c * mu + beta * mu^2) / rho_mean * grad_kernel diff --git a/src/schemes/structure/total_lagrangian_sph/rhs.jl b/src/schemes/structure/total_lagrangian_sph/rhs.jl index a0fb522b4d..192dcf00b6 100644 --- a/src/schemes/structure/total_lagrangian_sph/rhs.jl +++ b/src/schemes/structure/total_lagrangian_sph/rhs.jl @@ -56,14 +56,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/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 From 792194d0487838a65e58a12aeb588f4bc790f060 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 10 Apr 2026 17:53:20 +0200 Subject: [PATCH 12/12] Reformat --- src/schemes/fluid/viscosity.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index f3318f9433..b39ccb0ef4 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -123,7 +123,8 @@ end 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_viscosity = div_fast(m_b * alpha * c * mu + m_b * beta * mu^2, rho_mean) * + grad_kernel dv_particle[] += viscosity_correction * dv_viscosity end