Adaptive Particle Refinement: Particle refinement prototype

examples/fluid/particle_refinement_2d.jl

@@ -0,0 +1,86 @@
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+fluid_particle_spacing = 0.05
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 3
+
+# ==========================================================================================
+# ==== Experiment Setup
+gravity = 9.81
+tspan = (0.0, 1.0)
+
+# Boundary geometry and initial fluid particle positions
+initial_fluid_size = (2.0, 2.0)
+tank_size = (2.0, 2.0)
+
+fluid_density = 1000.0
+sound_speed = 10.0
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7, clip_negative_pressure=false)
+
+tank = RectangularTank(fluid_particle_spacing, initial_fluid_size, tank_size, fluid_density,
+                       n_layers=boundary_layers, spacing_ratio=2.0)
+
+sphere = SphereShape(0.5 * fluid_particle_spacing, 0.3, (1.0, 1.0),
+                     fluid_density, sphere_type=RoundSphere())
+# ==========================================================================================
+# ==== Fluid
+
+fluid = setdiff(tank.fluid, sphere)
+
+smoothing_length = 1.2 * fluid_particle_spacing
+smoothing_kernel = SchoenbergQuinticSplineKernel{2}()
+
+fluid_density_calculator = ContinuityDensity()
+pressure = 1000.0
+particle_refinement = ParticleRefinement(;
+                                         refinement_pattern=TrixiParticles.HexagonalSplitting(),
+                                         max_spacing_ratio=1.2)
+fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length,
+                                           sound_speed; viscosity=ViscosityAdami(; nu=1e-4),
+                                           particle_refinement=particle_refinement,
+                                           transport_velocity=TransportVelocityAdami(pressure),
+                                           acceleration=(0.0, -gravity))
+
+# ==========================================================================================
+# ==== Boundary
+
+# This is to set another boundary density calculation with `trixi_include`
+boundary_density_calculator = AdamiPressureExtrapolation()
+
+# This is to set wall viscosity with `trixi_include`
+viscosity_wall = nothing
+boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass,
+                                             boundary_density_calculator,
+                                             smoothing_kernel, smoothing_length,
+                                             viscosity=viscosity_wall)
+boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, movement=nothing)
+
+boundary_model_sphere = BoundaryModelDummyParticles(sphere.density, sphere.mass,
+                                                    boundary_density_calculator,
+                                                    smoothing_kernel, smoothing_length,
+                                                    viscosity=viscosity_wall)
+boundary_system_sphere = BoundarySPHSystem(sphere, boundary_model_sphere, movement=nothing)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(fluid_system, boundary_system, boundary_system_sphere)
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=50)
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+# This is to easily add a new callback with `trixi_include`
+extra_callback = nothing
+
+refinement_callback = ParticleRefinementCallback(interval=1)
+
+callbacks = CallbackSet(info_callback, saving_callback, extra_callback, UpdateCallback(),
+                        refinement_callback)
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks);

src/TrixiParticles.jl

@@ -22,7 +22,8 @@ using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @series
 using Random: seed!
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
-                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem
+                 get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem,
+                 RecursiveArrayTools
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
@@ -45,21 +46,23 @@ include("callbacks/callbacks.jl")
 include("general/general.jl")
 include("setups/setups.jl")
 include("schemes/schemes.jl")
+include("refinement/refinement.jl")
 
 # Note that `semidiscretization.jl` depends on the system types and has to be
 # included separately. `gpu.jl` in turn depends on the semidiscretization type.
 include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("visualization/write2vtk.jl")
 include("visualization/recipes_plots.jl")
