Validation hydrostatic water column

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
-version = "0.4.4-dev"
 authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"]
+version = "0.4.4-dev"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

examples/fsi/hydrostatic_water_column_2d.jl

@@ -0,0 +1,160 @@
+# ==========================================================================================
+# 2D Hydrostatic Water Column on an Elastic Plate (FSI)
+#
+# Case "Elastic plate under a hydrostatic water column" as described in
+# "A fluid–structure interaction model for free-surface flows and flexible structures
+# using smoothed particle hydrodynamics on a GPU" by J. O'Connor and B.D. Rogers
+# published in Journal of Fluids and Structures
+# https://doi.org/10.1016/j.jfluidstructs.2021.103312
+# ==========================================================================================
+
+using TrixiParticles
+using OrdinaryDiffEq
+using JSON
+
+# ==========================================================================================
+# ==== Options
+use_edac = true  # Use EDAC or WCSPH
+
+# ==========================================================================================
+# ==== Experiment Parameters
+# The paper uses 5 to 40 (3 for short runtime in CI)
+n_particles_plate_y = 5
+boundary_layers = 3
+spacing_ratio = 1
+gravity = 9.81
+# Recommended: run at least to 0.5 (paper runs to 1.0); using 0.3 for speed in CI
+tspan = (0.0, 0.5)
+
+initial_fluid_size = (1.0, 2.0)
+plate_size = (1.0, 0.05)
+
+fluid_density = 1000.0
+structure_density = 2700.0
+
+E = 67.5e9
+nu = 0.3
+sound_speed = 50.0
+
+# Particle spacings (structure and fluid share same spacing)
+structure_particle_spacing = plate_size[2] / (n_particles_plate_y - 1)
+fluid_particle_spacing = structure_particle_spacing
+boundary_particle_spacing = fluid_particle_spacing / spacing_ratio
+
+# Analytical solution (constant)
+D = E * plate_size[2]^3 / (12 * (1 - nu^2))
+analytical_value = -0.0026 * gravity *
+                   (fluid_density * initial_fluid_size[2] +
+                    structure_density * plate_size[2]) / D
+
+# ==========================================================================================
+# ==== Geometry Definitions: Tank, Beam, and Fixed Particles
+n_particles_plate_x = round(Int, plate_size[1] / structure_particle_spacing + 1)
+n_particles_per_dimension = (n_particles_plate_x, n_particles_plate_y)
+
+plate = RectangularShape(structure_particle_spacing, n_particles_per_dimension,
+                         (0.0, -plate_size[2]), density=structure_density,
+                         place_on_shell=true)
+
+left_wall = RectangularShape(structure_particle_spacing, (3, n_particles_plate_y),
+                             (-3 * structure_particle_spacing, -plate_size[2]),
+                             density=structure_density, place_on_shell=true)
+right_wall = RectangularShape(structure_particle_spacing, (3, n_particles_plate_y),
+                              (plate_size[1] + structure_particle_spacing,
+                               -plate_size[2]), density=structure_density,
+                              place_on_shell=true)
+fixed_particles = union(left_wall, right_wall)
+
+structure_geometry = union(plate, fixed_particles)
+
+# ==========================================================================================
+# ==== Smoothing Kernel, Boundary, and Related Quantities
+smoothing_kernel = WendlandC2Kernel{2}()
+smoothing_length_structure = sqrt(2) * structure_particle_spacing
+
+smoothing_length_fluid = sqrt(2) * fluid_particle_spacing
+
+hydrodynamic_densities = fluid_density * ones(size(structure_geometry.density))
+hydrodynamic_masses = hydrodynamic_densities *
+                      structure_particle_spacing^ndims(structure_geometry)
+boundary_density_calculator = AdamiPressureExtrapolation()
+
+# ==========================================================================================
+# ==== Run Simulations
+state_equation = use_edac ? nothing :
+                 StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7, clip_negative_pressure=false)
+
+tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, (plate_size[1], 3.0),
+                       min_coordinates=(0.0, fluid_particle_spacing / 2),
+                       fluid_density, n_layers=boundary_layers,
+                       spacing_ratio=spacing_ratio,
+                       faces=(true, true, false, false),
+                       acceleration=(0.0, -gravity),
+                       state_equation=state_equation)
+
+if use_edac
+    fluid_system = EntropicallyDampedSPHSystem(tank.fluid, smoothing_kernel,
+                                               smoothing_length_fluid, sound_speed,
+                                               acceleration=(0.0, -gravity),
+                                               correction=ShepardKernelCorrection(),
+                                               source_terms=SourceTermDamping(;
+                                                                              damping_coefficient=0.05))
+else
+    fluid_density_calculator = ContinuityDensity()
+    density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
+    # density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
+    fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
+                                               state_equation, smoothing_kernel,
+                                               smoothing_length_fluid,
+                                               density_diffusion=density_diffusion,
+                                               acceleration=(0.0, -gravity),
+                                               source_terms=SourceTermDamping(;
+                                                                              damping_coefficient=0.05))
+end
+
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
+                                             state_equation=state_equation,
+                                             boundary_density_calculator,
+                                             smoothing_kernel, smoothing_length_fluid)
+boundary_system = WallBoundarySystem(tank.boundary, boundary_model)
+boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densities,
+                                                       hydrodynamic_masses,
+                                                       state_equation=state_equation,
+                                                       boundary_density_calculator,
+                                                       smoothing_kernel,
+                                                       smoothing_length_structure)
+structure_system = TotalLagrangianSPHSystem(structure_geometry, smoothing_kernel,
+                                            smoothing_length_structure,
+                                            E, nu, boundary_model=boundary_model_structure,
+                                            n_clamped_particles=nparticles(fixed_particles),
+                                            acceleration=(0.0, -gravity))
+
+min_corner = min.(minimum(structure_geometry.coordinates, dims=2),
+                  minimum(tank.boundary.coordinates, dims=2),
+                  minimum(tank.fluid.coordinates, dims=2)) .- smoothing_length_fluid
+max_corner = max.(maximum(structure_geometry.coordinates, dims=2),
+                  maximum(tank.boundary.coordinates, dims=2),
+                  maximum(tank.fluid.coordinates, dims=2)) .+ smoothing_length_fluid
+
+cell_list = FullGridCellList(; min_corner, max_corner)
+neighborhood_search = GridNeighborhoodSearch{2}(; update_strategy=ParallelUpdate(),
+                                                cell_list)
+semi = Semidiscretization(structure_system, fluid_system, boundary_system,
+                          neighborhood_search=neighborhood_search,
+                          parallelization_backend=PolyesterBackend())
+ode = semidiscretize(semi, tspan)
+
+split_integration = SplitIntegrationCallback(CarpenterKennedy2N54(williamson_condition=false),
+                                             callback=StepsizeCallback(cfl=1.6),
+                                             dt=1.0) # dt will be overwritten by the stepsize cb
+
+# This can be overwritten with `trixi_include`
+extra_callback = nothing
+
+info_callback = InfoCallback(interval=100)
+saving_callback = SolutionSavingCallback(dt=0.1, prefix="")
+callbacks = CallbackSet(info_callback, saving_callback, split_integration, extra_callback)
+
+sol = solve(ode, RDPK3SpFSAL49(), dt=1e-8, reltol=1e-6, save_everystep=false,
+            callback=callbacks)

