make update! more abstract (#162)

* make update more compact

* move boundary models

* fix `SummationDensity` bug

* implement suggestions

* revert moving compute_density!()

* delete needless function

* add `update!` for `ContinuityDensity`

* fix bug

* add missing pressure

src/general/density_calculators.jl

@@ -49,3 +49,31 @@ end
 @inline function set_particle_density(particle, v, ::ContinuityDensity, system, density)
     v[end, particle] = density
 end
+
+function summation_density!(system, system_index, semi, u, u_ode, density;
+                            particles=each_moving_particle(system))
+    @unpack systems, neighborhood_searches = semi
+
+    density .= zero(eltype(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)
+
+        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,
+                              neighborhood_search,
+                              particles=particles) do particle, neighbor, pos_diff, distance
+            mass = hydrodynamic_mass(neighbor_system, neighbor)
+            density[particle] += mass * smoothing_kernel(system, distance)
+        end
+    end
+end

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -151,80 +151,31 @@ function create_cache(initial_density, ::AdamiPressureExtrapolation)
 end
 
 @inline function particle_density(v, model::BoundaryModelDummyParticles, system, particle)
-    return particle_density(v, model.density_calculator, system, particle)
+    return particle_density(v, model.density_calculator, model, particle)
 end
 
 # Note that the other density calculators are dispatched in `density_calculators.jl`
-@inline function particle_density(v, ::AdamiPressureExtrapolation, system, particle)
-    @unpack cache = system.boundary_model
+@inline function particle_density(v, ::AdamiPressureExtrapolation, boundary_model, particle)
+    @unpack cache = boundary_model
 
     return cache.density[particle]
 end
 
-@inline function update!(boundary_model::BoundaryModelDummyParticles,
-                         system, system_index, v, u, v_ode, u_ode, semi)
-    @unpack pressure, density_calculator = boundary_model
-    @unpack systems, neighborhood_searches = semi
-
-    pressure .= zero(eltype(pressure))
-
-    compute_quantities!(boundary_model, density_calculator,
-                        system, system_index, v, u, v_ode, u_ode, semi)
-
-    return boundary_model
-end
-
-function compute_quantities!(boundary_model, ::SummationDensity,
-                             system, system_index, v, u, v_ode, u_ode, semi)
-    @unpack systems, neighborhood_searches = semi
-    @unpack state_equation, pressure, cache = boundary_model
-    @unpack density = cache # Density is in the cache for SummationDensity
-
-    density .= zero(eltype(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)
-
-        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,
-                              neighborhood_search;
-                              particles=eachparticle(system)) do particle, neighbor,
-                                                                 pos_diff, distance
-            mass = hydrodynamic_mass(neighbor_system, neighbor)
-            density[particle] += mass * smoothing_kernel(boundary_model, distance)
-        end
-    end
-
-    for particle in eachparticle(system)
-        pressure[particle] = state_equation(particle_density(v, boundary_model, particle))
-    end
-end
-
-function compute_quantities!(boundary_model, ::ContinuityDensity,
-                             system, system_index, v, u, v_ode, u_ode, semi)
-    @unpack systems, neighborhood_searches = semi
-    @unpack pressure, state_equation = boundary_model
+@inline function update!(boundary_model::BoundaryModelDummyParticles, system, system_index,
+                         v, u, v_ode, u_ode, semi)
+    @unpack density_calculator = boundary_model
 
-    for particle in eachparticle(system)
-        pressure[particle] = state_equation(particle_density(v, boundary_model, particle))
-    end
+    update!(boundary_model, density_calculator, system, system_index, v, u,
+            v_ode, u_ode, semi)
 end
 
-function compute_quantities!(boundary_model, ::AdamiPressureExtrapolation,
-                             system, system_index, v, u, v_ode, u_ode, semi)
+@inline function update!(boundary_model, density_calculator::AdamiPressureExtrapolation,
+                         system, system_index, v, u, v_ode, u_ode, semi)
     @unpack systems, neighborhood_searches = semi
     @unpack pressure, state_equation, cache = boundary_model
-    @unpack density, volume = cache
+    @unpack volume, density = cache
 
+    pressure .= zero(eltype(pressure))
     density .= zero(eltype(density))
     volume .= zero(eltype(volume))
 
@@ -283,9 +234,29 @@ end
     end
 end
 
-@inline function adami_pressure_extrapolation!(boundary_model, system,
-                                               neighbor_system,
+@inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system,
                                                system_coords, neighbor_coords,
                                                v_neighbor_system, neighborhood_search)
     return boundary_model
 end
+
+@inline function update!(boundary_model, density_calculator::SummationDensity,
+                         system, system_index, v, u, v_ode, u_ode, semi)
+    @unpack cache, state_equation, pressure = boundary_model
+
+    summation_density!(system, system_index, semi, u, u_ode, cache.density,
+                       particles=eachparticle(system))
+
+    for particle in eachparticle(system)
+        pressure[particle] = state_equation(particle_density(v, boundary_model, particle))
+    end
+end
+
+@inline function update!(boundary_model, density_calculator::ContinuityDensity,
+                         system, system_index, v, u, v_ode, u_ode, semi)
+    @unpack state_equation, pressure = boundary_model
+
+    for particle in eachparticle(system)
+        pressure[particle] = state_equation(particle_density(v, boundary_model, particle))
+    end
+end

src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -67,15 +67,6 @@ function Base.show(io::IO, model::BoundaryModelMonaghanKajtar)
     print(io, "BoundaryModelMonaghanKajtar")
 end
 
-@inline function particle_density(v, model::BoundaryModelMonaghanKajtar, system,
-                                  particle)
-    @unpack hydrodynamic_mass, boundary_particle_spacing = model
-
-    # This model does not use any particle density. However, a mean density is used for
-    # `ArtificialViscosityMonaghan` in the fluid interaction.
-    return hydrodynamic_mass[particle] / boundary_particle_spacing^ndims(system)
-end
-
 @inline function boundary_particle_impact(particle, boundary_particle,
                                           boundary_model::BoundaryModelMonaghanKajtar,
                                           v_particle_system, v_boundary_system,
@@ -102,8 +93,16 @@ end
     return 1.77 / 32 * (1 + 5 / 2 * q + 2 * q^2) * (2 - q)^5
 end
 
-@inline function update!(boundary_model::BoundaryModelMonaghanKajtar,
-                         system, system_index, v, u, v_ode, u_ode, semi)
+@inline function particle_density(v, model::BoundaryModelMonaghanKajtar, system, particle)
+    @unpack hydrodynamic_mass, boundary_particle_spacing = model
+
+    # This model does not use any particle density. However, a mean density is used for
+    # `ArtificialViscosityMonaghan` in the fluid interaction.
+    return hydrodynamic_mass[particle] / boundary_particle_spacing^ndims(system)
+end
+
+@inline function update!(boundary_model::BoundaryModelMonaghanKajtar, system, system_index,
+                         v, u, v_ode, u_ode, semi)
     # Nothing to do in the update step
     return boundary_model
 end

src/schemes/boundary/system.jl

@@ -139,6 +139,11 @@ end
     return system.boundary_model.hydrodynamic_mass[particle]
 end
 
+@inline function smoothing_kernel(system::BoundarySPHSystem, distance)
+    @unpack smoothing_kernel, smoothing_length = system.boundary_model
+    return kernel(smoothing_kernel, distance, smoothing_length)
+end
+
 # This update depends on the computed quantities of the fluid system and therefore
 # has to be in `update_final!` after `update_quantities!`.
 function update_final!(system::BoundarySPHSystem, system_index, v, u, v_ode, u_ode, semi, t)

src/schemes/fluid/weakly_compressible_sph/system.jl

@@ -122,45 +122,21 @@ function update_quantities!(system::WeaklyCompressibleSPHSystem, system_index, v
                             v_ode, u_ode, semi, t)
     @unpack density_calculator = system
 
-    compute_quantities(v, u, density_calculator, system, system_index, u_ode, semi)
+    compute_density!(system, system_index, semi, u, u_ode, density_calculator)
 
-    return system
-end
-
-function compute_quantities(v, u, ::ContinuityDensity, system, system_index, u_ode, semi)
     compute_pressure!(system, v)
 end
 
-function compute_quantities(v, u, ::SummationDensity, system, system_index, u_ode, semi)
-    @unpack systems, neighborhood_searches = semi
+function compute_density!(system::WeaklyCompressibleSPHSystem, system_index, semi, u, u_ode,
+                          ::SummationDensity)
     @unpack cache = system
     @unpack density = cache # Density is in the cache for SummationDensity
 
-    density .= zero(eltype(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)
-
-        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,
-                              neighborhood_search) do particle, neighbor, pos_diff, distance
-            mass = hydrodynamic_mass(neighbor_system, neighbor)
-            density[particle] += mass * smoothing_kernel(system, distance)
-        end
-    end
-
-    compute_pressure!(system, v)
+    summation_density!(system, system_index, semi, u, u_ode, density)
 end
 
+compute_density!(system, system_index, semi, u, u_ode, ::ContinuityDensity) = system
+
 function compute_pressure!(system, v)
     @unpack state_equation, pressure = system