diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bcdce98ac..89026040c 100644 --- a/src/TrixiParticles.jl +++ b/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 diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..cf9df1b5e 100644 --- a/src/general/semidiscretization.jl +++ b/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 diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl new file mode 100644 index 000000000..bf886ac12 --- /dev/null +++ b/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 diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl new file mode 100644 index 000000000..eede0fc00 --- /dev/null +++ b/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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..94064555d 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/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)