src/callbacks/callbacks.jl

@@ -63,6 +63,11 @@ function set_callbacks_used!(semi, integrator)
     return semi
 end
 
+function set_callbacks_used!(p::NamedTuple{(:v_ode, :u_ode, :semi, :semi_split)},
+                             integrator)
+    set_callbacks_used!(p.semi_split, integrator)
+end
+
 include("info.jl")
 include("solution_saving.jl")
 include("density_reinit.jl")

src/general/custom_quantities.jl

@@ -18,6 +18,18 @@ function kinetic_energy(system, dv_ode, du_ode, v_ode, u_ode, semi, t)
     end
 end
 
+function kinetic_energy(system::AbstractStructureSystem,
+                        dv_ode, du_ode, v_ode, u_ode, semi, t)
+    v = wrap_v(v_ode, system, semi)
+    mass = system.mass
+    energy = zero(eltype(system))
+
+    return sum(each_active_particle(system)) do particle
+        v_i = current_velocity(v, system, particle)
+        energy += mass[particle] * dot(v_i, v_i) / 2
+    end
+end
+
 function kinetic_energy(system::AbstractBoundarySystem,
                         dv_ode, du_ode, v_ode, u_ode, semi, t)
     return zero(eltype(system))

src/schemes/schemes.jl

@@ -9,7 +9,7 @@ include("structure/discrete_element_method/discrete_element_method.jl")
 # Monaghan-Kajtar repulsive boundary particles require the `WallBoundarySystem`
 # and the `TotalLagrangianSPHSystem`.
 include("boundary/wall_boundary/monaghan_kajtar.jl")
