Skip to content

Commit 6faebb0

Browse files
efaulhaberefaulhaber
committed
Rewrite fluid interact kernel
1 parent 6bc3034 commit 6faebb0

1 file changed

Lines changed: 106 additions & 127 deletions

File tree

  • src/schemes/fluid/weakly_compressible_sph

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

Lines changed: 106 additions & 127 deletions
Original file line numberDiff line numberDiff line change
@@ -14,146 +14,125 @@ function interact!(dv, v_particle_system, u_particle_system,
1414

1515
system_coords = current_coordinates(u_particle_system, particle_system)
1616
neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system)
17-
18-
# In order to visualize quantities like pressure forces or viscosity forces, uncomment
19-
# the following code and the two other lines below that are marked as "debug example".
20-
# debug_array = zeros(ndims(particle_system), nparticles(particle_system))
21-
22-
# For `distance == 0`, the analytical gradient is zero, but the unsafe gradient
23-
# and the density diffusion divide by zero.
24-
# To account for rounding errors, we check if `distance` is almost zero.
25-
# Since the coordinates are in the order of the smoothing length `h`, `distance^2` is in
26-
# the order of `h^2`, so we need to check `distance < sqrt(eps(h^2))`.
27-
# Note that `sqrt(eps(h^2)) != eps(h)`.
28-
h = initial_smoothing_length(particle_system)
29-
almostzero = sqrt(eps(h^2))
30-
31-
# Loop over all pairs of particles and neighbors within the kernel cutoff
32-
foreach_point_neighbor(particle_system, neighbor_system,
33-
system_coords, neighbor_system_coords, semi;
34-
points=each_integrated_particle(particle_system)) do particle,
35-
neighbor,
36-
pos_diff,
37-
distance
38-
# Skip neighbors with the same position because the kernel gradient is zero.
39-
# Note that `return` only exits the closure, i.e., skips the current neighbor.
40-
skip_zero_distance(particle_system) && distance < almostzero && return
41-
42-
# Now that we know that `distance` is not zero, we can safely call the unsafe
43-
# version of the kernel gradient to avoid redundant zero checks.
44-
grad_kernel = smoothing_kernel_grad_unsafe(particle_system, pos_diff,
45-
distance, particle)
46-
47-
# `foreach_point_neighbor` makes sure that `particle` and `neighbor` are
48-
# in bounds of the respective system. For performance reasons, we use `@inbounds`
49-
# in this hot loop to avoid bounds checking when extracting particle quantities.
50-
rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
51-
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
52-
rho_mean = (rho_a + rho_b) / 2
53-
54-
v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
55-
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
56-
57-
# Determine correction factors.
58-
# This can be ignored, as these are all 1 when no correction is used.
59-
(viscosity_correction, pressure_correction,
60-
surface_tension_correction) = free_surface_correction(correction, particle_system,
61-
rho_mean)
62-
17+
neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, semi)
18+
19+
# For `distance == 0`, the analytical gradient is zero, but the unsafe gradient divides
20+
# by zero. To account for rounding errors, we check if `distance` is almost zero.
21+
# Since the coordinates are in the order of the compact support `c`, `distance^2` is in
22+
# the order of `c^2`, so we need to check `distance < sqrt(eps(c^2))`.
23+
# Note that `sqrt(eps(c^2)) != eps(c)`.
24+
compact_support_ = compact_support(particle_system, neighbor_system)
25+
almostzero = sqrt(eps(compact_support_^2))
26+
27+
@threaded semi for particle in each_integrated_particle(particle_system)
28+
# We are looping over the particles of `particle_system`, so it is guaranteed
29+
# that `particle` is in bounds of `particle_system`.
6330
m_a = @inbounds hydrodynamic_mass(particle_system, particle)
64-
m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
31+
p_a = @inbounds current_pressure(v_particle_system, particle_system, particle)
6532

