Restructure TVF (#864)

* Restructure TVF

* Remove dv term from newer paper

* Eliminate time step from spatial discretization

* Rename updates

* Fix EDAC

* Fix `current_density`

* Fix sign of background pressure

* Add reference

* Implement suggestions

* Remove `update_from_callback`

* Implement suggestions

* Fix tests

* Fix background pressure factor

* Improve comments

* Fix boundary bug

* Make N-Body system type-stable

* Fix packing

* Fix update callback

* Update src/schemes/fluid/transport_velocity.jl

Co-authored-by: Sven Berger <berger.sven@gmail.com>

---------

Co-authored-by: Sven Berger <berger.sven@gmail.com>

Author: efaulhaber

examples/n_body/n_body_system.jl

@@ -1,8 +1,8 @@
 using TrixiParticles
 using LinearAlgebra
 
-struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS}
-    initial_condition :: InitialCondition{ELTYPE}
+struct NBodySystem{NDIMS, ELTYPE <: Real, IC} <: TrixiParticles.System{NDIMS}
+    initial_condition :: IC
     mass              :: Array{ELTYPE, 1} # [particle]
     G                 :: ELTYPE
     buffer            :: Nothing
@@ -11,13 +11,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS}
         mass = copy(initial_condition.mass)
 
         new{size(initial_condition.coordinates, 1),
-            eltype(mass)}(initial_condition, mass, G, nothing)
+            eltype(mass), typeof(initial_condition)}(initial_condition, mass, G, nothing)
     end
 end
 
 TrixiParticles.timer_name(::NBodySystem) = "nbody"
 
-@inline Base.eltype(system::NBodySystem) = eltype(system.initial_condition.coordinates)
+@inline Base.eltype(system::NBodySystem{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE
 
 function TrixiParticles.write_u0!(u0, system::NBodySystem)
     u0 .= system.initial_condition.coordinates

src/callbacks/density_reinit.jl

@@ -69,7 +69,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr
         # Update systems to compute quantities like density and pressure.
         semi = integrator.p
         v_ode, u_ode = u.x
-        update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
+        update_systems_and_nhs(v_ode, u_ode, semi, t)
 
         # Apply the callback.
         cb(integrator)

src/callbacks/particle_shifting.jl

@@ -39,8 +39,7 @@ function particle_shifting!(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                     update_from_callback=true)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
         end
 
         @trixi_timeit timer() "particle shifting" foreach_system(semi) do system

src/callbacks/post_process.jl

@@ -240,8 +240,7 @@ function (pp::PostprocessCallback)(integrator)
         # still have the values from the last stage of the previous step if not updated here.
         @trixi_timeit timer() "update systems and nhs" begin
             # Don't create sub-timers here to avoid cluttering the timer output
-            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                     update_from_callback=true)
+            @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
         end
 
         foreach_system(semi) do system

src/callbacks/update.jl

@@ -57,7 +57,7 @@ function initial_update!(cb::UpdateCallback, u, t, integrator)
 
     # Tell systems that `UpdateCallback` is used
     foreach_system(semi) do system
-        update_callback_used!(system)
+        set_callback_flag!(system, true)
     end
 
     return cb(integrator)
@@ -78,7 +78,7 @@ function (update_callback!::UpdateCallback)(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)
+    update_systems_and_nhs(v_ode, u_ode, semi, t)
 
     # Update open boundaries first, since particles might be activated or deactivated
     @trixi_timeit timer() "update open boundary" foreach_system(semi) do system
@@ -89,6 +89,7 @@ function (update_callback!::UpdateCallback)(integrator)
         update_particle_packing(system, v_ode, u_ode, semi, integrator)
     end
 
+    # This is only used by the particle packing system and should be removed in the future
     @trixi_timeit timer() "update TVF" foreach_system(semi) do system
         update_transport_velocity!(system, v_ode, semi)
     end
@@ -141,5 +142,3 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "UpdateCallback", setup)
     end
 end
-
-update_callback_used!(system) = system

src/general/semidiscretization.jl

@@ -338,8 +338,8 @@ function semidiscretize(semi, tspan; reset_threads=true)
         # Initialize this system
         initialize!(system, semi_new)
 
-        # Only for systems requiring a mandatory callback
-        reset_callback_flag!(system)
+        # Only for systems requiring the use of the `UpdateCallback`
+        set_callback_flag!(system, false)
     end
 
     return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new)
@@ -373,8 +373,8 @@ function restart_with!(semi, sol; reset_threads=true)
 
         restart_with!(system, v, u)
 
-        # Only for systems requiring a mandatory callback
-        reset_callback_flag!(system)
+        # Only for systems requiring the use of the `UpdateCallback`
+        set_callback_flag!(system, false)
     end
 
     return semi
@@ -464,17 +464,33 @@ function drift!(du_ode, v_ode, u_ode, semi, t)
 end
 
 @inline function add_velocity!(du, v, particle, system)
+    # Generic fallback for all systems that don't define this function
     for i in 1:ndims(system)
-        du[i, particle] = v[i, particle]
+        @inbounds du[i, particle] = v[i, particle]
     end
 
     return du
 end
 
+# Solid wall boundary system doesn't integrate the particle positions
 @inline add_velocity!(du, v, particle, system::BoundarySPHSystem) = du
 
+@inline function add_velocity!(du, v, particle, system::FluidSystem)
+    # This is zero unless a transport velocity is used
+    delta_v_ = delta_v(system, particle)
+
+    for i in 1:ndims(system)
+        @inbounds du[i, particle] = v[i, particle] + delta_v_[i]
+    end
+
+    return du
+end
+
 function kick!(dv_ode, v_ode, u_ode, semi, t)
     @trixi_timeit timer() "kick!" begin
