Validate against dam break results from De Courcy (#934)

* fix

* improve plot

* improve plot

* change comment

* update

* update files

* format

* suppress warnings

* update with reference files on 1 thread

* format

* update files again with new Julia 1.11 release

* update error bounds

* review comments

* format

* Validate against dam break results from De Courcy

* Fix

* Update validation files

* Update JSON validation files

* Clean up timer output

* Update dam break validation files

* Explicitly list reference files

* Split plot file

* Update README and validation files

* Remove comment in dam break example file

* Add thread pinning comment

* Delete old CSV files

* Change example tests

* Reformat

* Fix example file

* Fix EDAC validation

* Add more validation JSON files

* Reset resolution in validation file

* Fix viscosity of dam break setups

* Shorten dam break validation test in CI

* Fix symplectic position verlet tests

* Only run validation tests for now

* Use x86 files

* Update x86 validation files again

* Fix tests and add high res files

* Re-enable tests

* Add last file to plot

* Reformat

* Fix tests

* Ignore warnings

* Reformat

* Fix GPU tests

* Fix comment

* Fix GPU tests

* Fix GPU tests

* Update examples/fluid/dam_break_2d.jl

Co-authored-by: Niklas Neher <[email protected]>

* Clarify README

* Fix comment about no-slip BC

---------

Co-authored-by: Sven Berger <[email protected]>
Co-authored-by: Niklas Neher <[email protected]>

Author: 3 people

.gitignore

@@ -12,7 +12,6 @@ docs/src/news.md
 docs/src/tutorials
 *_replaced.md
 run
-*.json
 !.zenodo.json
 *.png
 *.jpg

examples/fluid/dam_break_2d.jl

@@ -33,7 +33,7 @@ boundary_particle_spacing = fluid_particle_spacing / spacing_ratio
 # ==== Experiment Setup
 gravity = 9.81
 
-tspan = (0.0, 5.7 / sqrt(gravity))
+tspan = (0.0, 5.7 / sqrt(gravity / H))
 
 # Boundary geometry and initial fluid particle positions
 initial_fluid_size = (W, H)
@@ -50,29 +50,19 @@ tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fl
 
 # ==========================================================================================
 # ==== Fluid
-smoothing_length = 1.75 * fluid_particle_spacing
+smoothing_length = 2 * fluid_particle_spacing
 smoothing_kernel = WendlandC2Kernel{2}()
 
 fluid_density_calculator = ContinuityDensity()
+
 alpha = 0.02
 viscosity_fluid = ArtificialViscosityMonaghan(alpha=alpha, beta=0.0)
-# A typical formula to convert Artificial viscosity to a
-# kinematic viscosity is provided by Monaghan as
-# nu = alpha * smoothing_length * sound_speed/8
-
-# Alternatively a kinematic viscosity for water can be set
-# nu = 1.0e-6
 
-# This allows the use of a physical viscosity model like:
-# viscosity_fluid = ViscosityAdami(nu=nu)
-# or with additional dissipation through a Smagorinsky model
-# viscosity_fluid = ViscosityAdamiSGS(nu=nu)
-# For more details see the documentation "Viscosity model overview".
-
-# 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)
+# The density diffusion model by Molteni and Colagrossi shows unphysical effects at the
+# free surface in long-running simulations, but is significantly faster than the model
+# by Antuono. This simulation is short enough to use the faster model.
+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,
@@ -86,9 +76,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
 # ==== Boundary
 boundary_density_calculator = AdamiPressureExtrapolation()
 viscosity_wall = nothing
-# For a no-slip condition the corresponding wall viscosity without SGS can be set
-# viscosity_wall = ViscosityAdami(nu=nu)
-# viscosity_wall = ViscosityMorris(nu=nu)
+# For a no-slip boundary condition, define a wall viscosity:
+# viscosity_wall = viscosity_fluid
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              boundary_density_calculator,
@@ -114,15 +103,9 @@ info_callback = InfoCallback(interval=100)
 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
+extra_callback2 = nothing
 
 use_reinit = false
 density_reinit_cb = use_reinit ?
@@ -131,7 +114,7 @@ density_reinit_cb = use_reinit ?
 stepsize_callback = StepsizeCallback(cfl=0.9)
 
 callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback,
-                        density_reinit_cb, saving_paper)
+                        extra_callback2, density_reinit_cb)
 
 time_integration_scheme = CarpenterKennedy2N54(williamson_condition=false)
 sol = solve(ode, time_integration_scheme,

examples/fluid/dam_break_2d_gpu.jl

@@ -36,7 +36,7 @@ trixi_include(@__MODULE__,
               joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
               neighborhood_search=neighborhood_search,
               fluid_particle_spacing=fluid_particle_spacing,
-              tspan=tspan,
+              tspan=tspan, smoothing_length=smoothing_length,
               density_diffusion=density_diffusion,
               boundary_layers=boundary_layers, spacing_ratio=spacing_ratio,
               boundary_model=boundary_model,

src/TrixiParticles.jl

@@ -6,7 +6,7 @@ using Adapt: Adapt
 using Base: @propagate_inbounds
 using CSV: CSV
 using Dates
-using DataFrames: DataFrame
+using DataFrames: DataFrames, DataFrame
 using DelimitedFiles: DelimitedFiles
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect, PresetTimeCallback
 using FastPow: @fastpow

src/callbacks/post_process.jl

@@ -238,7 +238,11 @@ function (pp::PostprocessCallback)(integrator; from_initialize=false)
             # comes AFTER the `PostprocessCallback` in the `CallbackSet`.
             dv_ode, du_ode = zero(vu_ode).x
         else
-            dv_ode, du_ode = get_du(integrator).x
+            # Depending on the time integration scheme, this might call the RHS function
+            @trixi_timeit timer() "update du" begin
+                # Don't create sub-timers here to avoid cluttering the timer output
+                @notimeit timer() dv_ode, du_ode = get_du(integrator).x
+            end
         end
         v_ode, u_ode = vu_ode.x
         semi = integrator.p

test/examples/examples_fluid.jl

@@ -589,7 +589,7 @@
             @trixi_test_nowarn trixi_include(@__MODULE__,
                                              joinpath(examples_dir(), "fluid",
                                                       "dam_break_2d.jl"),
-                                             tspan=(0, 0.1), sol=nothing) [
+                                             tspan=(0, 0.3), sol=nothing) [
                 r"┌ Info: The desired tank length in y-direction .*\n",
                 r"└ New tank length in y-direction.*\n"]
 
@@ -607,7 +607,7 @@
             @trixi_test_nowarn trixi_include(@__MODULE__,
                                              joinpath(examples_dir(), "fluid",
                                                       "dam_break_2d.jl"),
-                                             tspan=(0, 0.1), sol=nothing,
+                                             tspan=(0, 0.3), sol=nothing,
                                              cfl=0.25) [
                 r"┌ Info: The desired tank length in y-direction .*\n",
                 r"└ New tank length in y-direction.*\n"]

test/examples/gpu.jl

@@ -505,19 +505,19 @@ end
                                           cut_off_bnd=false)
 
                 @test isapprox(Array(result.computed_density),
-                               Float32[500.33255, 893.09766, 997.7032, 1001.14355, 1001.234,
-                                       1001.0098, 1000.4352, 999.7572, 999.1139, 989.6319])
+                               Float32[500.96387, 859.06744, 989.0479, 1001.3735,
+                                       1001.30927, 1001.0831, 1000.7325, 1000.3575,
+                                       999.8011, 975.98553])
 
                 @test isapprox(Array(result.density),
-                               Float32[1002.3152, 1002.19653, 1001.99915, 1001.7685,
-                                       1001.5382,
-                                       1001.3093, 1001.0836, 1000.8649, 1000.635,
-                                       1000.4053])
+                               Float32[1002.1891, 1002.0748, 1001.8913, 1001.67126,
+                                       1001.4529, 1001.2382, 1001.0298, 1000.81964,
+                                       1000.594, 1000.3878])
 
                 @test isapprox(Array(result.pressure),
-                               Float32[5450.902, 5171.2856, 4706.551, 4163.9185, 3621.5042,
-                                       3082.6948, 2551.5725, 2036.1208, 1494.8608,
-                                       954.14355])
+                               Float32[5154.1177, 4885.1, 4452.6533, 3934.9075, 3420.5737,
+                                       2915.0933, 2424.0908, 1929.7888, 1398.8309,
+                                       913.2089])
             end
 
             @testset verbose=true "Plane" begin
