diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl index 751f60be0..43dadce78 100644 --- a/examples/fsi/dam_break_2d.jl +++ b/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)) # ========================================================================================== diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index aafa59083..bf5b73f39 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/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)) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 1c8f41bef..fcd0a92e6 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/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)) # ========================================================================================== diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl index cf5c9c4f5..037996f37 100644 --- a/examples/fsi/falling_water_column_2d.jl +++ b/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)) # ========================================================================================== diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl index fe71f1532..3c7fe4f64 100644 --- a/examples/solid/oscillating_beam_2d.jl +++ b/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 diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index d91d0f8b2..7b7a77dbe 100644 --- a/src/general/semidiscretization.jl +++ b/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 diff --git a/src/general/system.jl b/src/general/system.jl index a4cbee254..06c370f30 100644 --- a/src/general/system.jl +++ b/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 diff --git a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl b/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl index 0ea7970d5..85565e582 100644 --- a/src/schemes/boundary/monaghan_kajtar/monaghan_kajtar.jl +++ b/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()) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index df2d17cca..484413f36 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/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