Adaptive Particle Refinement: Resize semi v2

src/TrixiParticles.jl

@@ -51,6 +51,7 @@ include("general/semidiscretization.jl")
 include("general/gpu.jl")
 include("visualization/write2vtk.jl")
 include("visualization/recipes_plots.jl")
+include("refinement/resize.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition

src/general/semidiscretization.jl

@@ -73,14 +73,7 @@ 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))
+    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 +85,16 @@ function Semidiscretization(systems...;
     return Semidiscretization(systems, ranges_u, ranges_v, searches)
 end
 
+function ranges_vu(systems)
+    sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems]
+    ranges_u = [(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 = [(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)]
+
+    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

src/refinement/refinement.jl

@@ -0,0 +1,9 @@
+struct ParticleRefinement
+    n_particles_before_resize :: Ref{Int}
+    n_new_particles           :: Ref{Int}
+    delete_candidates         :: Vector{Bool}   # length(delete_candidates) == nparticles
+end
+
+function ParticleRefinement()
+    return ParticleRefinement(Ref(0), Ref(0), Bool[])
+end

src/refinement/resize.jl

@@ -0,0 +1,178 @@
+function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode)
+    # Resize all systems
+    foreach_system(semi) do system
+        capacity(system) > nparticles(system) && resize!(system, capacity(system))
+    end
+
+    resize!(v_ode, u_ode, _v_ode, _u_ode, semi)
+
+    return semi
+end
+
+function Base.deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode)
+    # Delete at specific indices
+    foreach_system(semi) do system
+        v = wrap_v(v_ode, system, semi)
+        u = wrap_u(u_ode, system, semi)
+        deleteat!(system, v, u)
+    end
+
+    resize!(v_ode, u_ode, _v_ode, _u_ode, semi)
+
+    return semi
+end
+
+function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization)
+    (; systems) = semi
+
+    copyto!(_v_ode, v_ode)
+    copyto!(_u_ode, u_ode)
+
+    # Get ranges after resizing the systems
+    ranges_v_new, ranges_u_new = ranges_vu(systems)
+
+    ranges_v_old = copy(semi.ranges_v)
+    ranges_u_old = copy(semi.ranges_u)
+
+    # Set ranges after resizing the systems
+    for i in 1:length(systems)
+        semi.ranges_v[i] = ranges_v_new[i]
+        semi.ranges_u[i] = ranges_u_new[i]
+    end
+
+    for i in eachindex(ranges_u_old)
+        length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i]))
+        for j in 0:(length_u - 1)
+            u_ode[ranges_u_new[i][1] + j] = _u_ode[ranges_u_old[i][1] + j]
+        end
+
+        length_v = min(length(ranges_v_old[i]), length(ranges_v_new[i]))
+        for j in 0:(length_v - 1)
+            v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j]
+        end
+    end
+
+    sizes_u = sum(u_nvariables(system) * n_moving_particles(system) for system in systems)
+    sizes_v = sum(v_nvariables(system) * n_moving_particles(system) for system in systems)
+
+    resize!(v_ode, sizes_v)
+    resize!(u_ode, sizes_u)
+
+    resize!(_v_ode, sizes_v)
+    resize!(_u_ode, sizes_u)
+
+    # TODO: Do the following in the callback
+    # resize!(integrator, (length(v_ode), length(u_ode)))
+
+    # # Tell OrdinaryDiffEq that u has been modified
+    # u_modified!(integrator, true)
+    return v_ode
+end
+
+Base.resize!(system::System, capacity_system) = system
+
+function Base.resize!(system::FluidSystem, capacity_system)
+    return resize!(system, system.particle_refinement, capacity_system)
+end
+
+Base.resize!(system, ::Nothing, capacity_system) = system
+
+function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int)
+    (; mass, pressure, cache, density_calculator) = system
+
+    refinement.n_particles_before_resize[] = nparticles(system)
+
+    resize!(mass, capacity_system)
+    resize!(pressure, capacity_system)
+    resize_density!(system, capacity_system, density_calculator)
+    # TODO
+    # resize_cache!(system, cache, n)
+end
+
+function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int)
+    (; mass, cache, density_calculator) = system
+
+    refinement.n_particles_before_resize[] = nparticles(system)
+
+    resize!(mass, capacity_system)
+    resize_density!(system, capacity_system, density_calculator)
+    # TODO
+    # resize_cache!(system, capacity_system)
+
+    return system
+end
+
+resize_density!(system, n::Int, ::SummationDensity) = resize!(system.cache.density, n)
+resize_density!(system, n::Int, ::ContinuityDensity) = system
+
+function resize_cache!(system, n::Int)
+    resize!(system.cache.smoothing_length, n)
+
+    return system
+end
+
+function resize_cache!(system::EntropicallyDampedSPHSystem, n)
+    resize!(system.cache.smoothing_length, n)
+    resize!(system.cache.pressure_average, n)
+    resize!(system.cache.neighbor_counter, n)
+
+    return system
+end
+
+Base.deleteat!(system::System, v, u) = system
+
+function Base.deleteat!(system::FluidSystem, v, u)
+    return deleteat!(system, system.particle_refinement, v, u)
+end
+
+Base.deleteat!(system::FluidSystem, ::Nothing, v, u) = system
+
+function Base.deleteat!(system::FluidSystem, refinement, v, u)
+    (; delete_candidates) = refinement
+
+    delete_counter = 0
+
+    for particle in eachparticle(system)
+        if !iszero(delete_candidates[particle])
+            # swap particles (keep -> delete)
+            dump_id = nparticles(system) - delete_counter
+
+            vel_keep = current_velocity(v, system, dump_id)
+            pos_keep = current_coords(u, system, dump_id)
+
+            mass_keep = hydrodynamic_mass(system, dump_id)
+            density_keep = particle_density(v, system, dump_id)
+            pressure_keep = particle_pressure(v, system, dump_id)
+            #TODO
+            # smoothing_length_keep = smoothing_length(system, dump_id)
+            # system.cache.smoothing_length[particle] = smoothing_length_keep
+
+            system.mass[particle] = mass_keep
+
+            set_particle_pressure!(v, system, particle, pressure_keep)
+
+            set_particle_density!(v, system, particle, density_keep)
+
+            for dim in 1:ndims(system)
+                v[dim, particle] = vel_keep[dim]
+                u[dim, particle] = pos_keep[dim]
+            end
+
+            delete_counter += 1
+        end
+    end
+
+    resize!(system, nparticles(system) - delete_counter)
+
+    return system
+end
+
+@inline capacity(system) = nparticles(system)
+
+@inline capacity(system::FluidSystem) = capacity(system, system.particle_refinement)
+
+@inline capacity(system, ::Nothing) = nparticles(system)
+
+@inline function capacity(system, particle_refinement)
+    return particle_refinement.n_new_particles[] + nparticles(system)
+end