@@ -527,22 +527,19 @@ end
 
                 result = interpolate_plane_2d(interpolation_start, interpolation_end,
                                               resolution, semi_new, semi_new.systems[1],
-                                              sol;
-                                              cut_off_bnd=false)
+                                              sol; cut_off_bnd=false)
 
                 @test isapprox(Array(result.computed_density),
-                               Float32[250.18625, 500.34482, 499.77225, 254.3632, 499.58026,
-                                       999.1413, 998.6351, 503.0122])
+                               Float32[250.47523, 500.98248, 500.49924, 253.77109,
+                                       500.0491, 1000.0818, 999.6527, 503.08704])
 
                 @test isapprox(Array(result.density),
-                               Float32[1002.34467, 1002.3365, 1001.5021, 999.7109,
-                                       1000.84863,
-                                       1000.8373, 1000.3423, 1000.20734])
+                               Float32[1002.2373, 1002.22516, 1001.4339, 999.3567,
+                                       1000.8318, 1000.80237, 1000.3408, 999.92017])
 
                 @test isapprox(Array(result.pressure),
-                               Float32[5520.0513, 5501.1846, 3536.2256, -680.5194,
-                                       1997.7814,
-                                       1971.0717, 805.8584, 488.4068])
+                               Float32[5267.867, 5239.2466, 3376.045, -1514.4237,
+                                       1958.2637, 1888.739, 802.3146, -187.81871])
             end
         end
     end

