WCSPH Pressure eq consistency

src/callbacks/info.jl

@@ -115,7 +115,7 @@ function initialize_info_callback(discrete_callback, u, t, integrator;
     show(io_context, MIME"text/plain"(), semi)
     println(io, "\n")
     systems = semi.systems
-    for system in systems
+    foreach_enumerate(systems) do (system_index, system)
         show(io_context, MIME"text/plain"(), system)
         println(io, "\n")
     end

src/general/semidiscretization.jl

@@ -31,6 +31,10 @@ struct Semidiscretization{S, RU, RV, NS, DC}
         ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i])
                          for i in eachindex(sizes_v))
 
+        # Check that the boundary systems are using a state equation if EDAC is not used.
+        # Other checks might be added here later.
+        check_configuration(systems)
+
         # Create (and initialize) a tuple of n neighborhood searches for each of the n systems
         # We will need one neighborhood search for each pair of systems.
         searches = Tuple(Tuple(create_neighborhood_search(system, neighbor,
@@ -479,3 +483,23 @@ function nhs_coords(system::BoundarySPHSystem,
     # Don't update
     return nothing
 end
+
+function check_configuration(systems)
+    foreach_enumerate(systems) do (system_index, system)
+        check_configuration(system, systems)
+    end
+end
+
+check_configuration(system, systems) = nothing
+
+function check_configuration(boundary_system::BoundarySPHSystem, systems)
+    (; boundary_model) = boundary_system
+
+    foreach_enumerate(systems) do (neighbor_index, neighbor)
+        if neighbor isa WeaklyCompressibleSPHSystem &&
+           boundary_model isa BoundaryModelDummyParticles &&
+           isnothing(boundary_model.state_equation)
+            throw(ArgumentError("`WeaklyCompressibleSPHSystem` cannot be used without setting a `state_equation` for all boundary systems"))
+        end
+    end
+end

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -108,7 +108,8 @@ end
                                        boundary_system,
                                        v_boundary_system,
                                        boundary_model::BoundaryModelDummyParticles, rho_a,
-                                       rho_b, pos_diff, distance, grad_kernel)
+                                       rho_b, pos_diff, distance, grad_kernel,
+                                       density_calculator)
     rho_a = particle_density(v_particle_system, particle_system, particle)
     rho_b = particle_density(v_boundary_system, boundary_system, boundary_particle)
 

src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -76,7 +76,8 @@ end
                                        v_particle_system, neighbor, neighbor_system,
                                        v_neighbor_system,
                                        boundary_model::BoundaryModelMonaghanKajtar, rho_a,
-                                       rho_b, pos_diff, distance, grad_kernel)
+                                       rho_b, pos_diff, distance, grad_kernel,
+                                       density_calculator)
     (; smoothing_length) = particle_system
     (; K, beta, boundary_particle_spacing) = boundary_model
 

src/schemes/fluid/weakly_compressible_sph/rhs.jl

@@ -41,8 +41,8 @@ function interact!(dv, v_particle_system, u_particle_system,
         dv_pressure = pressure_acceleration(pressure_correction, m_b, particle,
                                             particle_system, v_particle_system,
                                             neighbor, neighbor_system,
-                                            v_neighbor_system, rho_a, rho_b,
-                                            pos_diff, distance, grad_kernel)
+                                            v_neighbor_system, rho_a, rho_b, pos_diff,
+                                            distance, grad_kernel, density_calculator)
 
         dv_viscosity = viscosity_correction *
                        viscosity(particle_system, neighbor_system,
@@ -71,11 +71,28 @@ function interact!(dv, v_particle_system, u_particle_system,
     return dv
 end
 
+# As shown in "Variational and momentum preservation aspects of Smooth
+# Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation
+# this form has to be used with ContinuityDensity.
 @inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
                                        v_particle_system, neighbor,
                                        neighbor_system::WeaklyCompressibleSPHSystem,
                                        v_neighbor_system, rho_a, rho_b, pos_diff, distance,
-                                       grad_kernel)
+                                       grad_kernel, ::ContinuityDensity)
+    return (-m_b *
+            (particle_system.pressure[particle] + neighbor_system.pressure[neighbor]) /
+            (rho_a * rho_b) * grad_kernel) *
+           pressure_correction
+end
+
+# As shown in "Variational and momentum preservation aspects of Smooth
+# Particle Hydrodynamic formulations" by Bonet and Lok, 1999 for a consistent formulation
+# this form has to be used with SummationDensity.
+@inline function pressure_acceleration(pressure_correction, m_b, particle, particle_system,
+                                       v_particle_system, neighbor,
+                                       neighbor_system::WeaklyCompressibleSPHSystem,
+                                       v_neighbor_system, rho_a, rho_b, pos_diff, distance,
+                                       grad_kernel, ::SummationDensity)
     return (-m_b *
             (particle_system.pressure[particle] / rho_a^2 +
              neighbor_system.pressure[neighbor] / rho_b^2) * grad_kernel) *
@@ -87,13 +104,13 @@ end
                                        neighbor_system::Union{BoundarySPHSystem,
                                                               TotalLagrangianSPHSystem},
                                        v_neighbor_system, rho_a, rho_b, pos_diff, distance,
-                                       grad_kernel)
+                                       grad_kernel, density_calculator)
     (; boundary_model) = neighbor_system
 
     return pressure_acceleration(pressure_correction, m_b, particle, particle_system,
                                  v_particle_system, neighbor, neighbor_system,
                                  v_neighbor_system, boundary_model, rho_a, rho_b, pos_diff,
-                                 distance, grad_kernel)
+                                 distance, grad_kernel, density_calculator)
 end
 
 @inline function continuity_equation!(dv, density_calculator::ContinuityDensity,

src/schemes/solid/total_lagrangian_sph/rhs.jl

@@ -61,7 +61,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                    particle_system::TotalLagrangianSPHSystem,
                    neighbor_system::WeaklyCompressibleSPHSystem)
     (; boundary_model) = particle_system
-    (; state_equation, viscosity) = neighbor_system
+    (; density_calculator, state_equation, viscosity) = neighbor_system
     (; sound_speed) = state_equation
 
     system_coords = current_coordinates(u_particle_system, particle_system)
@@ -103,7 +103,7 @@ function interact!(dv, v_particle_system, u_particle_system,
                                             particle_system, v_particle_system,
                                             boundary_model, rho_a,
                                             rho_b, pos_diff, distance,
-                                            grad_kernel)
+                                            grad_kernel, density_calculator)
         dv_particle = dv_boundary + dv_viscosity
 
         for i in 1:ndims(particle_system)