+        # Check that the `UpdateCallback` is used if required
+        check_update_callback(semi)
+
         @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode)
 
         @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode,
@@ -490,8 +506,9 @@ function kick!(dv_ode, v_ode, u_ode, semi, t)
     return dv_ode
 end
 
-# Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces
-function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false)
+# Update the systems and neighborhood searches (NHS) for a simulation
+# before calling `interact!` to compute forces.
+function update_systems_and_nhs(v_ode, u_ode, semi, t)
     # First update step before updating the NHS
     # (for example for writing the current coordinates in the solid system)
     foreach_system(semi) do system
@@ -523,12 +540,21 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals
         update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     end
 
+    # This update depends on the computed quantities of the fluid system and therefore
+    # needs to be after `update_quantities!`.
+    foreach_system(semi) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+
+        update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
+    end
+
     # Final update step for all remaining systems
     foreach_system(semi) do system
         v = wrap_v(v_ode, system, semi)
         u = wrap_u(u_ode, system, semi)
 
-        update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback)
+        update_final!(system, v, u, v_ode, u_ode, semi, t)
     end
 end
 
@@ -858,6 +884,16 @@ function update!(neighborhood_search, x, y, semi; points_moving=(true, false),
                            parallelization_backend=semi.parallelization_backend)
 end
 
+function check_update_callback(semi)
+    foreach_system(semi) do system
+        # This check will be optimized away if the system does not require the callback
+        if requires_update_callback(system) && !update_callback_used(system)
+            system_name = system |> typeof |> nameof
+            throw(ArgumentError("`UpdateCallback` is required for `$system_name`"))
+        end
+    end
+end
+
 function check_configuration(systems,
                              nhs::Union{Nothing, PointNeighbors.AbstractNeighborhoodSearch})
     foreach_system(systems) do system

src/general/system.jl

@@ -126,7 +126,7 @@ end
                                  system_correction(system), system, particle)
 end
 
-# System update orders. This can be dispatched if needed.
+# System updates do nothing by default, but can be dispatched if needed
 function update_positions!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
@@ -139,21 +139,27 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false)
+function update_boundary_interpolation!(system, v, u, v_ode, u_ode, semi, t)
     return system
 end
 
-# Only for systems requiring a mandatory callback
-reset_callback_flag!(system) = system
+function update_final!(system, v, u, v_ode, u_ode, semi, t)
+    return system
+end
+
+# Only for systems requiring the use of the `UpdateCallback`
+@inline requires_update_callback(system) = false
+@inline update_callback_used(system) = false
+@inline set_callback_flag!(system, value) = system
 
-initial_smoothing_length(system) = smoothing_length(system, nothing)
+@inline initial_smoothing_length(system) = smoothing_length(system, nothing)
 
-function smoothing_length(system, particle)
+@inline function smoothing_length(system, particle)
     return system.smoothing_length
 end
 
-system_smoothing_kernel(system) = system.smoothing_kernel
-system_correction(system) = nothing
+@inline system_smoothing_kernel(system) = system.smoothing_kernel
+@inline system_correction(system) = nothing
 
 @inline particle_spacing(system, particle) = system.initial_condition.particle_spacing
 

src/io/write_vtk.jl

@@ -58,8 +58,7 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix
     # still have the values from the last stage of the previous step if not updated here.
     @trixi_timeit timer() "update systems" begin
         # Don't create sub-timers here to avoid cluttering the timer output
-        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t;
-                                                 update_from_callback=true)
+        @notimeit timer() update_systems_and_nhs(v_ode, u_ode, semi, t)
     end
 
     filenames = system_names(systems)

src/preprocessing/particle_packing/system.jl

@@ -215,14 +215,15 @@ end
     return ndims(system)
 end
 
-function reset_callback_flag!(system::ParticlePackingSystem)
-    system.update_callback_used[] = false
+@inline requires_update_callback(system::ParticlePackingSystem) = true
+@inline update_callback_used(system::ParticlePackingSystem) = system.update_callback_used[]
+
+function set_callback_flag!(system::ParticlePackingSystem, value)
+    system.update_callback_used[] = value
 
     return system
 end
 
-update_callback_used!(system::ParticlePackingSystem) = system.update_callback_used[] = true
-
 function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true)
     vtk["velocity"] = [advection_velocity(v, system, particle)
                        for particle in active_particles(system)]
@@ -289,15 +290,6 @@ function update_position!(u, system::ParticlePackingSystem, semi)
     return u
 end
 
-function update_final!(system::ParticlePackingSystem, v, u, v_ode, u_ode, semi, t;
-                       update_from_callback=false)
-    if !update_from_callback && !(system.update_callback_used[])
-        throw(ArgumentError("`UpdateCallback` is required when using `ParticlePackingSystem`"))
-    end
-
-    return system
-end
-
 # Skip for systems without `SignedDistanceField`
 constrain_particles_onto_surface!(u, system::ParticlePackingSystem{Nothing}, semi) = u
 
@@ -398,6 +390,10 @@ end
     return system
 end
 
+@inline function update_transport_velocity!(system, v_ode, semi)
+    return system
+end
+
 # Skip for fixed systems
 @inline add_velocity!(du, v, particle, system::ParticlePackingSystem{<:Any, true}) = du
 

src/schemes/boundary/dummy_particles/dummy_particles.jl

@@ -311,17 +311,16 @@ end
                                   system, v, u, v_ode, u_ode, semi)
     (; correction, density_calculator) = boundary_model
 
-    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
+    compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
 
+    # These are only computed when using corrections
+    compute_correction_values!(system, correction, u, v_ode, u_ode, semi)
     compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode,
                                         semi)
-
     # `kernel_correct_density!` only performed for `SummationDensity`
     kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction,
                             density_calculator)
 
-    compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi)
-
     return boundary_model
 end