diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index d6e0c2996..9f7e685d5 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -20,7 +20,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 @@ -50,6 +51,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("particle_refinement/particle_refinement.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition @@ -57,7 +59,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, OutFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, - PostprocessCallback, StepsizeCallback, UpdateCallback + PostprocessCallback, StepsizeCallback, UpdateCallback, ParticleRefinementCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, @@ -84,5 +86,7 @@ export kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, export interpolate_line, interpolate_point, interpolate_plane_3d, interpolate_plane_2d, interpolate_plane_2d_vtk export SurfaceTensionAkinci, CohesionForceAkinci +export ParticleRefinement, RefinementZone +export CubicSplitting, TriangularSplitting, HexagonalSplitting end # module diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 1aefdff72..ed573bbf7 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") include("update.jl") +include("refinement_callback.jl") diff --git a/src/callbacks/refinement_callback.jl b/src/callbacks/refinement_callback.jl new file mode 100644 index 000000000..d2b829d86 --- /dev/null +++ b/src/callbacks/refinement_callback.jl @@ -0,0 +1,145 @@ +mutable struct ParticleRefinementCallback{I} + interval :: I + ranges_u_cache :: Tuple + ranges_v_cache :: Tuple + nparticles_cache :: Tuple + eachparticle_cache :: Tuple +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 + + refinement_callback = ParticleRefinementCallback(interval, (), (), (), ()) + + if dt > 0 + # Add a `tstop` every `dt`, and save the final solution. + 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(refinement_callback, refinement_callback, + initialize=initial_refinement!, + save_positions=(false, false)) + end +end + +# initialize +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_refinement!(cb.affect!, u, t, integrator) +end + +function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator) + cb(integrator) +end + +# condition +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 integrator.stats.naccept % interval == 0 +end + +# affect +function (refinement_callback::ParticleRefinementCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # 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 + + refinement!(v_ode, u_ode, v_tmp, u_tmp, semi, refinement_callback, t) + + # 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 + + @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) + + 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), ")") +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 + refinement_cb = cb.affect! + setup = [ + "interval" => refinement_cb.interval, + ] + 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 + refinement_cb = cb.affect!.affect! + setup = [ + "dt" => refinement_cb.interval, + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..cf6c8c26c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -73,14 +73,9 @@ function Semidiscretization(systems...; # Other checks might be added here later. check_configuration(systems) - sizes_u = [u_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) - for i in eachindex(sizes_u)) - sizes_v = [v_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) - for i in eachindex(sizes_v)) + systems = create_child_systems(systems) + + ranges_v, ranges_u = ranges_vu(systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. @@ -92,6 +87,21 @@ function Semidiscretization(systems...; return Semidiscretization(systems, ranges_u, ranges_v, searches) end +function ranges_vu(systems) + 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)) + + # Each range is stored in a `Vector`, so we can mutate ranges in an immutable struct + return ranges_v, ranges_u +end + # Inline show function e.g. Semidiscretization(neighborhood_search=...) function Base.show(io::IO, semi::Semidiscretization) @nospecialize semi # reduce precompilation time @@ -370,7 +380,9 @@ end @inline function wrap_v(v_ode, system, semi) (; ranges_v) = semi - range = ranges_v[system_indices(system, semi)] + # 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) @@ -381,7 +393,9 @@ end @inline function wrap_u(u_ode, system, semi) (; ranges_u) = semi - range = ranges_u[system_indices(system, semi)] + # 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) 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/particle_refinement/refinement.jl b/src/particle_refinement/refinement.jl new file mode 100644 index 000000000..5ac228b27 --- /dev/null +++ b/src/particle_refinement/refinement.jl @@ -0,0 +1,317 @@ +include("refinement_pattern.jl") +include("refinement_criteria.jl") + +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 :: CNL + available_children :: Int + + # Depends on refinement pattern, particle spacing and parameters ϵ and α. + # 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_child_system()` at the beginning of the simulation + system_child::System + + # API --> parent system with `RL=0` + function ParticleRefinement(refinement_criteria...; + refinement_pattern=CubicSplitting(), + criteria_next_levels=[]) + ELTYPE = eltype(first(refinement_criteria)) + NDIMS = ndims(first(refinement_criteria)) + + return new{0, NDIMS, ELTYPE, typeof(refinement_pattern), + 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, 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), + typeof(criteria_next_levels)}([], [], refinement_pattern, + refinement_criteria, + criteria_next_levels, 0) + end +end + +@inline refinement_level(::ParticleRefinement{RL}) where {RL} = RL + +@inline child_set(system, particle_refinement) = Base.OneTo(nchilds(system, + particle_refinement)) + +@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_ = () + foreach_system(systems) do system + systems_ = (systems_..., create_child_system(system)...) + end + + return (systems..., systems_...) +end + +create_child_system(system) = () +function create_child_system(system::FluidSystem) + create_child_system(system, system.particle_refinement) +end + +create_child_system(system::FluidSystem, ::Nothing) = () + +function create_child_system(system::FluidSystem, + particle_refinement::ParticleRefinement{RL}) where {RL} + (; smoothing_length) = system + (; criteria_next_levels, refinement_pattern, refinement_criteria) = particle_refinement + + NDIMS = ndims(system) + + # Distribute values according to 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_ = smoothing_length * refinement_pattern.separation_parameter + 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 isempty(criteria_next_levels) + nothing + else + refinement_criteria = first(criteria_next_levels) + ParticleRefinement{level}(refinement_criteria, refinement_pattern, + criteria_next_levels[(level + 1):end]) + 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) + + particle_refinement.system_child = system_child + + return (system_child, + create_child_system(system_child, system_child.particle_refinement)...) +end + +# ==== Refinement +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, t) + end + + resize_and_copy!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) + + refine_particles!(callback, semi, v_ode, u_ode, _v_cache, _u_cache) +end + +check_refinement_criteria!(system, v_ode, u_ode, semi, t) = system + +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, t) = system + +function check_refinement_criteria!(system, particle_refinement::ParticleRefinement, + v_ode, u_ode, semi, t) + (; candidates, candidates_mass, refinement_criteria) = particle_refinement + + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + 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, t) + 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 + + return system +end + +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, _v_cache, _u_cache, callback, semi) + end +end + +refine_particles!(system, v_ode, u_ode, _v_cache, _u_cache, callback, semi) = system + +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 + +function refine_particles!(system::FluidSystem, ::Nothing, v_ode, u_ode, _v_cache, _u_cache, + callback, semi) + return system +end + +function refine_particles!(system_parent::FluidSystem, + particle_refinement::ParticleRefinement, + v_ode, u_ode, _v_cache, _u_cache, callback, semi) + (; 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) + + # Resized storage + v_child = wrap_v(v_ode, system_child, semi) + u_child = wrap_u(u_ode, system_child, semi) + + particle_refinement.available_children = length(candidates) * + nchilds(system_parent, particle_refinement) + + # Loop over all refinement candidates + mass_index = 1 + for particle_parent in candidates + mass_parent = candidates_mass[mass_index] + 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) + mass_index += 1 + end + 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_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 + + 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_children + + 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_children[particle_child] + for dim in 1:ndims(system_child) + u_child[dim, absolute_index] = child_coords[dim] + end + + volume = zero(eltype(system_child)) + 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 + + 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) + volume += kernel_weight + + v_b = current_velocity(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, absolute_index] += kernel_weight * v_b[dim] + end + + rho_a += kernel_weight * rho_b + p_a += kernel_weight * p_b + end + end + + if volume > eps() + for dim in 1:ndims(system_child) + v_child[dim, absolute_index] /= volume + end + + rho_a /= volume + p_a /= volume + + 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 +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/particle_refinement/refinement_criteria.jl b/src/particle_refinement/refinement_criteria.jl new file mode 100644 index 000000000..b777543c6 --- /dev/null +++ b/src/particle_refinement/refinement_criteria.jl @@ -0,0 +1,82 @@ +# Criteria of refinement: +# +# - fixed (/moving?) refinement zone +# - number of neighbors +# - problem specific criteria (e.g. high velocity gradient) + +abstract type RefinementCriteria{NDIMS, ELTYPE} end + +struct RefinementZone{NDIMS, ELTYPE, ZO} <: RefinementCriteria{NDIMS, ELTYPE} + zone_origin :: ZO + spanning_set :: Vector{SVector} + + 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 + else + 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_ = 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 + +@inline function (refinement_criterion::RefinementZone)(system, particle, + v, u, v_ode, u_ode, semi, t) + (; zone_origin, spanning_set) = refinement_criterion + 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] + # 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/particle_refinement/refinement_pattern.jl b/src/particle_refinement/refinement_pattern.jl new file mode 100644 index 000000000..2d94c2a0d --- /dev/null +++ b/src/particle_refinement/refinement_pattern.jl @@ -0,0 +1,136 @@ +function mass_distribution(system, refinement_pattern) + # TODO: + #= + 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} + separation_parameter :: ELTYPE + smoothing_ratio :: ELTYPE + center_particle :: Bool + + 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) + (; smoothing_length) = system + (; separation_parameter) = refinement_pattern + + direction_1 = [1 / sqrt(2), 1 / sqrt(2)] + direction_2 = [1 / sqrt(2), -1 / sqrt(2)] + direction_3 = -direction_1 + direction_4 = -direction_2 + + # Is it `particle_spacing * separation_parameter` or `smoothing_length * separation_parameter`? + 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]) + end + + return reinterpret(reshape, SVector{2, eltype(system)}, relative_position) +end + +struct TriangularSplitting{ELTYPE} + separation_parameter :: ELTYPE + smoothing_ratio :: ELTYPE + center_particle :: Bool + + function TriangularSplitting(; 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::TriangularSplitting) = 3 + + refinement_pattern.center_particle + +function relative_position_children(system::System{2}, + refinement_pattern::TriangularSplitting) + (; 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(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]) + end + + 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) + (; smoothing_length) = system + (; 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(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]) + end + + return reinterpret(reshape, SVector{2, eltype(system)}, relative_position) +end diff --git a/src/particle_refinement/resize.jl b/src/particle_refinement/resize.jl new file mode 100644 index 000000000..ccdafa32c --- /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 more than 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 diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index a3f47f611..2df65f593 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -376,6 +376,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 diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index cf73f5876..86c7c3849 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -42,7 +42,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC @@ -56,6 +56,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, pressure_acceleration_formulation :: PF source_terms :: ST buffer :: B + particle_refinement :: PR cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, @@ -65,6 +66,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), + particle_refinement=nothing, source_terms=nothing, buffer_size=nothing) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -97,11 +99,11 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(pressure_acceleration), typeof(source_terms), - typeof(buffer), + typeof(buffer), 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, - buffer, cache) + buffer, particle_refinement, cache) end end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 9e2289a8c..fec09f041 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -80,3 +80,66 @@ include("viscosity.jl") include("surface_tension.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) = particle + +@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) = particle + +@inline function set_particle_pressure(particle, v, system::EntropicallyDampedSPHSystem, + 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 diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 482261fc1..3534fb5a5 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,8 +44,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} +struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, + COR, PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -61,6 +61,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, source_terms :: ST surface_tension :: SRFT buffer :: B + particle_refinement :: PR cache :: C end @@ -74,6 +75,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), + particle_refinement=nothing, correction=nothing, source_terms=nothing, surface_tension=nothing) buffer = isnothing(buffer_size) ? nothing : @@ -121,8 +123,8 @@ function WeaklyCompressibleSPHSystem(initial_condition, smoothing_kernel, smoothing_length, acceleration_, viscosity, density_diffusion, correction, - pressure_acceleration, - source_terms, surface_tension, buffer, cache) + pressure_acceleration, source_terms, surface_tension, + buffer, particle_refinement, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -167,6 +169,7 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.surface_tension) print(io, ", ", system.acceleration) print(io, ", ", system.source_terms) + print(io, ", ", system.particle_refinement) print(io, ") with ", nparticles(system), " particles") end @@ -194,6 +197,10 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "surface tension", system.surface_tension) 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 diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 63d7628f4..034fe8bf0 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -107,36 +107,39 @@ function trixi2vtk(v_, u_, t, system_, periodic_box; output_directory="out", pre 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 @trixi_timeit timer() "write to vtk" vtk_grid(file, points, cells) do vtk - # dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem - write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) + if !isempty(points) + # dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem + write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) - # Store particle index - vtk["index"] = active_particles(system) - vtk["time"] = t + # Store particle index + vtk["index"] = active_particles(system) + vtk["time"] = t - if write_meta_data - vtk["solver_version"] = git_hash - vtk["julia_version"] = string(VERSION) - end + if write_meta_data + vtk["solver_version"] = 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 diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 11a2d7f56..e203b6caf 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -29,8 +29,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}(search_radius=0.2, eachpoint=1:2), 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 68f69aea8..0539d79fa 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}(), nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, [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 c344173b4..0ce7a3fa6 100644 --- a/test/unittest.jl +++ b/test/unittest.jl @@ -6,5 +6,6 @@ include("setups/setups.jl") include("systems/systems.jl") include("schemes/schemes.jl") + include("particle_refinement/particle_refinement.jl") include("preprocessing/preprocessing.jl") end;