test/validation/validation.jl

@@ -38,50 +38,66 @@
     end
 
     @trixi_testset "dam_break_2d" begin
+        # Use `SerialUpdate()` to obtain consistent results when using multiple
+        # threads and a shorter tspan to speed up CI tests.
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(validation_dir(), "dam_break_2d",
-                                                  "validation_dam_break_2d.jl")) [
+                                                  "validation_dam_break_2d.jl"),
+                                         update_strategy=SerialUpdate(),
+                                         tspan=(0.0, 4 / sqrt(9.81 / 0.6))) [
             r"┌ Info: The desired tank length in y-direction.*\n",
-            r"└ New tank length in y-direction is set to.*\n"
+            r"└ New tank length in y-direction is set to.*\n",
+            r"WARNING: Method definition max_x_coord.*\n",
+            r"WARNING: Method definition interpolated_pressure.*\n"
         ]
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
 
         if Sys.ARCH === :aarch64
             # MacOS ARM produces slightly different pressure values than x86.
             # Note that pressure values are in the order of 1e5.
-            @test isapprox(error_edac_P1, 0, atol=5e-6)
-            @test isapprox(error_edac_P2, 0, atol=4e-11)
-            @test isapprox(error_wcsph_P1, 0, atol=400.0)
-            @test isapprox(error_wcsph_P2, 0, atol=0.03)
+            @test isapprox(error_edac_P1, 0, atol=eps(1e5))
+            @test isapprox(error_edac_P2, 0, atol=eps(1e5))
+            @test isapprox(error_wcsph_P1, 0, atol=eps(1e5))
+            @test isapprox(error_wcsph_P2, 0, atol=eps(1e5))
         elseif VERSION < v"1.11"
-            # 1.10 produces slightly different pressure values than 1.11.
-            # This is most likely due to muladd and FMA instructions in the
-            # density diffusion update (inside the StaticArrays matrix-vector product).
-            # Note that pressure values are in the order of 1e5.
-            @test isapprox(error_edac_P1, 0, atol=eps())
-            @test isapprox(error_edac_P2, 0, atol=eps())
-            @test isapprox(error_wcsph_P1, 0, atol=8.0)
-            @test isapprox(error_wcsph_P2, 0, atol=5e-4)
+            # 1.10 produces slightly different pressure values than 1.11
+            @test isapprox(error_edac_P1, 0, atol=eps(1e5))
+            @test isapprox(error_edac_P2, 0, atol=eps(1e5))
+            @test isapprox(error_wcsph_P1, 0, atol=eps(1e5))
+            @test isapprox(error_wcsph_P2, 0, atol=eps(1e5))
         else
-            # Reference values are computed with 1.11
-            @test isapprox(error_edac_P1, 0, atol=5e-7)
-            @test isapprox(error_edac_P2, 0, atol=4e-11)
+            # Reference values are computed with 1.11 on x86
             @test isapprox(error_wcsph_P1, 0, atol=eps())
             @test isapprox(error_wcsph_P2, 0, atol=eps())
