Fix tests and add comments

Author: efaulhaber

src/schemes/fluid/implicit_incompressible_sph/rhs.jl

@@ -47,17 +47,14 @@ 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)
 
+        p_a = @inbounds current_pressure(v_particle_system, particle_system, particle)
         # The following call is equivalent to
-        #     `p_a = particle_pressure(v_particle_system, particle_system, particle)`
-        #     `p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)`
-        # Only when the neighbor system is a `WallBoundarySystem` or a `TotalLagrangianSPHSystem`
-        # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is
-        # the pressure of the fluid particle.
-        p_a,
-        p_b = @inbounds particle_neighbor_pressure(v_particle_system,
-                                                   v_neighbor_system,
-                                                   particle_system, neighbor_system,
-                                                   particle, neighbor)
+        #     `p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)`
+        # Only when the neighbor system is a `WallBoundarySystem`
+        # or a `TotalLagrangianSPHSystem` with the boundary model `PressureMirroring`,
+        # this will return `p_b = p_a`, which is the pressure of the fluid particle.
+        p_b = @inbounds neighbor_pressure(v_neighbor_system, neighbor_system,
+                                          neighbor, p_a)
 
         dv_pressure = pressure_acceleration(particle_system, neighbor_system,
                                             particle, neighbor,

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -147,19 +147,28 @@ end
     return velocity_and_density(v, density_calculator, system, particle)
 end
 
+@propagate_inbounds function velocity_and_density(v, system, particle)
+    # Call the default method below
+    return velocity_and_density(v, nothing, system, particle)
+end
+
+# Default method, which simply calls `current_velocity` and `current_density` separately.
 @propagate_inbounds function velocity_and_density(v, _, system, particle)
     v_particle = current_velocity(v, system, particle)
     rho_particle = current_density(v, system, particle)
 
     return v_particle, rho_particle
 end
 
+# Optimized version for WCSPH with `ContinuityDensity` in 3D on GPUs,
+# which combines the velocity and density load into one wide load.
 @inline function velocity_and_density(v::AbstractGPUArray, ::ContinuityDensity,
                                       ::WeaklyCompressibleSPHSystem{3}, particle)
     # Since `v` is stored as a 4 x N matrix, this aligned load extracts one column
     # of `v` corresponding to `particle`.
     # As opposed to `extract_svector`, this will translate to a single wide load instruction
     # on the GPU, which is faster than 4 separate loads.
+    # Note that this doesn't work for 2D because it requires a stride of 2^n.
     vrho_a = vloada(Vec{4, eltype(v)}, pointer(v, 4 * (particle - 1) + 1))
 
     # The column of `v` is ordered as (v_x, v_y, v_z, rho)