From ad5be9005404d425a84481daba673d4d9b14c234 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Mar 2024 16:57:40 +0100 Subject: [PATCH 01/20] first non-working prototype --- src/general/refinement.jl | 214 ++++++++++++++++++ src/general/semidiscretization.jl | 36 ++- src/general/system.jl | 2 + .../fluid/weakly_compressible_sph/system.jl | 17 +- 4 files changed, 251 insertions(+), 18 deletions(-) create mode 100644 src/general/refinement.jl diff --git a/src/general/refinement.jl b/src/general/refinement.jl new file mode 100644 index 000000000..562af9d3e --- /dev/null +++ b/src/general/refinement.jl @@ -0,0 +1,214 @@ +# Criteria of refinement: +# +# - fixed (/moving?) refinement zone +# - number of neighbors +# - problem specific criteria (e.g. high velocity gradient) + +include("refinement_pattern.jl") + +mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} + candidates :: Vector{ELTYPE} + refinement_levels :: Int + refinement_pattern :: RP + refinement_criteria :: RC + + # Depends on refinement pattern, particle spacing and parameters ϵ and α. + # Should be obtained prior to simulation in `create_system_child()` + rel_position_childs :: Vector{SVector{NDIMS, ELTYPE}} + mass_child :: ELTYPE + available_childs :: Int + + # It is essential to know the child system, which is empty at the beginning + # and will be created in `create_system_child()` + system_child::System + + # API --> parent system with `RL=0` + function ParticleRefinement(refinement_criteria...; + refinement_pattern=CubicSplitting(), + refinement_levels=1) + ELTYPE = eltype(refinement_criteria) + NDIMS = ndims(refinement_criteria[1]) + + return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), + typeof(refinement_criteria)}([], refinement_levels, refinement_pattern, + refinement_criteria) + end + + # Internal constructor for multiple refinement levels + function ParticleRefinement{RL}(refinement_criteria::Tuple, + refinement_pattern, refinement_levels) where {RL} + ELTYPE = eltype(refinement_criteria) + NDIMS = ndims(refinement_criteria[1]) + + return new{RL, NDIMS, ELTYPE, typeof(refinement_pattern), + typeof(refinement_criteria)}([], refinement_levels, refinement_pattern, + refinement_criteria) + end +end + +@inline Base.ndims(::ParticleRefinement{RL, NDIMS}) where {RL, NDIMS} = ndims + +@inline function child_set(system, particle_refinement) + (; available_childs) = particle_refinement + + start_iter = nparticles(system) + 1 - available_childs + end_iter = start_iter + nchilds(system, particle_refinement) + + return start_iter:end_iter +end + +# ==== Create child systems + +function create_system_childs(systems) + systems_ = () + foreach_system(systems) do system + systems_ = (systems_..., create_system_child(system, system.particle_refinement)...) + end + + return (systems..., systems_...) +end + +create_system_child(system, ::Nothing) = () + +function create_system_child(system::WeaklyCompressibleSPHSystem, + particle_refinement::ParticleRefinement{RL}) where {RL} + (; refinement_levels, refinement_pattern, refinement_criteria) = particle_refinement + (; density_calculator, state_equation, smoothing_kernel, + pressure_acceleration_formulation, viscosity, density_diffusion, + acceleration, correction, source_terms) = system + + NDIMS = ndims(system) + + # Distribute values according to refinement pattern + smoothing_length_ = smoothing_length_child(system, refinement_pattern) + + # Create "empty" `InitialCondition` for child system + particle_spacing_ = particle_spacing_child(system, refinement_pattern) + coordinates_ = zeros(NDIMS, 2) + velocity_ = similar(coordinates_) + density_ = system.initial_condition.density[1] + pressure_ = system.initial_condition.pressure[1] + mass_ = nothing + + empty_ic = InitialCondition{NDIMS}(coordinates_, velocity_, mass_, density_, pressure_, + particle_spacing_) + + # Let recursive dispatch handle multiple refinement levels + level = RL + 1 + particle_refinement_ = if level == refinement_levels + nothing + else + ParticleRefinement{level}(refinement_criteria, refinement_pattern, + refinement_levels) + end + + system_child = WeaklyCompressibleSPHSystem(empty_ic, density_calculator, state_equation, + smoothing_kernel, smoothing_length_; + pressure_acceleration=pressure_acceleration_formulation, + viscosity, density_diffusion, acceleration, + correction, source_terms, + particle_refinement=particle_refinement_) + + # Empty mass vector leads to `nparticles(system_child) = 0` + resize!(system_child.mass, 0) + + particle_refinement.system_child = system_child + + return (system_child, + create_system_child(system_child, system_child.particle_refinement)...) +end + +# ==== Refinement +function refinement!(system::FluidSystem, v, u, v_ode, u_ode, semi) + refinement!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) +end + +refinement!(system, ::Nothing, v, u, v_ode, u_ode, semi) = system + +function refinement!(system, particle_refinement::ParticleRefinement, + v, u, v_ode, u_ode, semi) + check_refinement_criteria!(system, particle_refinement, v, u, v_ode, u_ode, semi) + + refine_particles!(system, particle_refinement, v, u, v_ode, u_ode, semi) +end + +function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, + v, u, v_ode, u_ode, semi) + (; candidates, refinement_criteria) = particle_refinement + + !isempty(candidates) && resize!(candidates, 0) + + for particle in each_moving_particle(system) + for refinement_criterion in refinement_criteria + if particle != last(candidates) && + refinement_criterion(system, particle, v, u, v_ode, u_ode, semi) + push!(candidates, particle) + end + end + end + + return system +end + +function refine_particles!(system_parent, particle_refinement::ParticleRefinement, + v, u, v_ode, u_ode, semi) + (; candidates, system_child) = particle_refinement + + capacity_parent = nparticles(system_parent) - length(candidates) + capacity_child = nparticles(system_child) + + length(candidates) * nchilds(system_parent, particle_refinement) + + if !isempty(candidates) + resize!(system_child, capacity_child) + + particle_refinement.available_childs = capacity_child + + for particle_parent in candidates + particle_childs = child_set(system_child, particle_refinement) + + bear_childs!(system_child, system_parent, particle_childs, particle_parent, + particle_refinement, v_ode, u_ode, semi) + + particle_refinement.available_childs -= nchilds(system, particle_refinement) + end + + resize!(system_parent, capacity_parent) + end +end + +# 6 (8) unkowns in 2d (3D) need to be determined for each newly born child particle +# --> mass, position, velocity, smoothing length +# +# Reducing the dof by using a fixed regular refinement pattern +# (given: position and number of child particles) +function bear_childs!(system_child, system_parent, particle_childs, particle_parent, + particle_refinement, v_ode, u_ode, semi) + # spread child positions according to the refinement pattern + set_positions!(system_child, system_parent, particle_childs, particle_parent, + particle_refinement.refinement_pattern, semi, v_ode, u_ode) + + for particle_child in particle_childs + system_child.mass[particle_child] = particle_refinement.mass_child + + # Interpolate field variables + interpolate_particle_pressure!(system_child, system_parent, + particle_childs, particle_parent, semi, v_ode, u_ode) + interpolate_particle_density!(system_child, system_parent, + particle_childs, particle_parent, semi, v_ode, u_ode) + interpolate_current_velocity!(system_child, system_parent, + particle_childs, particle_parent, semi, v_ode, u_ode) + end + + return system_child +end + +function Base.resize!(system::WeaklyCompressibleSPHSystem, capacity) + (; mass, pressure, cache, density_calculator) = system + + resize!(mass, capacity) + resize!(pressure, capacity) + resize!(cache, capacity, density_calculator) +end + +Base.resize!(cache, capacity, ::SummationDensity) = resize!(cache.density, capacity) +Base.resize!(cache, capacity, ::ContinuityDensity) = cache diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 222b0d906..0e516b4df 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -63,14 +63,9 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear # 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)) + #systems = create_subsystems(systems) + + ranges_u, ranges_v = ranges_uv(systems) # Create (and initialize) a tuple of n neighborhood searches for each of the n systems # We will need one neighborhood search for each pair of systems. @@ -84,6 +79,19 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear return Semidiscretization(systems, ranges_u, ranges_v, searches) end +function ranges_uv(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)) + + return ranges_u, ranges_v +end + # Inline show function e.g. Semidiscretization(neighborhood_search=...) function Base.show(io::IO, semi::Semidiscretization) @nospecialize semi # reduce precompilation time @@ -315,7 +323,7 @@ end @inline function wrap_u(u_ode, system, semi) (; ranges_u) = semi - range = ranges_u[system_indices(system, semi)] + range = ranges_u[system_indices(system, semi)][1] @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system) @@ -329,7 +337,7 @@ end @inline function wrap_v(v_ode, system, semi) (; ranges_v) = semi - range = ranges_v[system_indices(system, semi)] + range = ranges_v[system_indices(system, semi)][1] @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) @@ -430,6 +438,14 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) update_final!(system, v, u, v_ode, u_ode, semi, t) end + + # Particle refinement + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + refinement!(system, v, u, v_ode, u_ode, semi) + end end function update_nhs(u_ode, semi) diff --git a/src/general/system.jl b/src/general/system.jl index c201bf3a5..69934765a 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -19,6 +19,8 @@ initialize!(system, neighborhood_search) = system @inline eachparticle(system) = Base.OneTo(nparticles(system)) @inline each_moving_particle(system) = Base.OneTo(n_moving_particles(system)) +@inline refinement!(system, v, u, v_ode, u_ode, semi) = system + # This should not be dispatched by system type. We always expect to get a column of `A`. @inline function extract_svector(A, system, i) extract_svector(A, Val(ndims(system)), i) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 8e98e3a13..bec18d05e 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -40,7 +40,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, - V, DD, COR, PF, ST, C} <: FluidSystem{NDIMS} + V, DD, COR, PF, ST, PR, C} <: FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] pressure :: Array{ELTYPE, 1} # [particle] @@ -54,6 +54,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, correction :: COR pressure_acceleration_formulation :: PF source_terms :: ST + particle_refinement :: PR cache :: C function WeaklyCompressibleSPHSystem(initial_condition, @@ -63,6 +64,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), + particle_refinement=nothing, correction=nothing, source_terms=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -100,13 +102,12 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, typeof(state_equation), typeof(smoothing_kernel), typeof(viscosity), typeof(density_diffusion), typeof(correction), typeof(pressure_acceleration), - typeof(source_terms), typeof(cache)}(initial_condition, mass, pressure, - density_calculator, state_equation, - smoothing_kernel, smoothing_length, - acceleration_, viscosity, - density_diffusion, correction, - pressure_acceleration, - source_terms, cache) + typeof(source_terms), typeof(particle_refinement), + typeof(cache)}(initial_condition, mass, pressure, density_calculator, + state_equation, smoothing_kernel, smoothing_length, + acceleration_, viscosity, density_diffusion, correction, + pressure_acceleration, source_terms, particle_refinement, + cache) end end From 51f89c2a372c24cc732880d63a2da5c41f267a56 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 14 Mar 2024 19:57:24 +0100 Subject: [PATCH 02/20] revise --- src/TrixiParticles.jl | 1 + src/general/density_calculators.jl | 15 ++++- src/general/refinement.jl | 97 ++++++++++++++++++++---------- src/general/refinement_pattern.jl | 58 ++++++++++++++++++ src/general/semidiscretization.jl | 7 --- 5 files changed, 138 insertions(+), 40 deletions(-) create mode 100644 src/general/refinement_pattern.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index e576cc37b..5ea61a880 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -38,6 +38,7 @@ include("schemes/schemes.jl") include("general/semidiscretization.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("general/refinement.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index ee89c92a6..fa0ba485b 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -35,6 +35,16 @@ end return v[end, particle] end +# *Note* that these functions are intended to internally set the density for buffer particles +# and density correction. It cannot be used to set up an initial condition, +# as the particle density depends on the particle positions. + +@inline set_particle_density(particle, v, ::SummationDensity, system, density) = system + +@inline function set_particle_density(particle, v, ::ContinuityDensity, system, density) + v[end, particle] = density +end + function summation_density!(system, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) @@ -52,7 +62,10 @@ function summation_density!(system, semi, u, u_ode, density; for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, nhs, particles=particles) do particle, neighbor, pos_diff, distance mass = hydrodynamic_mass(neighbor_system, neighbor) - density[particle] += mass * smoothing_kernel(system, distance) + + h_mean = (system.smoothing_length + neighbor_system.smoothing_length)/2 + #density[particle] += mass * smoothing_kernel(system, distance) + density[particle] += mass * kernel(system.smoothing_kernel, distance, h_mean) end end end diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 562af9d3e..2fed9ebd9 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -22,6 +22,10 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} # and will be created in `create_system_child()` system_child::System + # internal `resize!`able storage + _u :: Vector{ELTYPE} + _v :: Vector{ELTYPE} + # API --> parent system with `RL=0` function ParticleRefinement(refinement_criteria...; refinement_pattern=CubicSplitting(), @@ -48,14 +52,8 @@ end @inline Base.ndims(::ParticleRefinement{RL, NDIMS}) where {RL, NDIMS} = ndims -@inline function child_set(system, particle_refinement) - (; available_childs) = particle_refinement - - start_iter = nparticles(system) + 1 - available_childs - end_iter = start_iter + nchilds(system, particle_refinement) - - return start_iter:end_iter -end +@inline child_set(system, particle_refinement) = Base.OneTo(nchilds(system, + particle_refinement)) # ==== Create child systems @@ -151,28 +149,42 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem end function refine_particles!(system_parent, particle_refinement::ParticleRefinement, - v, u, v_ode, u_ode, semi) - (; candidates, system_child) = particle_refinement - - capacity_parent = nparticles(system_parent) - length(candidates) - capacity_child = nparticles(system_child) + - length(candidates) * nchilds(system_parent, particle_refinement) + v_parent, u_parent, v_ode, u_ode, semi) + (; candidates, system_child, _u, _v) = particle_refinement if !isempty(candidates) + n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) + capacity_parent = nparticles(system_parent) - length(candidates) + capacity_child = nparticles(system_child) + n_new_child + + particle_refinement.available_childs = n_new_child + + # Resize child system (extending) resize!(system_child, capacity_child) - particle_refinement.available_childs = capacity_child + # Resize internal storage for new child particles + resize!(_u, n_new_child) + resize!(_v, n_new_child) - for particle_parent in candidates - particle_childs = child_set(system_child, particle_refinement) + u_child = PtrArray(pointer(_u), + (StaticInt(u_nvariables(system_child)), n_new_child)) + v_child = PtrArray(pointer(_v), + (StaticInt(v_nvariables(system_child)), n_new_child)) - bear_childs!(system_child, system_parent, particle_childs, particle_parent, - particle_refinement, v_ode, u_ode, semi) + # Loop over all refinement candidates + for particle_parent in candidates + bear_childs!(system_child, system_parent, particle_parent, particle_refinement, + v_parent, u_parent, v_child, u_child, semi) particle_refinement.available_childs -= nchilds(system, particle_refinement) end - resize!(system_parent, capacity_parent) + # Resize parent system (reducing) + restructure_and_resize!(system_parent, capacity_parent) + + # Resize `v_ode` and `u_ode` + resize!(v_ode, u_ode, semi, system_parent, system_child, + capacity_parent, capacity_child) end end @@ -181,22 +193,37 @@ end # # Reducing the dof by using a fixed regular refinement pattern # (given: position and number of child particles) -function bear_childs!(system_child, system_parent, particle_childs, particle_parent, - particle_refinement, v_ode, u_ode, semi) - # spread child positions according to the refinement pattern - set_positions!(system_child, system_parent, particle_childs, particle_parent, - particle_refinement.refinement_pattern, semi, v_ode, u_ode) +function bear_childs!(system_child, system_parent, particle_parent, particle_refinement, + v_parent, u_parent, v_child, u_child, semi) + (; rel_position_childs, available_childs, candidates) = particle_refinement + n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) + + parent_coords = current_coords(u_parent, system_parent, particle_parent) + + # Loop over all child particles of parent particle + # The number of child particles depends on the refinement pattern + for particle_child in child_set(system_parent, particle_refinement) + absolute_index = particle_child + nparticles(system_child) - available_childs + relative_index = particle_child + n_new_child - available_childs - for particle_child in particle_childs - system_child.mass[particle_child] = particle_refinement.mass_child + system_child.mass[absolute_index] = particle_refinement.mass_child + + # spread child positions according to the refinement pattern + child_coords = parent_coords + rel_position_childs[particle_child] + for dim in 1:ndims(system_child) + u_child[relative_index, dim] = child_coords[dim] + end # Interpolate field variables interpolate_particle_pressure!(system_child, system_parent, - particle_childs, particle_parent, semi, v_ode, u_ode) + particle_childs, particle_parent, + v_parent, u_parent, v_child, u_child, semi) interpolate_particle_density!(system_child, system_parent, - particle_childs, particle_parent, semi, v_ode, u_ode) + particle_childs, particle_parent, + v_parent, u_parent, v_child, u_child, semi) interpolate_current_velocity!(system_child, system_parent, - particle_childs, particle_parent, semi, v_ode, u_ode) + particle_childs, particle_parent, + v_parent, u_parent, v_child, u_child, semi) end return system_child @@ -210,5 +237,11 @@ function Base.resize!(system::WeaklyCompressibleSPHSystem, capacity) resize!(cache, capacity, density_calculator) end -Base.resize!(cache, capacity, ::SummationDensity) = resize!(cache.density, capacity) -Base.resize!(cache, capacity, ::ContinuityDensity) = cache +resize!(cache, capacity, ::SummationDensity) = resize!(cache.density, capacity) +resize!(cache, capacity, ::ContinuityDensity) = cache + +function resize!(v_ode, u_ode, semi, system_parent, system_child, + capacity_parent, capacity_child) + (; particle_refinement) = system_parent + (; candidates) = particle_refinement +end diff --git a/src/general/refinement_pattern.jl b/src/general/refinement_pattern.jl new file mode 100644 index 000000000..12d395c0e --- /dev/null +++ b/src/general/refinement_pattern.jl @@ -0,0 +1,58 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + + function CubicSplitting(; epsilon=0.5, alpha=0.5) + new{typeof(epsilon)}(epsilon, alpha) + end +end + +function relative_positions(refinement_pattern::CubicSplitting, ::System{2}, + particle_spacing) + (; epsilon) = refinement_pattern + + direction_1 = normalize([1.0, 1.0]) + direction_2 = normalize([1.0, -1.0]) + direction_3 = -direction_1 + direction_4 = -direction_2 + + relative_position = hcat(particle_spacing * epsilon * direction_1, + particle_spacing * epsilon * direction_2, + particle_spacing * epsilon * direction_3, + particle_spacing * epsilon * direction_4) + + return reinterpret(reshape, SVector{2, typeof(epsilon)}, relative_position) +end + +# TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? +# See https://www.sciencedirect.com/science/article/pii/S0020740319317023 +@inline nchilds(system, refinement_pattern) = 2^ndims(system) + +@inline smoothing_length_child(system, refinement_pattern) = system.smoothing_length + +@inline mass_child(system, refinement_pattern) = system.mass + +@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing + +# ==== Refinement criteria +struct RefinementZone end # TODO + +@inline function (refinement_criterion::RefinementZone)(system, particle, + v, u, v_ode, u_ode, semi) + (; zone_origin, spanning_set) = refinement_criterion + particle_position = current_coords(u, system, particle) - zone_origin + + for dim in 1:ndims(system) + span_dim = spanning_set[dim] + # Checks whether the projection of the particle position + # falls within the range of the zone. + if !(0 <= dot(particle_position, span_dim) <= dot(span_dim, span_dim)) + + # Particle is not in refinement zone. + return false + end + end + + # Particle is in refinement zone. + return true +end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 0e516b4df..1da9c23ab 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -439,13 +439,6 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) update_final!(system, v, u, v_ode, u_ode, semi, t) end - # Particle refinement - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - - refinement!(system, v, u, v_ode, u_ode, semi) - end end function update_nhs(u_ode, semi) From d2ebf190dcaed1ab10b11dd53be75a9e89e019e4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 16 Mar 2024 12:08:35 +0100 Subject: [PATCH 03/20] revise resize --- src/callbacks/refinement_callback.jl | 131 ++++++++++++++++++ src/general/refinement.jl | 191 ++++++++++++++++++++++++--- 2 files changed, 306 insertions(+), 16 deletions(-) create mode 100644 src/callbacks/refinement_callback.jl diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl new file mode 100644 index 000000000..9783f4365 --- /dev/null +++ b/src/callbacks/refinement_callback.jl @@ -0,0 +1,131 @@ +mutable struct ParticleRefinementCallback{I} + interval :: I + n_candidates :: Vector{Int} + n_childs :: Vector{Int} + ranges_u_cache :: RU + ranges_v_cache :: RV + eachparticle_cache :: RP + + # internal `resize!`able storage + _u_ode::Vector{Float64} + _v_ode::Vector{Float64} +end + +function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) + if dt > 0 && interval !== -1 + throw(ArgumentError("Setting both interval and dt is not supported!")) + end + + # Update in intervals in terms of simulation time + if dt > 0 + interval = Float64(dt) + + # Update every time step (default) + elseif interval == -1 + interval = 1 + end + + update_callback! = ParticleRefinementCallback(interval) + + if dt > 0 && update + # Add a `tstop` every `dt`, and save the final solution. + return PeriodicCallback(update_callback!, dt, + initialize=initial_update!, + save_positions=(false, false)) + else + # The first one is the condition, the second the affect! + return DiscreteCallback(update_callback!, update_callback!, + initialize=initial_update!, + save_positions=(false, false)) + end +end + +# initialize +function initial_update!(cb, u, t, integrator) + # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`) + # or `cb.affect!.affect!` (with `PeriodicCallback`). + # Let recursive dispatch handle this. + + initial_update!(cb.affect!, u, t, integrator) +end + +function initial_update!(cb::ParticleRefinementCallback, u, t, integrator) + cb.update && cb(integrator) +end + +# condition +function (update_callback!::ParticleRefinementCallback)(u, t, integrator) + (; interval, update) = update_callback! + + # With error-based step size control, some steps can be rejected. Thus, + # `integrator.iter >= integrator.stats.naccept` + # (total #steps) (#accepted steps) + # We need to check the number of accepted steps since callbacks are not + # activated after a rejected step. + return update && (integrator.stats.naccept % interval == 0) +end + +# affect +function (update_callback!::ParticleRefinementCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # 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. + refinement!(v_ode, u_ode, semi, callback) + + update_systems_and_nhs(v_ode, u_ode, semi, t) + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(interval=", + (cb.affect!.update ? cb.affect!.interval : "-"), + ")") +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(dt=", cb.affect!.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect! + setup = [ + "interval" => update_cb.update ? update_cb.interval : "-", + "update" => update_cb.update ? "yes" : "no" + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect!.affect! + setup = [ + "dt" => update_cb.interval, + "update" => update_cb.update ? "yes" : "no" + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 2fed9ebd9..09d810402 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -22,10 +22,6 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} # and will be created in `create_system_child()` system_child::System - # internal `resize!`able storage - _u :: Vector{ELTYPE} - _v :: Vector{ELTYPE} - # API --> parent system with `RL=0` function ParticleRefinement(refinement_criteria...; refinement_pattern=CubicSplitting(), @@ -55,6 +51,9 @@ end @inline child_set(system, particle_refinement) = Base.OneTo(nchilds(system, particle_refinement)) +@inline ncandidates(::Nothing) = 0 +@inline ncandidates(particle_refinement) = particle_refinement.candidates + # ==== Create child systems function create_system_childs(systems) @@ -117,23 +116,29 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, end # ==== Refinement -function refinement!(system::FluidSystem, v, u, v_ode, u_ode, semi) - refinement!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) -end +function refinement!(v_ode, u_ode, semi, callback) + foreach_system(semi) do system + check_refinement_criteria!(system, v_ode, u_ode, semi) + end -refinement!(system, ::Nothing, v, u, v_ode, u_ode, semi) = system + resize_and_copy!(callback, semi, v_ode, u_ode) -function refinement!(system, particle_refinement::ParticleRefinement, - v, u, v_ode, u_ode, semi) - check_refinement_criteria!(system, particle_refinement, v, u, v_ode, u_ode, semi) + refine_particles!(callback, semi, v_ode, u_ode) +end - refine_particles!(system, particle_refinement, v, u, v_ode, u_ode, semi) +check_refinement_criteria!(system, v_ode, u_ode, semi) = system + +function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) + check_refinement_criteria!(system, system.particle_refinement, v_ode, u_ode, semi) end function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, - v, u, v_ode, u_ode, semi) + v_ode, u_ode, semi) (; candidates, refinement_criteria) = particle_refinement + v = wrap_v(v_ode, systm, semi) + u = wrap_u(u_ode, systm, semi) + !isempty(candidates) && resize!(candidates, 0) for particle in each_moving_particle(system) @@ -148,6 +153,56 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem return system end +function resize_and_copy!(callback, semi, v_ode, u_ode) + (; _v_ode, _u_ode, n_candidates, n_childs, + ranges_v_cache, ranges_u_cache, each_particle_cache) = callback + + # Get non-`resize!`d ranges + ranges_v_cache, ranges_u_cache = ranges_uv(systems) + each_particle_cache = Tuple(eachparticle(system) for system in semi.systems) + + # Resize internal storage + n_total_particles = sum(nparticles.(semi.systems)) + resize!(_v_ode, capacity, n_total_particles) + resize!(_u_ode, capacity, n_total_particles) + + # Count all candidates for refinement + n_candidates = 0 + n_childs = 0 + foreach_system(semi) do system + n_candidates += ncandidates(system.particle_refinement) + n_childs += ncandidates(system.particle_refinement) * + nchilds(system, system.refinement_pattern) + end + + capacity = n_total_particles - n_candidates + n_childs + + # Resize integrated values + resize!(v_ode, capacity) + resize!(u_ode, capacity) + + # Resize all systems + foreach_system(semi) do system + resize!(system, system.particle_refinement) + end + + # Set `resize!`d ranges + semi.ranges_v, semi.ranges_u = ranges_uv(systems) + + # Preserve non-changing values + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + _v = _wrap_v(_v_ode, system, callback) + _u = _wrap_u(_u_ode, system, callback) + + copy_values_v!(v, _v, system, system.particle_refinement, semi) + copy_values_u!(u, _u, system, system.particle_refinement, semi) + end + + return callback +end + function refine_particles!(system_parent, particle_refinement::ParticleRefinement, v_parent, u_parent, v_ode, u_ode, semi) (; candidates, system_child, _u, _v) = particle_refinement @@ -229,6 +284,24 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref return system_child end +@inline Base.resize!(system::System, ::Nothing) = system + +@inline function Base.resize!(system::System, particle_refinement::ParticleRefinement) + (; candidates, system_child) = particle_refinement + + if !isempty(candidates) + n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) + capacity_parent = nparticles(system_parent) - length(candidates) + capacity_child = nparticles(system_child) + n_new_child + + # Resize child system (extending) + resize!(system_child, capacity_child) + + # Resize parent system (reducing) + resize!(system_parent, capacity_parent) + end +end + function Base.resize!(system::WeaklyCompressibleSPHSystem, capacity) (; mass, pressure, cache, density_calculator) = system @@ -240,8 +313,94 @@ end resize!(cache, capacity, ::SummationDensity) = resize!(cache.density, capacity) resize!(cache, capacity, ::ContinuityDensity) = cache -function resize!(v_ode, u_ode, semi, system_parent, system_child, - capacity_parent, capacity_child) - (; particle_refinement) = system_parent +@inline function _wrap_u(_u_ode, system, callback) + (; ranges_u_cache) = callback + + range = ranges_u_cache[system_indices(system, semi)] + + @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system) + + # This is a non-allocating version of: + # return unsafe_wrap(Array{eltype(_u_ode), 2}, pointer(view(_u_ode, range)), + # (u_nvariables(system), n_moving_particles(system))) + return PtrArray(pointer(view(_u_ode, range)), + (StaticInt(u_nvariables(system)), n_moving_particles(system))) +end + +@inline function _wrap_v(_v_ode, system, callback) + (; ranges_v_cache) = callback + + range = ranges_v_cache[system_indices(system, semi)][1] + + @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) + + # This is a non-allocating version of: + # return unsafe_wrap(Array{eltype(_v_ode), 2}, pointer(view(_v_ode, range)), + # (v_nvariables(system), n_moving_particles(system))) + return PtrArray(pointer(view(_v_ode, range)), + (StaticInt(v_nvariables(system)), n_moving_particles(system))) +end + +# `v_new` >= `v_old` +function copy_values_v!(v_new, v_old, system, ::Nothing, semi, callback) + (; each_particle_cache) = callback + + for particle in each_particle_cache[system_indices(system, semi)] + for i in 1:v_nvariables(system) + v_new[i, particle] = v_old[i, particle] + end + end +end + +# `v_new` <= `v_old` +function copy_values_v!(v_new, v_old, system, particle_refinement::ParticleRefinement, + semi, callback) (; candidates) = particle_refinement + (; eachparticle_cache) = callback + + # Mark candidates + v_old[1, candidates] .= Inf + + # Copy only non-refined particles + new_particle_id = 1 + for particle in eachparticle_cache[system_indices(system, semi)] + if v_new[1, particle] < Inf + for i in 1:v_nvariables(system) + v_new[i, new_particle_id] = v_old[i, particle] + end + new_particle_id += 1 + end + end +end + +# `u_new` >= `u_old` +function copy_values_u!(u_new, u_old, system, ::Nothing, semi, callback) + (; each_particle_cache) = callback + + for particle in each_particle_cache[system_indices(system, semi)] + for i in 1:u_nuariables(system) + u_new[i, particle] = u_old[i, particle] + end + end +end + +# `u_new` <= `u_old` +function copy_values_u!(u_new, u_old, system, particle_refinement::ParticleRefinement, + semi, callback) + (; candidates) = particle_refinement + (; eachparticle_cache) = callback + + # Mark candidates + u_old[1, candidates] .= Inf + + # Copy only non-refined particles + new_particle_id = 1 + for particle in eachparticle_cache[system_indices(system, semi)] + if u_new[1, particle] < Inf + for i in 1:u_nuariables(system) + u_new[i, new_particle_id] = u_old[i, particle] + end + new_particle_id += 1 + end + end end From 9b1526865082807e05dd7500160106442fd494b7 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Sat, 16 Mar 2024 12:09:04 +0100 Subject: [PATCH 04/20] revise refinement --- src/callbacks/refinement_callback.jl | 56 +++++----- src/general/refinement.jl | 153 ++++++++++++++++----------- 2 files changed, 114 insertions(+), 95 deletions(-) diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl index 9783f4365..a4bd3ca0c 100644 --- a/src/callbacks/refinement_callback.jl +++ b/src/callbacks/refinement_callback.jl @@ -1,10 +1,10 @@ mutable struct ParticleRefinementCallback{I} interval :: I - n_candidates :: Vector{Int} - n_childs :: Vector{Int} - ranges_u_cache :: RU - ranges_v_cache :: RV - eachparticle_cache :: RP + n_candidates :: Int + n_childs :: Int + ranges_u_cache :: NTuple + ranges_v_cache :: NTuple + eachparticle_cache :: NTuple # internal `resize!`able storage _u_ode::Vector{Float64} @@ -25,55 +25,53 @@ function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) interval = 1 end - update_callback! = ParticleRefinementCallback(interval) + refinement_callback = ParticleRefinementCallback(interval, 0, 0, (), (), (), [0.0], [0.0]) - if dt > 0 && update + if dt > 0 # Add a `tstop` every `dt`, and save the final solution. - return PeriodicCallback(update_callback!, dt, - initialize=initial_update!, + return PeriodicCallback(refinement_callback, dt, + initialize=initial_refinement!, save_positions=(false, false)) else # The first one is the condition, the second the affect! - return DiscreteCallback(update_callback!, update_callback!, - initialize=initial_update!, + return DiscreteCallback(refinement_callback, refinement_callback, + initialize=initial_refinement!, save_positions=(false, false)) end end # initialize -function initial_update!(cb, u, t, integrator) +function initial_refinement!(cb, u, t, integrator) # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`) # or `cb.affect!.affect!` (with `PeriodicCallback`). # Let recursive dispatch handle this. - initial_update!(cb.affect!, u, t, integrator) + initial_refinement!(cb.affect!, u, t, integrator) end -function initial_update!(cb::ParticleRefinementCallback, u, t, integrator) - cb.update && cb(integrator) +function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator) + cb(integrator) end # condition -function (update_callback!::ParticleRefinementCallback)(u, t, integrator) - (; interval, update) = update_callback! +function (refinement_callback::ParticleRefinementCallback)(u, t, integrator) + (; interval) = refinement_callback # With error-based step size control, some steps can be rejected. Thus, # `integrator.iter >= integrator.stats.naccept` # (total #steps) (#accepted steps) # We need to check the number of accepted steps since callbacks are not # activated after a rejected step. - return update && (integrator.stats.naccept % interval == 0) + return integrator.stats.naccept % interval == 0 end # affect -function (update_callback!::ParticleRefinementCallback)(integrator) +function (refinement_callback::ParticleRefinementCallback)(integrator) t = integrator.t semi = integrator.p v_ode, u_ode = integrator.u.x - # 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. - refinement!(v_ode, u_ode, semi, callback) + refinement!(v_ode, u_ode, semi, refinement_callback) update_systems_and_nhs(v_ode, u_ode, semi, t) @@ -85,9 +83,7 @@ end function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) @nospecialize cb # reduce precompilation time - print(io, "ParticleRefinementCallback(interval=", - (cb.affect!.update ? cb.affect!.interval : "-"), - ")") + print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")") end function Base.show(io::IO, @@ -104,10 +100,9 @@ function Base.show(io::IO, ::MIME"text/plain", if get(io, :compact, false) show(io, cb) else - update_cb = cb.affect! + refinement_cb = cb.affect! setup = [ - "interval" => update_cb.update ? update_cb.interval : "-", - "update" => update_cb.update ? "yes" : "no" + "interval" => refinement_cb.interval, ] summary_box(io, "ParticleRefinementCallback", setup) end @@ -121,10 +116,9 @@ function Base.show(io::IO, ::MIME"text/plain", if get(io, :compact, false) show(io, cb) else - update_cb = cb.affect!.affect! + refinement_cb = cb.affect!.affect! setup = [ - "dt" => update_cb.interval, - "update" => update_cb.update ? "yes" : "no" + "dt" => refinement_cb.interval, ] summary_box(io, "ParticleRefinementCallback", setup) end diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 09d810402..e4291cb38 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -159,7 +159,7 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) # Get non-`resize!`d ranges ranges_v_cache, ranges_u_cache = ranges_uv(systems) - each_particle_cache = Tuple(eachparticle(system) for system in semi.systems) + each_particle_cache = Tuple(get_iterator(system) for system in semi.systems) # Resize internal storage n_total_particles = sum(nparticles.(semi.systems)) @@ -172,7 +172,7 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) foreach_system(semi) do system n_candidates += ncandidates(system.particle_refinement) n_childs += ncandidates(system.particle_refinement) * - nchilds(system, system.refinement_pattern) + nchilds(system, system.refinement_pattern) end capacity = n_total_particles - n_candidates + n_childs @@ -196,50 +196,48 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) _v = _wrap_v(_v_ode, system, callback) _u = _wrap_u(_u_ode, system, callback) - copy_values_v!(v, _v, system, system.particle_refinement, semi) - copy_values_u!(u, _u, system, system.particle_refinement, semi) + copy_values_v!(v, _v, system, system.particle_refinement, semi, callback) + copy_values_u!(u, _u, system, system.particle_refinement, semi, callback) end return callback end -function refine_particles!(system_parent, particle_refinement::ParticleRefinement, - v_parent, u_parent, v_ode, u_ode, semi) - (; candidates, system_child, _u, _v) = particle_refinement +refine_particles!(semi::Semidiscretization{Nothing}, v_ode, u_ode) = semi - if !isempty(candidates) - n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) - capacity_parent = nparticles(system_parent) - length(candidates) - capacity_child = nparticles(system_child) + n_new_child +function refine_particles!(callback, semi, v_ode, u_ode) - particle_refinement.available_childs = n_new_child + # Refine particles in all systems + foreach_system(semi) do system + refine_particles!(system, sysem.particle_refinement, v_ode, u_ode, callback, semi) + end +end - # Resize child system (extending) - resize!(system_child, capacity_child) +refine_particles!(system, ::Nothing, v, u, semi) = system - # Resize internal storage for new child particles - resize!(_u, n_new_child) - resize!(_v, n_new_child) +function refine_particles!(system_parent, particle_refinement::ParticleRefinement, + v_ode, u_ode, callback, semi) + (; _v_ode, _u_ode) = callback + (; candidates, system_child, available_childs) = particle_refinement - u_child = PtrArray(pointer(_u), - (StaticInt(u_nvariables(system_child)), n_new_child)) - v_child = PtrArray(pointer(_v), - (StaticInt(v_nvariables(system_child)), n_new_child)) + if !isempty(candidates) + # Old storage + v_parent = _wrap_v(_v_ode, system_parent, callback) + u_parent = _wrap_u(_u_ode, system_parent, callback) + + # Resized storage + v_child = wrap_v(v_ode, system_child, semi) + u_child = wrap_u(u_ode, system_child, semi) + + available_childs = length(candidates) * nchilds(system_parent, particle_refinement) # Loop over all refinement candidates for particle_parent in candidates bear_childs!(system_child, system_parent, particle_parent, particle_refinement, v_parent, u_parent, v_child, u_child, semi) - particle_refinement.available_childs -= nchilds(system, particle_refinement) + available_childs -= nchilds(system, particle_refinement) end - - # Resize parent system (reducing) - restructure_and_resize!(system_parent, capacity_parent) - - # Resize `v_ode` and `u_ode` - resize!(v_ode, u_ode, semi, system_parent, system_child, - capacity_parent, capacity_child) end end @@ -250,35 +248,62 @@ end # (given: position and number of child particles) function bear_childs!(system_child, system_parent, particle_parent, particle_refinement, v_parent, u_parent, v_child, u_child, semi) - (; rel_position_childs, available_childs, candidates) = particle_refinement - n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) + (; rel_position_childs, available_childs, mass_child) = particle_refinement + nhs = get_neighborhood_search(system_parent, system_parent, semi) parent_coords = current_coords(u_parent, system_parent, particle_parent) + system_child.mass .= mass_child + # Loop over all child particles of parent particle # The number of child particles depends on the refinement pattern for particle_child in child_set(system_parent, particle_refinement) absolute_index = particle_child + nparticles(system_child) - available_childs - relative_index = particle_child + n_new_child - available_childs - - system_child.mass[absolute_index] = particle_refinement.mass_child # spread child positions according to the refinement pattern child_coords = parent_coords + rel_position_childs[particle_child] for dim in 1:ndims(system_child) - u_child[relative_index, dim] = child_coords[dim] + u_child[absolute_index, dim] = child_coords[dim] + end + + volume = zero(eltype(system_child)) + p_a = zero(eltype(system_child)) + rho_a = zero(eltype(system_child)) + + for neighbor in eachneighbor(child_coords, nhs) + neighbor_coords = current_coords(u_parent, system_parent, neighbor) + pos_diff = child_coords - neighbor_coords + + distance2 = dot(pos_diff, pos_diff) + + if distance2 <= search_radius^2 + distance = sqrt(distance2) + kernel_weight = smoothing_kernel(system_parent, distance) + volume += kernel_weight + + v_b = current_velocity(v_parent, system_parent, neighbor) + p_b = particle_pressure(v_parent, system_parent, neighbor) + rho_b = particle_density(v_parent, system_parent, neighbor) + + for dim in 1:ndims(system_child) + v_child[dim, system_child] += kernel_weight * v_b[dim] + end + + rho_a += kernel_weight * rho_b + p_a += kernel_weight * p_b + end + end + + for dim in 1:ndims(system_child) + v[dim, system_child] ./ volume end - # Interpolate field variables - interpolate_particle_pressure!(system_child, system_parent, - particle_childs, particle_parent, - v_parent, u_parent, v_child, u_child, semi) - interpolate_particle_density!(system_child, system_parent, - particle_childs, particle_parent, - v_parent, u_parent, v_child, u_child, semi) - interpolate_current_velocity!(system_child, system_parent, - particle_childs, particle_parent, - v_parent, u_parent, v_child, u_child, semi) + rho_a /= volume + p_a /= volume + + set_particle_density(particle_child, v_child, system_child.density_calculator, + system_child, rho_a) + set_particle_pressure(particle_child, v_child, system_child, p_a) end return system_child @@ -330,7 +355,7 @@ end @inline function _wrap_v(_v_ode, system, callback) (; ranges_v_cache) = callback - range = ranges_v_cache[system_indices(system, semi)][1] + range = ranges_v_cache[system_indices(system, semi)] @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) @@ -355,21 +380,15 @@ end # `v_new` <= `v_old` function copy_values_v!(v_new, v_old, system, particle_refinement::ParticleRefinement, semi, callback) - (; candidates) = particle_refinement (; eachparticle_cache) = callback - # Mark candidates - v_old[1, candidates] .= Inf - # Copy only non-refined particles new_particle_id = 1 for particle in eachparticle_cache[system_indices(system, semi)] - if v_new[1, particle] < Inf - for i in 1:v_nvariables(system) - v_new[i, new_particle_id] = v_old[i, particle] - end - new_particle_id += 1 + for i in 1:v_nvariables(system) + v_new[i, new_particle_id] = v_old[i, particle] end + new_particle_id += 1 end end @@ -387,20 +406,26 @@ end # `u_new` <= `u_old` function copy_values_u!(u_new, u_old, system, particle_refinement::ParticleRefinement, semi, callback) - (; candidates) = particle_refinement (; eachparticle_cache) = callback - # Mark candidates - u_old[1, candidates] .= Inf - # Copy only non-refined particles new_particle_id = 1 for particle in eachparticle_cache[system_indices(system, semi)] - if u_new[1, particle] < Inf - for i in 1:u_nuariables(system) - u_new[i, new_particle_id] = u_old[i, particle] - end - new_particle_id += 1 + for i in 1:u_nuariables(system) + u_new[i, new_particle_id] = u_old[i, particle] end + new_particle_id += 1 end end + +@inline get_iterator(system) = get_iterator(system, system.particle_refinement) + +@inline get_iterator(system, ::Nothing) = eachparticle(system) + +@inline function get_iterator(system, particle_refinement::ParticleRefinement) + (; candidates) = particle_refinement + + # Filter candidates + #return Iterators.filter(i -> !(i in candidates), eachparticle(system)) + return setdiff(eachparticle(system), candidates) +end From cb28f4ef85c99e484bd85db8cc4e7a735ab24183 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Mar 2024 18:02:07 +0100 Subject: [PATCH 05/20] fix callback --- src/callbacks/callbacks.jl | 1 + src/callbacks/refinement_callback.jl | 15 +++++++-------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 08b28187a..fa121d1bc 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("solution_saving.jl") include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") +include("refinement_callback.jl") diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl index a4bd3ca0c..e342581a1 100644 --- a/src/callbacks/refinement_callback.jl +++ b/src/callbacks/refinement_callback.jl @@ -1,10 +1,9 @@ mutable struct ParticleRefinementCallback{I} interval :: I - n_candidates :: Int - n_childs :: Int - ranges_u_cache :: NTuple - ranges_v_cache :: NTuple - eachparticle_cache :: NTuple + ranges_u_cache :: Tuple + ranges_v_cache :: Tuple + nparticles_cache :: Tuple + eachparticle_cache :: Tuple # internal `resize!`able storage _u_ode::Vector{Float64} @@ -25,7 +24,7 @@ function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) interval = 1 end - refinement_callback = ParticleRefinementCallback(interval, 0, 0, (), (), (), [0.0], [0.0]) + refinement_callback = ParticleRefinementCallback(interval, (), (), (), (), [0.0], [0.0]) if dt > 0 # Add a `tstop` every `dt`, and save the final solution. @@ -102,7 +101,7 @@ function Base.show(io::IO, ::MIME"text/plain", else refinement_cb = cb.affect! setup = [ - "interval" => refinement_cb.interval, + "interval" => refinement_cb.interval ] summary_box(io, "ParticleRefinementCallback", setup) end @@ -118,7 +117,7 @@ function Base.show(io::IO, ::MIME"text/plain", else refinement_cb = cb.affect!.affect! setup = [ - "dt" => refinement_cb.interval, + "dt" => refinement_cb.interval ] summary_box(io, "ParticleRefinementCallback", setup) end From eed7bcc2a8c4bcd2b166e310135ea926a4977439 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Mar 2024 18:02:51 +0100 Subject: [PATCH 06/20] add setter --- src/general/density_calculators.jl | 10 ---------- src/general/refinement.jl | 14 ++++++++++++++ src/schemes/fluid/fluid.jl | 23 +++++++++++++++++++++++ 3 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index fa0ba485b..3b4f257b9 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -35,16 +35,6 @@ end return v[end, particle] end -# *Note* that these functions are intended to internally set the density for buffer particles -# and density correction. It cannot be used to set up an initial condition, -# as the particle density depends on the particle positions. - -@inline set_particle_density(particle, v, ::SummationDensity, system, density) = system - -@inline function set_particle_density(particle, v, ::ContinuityDensity, system, density) - v[end, particle] = density -end - function summation_density!(system, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) diff --git a/src/general/refinement.jl b/src/general/refinement.jl index e4291cb38..32f072bf2 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -429,3 +429,17 @@ end #return Iterators.filter(i -> !(i in candidates), eachparticle(system)) return setdiff(eachparticle(system), candidates) end + +@inline particle_pressure_parent(v, ::WeaklyCompressibleSPHSystem, particle) = 0.0 +@inline particle_pressure_parent(v, system::EntropicallyDampedSPHSystem, particle) = particle_pressure(v, + system, + particle) + +@inline particle_density_parent(v, system, particle) = particle_density_parent(v, system, + system.density_calculator, + particle) + +@inline particle_density_parent(v, system, ::SummationDensity, particle) = 0.0 +@inline particle_density_parent(v, system, ::ContinuityDensity, particle) = particle_density(v, + system, + particle) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 14f16cdf6..514791eec 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -79,3 +79,26 @@ include("pressure_acceleration.jl") include("viscosity.jl") include("weakly_compressible_sph/weakly_compressible_sph.jl") include("entropically_damped_sph/entropically_damped_sph.jl") + +# *Note* that these functions are intended to internally set the density for buffer particles +# and density correction. It cannot be used to set up an initial condition, +# as the particle density depends on the particle positions. + +@inline set_particle_density(particle, v, ::SummationDensity, system, density) = nothing + +@inline function set_particle_density(particle, v, ::ContinuityDensity, + system::WeaklyCompressibleSPHSystem, density) + v[end, particle] = density +end + +@inline function set_particle_density(particle, v, ::ContinuityDensity, + system::EntropicallyDampedSPHSystem, density) + v[end - 1, particle] = density +end + +@inline set_particle_pressure(particle, v, system, pressure) = nothing + +@inline function set_particle_pressure(particle, v, system::EntropicallyDampedSPHSystem, + pressure) + v[end, particle] = pressure +end From bf887ded4abe0ae169d8cb449dcbfa59922c0a2f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Mar 2024 18:07:43 +0100 Subject: [PATCH 07/20] fix resize stuff --- src/general/refinement.jl | 208 ++++++++++++++---------------- src/general/refinement_pattern.jl | 69 +++++++++- src/general/semidiscretization.jl | 8 +- 3 files changed, 160 insertions(+), 125 deletions(-) diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 32f072bf2..83c7c9ed6 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -7,16 +7,15 @@ include("refinement_pattern.jl") mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} - candidates :: Vector{ELTYPE} + candidates :: Vector{Int} refinement_levels :: Int refinement_pattern :: RP refinement_criteria :: RC + available_childs :: Int # Depends on refinement pattern, particle spacing and parameters ϵ and α. # Should be obtained prior to simulation in `create_system_child()` - rel_position_childs :: Vector{SVector{NDIMS, ELTYPE}} - mass_child :: ELTYPE - available_childs :: Int + rel_position_childs::Vector{SVector{NDIMS, ELTYPE}} # It is essential to know the child system, which is empty at the beginning # and will be created in `create_system_child()` @@ -26,23 +25,23 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} function ParticleRefinement(refinement_criteria...; refinement_pattern=CubicSplitting(), refinement_levels=1) - ELTYPE = eltype(refinement_criteria) + ELTYPE = eltype(refinement_criteria[1]) NDIMS = ndims(refinement_criteria[1]) return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], refinement_levels, refinement_pattern, - refinement_criteria) + typeof(refinement_criteria)}([0], refinement_levels, refinement_pattern, + refinement_criteria, 0) end # Internal constructor for multiple refinement levels function ParticleRefinement{RL}(refinement_criteria::Tuple, refinement_pattern, refinement_levels) where {RL} - ELTYPE = eltype(refinement_criteria) + ELTYPE = eltype(refinement_criteria[1]) NDIMS = ndims(refinement_criteria[1]) return new{RL, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], refinement_levels, refinement_pattern, - refinement_criteria) + typeof(refinement_criteria)}([0], refinement_levels, refinement_pattern, + refinement_criteria, 0) end end @@ -51,8 +50,8 @@ end @inline child_set(system, particle_refinement) = Base.OneTo(nchilds(system, particle_refinement)) -@inline ncandidates(::Nothing) = 0 -@inline ncandidates(particle_refinement) = particle_refinement.candidates +@inline nchilds(system, ::Nothing) = 0 +@inline nchilds(system, pr::ParticleRefinement) = nchilds(system, pr.refinement_pattern) # ==== Create child systems @@ -78,6 +77,7 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, # Distribute values according to refinement pattern smoothing_length_ = smoothing_length_child(system, refinement_pattern) + particle_refinement.rel_position_childs = refinement_pattern(system) # Create "empty" `InitialCondition` for child system particle_spacing_ = particle_spacing_child(system, refinement_pattern) @@ -107,7 +107,7 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, particle_refinement=particle_refinement_) # Empty mass vector leads to `nparticles(system_child) = 0` - resize!(system_child.mass, 0) + Base.resize!(system_child.mass, 0) particle_refinement.system_child = system_child @@ -132,18 +132,20 @@ function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) check_refinement_criteria!(system, system.particle_refinement, v_ode, u_ode, semi) end +check_refinement_criteria!(system, ::Nothing, v_ode, u_ode, semi) = system + function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, v_ode, u_ode, semi) (; candidates, refinement_criteria) = particle_refinement - v = wrap_v(v_ode, systm, semi) - u = wrap_u(u_ode, systm, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) - !isempty(candidates) && resize!(candidates, 0) + !isempty(candidates) && Base.resize!(candidates, 0) for particle in each_moving_particle(system) for refinement_criterion in refinement_criteria - if particle != last(candidates) && + if (isempty(candidates) || particle != last(candidates)) && refinement_criterion(system, particle, v, u, v_ode, u_ode, semi) push!(candidates, particle) end @@ -154,89 +156,87 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem end function resize_and_copy!(callback, semi, v_ode, u_ode) - (; _v_ode, _u_ode, n_candidates, n_childs, - ranges_v_cache, ranges_u_cache, each_particle_cache) = callback + (; _v_ode, _u_ode) = callback # Get non-`resize!`d ranges - ranges_v_cache, ranges_u_cache = ranges_uv(systems) - each_particle_cache = Tuple(get_iterator(system) for system in semi.systems) + callback.ranges_v_cache, callback.ranges_u_cache = ranges_vu(semi.systems) + callback.eachparticle_cache = Tuple(get_iterator(system) for system in semi.systems) + callback.nparticles_cache = Tuple(n_moving_particles(system) for system in semi.systems) # Resize internal storage - n_total_particles = sum(nparticles.(semi.systems)) - resize!(_v_ode, capacity, n_total_particles) - resize!(_u_ode, capacity, n_total_particles) - - # Count all candidates for refinement - n_candidates = 0 - n_childs = 0 - foreach_system(semi) do system - n_candidates += ncandidates(system.particle_refinement) - n_childs += ncandidates(system.particle_refinement) * - nchilds(system, system.refinement_pattern) - end + Base.resize!(_v_ode, length(v_ode)) + Base.resize!(_u_ode, length(u_ode)) - capacity = n_total_particles - n_candidates + n_childs - - # Resize integrated values - resize!(v_ode, capacity) - resize!(u_ode, capacity) + _v_ode .= v_ode + _u_ode .= u_ode # Resize all systems foreach_system(semi) do system - resize!(system, system.particle_refinement) + resize_system!(system, system.particle_refinement) end # Set `resize!`d ranges - semi.ranges_v, semi.ranges_u = ranges_uv(systems) + ranges_v_tmp, ranges_u_tmp = ranges_vu(semi.systems) + for i in 1:length(semi.systems) + semi.ranges_v[i][1] = ranges_v_tmp[i][1] + semi.ranges_u[i][1] = ranges_u_tmp[i][1] + end + + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in semi.systems) + sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in semi.systems) + + # Resize integrated values + Base.resize!(v_ode, sum(sizes_v)) + Base.resize!(u_ode, sum(sizes_u)) # Preserve non-changing values foreach_system(semi) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - _v = _wrap_v(_v_ode, system, callback) - _u = _wrap_u(_u_ode, system, callback) + _v = _wrap_v(_v_ode, system, semi, callback) + _u = _wrap_u(_u_ode, system, semi, callback) - copy_values_v!(v, _v, system, system.particle_refinement, semi, callback) - copy_values_u!(u, _u, system, system.particle_refinement, semi, callback) + copy_values_v!(v, _v, system, semi, callback) + copy_values_u!(u, _u, system, semi, callback) end return callback end -refine_particles!(semi::Semidiscretization{Nothing}, v_ode, u_ode) = semi - function refine_particles!(callback, semi, v_ode, u_ode) # Refine particles in all systems foreach_system(semi) do system - refine_particles!(system, sysem.particle_refinement, v_ode, u_ode, callback, semi) + refine_particles!(system, system.particle_refinement, v_ode, u_ode, callback, semi) end end -refine_particles!(system, ::Nothing, v, u, semi) = system +refine_particles!(system, ::Nothing, v_ode, u_ode, callback, semi) = system function refine_particles!(system_parent, particle_refinement::ParticleRefinement, v_ode, u_ode, callback, semi) (; _v_ode, _u_ode) = callback - (; candidates, system_child, available_childs) = particle_refinement + (; candidates, system_child) = particle_refinement if !isempty(candidates) # Old storage - v_parent = _wrap_v(_v_ode, system_parent, callback) - u_parent = _wrap_u(_u_ode, system_parent, callback) + v_parent = _wrap_v(_v_ode, system_parent, semi, callback) + u_parent = _wrap_u(_u_ode, system_parent, semi, callback) # Resized storage v_child = wrap_v(v_ode, system_child, semi) u_child = wrap_u(u_ode, system_child, semi) - available_childs = length(candidates) * nchilds(system_parent, particle_refinement) + particle_refinement.available_childs = length(candidates) * + nchilds(system_parent, particle_refinement) # Loop over all refinement candidates for particle_parent in candidates bear_childs!(system_child, system_parent, particle_parent, particle_refinement, v_parent, u_parent, v_child, u_child, semi) - available_childs -= nchilds(system, particle_refinement) + particle_refinement.available_childs -= nchilds(system_parent, + particle_refinement) end end end @@ -248,22 +248,26 @@ end # (given: position and number of child particles) function bear_childs!(system_child, system_parent, particle_parent, particle_refinement, v_parent, u_parent, v_child, u_child, semi) - (; rel_position_childs, available_childs, mass_child) = particle_refinement + (; rel_position_childs, available_childs, refinement_pattern) = particle_refinement nhs = get_neighborhood_search(system_parent, system_parent, semi) parent_coords = current_coords(u_parent, system_parent, particle_parent) - system_child.mass .= mass_child - # Loop over all child particles of parent particle # The number of child particles depends on the refinement pattern for particle_child in child_set(system_parent, particle_refinement) absolute_index = particle_child + nparticles(system_child) - available_childs + # TODO: Handle different masses. Problem: `particle_parent` does not have an + # mass-entry anymore, since `system_parent` was resized. + mass_ = system_parent.mass[1] + system_child.mass[absolute_index] = mass_child(system_parent, mass_, + refinement_pattern) + # spread child positions according to the refinement pattern child_coords = parent_coords + rel_position_childs[particle_child] for dim in 1:ndims(system_child) - u_child[absolute_index, dim] = child_coords[dim] + u_child[dim, absolute_index] = child_coords[dim] end volume = zero(eltype(system_child)) @@ -276,17 +280,17 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref distance2 = dot(pos_diff, pos_diff) - if distance2 <= search_radius^2 + if distance2 <= nhs.search_radius^2 distance = sqrt(distance2) kernel_weight = smoothing_kernel(system_parent, distance) volume += kernel_weight v_b = current_velocity(v_parent, system_parent, neighbor) - p_b = particle_pressure(v_parent, system_parent, neighbor) - rho_b = particle_density(v_parent, system_parent, neighbor) + p_b = particle_pressure_parent(v_parent, system_parent, neighbor) + rho_b = particle_density_parent(v_parent, system_parent, neighbor) for dim in 1:ndims(system_child) - v_child[dim, system_child] += kernel_weight * v_b[dim] + v_child[dim, absolute_index] += kernel_weight * v_b[dim] end rho_a += kernel_weight * rho_b @@ -295,7 +299,7 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref end for dim in 1:ndims(system_child) - v[dim, system_child] ./ volume + v_child[dim, absolute_index] ./ volume end rho_a /= volume @@ -309,77 +313,66 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref return system_child end -@inline Base.resize!(system::System, ::Nothing) = system +@inline resize_system!(system::System, ::Nothing) = system -@inline function Base.resize!(system::System, particle_refinement::ParticleRefinement) +@inline function resize_system!(system::System, particle_refinement::ParticleRefinement) (; candidates, system_child) = particle_refinement if !isempty(candidates) - n_new_child = length(candidates) * nchilds(system_parent, particle_refinement) - capacity_parent = nparticles(system_parent) - length(candidates) + n_new_child = length(candidates) * nchilds(system, particle_refinement) + capacity_parent = nparticles(system) - length(candidates) capacity_child = nparticles(system_child) + n_new_child # Resize child system (extending) - resize!(system_child, capacity_child) + resize_system!(system_child, capacity_child) # Resize parent system (reducing) - resize!(system_parent, capacity_parent) + resize_system!(system, capacity_parent) end end -function Base.resize!(system::WeaklyCompressibleSPHSystem, capacity) +function resize_system!(system::FluidSystem, capacity::Int) (; mass, pressure, cache, density_calculator) = system - resize!(mass, capacity) - resize!(pressure, capacity) + Base.resize!(mass, capacity) + Base.resize!(pressure, capacity) resize!(cache, capacity, density_calculator) end -resize!(cache, capacity, ::SummationDensity) = resize!(cache.density, capacity) -resize!(cache, capacity, ::ContinuityDensity) = cache +resize!(cache, capacity::Int, ::SummationDensity) = Base.resize!(cache.density, capacity) +resize!(cache, capacity::Int, ::ContinuityDensity) = cache -@inline function _wrap_u(_u_ode, system, callback) - (; ranges_u_cache) = callback +@inline function _wrap_u(_u_ode, system, semi, callback) + (; ranges_u_cache, nparticles_cache) = callback - range = ranges_u_cache[system_indices(system, semi)] + range = ranges_u_cache[system_indices(system, semi)][1] + n_particles = nparticles_cache[system_indices(system, semi)] - @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system) + @boundscheck @assert length(range) == u_nvariables(system) * n_particles # This is a non-allocating version of: # return unsafe_wrap(Array{eltype(_u_ode), 2}, pointer(view(_u_ode, range)), - # (u_nvariables(system), n_moving_particles(system))) + # (u_nvariables(system), n_particles)) return PtrArray(pointer(view(_u_ode, range)), - (StaticInt(u_nvariables(system)), n_moving_particles(system))) + (StaticInt(u_nvariables(system)), n_particles)) end -@inline function _wrap_v(_v_ode, system, callback) - (; ranges_v_cache) = callback +@inline function _wrap_v(_v_ode, system, semi, callback) + (; ranges_v_cache, nparticles_cache) = callback - range = ranges_v_cache[system_indices(system, semi)] + range = ranges_v_cache[system_indices(system, semi)][1] + n_particles = nparticles_cache[system_indices(system, semi)] - @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) + @boundscheck @assert length(range) == v_nvariables(system) * n_particles # This is a non-allocating version of: # return unsafe_wrap(Array{eltype(_v_ode), 2}, pointer(view(_v_ode, range)), - # (v_nvariables(system), n_moving_particles(system))) + # (v_nvariables(system), n_particles)) return PtrArray(pointer(view(_v_ode, range)), - (StaticInt(v_nvariables(system)), n_moving_particles(system))) -end - -# `v_new` >= `v_old` -function copy_values_v!(v_new, v_old, system, ::Nothing, semi, callback) - (; each_particle_cache) = callback - - for particle in each_particle_cache[system_indices(system, semi)] - for i in 1:v_nvariables(system) - v_new[i, particle] = v_old[i, particle] - end - end + (StaticInt(v_nvariables(system)), n_particles)) end -# `v_new` <= `v_old` -function copy_values_v!(v_new, v_old, system, particle_refinement::ParticleRefinement, - semi, callback) +function copy_values_v!(v_new, v_old, system, semi, callback) (; eachparticle_cache) = callback # Copy only non-refined particles @@ -392,26 +385,13 @@ function copy_values_v!(v_new, v_old, system, particle_refinement::ParticleRefin end end -# `u_new` >= `u_old` -function copy_values_u!(u_new, u_old, system, ::Nothing, semi, callback) - (; each_particle_cache) = callback - - for particle in each_particle_cache[system_indices(system, semi)] - for i in 1:u_nuariables(system) - u_new[i, particle] = u_old[i, particle] - end - end -end - -# `u_new` <= `u_old` -function copy_values_u!(u_new, u_old, system, particle_refinement::ParticleRefinement, - semi, callback) +function copy_values_u!(u_new, u_old, system, semi, callback) (; eachparticle_cache) = callback # Copy only non-refined particles new_particle_id = 1 for particle in eachparticle_cache[system_indices(system, semi)] - for i in 1:u_nuariables(system) + for i in 1:u_nvariables(system) u_new[i, new_particle_id] = u_old[i, particle] end new_particle_id += 1 diff --git a/src/general/refinement_pattern.jl b/src/general/refinement_pattern.jl index 12d395c0e..c77ab088c 100644 --- a/src/general/refinement_pattern.jl +++ b/src/general/refinement_pattern.jl @@ -1,3 +1,5 @@ +abstract type RefinementCriteria{NDIMS, ELTYPE} end + struct CubicSplitting{ELTYPE} epsilon :: ELTYPE alpha :: ELTYPE @@ -7,8 +9,9 @@ struct CubicSplitting{ELTYPE} end end -function relative_positions(refinement_pattern::CubicSplitting, ::System{2}, - particle_spacing) +function (refinement_pattern::CubicSplitting)(system::System{2}) + (; initial_condition) = system + (; particle_spacing) = initial_condition (; epsilon) = refinement_pattern direction_1 = normalize([1.0, 1.0]) @@ -26,16 +29,68 @@ end # TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? # See https://www.sciencedirect.com/science/article/pii/S0020740319317023 -@inline nchilds(system, refinement_pattern) = 2^ndims(system) +@inline nchilds(system, rp::CubicSplitting) = 2^ndims(system) -@inline smoothing_length_child(system, refinement_pattern) = system.smoothing_length +@inline mass_child(system, mass, rp::CubicSplitting) = mass / nchilds(system, rp) -@inline mass_child(system, refinement_pattern) = system.mass +@inline smoothing_length_child(system, refinement_pattern) = refinement_pattern.alpha * + system.smoothing_length -@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing +@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing * + refinement_pattern.epsilon # ==== Refinement criteria -struct RefinementZone end # TODO +struct RefinementZone{NDIMS, ELTYPE} <: RefinementCriteria{NDIMS, ELTYPE} + zone_origin :: SVector{NDIMS, ELTYPE} + spanning_set :: Vector{SVector} + + function RefinementZone(plane_points, zone_width) + NDIMS = length(plane_points) + ELTYPE = typeof(zone_width) + + # Vectors spanning the zone. + spanning_set = spanning_vectors(plane_points, zone_width) + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + zone_origin = SVector(plane_points[1]...) + + return new{NDIMS, ELTYPE}(zone_origin, spanning_set_) + end +end + +@inline Base.ndims(::RefinementCriteria{NDIMS}) where {NDIMS} = NDIMS +@inline Base.eltype(::RefinementCriteria{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE + +function spanning_vectors(plane_points, zone_width) + + # Convert to tuple + return spanning_vectors(tuple(plane_points...), zone_width) +end + +function spanning_vectors(plane_points::NTuple{2}, zone_width) + plane_size = plane_points[2] - plane_points[1] + + # Calculate normal vector of plane + b = Vector(normalize([-plane_size[2]; plane_size[1]]) * zone_width) + + return hcat(b, plane_size) +end + +function spanning_vectors(plane_points::NTuple{3}, zone_width) + # Vectors spanning the plane + edge1 = plane_points[2] - plane_points[1] + edge2 = plane_points[3] - plane_points[1] + + if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) + throw(ArgumentError("the provided points do not span a rectangular plane")) + end + + # Calculate normal vector of plane + c = Vector(normalize(cross(edge2, edge1)) * zone_width) + + return hcat(c, edge1, edge2) +end @inline function (refinement_criterion::RefinementZone)(system, particle, v, u, v_ode, u_ode, semi) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index f36ecca92..8da937db3 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -65,9 +65,9 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear # Other checks might be added here later. check_configuration(systems) - #systems = create_subsystems(systems) + systems = create_system_childs(systems) - ranges_u, ranges_v = ranges_uv(systems) + ranges_v, ranges_u = ranges_vu(systems) # Create (and initialize) a tuple of n neighborhood searches for each of the n systems # We will need one neighborhood search for each pair of systems. @@ -81,7 +81,7 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear return Semidiscretization(systems, ranges_u, ranges_v, searches) end -function ranges_uv(systems) +function ranges_vu(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])] @@ -91,7 +91,7 @@ function ranges_uv(systems) ranges_v = Tuple([(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i])] for i in eachindex(sizes_v)) - return ranges_u, ranges_v + return ranges_v, ranges_u end # Inline show function e.g. Semidiscretization(neighborhood_search=...) From fe4d47dad8128e8cdb37c43fd8da28a6efd29530 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 18 Mar 2024 20:25:47 +0100 Subject: [PATCH 08/20] fix dispatching --- src/general/refinement.jl | 52 ++++++++++++++++++++++------------ src/schemes/boundary/system.jl | 18 ++++++++++++ 2 files changed, 52 insertions(+), 18 deletions(-) diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 83c7c9ed6..983a8db45 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -58,12 +58,17 @@ end function create_system_childs(systems) systems_ = () foreach_system(systems) do system - systems_ = (systems_..., create_system_child(system, system.particle_refinement)...) + systems_ = (systems_..., create_system_child(system)...) end return (systems..., systems_...) end +create_system_child(system) = () +function create_system_child(system::FluidSystem) + create_system_child(system, system.particle_refinement) +end + create_system_child(system, ::Nothing) = () function create_system_child(system::WeaklyCompressibleSPHSystem, @@ -107,7 +112,7 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, particle_refinement=particle_refinement_) # Empty mass vector leads to `nparticles(system_child) = 0` - Base.resize!(system_child.mass, 0) + resize!(system_child.mass, 0) particle_refinement.system_child = system_child @@ -164,15 +169,15 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) callback.nparticles_cache = Tuple(n_moving_particles(system) for system in semi.systems) # Resize internal storage - Base.resize!(_v_ode, length(v_ode)) - Base.resize!(_u_ode, length(u_ode)) + resize!(_v_ode, length(v_ode)) + resize!(_u_ode, length(u_ode)) _v_ode .= v_ode _u_ode .= u_ode # Resize all systems foreach_system(semi) do system - resize_system!(system, system.particle_refinement) + resize_system!(system) end # Set `resize!`d ranges @@ -182,12 +187,12 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) semi.ranges_u[i][1] = ranges_u_tmp[i][1] end - sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in semi.systems) sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in semi.systems) + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in semi.systems) # Resize integrated values - Base.resize!(v_ode, sum(sizes_v)) - Base.resize!(u_ode, sum(sizes_u)) + resize!(v_ode, sum(sizes_v)) + resize!(u_ode, sum(sizes_u)) # Preserve non-changing values foreach_system(semi) do system @@ -207,10 +212,16 @@ function refine_particles!(callback, semi, v_ode, u_ode) # Refine particles in all systems foreach_system(semi) do system - refine_particles!(system, system.particle_refinement, v_ode, u_ode, callback, semi) + refine_particles!(system, v_ode, u_ode, callback, semi) end end +refine_particles!(system, v_ode, u_ode, callback, semi) = system + +function refine_particles!(system::FluidSystem, v_ode, u_ode, callback, semi) + refine_particles!(system, system.particle_refinement, v_ode, u_ode, callback, semi) +end + refine_particles!(system, ::Nothing, v_ode, u_ode, callback, semi) = system function refine_particles!(system_parent, particle_refinement::ParticleRefinement, @@ -313,9 +324,14 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref return system_child end -@inline resize_system!(system::System, ::Nothing) = system +@inline resize_system!(system) = system +@inline resize_system!(system::FluidSystem) = resize_system!(system, + system.particle_refinement) -@inline function resize_system!(system::System, particle_refinement::ParticleRefinement) +@inline resize_system!(system::FluidSystem, ::Nothing) = system + +@inline function resize_system!(system::FluidSystem, + particle_refinement::ParticleRefinement) (; candidates, system_child) = particle_refinement if !isempty(candidates) @@ -334,13 +350,13 @@ end function resize_system!(system::FluidSystem, capacity::Int) (; mass, pressure, cache, density_calculator) = system - Base.resize!(mass, capacity) - Base.resize!(pressure, capacity) - resize!(cache, capacity, density_calculator) + resize!(mass, capacity) + resize!(pressure, capacity) + resize_cache!(cache, capacity, density_calculator) end -resize!(cache, capacity::Int, ::SummationDensity) = Base.resize!(cache.density, capacity) -resize!(cache, capacity::Int, ::ContinuityDensity) = cache +resize_cache!(cache, capacity::Int, ::SummationDensity) = resize!(cache.density, capacity) +resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache @inline function _wrap_u(_u_ode, system, semi, callback) (; ranges_u_cache, nparticles_cache) = callback @@ -397,8 +413,8 @@ function copy_values_u!(u_new, u_old, system, semi, callback) new_particle_id += 1 end end - -@inline get_iterator(system) = get_iterator(system, system.particle_refinement) +@inline get_iterator(system) = eachparticle(system) +@inline get_iterator(system::FluidSystem) = get_iterator(system, system.particle_refinement) @inline get_iterator(system, ::Nothing) = eachparticle(system) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 6cb5b7bd0..35f6187f3 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -268,6 +268,24 @@ function write_v0!(v0, return v0 end +copy_values_v!(v_new, v_old, system::BoundarySPHSystem, semi, callback) = v_new +copy_values_u!(u_new, u_old, system::BoundarySPHSystem, semi, callback) = u_new + +function copy_values_v!(v_new, v_old, + system::BoundarySPHSystem{<:BoundaryModelDummyParticles{ContinuityDensity}}, + semi, callback) + (; eachparticle_cache) = callback + + # Copy only non-refined particles + new_particle_id = 1 + for particle in eachparticle_cache[system_indices(system, semi)] + for i in 1:v_nvariables(system) + v_new[i, new_particle_id] = v_old[i, particle] + end + new_particle_id += 1 + end +end + function restart_with!(system::BoundarySPHSystem, v, u) return system end From 92696c7a8b8f6996d5f66aa41d0859d2de8c4432 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 19 Mar 2024 16:28:48 +0100 Subject: [PATCH 09/20] first working prototype --- src/TrixiParticles.jl | 4 +- src/callbacks/refinement_callback.jl | 35 +++++++++--- src/general/refinement.jl | 79 +++++++++++++------------- src/general/semidiscretization.jl | 13 +++-- src/neighborhood_search/grid_nhs.jl | 12 +++- src/neighborhood_search/trivial_nhs.jl | 8 ++- src/schemes/fluid/fluid.jl | 4 +- src/visualization/write2vtk.jl | 44 +++++++------- 8 files changed, 121 insertions(+), 78 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 5ea61a880..eb045656e 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -16,7 +16,8 @@ using Polyester: Polyester, @batch using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @series using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, - get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem + get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, + RecursiveArrayTools @reexport using StaticArrays: SVector using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt @@ -71,5 +72,6 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk +export ParticleRefinement, RefinementZone, CubicSplitting end # module diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl index e342581a1..a571b73c8 100644 --- a/src/callbacks/refinement_callback.jl +++ b/src/callbacks/refinement_callback.jl @@ -4,10 +4,6 @@ mutable struct ParticleRefinementCallback{I} ranges_v_cache :: Tuple nparticles_cache :: Tuple eachparticle_cache :: Tuple - - # internal `resize!`able storage - _u_ode::Vector{Float64} - _v_ode::Vector{Float64} end function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) @@ -24,7 +20,7 @@ function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) interval = 1 end - refinement_callback = ParticleRefinementCallback(interval, (), (), (), (), [0.0], [0.0]) + refinement_callback = ParticleRefinementCallback(interval, (), (), (), ()) if dt > 0 # Add a `tstop` every `dt`, and save the final solution. @@ -70,9 +66,32 @@ function (refinement_callback::ParticleRefinementCallback)(integrator) semi = integrator.p v_ode, u_ode = integrator.u.x - refinement!(v_ode, u_ode, semi, refinement_callback) + # Update NHS + @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi) + + # Basically `get_tmp_cache(integrator)` to write into in order to be non-allocating + # https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Caches + v_tmp, u_tmp = integrator.cache.tmp.x + + v_tmp .= v_ode + u_tmp .= u_ode - update_systems_and_nhs(v_ode, u_ode, semi, t) + refinement!(v_ode, u_ode, v_tmp, u_tmp, semi, refinement_callback) + + # Resize neighborhood search + foreach_system(semi) do system + foreach_system(semi) do neighbor_system + search = get_neighborhood_search(system, neighbor_system, semi) + u_neighbor = wrap_u(u_ode, neighbor_system, semi) + + resize_nhs!(search, system, neighbor_system, u_neighbor) + end + end + + resize!(integrator, (length(v_ode), length(u_ode))) + + @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, + semi, t) # Tell OrdinaryDiffEq that u has been modified u_modified!(integrator, true) @@ -80,6 +99,8 @@ function (refinement_callback::ParticleRefinementCallback)(integrator) return integrator end +Base.resize!(a::RecursiveArrayTools.ArrayPartition, sizes::Tuple) = resize!.(a.x, sizes) + function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) @nospecialize cb # reduce precompilation time print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")") diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 983a8db45..72db358c3 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -121,14 +121,14 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, end # ==== Refinement -function refinement!(v_ode, u_ode, semi, callback) +function refinement!(v_ode, u_ode, _v_cache, _u_cache, semi, callback) foreach_system(semi) do system check_refinement_criteria!(system, v_ode, u_ode, semi) end - resize_and_copy!(callback, semi, v_ode, u_ode) + resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) - refine_particles!(callback, semi, v_ode, u_ode) + refine_particles!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) end check_refinement_criteria!(system, v_ode, u_ode, semi) = system @@ -160,21 +160,12 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem return system end -function resize_and_copy!(callback, semi, v_ode, u_ode) - (; _v_ode, _u_ode) = callback - +function resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) # Get non-`resize!`d ranges callback.ranges_v_cache, callback.ranges_u_cache = ranges_vu(semi.systems) callback.eachparticle_cache = Tuple(get_iterator(system) for system in semi.systems) callback.nparticles_cache = Tuple(n_moving_particles(system) for system in semi.systems) - # Resize internal storage - resize!(_v_ode, length(v_ode)) - resize!(_u_ode, length(u_ode)) - - _v_ode .= v_ode - _u_ode .= u_ode - # Resize all systems foreach_system(semi) do system resize_system!(system) @@ -198,8 +189,8 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) foreach_system(semi) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - _v = _wrap_v(_v_ode, system, semi, callback) - _u = _wrap_u(_u_ode, system, semi, callback) + _v = _wrap_v(_v_cache, system, semi, callback) + _u = _wrap_u(_u_cache, system, semi, callback) copy_values_v!(v, _v, system, semi, callback) copy_values_u!(u, _u, system, semi, callback) @@ -208,31 +199,36 @@ function resize_and_copy!(callback, semi, v_ode, u_ode) return callback end -function refine_particles!(callback, semi, v_ode, u_ode) +function refine_particles!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) # Refine particles in all systems foreach_system(semi) do system - refine_particles!(system, v_ode, u_ode, callback, semi) + refine_particles!(system, v_ode, u_ode, _v_cache, _u_cache, callback, semi) end end -refine_particles!(system, v_ode, u_ode, callback, semi) = system +refine_particles!(system, v_ode, u_ode, _v_cache, _u_cache, callback, semi) = system -function refine_particles!(system::FluidSystem, v_ode, u_ode, callback, semi) - refine_particles!(system, system.particle_refinement, v_ode, u_ode, callback, semi) +function refine_particles!(system::FluidSystem, v_ode, u_ode, _v_cache, _u_cache, + callback, semi) + refine_particles!(system, system.particle_refinement, v_ode, u_ode, _v_cache, _u_cache, + callback, semi) end -refine_particles!(system, ::Nothing, v_ode, u_ode, callback, semi) = system +function refine_particles!(system::FluidSystem, ::Nothing, v_ode, u_ode, _v_cache, _u_cache, + callback, semi) + return system +end -function refine_particles!(system_parent, particle_refinement::ParticleRefinement, - v_ode, u_ode, callback, semi) - (; _v_ode, _u_ode) = callback +function refine_particles!(system_parent::FluidSystem, + particle_refinement::ParticleRefinement, + v_ode, u_ode, _v_cache, _u_cache, callback, semi) (; candidates, system_child) = particle_refinement if !isempty(candidates) # Old storage - v_parent = _wrap_v(_v_ode, system_parent, semi, callback) - u_parent = _wrap_u(_u_ode, system_parent, semi, callback) + v_parent = _wrap_v(_v_cache, system_parent, semi, callback) + u_parent = _wrap_u(_u_cache, system_parent, semi, callback) # Resized storage v_child = wrap_v(v_ode, system_child, semi) @@ -309,16 +305,18 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref end end - for dim in 1:ndims(system_child) - v_child[dim, absolute_index] ./ volume - end + if volume > eps() + for dim in 1:ndims(system_child) + v_child[dim, absolute_index] /= volume + end - rho_a /= volume - p_a /= volume + rho_a /= volume + p_a /= volume - set_particle_density(particle_child, v_child, system_child.density_calculator, - system_child, rho_a) - set_particle_pressure(particle_child, v_child, system_child, p_a) + set_particle_density(absolute_index, v_child, system_child.density_calculator, + system_child, rho_a) + set_particle_pressure(absolute_index, v_child, system_child, p_a) + end end return system_child @@ -358,7 +356,7 @@ end resize_cache!(cache, capacity::Int, ::SummationDensity) = resize!(cache.density, capacity) resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache -@inline function _wrap_u(_u_ode, system, semi, callback) +@inline function _wrap_u(_u_cache, system, semi, callback) (; ranges_u_cache, nparticles_cache) = callback range = ranges_u_cache[system_indices(system, semi)][1] @@ -367,13 +365,13 @@ resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache @boundscheck @assert length(range) == u_nvariables(system) * n_particles # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(_u_ode), 2}, pointer(view(_u_ode, range)), + # return unsafe_wrap(Array{eltype(_u_cache), 2}, pointer(view(_u_cache, range)), # (u_nvariables(system), n_particles)) - return PtrArray(pointer(view(_u_ode, range)), + return PtrArray(pointer(view(_u_cache, range)), (StaticInt(u_nvariables(system)), n_particles)) end -@inline function _wrap_v(_v_ode, system, semi, callback) +@inline function _wrap_v(_v_cache, system, semi, callback) (; ranges_v_cache, nparticles_cache) = callback range = ranges_v_cache[system_indices(system, semi)][1] @@ -382,9 +380,9 @@ end @boundscheck @assert length(range) == v_nvariables(system) * n_particles # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(_v_ode), 2}, pointer(view(_v_ode, range)), + # return unsafe_wrap(Array{eltype(_v_cache), 2}, pointer(view(_v_cache, range)), # (v_nvariables(system), n_particles)) - return PtrArray(pointer(view(_v_ode, range)), + return PtrArray(pointer(view(_v_cache, range)), (StaticInt(v_nvariables(system)), n_particles)) end @@ -413,6 +411,7 @@ function copy_values_u!(u_new, u_old, system, semi, callback) new_particle_id += 1 end end + @inline get_iterator(system) = eachparticle(system) @inline get_iterator(system::FluidSystem) = get_iterator(system, system.particle_refinement) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 8da937db3..5a483289a 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -82,15 +82,16 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear end function ranges_vu(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)) + 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)) + return ranges_v, ranges_u end @@ -147,8 +148,8 @@ function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSear end @inline function compact_support(system, neighbor) - (; smoothing_kernel, smoothing_length) = system - return compact_support(smoothing_kernel, smoothing_length) + smoothing_length_avg = average_smoothing_length(system, neighbor) + return compact_support(system.smoothing_kernel, smoothing_length_avg) end @inline function compact_support(system::TotalLagrangianSPHSystem, diff --git a/src/neighborhood_search/grid_nhs.jl b/src/neighborhood_search/grid_nhs.jl index dbaed0180..470502392 100644 --- a/src/neighborhood_search/grid_nhs.jl +++ b/src/neighborhood_search/grid_nhs.jl @@ -87,7 +87,7 @@ since not sorting makes our implementation a lot faster (although less paralleli In: Computer Graphics Forum 30.1 (2011), pages 99–112. [doi: 10.1111/J.1467-8659.2010.01832.X](https://doi.org/10.1111/J.1467-8659.2010.01832.X) """ -struct GridNeighborhoodSearch{NDIMS, ELTYPE, PB} +mutable struct GridNeighborhoodSearch{NDIMS, ELTYPE, PB} hashtable :: Dict{NTuple{NDIMS, Int}, Vector{Int}} search_radius :: ELTYPE empty_vector :: Vector{Int} # Just an empty vector (used in `eachneighbor`) @@ -405,3 +405,13 @@ end function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, u) return nhs end + +function resize_nhs!(nhs::GridNeighborhoodSearch, system, neighbor_system, u_neighbor) + NDIMS = ndims(neighbor_system) + n_particles = nparticles(neighbor_system) + nhs.cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_particles, Threads.nthreads()) + + initialize!(nhs, nhs_coords(system, neighbor_system, u_neighbor)) + + return nhs +end diff --git a/src/neighborhood_search/trivial_nhs.jl b/src/neighborhood_search/trivial_nhs.jl index 3902f1d3f..0b7c14998 100644 --- a/src/neighborhood_search/trivial_nhs.jl +++ b/src/neighborhood_search/trivial_nhs.jl @@ -60,7 +60,7 @@ internal function `eachneighbor`. └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} +mutable struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} search_radius :: ELTYPE eachparticle :: EP periodic_box :: PB @@ -95,3 +95,9 @@ end @inline initialize!(search::TrivialNeighborhoodSearch, coords_fun) = search @inline update!(search::TrivialNeighborhoodSearch, coords_fun) = search @inline eachneighbor(coords, search::TrivialNeighborhoodSearch) = search.eachparticle + +function resize_nhs!(search::TrivialNeighborhoodSearch, system, neighbor, u_neighbor) + search.eachparticle = eachparticle(neighbor) + + return search +end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 514791eec..79a2aa073 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -84,7 +84,7 @@ include("entropically_damped_sph/entropically_damped_sph.jl") # and density correction. It cannot be used to set up an initial condition, # as the particle density depends on the particle positions. -@inline set_particle_density(particle, v, ::SummationDensity, system, density) = nothing +@inline set_particle_density(particle, v, ::SummationDensity, system, density) = particle @inline function set_particle_density(particle, v, ::ContinuityDensity, system::WeaklyCompressibleSPHSystem, density) @@ -96,7 +96,7 @@ end v[end - 1, particle] = density end -@inline set_particle_pressure(particle, v, system, pressure) = nothing +@inline set_particle_pressure(particle, v, system, pressure) = particle @inline function set_particle_pressure(particle, v, system::EntropicallyDampedSPHSystem, pressure) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 9db300269..02ffb1f92 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -88,35 +88,38 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix points = periodic_coords(current_coordinates(u, system), periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] - if abs(maximum(points)) > max_coordinates || abs(minimum(points)) > max_coordinates - println("Warning: At least one particle's absolute coordinates exceed `max_coordinates`" - * - " and have been clipped") - for i in eachindex(points) - points[i] = clamp(points[i], -max_coordinates, max_coordinates) + if !isempty(points) + if abs(maximum(points)) > max_coordinates || abs(minimum(points)) > max_coordinates + println("Warning: At least one particle's absolute coordinates exceed `max_coordinates`" + * + " and have been clipped") + for i in eachindex(points) + points[i] = clamp(points[i], -max_coordinates, max_coordinates) + end end end vtk_grid(file, points, cells) do vtk - write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) + if !isempty(points) + write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) - # Store particle index - vtk["index"] = eachparticle(system) - vtk["time"] = t + # Store particle index + vtk["index"] = eachparticle(system) + vtk["time"] = t - if write_meta_data - vtk["solver_version"] = get_git_hash() - vtk["julia_version"] = string(VERSION) - end + if write_meta_data + vtk["solver_version"] = get_git_hash() + vtk["julia_version"] = string(VERSION) + end - # Extract custom quantities for this system - for (key, quantity) in custom_quantities - value = custom_quantity(quantity, v, u, t, system) - if value !== nothing - vtk[string(key)] = value + # Extract custom quantities for this system + for (key, quantity) in custom_quantities + value = custom_quantity(quantity, v, u, t, system) + if value !== nothing + vtk[string(key)] = value + end end end - # Add to collection pvd[t] = vtk end @@ -178,6 +181,7 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) for particle in eachparticle(system)] vtk["pressure"] = [particle_pressure(v, system, particle) for particle in eachparticle(system)] + vtk["mass"] = system.mass if write_meta_data vtk["acceleration"] = system.acceleration From c9d0e736036c5e4cd199bd140e0523fd3d79a3d1 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 09:29:53 +0100 Subject: [PATCH 10/20] add edac and fix bugs --- src/TrixiParticles.jl | 2 +- src/callbacks/refinement_callback.jl | 8 +- src/general/refinement.jl | 84 +++++++++++++------ src/general/semidiscretization.jl | 1 - src/general/system.jl | 2 - .../fluid/entropically_damped_sph/system.jl | 10 ++- 6 files changed, 71 insertions(+), 36 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index eb045656e..2e8bddfac 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -46,7 +46,7 @@ export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, BoundarySPHSystem export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback + PostprocessCallback, StepsizeCallback, ParticleRefinementCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl index a571b73c8..e2e1edfab 100644 --- a/src/callbacks/refinement_callback.jl +++ b/src/callbacks/refinement_callback.jl @@ -88,11 +88,11 @@ function (refinement_callback::ParticleRefinementCallback)(integrator) end end - resize!(integrator, (length(v_ode), length(u_ode))) - @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, semi, t) + resize!(integrator, (length(v_ode), length(u_ode))) + # Tell OrdinaryDiffEq that u has been modified u_modified!(integrator, true) @@ -122,7 +122,7 @@ function Base.show(io::IO, ::MIME"text/plain", else refinement_cb = cb.affect! setup = [ - "interval" => refinement_cb.interval + "interval" => refinement_cb.interval, ] summary_box(io, "ParticleRefinementCallback", setup) end @@ -138,7 +138,7 @@ function Base.show(io::IO, ::MIME"text/plain", else refinement_cb = cb.affect!.affect! setup = [ - "dt" => refinement_cb.interval + "dt" => refinement_cb.interval, ] summary_box(io, "ParticleRefinementCallback", setup) end diff --git a/src/general/refinement.jl b/src/general/refinement.jl index 72db358c3..63e0481f5 100644 --- a/src/general/refinement.jl +++ b/src/general/refinement.jl @@ -7,11 +7,11 @@ include("refinement_pattern.jl") mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} - candidates :: Vector{Int} - refinement_levels :: Int - refinement_pattern :: RP - refinement_criteria :: RC - available_childs :: Int + candidates :: Vector{Int} + refinement_pattern :: RP + refinement_criteria :: RC + criteria_next_levels :: Vector{RC} + available_childs :: Int # Depends on refinement pattern, particle spacing and parameters ϵ and α. # Should be obtained prior to simulation in `create_system_child()` @@ -24,24 +24,26 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} # API --> parent system with `RL=0` function ParticleRefinement(refinement_criteria...; refinement_pattern=CubicSplitting(), - refinement_levels=1) + criteria_next_levels=[]) ELTYPE = eltype(refinement_criteria[1]) NDIMS = ndims(refinement_criteria[1]) return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([0], refinement_levels, refinement_pattern, - refinement_criteria, 0) + typeof(refinement_criteria)}([], refinement_pattern, + refinement_criteria, + criteria_next_levels, 0) end # Internal constructor for multiple refinement levels function ParticleRefinement{RL}(refinement_criteria::Tuple, - refinement_pattern, refinement_levels) where {RL} + refinement_pattern, criteria_next_levels) where {RL} ELTYPE = eltype(refinement_criteria[1]) NDIMS = ndims(refinement_criteria[1]) return new{RL, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([0], refinement_levels, refinement_pattern, - refinement_criteria, 0) + typeof(refinement_criteria)}([], refinement_pattern, + refinement_criteria, + criteria_next_levels, 0) end end @@ -71,12 +73,18 @@ end create_system_child(system, ::Nothing) = () -function create_system_child(system::WeaklyCompressibleSPHSystem, +function create_system_child(system::FluidSystem, particle_refinement::ParticleRefinement{RL}) where {RL} - (; refinement_levels, refinement_pattern, refinement_criteria) = particle_refinement - (; density_calculator, state_equation, smoothing_kernel, - pressure_acceleration_formulation, viscosity, density_diffusion, - acceleration, correction, source_terms) = system + (; criteria_next_levels, refinement_pattern, refinement_criteria) = particle_refinement + + if system isa WeaklyCompressibleSPHSystem + (; density_calculator, state_equation, smoothing_kernel, + pressure_acceleration_formulation, viscosity, density_diffusion, + acceleration, correction, source_terms) = system + else + (; density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, + pressure_acceleration_formulation, acceleration, correction, source_terms) = system + end NDIMS = ndims(system) @@ -97,19 +105,32 @@ function create_system_child(system::WeaklyCompressibleSPHSystem, # Let recursive dispatch handle multiple refinement levels level = RL + 1 - particle_refinement_ = if level == refinement_levels + particle_refinement_ = if isempty(criteria_next_levels) nothing else + refinement_criteria = first(criteria_next_levels) ParticleRefinement{level}(refinement_criteria, refinement_pattern, - refinement_levels) + criteria_next_levels[(level + 1):end]) end - system_child = WeaklyCompressibleSPHSystem(empty_ic, density_calculator, state_equation, - smoothing_kernel, smoothing_length_; - pressure_acceleration=pressure_acceleration_formulation, - viscosity, density_diffusion, acceleration, - correction, source_terms, - particle_refinement=particle_refinement_) + if system isa WeaklyCompressibleSPHSystem + system_child = WeaklyCompressibleSPHSystem(empty_ic, density_calculator, + state_equation, + smoothing_kernel, smoothing_length_; + pressure_acceleration=pressure_acceleration_formulation, + viscosity, density_diffusion, + acceleration, + correction, source_terms, + particle_refinement=particle_refinement_) + else + alpha = nu_edac * 8 / (smoothing_length_ * sound_speed) + system_child = EntropicallyDampedSPHSystem(empty_ic, smoothing_kernel, + smoothing_length_, sound_speed; + pressure_acceleration=pressure_acceleration_formulation, + density_calculator, alpha, viscosity, + acceleration, source_terms, + particle_refinement=particle_refinement_) + end # Empty mass vector leads to `nparticles(system_child) = 0` resize!(system_child.mass, 0) @@ -281,6 +302,10 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref p_a = zero(eltype(system_child)) rho_a = zero(eltype(system_child)) + for dim in 1:ndims(system_child) + v_child[dim, absolute_index] = zero(eltype(system_child)) + end + for neighbor in eachneighbor(child_coords, nhs) neighbor_coords = current_coords(u_parent, system_parent, neighbor) pos_diff = child_coords - neighbor_coords @@ -337,6 +362,8 @@ end capacity_parent = nparticles(system) - length(candidates) capacity_child = nparticles(system_child) + n_new_child + capacity_parent <= 0 && error("`RefinementCriteria` affects all particles") + # Resize child system (extending) resize_system!(system_child, capacity_child) @@ -345,7 +372,7 @@ end end end -function resize_system!(system::FluidSystem, capacity::Int) +function resize_system!(system::WeaklyCompressibleSPHSystem, capacity::Int) (; mass, pressure, cache, density_calculator) = system resize!(mass, capacity) @@ -353,6 +380,13 @@ function resize_system!(system::FluidSystem, capacity::Int) resize_cache!(cache, capacity, density_calculator) end +function resize_system!(system::EntropicallyDampedSPHSystem, capacity::Int) + (; mass, cache, density_calculator) = system + + resize!(mass, capacity) + resize_cache!(cache, capacity, density_calculator) +end + resize_cache!(cache, capacity::Int, ::SummationDensity) = resize!(cache.density, capacity) resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 5a483289a..3a81b6c96 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -441,7 +441,6 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) update_final!(system, v, u, v_ode, u_ode, semi, t) end - end function update_nhs(u_ode, semi) diff --git a/src/general/system.jl b/src/general/system.jl index 69934765a..c201bf3a5 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -19,8 +19,6 @@ initialize!(system, neighborhood_search) = system @inline eachparticle(system) = Base.OneTo(nparticles(system)) @inline each_moving_particle(system) = Base.OneTo(n_moving_particles(system)) -@inline refinement!(system, v, u, v_ode, u_ode, semi) = system - # This should not be dispatched by system type. We always expect to get a column of `A`. @inline function extract_svector(A, system, i) extract_svector(A, Val(ndims(system)), i) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index bfa07c051..0210f3344 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -40,7 +40,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, - PF, ST, C} <: FluidSystem{NDIMS} + PF, ST, PR, C} <: FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] density_calculator :: DC @@ -53,6 +53,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, correction :: Nothing pressure_acceleration_formulation :: PF source_terms :: ST + particle_refinement :: PR cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, @@ -62,6 +63,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), + particle_refinement=nothing, source_terms=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -87,10 +89,12 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, cache = create_cache_density(initial_condition, density_calculator) new{NDIMS, ELTYPE, typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity), typeof(pressure_acceleration), typeof(source_terms), + typeof(viscosity), typeof(pressure_acceleration), + typeof(source_terms), typeof(particle_refinement), typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, source_terms, cache) + nothing, pressure_acceleration, source_terms, + particle_refinement, cache) end end From f75c76391bf8972f79b5c72ad8cc4f63218d0cbb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 11:29:11 +0100 Subject: [PATCH 11/20] undo mean smoothing length --- src/general/density_calculators.jl | 5 +---- src/general/semidiscretization.jl | 13 +++++++++---- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 3b4f257b9..ee89c92a6 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -52,10 +52,7 @@ function summation_density!(system, semi, u, u_ode, density; for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, nhs, particles=particles) do particle, neighbor, pos_diff, distance mass = hydrodynamic_mass(neighbor_system, neighbor) - - h_mean = (system.smoothing_length + neighbor_system.smoothing_length)/2 - #density[particle] += mass * smoothing_kernel(system, distance) - density[particle] += mass * kernel(system.smoothing_kernel, distance, h_mean) + density[particle] += mass * smoothing_kernel(system, distance) end end end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 61814827f..fd7793455 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -97,6 +97,7 @@ function ranges_vu(systems) ranges_u = Tuple([(sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i])] for i in eachindex(sizes_u)) + # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct return ranges_v, ranges_u end @@ -155,8 +156,8 @@ function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSear end @inline function compact_support(system, neighbor) - smoothing_length_avg = average_smoothing_length(system, neighbor) - return compact_support(system.smoothing_kernel, smoothing_length_avg) + (; smoothing_kernel, smoothing_length) = system + return compact_support(smoothing_kernel, smoothing_length) end @inline function compact_support(system::TotalLagrangianSPHSystem, @@ -333,7 +334,9 @@ end @inline function wrap_u(u_ode, system, semi) (; ranges_u) = semi - range = ranges_u[system_indices(system, semi)][1] + # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct. + # The `Vector` has length 1, thus, take the first element. + range = first(ranges_u[system_indices(system, semi)]) @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system) @@ -347,7 +350,9 @@ end @inline function wrap_v(v_ode, system, semi) (; ranges_v) = semi - range = ranges_v[system_indices(system, semi)][1] + # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct. + # The `Vector` has length 1, thus, take the first element. + range = first(ranges_v[system_indices(system, semi)]) @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) From 16412031a0cef5f8296eb7e17178e54f9713b459 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 11:44:25 +0100 Subject: [PATCH 12/20] move to separate folder --- src/TrixiParticles.jl | 2 +- .../particle_refinement.jl | 1 + .../refinement.jl | 7 +-- .../refinement_criteria.jl} | 46 +++---------------- src/particle_refinement/refinement_pattern.jl | 38 +++++++++++++++ 5 files changed, 47 insertions(+), 47 deletions(-) create mode 100644 src/particle_refinement/particle_refinement.jl rename src/{general => particle_refinement}/refinement.jl (99%) rename src/{general/refinement_pattern.jl => particle_refinement/refinement_criteria.jl} (60%) create mode 100644 src/particle_refinement/refinement_pattern.jl diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 2121ecb4b..cd415c2b2 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -39,7 +39,7 @@ include("schemes/schemes.jl") include("general/semidiscretization.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") -include("general/refinement.jl") +include("particle_refinement/particle_refinement.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition diff --git a/src/particle_refinement/particle_refinement.jl b/src/particle_refinement/particle_refinement.jl new file mode 100644 index 000000000..3a6a7e965 --- /dev/null +++ b/src/particle_refinement/particle_refinement.jl @@ -0,0 +1 @@ +include("refinement.jl") diff --git a/src/general/refinement.jl b/src/particle_refinement/refinement.jl similarity index 99% rename from src/general/refinement.jl rename to src/particle_refinement/refinement.jl index 63e0481f5..96b0ef964 100644 --- a/src/general/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -1,10 +1,5 @@ -# Criteria of refinement: -# -# - fixed (/moving?) refinement zone -# - number of neighbors -# - problem specific criteria (e.g. high velocity gradient) - include("refinement_pattern.jl") +include("refinement_criteria.jl") mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} candidates :: Vector{Int} diff --git a/src/general/refinement_pattern.jl b/src/particle_refinement/refinement_criteria.jl similarity index 60% rename from src/general/refinement_pattern.jl rename to src/particle_refinement/refinement_criteria.jl index c77ab088c..3f28635c1 100644 --- a/src/general/refinement_pattern.jl +++ b/src/particle_refinement/refinement_criteria.jl @@ -1,45 +1,11 @@ -abstract type RefinementCriteria{NDIMS, ELTYPE} end - -struct CubicSplitting{ELTYPE} - epsilon :: ELTYPE - alpha :: ELTYPE - - function CubicSplitting(; epsilon=0.5, alpha=0.5) - new{typeof(epsilon)}(epsilon, alpha) - end -end - -function (refinement_pattern::CubicSplitting)(system::System{2}) - (; initial_condition) = system - (; particle_spacing) = initial_condition - (; epsilon) = refinement_pattern - - direction_1 = normalize([1.0, 1.0]) - direction_2 = normalize([1.0, -1.0]) - direction_3 = -direction_1 - direction_4 = -direction_2 +# Criteria of refinement: +# +# - fixed (/moving?) refinement zone +# - number of neighbors +# - problem specific criteria (e.g. high velocity gradient) - relative_position = hcat(particle_spacing * epsilon * direction_1, - particle_spacing * epsilon * direction_2, - particle_spacing * epsilon * direction_3, - particle_spacing * epsilon * direction_4) - - return reinterpret(reshape, SVector{2, typeof(epsilon)}, relative_position) -end - -# TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? -# See https://www.sciencedirect.com/science/article/pii/S0020740319317023 -@inline nchilds(system, rp::CubicSplitting) = 2^ndims(system) - -@inline mass_child(system, mass, rp::CubicSplitting) = mass / nchilds(system, rp) - -@inline smoothing_length_child(system, refinement_pattern) = refinement_pattern.alpha * - system.smoothing_length - -@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing * - refinement_pattern.epsilon +abstract type RefinementCriteria{NDIMS, ELTYPE} end -# ==== Refinement criteria struct RefinementZone{NDIMS, ELTYPE} <: RefinementCriteria{NDIMS, ELTYPE} zone_origin :: SVector{NDIMS, ELTYPE} spanning_set :: Vector{SVector} diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl new file mode 100644 index 000000000..d65d293df --- /dev/null +++ b/src/particle_refinement/refinement_pattern.jl @@ -0,0 +1,38 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + + function CubicSplitting(; epsilon=0.5, alpha=0.5) + new{typeof(epsilon)}(epsilon, alpha) + end +end + +function (refinement_pattern::CubicSplitting)(system::System{2}) + (; initial_condition) = system + (; particle_spacing) = initial_condition + (; epsilon) = refinement_pattern + + direction_1 = normalize([1.0, 1.0]) + direction_2 = normalize([1.0, -1.0]) + direction_3 = -direction_1 + direction_4 = -direction_2 + + relative_position = hcat(particle_spacing * epsilon * direction_1, + particle_spacing * epsilon * direction_2, + particle_spacing * epsilon * direction_3, + particle_spacing * epsilon * direction_4) + + return reinterpret(reshape, SVector{2, typeof(epsilon)}, relative_position) +end + +# TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? +# See https://www.sciencedirect.com/science/article/pii/S0020740319317023 +@inline nchilds(system, rp::CubicSplitting) = 2^ndims(system) + +@inline mass_child(system, mass, rp::CubicSplitting) = mass / nchilds(system, rp) + +@inline smoothing_length_child(system, refinement_pattern) = refinement_pattern.alpha * + system.smoothing_length + +@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing * + refinement_pattern.epsilon From 76e77465644a7bc40e55d93b2004651992dc42c1 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 12:05:38 +0100 Subject: [PATCH 13/20] rename --- src/general/semidiscretization.jl | 2 +- src/particle_refinement/refinement.jl | 98 ++++++++++--------- src/particle_refinement/refinement_pattern.jl | 9 +- 3 files changed, 60 insertions(+), 49 deletions(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index fd7793455..8b27a05a0 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -69,7 +69,7 @@ function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSear # Other checks might be added here later. check_configuration(systems) - systems = create_system_childs(systems) + systems = create_child_systems(systems) ranges_v, ranges_u = ranges_vu(systems) diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 96b0ef964..9961b8f00 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -3,17 +3,19 @@ include("refinement_criteria.jl") mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} candidates :: Vector{Int} + candidates_mass :: Vector{ELTYPE} refinement_pattern :: RP refinement_criteria :: RC criteria_next_levels :: Vector{RC} - available_childs :: Int + available_children :: Int # Depends on refinement pattern, particle spacing and parameters ϵ and α. - # Should be obtained prior to simulation in `create_system_child()` - rel_position_childs::Vector{SVector{NDIMS, ELTYPE}} + # Should be obtained prior to simulation in `create_child_system()` + rel_position_children :: Vector{SVector{NDIMS, ELTYPE}} + mass_ratio :: Vector{ELTYPE} # It is essential to know the child system, which is empty at the beginning - # and will be created in `create_system_child()` + # and will be created in `create_child_system()` system_child::System # API --> parent system with `RL=0` @@ -24,7 +26,7 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} NDIMS = ndims(refinement_criteria[1]) return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], refinement_pattern, + typeof(refinement_criteria)}([], [], refinement_pattern, refinement_criteria, criteria_next_levels, 0) end @@ -36,7 +38,7 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} NDIMS = ndims(refinement_criteria[1]) return new{RL, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], refinement_pattern, + typeof(refinement_criteria)}([], [], refinement_pattern, refinement_criteria, criteria_next_levels, 0) end @@ -52,40 +54,33 @@ end # ==== Create child systems -function create_system_childs(systems) +function create_child_systems(systems) systems_ = () foreach_system(systems) do system - systems_ = (systems_..., create_system_child(system)...) + systems_ = (systems_..., create_child_system(system)...) end return (systems..., systems_...) end -create_system_child(system) = () -function create_system_child(system::FluidSystem) - create_system_child(system, system.particle_refinement) +create_child_system(system) = () +function create_child_system(system::FluidSystem) + create_child_system(system, system.particle_refinement) end -create_system_child(system, ::Nothing) = () +create_child_system(system::FluidSystem, ::Nothing) = () -function create_system_child(system::FluidSystem, +function create_child_system(system::FluidSystem, particle_refinement::ParticleRefinement{RL}) where {RL} (; criteria_next_levels, refinement_pattern, refinement_criteria) = particle_refinement - if system isa WeaklyCompressibleSPHSystem - (; density_calculator, state_equation, smoothing_kernel, - pressure_acceleration_formulation, viscosity, density_diffusion, - acceleration, correction, source_terms) = system - else - (; density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, - pressure_acceleration_formulation, acceleration, correction, source_terms) = system - end - NDIMS = ndims(system) # Distribute values according to refinement pattern smoothing_length_ = smoothing_length_child(system, refinement_pattern) - particle_refinement.rel_position_childs = refinement_pattern(system) + particle_refinement.rel_position_children = relative_position_children(system, + refinement_pattern) + particle_refinement.mass_ratio = mass_distribution(system, refinement_pattern) # Create "empty" `InitialCondition` for child system particle_spacing_ = particle_spacing_child(system, refinement_pattern) @@ -109,15 +104,21 @@ function create_system_child(system::FluidSystem, end if system isa WeaklyCompressibleSPHSystem + (; density_calculator, state_equation, smoothing_kernel, + pressure_acceleration_formulation, viscosity, density_diffusion, + acceleration, correction, source_terms) = system + system_child = WeaklyCompressibleSPHSystem(empty_ic, density_calculator, - state_equation, - smoothing_kernel, smoothing_length_; + state_equation, smoothing_kernel, + smoothing_length_; pressure_acceleration=pressure_acceleration_formulation, viscosity, density_diffusion, - acceleration, - correction, source_terms, + acceleration, correction, source_terms, particle_refinement=particle_refinement_) else + (; density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, + pressure_acceleration_formulation, acceleration, correction, source_terms) = system + alpha = nu_edac * 8 / (smoothing_length_ * sound_speed) system_child = EntropicallyDampedSPHSystem(empty_ic, smoothing_kernel, smoothing_length_, sound_speed; @@ -133,7 +134,7 @@ function create_system_child(system::FluidSystem, particle_refinement.system_child = system_child return (system_child, - create_system_child(system_child, system_child.particle_refinement)...) + create_child_system(system_child, system_child.particle_refinement)...) end # ==== Refinement @@ -157,18 +158,22 @@ check_refinement_criteria!(system, ::Nothing, v_ode, u_ode, semi) = system function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, v_ode, u_ode, semi) - (; candidates, refinement_criteria) = particle_refinement + (; candidates, candidates_mass, refinement_criteria) = particle_refinement v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - !isempty(candidates) && Base.resize!(candidates, 0) + Base.resize!(candidates, 0) + Base.resize!(candidates_mass, 0) for particle in each_moving_particle(system) for refinement_criterion in refinement_criteria if (isempty(candidates) || particle != last(candidates)) && refinement_criterion(system, particle, v, u, v_ode, u_ode, semi) push!(candidates, particle) + # Store mass of candidate, since we lose the mass of the particle + # when resizing the systems + push!(candidates_mass, system.mass[particle]) end end end @@ -239,7 +244,7 @@ end function refine_particles!(system_parent::FluidSystem, particle_refinement::ParticleRefinement, v_ode, u_ode, _v_cache, _u_cache, callback, semi) - (; candidates, system_child) = particle_refinement + (; candidates, candidates_mass, system_child) = particle_refinement if !isempty(candidates) # Old storage @@ -250,16 +255,19 @@ function refine_particles!(system_parent::FluidSystem, v_child = wrap_v(v_ode, system_child, semi) u_child = wrap_u(u_ode, system_child, semi) - particle_refinement.available_childs = length(candidates) * - nchilds(system_parent, particle_refinement) + particle_refinement.available_children = length(candidates) * + nchilds(system_parent, particle_refinement) # Loop over all refinement candidates + mass_index = 1 for particle_parent in candidates - bear_childs!(system_child, system_parent, particle_parent, particle_refinement, - v_parent, u_parent, v_child, u_child, semi) + mass_parent = candidates_mass[mass_index] + bear_childs!(system_child, system_parent, particle_parent, mass_parent, + particle_refinement, v_parent, u_parent, v_child, u_child, semi) - particle_refinement.available_childs -= nchilds(system_parent, - particle_refinement) + particle_refinement.available_children -= nchilds(system_parent, + particle_refinement) + mass_index += 1 end end end @@ -269,9 +277,9 @@ end # # Reducing the dof by using a fixed regular refinement pattern # (given: position and number of child particles) -function bear_childs!(system_child, system_parent, particle_parent, particle_refinement, - v_parent, u_parent, v_child, u_child, semi) - (; rel_position_childs, available_childs, refinement_pattern) = particle_refinement +function bear_childs!(system_child, system_parent, particle_parent, mass_parent, + particle_refinement, v_parent, u_parent, v_child, u_child, semi) + (; rel_position_children, available_children, mass_ratio) = particle_refinement nhs = get_neighborhood_search(system_parent, system_parent, semi) parent_coords = current_coords(u_parent, system_parent, particle_parent) @@ -279,16 +287,12 @@ function bear_childs!(system_child, system_parent, particle_parent, particle_ref # Loop over all child particles of parent particle # The number of child particles depends on the refinement pattern for particle_child in child_set(system_parent, particle_refinement) - absolute_index = particle_child + nparticles(system_child) - available_childs + absolute_index = particle_child + nparticles(system_child) - available_children - # TODO: Handle different masses. Problem: `particle_parent` does not have an - # mass-entry anymore, since `system_parent` was resized. - mass_ = system_parent.mass[1] - system_child.mass[absolute_index] = mass_child(system_parent, mass_, - refinement_pattern) + system_child.mass[absolute_index] = mass_parent * mass_ratio[particle_child] # spread child positions according to the refinement pattern - child_coords = parent_coords + rel_position_childs[particle_child] + child_coords = parent_coords + rel_position_children[particle_child] for dim in 1:ndims(system_child) u_child[dim, absolute_index] = child_coords[dim] end diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl index d65d293df..00c48a9fb 100644 --- a/src/particle_refinement/refinement_pattern.jl +++ b/src/particle_refinement/refinement_pattern.jl @@ -7,7 +7,8 @@ struct CubicSplitting{ELTYPE} end end -function (refinement_pattern::CubicSplitting)(system::System{2}) +function relative_position_children(system::System{2}, + refinement_pattern::CubicSplitting) (; initial_condition) = system (; particle_spacing) = initial_condition (; epsilon) = refinement_pattern @@ -25,6 +26,12 @@ function (refinement_pattern::CubicSplitting)(system::System{2}) return reinterpret(reshape, SVector{2, typeof(epsilon)}, relative_position) end +function mass_distribution(system::System{2}, refinement_pattern::CubicSplitting) + lambda = 1 / nchilds(system, refinement_pattern) + + return lambda .* ones(nchilds(system, refinement_pattern)) +end + # TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? # See https://www.sciencedirect.com/science/article/pii/S0020740319317023 @inline nchilds(system, rp::CubicSplitting) = 2^ndims(system) From b9a1b917899c0437391129ed6d8cdc71b265d9fb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 15:35:07 +0100 Subject: [PATCH 14/20] add `copy_system` --- src/particle_refinement/refinement.jl | 28 ++----------------- src/schemes/fluid/fluid.jl | 40 +++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 25 deletions(-) diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 9961b8f00..94ce79768 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -53,7 +53,6 @@ end @inline nchilds(system, pr::ParticleRefinement) = nchilds(system, pr.refinement_pattern) # ==== Create child systems - function create_child_systems(systems) systems_ = () foreach_system(systems) do system @@ -103,30 +102,9 @@ function create_child_system(system::FluidSystem, criteria_next_levels[(level + 1):end]) end - if system isa WeaklyCompressibleSPHSystem - (; density_calculator, state_equation, smoothing_kernel, - pressure_acceleration_formulation, viscosity, density_diffusion, - acceleration, correction, source_terms) = system - - system_child = WeaklyCompressibleSPHSystem(empty_ic, density_calculator, - state_equation, smoothing_kernel, - smoothing_length_; - pressure_acceleration=pressure_acceleration_formulation, - viscosity, density_diffusion, - acceleration, correction, source_terms, - particle_refinement=particle_refinement_) - else - (; density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, - pressure_acceleration_formulation, acceleration, correction, source_terms) = system - - alpha = nu_edac * 8 / (smoothing_length_ * sound_speed) - system_child = EntropicallyDampedSPHSystem(empty_ic, smoothing_kernel, - smoothing_length_, sound_speed; - pressure_acceleration=pressure_acceleration_formulation, - density_calculator, alpha, viscosity, - acceleration, source_terms, - particle_refinement=particle_refinement_) - end + system_child = copy_system(system; initial_condition=empty_ic, + smoothing_length=smoothing_length_, + particle_refinement=particle_refinement_) # Empty mass vector leads to `nparticles(system_child) = 0` resize!(system_child.mass, 0) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 831fc372d..fff986927 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -101,3 +101,43 @@ end pressure) v[end, particle] = pressure end + +function copy_system(system::WeaklyCompressibleSPHSystem; + initial_condition=system.initial_condition, + density_calculator=system.density_calculator, + state_equation=system.state_equation, + smoothing_kernel=system.smoothing_kernel, + smoothing_length=system.smoothing_length, + pressure_acceleration=system.pressure_acceleration_formulation, + viscosity=system.viscosity, + density_diffusion=system.density_diffusion, + acceleration=system.acceleration, + particle_refinement=system.particle_refinement, + correction=system.correction, + source_terms=system.source_terms) + return WeaklyCompressibleSPHSystem(initial_condition, + density_calculator, state_equation, + smoothing_kernel, smoothing_length; + pressure_acceleration, viscosity, density_diffusion, + acceleration, particle_refinement, correction, + source_terms) +end + +function copy_system(system::EntropicallyDampedSPHSystem; + initial_condition=system.initial_condition, + smoothing_kernel=system.smoothing_kernel, + smoothing_length=system.smoothing_length, + sound_speed=system.sound_speed, + pressure_acceleration=system.pressure_acceleration_formulation, + density_calculator=system.density_calculator, + alpha=0.5, + viscosity=system.viscosity, + acceleration=system.acceleration, + particle_refinement=system.particle_refinement, + source_terms=system.source_terms) + return EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, sound_speed; + pressure_acceleration, density_calculator, alpha, + viscosity, + acceleration, particle_refinement, source_terms) +end From b2b0a22f67cd44671b0574949537787df2c0066c Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 18:50:42 +0100 Subject: [PATCH 15/20] add more refinement pattern --- src/TrixiParticles.jl | 3 +- src/particle_refinement/refinement.jl | 5 +- src/particle_refinement/refinement_pattern.jl | 138 ++++++++++++++---- 3 files changed, 118 insertions(+), 28 deletions(-) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index cd415c2b2..6725414f9 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -72,6 +72,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk -export ParticleRefinement, RefinementZone, CubicSplitting +export ParticleRefinement, RefinementZone +export CubicSplitting, TriangularSplitting, HexagonalSplitting end # module diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 94ce79768..901d6d226 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -71,18 +71,19 @@ create_child_system(system::FluidSystem, ::Nothing) = () function create_child_system(system::FluidSystem, particle_refinement::ParticleRefinement{RL}) where {RL} + (; particle_spacing) = system.initial_condition (; criteria_next_levels, refinement_pattern, refinement_criteria) = particle_refinement NDIMS = ndims(system) # Distribute values according to refinement pattern - smoothing_length_ = smoothing_length_child(system, refinement_pattern) + smoothing_length_ = refinement_pattern.smoothing_ratio * system.smoothing_length particle_refinement.rel_position_children = relative_position_children(system, refinement_pattern) particle_refinement.mass_ratio = mass_distribution(system, refinement_pattern) # Create "empty" `InitialCondition` for child system - particle_spacing_ = particle_spacing_child(system, refinement_pattern) + particle_spacing_ = particle_spacing * refinement_pattern.separation_parameter #TODO: Is this really the new particle spacing? coordinates_ = zeros(NDIMS, 2) velocity_ = similar(coordinates_) density_ = system.initial_condition.density[1] diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl index 00c48a9fb..5946874cf 100644 --- a/src/particle_refinement/refinement_pattern.jl +++ b/src/particle_refinement/refinement_pattern.jl @@ -1,45 +1,133 @@ +function mass_distribution(system, refinement_pattern) + if refinement_pattern.center_particle + # solve minimisation problem + + # For `HexagonalSplitting`, `separation_parameter=0.4` and `smoothing_ratio=0.9`: + lambdas = fill(0.1369, (6,)) + push!(lambdas, 0.1787) # central particle + + if isapprox(sum(lambdas), 1.0, rtol=1e-4) + return SVector(lambdas...) + end + + error("no mass conservation") + else + lambda = 1 / nchilds(system, refinement_pattern) + + return fill(lambda, SVector{nchilds(system, refinement_pattern), eltype(system)}) + end +end + struct CubicSplitting{ELTYPE} - epsilon :: ELTYPE - alpha :: ELTYPE + separation_parameter :: ELTYPE + smoothing_ratio :: ELTYPE + center_particle :: Bool - function CubicSplitting(; epsilon=0.5, alpha=0.5) - new{typeof(epsilon)}(epsilon, alpha) + function CubicSplitting(; separation_parameter=0.5, smoothing_ratio=0.5, + center_particle=true) + new{typeof(separation_parameter)}(separation_parameter, smoothing_ratio, + center_particle) end end +@inline nchilds(system, refinement_pattern::CubicSplitting) = (2^ndims(system) + + refinement_pattern.center_particle) + function relative_position_children(system::System{2}, refinement_pattern::CubicSplitting) - (; initial_condition) = system - (; particle_spacing) = initial_condition - (; epsilon) = refinement_pattern + (; particle_spacing) = system.initial_condition + (; separation_parameter) = refinement_pattern - direction_1 = normalize([1.0, 1.0]) - direction_2 = normalize([1.0, -1.0]) + direction_1 = [1 / sqrt(2), 1 / sqrt(2)] + direction_2 = [1 / sqrt(2), -1 / sqrt(2)] direction_3 = -direction_1 direction_4 = -direction_2 - relative_position = hcat(particle_spacing * epsilon * direction_1, - particle_spacing * epsilon * direction_2, - particle_spacing * epsilon * direction_3, - particle_spacing * epsilon * direction_4) + # Is it `particle_spacing * separation_parameter` or `smoothing_length * separation_parameter`? + relative_position = hcat(particle_spacing * separation_parameter * direction_1, + particle_spacing * separation_parameter * direction_2, + particle_spacing * separation_parameter * direction_3, + particle_spacing * separation_parameter * direction_4) - return reinterpret(reshape, SVector{2, typeof(epsilon)}, relative_position) + if refinement_pattern.center_particle + relative_position = hcat(relative_position, [0.0, 0.0]) + end + + return reinterpret(reshape, SVector{2, eltype(system)}, relative_position) end -function mass_distribution(system::System{2}, refinement_pattern::CubicSplitting) - lambda = 1 / nchilds(system, refinement_pattern) +struct TriangularSplitting{ELTYPE} + separation_parameter :: ELTYPE + smoothing_ratio :: ELTYPE + center_particle :: Bool - return lambda .* ones(nchilds(system, refinement_pattern)) + function TriangularSplitting(; separation_parameter=0.5, smoothing_ratio=0.5, + center_particle=false) + new{typeof(separation_parameter)}(separation_parameter, smoothing_ratio, + center_particle) + end end -# TODO: Clarify refinement pattern. Cubic splitting? Triangular or hexagonal? -# See https://www.sciencedirect.com/science/article/pii/S0020740319317023 -@inline nchilds(system, rp::CubicSplitting) = 2^ndims(system) +@inline nchilds(system::System{2}, refinement_pattern::TriangularSplitting) = 3 + + refinement_pattern.center_particle -@inline mass_child(system, mass, rp::CubicSplitting) = mass / nchilds(system, rp) +function relative_position_children(system::System{2}, + refinement_pattern::TriangularSplitting) + (; particle_spacing) = system.initial_condition + (; separation_parameter) = refinement_pattern -@inline smoothing_length_child(system, refinement_pattern) = refinement_pattern.alpha * - system.smoothing_length + direction_1 = [0.0, 1.0] + direction_2 = [-sqrt(3) / 2, -0.5] + direction_3 = [sqrt(3) / 2, -0.5] + + relative_position = hcat(particle_spacing * separation_parameter * direction_1, + particle_spacing * separation_parameter * direction_2, + particle_spacing * separation_parameter * direction_3) + + if refinement_pattern.center_particle + relative_position = hcat(relative_position, [0.0, 0.0]) + end -@inline particle_spacing_child(system, refinement_pattern) = system.initial_condition.particle_spacing * - refinement_pattern.epsilon + return reinterpret(reshape, SVector{2, eltype(system)}, relative_position) +end + +struct HexagonalSplitting{ELTYPE} + separation_parameter :: ELTYPE + smoothing_ratio :: ELTYPE + center_particle :: Bool + + function HexagonalSplitting(; separation_parameter=0.5, smoothing_ratio=0.5, + center_particle=true) + new{typeof(separation_parameter)}(separation_parameter, smoothing_ratio, + center_particle) + end +end + +@inline nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) = 6 + + refinement_pattern.center_particle + +function relative_position_children(system::System{2}, + refinement_pattern::HexagonalSplitting) + (; particle_spacing) = system.initial_condition + (; separation_parameter) = refinement_pattern + + direction_1 = [1.0, 0.0] + direction_2 = [-1.0, 0.0] + direction_3 = [0.5, sqrt(3) / 2] + direction_4 = [0.5, -sqrt(3) / 2] + direction_5 = [-0.5, sqrt(3) / 2] + direction_6 = [-0.5, -sqrt(3) / 2] + + relative_position = hcat(particle_spacing * separation_parameter * direction_1, + particle_spacing * separation_parameter * direction_2, + particle_spacing * separation_parameter * direction_3, + particle_spacing * separation_parameter * direction_4, + particle_spacing * separation_parameter * direction_5, + particle_spacing * separation_parameter * direction_6) + + if refinement_pattern.center_particle + relative_position = hcat(relative_position, [0.0, 0.0]) + end + + return reinterpret(reshape, SVector{2, eltype(system)}, relative_position) +end From 74c70c9c1b6f4959bd257d5b50fb72a5a52051b6 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 20 Mar 2024 23:19:56 +0100 Subject: [PATCH 16/20] revise refinement criteria --- src/callbacks/refinement_callback.jl | 2 +- src/particle_refinement/refinement.jl | 16 ++--- .../refinement_criteria.jl | 70 +++++++++++-------- 3 files changed, 49 insertions(+), 39 deletions(-) diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl index e2e1edfab..d2b829d86 100644 --- a/src/callbacks/refinement_callback.jl +++ b/src/callbacks/refinement_callback.jl @@ -76,7 +76,7 @@ function (refinement_callback::ParticleRefinementCallback)(integrator) v_tmp .= v_ode u_tmp .= u_ode - refinement!(v_ode, u_ode, v_tmp, u_tmp, semi, refinement_callback) + refinement!(v_ode, u_ode, v_tmp, u_tmp, semi, refinement_callback, t) # Resize neighborhood search foreach_system(semi) do system diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 901d6d226..4f0b24434 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -117,9 +117,9 @@ function create_child_system(system::FluidSystem, end # ==== Refinement -function refinement!(v_ode, u_ode, _v_cache, _u_cache, semi, callback) +function refinement!(v_ode, u_ode, _v_cache, _u_cache, semi, callback, t) foreach_system(semi) do system - check_refinement_criteria!(system, v_ode, u_ode, semi) + check_refinement_criteria!(system, v_ode, u_ode, semi, t) end resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) @@ -127,16 +127,16 @@ function refinement!(v_ode, u_ode, _v_cache, _u_cache, semi, callback) refine_particles!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) end -check_refinement_criteria!(system, v_ode, u_ode, semi) = system +check_refinement_criteria!(system, v_ode, u_ode, semi, t) = system -function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) - check_refinement_criteria!(system, system.particle_refinement, v_ode, u_ode, semi) +function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi, t) + check_refinement_criteria!(system, system.particle_refinement, v_ode, u_ode, semi, t) end -check_refinement_criteria!(system, ::Nothing, v_ode, u_ode, semi) = system +check_refinement_criteria!(system, ::Nothing, v_ode, u_ode, semi, t) = system function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, - v_ode, u_ode, semi) + v_ode, u_ode, semi, t) (; candidates, candidates_mass, refinement_criteria) = particle_refinement v = wrap_v(v_ode, system, semi) @@ -148,7 +148,7 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem for particle in each_moving_particle(system) for refinement_criterion in refinement_criteria if (isempty(candidates) || particle != last(candidates)) && - refinement_criterion(system, particle, v, u, v_ode, u_ode, semi) + refinement_criterion(system, particle, v, u, v_ode, u_ode, semi, t) push!(candidates, particle) # Store mass of candidate, since we lose the mass of the particle # when resizing the systems diff --git a/src/particle_refinement/refinement_criteria.jl b/src/particle_refinement/refinement_criteria.jl index 3f28635c1..392d9c495 100644 --- a/src/particle_refinement/refinement_criteria.jl +++ b/src/particle_refinement/refinement_criteria.jl @@ -6,62 +6,72 @@ abstract type RefinementCriteria{NDIMS, ELTYPE} end -struct RefinementZone{NDIMS, ELTYPE} <: RefinementCriteria{NDIMS, ELTYPE} - zone_origin :: SVector{NDIMS, ELTYPE} +struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE} + zone_origin :: ZO spanning_set :: Vector{SVector} - function RefinementZone(plane_points, zone_width) - NDIMS = length(plane_points) - ELTYPE = typeof(zone_width) + function RefinementZone(edge_lengths; + zone_origin=ntuple(_ -> 0.0, length(edge_lengths)), + rotation=nothing) # TODO + NDIMS = length(edge_lengths) + ELTYPE = eltype(edge_lengths) - # Vectors spanning the zone. - spanning_set = spanning_vectors(plane_points, zone_width) - - spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + if zone_origin isa Function + zone_origin_function = zone_origin + else + zone_origin_function = (v, u, v_ode, u_ode, t, system, semi) -> SVector(zone_origin...) + end - zone_origin = SVector(plane_points[1]...) + # Vectors spanning the zone. + spanning_set = spanning_vectors(edge_lengths, rotation) - return new{NDIMS, ELTYPE}(zone_origin, spanning_set_) + return new{NDIMS, ELTYPE, + typeof(zone_origin_function)}(zone_origin_function, spanning_set) end end @inline Base.ndims(::RefinementCriteria{NDIMS}) where {NDIMS} = NDIMS @inline Base.eltype(::RefinementCriteria{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE -function spanning_vectors(plane_points, zone_width) +function spanning_vectors(edge_lengths, rotation) # Convert to tuple - return spanning_vectors(tuple(plane_points...), zone_width) + return spanning_vectors(tuple(edge_lengths...), rotation) end -function spanning_vectors(plane_points::NTuple{2}, zone_width) - plane_size = plane_points[2] - plane_points[1] +function spanning_vectors(edge_lengths::NTuple{2}, rotation) + ELTYPE = eltype(edge_lengths) + vec1 = [edge_lengths[1], 0.0] + vec2 = [0.0, edge_lengths[2]] - # Calculate normal vector of plane - b = Vector(normalize([-plane_size[2]; plane_size[1]]) * zone_width) + if !isnothing(rotation) + # rotate vecs + return reinterpret(reshape, SVector{2, ELTYPE}, hcat(vec1, vec2)) + end - return hcat(b, plane_size) + return reinterpret(reshape, SVector{2, ELTYPE}, hcat(vec1, vec2)) end -function spanning_vectors(plane_points::NTuple{3}, zone_width) - # Vectors spanning the plane - edge1 = plane_points[2] - plane_points[1] - edge2 = plane_points[3] - plane_points[1] +function spanning_vectors(plane_points::NTuple{3}, rotation) + ELTYPE = eltype(edge_lengths) - if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) - throw(ArgumentError("the provided points do not span a rectangular plane")) - end + vec1 = [edge_lengths[1], 0.0, 0.0] + vec2 = [0.0, edge_lengths[2], 0.0] + vec3 = [0.0, 0.0, edge_lengths[3]] - # Calculate normal vector of plane - c = Vector(normalize(cross(edge2, edge1)) * zone_width) + if !isnothing(rotation) + # rotate vecs + return reinterpret(reshape, SVector{3, ELTYPE}, hcat(vec1, vec2, vec3)) + end - return hcat(c, edge1, edge2) + return reinterpret(reshape, SVector{3, ELTYPE}, hcat(vec1, vec2, vec3)) end @inline function (refinement_criterion::RefinementZone)(system, particle, - v, u, v_ode, u_ode, semi) + v, u, v_ode, u_ode, semi, t) (; zone_origin, spanning_set) = refinement_criterion - particle_position = current_coords(u, system, particle) - zone_origin + particle_position = current_coords(u, system, particle) - + zone_origin(v, u, v_ode, u_ode, t, system, semi) for dim in 1:ndims(system) span_dim = spanning_set[dim] From 3dbe9f03de7c3ed04494dcb16dff6c0ccbba7de8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 21 Mar 2024 09:53:18 +0100 Subject: [PATCH 17/20] add comments --- src/particle_refinement/refinement.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 4f0b24434..ceceea587 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -15,7 +15,7 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} mass_ratio :: Vector{ELTYPE} # It is essential to know the child system, which is empty at the beginning - # and will be created in `create_child_system()` + # and will be created in `create_child_system()` at the beginning of the simulation system_child::System # API --> parent system with `RL=0` @@ -290,6 +290,12 @@ function bear_childs!(system_child, system_parent, particle_parent, mass_parent, distance2 = dot(pos_diff, pos_diff) + # TODO: Check the following statement + # + # For the Navier–Stokes equations Feldman and Bonet showed that + # the only way to conserve both total momentum and energy is to define + # the velocities of the daughter particles `v_child` equal to the + # the velocity of the original parent particle therefore: `v_child = v_parent` if distance2 <= nhs.search_radius^2 distance = sqrt(distance2) kernel_weight = smoothing_kernel(system_parent, distance) From fb9f0eaac47aef92a397059a1c65cc20ded0efb3 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 25 Mar 2024 18:17:28 +0100 Subject: [PATCH 18/20] separation parameter: smoothing length --- src/particle_refinement/refinement.jl | 41 +++++----- .../refinement_criteria.jl | 74 +++++++++---------- src/particle_refinement/refinement_pattern.jl | 34 ++++----- .../fluid/weakly_compressible_sph/system.jl | 4 + 4 files changed, 79 insertions(+), 74 deletions(-) diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index ceceea587..c26436124 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -1,12 +1,12 @@ include("refinement_pattern.jl") include("refinement_criteria.jl") -mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} +mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC, CNL} candidates :: Vector{Int} candidates_mass :: Vector{ELTYPE} refinement_pattern :: RP refinement_criteria :: RC - criteria_next_levels :: Vector{RC} + criteria_next_levels :: CNL available_children :: Int # Depends on refinement pattern, particle spacing and parameters ϵ and α. @@ -22,29 +22,36 @@ mutable struct ParticleRefinement{RL, NDIMS, ELTYPE, RP, RC} function ParticleRefinement(refinement_criteria...; refinement_pattern=CubicSplitting(), criteria_next_levels=[]) - ELTYPE = eltype(refinement_criteria[1]) - NDIMS = ndims(refinement_criteria[1]) + ELTYPE = eltype(first(refinement_criteria)) + NDIMS = ndims(first(refinement_criteria)) return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], [], refinement_pattern, - refinement_criteria, - criteria_next_levels, 0) + typeof(refinement_criteria), + typeof(criteria_next_levels)}([], [], refinement_pattern, + refinement_criteria, + criteria_next_levels, 0) end # Internal constructor for multiple refinement levels - function ParticleRefinement{RL}(refinement_criteria::Tuple, - refinement_pattern, criteria_next_levels) where {RL} - ELTYPE = eltype(refinement_criteria[1]) - NDIMS = ndims(refinement_criteria[1]) + function ParticleRefinement{RL}(refinement_criteria, refinement_pattern, + criteria_next_levels) where {RL} + if refinement_criteria isa Tuple + ELTYPE = eltype(first(refinement_criteria)) + NDIMS = ndims(first(refinement_criteria)) + else + ELTYPE = eltype(refinement_criteria) + NDIMS = ndims(refinement_criteria) + end return new{RL, NDIMS, ELTYPE, typeof(refinement_pattern), - typeof(refinement_criteria)}([], [], refinement_pattern, - refinement_criteria, - criteria_next_levels, 0) + typeof(refinement_criteria), + typeof(criteria_next_levels)}([], [], refinement_pattern, + refinement_criteria, + criteria_next_levels, 0) end end -@inline Base.ndims(::ParticleRefinement{RL, NDIMS}) where {RL, NDIMS} = ndims +@inline refinement_level(::ParticleRefinement{RL}) where {RL} = RL @inline child_set(system, particle_refinement) = Base.OneTo(nchilds(system, particle_refinement)) @@ -71,7 +78,7 @@ create_child_system(system::FluidSystem, ::Nothing) = () function create_child_system(system::FluidSystem, particle_refinement::ParticleRefinement{RL}) where {RL} - (; particle_spacing) = system.initial_condition + (; smoothing_length) = system (; criteria_next_levels, refinement_pattern, refinement_criteria) = particle_refinement NDIMS = ndims(system) @@ -83,7 +90,7 @@ function create_child_system(system::FluidSystem, particle_refinement.mass_ratio = mass_distribution(system, refinement_pattern) # Create "empty" `InitialCondition` for child system - particle_spacing_ = particle_spacing * refinement_pattern.separation_parameter #TODO: Is this really the new particle spacing? + particle_spacing_ = smoothing_length * refinement_pattern.separation_parameter coordinates_ = zeros(NDIMS, 2) velocity_ = similar(coordinates_) density_ = system.initial_condition.density[1] diff --git a/src/particle_refinement/refinement_criteria.jl b/src/particle_refinement/refinement_criteria.jl index 392d9c495..54ea2311c 100644 --- a/src/particle_refinement/refinement_criteria.jl +++ b/src/particle_refinement/refinement_criteria.jl @@ -10,11 +10,25 @@ struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE} zone_origin :: ZO spanning_set :: Vector{SVector} - function RefinementZone(edge_lengths; - zone_origin=ntuple(_ -> 0.0, length(edge_lengths)), - rotation=nothing) # TODO - NDIMS = length(edge_lengths) - ELTYPE = eltype(edge_lengths) + function RefinementZone(; edge_length_x=0.0, edge_length_y=0.0, edge_length_z=nothing, + zone_origin, rotation=nothing) # TODO + if isnothing(edge_length_z) + NDIMS = 2 + elseif edge_length_z < eps() + throw(ArgumentError("`edge_length_z` must be either `nothing` for a 2D problem or greater than zero for a 3D problem")) + else + NDIMS = 3 + end + + ELTYPE = eltype(edge_length_x) + + if edge_length_x * edge_length_y < eps() + throw(ArgumentError("edge lengths must be greater than zero")) + end + + if length(zone_origin) != NDIMS + throw(ArgumentError("`zone_origin` must be a `Vector` of size $NDIMS for a $NDIMS-D Problem")) + end if zone_origin isa Function zone_origin_function = zone_origin @@ -22,51 +36,31 @@ struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE} zone_origin_function = (v, u, v_ode, u_ode, t, system, semi) -> SVector(zone_origin...) end + edge_lengths = if NDIMS == 2 + [edge_length_x, edge_length_y] + else + [edge_length_x, edge_length_y, edge_length_z] + end + # Vectors spanning the zone. - spanning_set = spanning_vectors(edge_lengths, rotation) + spanning_set_ = I(NDIMS) .* edge_lengths' + + spanning_set = if !isnothing(rotation) + # rotate vecs + reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set_) + else + reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set_) + end return new{NDIMS, ELTYPE, typeof(zone_origin_function)}(zone_origin_function, spanning_set) end + end @inline Base.ndims(::RefinementCriteria{NDIMS}) where {NDIMS} = NDIMS @inline Base.eltype(::RefinementCriteria{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE -function spanning_vectors(edge_lengths, rotation) - - # Convert to tuple - return spanning_vectors(tuple(edge_lengths...), rotation) -end - -function spanning_vectors(edge_lengths::NTuple{2}, rotation) - ELTYPE = eltype(edge_lengths) - vec1 = [edge_lengths[1], 0.0] - vec2 = [0.0, edge_lengths[2]] - - if !isnothing(rotation) - # rotate vecs - return reinterpret(reshape, SVector{2, ELTYPE}, hcat(vec1, vec2)) - end - - return reinterpret(reshape, SVector{2, ELTYPE}, hcat(vec1, vec2)) -end - -function spanning_vectors(plane_points::NTuple{3}, rotation) - ELTYPE = eltype(edge_lengths) - - vec1 = [edge_lengths[1], 0.0, 0.0] - vec2 = [0.0, edge_lengths[2], 0.0] - vec3 = [0.0, 0.0, edge_lengths[3]] - - if !isnothing(rotation) - # rotate vecs - return reinterpret(reshape, SVector{3, ELTYPE}, hcat(vec1, vec2, vec3)) - end - - return reinterpret(reshape, SVector{3, ELTYPE}, hcat(vec1, vec2, vec3)) -end - @inline function (refinement_criterion::RefinementZone)(system, particle, v, u, v_ode, u_ode, semi, t) (; zone_origin, spanning_set) = refinement_criterion diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl index 5946874cf..f7024bc3d 100644 --- a/src/particle_refinement/refinement_pattern.jl +++ b/src/particle_refinement/refinement_pattern.jl @@ -35,7 +35,7 @@ end function relative_position_children(system::System{2}, refinement_pattern::CubicSplitting) - (; particle_spacing) = system.initial_condition + (; smoothing_length) = system (; separation_parameter) = refinement_pattern direction_1 = [1 / sqrt(2), 1 / sqrt(2)] @@ -44,10 +44,10 @@ function relative_position_children(system::System{2}, direction_4 = -direction_2 # Is it `particle_spacing * separation_parameter` or `smoothing_length * separation_parameter`? - relative_position = hcat(particle_spacing * separation_parameter * direction_1, - particle_spacing * separation_parameter * direction_2, - particle_spacing * separation_parameter * direction_3, - particle_spacing * separation_parameter * direction_4) + relative_position = hcat(smoothing_length * separation_parameter * direction_1, + smoothing_length * separation_parameter * direction_2, + smoothing_length * separation_parameter * direction_3, + smoothing_length * separation_parameter * direction_4) if refinement_pattern.center_particle relative_position = hcat(relative_position, [0.0, 0.0]) @@ -62,7 +62,7 @@ struct TriangularSplitting{ELTYPE} center_particle :: Bool function TriangularSplitting(; separation_parameter=0.5, smoothing_ratio=0.5, - center_particle=false) + center_particle=true) new{typeof(separation_parameter)}(separation_parameter, smoothing_ratio, center_particle) end @@ -73,16 +73,16 @@ end function relative_position_children(system::System{2}, refinement_pattern::TriangularSplitting) - (; particle_spacing) = system.initial_condition + (; smoothing_length) = system (; separation_parameter) = refinement_pattern direction_1 = [0.0, 1.0] direction_2 = [-sqrt(3) / 2, -0.5] direction_3 = [sqrt(3) / 2, -0.5] - relative_position = hcat(particle_spacing * separation_parameter * direction_1, - particle_spacing * separation_parameter * direction_2, - particle_spacing * separation_parameter * direction_3) + relative_position = hcat(smoothing_length * separation_parameter * direction_1, + smoothing_length * separation_parameter * direction_2, + smoothing_length * separation_parameter * direction_3) if refinement_pattern.center_particle relative_position = hcat(relative_position, [0.0, 0.0]) @@ -108,7 +108,7 @@ end function relative_position_children(system::System{2}, refinement_pattern::HexagonalSplitting) - (; particle_spacing) = system.initial_condition + (; smoothing_length) = system (; separation_parameter) = refinement_pattern direction_1 = [1.0, 0.0] @@ -118,12 +118,12 @@ function relative_position_children(system::System{2}, direction_5 = [-0.5, sqrt(3) / 2] direction_6 = [-0.5, -sqrt(3) / 2] - relative_position = hcat(particle_spacing * separation_parameter * direction_1, - particle_spacing * separation_parameter * direction_2, - particle_spacing * separation_parameter * direction_3, - particle_spacing * separation_parameter * direction_4, - particle_spacing * separation_parameter * direction_5, - particle_spacing * separation_parameter * direction_6) + relative_position = hcat(smoothing_length * separation_parameter * direction_1, + smoothing_length * separation_parameter * direction_2, + smoothing_length * separation_parameter * direction_3, + smoothing_length * separation_parameter * direction_4, + smoothing_length * separation_parameter * direction_5, + smoothing_length * separation_parameter * direction_6) if refinement_pattern.center_particle relative_position = hcat(relative_position, [0.0, 0.0]) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index bec18d05e..28c07fd57 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -147,6 +147,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.density_diffusion) print(io, ", ", system.acceleration) print(io, ", ", system.source_terms) + print(io, ", ", system.particle_refinement) print(io, ") with ", nparticles(system), " particles") end @@ -168,6 +169,9 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "density diffusion", system.density_diffusion) summary_line(io, "acceleration", system.acceleration) summary_line(io, "source terms", system.source_terms |> typeof |> nameof) + if !isnothing(system.particle_refinement) + summary_line(io, "refinement level", refinement_level(system.particle_refinement)) + end summary_footer(io) end end From 5f19582fe96a1f9ddabd2a86c27b28b1dbb648a9 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 26 Mar 2024 15:15:59 +0100 Subject: [PATCH 19/20] move `resize` stuff in separate file --- src/particle_refinement/refinement.jl | 153 +------------------------ src/particle_refinement/resize.jl | 154 ++++++++++++++++++++++++++ 2 files changed, 156 insertions(+), 151 deletions(-) create mode 100644 src/particle_refinement/resize.jl diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index c26436124..0fd687b51 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -59,6 +59,8 @@ end @inline nchilds(system, ::Nothing) = 0 @inline nchilds(system, pr::ParticleRefinement) = nchilds(system, pr.refinement_pattern) +include("resize.jl") + # ==== Create child systems function create_child_systems(systems) systems_ = () @@ -167,45 +169,6 @@ function check_refinement_criteria!(system, particle_refinement::ParticleRefinem return system end -function resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) - # Get non-`resize!`d ranges - callback.ranges_v_cache, callback.ranges_u_cache = ranges_vu(semi.systems) - callback.eachparticle_cache = Tuple(get_iterator(system) for system in semi.systems) - callback.nparticles_cache = Tuple(n_moving_particles(system) for system in semi.systems) - - # Resize all systems - foreach_system(semi) do system - resize_system!(system) - end - - # Set `resize!`d ranges - ranges_v_tmp, ranges_u_tmp = ranges_vu(semi.systems) - for i in 1:length(semi.systems) - semi.ranges_v[i][1] = ranges_v_tmp[i][1] - semi.ranges_u[i][1] = ranges_u_tmp[i][1] - end - - sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in semi.systems) - sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in semi.systems) - - # Resize integrated values - resize!(v_ode, sum(sizes_v)) - resize!(u_ode, sum(sizes_u)) - - # Preserve non-changing values - foreach_system(semi) do system - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) - _v = _wrap_v(_v_cache, system, semi, callback) - _u = _wrap_u(_u_cache, system, semi, callback) - - copy_values_v!(v, _v, system, semi, callback) - copy_values_u!(u, _u, system, semi, callback) - end - - return callback -end - function refine_particles!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) # Refine particles in all systems @@ -338,118 +301,6 @@ function bear_childs!(system_child, system_parent, particle_parent, mass_parent, return system_child end -@inline resize_system!(system) = system -@inline resize_system!(system::FluidSystem) = resize_system!(system, - system.particle_refinement) - -@inline resize_system!(system::FluidSystem, ::Nothing) = system - -@inline function resize_system!(system::FluidSystem, - particle_refinement::ParticleRefinement) - (; candidates, system_child) = particle_refinement - - if !isempty(candidates) - n_new_child = length(candidates) * nchilds(system, particle_refinement) - capacity_parent = nparticles(system) - length(candidates) - capacity_child = nparticles(system_child) + n_new_child - - capacity_parent <= 0 && error("`RefinementCriteria` affects all particles") - - # Resize child system (extending) - resize_system!(system_child, capacity_child) - - # Resize parent system (reducing) - resize_system!(system, capacity_parent) - end -end - -function resize_system!(system::WeaklyCompressibleSPHSystem, capacity::Int) - (; mass, pressure, cache, density_calculator) = system - - resize!(mass, capacity) - resize!(pressure, capacity) - resize_cache!(cache, capacity, density_calculator) -end - -function resize_system!(system::EntropicallyDampedSPHSystem, capacity::Int) - (; mass, cache, density_calculator) = system - - resize!(mass, capacity) - resize_cache!(cache, capacity, density_calculator) -end - -resize_cache!(cache, capacity::Int, ::SummationDensity) = resize!(cache.density, capacity) -resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache - -@inline function _wrap_u(_u_cache, system, semi, callback) - (; ranges_u_cache, nparticles_cache) = callback - - range = ranges_u_cache[system_indices(system, semi)][1] - n_particles = nparticles_cache[system_indices(system, semi)] - - @boundscheck @assert length(range) == u_nvariables(system) * n_particles - - # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(_u_cache), 2}, pointer(view(_u_cache, range)), - # (u_nvariables(system), n_particles)) - return PtrArray(pointer(view(_u_cache, range)), - (StaticInt(u_nvariables(system)), n_particles)) -end - -@inline function _wrap_v(_v_cache, system, semi, callback) - (; ranges_v_cache, nparticles_cache) = callback - - range = ranges_v_cache[system_indices(system, semi)][1] - n_particles = nparticles_cache[system_indices(system, semi)] - - @boundscheck @assert length(range) == v_nvariables(system) * n_particles - - # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(_v_cache), 2}, pointer(view(_v_cache, range)), - # (v_nvariables(system), n_particles)) - return PtrArray(pointer(view(_v_cache, range)), - (StaticInt(v_nvariables(system)), n_particles)) -end - -function copy_values_v!(v_new, v_old, system, semi, callback) - (; eachparticle_cache) = callback - - # Copy only non-refined particles - new_particle_id = 1 - for particle in eachparticle_cache[system_indices(system, semi)] - for i in 1:v_nvariables(system) - v_new[i, new_particle_id] = v_old[i, particle] - end - new_particle_id += 1 - end -end - -function copy_values_u!(u_new, u_old, system, semi, callback) - (; eachparticle_cache) = callback - - # Copy only non-refined particles - new_particle_id = 1 - for particle in eachparticle_cache[system_indices(system, semi)] - for i in 1:u_nvariables(system) - u_new[i, new_particle_id] = u_old[i, particle] - end - new_particle_id += 1 - end -end - -@inline get_iterator(system) = eachparticle(system) -@inline get_iterator(system::FluidSystem) = get_iterator(system, system.particle_refinement) - -@inline get_iterator(system, ::Nothing) = eachparticle(system) - -@inline function get_iterator(system, particle_refinement::ParticleRefinement) - (; candidates) = particle_refinement - - # Filter candidates - #return Iterators.filter(i -> !(i in candidates), eachparticle(system)) - return setdiff(eachparticle(system), candidates) -end - @inline particle_pressure_parent(v, ::WeaklyCompressibleSPHSystem, particle) = 0.0 @inline particle_pressure_parent(v, system::EntropicallyDampedSPHSystem, particle) = particle_pressure(v, system, diff --git a/src/particle_refinement/resize.jl b/src/particle_refinement/resize.jl new file mode 100644 index 000000000..d7db857c7 --- /dev/null +++ b/src/particle_refinement/resize.jl @@ -0,0 +1,154 @@ +function resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) + # Get non-`resize!`d ranges + callback.ranges_v_cache, callback.ranges_u_cache = ranges_vu(semi.systems) + callback.eachparticle_cache = Tuple(get_iterator(system) for system in semi.systems) + callback.nparticles_cache = Tuple(n_moving_particles(system) for system in semi.systems) + + # Resize all systems + foreach_system(semi) do system + resize_system!(system) + end + + # Set `resize!`d ranges + ranges_v_tmp, ranges_u_tmp = ranges_vu(semi.systems) + for i in 1:length(semi.systems) + semi.ranges_v[i][1] = ranges_v_tmp[i][1] + semi.ranges_u[i][1] = ranges_u_tmp[i][1] + end + + sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in semi.systems) + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in semi.systems) + + # Resize integrated values + resize!(v_ode, sum(sizes_v)) + resize!(u_ode, sum(sizes_u)) + + # Preserve non-changing values + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + _v = _wrap_v(_v_cache, system, semi, callback) + _u = _wrap_u(_u_cache, system, semi, callback) + + copy_values_v!(v, _v, system, semi, callback) + copy_values_u!(u, _u, system, semi, callback) + end + + return callback +end + +@inline resize_system!(system) = system +@inline resize_system!(system::FluidSystem) = resize_system!(system, + system.particle_refinement) + +@inline resize_system!(system::FluidSystem, ::Nothing) = system + +@inline function resize_system!(system::FluidSystem, + particle_refinement::ParticleRefinement) + (; system_child, candidates) = particle_refinement + + candidates_refine = length(candidates) + + if candidates_refine != 0 || candidates_coarsen != 0 + n_new_child = candidates_refine * nchilds(system, particle_refinement) + capacity_parent = nparticles(system) - candidates_refine + capacity_child = nparticles(system_child) + n_new_child + + capacity_parent <= 0 && error("`RefinementCriteria` affects all particles") + + # Resize child system (extending) + resize_system!(system_child, capacity_child) + + # Resize parent system (reducing) + resize_system!(system, capacity_parent) + end +end + +function resize_system!(system::WeaklyCompressibleSPHSystem, capacity::Int) + (; mass, pressure, cache, density_calculator) = system + + resize!(mass, capacity) + resize!(pressure, capacity) + resize_cache!(cache, capacity, density_calculator) +end + +function resize_system!(system::EntropicallyDampedSPHSystem, capacity::Int) + (; mass, cache, density_calculator) = system + + resize!(mass, capacity) + resize_cache!(cache, capacity, density_calculator) +end + +resize_cache!(cache, capacity::Int, ::SummationDensity) = resize!(cache.density, capacity) +resize_cache!(cache, capacity::Int, ::ContinuityDensity) = cache + +@inline function _wrap_u(_u_cache, system, semi, callback) + (; ranges_u_cache, nparticles_cache) = callback + + range = ranges_u_cache[system_indices(system, semi)][1] + n_particles = nparticles_cache[system_indices(system, semi)] + + @boundscheck @assert length(range) == u_nvariables(system) * n_particles + + # This is a non-allocating version of: + # return unsafe_wrap(Array{eltype(_u_cache), 2}, pointer(view(_u_cache, range)), + # (u_nvariables(system), n_particles)) + return PtrArray(pointer(view(_u_cache, range)), + (StaticInt(u_nvariables(system)), n_particles)) +end + +@inline function _wrap_v(_v_cache, system, semi, callback) + (; ranges_v_cache, nparticles_cache) = callback + + range = ranges_v_cache[system_indices(system, semi)][1] + n_particles = nparticles_cache[system_indices(system, semi)] + + @boundscheck @assert length(range) == v_nvariables(system) * n_particles + + # This is a non-allocating version of: + # return unsafe_wrap(Array{eltype(_v_cache), 2}, pointer(view(_v_cache, range)), + # (v_nvariables(system), n_particles)) + return PtrArray(pointer(view(_v_cache, range)), + (StaticInt(v_nvariables(system)), n_particles)) +end + +function copy_values_v!(v_new, v_old, system, semi, callback) + (; eachparticle_cache) = callback + + # Copy only unrefined particles + new_particle_id = 1 + for particle in eachparticle_cache[system_indices(system, semi)] + for i in 1:v_nvariables(system) + v_new[i, new_particle_id] = v_old[i, particle] + end + new_particle_id += 1 + end +end + +function copy_values_u!(u_new, u_old, system, semi, callback) + (; eachparticle_cache) = callback + + # Copy only unrefined particles + new_particle_id = 1 + for particle in eachparticle_cache[system_indices(system, semi)] + for i in 1:u_nvariables(system) + u_new[i, new_particle_id] = u_old[i, particle] + end + new_particle_id += 1 + end +end + +@inline get_iterator(system) = eachparticle(system) +@inline get_iterator(system::FluidSystem) = get_iterator(system, system.particle_refinement) + +@inline get_iterator(system::FluidSystem, ::Nothing) = eachparticle(system) + +@inline function get_iterator(system::FluidSystem, + particle_refinement::ParticleRefinement) + (; candidates) = particle_refinement + + # Filter candidates + # Uncomment for benchmark + # return Iterators.filter(i -> !(i in candidates), eachparticle(system)) + return setdiff(eachparticle(system), candidates) +end From e5a890bdc2b69bfdbc3a6d3d88b15c72ea6384bb Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 26 Mar 2024 17:48:15 +0100 Subject: [PATCH 20/20] add tests --- src/particle_refinement/refinement.jl | 11 +- .../refinement_criteria.jl | 1 - src/particle_refinement/refinement_pattern.jl | 9 +- src/particle_refinement/resize.jl | 2 +- .../fluid/weakly_compressible_sph/system.jl | 3 +- test/general/semidiscretization.jl | 4 +- .../particle_refinement.jl | 2 + test/particle_refinement/refinement.jl | 223 ++++++++++++++++++ test/particle_refinement/refinement_resize.jl | 155 ++++++++++++ test/systems/wcsph_system.jl | 2 +- test/unittest.jl | 1 + 11 files changed, 399 insertions(+), 14 deletions(-) create mode 100644 test/particle_refinement/particle_refinement.jl create mode 100644 test/particle_refinement/refinement.jl create mode 100644 test/particle_refinement/refinement_resize.jl diff --git a/src/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl index 0fd687b51..5ac228b27 100644 --- a/src/particle_refinement/refinement.jl +++ b/src/particle_refinement/refinement.jl @@ -196,6 +196,8 @@ function refine_particles!(system_parent::FluidSystem, (; candidates, candidates_mass, system_child) = particle_refinement if !isempty(candidates) + nhs = get_neighborhood_search(system_parent, semi) + # Old storage v_parent = _wrap_v(_v_cache, system_parent, semi, callback) u_parent = _wrap_u(_u_cache, system_parent, semi, callback) @@ -211,8 +213,8 @@ function refine_particles!(system_parent::FluidSystem, mass_index = 1 for particle_parent in candidates mass_parent = candidates_mass[mass_index] - bear_childs!(system_child, system_parent, particle_parent, mass_parent, - particle_refinement, v_parent, u_parent, v_child, u_child, semi) + bear_children!(system_child, system_parent, particle_parent, mass_parent, nhs, + particle_refinement, v_parent, u_parent, v_child, u_child) particle_refinement.available_children -= nchilds(system_parent, particle_refinement) @@ -226,11 +228,10 @@ end # # Reducing the dof by using a fixed regular refinement pattern # (given: position and number of child particles) -function bear_childs!(system_child, system_parent, particle_parent, mass_parent, - particle_refinement, v_parent, u_parent, v_child, u_child, semi) +function bear_children!(system_child, system_parent, particle_parent, mass_parent, nhs, + particle_refinement, v_parent, u_parent, v_child, u_child) (; rel_position_children, available_children, mass_ratio) = particle_refinement - nhs = get_neighborhood_search(system_parent, system_parent, semi) parent_coords = current_coords(u_parent, system_parent, particle_parent) # Loop over all child particles of parent particle diff --git a/src/particle_refinement/refinement_criteria.jl b/src/particle_refinement/refinement_criteria.jl index 54ea2311c..b777543c6 100644 --- a/src/particle_refinement/refinement_criteria.jl +++ b/src/particle_refinement/refinement_criteria.jl @@ -55,7 +55,6 @@ struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE} return new{NDIMS, ELTYPE, typeof(zone_origin_function)}(zone_origin_function, spanning_set) end - end @inline Base.ndims(::RefinementCriteria{NDIMS}) where {NDIMS} = NDIMS diff --git a/src/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl index f7024bc3d..2d94c2a0d 100644 --- a/src/particle_refinement/refinement_pattern.jl +++ b/src/particle_refinement/refinement_pattern.jl @@ -1,4 +1,6 @@ function mass_distribution(system, refinement_pattern) + # TODO: + #= if refinement_pattern.center_particle # solve minimisation problem @@ -12,10 +14,11 @@ function mass_distribution(system, refinement_pattern) error("no mass conservation") else - lambda = 1 / nchilds(system, refinement_pattern) + =# + lambda = 1 / nchilds(system, refinement_pattern) - return fill(lambda, SVector{nchilds(system, refinement_pattern), eltype(system)}) - end + return fill(lambda, SVector{nchilds(system, refinement_pattern), eltype(system)}) + #end end struct CubicSplitting{ELTYPE} diff --git a/src/particle_refinement/resize.jl b/src/particle_refinement/resize.jl index d7db857c7..ccdafa32c 100644 --- a/src/particle_refinement/resize.jl +++ b/src/particle_refinement/resize.jl @@ -54,7 +54,7 @@ end capacity_parent = nparticles(system) - candidates_refine capacity_child = nparticles(system_child) + n_new_child - capacity_parent <= 0 && error("`RefinementCriteria` affects all particles") + capacity_parent < 0 && error("`RefinementCriteria` affects more than all particles") # Resize child system (extending) resize_system!(system_child, capacity_child) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 28c07fd57..d58bdd257 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -170,7 +170,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "acceleration", system.acceleration) summary_line(io, "source terms", system.source_terms |> typeof |> nameof) if !isnothing(system.particle_refinement) - summary_line(io, "refinement level", refinement_level(system.particle_refinement)) + summary_line(io, "refinement level", + refinement_level(system.particle_refinement)) end summary_footer(io) end diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index ac0f34300..a0be9f1b1 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -23,8 +23,8 @@ semi = Semidiscretization(system1, system2, neighborhood_search=nothing) # Verification - @test semi.ranges_u == (1:6, 7:18) - @test semi.ranges_v == (1:6, 7:12) + @test semi.ranges_u == ([1:6], [7:18]) + @test semi.ranges_v == ([1:6], [7:12]) nhs = ((TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(2)), TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(3))), diff --git a/test/particle_refinement/particle_refinement.jl b/test/particle_refinement/particle_refinement.jl new file mode 100644 index 000000000..899e54e30 --- /dev/null +++ b/test/particle_refinement/particle_refinement.jl @@ -0,0 +1,2 @@ +include("refinement.jl") +include("refinement_resize.jl") diff --git a/test/particle_refinement/refinement.jl b/test/particle_refinement/refinement.jl new file mode 100644 index 000000000..c42ae1478 --- /dev/null +++ b/test/particle_refinement/refinement.jl @@ -0,0 +1,223 @@ +@trixi_testset "Create System" begin + particle_spacing = 0.1 + smoothing_length = 1.2 * particle_spacing + + refinement_criterion_1 = RefinementZone(edge_length_x=1.0, edge_length_y=0.5, + zone_origin=(0.0, 0.0)) + refinement_criterion_2 = RefinementZone(edge_length_x=0.5, edge_length_y=0.75, + zone_origin=(0.5, 1.0)) + + ic = InitialCondition(; particle_spacing, coordinates=ones(2, 2), density=[1.0, 1.0]) + + @testset "Single Criteria" begin + particle_refinement = ParticleRefinement(refinement_criterion_1) + + system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + smoothing_length; + particle_refinement) + + systems = TrixiParticles.create_child_systems((system_parent,)) + + system_child = systems[2] + + particle_spacing_child = system_child.initial_condition.particle_spacing + factor = particle_refinement.refinement_pattern.separation_parameter + + @test length(systems) == 2 + @test nparticles(system_parent) == 2 + @test nparticles(system_child) == 0 + @test system_parent.particle_refinement.system_child == system_child + @test TrixiParticles.refinement_level(system_parent.particle_refinement) == 0 + @test system_child.particle_refinement isa Nothing + @test particle_spacing_child == factor * smoothing_length + # TODO: Test mass distribution + end + + @testset "Multilevel Refinement" begin + particle_refinement = ParticleRefinement(refinement_criterion_1, + criteria_next_levels=[ + refinement_criterion_1, + ]) + + system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + smoothing_length; + particle_refinement) + + systems = TrixiParticles.create_child_systems((system_parent,)) + + system_child_1 = systems[2] + system_child_2 = systems[3] + + particle_spacing_child_1 = system_child_1.initial_condition.particle_spacing + particle_spacing_child_2 = system_child_2.initial_condition.particle_spacing + factor1 = particle_refinement.refinement_pattern.separation_parameter + factor2 = factor1^2 + + @test length(systems) == 3 + @test nparticles(system_parent) == 2 + @test nparticles(system_child_1) == 0 + @test nparticles(system_child_2) == 0 + @test system_parent.particle_refinement.system_child == system_child_1 + @test system_child_1.particle_refinement.system_child == system_child_2 + @test TrixiParticles.refinement_level(system_parent.particle_refinement) == 0 + @test TrixiParticles.refinement_level(system_child_1.particle_refinement) == 1 + @test system_child_2.particle_refinement isa Nothing + @test particle_spacing_child_1 == factor1 * smoothing_length + @test particle_spacing_child_2 == factor2 * smoothing_length + # TODO: Test mass distribution + end + + @testset "Multiple Criteria" begin + particle_refinement = ParticleRefinement(refinement_criterion_1, + refinement_criterion_2) + + system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + smoothing_length; + particle_refinement) + + @test length(system_parent.particle_refinement.refinement_criteria) == 2 + + (; refinement_criteria) = system_parent.particle_refinement + @test refinement_criteria[1].zone_origin(0, 0, 0, 0, 0, 0, 0) == [0.0, 0.0] + @test refinement_criteria[2].zone_origin(0, 0, 0, 0, 0, 0, 0) == [0.5, 1.0] + end +end + +@trixi_testset "Refinement Pattern" begin + struct MockSystem <: TrixiParticles.System{2} + smoothing_length::Any + end + + Base.eltype(::MockSystem) = Float64 + + mock_system = MockSystem(1.0) + + refinement_patterns = [ + TriangularSplitting(), + CubicSplitting(), + HexagonalSplitting(), + TriangularSplitting(; center_particle=false), + CubicSplitting(; center_particle=false), + HexagonalSplitting(; center_particle=false), + ] + + expected_positions = [ + [[0.0, 0.5], [-0.4330127018922193, -0.25], [0.4330127018922193, -0.25], [0.0, 0.0]], + [ + [0.35355339059327373, 0.35355339059327373], + [0.35355339059327373, -0.35355339059327373], + [-0.35355339059327373, -0.35355339059327373], + [-0.35355339059327373, 0.35355339059327373], + [0.0, 0.0], + ], + [ + [0.5, 0.0], + [-0.5, 0.0], + [0.25, 0.4330127018922193], + [0.25, -0.4330127018922193], + [-0.25, 0.4330127018922193], + [-0.25, -0.4330127018922193], + [0.0, 0.0], + ], + [[0.0, 0.5], [-0.4330127018922193, -0.25], [0.4330127018922193, -0.25]], + [ + [0.35355339059327373, 0.35355339059327373], + [0.35355339059327373, -0.35355339059327373], + [-0.35355339059327373, -0.35355339059327373], + [-0.35355339059327373, 0.35355339059327373], + ], + [ + [0.5, 0.0], + [-0.5, 0.0], + [0.25, 0.4330127018922193], + [0.25, -0.4330127018922193], + [-0.25, 0.4330127018922193], + [-0.25, -0.4330127018922193], + ], + ] + + @testset "$(refinement_patterns[i])" for i in eachindex(refinement_patterns) + positions = TrixiParticles.relative_position_children(mock_system, + refinement_patterns[i]) + + @test expected_positions[i] ≈ positions + end +end + +@trixi_testset "Refine Particle" begin + nx = 5 + n_particles = nx^2 + particle_parent = 12 + particle_spacing = 0.1 + smoothing_length = 1.2 * particle_spacing + + nhs = TrixiParticles.TrivialNeighborhoodSearch{2}(2smoothing_length, 1:n_particles) + + ic = RectangularShape(particle_spacing, (nx, nx), (0.0, 0.0), velocity=ones(2), + mass=1.0, density=2.0) + + system_child = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + smoothing_length) + + refinement_criterion = RefinementZone(edge_length_x=Inf, edge_length_y=Inf, + zone_origin=(0.0, 0.0)) + + refinement_patterns = [ + TriangularSplitting(), + CubicSplitting(), + HexagonalSplitting(), + TriangularSplitting(; center_particle=false), + CubicSplitting(; center_particle=false), + HexagonalSplitting(; center_particle=false), + ] + + @testset "$refinement_pattern" for refinement_pattern in refinement_patterns + particle_refinement = ParticleRefinement(refinement_criterion; refinement_pattern) + + system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + smoothing_length; + particle_refinement) + + n_children = TrixiParticles.nchilds(system_parent, particle_refinement) + + resize!(system_child.mass, n_children) + + relative_positions = TrixiParticles.relative_position_children(system_parent, + refinement_pattern) + mass_ratios = TrixiParticles.mass_distribution(system_parent, + refinement_pattern) + + particle_refinement.rel_position_children = relative_positions + particle_refinement.available_children = n_children + particle_refinement.mass_ratio = mass_ratios + + v_parent = vcat(ic.velocity, ic.density') + u_parent = copy(ic.coordinates) + + v_child = Array{Float64, 2}(undef, 3, n_children) + u_child = Array{Float64, 2}(undef, 2, n_children) + + mass_parent = system_parent.mass[particle_parent] + + TrixiParticles.bear_children!(system_child, system_parent, particle_parent, + mass_parent, nhs, particle_refinement, + v_parent, u_parent, v_child, u_child) + + parent_pos = u_parent[:, particle_parent] + for child in 1:n_children + @test u_child[:, child] ≈ parent_pos .+ relative_positions[child] + @test v_child[:, child] == [1.0, 1.0, 2.0] + @test system_child.mass[child] == mass_parent * mass_ratios[child] + end + end +end diff --git a/test/particle_refinement/refinement_resize.jl b/test/particle_refinement/refinement_resize.jl new file mode 100644 index 000000000..d0f190908 --- /dev/null +++ b/test/particle_refinement/refinement_resize.jl @@ -0,0 +1,155 @@ +@trixi_testset "Resize and Copy" begin + n_particles = 10 + dict_ = Dict(12 => "Filled", 0 => "Empty") + + @testset "$(dict_[n_particles_child]) Child System" for n_particles_child in [0, 12] + refinement_patterns = [ + TriangularSplitting(), + CubicSplitting(), + HexagonalSplitting(), + ] + + candidates = [[1, 10], [1, 3, 5, 7, 9], [2, 4, 6, 8], collect(1:n_particles)] + + old_ranges_v = ([1:(3n_particles)], + [(3n_particles + 1):(3n_particles + 3n_particles_child)]) + old_ranges_u = ([1:(2n_particles)], + [(2n_particles + 1):(2n_particles + 2n_particles_child)]) + + old_nparticles = (n_particles, n_particles_child) + + ic = InitialCondition(; particle_spacing=0.1, coordinates=ones(2, n_particles), + density=ones(n_particles)) + + refinement_criterion = RefinementZone(edge_length_x=Inf, edge_length_y=Inf, + zone_origin=(0.0, 0.0)) + + refinement_callback = ParticleRefinementCallback(interval=1) + callback = refinement_callback.affect! + + @testset "Refinement Pattern $refinement_pattern" for refinement_pattern in refinement_patterns + @testset "Refinement Candidates $j" for j in eachindex(candidates) + particle_refinement = ParticleRefinement(refinement_criterion; + refinement_pattern=refinement_pattern) + + system_parent = WeaklyCompressibleSPHSystem(ic, ContinuityDensity(), + nothing, + SchoenbergCubicSplineKernel{2}(), + 1.0; particle_refinement) + + semi = Semidiscretization(system_parent) + + # Resize emtpy child system + resize!(semi.systems[2].mass, n_particles_child) + + # Add candidates + system_parent.particle_refinement.candidates = candidates[j] + + # Create vectors filled with the corresponding particle index + v_ode_parent = reshape(stack([i * ones(Int, 3) for i in 1:n_particles]), + (3 * n_particles,)) + u_ode_parent = reshape(stack([i * ones(Int, 2) for i in 1:n_particles]), + (2 * n_particles,)) + + if n_particles_child > 0 + v_ode_child = reshape(stack([i * ones(Int, 3) + for i in 1:n_particles_child]), + (3 * n_particles_child,)) + u_ode_child = reshape(stack([i * ones(Int, 2) + for i in 1:n_particles_child]), + (2 * n_particles_child,)) + + v_ode = vcat(v_ode_parent, v_ode_child) + u_ode = vcat(u_ode_parent, u_ode_child) + else + v_ode = v_ode_parent + u_ode = u_ode_parent + end + + _v_cache = copy(v_ode) + _u_cache = copy(u_ode) + + TrixiParticles.resize_and_copy!(callback, semi, v_ode, u_ode, + _v_cache, _u_cache) + + # Iterator without candidates + eachparticle_excluding_candidates = (setdiff(1:n_particles, candidates[j]), + Base.OneTo(n_particles_child)) + + n_candidates = length(candidates[j]) + n_children = TrixiParticles.nchilds(system_parent, particle_refinement) + n_total_particles = n_particles + n_particles_child - n_candidates + + n_candidates * n_children + + # Ranges after resizing + new_ranges_v = ([1:(3nparticles(system_parent))], + [(3nparticles(system_parent) + 1):(3n_total_particles)]) + new_ranges_u = ([1:(2nparticles(system_parent))], + [(2nparticles(system_parent) + 1):(2n_total_particles)]) + + @testset "Iterators And Ranges" begin + @test callback.nparticles_cache == old_nparticles + @test callback.ranges_v_cache == old_ranges_v + @test callback.ranges_u_cache == old_ranges_u + @test callback.eachparticle_cache == eachparticle_excluding_candidates + + @test nparticles(system_parent) == n_particles - n_candidates == + length(eachparticle_excluding_candidates[1]) + @test nparticles(semi.systems[2]) == + n_candidates * n_children + n_particles_child + + @test semi.ranges_v == new_ranges_v + @test semi.ranges_u == new_ranges_u + end + + @testset "Resized Integrator-Array" begin + # Parent system + v_parent = TrixiParticles.wrap_v(v_ode, system_parent, semi) + + @test size(v_parent, 2) == length(eachparticle_excluding_candidates[1]) + + # Test `copy_values_v!` + for particle in 1:nparticles(system_parent) + for dim in 1:3 + @test v_parent[dim, particle] == + eachparticle_excluding_candidates[1][particle] + end + end + + u_parent = TrixiParticles.wrap_u(u_ode, system_parent, semi) + + @test size(u_parent, 2) == length(eachparticle_excluding_candidates[1]) + + # Test `copy_values_u!` + for particle in 1:nparticles(system_parent) + for dim in 1:ndims(system_parent) + @test u_parent[dim, particle] == + eachparticle_excluding_candidates[1][particle] + end + end + + # Child system + v_child = TrixiParticles.wrap_v(v_ode, semi.systems[2], semi) + @test size(v_child, 2) == n_candidates * n_children + n_particles_child + + # Test `copy_values_v!` + for particle in 1:n_particles_child + for dim in 1:3 + @test v_child[dim, particle] == particle + end + end + + u_child = TrixiParticles.wrap_u(u_ode, semi.systems[2], semi) + @test size(u_child, 2) == n_candidates * n_children + n_particles_child + + # Test `copy_values_u!` + for particle in 1:n_particles_child + for dim in 1:2 + @test u_child[dim, particle] == particle + end + end + end + end + end + end +end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 0749261bd..83bffd45b 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -193,7 +193,7 @@ smoothing_length, density_diffusion=density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), [0.0, 0.0], nothing, nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ diff --git a/test/unittest.jl b/test/unittest.jl index 01308bd44..06878d714 100644 --- a/test/unittest.jl +++ b/test/unittest.jl @@ -7,4 +7,5 @@ include("setups/setups.jl") include("systems/systems.jl") include("schemes/schemes.jl") + include("particle_refinement/particle_refinement.jl") end;