Skip to content

Commit 425a7b4

Browse files
efaulhaberefaulhaber
authored
Make viscosity inplace and avoid duplicate velocity loads (#1125)
* Avoid duplicate velocity loads for viscosity calculation * Make viscosity inplace * Rename `dv_viscosity` and add free surface correction * Optimize viscosity * Fix dispatch * Remove duplicate function * Fix tests * Propagate inbounds into `viscous_velocity` functions * Reformat * Implement copilot suggestions * Fix TLSPH viscosity * Reformat
1 parent ca35699 commit 425a7b4

16 files changed

Lines changed: 401 additions & 272 deletions

File tree

‎src/general/interpolation.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -593,7 +593,9 @@ end
593593
other_density[point] += m_b * W_ab
594594

595595
if include_wall_velocity
596-
velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor)
596+
velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor)
597+
velocity_neighbor = viscous_velocity(v, neighbor_system, neighbor,
598+
velocity_neighbor_)
597599
for i in axes(velocity_neighbor, 1)
598600
cache.velocity[i, point] += velocity_neighbor[i] * volume_b * W_ab
599601
end

‎src/schemes/boundary/open_boundary/rhs.jl‎

Lines changed: 14 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
4141
rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
4242
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
4343

44+
v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
45+
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
46+
4447
m_a = @inbounds hydrodynamic_mass(particle_system, particle)
4548
m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
4649

@@ -58,15 +61,17 @@ function interact!(dv, v_particle_system, u_particle_system,
5861
dv_pressure_boundary = 2 * p_boundary * (m_b / (rho_a * rho_b)) * grad_kernel
5962

6063
# Propagate `@inbounds` to the viscosity function, which accesses particle data
61-
dv_viscosity_ = @inbounds dv_viscosity(viscosity_model(fluid_system,
62-
neighbor_system),
63-
particle_system, neighbor_system,
64-
v_particle_system, v_neighbor_system,
65-
particle, neighbor, pos_diff, distance,
66-
sound_speed, m_a, m_b, rho_a, rho_b,
67-
grad_kernel)
68-
69-
dv_particle = dv_pressure + dv_viscosity_ + dv_pressure_boundary
64+
dv_viscosity_ = Ref(zero(pos_diff))
65+
@inbounds dv_viscosity!(dv_viscosity_,
66+
viscosity_model(fluid_system,
67+
neighbor_system),
68+
particle_system, neighbor_system,
69+
v_particle_system, v_neighbor_system,
70+
particle, neighbor, pos_diff, distance,
71+
sound_speed, m_a, m_b, rho_a, rho_b,
72+
v_a, v_b, grad_kernel)
73+
74+
dv_particle = dv_pressure + dv_viscosity_[] + dv_pressure_boundary
7075

7176
for i in 1:ndims(particle_system)
7277
@inbounds dv[i, particle] += dv_particle[i]

‎src/schemes/boundary/open_boundary/system.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -654,8 +654,10 @@ function interpolate_velocity!(system::OpenBoundarySystem, boundary_zone,
654654
W_ab = kernel(smoothing_kernel, distance, smoothing_length)
655655
@inbounds shepard_coefficient[point] += volume_b * W_ab
656656

657+
velocity_neighbor_ = @inbounds current_velocity(v_neighbor, neighbor_system,
658+
neighbor)
657659
velocity_neighbor = @inbounds viscous_velocity(v_neighbor, neighbor_system,
658-
neighbor)
660+
neighbor, velocity_neighbor_)
659661
for i in axes(velocity_neighbor, 1)
660662
@inbounds sample_velocity[i,
661663
point] += velocity_neighbor[i] * volume_b * W_ab

‎src/schemes/boundary/wall_boundary/system.jl‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -145,15 +145,20 @@ end
145145
return zero(SVector{ndims(system), eltype(system)})
146146
end
147147

148-
@inline function viscous_velocity(v, system::WallBoundarySystem, particle)
149-
return viscous_velocity(v, system.boundary_model.viscosity, system, particle)
148+
@propagate_inbounds function viscous_velocity(v, system::WallBoundarySystem,
149+
particle, v_particle)
150+
return viscous_velocity(v, system.boundary_model.viscosity, system,
151+
particle, v_particle)
150152
end
151153

152-
@inline function viscous_velocity(v, ::Nothing, system, particle)
153-
return current_velocity(v, system, particle)
154+
@inline function viscous_velocity(v, ::Nothing, system, particle, v_particle)
155+
# Regular particle velocity is used for the viscosity calculation by default
156+
return v_particle
154157
end
155158

156-
@inline function viscous_velocity(v, viscosity, system, particle)
159+
@propagate_inbounds function viscous_velocity(v, viscosity, system, particle, v_particle)
160+
# Wall velocity in the viscosity calculation contains the physical wall velocity
161+
# and an interpolated velocity when a wall viscosity (no-slip BC) is used.
157162
return extract_svector(system.boundary_model.cache.wall_velocity, system, particle)
158163
end
159164

‎src/schemes/fluid/entropically_damped_sph/rhs.jl‎

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,9 @@ function interact!(dv, v_particle_system, u_particle_system,
4242
rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
4343
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
4444

45+
v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
46+
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
47+
4548
p_a = @inbounds current_pressure(v_particle_system, particle_system, particle)
4649
p_b = @inbounds current_pressure(v_neighbor_system, neighbor_system, neighbor)
4750

@@ -62,13 +65,12 @@ function interact!(dv, v_particle_system, u_particle_system,
6265
rho_b, pos_diff, distance, grad_kernel,
6366
correction)
6467

65-
dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
66-
v_particle_system, v_neighbor_system,
67-
particle, neighbor, pos_diff, distance,
68-
sound_speed, m_a, m_b, rho_a, rho_b,
69-
grad_kernel)
70-
71-
dv_particle = Ref(dv_pressure + dv_viscosity_)
68+
dv_particle = Ref(dv_pressure)
69+
@inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
70+
v_particle_system, v_neighbor_system,
71+
particle, neighbor, pos_diff, distance,
72+
sound_speed, m_a, m_b, rho_a, rho_b,
73+
v_a, v_b, grad_kernel)
7274

