Implement the Particle Shifting Technique (PST)

src/TrixiParticles.jl

@@ -41,10 +41,11 @@ using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 # `util.jl` needs to be first because of the macros `@trixi_timeit` and `@threaded`
 include("util.jl")
 include("general/system.jl")
-include("callbacks/callbacks.jl")
 include("general/general.jl")
 include("setups/setups.jl")
 include("schemes/schemes.jl")
+# `callbacks.jl` requires the system types to be defined
+include("callbacks/callbacks.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.
@@ -60,7 +61,8 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
        BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem
 export BoundaryZone, InFlow, OutFlow, BidirectionalFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
-       PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback
+       PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback,
+       ParticleShiftingCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller, TransportVelocityAdami
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,

src/callbacks/callbacks.jl

@@ -32,3 +32,4 @@ include("post_process.jl")
 include("stepsize.jl")
 include("update.jl")
 include("steady_state_reached.jl")
+include("particle_shifting.jl")

src/callbacks/particle_shifting.jl

@@ -0,0 +1,123 @@
+@doc raw"""
+    ParticleShiftingCallback()
+
+Callback to apply the Particle Shifting Technique by [Sun et al. (2017)](@cite Sun2017).
+Following the original paper, the callback is applied in every time step and not
+in every stage of a multi-stage time integration method to reduce the computational
+cost and improve the stability of the scheme.
+
+## References
+[Sun2017](@cite)
+"""
+function ParticleShiftingCallback()
+    # The first one is the `condition`, the second the `affect!`
+    return DiscreteCallback((particle_shifting_condition), particle_shifting!,
+                            save_positions=(false, false))
+end
+
+# `condition`
+function particle_shifting_condition(u, t, integrator)
+    return true
+end
+
+# `affect!`
+function particle_shifting!(integrator)
+    t = integrator.t
+    semi = integrator.p
+    v_ode, u_ode = integrator.u.x
+    dt = integrator.dt
+    # Internal cache vector, which is safe to use as temporary array
+    u_cache = first(get_tmp_cache(integrator))
+
+    # Update quantities that are stored in the systems. These quantities (e.g. pressure)
+    # still have the values from the last stage of the previous step if not updated here.
+    update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
+
+    @trixi_timeit timer() "particle shifting" foreach_system(semi) do system
+        u = wrap_u(u_ode, system, semi)
+        v = wrap_v(v_ode, system, semi)
+        particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt)
+    end
+
+    # Tell OrdinaryDiffEq that `u` has been modified
+    u_modified!(integrator, true)
+
+    return integrator
+end
+
+function particle_shifting!(u, v, system, v_ode, u_ode, semi, u_cache, dt)
+    return u
+end
+
+function particle_shifting!(u, v, system::WeaklyCompressibleSPHSystem, v_ode, u_ode, semi,
+                            u_cache, dt)
+    # Wrap the cache vector to an NDIMS x NPARTICLES matrix.
+    # We need this buffer because we cannot safely update `u` while iterating over it.
+    delta_r = wrap_u(u_cache, system, semi)
+    set_zero!(delta_r)
+
+    v_max = maximum(particle -> norm(current_velocity(v, system, particle)),
+                    eachparticle(system))
+
+    # TODO this needs to be adapted to multi-resolution.
+    # Section 3.2 explains what else needs to be changed.
+    Wdx = smoothing_kernel(system, particle_spacing(system, 1), 1)
+    h = smoothing_length(system, 1)
+
+    foreach_system(semi) do neighbor_system
+        u_neighbor = wrap_u(u_ode, neighbor_system, semi)
+        v_neighbor = wrap_v(v_ode, neighbor_system, semi)
+
+        system_coords = current_coordinates(u, system)
+        neighbor_coords = current_coordinates(u_neighbor, neighbor_system)
+
+        foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords,
+                               semi) do particle, neighbor, pos_diff, distance
+            m_b = hydrodynamic_mass(neighbor_system, neighbor)
+            rho_a = particle_density(v, system, particle)
+            rho_b = particle_density(v_neighbor, neighbor_system, neighbor)
+
+            kernel = smoothing_kernel(system, distance, particle)
+            grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle)
+
+            # According to p. 29 below Eq. 9
+            R = 0.2
+            n = 4
+
+            # Eq. 7 in Sun et al. (2017).
+            # CFL * Ma can be rewritten as Δt * v_max / h (see p. 29, right above Eq. 9).
+            delta_r_ = -dt * v_max * 2 * h * (1 + R * (kernel / Wdx)^n) *
+                       m_b / (rho_a + rho_b) * grad_kernel
+
+            # Write into the buffer
+            for i in eachindex(delta_r_)
+                @inbounds delta_r[i, particle] += delta_r_[i]
+            end
+        end
+    end
+
+    # Add δ_r from the buffer to the current coordinates
+    @threaded semi for particle in eachparticle(system)
+        for i in 1:ndims(system)
+            @inbounds u[i, particle] += delta_r[i, particle]
+        end
+    end
+
+    return u
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, typeof(particle_shifting!)})
+    @nospecialize cb # reduce precompilation time
+    print(io, "ParticleShiftingCallback()")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   cb::DiscreteCallback{<:Any, typeof(particle_shifting!)})
+    @nospecialize cb # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, cb)
+    else
+        summary_box(io, "ParticleShiftingCallback")
+    end
+end