+            # Why are these errors not zero?
+            @test isapprox(error_edac_P1, 0, atol=eps(1e5))
+            @test isapprox(error_edac_P2, 0, atol=eps(1e5))
         end
 
         # Ignore method redefinitions from duplicate `include("../validation_util.jl")`
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(validation_dir(), "dam_break_2d",
-                                                  "plot_dam_break_results.jl")) [
+                                                  "plot_pressure_sensors.jl")) [
             r"WARNING: Method definition linear_interpolation.*\n",
             r"WARNING: Method definition interpolated_mse.*\n",
-            r"WARNING: Method definition extract_number_from_filename.*\n",
-            r"WARNING: Method definition extract_resolution_from_filename.*\n"
+            r"WARNING: Method definition extract_number_from_filename.*\n"
+        ]
+        # Verify number of plots
+        @test length(axs_edac[1].scene.plots) >= 2
+        @test length(axs_wcsph[1].scene.plots) >= 2
+
+        # Ignore method redefinitions from duplicate `include("../validation_util.jl")`
+        @trixi_test_nowarn trixi_include(@__MODULE__,
+                                         joinpath(validation_dir(), "dam_break_2d",
+                                                  "plot_surge_front.jl")) [
+            r"WARNING: Method definition linear_interpolation.*\n",
+            r"WARNING: Method definition interpolated_mse.*\n",
+            r"WARNING: Method definition extract_number_from_filename.*\n"
         ]
         # Verify number of plots
         @test length(axs_edac[1].scene.plots) >= 2
+        @test length(axs_wcsph[1].scene.plots) >= 2
     end
 
     @trixi_testset "TGV_2D" begin

validation/dam_break_2d/README.md

@@ -1,8 +1,33 @@
+The files in this folder provide a 2D dam break validation case for TrixiParticles.jl
+based on the following references:
+1. J. J. De Courcy, T. C. S. Rendall, L. Constantin, B. Titurus, J. E. Cooper.
+   "Incompressible δ-SPH via artificial compressibility".
+   In: Computer Methods in Applied Mechanics and Engineering, Volume 420 (2024),
+   https://doi.org/10.1016/j.cma.2023.116700
+2. S. Marrone, M. Antuono, A. Colagrossi, G. Colicchio, D. le Touzé, G. Graziani.
+   "δ-SPH model for simulating violent impact flows".
+   In: Computer Methods in Applied Mechanics and Engineering, Volume 200, Issues 13–16 (2011), pages 1526–1542.
+   https://doi.org/10.1016/J.CMA.2010.12.016
+3. J. C. Martin, W. J. Moyce, William George Penney, A. T. Price, C. K. Thornhill.
+   "Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane", Phil. Tran. R. Soc. London, Volume 244, Issue 882 (1952).
+   https://doi.org/10.1098/rsta.1952.0006
+
 The following files are provided here:
 
-1) Validation simulation: validation_dam_break_2d.jl
-2) Comparison with TrixiParticles.jl and literature reference: plot_dam_break_results.jl
-3) TrixiParticles.jl reference files: validation_reference_[0015, 00075, 0001875].j (0025 is the default CI resolution. Note that resolution 0015 takes around 2-5 minutes, 00075 around 5-10 minutes, 0001875 around 2-4 hours).
-4) exp_surge_front.csv was manually extracted from Martin and Moyce, "Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane", Phil. Tran. R. Soc. London, 1952, doi: https://doi.org/10.1098/rsta.1952.0006
-5) exp_pressure_sensor_P*.csv was manually extracted from  B. Buchner, "Green Water on Ship-type Offshore Structures", Ph.D. Thesis, Delft University of Technology, 2002.
-6) sim_pressure_sensor_P*.csv was manually extracted from  S. Marrone et al., "δ-SPH model for simulating violent impact flows", Computer Methods in Applied Mechanics and Engineering, 2011, doi: https://doi.org/10.1016/j.cma.2010.12.016
+1. `setup_marrone_2011.jl`: Setup file for the dam break simulation based on Marrone et al. (2011).
+   This can be used to reproduce the images of the breaking wave in Marrone et al. (2011).
+2. `validation_dam_break_2d.jl`: Script to run the dam break simulation
+   based on De Courcy et al. (2024) and extract pressure sensor data.
+   **Note**: This paper used particle shifting, which is not yet implemented
+   in TrixiParticles.jl for free surface simulations.
+3. `plot_pressure_sensors.jl`: Script to plot the pressure sensor data and compare with
+   De Courcy et al. (2024). This script can plot both the reference data included
+   in this folder and simulation results from the `out` folder (if available).
+4. `plot_surge_front.jl`: Script to plot the surge front position and compare with
+   Martin et al. (1952). This script can plot both the reference data included
+   in this folder and simulation results from the `out` folder (if available).
+
+Note that the reference files for resolutions 40 and 80 were obtained with `SerialUpdate()`
+and `--check-bounds=yes` on an x86 CPU.
+The reference file for resolution 400 was obtained with `update_strategy=nothing`
+and `--check-bounds=auto`, also on an x86 CPU.