7375
# Extra terms in the momentum equation when using a shifting technique
7476
@inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),

‎src/schemes/fluid/implicit_incompressible_sph/rhs.jl‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,9 @@ function interact!(dv, v_particle_system, u_particle_system,
4141
rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
4242
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
4343

44+
v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
45+
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
46+
4447
m_a = @inbounds hydrodynamic_mass(particle_system, particle)
4548
m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
4649

@@ -62,14 +65,15 @@ function interact!(dv, v_particle_system, u_particle_system,
6265
distance, grad_kernel, nothing)
6366

6467
# Propagate `@inbounds` to the viscosity function, which accesses particle data
65-
dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
66-
v_particle_system, v_neighbor_system,
67-
particle, neighbor, pos_diff, distance,
68-
sound_speed, m_a, m_b, rho_a, rho_b,
69-
grad_kernel)
68+
dv_viscosity_ = Ref(zero(pos_diff))
69+
@inbounds dv_viscosity!(dv_viscosity_, particle_system, neighbor_system,
70+
v_particle_system, v_neighbor_system,
71+
particle, neighbor, pos_diff, distance,
72+
sound_speed, m_a, m_b, rho_a, rho_b,
73+
v_a, v_b, grad_kernel)
7074

7175
for i in 1:ndims(particle_system)
72-
@inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i]
76+
@inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[][i]
7377
end
7478
end
7579
return dv

‎src/schemes/fluid/implicit_incompressible_sph/system.jl‎

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -281,16 +281,21 @@ function calculate_predicted_velocity_and_d_ii_values!(system::ImplicitIncompres
281281
rho_a = @inbounds current_density(v_particle_system, system, particle)
282282
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
283283

284+
v_a = @inbounds current_velocity(v_particle_system, system, particle)
285+
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
286+
284287
grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
285288

286-
dv_viscosity_ = @inbounds dv_viscosity(system, neighbor_system,
287-
v_particle_system, v_neighbor_system,
288-
particle, neighbor, pos_diff, distance,
289-
sound_speed, m_a, m_b, rho_a, rho_b,
290-
grad_kernel)
289+
dv_viscosity_ = Ref(zero(pos_diff))
290+
@inbounds dv_viscosity!(dv_viscosity_, system, neighbor_system,
291+
v_particle_system, v_neighbor_system,
292+
particle, neighbor, pos_diff, distance,
293+
sound_speed, m_a, m_b, rho_a, rho_b,
294+
v_a, v_b, grad_kernel)
291295
# Add all other non-pressure forces
292296
for i in 1:ndims(system)
293-
@inbounds advection_velocity[i, particle] += time_step * dv_viscosity_[i]
297+
@inbounds advection_velocity[i,
298+
particle] += time_step * dv_viscosity_[][i]
294299
end
295300
# Calculate d_ii with eq. 9 in Ihmsen et al. (2013)
296301
for i in 1:ndims(system)

0 commit comments

Comments
 (0)