src/general/system.jl

@@ -142,8 +142,7 @@ end
 system_smoothing_kernel(system) = system.smoothing_kernel
 system_correction(system) = nothing
 
-# TODO
-# @inline particle_spacing(system, particle) = system.initial_condition.particle_spacing
+@inline particle_spacing(system, particle) = system.initial_condition.particle_spacing
 
 # Assuming a constant particle spacing one can calculate the number of neighbors within the
 # compact support for an undisturbed particle distribution.

src/schemes/fluid/fluid.jl

@@ -45,17 +45,16 @@ function initial_smoothing_length(system, refinement)
            system.initial_condition.particle_spacing
 end
 
-# TODO
-# @inline function particle_spacing(system::FluidSystem, particle)
-#     return particle_spacing(system, system.particle_refinement, particle)
-# end
+@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 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
+@inline function particle_spacing(system, refinement, particle)
+    (; smoothing_length_factor) = system.cache
+    return smoothing_length(system, particle) / smoothing_length_factor
+end
 
 function write_u0!(u0, system::FluidSystem)
     (; initial_condition) = system

src/visualization/recipes_plots.jl

@@ -39,7 +39,10 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization, particle_
                                xlims=(-Inf, Inf), ylims=(-Inf, Inf))
     systems_data = map(enumerate(semi.systems)) do (i, system)
         u = wrap_u(u_ode, system, semi)
-        coordinates = active_coordinates(u, system)
+        periodic_box = get_neighborhood_search(system, semi).periodic_box
+        coordinates = PointNeighbors.periodic_coords(active_coordinates(u, system),
+                                                     periodic_box)
+
         x = collect(coordinates[1, :])
         y = collect(coordinates[2, :])
 

src/visualization/write2vtk.jl

@@ -124,9 +124,8 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre
         vtk["index"] = active_particles(system)
         vtk["time"] = t
 
-        # TODO
-        # vtk["particle_spacing"] = [particle_spacing(system, particle)
-        #                            for particle in active_particles(system)]
+        vtk["particle_spacing"] = [particle_spacing(system, particle)
+                                   for particle in active_particles(system)]
 
         if write_meta_data
             vtk["solver_version"] = git_hash
@@ -372,9 +371,8 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d
         vtk["lame_lambda"] = system.lame_lambda
         vtk["lame_mu"] = system.lame_mu
         vtk["smoothing_kernel"] = type2string(system.smoothing_kernel)
-        # TODO
-        # vtk["smoothing_length_factor"] = initial_smoothing_length(system) /
-        #                                  particle_spacing(system, 1)
+        vtk["smoothing_length_factor"] = initial_smoothing_length(system) /
+                                         particle_spacing(system, 1)
     end
 
     write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)