Add particle refinement

src/TrixiParticles.jl

@@ -20,7 +20,8 @@ using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @series
 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
@@ -50,14 +51,15 @@ include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("visualization/write2vtk.jl")
 include("visualization/recipes_plots.jl")
+include("particle_refinement/particle_refinement.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
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,
@@ -84,5 +86,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure,
 export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d,
        interpolate_plane_2d_vtk
 export SurfaceTensionAkinci, CohesionForceAkinci
+export ParticleRefinement, RefinementZone
+export CubicSplitting, TriangularSplitting, HexagonalSplitting
 
 end # module

src/callbacks/callbacks.jl

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

src/callbacks/refinement_callback.jl

@@ -0,0 +1,145 @@
+mutable struct ParticleRefinementCallback{I}
+    interval           :: I
+    ranges_u_cache     :: Tuple
+    ranges_v_cache     :: Tuple
+    nparticles_cache   :: Tuple
+    eachparticle_cache :: Tuple
+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)
+    cb(integrator)
+end
+
+# condition
+function (refinement_callback::ParticleRefinementCallback)(u, t, integrator)
+    (; interval) = refinement_callback
+
+    # With error-based step size control, some steps can be rejected. Thus,
+    #   `integrator.iter >= integrator.stats.naccept`
+    #    (total #steps)       (#accepted steps)
+    # We need to check the number of accepted steps since callbacks are not
+    # activated after a rejected step.
+    return integrator.stats.naccept % interval == 0
+end
+
+# affect
+function (refinement_callback::ParticleRefinementCallback)(integrator)
+    t = integrator.t
+    semi = integrator.p
+    v_ode, u_ode = integrator.u.x
+
+    # Update NHS
+    @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi)
+
+    # 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!(v_ode, u_ode, v_tmp, u_tmp, semi, refinement_callback, t)
+
+    # Resize neighborhood search
+    foreach_system(semi) do system
+        foreach_system(semi) do neighbor_system
+            search = get_neighborhood_search(system, neighbor_system, semi)
+            u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+
+            resize_nhs!(search, system, neighbor_system, u_neighbor)
+        end
+    end
+
+    @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
+                                                                          semi, 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/semidiscretization.jl

@@ -73,14 +73,9 @@ 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))
+    systems = create_child_systems(systems)
+
+    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 +87,21 @@ function Semidiscretization(systems...;
     return Semidiscretization(systems, ranges_u, ranges_v, searches)
 end
 
+function ranges_vu(systems)
+    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))
+
+    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))
+
+    # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct
+    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
@@ -370,7 +380,9 @@ end
 @inline function wrap_v(v_ode, system, semi)
     (; ranges_v) = semi
 
-    range = ranges_v[system_indices(system, semi)]
+    # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct.
+    # The `Vector` has length 1, thus, take the first element.
+    range = first(ranges_v[system_indices(system, semi)])
 
     @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system)
 
@@ -381,7 +393,9 @@ end
 @inline function wrap_u(u_ode, system, semi)
     (; ranges_u) = semi
 
-    range = ranges_u[system_indices(system, semi)]
+    # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct.
+    # The `Vector` has length 1, thus, take the first element.
+    range = first(ranges_u[system_indices(system, semi)])
 
     @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system)
 

src/particle_refinement/particle_refinement.jl

@@ -0,0 +1 @@
+include("refinement.jl")