src/schemes/fluid/entropically_damped_sph/system.jl

@@ -44,7 +44,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more
                     gravity-like source terms.
 """
 struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
-                                   PF, ST, B, C} <: FluidSystem{NDIMS, IC}
+                                   PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC}
     initial_condition                 :: IC
     mass                              :: M # Vector{ELTYPE}: [particle]
     density_calculator                :: DC
@@ -59,55 +59,55 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV,
     transport_velocity                :: TV
     source_terms                      :: ST
     buffer                            :: B
+    particle_refinement               :: PR
     cache                             :: C
+end
 
-    function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
-                                         smoothing_length, sound_speed;
-                                         pressure_acceleration=inter_particle_averaged_pressure,
-                                         density_calculator=SummationDensity(),
-                                         transport_velocity=nothing,
-                                         alpha=0.5, viscosity=nothing,
-                                         acceleration=ntuple(_ -> 0.0,
-                                                             ndims(smoothing_kernel)),
-                                         source_terms=nothing, buffer_size=nothing)
-        buffer = isnothing(buffer_size) ? nothing :
-                 SystemBuffer(nparticles(initial_condition), buffer_size)
+function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel,
+                                     smoothing_length, sound_speed;
+                                     pressure_acceleration=inter_particle_averaged_pressure,
+                                     density_calculator=SummationDensity(),
+                                     transport_velocity=nothing,
+                                     alpha=0.5, viscosity=nothing,
+                                     acceleration=ntuple(_ -> 0.0,
+                                                         ndims(smoothing_kernel)),
+                                     particle_refinement=nothing,
+                                     source_terms=nothing, buffer_size=nothing)
+    buffer = isnothing(buffer_size) ? nothing :
+             SystemBuffer(nparticles(initial_condition), buffer_size)
 
-        initial_condition = allocate_buffer(initial_condition, buffer)
+    initial_condition = allocate_buffer(initial_condition, buffer)
 
-        NDIMS = ndims(initial_condition)
-        ELTYPE = eltype(initial_condition)
+    NDIMS = ndims(initial_condition)
+    ELTYPE = eltype(initial_condition)
 
-        mass = copy(initial_condition.mass)
+    mass = copy(initial_condition.mass)
 
-        if ndims(smoothing_kernel) != NDIMS
-            throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
-        end
+    if ndims(smoothing_kernel) != NDIMS
+        throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem"))
+    end
 
-        acceleration_ = SVector(acceleration...)
-        if length(acceleration_) != NDIMS
-            throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
-        end
+    acceleration_ = SVector(acceleration...)
+    if length(acceleration_) != NDIMS
+        throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem"))
+    end
 
-        pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
-                                                                         density_calculator,
-                                                                         NDIMS, ELTYPE,
-                                                                         nothing)
+    pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration,
+                                                                     density_calculator,
+                                                                     NDIMS, ELTYPE,
+                                                                     nothing)
 
-        nu_edac = (alpha * smoothing_length * sound_speed) / 8
+    nu_edac = (alpha * smoothing_length * sound_speed) /
+              (2 * ndims(initial_condition) + 4)
 
-        cache = create_cache_density(initial_condition, density_calculator)
-        cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
+    cache = create_cache_density(initial_condition, density_calculator)
+    cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...)
 
-        new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass),
-            typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity),
-            typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms),
-            typeof(buffer),
-            typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel,
-                           smoothing_length, sound_speed, viscosity, nu_edac, acceleration_,
-                           nothing, pressure_acceleration, transport_velocity, source_terms,
-                           buffer, cache)
-    end
+    return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator,
+                                       smoothing_kernel, smoothing_length, sound_speed,
+                                       viscosity, nu_edac, acceleration_, nothing,
+                                       pressure_acceleration, transport_velocity,
+                                       source_terms, buffer, particle_refinement, cache)
 end
 
 function Base.show(io::IO, system::EntropicallyDampedSPHSystem)