+include("refinement/resize.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
        BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow,
        OutFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
-       PostprocessCallback, StepsizeCallback, UpdateCallback
+       PostprocessCallback, StepsizeCallback, UpdateCallback, ParticleRefinementCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
@@ -72,7 +75,7 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr
 export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation,
        PressureMirroring, PressureZeroing, BoundaryModelLastiwka,
        BernoulliPressureExtrapolation
-
+export ParticleRefinement, SpatialRefinementCriterion
 export BoundaryMovement
 export examples_dir, validation_dir, trixi_include
 export trixi2vtk

src/callbacks/callbacks.jl

@@ -31,3 +31,4 @@ include("density_reinit.jl")
 include("post_process.jl")
 include("stepsize.jl")
 include("update.jl")
+include("refinement.jl")

src/callbacks/refinement.jl

@@ -0,0 +1,126 @@
+struct ParticleRefinementCallback{I}
+    interval::I
+end
+
+function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0)
+    if dt > 0 && interval !== -1
+        throw(ArgumentError("setting both `interval` and `dt` is not supported"))
+    end
+
+    # Update in intervals in terms of simulation time
+    if dt > 0
+        interval = Float64(dt)
+
+        # Update every time step (default)
+    elseif interval == -1
+        interval = 1
+    end
+
+    refinement_callback = ParticleRefinementCallback(interval)
+
+    if dt > 0
+        # Add a `tstop` every `dt`, and save the final solution.
+        return PeriodicCallback(refinement_callback, dt,
+                                initialize=initial_refinement!,
+                                save_positions=(false, false))
+    else
+        # The first one is the `condition`, the second the `affect!`
+        return DiscreteCallback(refinement_callback, refinement_callback,
+                                initialize=initial_refinement!,
+                                save_positions=(false, false))
+    end
+end
+
+# initialize
+function initial_refinement!(cb, u, t, integrator)
+    # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`)
+    # or `cb.affect!.affect!` (with `PeriodicCallback`).
+    # Let recursive dispatch handle this.
+
+    initial_refinement!(cb.affect!, u, t, integrator)
+end
+
+function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator)
+    semi = integrator.p
+
+    foreach_system(semi) do system
+        resize_refinement!(system)
+    end
+
+    cb(integrator)
+end
+
+# condition
+function (refinement_callback::ParticleRefinementCallback)(u, t, integrator)
+    (; interval) = refinement_callback
+
+    return condition_integrator_interval(integrator, interval)
+end
+
+# affect
+function (refinement_callback::ParticleRefinementCallback)(integrator)
+    t = integrator.t
+    semi = integrator.p
+    v_ode, u_ode = integrator.u.x
+
+    # Basically `get_tmp_cache(integrator)` to write into in order to be non-allocating
+    # https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Caches
+    v_tmp, u_tmp = integrator.cache.tmp.x
+
+    v_tmp .= v_ode
+    u_tmp .= u_ode
+
+    refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t)
+
+    resize!(integrator, (length(v_ode), length(u_ode)))
+
+    # Tell OrdinaryDiffEq that u has been modified
+    u_modified!(integrator, true)
+
+    return integrator
+end
+
+Base.resize!(a::RecursiveArrayTools.ArrayPartition, sizes::Tuple) = resize!.(a.x, sizes)
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback})
+    @nospecialize cb # reduce precompilation time
+    print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")")
+end
+
+function Base.show(io::IO,
+                   cb::DiscreteCallback{<:Any,
+                                        <:PeriodicCallbackAffect{<:ParticleRefinementCallback}})
+    @nospecialize cb # reduce precompilation time
+    print(io, "ParticleRefinementCallback(dt=", cb.affect!.affect!.interval, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        refinement_cb = cb.affect!
+        setup = [
+            "interval" => refinement_cb.interval
+        ]
+        summary_box(io, "ParticleRefinementCallback", setup)
+    end
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any,
+                                        <:PeriodicCallbackAffect{<:ParticleRefinementCallback}})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        refinement_cb = cb.affect!.affect!
+        setup = [
+            "dt" => refinement_cb.interval
+        ]
+        summary_box(io, "ParticleRefinementCallback", setup)
+    end
+end

src/general/corrections.jl

@@ -220,7 +220,7 @@ function compute_correction_values!(system,
             kernel_correction_coefficient[particle] += volume *
                                                        smoothing_kernel(system, distance)
             if distance > sqrt(eps())
-                tmp = volume * smoothing_kernel_grad(system, pos_diff, distance)
+                tmp = volume * smoothing_kernel_grad(system, pos_diff, distance, particle)
                 for i in axes(dw_gamma, 1)
                     dw_gamma[i, particle] += tmp[i]
                 end
@@ -311,7 +311,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search,
                                                            pos_diff, distance
         volume = mass[neighbor] / density_fun(neighbor)
 
-        grad_kernel = smoothing_kernel_grad(system, pos_diff, distance)
+        grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
 
         iszero(grad_kernel) && return
 
@@ -349,7 +349,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system,
 
             function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance,
                                          smoothing_length, system, particle)
