Merge branch 'main' into dev

Author: efaulhaber

src/callbacks/update.jl

@@ -55,10 +55,8 @@ end
 function initial_update!(cb::UpdateCallback, u, t, integrator)
     semi = integrator.p
 
-    # Tell systems that `UpdateCallback` is used
-    foreach_system(semi) do system
-        set_callback_flag!(system, true)
-    end
+    # Tell the semidiscretization that the `UpdateCallback` is used
+    semi.update_callback_used[] = true
 
     return cb(integrator)
 end

src/general/semidiscretization.jl

@@ -47,23 +47,26 @@ semi = Semidiscretization(fluid_system, boundary_system,
 └──────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```
 """
-struct Semidiscretization{BACKEND, S, RU, RV, NS}
+struct Semidiscretization{BACKEND, S, RU, RV, NS, UCU}
     systems                 :: S
     ranges_u                :: RU
     ranges_v                :: RV
     neighborhood_searches   :: NS
     parallelization_backend :: BACKEND
+    update_callback_used    :: UCU
 
     # Dispatch at `systems` to distinguish this constructor from the one below when
     # 4 systems are passed.
     # This is an internal constructor only used in `test/count_allocations.jl`
     # and by Adapt.jl.
     function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches,
-                                parallelization_backend::PointNeighbors.ParallelizationBackend)
+                                parallelization_backend::PointNeighbors.ParallelizationBackend,
+                                update_callback_used)
         new{typeof(parallelization_backend), typeof(systems), typeof(ranges_u),
-            typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v,
-                                                             neighborhood_searches,
-                                                             parallelization_backend)
+            typeof(ranges_v), typeof(neighborhood_searches),
+            typeof(update_callback_used)}(systems, ranges_u, ranges_v,
+                                          neighborhood_searches, parallelization_backend,
+                                          update_callback_used)
     end
 end
 
@@ -92,8 +95,13 @@ function Semidiscretization(systems::Union{System, Nothing}...;
                            for neighbor in systems)
                      for system in systems)
 
+    # These will be set to true inside the `UpdateCallback`.
+    # Some techniques require the use of this callback, and this flag can be used
+    # to determine if the callback is used in a simulation.
+    update_callback_used = Ref(false)
+
     return Semidiscretization(systems, ranges_u, ranges_v, searches,
-                              parallelization_backend)
+                              parallelization_backend, update_callback_used)
 end
 
 # Inline show function e.g. Semidiscretization(neighborhood_search=...)
@@ -328,7 +336,8 @@ function semidiscretize(semi, tspan; reset_threads=true)
         semi_new = Semidiscretization(set_system_links.(semi_.systems, Ref(semi_)),
                                       semi_.ranges_u, semi_.ranges_v,
                                       semi_.neighborhood_searches,
-                                      semi_.parallelization_backend)
+                                      semi_.parallelization_backend,
+                                      semi_.update_callback_used)
     else
         semi_new = semi
     end
@@ -337,11 +346,11 @@ function semidiscretize(semi, tspan; reset_threads=true)
     foreach_system(semi_new) do system
         # Initialize this system
         initialize!(system, semi_new)
-
-        # Only for systems requiring the use of the `UpdateCallback`
-        set_callback_flag!(system, false)
     end
 
+    # Reset callback flag that will be set by the `UpdateCallback`
+    semi_new.update_callback_used[] = false
+
     return DynamicalODEProblem(kick!, drift!, v0_ode, u0_ode, tspan, semi_new)
 end
 
@@ -372,11 +381,11 @@ function restart_with!(semi, sol; reset_threads=true)
         u = wrap_u(sol.u[end].x[2], system, semi)
 
         restart_with!(system, v, u)
-
-        # Only for systems requiring the use of the `UpdateCallback`
-        set_callback_flag!(system, false)
     end
 
+    # Reset callback flag that will be set by the `UpdateCallback`
+    semi.update_callback_used[] = false
+
     return semi
 end
 
@@ -887,7 +896,7 @@ 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)
+        if requires_update_callback(system) && !semi.update_callback_used[]
             system_name = system |> typeof |> nameof
             throw(ArgumentError("`UpdateCallback` is required for `$system_name`"))
         end
@@ -1007,6 +1016,5 @@ function set_system_links(system::OpenBoundarySPHSystem, semi)
                                  system.reference_pressure,
                                  system.reference_density,
                                  system.buffer,
-                                 system.update_callback_used,
                                  system.cache)
 end

src/general/system.jl

@@ -149,8 +149,6 @@ 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
 
 @inline initial_smoothing_length(system) = smoothing_length(system, nothing)
 

src/preprocessing/particle_packing/system.jl

@@ -56,7 +56,7 @@ For more information on the methods, see [particle packing](@ref particle_packin
                               Recommended values are `0.8` or `0.9`.
 """
 struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
-                             IC, M, D, K, N, SD, UCU} <: FluidSystem{NDIMS}
+                             IC, M, D, K, N, SD} <: FluidSystem{NDIMS}
     initial_condition              :: IC
     advection_velocity             :: AV
     mass                           :: M
@@ -73,7 +73,6 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
     signed_distances               :: SD # Only for visualization
     particle_refinement            :: PR
     buffer                         :: Nothing
-    update_callback_used           :: UCU
     cache                          :: C
 
     # This constructor is necessary for Adapt.jl to work with this struct.
@@ -84,22 +83,20 @@ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV,
                                    signed_distance_field,
                                    is_boundary, shift_length, neighborhood_search,
                                    signed_distances, particle_refinement, buffer,
-                                   update_callback_used, fixed_system, cache,
-                                   advection_velocity)
+                                   fixed_system, cache, advection_velocity)
         return new{typeof(signed_distance_field), fixed_system, ndims(smoothing_kernel),
                    eltype(density), typeof(particle_refinement), typeof(cache),
                    typeof(advection_velocity), typeof(initial_condition), typeof(mass),
                    typeof(density), typeof(smoothing_kernel), typeof(neighborhood_search),
-                   typeof(signed_distances),
-                   typeof(update_callback_used)}(initial_condition, advection_velocity,
-                                                 mass, density, particle_spacing,
-                                                 smoothing_kernel,
-                                                 smoothing_length_interpolation,
-                                                 background_pressure, place_on_shell,
-                                                 signed_distance_field, is_boundary,
-                                                 shift_length, neighborhood_search,
-                                                 signed_distances, particle_refinement,
-                                                 buffer, update_callback_used, cache)
+                   typeof(signed_distances)}(initial_condition, advection_velocity,
+                                             mass, density, particle_spacing,
+                                             smoothing_kernel,
+                                             smoothing_length_interpolation,
+                                             background_pressure, place_on_shell,
+                                             signed_distance_field, is_boundary,
+                                             shift_length, neighborhood_search,
+                                             signed_distances, particle_refinement,
+                                             buffer, cache)
     end
 end
 
@@ -167,8 +164,7 @@ function ParticlePackingSystem(shape::InitialCondition;
                                  background_pressure, place_on_shell, signed_distance_field,
                                  is_boundary, shift_length, nhs,
                                  fill(zero(ELTYPE), nparticles(shape)), particle_refinement,
-                                 nothing, Ref(false), fixed_system, cache,
-                                 advection_velocity)
+                                 nothing, fixed_system, cache, advection_velocity)
 end
 
 function Base.show(io::IO, system::ParticlePackingSystem)
@@ -219,13 +215,6 @@ end
 end
 
 @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
 
 function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true)
     vtk["velocity"] = [advection_velocity(v, system, particle)

src/schemes/boundary/open_boundary/system.jl

@@ -37,47 +37,44 @@ Open boundary system for in- and outflow particles.
     It is GPU-compatible (e.g., with CUDA.jl and AMDGPU.jl), but currently **not** supported with Metal.jl.
 """
 struct OpenBoundarySPHSystem{BM, ELTYPE, NDIMS, IC, FS, FSI, ARRAY1D, BC, FC, BZ, RV,
-                             RP, RD, B, UCU, C} <: System{NDIMS}
-    boundary_model       :: BM
-    initial_condition    :: IC
-    fluid_system         :: FS
-    fluid_system_index   :: FSI
-    smoothing_length     :: ELTYPE
-    mass                 :: ARRAY1D # Array{ELTYPE, 1}: [particle]
-    density              :: ARRAY1D # Array{ELTYPE, 1}: [particle]
-    volume               :: ARRAY1D # Array{ELTYPE, 1}: [particle]
-    pressure             :: ARRAY1D # Array{ELTYPE, 1}: [particle]
-    boundary_candidates  :: BC      # Array{UInt32, 1}: [particle]
-    fluid_candidates     :: FC      # Array{UInt32, 1}: [particle]
-    boundary_zone        :: BZ
-    reference_velocity   :: RV
-    reference_pressure   :: RP
-    reference_density    :: RD
-    buffer               :: B
-    update_callback_used :: UCU
-    cache                :: C
+                             RP, RD, B, C} <: System{NDIMS}
+    boundary_model      :: BM
+    initial_condition   :: IC
+    fluid_system        :: FS
+    fluid_system_index  :: FSI
+    smoothing_length    :: ELTYPE
+    mass                :: ARRAY1D # Array{ELTYPE, 1}: [particle]
+    density             :: ARRAY1D # Array{ELTYPE, 1}: [particle]
+    volume              :: ARRAY1D # Array{ELTYPE, 1}: [particle]
+    pressure            :: ARRAY1D # Array{ELTYPE, 1}: [particle]
+    boundary_candidates :: BC      # Array{UInt32, 1}: [particle]
+    fluid_candidates    :: FC      # Array{UInt32, 1}: [particle]
+    boundary_zone       :: BZ
+    reference_velocity  :: RV
+    reference_pressure  :: RP
+    reference_density   :: RD
+    buffer              :: B
+    cache               :: C
 end
 
 function OpenBoundarySPHSystem(boundary_model, initial_condition, fluid_system,
                                fluid_system_index, smoothing_length, mass, density, volume,
                                pressure, boundary_candidates, fluid_candidates,
                                boundary_zone, reference_velocity,
-                               reference_pressure, reference_density, buffer,
-                               update_callback_used, cache)
+                               reference_pressure, reference_density, buffer, cache)
     OpenBoundarySPHSystem{typeof(boundary_model), eltype(mass), ndims(initial_condition),
                           typeof(initial_condition), typeof(fluid_system),
                           typeof(fluid_system_index), typeof(mass),
                           typeof(boundary_candidates), typeof(fluid_candidates),
                           typeof(boundary_zone), typeof(reference_velocity),
                           typeof(reference_pressure), typeof(reference_density),
-                          typeof(buffer), typeof(update_callback_used),
+                          typeof(buffer),
                           typeof(cache)}(boundary_model, initial_condition, fluid_system,
                                          fluid_system_index, smoothing_length, mass,
                                          density, volume, pressure, boundary_candidates,
                                          fluid_candidates, boundary_zone,
                                          reference_velocity, reference_pressure,
-                                         reference_density, buffer, update_callback_used,
-                                         cache)
+                                         reference_density, buffer, cache)
 end
 
 function OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
@@ -149,8 +146,6 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
                                        reference_density, reference_velocity,
                                        reference_pressure)
 
-    # These will be set later
-    update_callback_used = Ref(false)
     fluid_system_index = Ref(0)
 
     smoothing_length = initial_smoothing_length(fluid_system)
@@ -162,8 +157,7 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone;
                                  fluid_system_index, smoothing_length, mass, density,
                                  volume, pressure, boundary_candidates, fluid_candidates,
                                  boundary_zone, reference_velocity_,
-                                 reference_pressure_, reference_density_, buffer,
-                                 update_callback_used, cache)
+                                 reference_pressure_, reference_density_, buffer, cache)
 end
 
 function create_cache_open_boundary(boundary_model, initial_condition,
@@ -226,14 +220,8 @@ end
     return ELTYPE
 end
 
+# The `UpdateCallback` is required to update particle positions between time steps
 @inline requires_update_callback(system::OpenBoundarySPHSystem) = true
-@inline update_callback_used(system::OpenBoundarySPHSystem) = system.update_callback_used[]
-
-function set_callback_flag!(system::OpenBoundarySPHSystem, value)
-    system.update_callback_used[] = value
-
-    return system
-end
 
 function corresponding_fluid_system(system::OpenBoundarySPHSystem, semi)
     return system.fluid_system

src/schemes/fluid/entropically_damped_sph/rhs.jl

@@ -27,10 +27,12 @@ function interact!(dv, v_particle_system, u_particle_system,
         p_a = current_pressure(v_particle_system, particle_system, particle)
         p_b = current_pressure(v_neighbor_system, neighbor_system, neighbor)
 
-        # This technique is for a more robust `pressure_acceleration` but only with TVF.
-        # It results only in significant improvement for EDAC and not for WCSPH.
-        # See Ramachandran (2019) p. 582
-        # Note that the return value is zero when not using EDAC with TVF.
+        # This technique by Basa et al. 2017 (10.1002/fld.1927) aims to reduce numerical
+        # errors due to large pressures by subtracting the average pressure of neighboring
+        # particles.
+        # It results in significant improvement for EDAC, especially with TVF,
+        # but not for WCSPH, according to Ramachandran & Puri (2019), Section 3.2.
+        # Note that the return value is zero when not using average pressure reduction.
         p_avg = average_pressure(particle_system, particle)
 
         m_a = hydrodynamic_mass(particle_system, particle)