diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 66b3d20f9e..1611719a03 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -77,7 +77,7 @@ export PenaltyForceGanzenmueller, TransportVelocityAdami, ParticleShiftingTechni ContinuityEquationTermSun2019, MomentumEquationTermSun2019 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, SchoenbergQuinticSplineKernel, GaussianKernel, WendlandC2Kernel, WendlandC4Kernel, - WendlandC6Kernel, SpikyKernel, Poly6Kernel + WendlandC6Kernel, SpikyKernel, Poly6Kernel, MarroneMLSKernel export StateEquationCole, StateEquationIdealGas export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris, ViscosityAdamiSGS, ViscosityMorrisSGS @@ -86,7 +86,7 @@ export tensile_instability_control export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BoundaryModelCharacteristicsLastiwka, BoundaryModelMirroringTafuni, BoundaryModelDynamicalPressureZhang, - BernoulliPressureExtrapolation, PressureBoundaries + BernoulliPressureExtrapolation, PressureBoundaries, MarronePressureExtrapolation export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring export HertzContactModel, LinearContactModel export PrescribedMotion, OscillatingMotion2D diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 4eb1af8f00..549e90b5f0 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -96,27 +96,37 @@ struct InitialCondition{ELTYPE, MATRIX, VECTOR} mass :: VECTOR # Array{ELTYPE, 1} density :: VECTOR # Array{ELTYPE, 1} pressure :: VECTOR # Array{ELTYPE, 1} + normals :: Union{Nothing, MATRIX} end # The default constructor needs to be accessible for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. function InitialCondition(; coordinates, density, velocity=zeros(size(coordinates, 1)), - mass=nothing, pressure=0.0, particle_spacing=-1.0) + mass=nothing, pressure=0.0, particle_spacing=-1.0, + normals=nothing) NDIMS = size(coordinates, 1) return InitialCondition{NDIMS}(coordinates, velocity, mass, density, - pressure, particle_spacing) + pressure, particle_spacing, normals) end -# Function barrier to make `NDIMS` static and therefore SVectors type-stable +# Wrapper function for outer constructor without `normals`-keyword function InitialCondition{NDIMS}(coordinates, velocity, mass, density, pressure, particle_spacing) where {NDIMS} + return InitialCondition{NDIMS}(coordinates, velocity, mass, density, + pressure, particle_spacing, nothing) +end + +# Function barrier to make `NDIMS` static and therefore SVectors type-stable +function InitialCondition{NDIMS}(coordinates, velocity, mass, density, + pressure, particle_spacing, normals) where {NDIMS} ELTYPE = eltype(coordinates) n_particles = size(coordinates, 2) if n_particles == 0 return InitialCondition(particle_spacing, coordinates, zeros(ELTYPE, NDIMS, 0), - zeros(ELTYPE, 0), zeros(ELTYPE, 0), zeros(ELTYPE, 0)) + zeros(ELTYPE, 0), zeros(ELTYPE, 0), zeros(ELTYPE, 0), + normals) end # SVector of coordinates to pass to functions. @@ -187,8 +197,15 @@ function InitialCondition{NDIMS}(coordinates, velocity, mass, density, masses = mass_fun.(coordinates_svector) end + if normals isa AbstractMatrix + if size(coordinates) != size(normals) + throw(ArgumentError("`coordinates` and `normals` must be of the same size")) + end + end + return InitialCondition(particle_spacing, coordinates, ELTYPE.(velocities), - ELTYPE.(masses), ELTYPE.(densities), ELTYPE.(pressures)) + ELTYPE.(masses), ELTYPE.(densities), ELTYPE.(pressures), + normals) end function Base.show(io::IO, ic::InitialCondition) diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index 192eb5fee1..270a966a3a 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -64,9 +64,12 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, ELTYPE = eltype(smoothing_length) n_particles = length(initial_density) - cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., - create_cache_model(initial_density, density_calculator, NDIMS)..., - create_cache_model(correction, initial_density, NDIMS, n_particles)...) + cache = (; create_cache_boundary_velocity(viscosity, n_particles, NDIMS)..., + create_cache_boundary_density(initial_density, density_calculator, NDIMS)..., + create_cache_boundary_correction(correction, initial_density, NDIMS, + n_particles)..., + create_cache_boundary_interpolation(density_calculator, n_particles, NDIMS, + eltype(initial_density))...) # If the `reference_density_spacing` is set calculate the `ideal_neighbor_count` if reference_particle_spacing > 0 @@ -200,39 +203,63 @@ struct PressureBoundaries{ELTYPE} return new{eltype(time_step)}(time_step, omega) end end -@inline create_cache_model(correction, density, NDIMS, nparticles) = (;) -function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particles) +@doc raw""" + MarronePressureExtrapolation() + +`density_calculator` for `BoundaryModelDummyParticles`. + +!!! note + This boundary model was orignally proposed in [𝛿-SPH Model for Simulating Violent Impact Flows](https://www.researchgate.net/publication/241077909_-SPH_model_for_simulating_violent_impact_flows). +""" +struct MarronePressureExtrapolation + allow_loop_flipping::Bool + + function MarronePressureExtrapolation(; allow_loop_flipping=true) + return new(allow_loop_flipping) + end +end + +@inline create_cache_boundary_correction(correction, density, NDIMS, nparticles) = (;) + +function create_cache_boundary_correction(::ShepardKernelCorrection, density, NDIMS, + n_particles) return (; kernel_correction_coefficient=similar(density)) end -function create_cache_model(::KernelCorrection, density, NDIMS, n_particles) +function create_cache_boundary_correction(::KernelCorrection, density, NDIMS, n_particles) dw_gamma = Array{Float64}(undef, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density), dw_gamma) end -function create_cache_model(::Union{GradientCorrection, BlendedGradientCorrection}, density, - NDIMS, n_particles) +function create_cache_boundary_correction(::Union{GradientCorrection, + BlendedGradientCorrection}, + density, + NDIMS, n_particles) correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) return (; correction_matrix) end -function create_cache_model(::MixedKernelGradientCorrection, density, NDIMS, n_particles) +function create_cache_boundary_correction(::MixedKernelGradientCorrection, density, NDIMS, + n_particles) dw_gamma = Array{Float64}(undef, NDIMS, n_particles) correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end -function create_cache_model(initial_density, - ::Union{SummationDensity, PressureMirroring, PressureZeroing}, - NDIMS) +function create_cache_boundary_density(initial_density, + ::Union{SummationDensity, PressureMirroring, + PressureZeroing}, NDIMS) density = copy(initial_density) return (; density) end -function create_cache_model(initial_density, - density_calculator::PressureBoundaries, NDIMS) +@inline create_cache_boundary_density(initial_density, + ::ContinuityDensity, NDIMS) = (; initial_density) + +function create_cache_boundary_density(initial_density, + density_calculator::PressureBoundaries, NDIMS) ELTYPE = eltype(initial_density) n_particles = size(initial_density, 1) @@ -253,22 +280,19 @@ function create_cache_model(initial_density, omega, time_step, density_error) end -function create_cache_model(initial_density, ::ContinuityDensity, NDIMS) - return (; initial_density) -end - -function create_cache_model(initial_density, - ::Union{AdamiPressureExtrapolation, - BernoulliPressureExtrapolation}, NDIMS) +function create_cache_boundary_density(initial_density, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation, + MarronePressureExtrapolation}, NDIMS) density = copy(initial_density) volume = similar(initial_density) return (; density, volume) end -@inline create_cache_model(viscosity::Nothing, n_particles, n_dims) = (;) +@inline create_cache_boundary_velocity(viscosity::Nothing, n_particles, n_dims) = (;) -function create_cache_model(viscosity, n_particles, n_dims) +function create_cache_boundary_velocity(viscosity, n_particles, n_dims) ELTYPE = eltype(viscosity.epsilon) wall_velocity = zeros(ELTYPE, n_dims, n_particles) @@ -276,17 +300,39 @@ function create_cache_model(viscosity, n_particles, n_dims) return (; wall_velocity) end -@inline reset_cache!(cache, viscosity) = set_zero!(cache.volume) +@inline create_cache_boundary_interpolation(density_calculator, n_particles, n_dims, + ELTYPE) = (;) -function reset_cache!(cache, viscosity::ViscosityAdami) - (; volume, wall_velocity) = cache +function create_cache_boundary_interpolation(density_calculator::MarronePressureExtrapolation, + n_particles, n_dims, ELTYPE) + # TODO: `wall_velocity` should be created in `create_cache_velocity` + wall_velocity = zeros(ELTYPE, n_dims, n_particles) + interpolation_coords = fill(zero(SVector{n_dims, ELTYPE}), n_particles) + _pressure = zeros(ELTYPE, n_particles) + + return (; interpolation_coords, wall_velocity, _pressure) +end + +function reset_cache!(model::BoundaryModelDummyParticles) + (; cache, viscosity, density_calculator) = model + (; volume) = cache set_zero!(volume) - set_zero!(wall_velocity) + reset_cache!(cache, viscosity) + reset_cache!(cache, density_calculator) return cache end +@inline reset_cache!(_, _) = nothing + +@inline reset_cache!(cache, viscosity::ViscosityAdami) = set_zero!(cache.wall_velocity) + +function reset_cache!(cache, density_calculator::MarronePressureExtrapolation) + set_zero!(cache.wall_velocity) + set_zero!(cache._pressure) +end + function Base.show(io::IO, model::BoundaryModelDummyParticles) @nospecialize model # reduce precompilation time @@ -320,7 +366,7 @@ end ::Union{SummationDensity, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BernoulliPressureExtrapolation, - PressureBoundaries}, + PressureBoundaries, MarronePressureExtrapolation}, model::BoundaryModelDummyParticles) # When using `SummationDensity`, the density is stored in the cache return model.cache.density @@ -348,10 +394,11 @@ end function compute_density!(boundary_model, ::Union{ContinuityDensity, AdamiPressureExtrapolation, BernoulliPressureExtrapolation, + MarronePressureExtrapolation, PressureMirroring, PressureZeroing}, system, v, u, v_ode, u_ode, semi) # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`. - # For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. + # For `AdamiPressureExtrapolation`, `MarronePressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. return boundary_model end @@ -441,7 +488,8 @@ end function compute_pressure!(boundary_model, ::Union{AdamiPressureExtrapolation, - BernoulliPressureExtrapolation}, + BernoulliPressureExtrapolation, + MarronePressureExtrapolation}, system, v, u, v_ode, u_ode, semi) (; pressure, cache, viscosity) = boundary_model (; allow_loop_flipping) = boundary_model.density_calculator @@ -449,7 +497,7 @@ function compute_pressure!(boundary_model, set_zero!(pressure) # Set `volume` to zero. For `ViscosityAdami` the `wall_velocity` is also set to zero. - reset_cache!(cache, viscosity) + reset_cache!(boundary_model) system_coords = current_coordinates(u, system) diff --git a/src/schemes/boundary/wall_boundary/marrone.jl b/src/schemes/boundary/wall_boundary/marrone.jl new file mode 100644 index 0000000000..f0c56febf0 --- /dev/null +++ b/src/schemes/boundary/wall_boundary/marrone.jl @@ -0,0 +1,184 @@ +@doc raw""" + MarroneMLSKernel{NDIMS}() + +The Moving Least-Squares Kernel by Marrone et al. is used to compute the pressure of dummy particles for `MarronePressureExtrapolation`. +""" + +struct MarroneMLSKernel{NDIMS} <: AbstractSmoothingKernel{NDIMS} + inner_kernel :: AbstractSmoothingKernel + basis :: Array{Float64, 3} + momentum :: Array{Float64, 3} +end + +function MarroneMLSKernel(inner_kernel::AbstractSmoothingKernel{NDIMS}, + n_boundary_particles, n_fluid_particles) where {NDIMS} + basis = zeros(n_boundary_particles, n_fluid_particles, NDIMS+1) # Big sparse tensor + momentum = zeros(n_boundary_particles, NDIMS+1, NDIMS+1) + + return MarroneMLSKernel{NDIMS}(inner_kernel, basis, momentum) +end + +# Compute the Marrone MLS-Kernel for the points with indices i and j +@inline function boundary_kernel_marrone(smoothing_kernel::MarroneMLSKernel, i, j, distance, + smoothing_length) + (; inner_kernel, basis, momentum) = smoothing_kernel + ndims = size(momentum, 2) + canonical_vector = ndims == 3 ? [1.0, 0.0, 0.0] : [1.0, 0.0, 0.0, 0.0] + + kernel_weight = kernel(inner_kernel, distance, smoothing_length) + M_inv = momentum[i, :, :] + + return dot((M_inv * canonical_vector), (basis[i, j, :] * kernel_weight)) +end + +function compute_basis_marrone(kernel::MarroneMLSKernel, system, neighbor_system, + system_coords, neighbor_coords, semi) + (; basis) = kernel + basis .= 0.0 + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; + points=eachparticle(system)) do particle, neighbor, + pos_diff, distance + basis[particle, neighbor, :] = [1; -pos_diff] # We use -pos_diff since we need to calculate coords_neighbor - coords_particle + end +end + +function compute_momentum_marrone(kernel::MarroneMLSKernel, system, neighbor_system, + system_coords, neighbor_coords, v_neighbor_system, semi, + smoothing_length) + (; inner_kernel, basis, momentum) = kernel + momentum .= 0.0 + + foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; + points=eachparticle(system)) do particle, neighbor, + pos_diff, distance + kernel_weight = TrixiParticles.kernel(inner_kernel, distance, smoothing_length) + density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) + volume_neighbor = density_neighbor != 0 ? + hydrodynamic_mass(neighbor_system, neighbor) / density_neighbor : + 0 + momentum[particle, :, + :] += basis[particle, neighbor, :] .* + basis[particle, neighbor, :]' .* + kernel_weight .* volume_neighbor + end + + n = size(momentum, 2) + for particle in eachparticle(system) + if abs(det(momentum[particle, :, :])) < 1.0f-9 + momentum[particle, :, + :] = Matrix{eltype(momentum[particle, :, :])}(I, n, n) # Create identity matrix of fitting size and type + else + momentum[particle, :, :] = inv(momentum[particle, :, :]) + end + end +end + +# I dont know if this is correct... +@inline compact_support(kernel::MarroneMLSKernel, + h) = compact_support(kernel.inner_kernel, h) + +struct DummyKernel{NDIMS} <: AbstractSmoothingKernel{NDIMS} end + +function kernel(kernel::DummyKernel, r::Real, h) + return 1 +end + +@inline function boundary_pressure_extrapolation!(parallel::Val{true}, + boundary_model::BoundaryModelDummyParticles{MarronePressureExtrapolation}, + system, + neighbor_system::AbstractFluidSystem, + system_coords, neighbor_coords, v, + v_neighbor_system, + semi) + (; pressure, cache, viscosity, density_calculator, smoothing_kernel, + smoothing_length) = boundary_model + (; interpolation_coords, _pressure) = cache + + # Interpret Vector of SVectors as Matrix instead to work with the neighborhood search + interpolation_coords_view = reinterpret(reshape, eltype(eltype(interpolation_coords)), + interpolation_coords) + + compute_basis_marrone(smoothing_kernel, system, neighbor_system, + interpolation_coords_view, + neighbor_coords, semi) + compute_momentum_marrone(smoothing_kernel, system, neighbor_system, + interpolation_coords_view, + neighbor_coords, v_neighbor_system, semi, + smoothing_length) + + # Loop over all pairs of interpolation points and fluid particles within the kernel cutoff + foreach_point_neighbor(system, neighbor_system, interpolation_coords_view, + neighbor_coords, + semi) do particle, neighbor, + pos_diff, distance + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, _pressure) + end + + # Copy the updated pressure values from the buffer + pressure .= _pressure +end + +@inline function boundary_pressure_inner!(boundary_model, + boundary_density_calculator::MarronePressureExtrapolation, + system, neighbor_system::AbstractFluidSystem, v, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure) + (; smoothing_kernel, smoothing_length) = boundary_model + kernel_weight = boundary_kernel_marrone(smoothing_kernel, particle, neighbor, distance, + smoothing_length) + # Update the pressure + neighbor_pressure = current_pressure(v_neighbor_system, neighbor_system, neighbor) + neighbor_density = current_density(v_neighbor_system, neighbor_system, neighbor) + neighbor_volume = neighbor_density != 0 ? + hydrodynamic_mass(neighbor_system, neighbor) / neighbor_density : 0 + + pressure[particle] += neighbor_pressure * kernel_weight * neighbor_volume + + # Update the boundary particle velocity + # TODO: This method takes an additional parameter `neighbor_volume`, maybe unify all methods of this function + # to accept the same arguments + compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, + kernel_weight, particle, neighbor, neighbor_volume) +end + +# Change the dispatched type of `viscosity` so that it also accepts Marrone +function compute_smoothed_velocity!(cache, viscosity, neighbor_system, + v_neighbor_system, kernel_weight, particle, neighbor, + neighbor_volume) + neighbor_velocity = current_velocity(v_neighbor_system, neighbor_system, neighbor) + + # CHECK: is the neighbor_volume term necessary? + for dim in eachindex(neighbor_velocity) + @inbounds cache.wall_velocity[dim, + particle] += neighbor_velocity[dim] * kernel_weight * + neighbor_volume + end + + return cache +end + +# See `dummy_particles.jl` +function compute_boundary_density!(boundary_model::BoundaryModelDummyParticles{MarronePressureExtrapolation}, + system, system_coords, particle) + (; pressure, state_equation, cache, viscosity) = boundary_model + (; volume, density) = cache + + # The summation is only over fluid particles, thus the volume stays zero when a boundary + # particle isn't surrounded by fluid particles. + # Check the volume to avoid NaNs in pressure and velocity. + particle_volume = volume[particle] + if @inbounds particle_volume > eps() + # To impose no-slip condition + compute_wall_velocity!(viscosity, system, system_coords, particle) + end + + # Limit pressure to be non-negative to avoid attractive forces between fluid and + # boundary particles at free surfaces (sticking artifacts). + @inbounds pressure[particle] = max(pressure[particle], 0) + + # Apply inverse state equation to compute density (not used with EDAC) + # CHECK: this computes the density based on the updated pressure, is this correct? + inverse_state_equation!(density, state_equation, pressure, particle) +end diff --git a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl index 1df9a38cf9..80d4d8035d 100644 --- a/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl +++ b/src/schemes/boundary/wall_boundary/monaghan_kajtar.jl @@ -64,10 +64,10 @@ end return K / beta^(ndims(particle_system) - 1) * pos_diff / (distance * distance_from_singularity) * - boundary_kernel(distance, smoothing_length(particle_system, particle)) + boundary_kernel_monaghan(distance, smoothing_length(particle_system, particle)) end -@fastpow @inline function boundary_kernel(r, h) +@fastpow @inline function boundary_kernel_monaghan(r, h) q = r / h # TODO The neighborhood search fluid->boundary should use this search distance diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index 1226d8f8bc..7a8cfac11d 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -41,6 +41,17 @@ function WallBoundarySystem(initial_condition, model; prescribed_motion=nothing, adhesion_coefficient=0.0, color_value=0) coordinates = copy(initial_condition.coordinates) + # Compute the normal vectors for the boundary model + # We need to do this here, because the normals are stored in `initial_condition` + if !isnothing(initial_condition.normals) + (; coordinates, normals) = initial_condition + coords_s = reinterpret(reshape, SVector{ndims(model), eltype(coordinates)}, + coordinates) + norms_s = reinterpret(reshape, SVector{ndims(model), eltype(normals)}, normals) + interpolation_coords = coords_s .- (2 .* norms_s) # The normals point out, i.e. from the fluid to the boundary + model.cache.interpolation_coords .= interpolation_coords + end + ismoving = Ref(!isnothing(prescribed_motion)) initialize_prescribed_motion!(prescribed_motion, initial_condition) diff --git a/src/schemes/boundary/wall_boundary/wall_boundary.jl b/src/schemes/boundary/wall_boundary/wall_boundary.jl index f20a5dbf20..fe7ac12729 100644 --- a/src/schemes/boundary/wall_boundary/wall_boundary.jl +++ b/src/schemes/boundary/wall_boundary/wall_boundary.jl @@ -1,4 +1,5 @@ include("dummy_particles.jl") include("system.jl") +include("marrone.jl") # Monaghan-Kajtar repulsive boundary particles require the `WallBoundarySystem` # and the `TotalLagrangianSPHSystem` and are therefore included later. diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index 78c429592d..55a367fad2 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -9,7 +9,7 @@ include("structure/discrete_element_method/discrete_element_method.jl") # Monaghan-Kajtar repulsive boundary particles require the `WallBoundarySystem` # and the `TotalLagrangianSPHSystem`. include("boundary/wall_boundary/monaghan_kajtar.jl") -# Implicit incompressible SPH requires the `BoundarySPHSystem` +# Implicit incompressible SPH requires the `WallBoundarySystem` include("fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl") # Include rhs for all schemes diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index ce420fbdf3..20fe630ff8 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -100,7 +100,8 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} boundary_density=fluid_density, n_layers=1, spacing_ratio=1, min_coordinates=zeros(length(fluid_size)), - faces=Tuple(trues(2 * length(fluid_size)))) + faces=Tuple(trues(2 * length(fluid_size))), + normal=false) NDIMS = length(fluid_size) ELTYPE = eltype(particle_spacing) fluid_size_ = Tuple(ELTYPE.(fluid_size)) @@ -136,10 +137,11 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} boundary_spacing = particle_spacing / spacing_ratio boundary_coordinates, - face_indices = initialize_boundaries(boundary_spacing, - tank_size_, - n_boundaries_per_dim, - n_layers, faces) + face_indices, + corner_indices = initialize_boundaries(boundary_spacing, + tank_size_, + n_boundaries_per_dim, + n_layers, faces) boundary_masses = boundary_density * boundary_spacing^NDIMS * ones(ELTYPE, size(boundary_coordinates, 2)) @@ -151,10 +153,14 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} particle_spacing, n_particles_per_dim) + normals = normal == false ? nothing : + compute_normals(boundary_coordinates, boundary_spacing, + face_indices, corner_indices, faces) + boundary = InitialCondition(coordinates=boundary_coordinates, velocity=boundary_velocities, mass=boundary_masses, density=boundary_densities, - particle_spacing=boundary_spacing) + particle_spacing=boundary_spacing, normals=normals) # Move the tank corner in the negative coordinate directions to the desired position boundary.coordinates .+= min_coordinates @@ -187,6 +193,156 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} end end +function compute_normals(boundary_coordinates, boundary_spacing, + face_indices, corner_indices, faces) + _compute_normals(boundary_coordinates, boundary_spacing, face_indices, corner_indices, + faces, + Val(size(boundary_coordinates, 1))) +end + +# 2D +function _compute_normals(boundary_coordinates, boundary_spacing, + face_indices, corner_indices, faces, ::Val{2}) + normals = zeros(size(boundary_coordinates)) + face_indices = Tuple(vec(x) for x in face_indices) + offset = boundary_spacing / 2 + + #### Left boundary + if faces[1] + left_boundary = maximum(boundary_coordinates[1, face_indices[1]]) + offset + for idx in face_indices[1] + normals[1, idx] = -abs(boundary_coordinates[1, idx] - left_boundary) + end + end + + #### Right boundary + if faces[2] + right_boundary = minimum(boundary_coordinates[1, face_indices[2]]) - offset + for idx in face_indices[2] + normals[1, idx] = abs(boundary_coordinates[1, idx] - right_boundary) + end + end + + #### Bottom boundary + if faces[3] + bottom_boundary = maximum(boundary_coordinates[2, face_indices[3]]) + offset + for idx in face_indices[3] + normals[2, idx] = -abs(boundary_coordinates[2, idx] - bottom_boundary) + end + end + + #### Top boundary + if faces[4] + top_boundary = minimum(boundary_coordinates[2, face_indices[4]]) - offset + for idx in face_indices[4] + normals[2, idx] = abs(boundary_coordinates[2, idx] - top_boundary) + end + end + + # Bottom left corner + if faces[1] && faces[3] + boundary_corner_point = [maximum(boundary_coordinates[1, corner_indices[1]]) + maximum(boundary_coordinates[2, corner_indices[1]])] + corner_point = boundary_corner_point + [offset; offset] + for idx in corner_indices[1] + normals[:, idx] = boundary_coordinates[:, idx] - corner_point + end + end + + # Top left corner + if faces[1] && faces[4] + boundary_corner_point = [maximum(boundary_coordinates[1, corner_indices[2]]) + minimum(boundary_coordinates[2, corner_indices[2]])] + corner_point = boundary_corner_point + [offset; -offset] + for idx in corner_indices[2] + normals[:, idx] = boundary_coordinates[:, idx] - corner_point + end + end + + # Bottom right corner + if faces[2] && faces[3] + boundary_corner_point = [minimum(boundary_coordinates[1, corner_indices[3]]) + maximum(boundary_coordinates[2, corner_indices[3]])] + corner_point = boundary_corner_point + [-offset; offset] + for idx in corner_indices[3] + normals[:, idx] = boundary_coordinates[:, idx] - corner_point + end + end + + # Top right corner + if faces[2] && faces[4] + boundary_corner_point = [minimum(boundary_coordinates[1, corner_indices[4]]) + minimum(boundary_coordinates[2, corner_indices[4]])] + corner_point = boundary_corner_point + [-offset; -offset] + for idx in corner_indices[4] + normals[:, idx] = boundary_coordinates[:, idx] - corner_point + end + end + + return normals +end + +# 3D +# Note: havent properly tested this yet +function _compute_normals(boundary, fluid, face_indices, corner_indices, faces, ::Val{3}) + (; coordinates) = boundary + normals = zeros(size(coordinates)) + face_indices = Tuple(vec(x) for x in face_indices) + offset = (boundary.particle_spacing + fluid.particle_spacing) / 2 + + #### +x boundary + if faces[1] + x_neg_boundary = maximum(coordinates[1, face_indices[1]]) + offset + for idx in face_indices[1] + normals[1, idx] = abs(coordinates[1, idx] - x_neg_boundary) + end + end + + #### -x boundary + if faces[2] + x_pos_boundary = minimum(coordinates[1, face_indices[2]]) - offset + for idx in face_indices[2] + normals[1, idx] = -abs(coordinates[1, idx] - x_pos_boundary) + end + end + + #### +y boundary + if faces[3] + y_neg_boundary = maximum(coordinates[2, face_indices[3]]) + offset + for idx in face_indices[3] + normals[2, idx] = abs(coordinates[2, idx] - y_neg_boundary) + end + end + + #### -y boundary + if faces[4] + y_pos_boundary = minimum(coordinates[2, face_indices[4]]) - offset + for idx in face_indices[4] + normals[2, idx] = -abs(coordinates[2, idx] - y_pos_boundary) + end + end + + #### +z boundary + if faces[5] + z_neg_boundary = maximum(coordinates[3, face_indices[5]]) + offset + for idx in face_indices[5] + normals[3, idx] = abs(coordinates[3, idx] - z_neg_boundary) + end + end + + #### -z boundary + if faces[6] + z_pos_boundary = minimum(coordinates[3, face_indices[6]]) - offset + for idx in face_indices[6] + normals[3, idx] = -abs(coordinates[3, idx] - z_pos_boundary) + end + end + + # TODO: edges + + return normals +end + function round_n_particles(size, spacing, type) n_particles = round(Int, size / spacing) @@ -320,6 +476,10 @@ function initialize_boundaries(particle_spacing, tank_size::NTuple{2}, face_indices_2 = Array{Int, 2}(undef, n_layers, n_particles_y) face_indices_3 = Array{Int, 2}(undef, n_layers, n_particles_x) face_indices_4 = Array{Int, 2}(undef, n_layers, n_particles_x) + corner_indices_1 = Array{Int, 2}(undef, n_layers, n_layers) + corner_indices_2 = Array{Int, 2}(undef, n_layers, n_layers) + corner_indices_3 = Array{Int, 2}(undef, n_layers, n_layers) + corner_indices_4 = Array{Int, 2}(undef, n_layers, n_layers) # Create empty array to extend later depending on faces and corners to build boundary_coordinates = Array{typeof(particle_spacing), 2}(undef, 2, 0) @@ -409,6 +569,13 @@ function initialize_boundaries(particle_spacing, tank_size::NTuple{2}, (n_layers, n_layers), (layer_offset, layer_offset)) boundary_coordinates = hcat(boundary_coordinates, bottom_left_corner) + + # store the indices of each particle + particles_per_layer = n_layers + for i in 1:n_layers + corner_indices_1[i, :] = collect((index + 1):(particles_per_layer + index)) + index += particles_per_layer + end end # Top left @@ -417,6 +584,13 @@ function initialize_boundaries(particle_spacing, tank_size::NTuple{2}, (n_layers, n_layers), (layer_offset, tank_size[2])) boundary_coordinates = hcat(boundary_coordinates, top_left_corner) + + # store the indices of each particle + particles_per_layer = n_layers + for i in 1:n_layers + corner_indices_2[i, :] = collect((index + 1):(particles_per_layer + index)) + index += particles_per_layer + end end # Bottom right @@ -425,6 +599,13 @@ function initialize_boundaries(particle_spacing, tank_size::NTuple{2}, (n_layers, n_layers), (tank_size[1], layer_offset)) boundary_coordinates = hcat(boundary_coordinates, bottom_right_corner) + + # store the indices of each particle + particles_per_layer = n_layers + for i in 1:n_layers + corner_indices_3[i, :] = collect((index + 1):(particles_per_layer + index)) + index += particles_per_layer + end end # Top right @@ -433,10 +614,18 @@ function initialize_boundaries(particle_spacing, tank_size::NTuple{2}, (n_layers, n_layers), (tank_size[1], tank_size[2])) boundary_coordinates = hcat(boundary_coordinates, top_right_corner) + + # store the indices of each particle + particles_per_layer = n_layers + for i in 1:n_layers + corner_indices_4[i, :] = collect((index + 1):(particles_per_layer + index)) + index += particles_per_layer + end end return boundary_coordinates, - (face_indices_1, face_indices_2, face_indices_3, face_indices_4) + (face_indices_1, face_indices_2, face_indices_3, face_indices_4), + (corner_indices_1, corner_indices_2, corner_indices_3, corner_indices_4) end # 3D diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index 3c6a7466d5..5668413809 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -91,7 +91,8 @@ SphereShape(0.1, 0.5, (0.2, 0.4, 0.3), 1000.0, sphere_type=RoundSphere()) function SphereShape(particle_spacing, radius, center_position, density; sphere_type=VoxelSphere(), n_layers=-1, layer_outwards=false, cutout_min=(0.0, 0.0), cutout_max=(0.0, 0.0), place_on_shell=false, - velocity=zeros(length(center_position)), mass=nothing, pressure=0) + velocity=zeros(length(center_position)), mass=nothing, pressure=0, + normal=false) if particle_spacing < eps() throw(ArgumentError("`particle_spacing` needs to be positive and larger than $(eps())")) end @@ -117,8 +118,12 @@ function SphereShape(particle_spacing, radius, center_position, density; particles_not_in_cutout = map(!in_cutout, axes(coordinates, 2)) coordinates = coordinates[:, particles_not_in_cutout] + if normal + normals = compute_normals(coordinates, center_position, radius) + end + return InitialCondition(; coordinates, velocity, mass, density, pressure, - particle_spacing) + particle_spacing, normals) end """ @@ -423,3 +428,47 @@ function round_sphere(sphere, particle_spacing, radius, center::SVector{3}) return particle_coords end + +# Copy to REPL and run +# function plot_coords(boundary::Matrix{T}, normals=nothing) where {T} +# if size(boundary)[1] == 2 +# x_b, y_b = eachrow(boundary) + +# Plots.plot(x_b, y_b, seriestype=:scatter, color=:blue, label="Boundary") + +# if normals !== nothing +# u, v = eachrow(normals) +# Plots.quiver!(x_b, y_b, quiver=(u, v), aspect_ratio=1, label="Normals") +# end + +# elseif size(boundary)[1] == 3 +# x_b, y_b, z_b = eachrow(boundary) + +# Plots.plot(x_b, y_b, z_b, seriestype=:scatter, color=:blue, label="Boundary") +# end +# end + +function compute_normals(coordinates::Matrix{T}, center_position, + radius::T) where {T} + normals = zeros(size(coordinates)) + center_position = collect(center_position) + + for idx in 1:size(coordinates, 2) + coord = coordinates[:, idx] + + # Project the point at coord on the circle + diff = center_position - coord + + # Check for division-by-zero + diff_norm = norm(diff) + if iszero(diff_norm) + normal = zeros(T, 2) + else + proj = center_position + radius * (diff / diff_norm) + normal = proj - coord + end + normals[:, idx] = normal + end + + return normals +end diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 4f46abe701..d6f3f0e3d8 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -106,8 +106,7 @@ volume = boundary_model.cache.volume # Reset cache and perform pressure extrapolation - TrixiParticles.reset_cache!(boundary_model.cache, - boundary_model.viscosity) + TrixiParticles.reset_cache!(boundary_model) TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, boundary_system, @@ -158,8 +157,7 @@ end # Reset cache and perform pressure extrapolation - TrixiParticles.reset_cache!(boundary_model.cache, - boundary_model.viscosity) + TrixiParticles.reset_cache!(boundary_model) TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, boundary_system, @@ -252,8 +250,7 @@ TrixiParticles.compute_pressure!(fluid_system1, v_fluid, semi) TrixiParticles.set_zero!(boundary_model.pressure) - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - viscosity) + TrixiParticles.reset_cache!(boundary_model) TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, boundary_system, @@ -291,8 +288,7 @@ TrixiParticles.compute_pressure!(fluid_system2, v_fluid, semi) TrixiParticles.set_zero!(boundary_model.pressure) - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - viscosity) + TrixiParticles.reset_cache!(boundary_model) TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, boundary_system, @@ -334,8 +330,7 @@ TrixiParticles.compute_pressure!(fluid_system3, v_fluid, semi) TrixiParticles.set_zero!(boundary_model.pressure) - TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, - viscosity) + TrixiParticles.reset_cache!(boundary_model) TrixiParticles.boundary_pressure_extrapolation!(flipped_condition, boundary_model, diff --git a/test/schemes/boundary/marrone/marrone.jl b/test/schemes/boundary/marrone/marrone.jl new file mode 100644 index 0000000000..5de6022d25 --- /dev/null +++ b/test/schemes/boundary/marrone/marrone.jl @@ -0,0 +1,447 @@ +@testset verbose=true "Dummy Particles with `MarronePressureExtrapolation`" begin + @testset "Compute Boundary Normal Vectors" begin + particle_spacing = 1.0 + n_particles = 2 + n_layers = 1 + width = particle_spacing * n_particles + height = particle_spacing * n_particles + density = 257 + + tank = RectangularTank(particle_spacing, (width, height), (width, height), + density, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + + (; normals) = tank.boundary + normals_reference = [[-0.5 -0.5 0.5 0.5 0.0 0.0 -0.5 0.5] + [0.0 0.0 0.0 0.0 -0.5 -0.5 -0.5 -0.5]] + + @test normals == normals_reference + end + + @testset "`MarroneMLSKernel`" begin + particle_spacing = 1.0 + n_particles = 2 + n_layers = 1 + width = particle_spacing * n_particles + height = particle_spacing * n_particles + density = 257 + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = particle_spacing * 1.1 + state_equation = StateEquationCole(sound_speed=10, reference_density=257, + exponent=7) + + tank1 = RectangularTank(particle_spacing, (width, height), (width, height), + density, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + n_boundary_particles = size(tank1.boundary.coordinates, 2) + n_fluid_particles = size(tank1.fluid.coordinates, 2) + + mls_kernel = MarroneMLSKernel(smoothing_kernel, n_boundary_particles, + n_fluid_particles) + boundary_model = BoundaryModelDummyParticles(tank1.boundary.density, + tank1.boundary.mass, + state_equation=state_equation, + MarronePressureExtrapolation(), + mls_kernel, smoothing_length) + boundary_system = WallBoundarySystem(tank1.boundary, boundary_model) + viscosity = boundary_system.boundary_model.viscosity + semi = DummySemidiscretization() + fluid_system1 = WeaklyCompressibleSPHSystem(tank1.fluid, SummationDensity(), + state_equation, + smoothing_kernel, smoothing_length) + fluid_system1.cache.density .= tank1.fluid.density + v_fluid = zeros(2, TrixiParticles.nparticles(fluid_system1)) + + TrixiParticles.compute_pressure!(fluid_system1, v_fluid, semi) + + TrixiParticles.set_zero!(boundary_model.pressure) + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + viscosity) + + boundary_coords = tank1.boundary.coordinates + fluid_coords = tank1.fluid.coordinates + + expected_basis = zeros(n_boundary_particles, n_fluid_particles, 3) + @testset "Compute Marrone Basis" begin + (; basis) = mls_kernel + expected_basis[:, 1, + :] = [1.0 1.0 0.0 + 1.0 1.0 -1.0 + 1.0 -2.0 0.0 + 0.0 0.0 0.0 + 1.0 0.0 1.0 + 1.0 -1.0 1.0 + 1.0 1.0 1.0 + 0.0 0.0 0.0] + expected_basis[:, 2, + :] = [1.0 2.0 0.0 + 0.0 0.0 0.0 + 1.0 -1.0 0.0 + 1.0 -1.0 -1.0 + 1.0 1.0 1.0 + 1.0 0.0 1.0 + 0.0 0.0 0.0 + 1.0 -1.0 1.0] + expected_basis[:, 3, + :] = [1.0 1.0 1.0 + 1.0 1.0 0.0 + 0.0 0.0 0.0 + 1.0 -2.0 0.0 + 1.0 0.0 2.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0] + expected_basis[:, 4, + :] = [0.0 0.0 0.0 + 1.0 2.0 0.0 + 1.0 -1.0 1.0 + 1.0 -1.0 0.0 + 0.0 0.0 0.0 + 1.0 0.0 2.0 + 0.0 0.0 0.0 + 0.0 0.0 0.0] + + TrixiParticles.compute_basis_marrone(mls_kernel, boundary_system, fluid_system1, + boundary_coords, fluid_coords, semi) + @test basis == expected_basis + end + + expected_momentum = zeros(n_boundary_particles, 3, 3) + @testset "Compute Marrone Momentum" begin + (; inner_kernel, momentum) = mls_kernel + + # The momentum is computed as M_i = ∑ⱼ b_ij ⊗ b_ij' ⋅ kernel_weight_j ⋅ volume_j + # The volume of all fluid particles is 1.0, so we leave this factor out + # We need only consider the neighbors in the compact support of the inner kernel + for (i, fluid_indices) in Dict(1 => [1, 2, 3], 2 => [1, 3, 4], 5 => [1, 2, 3]) + for j in fluid_indices + diff = boundary_coords[:, i] - fluid_coords[:, j] + distance = sqrt(dot(diff, diff)) + expected_momentum[i, :, + :] += expected_basis[i, j, :] .* + expected_basis[i, j, :]' * + TrixiParticles.kernel(inner_kernel, distance, + smoothing_length) + end + + if abs(det(expected_momentum[i, :, :])) < 1.0f-9 + expected_momentum[i, :, + :] = Matrix{Float64}(I, 3, 3) + else + expected_momentum[i, :, :] = inv(expected_momentum[i, :, :]) + end + end + + TrixiParticles.compute_momentum_marrone(mls_kernel, boundary_system, + fluid_system1, boundary_coords, + fluid_coords, v_fluid, semi, + smoothing_length) + + @test expected_momentum[1, :, :] == momentum[1, :, :] + @test expected_momentum[2, :, :] == momentum[2, :, :] + @test expected_momentum[5, :, :] == momentum[5, :, :] + @test Matrix{Float64}(I, 3, 3) == momentum[7, :, :] + end + + @testset "Pressure Extrapolation" begin + particle_spacing = 1.0 + n_particles = 10 + n_layers = 4 + width = particle_spacing * n_particles + height = particle_spacing * n_particles + density = 257 + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 3 * particle_spacing + state_equation = StateEquationCole(sound_speed=10, reference_density=257, + exponent=7) + + tank1 = RectangularTank(particle_spacing, (width, height), (width, height), + density, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + n_boundary_particles = size(tank1.boundary.coordinates, 2) + n_fluid_particles = size(tank1.fluid.coordinates, 2) + + mls_kernel = MarroneMLSKernel(smoothing_kernel, n_boundary_particles, + n_fluid_particles) + + boundary_model = BoundaryModelDummyParticles(tank1.boundary.density, + tank1.boundary.mass, + state_equation=state_equation, + MarronePressureExtrapolation(), + mls_kernel, smoothing_length) + + boundary_system = WallBoundarySystem(tank1.boundary, boundary_model) + viscosity = boundary_system.boundary_model.viscosity + + semi = DummySemidiscretization() + + # In this testset, we verify that pressures in a static tank are extrapolated correctly. + # Use constant density equal to the reference density of the state equation, + # so that the pressure is constant zero. Then we test that the extrapolation also yields zero. + @testset "Constant Zero Pressure" begin + fluid_system1 = WeaklyCompressibleSPHSystem(tank1.fluid, SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length) + fluid_system1.cache.density .= tank1.fluid.density + v_fluid = zeros(2, TrixiParticles.nparticles(fluid_system1)) + + TrixiParticles.compute_pressure!(fluid_system1, v_fluid, semi) + + TrixiParticles.set_zero!(boundary_model.pressure) + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + viscosity) + + TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, + boundary_system, + fluid_system1, + tank1.boundary.coordinates, + tank1.fluid.coordinates, + v_fluid, + v_fluid, + semi) + + @test all(boundary_system.boundary_model.pressure .== 0.0) + @test all(fluid_system1.pressure .== 0.0) + end + + # Test whether the pressure is constant if the density of the state equation + # and in the tank are not the same. + # Then we test that the extrapolation yields some constant value. + @testset "Constant Non-Zero Pressure" begin + density = 260 + tank2 = RectangularTank(particle_spacing, (width, height), (width, height), + density, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + + fluid_system2 = WeaklyCompressibleSPHSystem(tank2.fluid, SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length) + + fluid_system2.cache.density .= tank2.fluid.density + v_fluid = zeros(2, TrixiParticles.nparticles(fluid_system2)) + TrixiParticles.compute_pressure!(fluid_system2, v_fluid, semi) + + TrixiParticles.set_zero!(boundary_model.pressure) + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + viscosity) + + TrixiParticles.boundary_pressure_extrapolation!(Val(true), boundary_model, + boundary_system, + fluid_system2, + tank2.boundary.coordinates, + tank2.fluid.coordinates, + v_fluid, + v_fluid, + semi) + + # Test that pressure of the fluid is indeed constant + @test all(isapprox.(fluid_system2.pressure, fluid_system2.pressure[1])) + # Test that boundary pressure equals fluid pressure + @test all(isapprox.(boundary_system.boundary_model.pressure, + fluid_system2.pressure[1], atol=1.0e-10)) + end + + # In this test, we initialize a fluid with a hydrostatic pressure gradient + # and check that this gradient is extrapolated correctly. + @testset "Hydrostatic Pressure Gradient" begin + particle_spacing = 0.1 + n_particles = 10 + n_layers = 2 + width = particle_spacing * n_particles + height = particle_spacing * n_particles + density = 257 + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 2 * particle_spacing + state_equation = StateEquationCole(sound_speed=10, reference_density=257, + exponent=7) + + tank1 = RectangularTank(particle_spacing, (width, height), (width, height), + density, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + n_boundary_particles = size(tank1.boundary.coordinates, 2) + n_fluid_particles = size(tank1.fluid.coordinates, 2) + + mls_kernel = MarroneMLSKernel(smoothing_kernel, n_boundary_particles, + n_fluid_particles) + + boundary_model = BoundaryModelDummyParticles(tank1.boundary.density, + tank1.boundary.mass, + state_equation=state_equation, + MarronePressureExtrapolation(), + mls_kernel, smoothing_length) + + boundary_system = WallBoundarySystem(tank1.boundary, boundary_model) + viscosity = boundary_system.boundary_model.viscosity + + semi = DummySemidiscretization() + + tank3 = RectangularTank(particle_spacing, (width, height), (width, height), + density, acceleration=[0.0, -9.81], + state_equation=state_equation, n_layers=n_layers, + faces=(true, true, true, false), normal=true) + + fluid_system3 = WeaklyCompressibleSPHSystem(tank3.fluid, SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length, + acceleration=[0.0, -9.81]) + fluid_system3.cache.density .= tank3.fluid.density + v_fluid = zeros(2, TrixiParticles.nparticles(fluid_system3)) + TrixiParticles.compute_pressure!(fluid_system3, v_fluid, semi) + + TrixiParticles.set_zero!(boundary_model.pressure) + TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, + viscosity) + + TrixiParticles.boundary_pressure_extrapolation!(Val(true), + boundary_model, + boundary_system, + fluid_system3, + tank3.boundary.coordinates, + tank3.fluid.coordinates, + v_fluid, + v_fluid, + semi) + + width_reference = particle_spacing * (n_particles + 2 * n_layers) + height_reference = particle_spacing * (n_particles + n_layers) + + # Define another tank without a boundary, where the fluid has the same size + # as fluid plus boundary in the other tank. + # The pressure gradient of this fluid should be the same as the extrapolated pressure + # of the boundary in the first tank. + tank_reference = RectangularTank(particle_spacing, + (width_reference, height_reference), + (width_reference, height_reference), + density, acceleration=[0.0, -9.81], + state_equation=state_equation, n_layers=0, + faces=(true, true, true, false)) + + # Because it is a pain to deal with the linear indices of the pressure arrays, + # we convert the matrices to Cartesian indices based on the coordinates. + function set_pressure!(pressure, coordinates, offset, system, + system_pressure) + for particle in TrixiParticles.eachparticle(system) + # Coordinates as integer indices + coords = coordinates[:, particle] ./ particle_spacing + # Move bottom left corner to (1, 1) + coords .+= offset + # Round to integer index + index = round.(Int, coords) + pressure[index...] = system_pressure[particle] + end + end + + # Set up the combined pressure matrix to store the pressure values of fluid + # and boundary together. + pressure = zeros(n_particles + 2 * n_layers, n_particles + n_layers) + + # The fluid starts at -0.5 * particle_spacing from (0, 0), + # so the boundary starts at -(n_layers + 0.5) * particle_spacing + set_pressure!(pressure, boundary_system.coordinates, n_layers + 0.5, + boundary_system, boundary_system.boundary_model.pressure) + + # The fluid starts at -0.5 * particle_spacing from (0, 0), + # so the boundary starts at -(n_layers + 0.5) * particle_spacing + set_pressure!(pressure, tank3.fluid.coordinates, n_layers + 0.5, + fluid_system3, fluid_system3.pressure) + pressure_reference = similar(pressure) + + # The fluid starts at -0.5 * particle_spacing from (0, 0) + set_pressure!(pressure_reference, tank_reference.fluid.coordinates, 0.5, + tank_reference.fluid, tank_reference.fluid.pressure) + + # TODO: test failing, maximum difference between approximation and correct solution is + # maximum(abs.(pressure - pressure-reference)) -> 813.51 + # @test all(isapprox.(pressure, pressure_reference, atol=4.0)) + end + + @testset "Numerical Consistency" begin + mls_kernel = MarroneMLSKernel(smoothing_kernel, n_fluid_particles, + n_fluid_particles) + fluid_coords = tank1.fluid.coordinates + fluid_system1 = WeaklyCompressibleSPHSystem(tank1.fluid, SummationDensity(), + state_equation, + smoothing_kernel, + smoothing_length) + fluid_system1.cache.density .= tank1.fluid.density + + TrixiParticles.compute_basis_marrone(mls_kernel, fluid_system1, + fluid_system1, fluid_coords, + fluid_coords, semi) + TrixiParticles.compute_momentum_marrone(mls_kernel, fluid_system1, + fluid_system1, + fluid_coords, + fluid_coords, v_fluid, semi, + smoothing_length) + + # We test that the `MarroneMLSKernel` correctly computes the + # first derivative of a constant function. + @testset "Zeroth Order Consistency" begin + zero_order_approx = zeros(n_fluid_particles) + constant = 3.0 + TrixiParticles.foreach_point_neighbor(fluid_system1, fluid_system1, + fluid_coords, fluid_coords, + semi) do particle, neighbor, + pos_diff, distance + neighbor_density = TrixiParticles.current_density(v_fluid, + fluid_system1, + neighbor) + neighbor_volume = neighbor_density != 0 ? + TrixiParticles.hydrodynamic_mass(fluid_system1, + neighbor) / + neighbor_density : 0 + + zero_order_approx[particle] += TrixiParticles.boundary_kernel_marrone(mls_kernel, + particle, + neighbor, + distance, + smoothing_length) * + constant * + neighbor_volume + end + @test all(isapprox.(zero_order_approx, constant, atol=1.0e-10)) + end + + # We test that the `MarroneMLSKernel` correctly computes the + # first derivative of a linear function. + @testset "First Order Consistency" begin + first_order_approx = zeros(n_fluid_particles) + a = [1, 2] + b = 3 + f(x) = dot(a, x) + b + + linear_mapping = [f(fluid_coords[:, particle]) + for particle in 1:n_fluid_particles] + TrixiParticles.foreach_point_neighbor(fluid_system1, fluid_system1, + fluid_coords, fluid_coords, + semi) do particle, neighbor, + pos_diff, distance + neighbor_density = TrixiParticles.current_density(v_fluid, + fluid_system1, + neighbor) + neighbor_volume = neighbor_density != 0 ? + TrixiParticles.hydrodynamic_mass(fluid_system1, + neighbor) / + neighbor_density : 0 + neighbor_val = f(fluid_coords[:, neighbor]) + + first_order_approx[particle] += TrixiParticles.boundary_kernel_marrone(mls_kernel, + particle, + neighbor, + distance, + smoothing_length) * + neighbor_val * + neighbor_volume + end + @test all(isapprox.(first_order_approx, linear_mapping, atol=1.0e-10)) + end + end + end + end +end diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index bb6e9d5483..baab34bf41 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -1,5 +1,6 @@ include("structure/total_lagrangian_sph/total_lagrangian_sph.jl") include("boundary/dummy_particles/dummy_particles.jl") +include("boundary/marrone/marrone.jl") include("boundary/monaghan_kajtar/monaghan_kajtar.jl") include("boundary/open_boundary/open_boundary.jl") include("boundary/prescribed_motion.jl")