-# Implicit incompressible SPH requires the `BoundarySPHSystem`
+# Implicit incompressible SPH requires the `WallBoundarySystem`
 include("fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl")
 
 # Include rhs for all schemes

test/examples/examples.jl

@@ -221,6 +221,46 @@
                 @test count_rhs_allocations(sol, semi) == 0
             end
         end
+
+        @trixi_testset "fsi/hydrostatic_water_column_2d.jl" begin
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(examples_dir(), "fsi",
+                                                      "hydrostatic_water_column_2d.jl"),
+                                             tspan=(0.0, 0.1), n_particles_plate_y=3) [
+                r"\[ Info: To create the self-interaction neighborhood search.*\n",
+                r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n",
+                r"│   caller = ip:0x0\n",
+                r"└ @ Core :-1\n"
+            ]
+            @test sol.retcode == ReturnCode.Success
+            if VERSION < v"1.12"
+                # Older Julia versions produce allocations because `get_neighborhood_search`
+                # is not type-stable with TLSPH.
+                @test count_rhs_allocations(sol, semi) < 500
+            else
+                @test count_rhs_allocations(sol, semi) == 0
+            end
+        end
+        @trixi_testset "fsi/hydrostatic_water_column_2d.jl with EDAC" begin
+            @trixi_test_nowarn trixi_include(@__MODULE__,
+                                             joinpath(examples_dir(), "fsi",
+                                                      "hydrostatic_water_column_2d.jl"),
+                                             tspan=(0.0, 0.1), n_particles_plate_y=3,
+                                             use_edac=true) [
+                r"\[ Info: To create the self-interaction neighborhood search.*\n",
+                r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n",
+                r"│   caller = ip:0x0\n",
+                r"└ @ Core :-1\n"
+            ]
+            @test sol.retcode == ReturnCode.Success
+            if VERSION < v"1.12"
+                # Older Julia versions produce allocations because `get_neighborhood_search`
+                # is not type-stable with TLSPH.
+                @test count_rhs_allocations(sol, semi) < 500
+            else
+                @test count_rhs_allocations(sol, semi) == 0
+            end
+        end
     end
 
     @testset verbose=true "N-Body" begin

test/systems/tlsph_system.jl

@@ -383,7 +383,8 @@
         for particle in TrixiParticles.eachparticle(system)
             system.deformation_grad[:, :, particle] = [1.0 0.2; 0.2 1.0]
             system.pk1_rho2[:, :,
-                            particle] = [1.0 0.5; 0.5 1.0] / material_densities[particle]^2
+                            particle] = [1.0 0.5; 0.5 1.0] /
+                                        material_densities[particle]^2
         end
 
         von_mises_stress = TrixiParticles.von_mises_stress(system)

test/validation/validation.jl

@@ -92,6 +92,24 @@
         @test length(axs_wcsph[1].scene.plots) >= 2
     end
 
+    @trixi_testset "hydrostatic_water_column_2d" begin
+        @trixi_test_nowarn trixi_include(@__MODULE__,
+                                         joinpath(validation_dir(),
+                                                  "hydrostatic_water_column_2d",
+                                                  "validation.jl"), tspan=(0.0, 0.3),
+                                         n_particles_plate_y=3) [
+            r"┌ Info: The desired tank length in y-direction.*\n",
+            r"└ New tank length in y-direction is set to.*\n",
+            r"┌ Warning: keyword `n_clamped_particles` is deprecated.*\n",
+            r"│   caller = ip:0x0\n",
+            r"└ @ Core :-1\n",
+            r"\[ Info: To create the self-interaction neighborhood search.*\n"
+        ]
+
+        # We compare the relative error to the analytical solution
+        @test isapprox(errors[:edac][2], 0.0, atol=0.033)
+        @test isapprox(errors[:wcsph][2], 0.0, atol=0.045)
+    end
     @trixi_testset "TGV_2D" begin
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(validation_dir(),