acceleration for tlsph AdamiPressureExtrapolation

and make the `boundary_model` optional for tlsph

examples/fsi/dam_break_2d.jl

@@ -125,7 +125,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
                                         E, nu,
                                         n_fixed_particles=n_particles_x,
                                         acceleration=(0.0, gravity),
-                                        boundary_model_solid,
+                                        boundary_model=boundary_model_solid,
                                         penalty_force=PenaltyForceGanzenmueller(alpha=0.01))
 
 # ==========================================================================================

examples/fsi/dam_break_gate_2d.jl

@@ -154,7 +154,7 @@ boundary_system_gate = BoundarySPHSystem(gate, boundary_model_gate,
 
 solid_system = TotalLagrangianSPHSystem(solid,
                                         solid_smoothing_kernel, solid_smoothing_length,
-                                        E, nu, boundary_model_solid,
+                                        E, nu, boundary_model=boundary_model_solid,
                                         n_fixed_particles=n_particles_x,
                                         acceleration=(0.0, gravity))
 

examples/fsi/falling_spheres_2d.jl

@@ -106,14 +106,14 @@ solid_system_1 = TotalLagrangianSPHSystem(sphere_1,
                                           solid_smoothing_kernel, solid_smoothing_length,
                                           E1, nu,
                                           acceleration=(0.0, gravity),
-                                          solid_boundary_model_1,
+                                          boundary_model=solid_boundary_model_1,
                                           penalty_force=PenaltyForceGanzenmueller(alpha=0.3))
 
 solid_system_2 = TotalLagrangianSPHSystem(sphere_2,
                                           solid_smoothing_kernel, solid_smoothing_length,
                                           E2, nu,
                                           acceleration=(0.0, gravity),
-                                          solid_boundary_model_2,
+                                          boundary_model=solid_boundary_model_2,
                                           penalty_force=PenaltyForceGanzenmueller(alpha=0.3))
 
 # ==========================================================================================

examples/fsi/falling_water_column_2d.jl

@@ -93,7 +93,8 @@ solid_system = TotalLagrangianSPHSystem(solid,
                                         smoothing_kernel, smoothing_length,
                                         E, nu,
                                         n_fixed_particles=nparticles(fixed_particles),
-                                        acceleration=(0.0, gravity), boundary_model,
+                                        acceleration=(0.0, gravity),
+                                        boundary_model=boundary_model,
                                         penalty_force=PenaltyForceGanzenmueller(alpha=0.1))
 
 # ==========================================================================================

examples/solid/oscillating_beam_2d.jl

@@ -50,8 +50,7 @@ solid_system = TotalLagrangianSPHSystem(solid,
                                         smoothing_kernel, smoothing_length,
                                         E, nu,
                                         n_fixed_particles=nparticles(fixed_particles),
-                                        acceleration=(0.0, acceleration),
-                                        nothing) # No boundary model
+                                        acceleration=(0.0, acceleration))
 
 # ==========================================================================================
 # ==== Simulation

src/general/semidiscretization.jl

@@ -412,6 +412,9 @@ function system_interaction!(dv_ode, v_ode, u_ode, semi)
             @trixi_timeit timer() timer_str begin
                 interact!(dv, v_system, u_system, v_neighbor, u_neighbor,
                           neighborhood_search, system, neighbor)
+                # Only for `TotalLagrangianSPHSystem` with boundary model
+                # `AdamiPressureExtrapolation`.
+                copy_dv!(system, dv)
             end
         end
     end
@@ -492,6 +495,14 @@ end
 
 check_configuration(system, systems) = nothing
 
+function check_configuration(system::TotalLagrangianSPHSystem, systems)
+    foreach_enumerate(systems) do (neighbor_index, neighbor)
+        if !isa(neighbor, TotalLagrangianSPHSystem) && isnothing(system.boundary_model)
+            throw(ArgumentError("Please specify a boundary model for `TotalLagrangianSPHSystem` when simulating a fluid-structure interaction."))
+        end
+    end
+end
+
 function check_configuration(boundary_system::BoundarySPHSystem, systems)
     (; boundary_model) = boundary_system
 

src/general/system.jl

@@ -61,11 +61,6 @@ end
 
 @inline current_velocity(v, system, particle) = extract_svector(v, system, particle)
 
-@inline function current_acceleration(system, particle)
-    # TODO: Return `dv` of solid particles
-    return SVector(ntuple(_ -> 0.0, Val(ndims(system))))
-end
-
 @inline function smoothing_kernel(system, distance)
     (; smoothing_kernel, smoothing_length) = system
     return kernel(smoothing_kernel, distance, smoothing_length)
@@ -101,3 +96,5 @@ end
 function update_final!(system, system_index, v, u, v_ode, u_ode, semi, t)
     return system
 end
+
+copy_dv!(system, dv) = system

src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl

@@ -56,6 +56,7 @@ struct BoundaryModelMonaghanKajtar{ELTYPE <: Real, V}
     boundary_particle_spacing :: ELTYPE
     hydrodynamic_mass         :: Vector{ELTYPE}
     viscosity                 :: V
