rework moving boundaries and support for BoundaryModelDummyParticles

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Morton = "2a6d852e-3fac-5a38-885c-fe708af2d09e"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"

examples/fluid/accelerated_tank_2d.jl

@@ -0,0 +1,91 @@
+# This setup is identical to `rectangular_tank_2d.jl`, except that now there is no gravity, and
+# the tank is accelerated upwards instead.
+# Note that the two setups are physically identical, but produce different numerical errors.
+using TrixiParticles
+using OrdinaryDiffEq
+
+gravity = 0.0
+
+# ==========================================================================================
+# ==== Fluid
+
+particle_spacing = 0.02
+
+# Ratio of fluid particle spacing to boundary particle spacing
+beta = 1
+boundary_layers = 3
+
+water_width = 2.0
+water_height = 0.9
+water_density = 1000.0
+
+tank_width = 2.0
+tank_height = 1.0
+
+sound_speed = 10 * sqrt(9.81 * water_height)
+state_equation = StateEquationCole(sound_speed, 7, water_density, 100000.0,
+                                   background_pressure=100000.0)
+
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+
+viscosity = ArtificialViscosityMonaghan(0.02, 0.0)
+
+tank = RectangularTank(particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta)
+
+# ==========================================================================================
+# ==== Boundary models
+
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density,
+                                             tank.boundary.mass, state_equation,
+                                             AdamiPressureExtrapolation(),
+                                             smoothing_kernel, smoothing_length)
+
+f_y(t) = 0.5 * 9.81 * t^2
+f_x(t) = 0.0
+
+is_moving(t) = true
+
+movement = BoundaryMovement((f_x, f_y), is_moving)
+
+# ==========================================================================================
+# ==== Systems
+
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, SummationDensity(), state_equation,
+                                           smoothing_kernel, smoothing_length,
+                                           viscosity=viscosity,
+                                           acceleration=(0.0, gravity))
+
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model,
+                                    movement=movement)
+
+# ==========================================================================================
+# ==== Simulation
+
+semi = Semidiscretization(fluid_system, boundary_system,
+                          #damping_coefficient=1e-5,
+                          neighborhood_search=SpatialHashingSearch)
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=10)
+saving_callback = SolutionSavingCallback(dt=0.02)
+
+callbacks = CallbackSet(info_callback, saving_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control.
+# Enable threading of the RK method for better performance on multiple threads.
+# Limiting of the maximum stepsize is necessary to prevent crashing.
+# When particles are approaching a wall in a uniform way, they can be advanced
+# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
+# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces
+# become extremely large when fluid particles are very close to boundary particles,
+# and the time integration method interprets this as an instability.
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

examples/fluid/dam_break_2d.jl

@@ -68,7 +68,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/dam_break_2d_corrections.jl

@@ -60,7 +60,7 @@ density_calculator_dict = Dict("no_correction" => SummationDensity(),
                                "kernel_gradient_summation_correction" => SummationDensity(),
                                "kernel_gradient_continuity_correction" => ContinuityDensity())
 
-boundary_system = BoundarySPHSystem(setup.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(setup.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/dam_break_3d.jl

@@ -61,7 +61,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity, 0.0))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/falling_water_column_2d.jl

@@ -57,7 +57,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

examples/fluid/moving_wall_2d.jl

@@ -0,0 +1,108 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+gravity = -9.81
+
+# ==========================================================================================
+# ==== Fluid
+
+particle_spacing = 0.03
+# Ratio of fluid particle spacing to boundary particle spacing
+beta_wall = 1
+beta_tank = 1
+boundary_layers = 3
+boundary_layers_wall = 3
+
+water_width = 1.0
+water_height = 0.8
+water_density = 1000.0
+
+tank_width = 4.0
+tank_height = 1.0
+
+sound_speed = 10 * sqrt(9.81 * water_height)
+
+state_equation = StateEquationCole(sound_speed, 7, water_density, 100_000.0,
+                                   background_pressure=100_000.0)
+
+viscosity = ArtificialViscosityMonaghan(0.1, 0.0)
+
+smoothing_length = 1.2 * particle_spacing
+smoothing_kernel = SchoenbergCubicSplineKernel{2}()
+
+tank = RectangularTank(particle_spacing, (water_width, water_height),
+                       (tank_width, tank_height), water_density,
+                       n_layers=boundary_layers, spacing_ratio=beta_tank)
+
+# ==========================================================================================
+# ==== Boundary
+
+boundary_particle_spacing = particle_spacing / beta_wall
+
+# Move right boundary
+wall_position = (tank.n_particles_per_dimension[1] + 1) * particle_spacing
+n_wall_particles_y = size(tank.face_indices[2], 2) * beta_wall
+
+wall = RectangularShape(boundary_particle_spacing,
+                        (boundary_layers_wall, n_wall_particles_y),
+                        (wall_position, boundary_particle_spacing), water_density)
+
+f_y(t) = 0.0
+f_x(t) = 0.5t^2
+
+is_moving(t) = t < 1.5
+
+movement = BoundaryMovement((f_x, f_y), is_moving)
+
+# ==========================================================================================
+# ==== Boundary models
+
+boundary_model_tank = BoundaryModelDummyParticles(tank.boundary.density,
+                                                  tank.boundary.mass, state_equation,
+                                                  AdamiPressureExtrapolation(),
+                                                  smoothing_kernel, smoothing_length)
+
+boundary_model_wall = BoundaryModelDummyParticles(wall.density, wall.mass, state_equation,
+                                                  AdamiPressureExtrapolation(),
+                                                  smoothing_kernel, smoothing_length)
+
+# ==========================================================================================
+# ==== Systems
+
+fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, SummationDensity(), state_equation,
+                                           smoothing_kernel, smoothing_length,
+                                           viscosity=viscosity, acceleration=(0.0, gravity))
+
+boundary_system_tank = BoundarySPHSystem(tank.boundary, boundary_model_tank)
+
+boundary_system_wall = BoundarySPHSystem(wall, boundary_model_wall,
+                                         movement=movement)
+
+# ==========================================================================================
+# ==== Simulation
+
+tspan = (0.0, 2.0)
+
+semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_wall,
+                          neighborhood_search=SpatialHashingSearch)
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+saving_callback = SolutionSavingCallback(dt=0.02)
+callbacks = CallbackSet(info_callback, saving_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control.
+# Enable threading of the RK method for better performance on multiple threads.
+# Limiting of the maximum stepsize is necessary to prevent crashing.
+# When particles are approaching a wall in a uniform way, they can be advanced
+# with large time steps. Close to the wall, the stepsize has to be reduced drastically.
+# Sometimes, the method fails to do so with Monaghan-Kajtar BC because forces
+# become extremely large when fluid particles are very close to boundary particles,
+# and the time integration method interprets this as an instability.
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

examples/fluid/rectangular_tank_2d.jl

@@ -53,7 +53,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation

examples/fsi/bending_beam_2d.jl

@@ -130,7 +130,7 @@ fluid_system = WeaklyCompressibleSPHSystem(fluid,
                                            acceleration=(0.0, gravity))
 
 # Replace the elastic solid by a boundary for testing purposes.
-# rigid_solid_system = BoundarySPHSystem(solid.coordinates, boundary_model_solid)
+# rigid_solid_system = BoundarySPHSystem(solid, boundary_model_solid)
 
 # ==========================================================================================
 # ==== Simulation

examples/fsi/dam_break_2d.jl

@@ -116,7 +116,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 solid_system = TotalLagrangianSPHSystem(solid,
                                         solid_smoothing_kernel, solid_smoothing_length,

examples/fsi/dam_break_gate_2d.jl

@@ -50,11 +50,13 @@ gate_position = (tank.n_particles_per_dimension[1] + 1) * fluid_particle_spacing
 gate = RectangularShape(fluid_particle_spacing / beta_gate,
                         (gate_layers,
                          round(Int, gate_height / fluid_particle_spacing * beta_gate)),
-                        (gate_position, fluid_particle_spacing / beta_gate),
-                        water_density)
+                        (gate_position, fluid_particle_spacing / beta_gate), water_density)
 
-# No moving boundaries for the relaxing step
-movement_function(coordinates, t) = false
+f_x(t) = 0.0
+f_y(t) = -285.115t^3 + 72.305t^2 + 0.1463t
+is_moving(t) = false # No moving boundaries for the relaxing step
+
+movement = BoundaryMovement((f_x, f_y), gate.coordinates, is_moving)
 
 # ==========================================================================================
 # ==== Solid
@@ -127,15 +129,12 @@ boundary_model_solid = BoundaryModelDummyParticles(hydrodynamic_densites,
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), state_equation,
                                            smoothing_kernel, smoothing_length,
-                                           viscosity=viscosity,
-                                           acceleration=(0.0, gravity))
+                                           viscosity=viscosity, acceleration=(0.0, gravity))
 
