Skip to content

Commit 87c9034

Browse files
efaulhaberefaulhaber
committed
Make viscosity inplace
1 parent 7578733 commit 87c9034

2 files changed

Lines changed: 42 additions & 37 deletions

File tree

‎src/schemes/fluid/viscosity.jl‎

Lines changed: 17 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -6,31 +6,33 @@ end
66
# Unpack the neighboring systems viscosity to dispatch on the viscosity type.
77
# This function is only necessary to allow `nothing` as viscosity.
88
# Otherwise, we could just apply the viscosity as a function directly.
9-
@propagate_inbounds function dv_viscosity(particle_system, neighbor_system,
9+
@propagate_inbounds function dv_viscosity(dv_particle, particle_system, neighbor_system,
1010
v_particle_system, v_neighbor_system,
1111
particle, neighbor, pos_diff, distance,
1212
sound_speed, m_a, m_b, rho_a, rho_b,
1313
v_a, v_b, grad_kernel)
1414
viscosity = viscosity_model(particle_system, neighbor_system)
1515

16-
return dv_viscosity(viscosity, particle_system, neighbor_system,
16+
return dv_viscosity(dv_particle, viscosity, particle_system, neighbor_system,
1717
v_particle_system, v_neighbor_system,
1818
particle, neighbor, pos_diff, distance,
1919
sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
2020
end
2121

22-
@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system,
22+
@propagate_inbounds function dv_viscosity(dv_particle, viscosity, particle_system,
23+
neighbor_system,
2324
v_particle_system, v_neighbor_system,
2425
particle, neighbor, pos_diff, distance,
2526
sound_speed, m_a, m_b, rho_a, rho_b,
2627
v_a, v_b, grad_kernel)
27-
return viscosity(particle_system, neighbor_system,
28+
return viscosity(dv_particle, particle_system, neighbor_system,
2829
v_particle_system, v_neighbor_system,
2930
particle, neighbor, pos_diff, distance,
3031
sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
3132
end
3233

33-
@inline function dv_viscosity(viscosity::Nothing, particle_system, neighbor_system,
34+
@inline function dv_viscosity(dv_particle, viscosity::Nothing, particle_system,
35+
neighbor_system,
3436
v_particle_system, v_neighbor_system,
3537
particle, neighbor, pos_diff, distance,
3638
sound_speed, m_a, m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
@@ -89,7 +91,8 @@ end
8991
end
9092

9193
@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan,
92-
ViscosityMorris})(particle_system,
94+
ViscosityMorris})(dv_particle,
95+
particle_system,
9396
neighbor_system,
9497
v_particle_system,
9598
v_neighbor_system,
@@ -115,15 +118,14 @@ end
115118
viscosity_model(particle_system, neighbor_system),
116119
smoothing_length_neighbor, sound_speed)
117120

118-
pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
119-
smoothing_length_average, grad_kernel, nu_a, nu_b)
120-
121-
return m_b * pi_ab
121+
viscosity(dv_particle, sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b,
122+
smoothing_length_average, grad_kernel, nu_a, nu_b, m_b)
122123
end
123124

124-
@inline function (viscosity::ArtificialViscosityMonaghan)(c, v_diff, pos_diff, distance,
125+
@inline function (viscosity::ArtificialViscosityMonaghan)(dv_particle, c, v_diff, pos_diff,
126+
distance,
125127
rho_mean, rho_a, rho_b, h,
126-
grad_kernel, nu_a, nu_b)
128+
grad_kernel, nu_a, nu_b, m_b)
127129
(; alpha, beta, epsilon) = viscosity
128130

129131
# v_ab ⋅ r_ab
@@ -137,10 +139,11 @@ end
137139
# Since this is one of the most performance critical functions, using fast divisions
138140
# here gives a significant speedup on GPUs.
139141
mu = div_fast(h * vr, distance^2 + epsilon * h^2)
140-
return div_fast(alpha * c * mu + beta * mu^2, rho_mean) * grad_kernel
142+
dv_particle[] += m_b * div_fast(alpha * c * mu + beta * mu^2, rho_mean) *
143+
grad_kernel
141144
end
142145

143-
return zero(v_diff)
146+
return dv_particle
144147
end
145148

146149
@inline function (viscosity::ViscosityMorris)(c, v_diff, pos_diff, distance, rho_mean,

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

Lines changed: 25 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -80,42 +80,44 @@ function interact!(dv, v_particle_system, u_particle_system,
8080
distance, grad_kernel, correction)
8181

8282
# Propagate `@inbounds` to the viscosity function, which accesses particle data
83-
dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system,
84-
v_particle_system, v_neighbor_system,
85-
particle, neighbor, pos_diff, distance,
86-
sound_speed, m_a, m_b, rho_a, rho_b,
87-
v_a, v_b, grad_kernel)
83+
@inbounds dv_viscosity(dv_particle, particle_system, neighbor_system,
84+
v_particle_system, v_neighbor_system,
85+
particle, neighbor, pos_diff, distance,
86+
sound_speed, m_a, m_b, rho_a, rho_b,
87+
v_a, v_b, grad_kernel)
8888

8989
# Extra terms in the momentum equation when using a shifting technique
90-
dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system),
91-
particle_system, neighbor_system,
92-
v_particle_system, v_neighbor_system,
93-
particle, neighbor, m_a, m_b, rho_a, rho_b,
94-
pos_diff, distance, grad_kernel, correction)
90+
# dv_tvf = @inbounds dv_shifting(shifting_technique(particle_system),
91+
# particle_system, neighbor_system,
92+
# v_particle_system, v_neighbor_system,
93+
# particle, neighbor, m_a, m_b, rho_a, rho_b,
94+
# pos_diff, distance, grad_kernel, correction)
9595

96-
dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b,
97-
particle_system, neighbor_system,
98-
particle, neighbor, pos_diff, distance,
99-
rho_a, rho_b, grad_kernel)
96+
# dv_surface_tension = surface_tension_force(surface_tension_a, surface_tension_b,
97+
# particle_system, neighbor_system,
98+
# particle, neighbor, pos_diff, distance,
99+
# rho_a, rho_b, grad_kernel)
100100