+    density_calculator        :: Nothing
 
     function BoundaryModelMonaghanKajtar(K, beta, boundary_particle_spacing, mass;
                                          viscosity=NoViscosity())

src/schemes/solid/total_lagrangian_sph/system.jl

@@ -75,7 +75,7 @@ The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGa
   In: International Journal for Numerical Methods in Engineering 48 (2000), pages 1359–1400.
   [doi: 10.1002/1097-0207](https://doi.org/10.1002/1097-0207)
 """
-struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIMS}
+struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF, C} <: System{NDIMS}
     initial_condition   :: InitialCondition{ELTYPE}
     initial_coordinates :: Array{ELTYPE, 2} # [dimension, particle]
     current_coordinates :: Array{ELTYPE, 2} # [dimension, particle]
@@ -94,11 +94,12 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM
     acceleration        :: SVector{NDIMS, ELTYPE}
     boundary_model      :: BM
     penalty_force       :: PF
+    cache               :: C
 
     function TotalLagrangianSPHSystem(initial_condition,
                                       smoothing_kernel, smoothing_length,
-                                      young_modulus, poisson_ratio, boundary_model;
-                                      n_fixed_particles=0,
+                                      young_modulus, poisson_ratio;
+                                      boundary_model=nothing, n_fixed_particles=0,
                                       acceleration=ntuple(_ -> 0.0,
                                                           ndims(smoothing_kernel)),
                                       penalty_force=nothing)
@@ -130,20 +131,33 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: System{NDIM
                       ((1 + poisson_ratio) * (1 - 2 * poisson_ratio))
         lame_mu = 0.5 * young_modulus / (1 + poisson_ratio)
 
-        return new{typeof(boundary_model),
-                   NDIMS, ELTYPE,
-                   typeof(smoothing_kernel),
-                   typeof(penalty_force)}(initial_condition, initial_coordinates,
-                                          current_coordinates, mass,
-                                          correction_matrix, pk1_corrected,
-                                          deformation_grad, material_density,
-                                          n_moving_particles, young_modulus, poisson_ratio,
-                                          lame_lambda, lame_mu,
-                                          smoothing_kernel, smoothing_length,
-                                          acceleration_, boundary_model, penalty_force)
+        cache = create_solid_cache(boundary_model, initial_condition)
+
+        return new{typeof(boundary_model), NDIMS, ELTYPE,
+                   typeof(smoothing_kernel), typeof(penalty_force),
+                   typeof(cache)}(initial_condition, initial_coordinates,
+                                  current_coordinates, mass, correction_matrix,
+                                  pk1_corrected, deformation_grad, material_density,
+                                  n_moving_particles, young_modulus, poisson_ratio,
+                                  lame_lambda, lame_mu, smoothing_kernel, smoothing_length,
+                                  acceleration_, boundary_model, penalty_force, cache)
     end
 end
 
+create_solid_cache(::Nothing, initial_condition) = (;)
+
+function create_solid_cache(boundary_model, initial_condition)
+    create_solid_cache(boundary_model, boundary_model.density_calculator, initial_condition)
+end
+
+create_solid_cache(boundary_model, density_calculator, initial_condition) = (;)
+
+function create_solid_cache(boundary_model, ::AdamiPressureExtrapolation, initial_condition)
+    acceleration = similar(initial_condition.velocity)
+
+    return (; acceleration)
+end
+
 function Base.show(io::IO, system::TotalLagrangianSPHSystem)
     @nospecialize system # reduce precompilation time
 
@@ -225,6 +239,11 @@ end
     return extract_svector(v, system, particle)
 end
 
+# Acceleration for `AdamiPressureExtrapolation`
+@inline function current_acceleration(system::TotalLagrangianSPHSystem, particle)
+    return extract_svector(system.cache.acceleration, system, particle)
+end
+
 @inline function viscous_velocity(v, system::TotalLagrangianSPHSystem, particle)
     return extract_svector(system.boundary_model.cache.wall_velocity, system, particle)
 end
@@ -429,3 +448,23 @@ end
 function viscosity_model(system::TotalLagrangianSPHSystem)
     return system.boundary_model.viscosity
 end
+
+function copy_dv!(system::TotalLagrangianSPHSystem, dv)
+    copy_dv!(system, system.boundary_model, dv)
+end
+
+copy_dv!(system, ::Nothing, dv) = system
+
+copy_dv!(system, model, dv) = copy_dv!(system, model, model.density_calculator, dv)
+
+copy_dv!(system, boundary_model, density_calculator, dv) = system
+
+function copy_dv!(system, boundary_model, ::AdamiPressureExtrapolation, dv)
+    for particle in each_moving_particle(system)
+        for dim in 1:ndims(system)
+            system.cache.acceleration[dim, particle] = dv[dim, particle]
+        end
+    end
+
+    return system
+end