diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index 56eb4ae12d..b8216db0e6 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -350,18 +350,13 @@ function write2vtk!(vtk, v, u, t, system::AbstractFluidSystem) rho_b = current_density(v, system, neighbor) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) - surface_tension[1:ndims(system), - particle] .+= surface_tension_force(surface_tension_a, - surface_tension_b, - system, - system, - particle, - neighbor, - pos_diff, - distance, - rho_a, - rho_b, - grad_kernel) + dv_surface_tension = Ref(zero(pos_diff)) + surface_tension_force!(dv_surface_tension, + surface_tension_a, surface_tension_b, + system, system, particle, neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel, 1) + + surface_tension[1:ndims(system), particle] .+= dv_surface_tension[] end vtk["surface_tension"] = surface_tension diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index cea18484d4..b91279702b 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -52,26 +52,27 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Extra terms in the momentum equation when using a shifting technique - dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), - particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, m_a, m_b, rho_a, rho_b, - pos_diff, distance, grad_kernel, correction) - - dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) + dv_particle = Ref(dv_pressure + dv_viscosity_) - 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_adhesion + # Extra terms in the momentum equation when using a shifting technique + @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, m_a, m_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) + + @inbounds surface_tension_force!(dv_particle, surface_tension_a, + surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel, 1) + + @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system, + neighbor_system, + particle, neighbor, pos_diff, distance) for i in 1:ndims(particle_system) - @inbounds dv[i, particle] += dv_particle[i] + @inbounds dv[i, particle] += dv_particle[][i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index ceaff12a21..ed9644b471 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -41,11 +41,11 @@ function update_shifting!(system, shifting, v, u, v_ode, u_ode, semi) end # Additional term in the momentum equation due to the shifting technique -@inline function dv_shifting(shifting, system, neighbor_system, - v_system, v_neighbor_system, particle, neighbor, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - return zero(grad_kernel) +@inline function dv_shifting!(dv_particle, shifting, system, neighbor_system, + v_system, v_neighbor_system, particle, neighbor, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + return dv_particle end # Additional term(s) in the continuity equation due to the shifting technique @@ -324,22 +324,23 @@ See [`ParticleShiftingTechnique`](@ref). struct MomentumEquationTermSun2019 end # Additional term in the momentum equation due to the shifting technique -@propagate_inbounds function dv_shifting(shifting::ParticleShiftingTechnique, system, - neighbor_system, - v_system, v_neighbor_system, particle, neighbor, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) - return dv_shifting(shifting.momentum_equation_term, system, neighbor_system, - v_system, v_neighbor_system, particle, neighbor, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) +@propagate_inbounds function dv_shifting!(dv_particle, + shifting::ParticleShiftingTechnique, + system, neighbor_system, + v_system, v_neighbor_system, particle, neighbor, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) + return dv_shifting!(dv_particle, shifting.momentum_equation_term, system, + neighbor_system, v_system, v_neighbor_system, + particle, neighbor, m_a, m_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) end -@propagate_inbounds function dv_shifting(::MomentumEquationTermSun2019, - system, neighbor_system, - v_system, v_neighbor_system, - particle, neighbor, m_a, m_b, rho_a, rho_b, - pos_diff, distance, grad_kernel, correction) +@propagate_inbounds function dv_shifting!(dv_particle, ::MomentumEquationTermSun2019, + system, neighbor_system, + v_system, v_neighbor_system, + particle, neighbor, m_a, m_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) delta_v_a = delta_v(system, particle) delta_v_b = delta_v(neighbor_system, neighbor) @@ -347,8 +348,11 @@ end v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) tensor_product = v_a * delta_v_a' + v_b * delta_v_b' - return m_b / rho_b * - (tensor_product * grad_kernel + v_a * dot(delta_v_a - delta_v_b, grad_kernel)) + dv_particle[] += m_b / rho_b * + (tensor_product * grad_kernel + + v_a * dot(delta_v_a - delta_v_b, grad_kernel)) + + return dv_particle end # `ParticleShiftingTechnique{<:Any, <:Any, true}` means `modify_continuity_equation=true` @@ -568,10 +572,11 @@ struct TransportVelocityAdami{modify_continuity_equation, T <: Real} <: end end -@propagate_inbounds function dv_shifting(::TransportVelocityAdami, system, neighbor_system, - v_system, v_neighbor_system, particle, neighbor, - m_a, m_b, rho_a, rho_b, pos_diff, distance, - grad_kernel, correction) +@propagate_inbounds function dv_shifting!(dv_particle, ::TransportVelocityAdami, + system, neighbor_system, + v_system, v_neighbor_system, particle, neighbor, + m_a, m_b, rho_a, rho_b, pos_diff, distance, + grad_kernel, correction) v_a = current_velocity(v_system, system, particle) delta_v_a = delta_v(system, particle) @@ -588,9 +593,11 @@ end # m_b * (A_a + A_b) / (ρ_a * ρ_b) * ∇W_ab. # In order to obtain this, we pass `p_a = A_a` and `p_b = A_b` to the # `pressure_acceleration` function. - return pressure_acceleration(system, neighbor_system, particle, neighbor, - m_a, m_b, A_a, A_b, rho_a, rho_b, pos_diff, - distance, grad_kernel, correction) + dv_particle[] += pressure_acceleration(system, neighbor_system, particle, neighbor, + m_a, m_b, A_a, A_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) + + return dv_particle end # The function above misuses the pressure acceleration function by passing a Matrix as `p_a`. diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 7a0232c51c..5656e95e12 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -162,66 +162,85 @@ end end # Skip -@inline function surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, particle, neighbor, - pos_diff, distance, rho_a, rho_b, grad_kernel) - return zero(pos_diff) +@inline function surface_tension_force!(dv_particle, surface_tension_a, surface_tension_b, + particle_system, neighbor_system, particle, + neighbor, pos_diff, distance, rho_a, rho_b, + grad_kernel, surface_tension_correction) + return dv_particle end -@inline function surface_tension_force(surface_tension_a::CohesionForceAkinci, - surface_tension_b::CohesionForceAkinci, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) +@inline function surface_tension_force!(dv_particle, + surface_tension_a::CohesionForceAkinci, + surface_tension_b::CohesionForceAkinci, + particle_system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel, + surface_tension_correction) + (; smoothing_kernel) = particle_system + # No cohesion with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) support_radius = compact_support(smoothing_kernel, smoothing_length(particle_system, particle)) - return cohesion_force_akinci(surface_tension_a, support_radius, m_b, pos_diff, distance) + dv_particle[] += surface_tension_correction * + cohesion_force_akinci(surface_tension_a, support_radius, m_b, + pos_diff, distance) + + return dv_particle end -@inline function surface_tension_force(surface_tension_a::SurfaceTensionAkinci, - surface_tension_b::SurfaceTensionAkinci, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, particle, - neighbor, - pos_diff, distance, rho_a, rho_b, grad_kernel) +@inline function surface_tension_force!(dv_particle, + surface_tension_a::SurfaceTensionAkinci, + surface_tension_b::SurfaceTensionAkinci, + particle_system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, particle, + neighbor, + pos_diff, distance, rho_a, rho_b, grad_kernel, + surface_tension_correction) (; smoothing_kernel) = particle_system (; surface_tension_coefficient) = surface_tension_a smoothing_length_ = smoothing_length(particle_system, particle) # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) n_a = surface_normal(particle_system, particle) n_b = surface_normal(neighbor_system, neighbor) support_radius = compact_support(smoothing_kernel, smoothing_length_) - return cohesion_force_akinci(surface_tension_a, support_radius, m_b, - pos_diff, distance) .- - (surface_tension_coefficient * (n_a - n_b) * smoothing_length_) -end + dv_particle[] += surface_tension_correction * + cohesion_force_akinci(surface_tension_a, support_radius, m_b, + pos_diff, distance) + dv_particle[] -= surface_tension_correction * surface_tension_coefficient * + (n_a - n_b) * smoothing_length_ -@inline function surface_tension_force(surface_tension_a::SurfaceTensionMorris, - surface_tension_b::SurfaceTensionMorris, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) + return dv_particle +end +@inline function surface_tension_force!(dv_particle, + surface_tension_a::SurfaceTensionMorris, + surface_tension_b::SurfaceTensionMorris, + particle_system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel, + surface_tension_correction) (; surface_tension_coefficient) = surface_tension_a # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle n_a = surface_normal(particle_system, particle) curvature_a = curvature(particle_system, particle) - return -surface_tension_coefficient / rho_a * curvature_a * n_a + dv_particle[] -= surface_tension_correction * surface_tension_coefficient / rho_a * + curvature_a * n_a + + return dv_particle end function compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) @@ -277,46 +296,55 @@ function compute_surface_delta_function!(system, ::SurfaceTensionMomentumMorris, return system end -@inline function surface_tension_force(surface_tension_a::SurfaceTensionMomentumMorris, - surface_tension_b::SurfaceTensionMomentumMorris, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractFluidSystem, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) +@inline function surface_tension_force!(dv_particle, + surface_tension_a::SurfaceTensionMomentumMorris, + surface_tension_b::SurfaceTensionMomentumMorris, + particle_system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel, + surface_tension_correction) (; surface_tension_coefficient) = surface_tension_a # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle S_a = stress_tensor(particle_system, particle) S_b = stress_tensor(neighbor_system, neighbor) m_b = hydrodynamic_mass(neighbor_system, neighbor) - return surface_tension_coefficient * m_b * (S_a + S_b) / (rho_a * rho_b) * grad_kernel + dv_particle[] += surface_tension_correction * surface_tension_coefficient * m_b * + (S_a + S_b) / (rho_a * rho_b) * grad_kernel + + return dv_particle end -@inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, - particle_system::AbstractFluidSystem, - neighbor_system::AbstractBoundarySystem, particle, neighbor, - pos_diff, distance) +@inline function adhesion_force!(dv_particle, + surface_tension::AkinciTypeSurfaceTension, + particle_system::AbstractFluidSystem, + neighbor_system::AbstractBoundarySystem, particle, + neighbor, + pos_diff, distance) (; adhesion_coefficient) = neighbor_system # No adhesion with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle # No reason to calculate the adhesion force if adhesion coefficient is near zero - abs(adhesion_coefficient) < eps() && return zero(pos_diff) + abs(adhesion_coefficient) < eps() && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) support_radius = compact_support(particle_system.smoothing_kernel, smoothing_length(particle_system, particle)) - return adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, - adhesion_coefficient) + dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, + distance, adhesion_coefficient) + + return dv_particle end -@inline function adhesion_force(surface_tension, particle_system, neighbor_system, particle, - neighbor, pos_diff, distance) - return zero(pos_diff) +@inline function adhesion_force!(dv_particle, surface_tension, particle_system, + neighbor_system, particle, neighbor, pos_diff, distance) + return dv_particle end diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index acb0f06894..f70891eae3 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -70,27 +70,28 @@ function interact!(dv, v_particle_system, u_particle_system, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - # Extra terms in the momentum equation when using a shifting technique - dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system), - particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, m_a, m_b, rho_a, rho_b, - pos_diff, distance, grad_kernel, correction) - - dv_surface_tension = surface_tension_correction * - surface_tension_force(surface_tension_a, surface_tension_b, - particle_system, neighbor_system, - particle, neighbor, pos_diff, distance, - rho_a, rho_b, grad_kernel) + dv_particle = Ref(dv_pressure + dv_viscosity_) - 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_adhesion + # Extra terms in the momentum equation when using a shifting technique + @inbounds dv_shifting!(dv_particle, shifting_technique(particle_system), + particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, m_a, m_b, rho_a, rho_b, + pos_diff, distance, grad_kernel, correction) + + @inbounds surface_tension_force!(dv_particle, + surface_tension_a, surface_tension_b, + particle_system, neighbor_system, + particle, neighbor, pos_diff, distance, + rho_a, rho_b, grad_kernel, + surface_tension_correction) + + @inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system, + neighbor_system, + particle, neighbor, pos_diff, distance) for i in 1:ndims(particle_system) - @inbounds dv[i, particle] += dv_particle[i] + @inbounds dv[i, particle] += dv_particle[][i] # Debug example # debug_array[i, particle] += dv_pressure[i] end diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index c9281aa5f9..20a726b7b1 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -282,23 +282,26 @@ function calc_normal!(system::AbstractFluidSystem, semi, surface_normal_method) end -@inline function adhesion_force(surface_tension::AkinciTypeSurfaceTension, - particle_system::AbstractFluidSystem, - neighbor_system::RigidBodySystem, - particle, neighbor, pos_diff, distance) +@inline function adhesion_force!(dv_particle, + surface_tension::AkinciTypeSurfaceTension, + particle_system::AbstractFluidSystem, + neighbor_system::RigidBodySystem, + particle, neighbor, pos_diff, distance) (; adhesion_coefficient) = neighbor_system # No adhesion with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return zero(pos_diff) + distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle - abs(adhesion_coefficient) < eps() && return zero(pos_diff) + abs(adhesion_coefficient) < eps() && return dv_particle m_b = hydrodynamic_mass(neighbor_system, neighbor) support_radius = compact_support(system_smoothing_kernel(particle_system), smoothing_length(particle_system, particle)) - return adhesion_force_akinci(surface_tension, support_radius, m_b, pos_diff, distance, - adhesion_coefficient) + dv_particle[] += adhesion_force_akinci(surface_tension, support_radius, m_b, + pos_diff, distance, adhesion_coefficient) + + return dv_particle end function write_u0!(u0, system::RigidBodySystem) diff --git a/src/schemes/structure/structure.jl b/src/schemes/structure/structure.jl index 1ac304d79f..5fe6d08f1d 100644 --- a/src/schemes/structure/structure.jl +++ b/src/schemes/structure/structure.jl @@ -71,12 +71,11 @@ function interact_structure_fluid!(dv, v_particle_system, neighbor, particle, pos_diff, distance, sound_speed, m_b, m_a, rho_b, rho_a, grad_kernel) - dv_adhesion = adhesion_force(surface_tension, neighbor_system, particle_system, - neighbor, particle, pos_diff, distance) + dv_particle = Ref(dv_boundary + dv_viscosity_) + adhesion_force!(dv_particle, surface_tension, neighbor_system, particle_system, + neighbor, particle, pos_diff, distance) - dv_fs = dv_boundary + dv_viscosity_ + dv_adhesion - - accumulate_structure_fluid_pair!(dv, dv_fs, particle_system, particle, m_b) + accumulate_structure_fluid_pair!(dv, dv_particle[], particle_system, particle, m_b) continuity_equation!(dv, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance,