Mixed kernel gradient (#261)

* add gradient correction stuff * impl * up * finish impl * update * fix testcases * fix * fix tests * format * working * format * format * fix * working trixi_include * fix * cleanup * format * fix * fix tests * fix * ups * typo * up * fix * format * update fix errors * fix * update * up * fix tests * fix tests * np * remove comment * revert some changes * up * up * update * format * cleanup * revert * revert * implement correction on wall * fix tests * fix pressure_acceleration for summation * format * revert * fix * fixes * fix * fix pressure_acceleration dispatches * typo * fix * rename * some review fixes * review * more fixes * add solid system type * review * fix tlsph call * fix comment * review * review * format * review * move * fix doc * review fix * more review * review fix * format * change parameter list * reorder params * format * fix * fix * fix * remove dispatched functions * format * fix comments * fix * remove system * format * change parameter list * review * unused variable * break * edit todo * remove breaks * format * rename grad_kernel to W_a
trixi-framework · Jan 5, 2024 · f789ffb · f789ffb
1 parent 88d287d
commit f789ffb
Show file tree

Hide file tree

Showing 19 changed files with 653 additions and 205 deletions.
diff --git a/examples/fluid/dam_break_2d.jl b/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);
diff --git a/examples/fluid/dam_break_2d_corrections.jl b/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
diff --git a/src/TrixiParticles.jl b/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