Adaptive Particle Refinement: Update particle spacing

src/TrixiParticles.jl

@@ -45,6 +45,7 @@ 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.
@@ -72,7 +73,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,129 @@
+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_update!,
+                                save_positions=(false, false))
+    else
+        # The first one is the `condition`, the second the `affect!`
+        return DiscreteCallback(refinement_callback, refinement_callback,
+                                initialize=initial_update!,
+                                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
+
+    # 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!(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

@@ -134,8 +134,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 +407,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.

src/refinement/refinement.jl

@@ -0,0 +1,133 @@
+include("refinement_criteria.jl")
+
+struct ParticleRefinement{RC, ELTYPE}
+    refinement_criteria :: RC
+    max_spacing_ratio   :: ELTYPE
+    mass_ref            :: Vector{ELTYPE} # length(mass_ref) == nparticles
+end
+
+function ParticleRefinement(; max_spacing_ratio,
+                            refinement_criteria=SpatialRefinementCriterion())
+    mass_ref = Vector{eltype(max_spacing_ratio)}()
+
+    if refinement_criteria isa Tuple
+        refinement_criteria = (refinement_criteria,)
+    end
+    return ParticleRefinement(refinement_criteria, max_spacing_ratio, mass_ref)
+end
+
+resize_refinement!(system) = system
+
+function resize_refinement!(system::FluidSystem)
+    resize_refinement!(system, system.particle_refinement)
+end
+
+resize_refinement!(system, ::Nothing) = system
+
+function resize_refinement!(system, particle_refinement)
+    resize!(particle_refinement.mass_ref, nparticles(system))
+
+    return system
+end
+
+function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t)
+    check_refinement_criteria!(semi, v_ode, u_ode)
+
+    # Update the spacing of particles (Algorthm 1)
+    update_particle_spacing(semi, u_ode)
+
+    # Split the particles (Algorithm 2)
+
+    # Merge the particles (Algorithm 3)
+
+    # Shift the particles
+
+    # Correct the particles
+
+    # Update smoothing lengths
+
+    # Resize neighborhood search
+
+    return semi
+end
+
+function update_particle_spacing(semi, u_ode)
+    foreach_system(semi) do system
+        u = wrap_u(u_ode, system, semi)
+        update_particle_spacing(system, u, semi)
+    end
+end
+
+# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl`
+@inline update_particle_spacing(system, u, semi) = systemerror
+
+@inline function update_particle_spacing(system::FluidSystem, u, semi)
+    update_particle_spacing(system, system.particle_refinement, u, semi)
+end
+
+@inline update_particle_spacing(system, ::Nothing, u, semi) = system
+
+@inline function update_particle_spacing(system::FluidSystem, particle_refinement, u, semi)
+    (; smoothing_length, smoothing_length_factor) = system.cache
+    (; mass_ref) = particle_refinement
+
+    system_coords = current_coordinates(u, system)
+
+    for particle in eachparticle(system)
+        dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, system_coords, particle)
+
+        if dp_max / dp_min < max_spacing_ratio^3
+            new_spacing = min(dp_max, max_spacing_ratio * dp_min)
+        else
+            new_spacing = dp_avg
+        end
+
+        smoothing_length[particle] = smoothing_length_factor * new_spacing
+        mass_ref[particle] = system.density[particle] * new_spacing^(ndims(system))
+    end
+
+    return system
+end
+
+@inline function min_max_avg_spacing(system, semi, system_coords, particle)
+    dp_min = Inf
+    dp_max = zero(eltype(system))
+    dp_avg = zero(eltype(system))
+    counter_neighbors = 0
+
+    foreach_system(semi) do neighbor_system
+        neighborhood_search = get_neighborhood_search(particle_system, neighbor_system,
+                                                      semi)
+
+        u_neighbor_system = wrap_u(u_ode, neighbor, semi)
+        neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)
+
+        PointNeighbors.foreach_neighbor(system_coords, neighbor_coords, neighborhood_search,
+                                        particle) do particle, neighbor, pos_diff, distance
+            dp_neighbor = particle_spacing(neighbor_system, neighbor)
+
+            dp_min = min(dp_min, dp_neighbor)
+            dp_max = max(dp_max, dp_neighbor)
+            dp_avg += dp_neighbor
+
+            counter_neighbors += 1
+        end
+    end
+
+    dp_avg / counter_neighbors
+
+    return dp_min, dp_max, dp_avg
+end
+
+@inline particle_spacing(system, particle) = system.initial_condition.particle_spacing
+
+@inline function particle_spacing(system::FluidSystem, particle)
+    return particle_spacing(system, system.particle_refinement, particle)
+end
+
+@inline particle_spacing(system, ::Nothing, _) = system.initial_condition.particle_spacing
+
+@inline function particle_spacing(system, refinement, particle)
+    (; smoothing_length_factor) = system.cache
+    return smoothing_length(system, particle) / smoothing_length_factor
+end