-boundary_system_tank = BoundarySPHSystem(tank.boundary.coordinates,
-                                         boundary_model_tank)
+boundary_system_tank = BoundarySPHSystem(tank.boundary, boundary_model_tank)
 
-boundary_system_gate = BoundarySPHSystem(gate.coordinates,
-                                         boundary_model_gate,
-                                         movement_function=movement_function)
+boundary_system_gate = BoundarySPHSystem(gate, boundary_model_gate,
+                                         movement=movement)
 
 solid_system = TotalLagrangianSPHSystem(solid,
                                         solid_smoothing_kernel, solid_smoothing_length,
@@ -173,20 +172,7 @@ sol = solve(ode, RDPK3SpFSAL49(),
 # Run full simulation
 tspan = (0.0, 1.0)
 
-function movement_function(coordinates, t)
-    if t < 0.1
-        particle_spacing = coordinates[2, 2] - coordinates[2, 1]
-        f(t) = -285.115t^3 + 72.305t^2 + 0.1463t + particle_spacing
-        pos_1 = coordinates[2, 1]
-        pos_2 = f(t)
-        diff_pos = pos_2 - pos_1
-        coordinates[2, :] .+= diff_pos
-
-        return true
-    end
-
-    return false
-end
+is_moving(t) = t < 0.1
 
 # Use solution of the relaxing step as initial coordinates
 restart_with!(semi, sol)

examples/fsi/falling_spheres_2d.jl

@@ -96,7 +96,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, ContinuityDensity(), stat
                                            viscosity=viscosity,
                                            acceleration=(0.0, gravity))
 
-boundary_system = BoundarySPHSystem(tank.boundary.coordinates, boundary_model)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 solid_system_1 = TotalLagrangianSPHSystem(sphere_1,
                                           solid_smoothing_kernel, solid_smoothing_length,

src/TrixiParticles.jl

@@ -17,6 +17,7 @@ using ThreadingUtilities
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
 @reexport using SimpleUnPack: @unpack
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes
+using ForwardDiff
 
 # util needs to be first because of macro @trixi_timeit
 include("util.jl")
@@ -40,6 +41,7 @@ export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
 export StateEquationIdealGas, StateEquationCole
 export ArtificialViscosityMonaghan, ViscosityAdami
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation
+export BoundaryMovement
 export SpatialHashingSearch
 export examples_dir, trixi_include
 export trixi2vtk

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -296,17 +296,19 @@ end
                           neighborhood_search;
                           particles=eachparticle(system)) do particle, neighbor,
                                                              pos_diff, distance
-        density_neighbor = particle_density(v_neighbor_system, neighbor_system,
-                                            neighbor)
+        density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor)
+
+        resulting_acc = neighbor_system.acceleration -
+                        current_acceleration(system, particle)
 
         kernel_weight = smoothing_kernel(boundary_model, distance)
 
-        # TODO moving boundaries
         pressure[particle] += (neighbor_system.pressure[neighbor] +
-                               dot(neighbor_system.acceleration,
-                                   density_neighbor * pos_diff)) * kernel_weight
+                               dot(resulting_acc, density_neighbor * pos_diff)) *
+                              kernel_weight
 
         cache.volume[particle] += kernel_weight
+
         compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system,
                                    kernel_weight, particle, neighbor)
     end