Merge branch 'main' into oriented-bbox

Author: LasNikas

.buildkite/pipeline.yml

@@ -2,7 +2,7 @@ steps:
   - label: "CUDA"
     plugins:
       - JuliaCI/julia#v1:
-          version: "1"
+          version: "1.11"
     agents:
       queue: "juliagpu"
       cuda: "*"
@@ -16,7 +16,7 @@ steps:
   - label: "AMDGPU"
     plugins:
       - JuliaCI/julia#v1:
-          version: "1"
+          version: "1.11"
     agents:
       queue: "juliagpu"
       rocm: "*"
@@ -30,7 +30,7 @@ steps:
   - label: "Metal"
     plugins:
       - JuliaCI/julia#v1:
-          version: "1"
+          version: "1.11"
     agents:
       queue: "juliaecosystem"
       os: "macos"

.github/workflows/ci.yml

@@ -43,7 +43,7 @@ jobs:
       matrix:
         version:
           - '1.10'
-          - '1'
+          - '1.11'
         os:
           - ubuntu-latest
         include:
@@ -53,9 +53,9 @@ jobs:
             os: windows-latest
           - version: '1.10'
             os: macos-14
-          - version: '1'
+          - version: '1.11'
             os: windows-latest
-          - version: '1'
+          - version: '1.11'
             os: macos-14
 
     steps:

.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

@@ -404,7 +404,7 @@
     end
 
     @trixi_testset "fluid/pipe_flow_3d.jl" begin
-        @trixi_test_nowarn trixi_include(@__MODULE__,
+        @trixi_test_nowarn trixi_include(@__MODULE__, tspan=(0, 0.1),
                                          joinpath(examples_dir(), "fluid",
                                                   "pipe_flow_3d.jl"))
         @test sol.retcode == ReturnCode.Success
@@ -415,7 +415,7 @@
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(examples_dir(), "fluid",
                                                   "poiseuille_flow_2d.jl"),
-                                         tspan=(0.0, 0.1))
+                                         tspan=(0.0, 0.02))
         @test fluid_system isa WeaklyCompressibleSPHSystem
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
@@ -425,7 +425,7 @@
         @trixi_test_nowarn trixi_include(@__MODULE__,
                                          joinpath(examples_dir(), "fluid",
                                                   "poiseuille_flow_2d.jl"),
-                                         tspan=(0.0, 0.1), wcsph=false)
+                                         tspan=(0.0, 0.02), wcsph=false)
         @test fluid_system isa EntropicallyDampedSPHSystem
         @test sol.retcode == ReturnCode.Success
         @test count_rhs_allocations(sol, semi) == 0
@@ -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