101-
dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
102-
particle, neighbor, pos_diff, distance)
101+
# dv_adhesion = adhesion_force(surface_tension_a, particle_system, neighbor_system,
102+
# particle, neighbor, pos_diff, distance)
103103

104104
# Determine correction factors.
105105
# This can usually be ignored, as these are all 1 when no correction is used.
106106
(viscosity_correction, pressure_correction,
107-
surface_tension_correction) = free_surface_correction(correction, particle_system,
108-
rho_mean)
107+
surface_tension_correction) = free_surface_correction(correction,
108+
particle_system,
109+
rho_a, rho_b)
109110

110111
# Accumulate contributions over all neighbors
111-
dv_particle[] += dv_pressure * pressure_correction +
112-
dv_viscosity_ * viscosity_correction +
113-
dv_tvf + dv_adhesion +
114-
dv_surface_tension * surface_tension_correction
112+
dv_particle[] += dv_pressure * pressure_correction
113+
# dv_viscosity_ * viscosity_correction +
114+
# dv_tvf + dv_adhesion +
115+
# dv_surface_tension * surface_tension_correction
115116

116117
# TODO If variable smoothing_length is used, this should use the neighbor smoothing length
117118
# Propagate `@inbounds` to the continuity equation, which accesses particle data
118-
@inbounds continuity_equation!(drho_particle, density_calculator, particle_system,
119+
@inbounds continuity_equation!(drho_particle, density_calculator,
120+
particle_system,
119121
neighbor_system, v_particle_system,
120122
v_neighbor_system, particle, neighbor,
121123
pos_diff, distance, m_b, rho_a, rho_b, vdiff,

0 commit comments

Comments
 (0)