6633
v_a = @inbounds current_velocity(v_particle_system, particle_system, particle)
67-
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
68-
69-
# The following call is equivalent to
70-
# `p_a = current_pressure(v_particle_system, particle_system, particle)`
71-
# `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
72-
# Only when the neighbor system is a `WallBoundarySystem` or a `TotalLagrangianSPHSystem`
73-
# with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is
74-
# the pressure of the fluid particle.
75-
p_a,
76-
p_b = @inbounds particle_neighbor_pressure(v_particle_system,
77-
v_neighbor_system,
78-
particle_system, neighbor_system,
79-
particle, neighbor)
80-
81-
dv_pressure = pressure_correction *
82-
pressure_acceleration(particle_system, neighbor_system,
83-
particle, neighbor,
84-
m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
85-
distance, grad_kernel, correction)
86-
87-
dv_particle = Ref(dv_pressure)
88-
89-
# Propagate `@inbounds` to the viscosity function, which accesses particle data
90-
@inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
91-
v_particle_system, v_neighbor_system,
92-
particle, neighbor, pos_diff, distance,
93-
sound_speed, m_a, m_b, rho_a, rho_b,
94-
v_a, v_b, grad_kernel, viscosity_correction)
95-
96-
# Extra terms in the momentum equation when using a shifting technique
97-
@inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
98-
particle_system, neighbor_system,
99-
v_particle_system, v_neighbor_system,
100-
particle, neighbor, m_a, m_b, rho_a, rho_b,
101-
pos_diff, distance, grad_kernel, correction)
102-
103-
@inbounds surface_tension_force!(dv_particle,
104-
surface_tension_a, surface_tension_b,
105-
particle_system, neighbor_system,
106-
particle, neighbor, pos_diff, distance,
107-
rho_a, rho_b, grad_kernel,
108-
surface_tension_correction)
109-
110-
@inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system,
111-
neighbor_system,
112-
particle, neighbor, pos_diff, distance)
113-
114-
for i in 1:ndims(particle_system)
115-
@inbounds dv[i, particle] += dv_particle[][i]
116-
# Debug example
117-
# debug_array[i, particle] += dv_pressure[i]
118-
end
34+
rho_a = @inbounds current_density(v_particle_system, particle_system, particle)
11935

36+
# Accumulate the RHS contributions over all neighbors before writing to `dv`,
37+
# to reduce the number of memory writes.
38+
# Note that we need a `Ref` in order to be able to update these variables
39+
# inside the closure in the `foreach_neighbor` loop.
40+
dv_particle = Ref(zero(v_a))
12041
drho_particle = Ref(zero(rho_a))
12142