-                return smoothing_kernel_grad(system, pos_diff, distance)
+                return smoothing_kernel_grad(system, pos_diff, distance, particle)
             end
 
             # Compute gradient of corrected kernel

src/general/density_calculators.jl

@@ -66,7 +66,7 @@ function summation_density!(system, semi, u, u_ode, density;
                                points=particles) do particle, neighbor,
                                                     pos_diff, distance
             mass = hydrodynamic_mass(neighbor_system, neighbor)
-            density[particle] += mass * smoothing_kernel(system, distance)
+            density[particle] += mass * smoothing_kernel(system, distance, particle)
         end
     end
 end

src/general/semidiscretization.jl

@@ -73,14 +73,7 @@ function Semidiscretization(systems...;
     # Other checks might be added here later.
     check_configuration(systems)
 
-    sizes_u = [u_nvariables(system) * n_moving_particles(system)
-               for system in systems]
-    ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i])
-                     for i in eachindex(sizes_u))
-    sizes_v = [v_nvariables(system) * n_moving_particles(system)
-               for system in systems]
-    ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i])
-                     for i in eachindex(sizes_v))
+    ranges_v, ranges_u = ranges_vu(systems)
 
     # Create a tuple of n neighborhood searches for each of the n systems.
     # We will need one neighborhood search for each pair of systems.
@@ -92,6 +85,16 @@ function Semidiscretization(systems...;
     return Semidiscretization(systems, ranges_u, ranges_v, searches)
 end
 
+function ranges_vu(systems)
+    sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems]
+    ranges_u = [(sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) for i in eachindex(sizes_u)]
+
+    sizes_v = [v_nvariables(system) * n_moving_particles(system) for system in systems]
+    ranges_v = [(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)]
+
+    return ranges_v, ranges_u
+end
+
 # Inline show function e.g. Semidiscretization(neighborhood_search=...)
 function Base.show(io::IO, semi::Semidiscretization)
     @nospecialize semi # reduce precompilation time
@@ -134,8 +137,8 @@ function create_neighborhood_search(neighborhood_search, system, neighbor)
 end
 
 @inline function compact_support(system, neighbor)
-    (; smoothing_kernel, smoothing_length) = system
-    return compact_support(smoothing_kernel, smoothing_length)
+    (; smoothing_kernel) = system
+    return compact_support(smoothing_kernel, smoothing_length(system, 1)) # TODO
 end
 
 @inline function compact_support(system::OpenBoundarySPHSystem, neighbor)
@@ -407,7 +410,9 @@ end
 function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization)
     (; systems) = semi
 
-    return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems)
+    return 1.0
+    # TODO
+    # return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems)
 end
 
 function drift!(du_ode, v_ode, u_ode, semi, t)

src/general/system.jl

@@ -106,30 +106,23 @@ end
 
 @inline set_particle_pressure!(v, system, particle, pressure) = v
 
-@inline function smoothing_kernel(system, distance)
-    (; smoothing_kernel, smoothing_length) = system
-    return kernel(smoothing_kernel, distance, smoothing_length)
+@inline function smoothing_kernel(system, distance, particle)
+    (; smoothing_kernel) = system
+    return kernel(smoothing_kernel, distance, smoothing_length(system, particle))
 end
 
-@inline function smoothing_kernel_deriv(system, distance)
-    (; smoothing_kernel, smoothing_length) = system
-    return kernel_deriv(smoothing_kernel, distance, smoothing_length)
-end
-
-@inline function smoothing_kernel_grad(system, pos_diff, distance)
-    return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length)
-end
-
-@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance)
+@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle)
     (; smoothing_kernel, smoothing_length) = system.boundary_model
 
     return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length)
 end
 
 @inline function smoothing_kernel_grad(system, pos_diff, distance, particle)
-    return corrected_kernel_grad(system.smoothing_kernel, pos_diff, distance,
-                                 system.smoothing_length, system.correction, system,
-                                 particle)
+    (; smoothing_kernel) = system
+
+    return corrected_kernel_grad(smoothing_kernel, pos_diff, distance,
+                                 smoothing_length(system, particle),
+                                 system.correction, system, particle)
 end
 
 # System update orders. This can be dispatched if needed.