Add basic unit test for density

src/general/density_calculators.jl

@@ -59,19 +59,18 @@ function summation_density!(system, system_index, semi, u, u_ode, density;
     # Use all other systems for the density summation
     @trixi_timeit timer() "compute density" foreach_enumerate(systems) do (neighbor_system_index,
                                                                            neighbor_system)
-        u_neighbor_system = wrap_u(u_ode, neighbor_system_index,
-                                   neighbor_system, semi)
+        u_neighbor_system = wrap_u(u_ode, neighbor_system_index, neighbor_system, semi)
 
         system_coords = current_coordinates(u, system)
         neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
 
         neighborhood_search = neighborhood_searches[system_index][neighbor_system_index]
 
         # Loop over all pairs of particles and neighbors within the kernel cutoff.
-        for_particle_neighbor(system, neighbor_system,
-                              system_coords, neighbor_coords,
+        for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords,
                               neighborhood_search,
-                              particles=particles) do particle, neighbor, pos_diff, distance
+                              particles=particles) do particle,
+                                                      neighbor, pos_diff, distance
             mass = hydrodynamic_mass(neighbor_system, neighbor)
             density[particle] += mass * smoothing_kernel(system, distance)
         end

src/general/initial_condition.jl

@@ -66,24 +66,28 @@ struct InitialCondition{ELTYPE}
 
     function InitialCondition(coordinates, velocities, masses, densities; pressure=0.0,
                               particle_spacing=-1.0)
+        ELTYPE = eltype(coordinates)
+
         if size(coordinates) != size(velocities)
             throw(ArgumentError("`coordinates` and `velocities` must be of the same size"))
         end
 
         if !(size(coordinates, 2) == length(masses) == length(densities))
-            throw(ArgumentError("the following must hold: " *
-                                "`size(coordinates, 2) == length(masses) == length(densities)`"))
+            throw(ArgumentError("Expected: size(coordinates, 2) == length(masses) == length(densities)\n" *
+                                "Got: size(coordinates, 2) = $(size(coordinates, 2)), " *
+                                "length(masses) = $(length(masses)), " *
+                                "length(densities) = $(length(densities))"))
         end
 
         if pressure isa Number
-            pressure = pressure * ones(length(masses))
+            pressure = pressure * ones(ELTYPE, length(masses))
         elseif length(pressure) != length(masses)
             throw(ArgumentError("`pressure` must either be a scalar or a vector of the " *
                                 "same length as `masses`"))
         end
 
-        return new{eltype(coordinates)}(particle_spacing, coordinates, velocities, masses,
-                                        densities, pressure)
+        return new{ELTYPE}(particle_spacing, coordinates, velocities, masses,
+                           densities, pressure)
     end
 end
 

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -101,6 +101,7 @@ function create_cache(n_particles, ELTYPE, ::SummationDensity)
 end
 
 function create_cache(n_particles, ELTYPE, ::ContinuityDensity)
+    # Density in this case is added to the end of 'v' and allocated by modifying 'v_nvariables'.
     return (;)
 end
 

test/general/density_calculator.jl

@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+
+# setup a single particle and calculate its density
+@testset verbose=true "DensityCalculators" begin
+    @testset verbose=true "SummationDensity" begin
+        water_density = 1000.0
+
+        # use reshape to create a matrix
+        init = InitialCondition(reshape([0.0, 0.0], 2, 1), reshape([0.0, 0.0], 2, 1),
+                                [water_density], [water_density])
+
+        smoothing_length = 1.0
+        smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+
+        state_equation = StateEquationCole(10, 7, water_density, 100000.0,
+                                           background_pressure=100000.0,
+                                           clip_negative_pressure=false)
+        viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
+
+        fluid_system = WeaklyCompressibleSPHSystem(init, SummationDensity(), state_equation,
+                                                   smoothing_kernel, smoothing_length,
+                                                   viscosity=viscosity,
+                                                   acceleration=(0.0, 0.0))
+
+        (; cache) = fluid_system
+        (; density) = cache # Density is in the cache for SummationDensity
+
+        semi = Semidiscretization(fluid_system, neighborhood_search=GridNeighborhoodSearch,
+                                  damping_coefficient=1e-5)
+
+        tspan = (0.0, 0.01)
+        ode = semidiscretize(semi, tspan)
+
+        sol = solve(ode, RDPK3SpFSAL49(),
+                    abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+                    reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+                    dtmax=1e-2, # Limit stepsize to prevent crashing
+                    save_everystep=false)
+
+        @test density[1] ===
+              water_density * TrixiParticles.kernel(smoothing_kernel, 0.0, smoothing_length)
+    end
+end

test/general/general.jl

@@ -1,4 +1,5 @@
 # `util.jl` is added here to avoid confusion with `test_util.jl`
 include("util.jl")
 include("initial_condition.jl")
+include("density_calculator.jl")
 include("semidiscretization.jl")