122-
# TODO If variable smoothing_length is used, this should use the neighbor smoothing length
123-
# Propagate `@inbounds` to the continuity equation, which accesses particle data
124-
@inbounds continuity_equation!(drho_particle, density_calculator,
125-
particle_system, neighbor_system,
126-
v_particle_system, v_neighbor_system,
127-
particle, neighbor, pos_diff, distance,
128-
m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
43+
# Loop over all neighbors within the kernel cutoff
44+
@inbounds PointNeighbors.foreach_neighbor(system_coords, neighbor_system_coords,
45+
neighborhood_search,
46+
particle) do particle, neighbor,
47+
pos_diff, distance
48+
# Skip neighbors with the same position because the kernel gradient is zero.
49+
# Note that `return` only exits the closure, i.e., skips the current neighbor.
50+
skip_zero_distance(particle_system) && distance < almostzero && return
51+
52+
# Now that we know that `distance` is not zero, we can safely call the unsafe
53+
# version of the kernel gradient to avoid redundant zero checks.
54+
grad_kernel = smoothing_kernel_grad_unsafe(particle_system, pos_diff,
55+
distance, particle)
56+
57+
# `foreach_neighbor` makes sure that `neighbor` is in bounds of `neighbor_system`
58+
m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor)
59+
v_b = @inbounds current_velocity(v_neighbor_system, neighbor_system, neighbor)
60+
rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor)
61+
rho_mean = (rho_a + rho_b) / 2
62+
63+
# The following call is equivalent to
64+
# `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
65+
# Only when the neighbor system is a `WallBoundarySystem`
66+
# or a `TotalLagrangianSPHSystem` with the boundary model `PressureMirroring`,
67+
# this will return `p_b = p_a`, which is the pressure of the fluid particle.
68+
p_b = @inbounds neighbor_pressure(v_neighbor_system, neighbor_system,
69+
neighbor, p_a)
70+
71+
# Determine correction factors.
72+
# This can usually be ignored, as these are all 1 when no correction is used.
73+
(viscosity_correction, pressure_correction,
74+
surface_tension_correction) = free_surface_correction(correction,
75+
particle_system,
76+
rho_mean)
77+
78+
# For `ContinuityDensity` without correction, this is equivalent to
79+
# dv_pressure = -m_b * (p_a + p_b) / (rho_a * rho_b) * grad_kernel
80+
dv_pressure = pressure_acceleration(particle_system, neighbor_system,
81+
particle, neighbor,
82+
m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff,
83+
distance, grad_kernel, correction)
84+
dv_particle[] += dv_pressure * pressure_correction
85+
86+
# Propagate `@inbounds` to the viscosity function, which accesses particle data
87+
@inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system,
88+
v_particle_system, v_neighbor_system,
89+
particle, neighbor, pos_diff, distance,
90+
sound_speed, m_a, m_b, rho_a, rho_b,
91+
v_a, v_b, grad_kernel, viscosity_correction)
92+
93+
# Extra terms in the momentum equation when using a shifting technique
94+
@inbounds dv_shifting!(dv_particle, shifting_technique(particle_system),
95+
particle_system, neighbor_system,
96+
v_particle_system, v_neighbor_system,
97+
particle, neighbor, m_a, m_b, rho_a, rho_b,
98+
pos_diff, distance, grad_kernel, correction)
99+
100+
@inbounds surface_tension_force!(dv_particle,
101+
surface_tension_a, surface_tension_b,
102+
particle_system, neighbor_system,
103+
particle, neighbor, pos_diff, distance,
104+
rho_a, rho_b, grad_kernel,
105+
surface_tension_correction)
106+
107+
@inbounds adhesion_force!(dv_particle, surface_tension_a, particle_system,
108+
neighbor_system,
109+
particle, neighbor, pos_diff, distance)
110+
111+
# TODO If variable smoothing_length is used, this should use the neighbor smoothing length
112+
# Propagate `@inbounds` to the continuity equation, which accesses particle data
113+
@inbounds continuity_equation!(drho_particle, density_calculator,
114+
particle_system, neighbor_system,
115+
v_particle_system, v_neighbor_system,
116+
particle, neighbor, pos_diff, distance,
117+
m_b, rho_a, rho_b, v_a, v_b, grad_kernel)
118+
end
129119

120+
for i in eachindex(dv_particle[])
121+
@inbounds dv[i, particle] += dv_particle[][i]
122+
end
130123
@inbounds write_drho_particle!(dv, density_calculator, drho_particle, particle)
131124
end
132-
# Debug example
133-
# periodic_box = neighborhood_search.periodic_box
134-
# Note: this saves a file in every stage of the integrator
135-
# if !@isdefined iter; iter = 0; end
136-
# TODO: This call should use public API. This requires some additional changes to simplify the calls.
137-
# trixi2vtk(v_particle_system, u_particle_system, -1.0, particle_system, periodic_box, debug=debug_array, prefix="debug", iter=iter += 1)
138125

139126
return dv
140127
end
141128

142-
@propagate_inbounds function particle_neighbor_pressure(v_particle_system,
143-
v_neighbor_system,
144-
particle_system, neighbor_system,
145-
particle, neighbor)
146-
p_a = current_pressure(v_particle_system, particle_system, particle)
147-
p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)
148-
149-
return p_a, p_b
129+
@propagate_inbounds function neighbor_pressure(v_neighbor_system, neighbor_system,
130+
neighbor, p_a)
131+
return current_pressure(v_neighbor_system, neighbor_system, neighbor)
150132
end
151133

152-
@inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system,
153-
particle_system,
154-
neighbor_system::WallBoundarySystem{<:BoundaryModelDummyParticles{PressureMirroring}},
155-
particle, neighbor)
156-
p_a = current_pressure(v_particle_system, particle_system, particle)
157-
158-
return p_a, p_a
134+
@inline function neighbor_pressure(v_neighbor_system,
135+
neighbor_system::WallBoundarySystem{<:BoundaryModelDummyParticles{PressureMirroring}},
136+
neighbor, p_a)
137+
return p_a
159138
end

0 commit comments

Comments
 (0)