Merge branch 'main' into fix-visualization

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/[email protected].23
+        uses: crate-ci/[email protected].26

.gitignore

@@ -10,3 +10,10 @@ docs/src/license.md
 docs/src/code_of_conduct.md
 docs/src/news.md
 run
+*.json
+*.png
+*.gif
+*.svg
+*.vtu
+*.pvd
+*.mp4

examples/fluid/dam_break_2d.jl

@@ -13,7 +13,7 @@ using OrdinaryDiffEq
 fluid_particle_spacing = 0.02
 
 # Change spacing ratio to 3 and boundary layers to 1 when using Monaghan-Kajtar boundary model
-boundary_layers = 4
+boundary_layers = 3
 spacing_ratio = 1
 
 boundary_particle_spacing = fluid_particle_spacing / spacing_ratio
@@ -31,7 +31,8 @@ fluid_density = 1000.0
 atmospheric_pressure = 100000.0
 sound_speed = 20 * sqrt(gravity * initial_fluid_size[2])
 state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure,
-                                   background_pressure=atmospheric_pressure)
+                                   background_pressure=atmospheric_pressure,
+                                   clip_negative_pressure=false)
 
 tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
                        n_layers=boundary_layers, spacing_ratio=spacing_ratio,
@@ -45,6 +46,7 @@ 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)
 
 fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator,
                                            state_equation, smoothing_kernel,
@@ -58,7 +60,8 @@ boundary_density_calculator = AdamiPressureExtrapolation()
 boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
                                              state_equation=state_equation,
                                              boundary_density_calculator,
-                                             smoothing_kernel, smoothing_length)
+                                             smoothing_kernel, smoothing_length,
+                                             correction=nothing)
 
 boundary_system = BoundarySPHSystem(tank.boundary, boundary_model)
 
@@ -68,7 +71,7 @@ semi = Semidiscretization(fluid_system, boundary_system,
                           neighborhood_search=GridNeighborhoodSearch)
 ode = semidiscretize(semi, tspan)
 
-info_callback = InfoCallback(interval=100)
+info_callback = InfoCallback(interval=250)
 saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
 
 use_reinit = false
@@ -87,5 +90,5 @@ callbacks = CallbackSet(info_callback, saving_callback, density_reinit_cb)
 sol = solve(ode, RDPK3SpFSAL49(),
             abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
             reltol=1e-5, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
-            dtmax=1e-2, # Limit stepsize to prevent crashing
+            dtmax=1e-3, # Limit stepsize to prevent crashing
             save_everystep=false, callback=callbacks);

examples/fluid/dam_break_2d_corrections.jl

@@ -3,19 +3,35 @@ using TrixiParticles
 fluid_density = 1000.0
 
 particle_spacing = 0.05
-smoothing_length = 1.15 * particle_spacing
 
-boundary_density_calculator = SummationDensity()
-
-gravity = 9.81
-tspan = (0.0, 5.7 / sqrt(gravity))
+tspan = (0.0, 5.7 / sqrt(9.81))
 
 correction_dict = Dict(
-    "no_correction" => Nothing(),
+    "no_correction" => nothing,
     "shepard_kernel_correction" => ShepardKernelCorrection(),
     "akinci_free_surf_correction" => AkinciFreeSurfaceCorrection(fluid_density),
-    "kernel_gradient_summation_correction" => KernelGradientCorrection(),
-    "kernel_gradient_continuity_correction" => KernelGradientCorrection(),
+    "kernel_gradient_summation_correction" => KernelCorrection(),
+    "kernel_gradient_continuity_correction" => KernelCorrection(),
+    "blended_gradient_summation_correction" => BlendedGradientCorrection(0.5),
+    "blended_gradient_continuity_correction" => BlendedGradientCorrection(0.2),
+    "gradient_summation_correction" => GradientCorrection(),
+    "mixed_kernel_gradient_summation_correction" => MixedKernelGradientCorrection(),
+    # "gradient_continuity_correction" => GradientCorrection(),
+    # "mixed_kernel_gradient_continuity_correction" => MixedKernelGradientCorrection(),
+)
+
+smoothing_length_dict = Dict(
+    "no_correction" => 3.0 * particle_spacing,
+    "shepard_kernel_correction" => 3.0 * particle_spacing,
+    "akinci_free_surf_correction" => 3.0 * particle_spacing,
+    "kernel_gradient_summation_correction" => 4.0 * particle_spacing,
+    "kernel_gradient_continuity_correction" => 3.5 * particle_spacing,
+    "blended_gradient_summation_correction" => 3.0 * particle_spacing,
+    "blended_gradient_continuity_correction" => 4.0 * particle_spacing,
+    "gradient_summation_correction" => 3.5 * particle_spacing,
+    "mixed_kernel_gradient_summation_correction" => 3.5 * particle_spacing,
+    "gradient_continuity_correction" => 4.5 * particle_spacing,
+    "mixed_kernel_gradient_continuity_correction" => 4.0 * particle_spacing,
 )
 
 density_calculator_dict = Dict(
@@ -24,33 +40,55 @@ density_calculator_dict = Dict(
     "akinci_free_surf_correction" => SummationDensity(),
     "kernel_gradient_summation_correction" => SummationDensity(),
     "kernel_gradient_continuity_correction" => ContinuityDensity(),
+    "blended_gradient_summation_correction" => SummationDensity(),
+    "blended_gradient_continuity_correction" => ContinuityDensity(),
+    "gradient_summation_correction" => SummationDensity(),
+    "gradient_continuity_correction" => ContinuityDensity(),
+    "mixed_kernel_gradient_summation_correction" => SummationDensity(),
+    "mixed_kernel_gradient_continuity_correction" => ContinuityDensity(),
+)
+
+smoothing_kernel_dict = Dict(
+    "no_correction" => WendlandC2Kernel{2}(),
+    "shepard_kernel_correction" => WendlandC2Kernel{2}(),
+    "akinci_free_surf_correction" => WendlandC2Kernel{2}(),
+    "kernel_gradient_summation_correction" => WendlandC6Kernel{2}(),
+    "kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(),
+    "blended_gradient_summation_correction" => WendlandC2Kernel{2}(),
+    "blended_gradient_continuity_correction" => WendlandC6Kernel{2}(),
+    "gradient_summation_correction" => WendlandC6Kernel{2}(),
+    "gradient_continuity_correction" => WendlandC6Kernel{2}(),
+    "mixed_kernel_gradient_summation_correction" => WendlandC6Kernel{2}(),
+    "mixed_kernel_gradient_continuity_correction" => WendlandC6Kernel{2}(),
 )
 
 trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
-              fluid_particle_spacing=particle_spacing, smoothing_length=smoothing_length,
+              fluid_particle_spacing=particle_spacing,
+              smoothing_length=3.0 * particle_spacing,
               boundary_density_calculator=ContinuityDensity(),
               fluid_density_calculator=ContinuityDensity(),
-              correction=Nothing(), use_reinit=true,
-              prefix="continuity_reinit", tspan=tspan)
-
-# Clip negative pressure to be able to use `SummationDensity`
-state_equation = StateEquationCole(sound_speed, 7, fluid_density, atmospheric_pressure,
-                                   background_pressure=atmospheric_pressure,
-                                   clip_negative_pressure=true)
+              correction=nothing, use_reinit=true,
+              prefix="continuity_reinit", tspan=tspan,
+              fluid_density=fluid_density, density_diffusion=nothing)
 
 for correction_name in keys(correction_dict)
     local fluid_density_calculator = density_calculator_dict[correction_name]
     local correction = correction_dict[correction_name]
