The grand reopening Part 1: Validation setup for dam break

examples/fluid/dam_break_2d.jl

@@ -8,49 +8,57 @@
 using TrixiParticles
 using OrdinaryDiffEq
 
+# Size parameters
+H = 0.6
+W = 2 * H
+
 # ==========================================================================================
 # ==== Resolution
-fluid_particle_spacing = 0.02
+fluid_particle_spacing = H / 40
 
 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model
-boundary_layers = 3
+boundary_layers = 4
 spacing_ratio = 1
 
 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio
 
 # ==========================================================================================
 # ==== Experiment Setup
 gravity = 9.81
+
 tspan = (0.0, 5.7 / sqrt(gravity))
 
 # Boundary geometry and initial fluid particle positions
-initial_fluid_size = (2.0, 1.0)
-tank_size = (floor(5.366 / boundary_particle_spacing) * boundary_particle_spacing, 4.0)
+initial_fluid_size = (W, H)
+tank_size = (floor(5.366 * H / boundary_particle_spacing) * boundary_particle_spacing, 4.0)
 
 fluid_density = 1000.0
-sound_speed = 20 * sqrt(gravity * initial_fluid_size[2])
+sound_speed = 20 * sqrt(gravity * H)
 state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
-                                   exponent=7, clip_negative_pressure=false)
+                                   exponent=1, clip_negative_pressure=false)
 
 tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                        n_layers=boundary_layers, spacing_ratio=spacing_ratio,
                        acceleration=(0.0, -gravity), state_equation=state_equation)
 
 # ==========================================================================================
 # ==== Fluid
-smoothing_length = 3.0 * fluid_particle_spacing
+smoothing_length = 3.5 * fluid_particle_spacing
 smoothing_kernel = WendlandC2Kernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
 viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0)
-density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1)
-# density_diffusion = DensityDiffusionAntuono(tank.fluid, delta=0.1)
+# Alternatively the density diffusion model by Molteni & Colagrossi can be used,
+# which will run faster.
+# 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, viscosity=viscosity,
                                            density_diffusion=density_diffusion,
-                                           acceleration=(0.0, -gravity), correction=nothing)
+                                           acceleration=(0.0, -gravity),
+                                           correction=nothing)
 
 # ==========================================================================================
 # ==== Boundary
@@ -65,18 +73,32 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system)
+semi = Semidiscretization(fluid_system, boundary_system, threaded_nhs_update=true)
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=100)
-saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+solution_prefix = ""
+saving_callback = SolutionSavingCallback(dt=0.02, prefix=solution_prefix)
+
+# Save at certain timepoints which allows comparison to the results of Marrone et al.,
+# i.e. (1.5, 2.36, 3.0, 5.7, 6.45).
+# Please note that the images in Marrone et al. are obtained at a particle_spacing = H/320,
+# which takes between 2 and 4 hours.
+saving_paper = SolutionSavingCallback(save_times=[0.0, 0.371, 0.584, 0.743, 1.411, 1.597],
+                                      prefix="marrone_times")
+
+# This can be overwritten with `trixi_include`
+extra_callback = nothing
 
 use_reinit = false
-density_reinit_cb = use_reinit ? DensityReinitializationCallback(semi.systems[1], dt=0.01) :
+density_reinit_cb = use_reinit ?
+                    DensityReinitializationCallback(semi.systems[1], interval=10) :
                     nothing
-stepsize_callback = StepsizeCallback(cfl=1.1)
+stepsize_callback = StepsizeCallback(cfl=0.9)
 
-callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback)
+callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback,
+                        density_reinit_cb, saving_paper)
 
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
             dt=1.0, # This is overwritten by the stepsize callback

src/TrixiParticles.jl

@@ -5,7 +5,7 @@ using Reexport: @reexport
 using CSV: CSV
 using Dates
 using DataFrames: DataFrame
-using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
+using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect, PresetTimeCallback
 using FastPow: @fastpow
 using ForwardDiff: ForwardDiff
 using JSON: JSON

src/callbacks/solution_saving.jl

@@ -1,8 +1,9 @@
 @doc raw"""
-    SolutionSavingCallback(; interval::Integer=0, dt=0.0,
+    SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_times=Array{Float64, 1}([]),
                            save_initial_solution=true, save_final_solution=true,
-                           output_directory="out", append_timestamp=false,
-                           max_coordinates=2^15, custom_quantities...)
+                           output_directory="out", append_timestamp=false, max_coordinates=2^15,
+                           custom_quantities...)
+
 
 Callback to save the current numerical solution in VTK format in regular intervals.
 Either pass `interval` to save every `interval` time steps,
@@ -19,6 +20,7 @@ To ignore a custom quantity for a specific system, return `nothing`.
 - `dt`:                         Save the solution in regular intervals of `dt` in terms
                                 of integration time by adding additional `tstops`
                                 (note that this may change the solution).
+- `save_times=[]`               List of times at which to save a solution.
 - `save_initial_solution=true`: Save the initial solution.
 - `save_final_solution=true`:   Save the final solution.
 - `output_directory="out"`:     Directory to save the VTK files.
@@ -64,6 +66,7 @@ saving_callback = SolutionSavingCallback(dt=0.1, my_custom_quantity=kinetic_ener
 """
 mutable struct SolutionSavingCallback{I, CQ}
     interval              :: I
+    save_times            :: Array{Float64, 1}
     save_initial_solution :: Bool
     save_final_solution   :: Bool
     write_meta_data       :: Bool
@@ -76,13 +79,14 @@ mutable struct SolutionSavingCallback{I, CQ}
 end
 
 function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
-                                save_initial_solution=true,
-                                save_final_solution=true,
+                                save_times=Array{Float64, 1}([]),
+                                save_initial_solution=true, save_final_solution=true,
                                 output_directory="out", append_timestamp=false,
                                 prefix="", verbose=false, write_meta_data=true,
                                 max_coordinates=Float64(2^15), custom_quantities...)
-    if dt > 0 && interval > 0
-        throw(ArgumentError("Setting both interval and dt is not supported!"))
+    if (dt > 0 && interval > 0) || (length(save_times) > 0 && (dt > 0 || interval > 0))
+        throw(ArgumentError("Setting multiple save times for the same solution " *
+                            "callback is not possible. Use either `dt`, `interval` or `save_times`."))
     end
 
     if dt > 0
@@ -93,13 +97,15 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
         output_directory *= string("_", Dates.format(now(), "YY-mm-ddTHHMMSS"))
     end
 
-    solution_callback = SolutionSavingCallback(interval,
+    solution_callback = SolutionSavingCallback(interval, save_times,
                                                save_initial_solution, save_final_solution,
                                                write_meta_data, verbose, output_directory,
                                                prefix, max_coordinates, custom_quantities,
                                                -1)
 
-    if dt > 0
+    if length(save_times) > 0
+        return PresetTimeCallback(save_times, solution_callback)
+    elseif dt > 0
         # Add a `tstop` every `dt`, and save the final solution.
         return PeriodicCallback(solution_callback, dt,
                                 initialize=initialize_save_cb!,