+    local smoothing_kernel = smoothing_kernel_dict[correction_name]
+    local smoothing_length = smoothing_length_dict[correction_name]
 
     println("="^100)
     println("fluid/dam_break_2d.jl with ", correction_name)
 
     trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
                   fluid_particle_spacing=particle_spacing,
                   smoothing_length=smoothing_length,
-                  boundary_density_calculator=boundary_density_calculator,
+                  boundary_density_calculator=SummationDensity(),
                   fluid_density_calculator=fluid_density_calculator,
                   correction=correction, use_reinit=false,
-                  state_equation=state_equation,
-                  prefix="$(correction_name)", tspan=tspan)
+                  clip_negative_pressure=(fluid_density_calculator isa SummationDensity),
+                  smoothing_kernel=smoothing_kernel,
+                  prefix="$(correction_name)", tspan=tspan,
+                  fluid_density=fluid_density, density_diffusion=Nothing(),
+                  boundary_layers=5)
 end

examples/n_body/n_body_benchmark_trixi.jl

@@ -101,7 +101,14 @@ function symplectic_euler!(velocities, coordinates, semi)
 end
 
 # One RHS evaluation is so fast that timers make it multiple times slower.
-TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+# Disabling timers throws a warning, which we suppress here in order to make the tests pass.
+# For some reason, this only works with a file and not with devnull. See issue #332.
+filename = tempname()
+open(filename, "w") do f
+    redirect_stderr(f) do
+        TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+    end
+end
 
 @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi))
 
@@ -112,5 +119,9 @@ end
 
 @printf("%.9f\n", energy(velocities, coordinates, particle_system, semi))
 
-# Enable timers again.
-TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+# Enable timers again
+open(filename, "w") do f
+    redirect_stderr(f) do
+        TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+    end
+end

examples/n_body/n_body_solar_system.jl

@@ -44,7 +44,14 @@ saving_callback = SolutionSavingCallback(dt=10day)
 callbacks = CallbackSet(info_callback, saving_callback)
 
 # One RHS evaluation is so fast that timers make it multiple times slower.
-TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+# Disabling timers throws a warning, which we suppress here in order to make the tests pass.
+# For some reason, this only works with a file and not with devnull. See issue #332.
+filename = tempname()
+open(filename, "w") do f
+    redirect_stderr(f) do
+        TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+    end
+end
 
 sol = solve(ode, SymplecticEuler(),
             dt=1.0e5,
@@ -53,5 +60,9 @@ sol = solve(ode, SymplecticEuler(),
 @printf("%.9e\n", energy(ode.u0.x..., particle_system, semi))
 @printf("%.9e\n", energy(sol.u[end].x..., particle_system, semi))
 
-# Enable timers again.
-TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+# Enable timers again
+open(filename, "w") do f
+    redirect_stderr(f) do
+        TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+    end
+end

src/TrixiParticles.jl

@@ -5,7 +5,7 @@ using Reexport: @reexport
 using Dates
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 using FastPow: @fastpow
-using LinearAlgebra: norm, dot, I, tr
+using LinearAlgebra: norm, dot, I, tr, inv, pinv, det
 using Morton: cartesian2morton
 using MuladdMacro: @muladd
 using Polyester: Polyester, @batch
@@ -55,7 +55,8 @@ export examples_dir, trixi_include
 export trixi2vtk
 export RectangularTank, RectangularShape, SphereShape
 export VoxelSphere, RoundSphere, reset_wall!
-export ShepardKernelCorrection, KernelGradientCorrection, AkinciFreeSurfaceCorrection
+export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection,
+       GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection
 export nparticles
 
 end # module