diff --git a/NEWS.md b/NEWS.md index ffe272e8a..5d7a295d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,8 @@ # Changelog TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) -used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. +used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. ## Version 0.2.x @@ -14,11 +14,15 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between ### Deprecated +## Version 0.1.3 + +### Added +Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC. ## Version 0.1.2 ### Added -A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH +A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH. ## Version 0.1.1 diff --git a/README.md b/README.md index b40e68fbc..76b288ed2 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ Its features include: ## Features - Incompressible Navier-Stokes - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC) - - Models: Surface Tension + - Models: Surface Tension, Open Boundaries - Solid-body mechanics - Methods: Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM) - Fluid-Structure Interaction diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 8b374a71e..ae4b9b676 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -234,3 +234,86 @@ Pages = [joinpath("schemes", "boundary", "monaghan_kajtar", "monaghan_kajtar.jl" - Alireza Valizadeh, Joseph J. Monaghan. "A study of solid wall models for weakly compressible SPH." In: Journal of Computational Physics 300 (2015), pages 5–19. [doi: 10.1016/J.JCP.2015.07.033](https://doi.org/10.1016/J.JCP.2015.07.033) + +# [Open Boundaries](@id open_boundary) + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "boundary", "open_boundary", "system.jl")] +``` + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "boundary", "open_boundary", "boundary_zones.jl")] +``` + +### [Method of characteristics](@id method_of_characteristics) + +The difficulty in non-reflecting boundary conditions, also called open boundaries, is to determine +the appropriate boundary values of the exact characteristics of the Euler equations. +Assuming the flow near the boundaries is normal to the boundary +and free of shock waves and significant viscous effects, it can be shown that three characteristic variables exist: + +- ``J_1``, associated with convection of entropy and propagates at flow velocity, +- ``J_2``, downstream-running characteristics, +- ``J_3``, upstream-running characteristics. + +Giles (1990) derived those variables based on a linearized set of governing equations: +```math +J_1 = -c_s^2 (\rho - \rho_{\text{ref}}) + (p - p_{\text{ref}}) +``` +```math +J_2 = \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}}) +``` +```math +J_3 = - \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}}) +``` +where the subscript "ref" denotes the reference flow near the boundaries, which can be prescribed. + +Specifying the reference variables is **not** equivalent to prescription of ``\rho``, ``v`` and ``p`` +directly, since the perturbation from the reference flow is allowed. + +Lastiwka et al. (2009) applied the method of characteristic to SPH and determined the number of variables that should be +**prescribed** at the boundary and the number which should be **propagated** from the fluid domain to the boundary: + +- For an **inflow** boundary: + - Prescribe *downstream*-running characteristics ``J_1`` and ``J_2`` + - Transmit ``J_3`` from the fluid domain (allow ``J_3`` to propagate upstream to the boundary). + +- For an **outflow** boundary: + - Prescribe *upstream*-running characteristic ``J_3`` + - Transmit ``J_1`` and ``J_2`` from the fluid domain. + +Prescribing is done by simply setting the characteristics to zero. To transmit the characteristics from the fluid +domain, or in other words, to carry the information of the fluid to the boundaries, Negi et al. (2020) use a Shepard Interpolation +```math +f_i = \frac{\sum_j^N f_j W_{ij}}{\sum_j^N W_{ij}}, +``` +where the ``i``-th particle is a boundary particle, ``f`` is either ``J_1``, ``J_2`` or ``J_3`` and ``N`` is the set of +neighboring fluid particles. + +To express pressure ``p``, density ``\rho`` and velocity ``v`` as functions of the characteristic variables, the system of equations +from the characteristic variables is inverted and gives +```math + \rho - \rho_{\text{ref}} = \frac{1}{c_s^2} \left( -J_1 + \frac{1}{2} J_2 + \frac{1}{2} J_3 \right), +``` +```math +u - u_{\text{ref}}= \frac{1}{2\rho c_s} \left( J_2 - J_3 \right), +``` +```math +p - p_{\text{ref}} = \frac{1}{2} \left( J_2 + J_3 \right). +``` +With ``J_1``, ``J_2`` and ``J_3`` determined, we can easily solve for the actual variables for each particle. + +### References +- M. B. Giles. "Nonreflecting boundary conditions for Euler equation calculations". + In: AIAA Journal, 28.12 pages 2050--2058. + [doi: 10.2514/3.10521](https://doi.org/10.2514/3.10521) +- M. Lastiwka, M. Basa, N. J. Quinlan. + "Permeable and non-reflecting boundary conditions in SPH". + In: International Journal for Numerical Methods in Fluids 61, (2009), pages 709--724. + [doi: 10.1002/fld.1971](https://doi.org/10.1002/fld.1971) +- P. Negi, P. Ramachandran, A. Haftu. + "An improved non-reflecting outlet boundary condition for weakly-compressible SPH". + In: Computer Methods in Applied Mechanics and Engineering 367, (2020), pages 113--119. + [doi: 10.1016/j.cma.2020.113119](https://doi.org/10.1016/j.cma.2020.113119) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl new file mode 100644 index 000000000..3d642696e --- /dev/null +++ b/examples/fluid/pipe_flow_2d.jl @@ -0,0 +1,131 @@ +# 2D channel flow simulation with open boundaries. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.4) + +flow_direction = [1.0, 0.0] +reynolds_number = 100 +const prescribed_velocity = 2.0 + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2]) + +fluid_density = 1000.0 + +# For this particular example, it is necessary to have a background pressure. +# Otherwise the suction at the outflow is to big and the simulation becomes unstable. +pressure = 1000.0 + +sound_speed = 10 * prescribed_velocity + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=pressure) + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + faces=(false, false, true, true)) + +# Shift pipe walls in negative x-direction for the inflow +pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers + +n_buffer_particles = 4 * pipe.n_particles_per_dimension[2] + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + +# Alternatively the WCSPH scheme can be used +# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +function velocity_function(pos, t) + # Use this for a time-dependent inflow velocity + # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) + + return SVector(prescribed_velocity, 0.0) +end + +inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing) + +open_boundary_in = OpenBoundarySPHSystem(inflow; sound_speed, fluid_system, + buffer_size=n_buffer_particles, + reference_pressure=pressure, + reference_velocity=velocity_function) + +outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), + flow_direction, open_boundary_layers, density=fluid_density, + particle_spacing) + +open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system, + buffer_size=n_buffer_particles, + reference_pressure=pressure, + reference_velocity=velocity_function) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, + boundary_system) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration) + dtmax=1e-2, # Limit stepsize to prevent crashing + save_everystep=false, callback=callbacks); diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index caf40a62e..c7226357a 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -5,12 +5,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothin initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] G :: ELTYPE + buffer :: Nothing function NBodySystem(initial_condition, G) mass = copy(initial_condition.mass) new{size(initial_condition.coordinates, 1), - eltype(mass)}(initial_condition, mass, G) + eltype(mass)}(initial_condition, mass, G, nothing) end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 2fe865d10..aa97ab720 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -44,7 +44,8 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, DEMSystem, BoundaryDEMSystem + BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, + OutFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback export ContinuityDensity, SummationDensity diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 5c8a7b5e1..74c6fe5c9 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -74,7 +74,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr # Update systems to compute quantities like density and pressure. semi = integrator.p v_ode, u_ode = u.x - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Apply the callback. cb(integrator) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index d68168258..9c39f3eb4 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -236,7 +236,7 @@ function (pp::PostprocessCallback)(integrator) new_data = false # Update systems to compute quantities like density and pressure - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) foreach_system(semi) do system if system isa BoundarySystem && pp.exclude_boundary diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 04552d57e..468a72a1c 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -139,7 +139,7 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in # Update systems to compute quantities like density and pressure semi = integrator.p v_ode, u_ode = u.x - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Apply the callback solution_callback(integrator) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index ca3b5f114..3595d5c36 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -52,7 +52,16 @@ function initial_update!(cb, u, t, integrator) initial_update!(cb.affect!, u, t, integrator) end -initial_update!(cb::UpdateCallback, u, t, integrator) = cb(integrator) +function initial_update!(cb::UpdateCallback, u, t, integrator) + semi = integrator.p + + # Tell systems that `UpdateCallback` is used + foreach_system(semi) do system + update_callback_used!(system) + end + + return cb(integrator) +end # `condition` function (update_callback!::UpdateCallback)(u, t, integrator) @@ -69,16 +78,12 @@ function (update_callback!::UpdateCallback)(integrator) # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Other updates might be added here later (e.g. Transport Velocity Formulation). - # @trixi_timeit timer() "update open boundary" foreach_system(semi) do system - # update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) - # end - # - # @trixi_timeit timer() "update TVF" foreach_system(semi) do system - # update_transport_velocity_eachstep!(system, v_ode, u_ode, semi, t) - # end + @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + end # Tell OrdinaryDiffEq that `u` has been modified u_modified!(integrator, true) @@ -128,3 +133,5 @@ function Base.show(io::IO, ::MIME"text/plain", summary_box(io, "UpdateCallback", setup) end end + +update_callback_used!(system) = system diff --git a/src/general/buffer.jl b/src/general/buffer.jl new file mode 100644 index 000000000..b12cf05f3 --- /dev/null +++ b/src/general/buffer.jl @@ -0,0 +1,87 @@ +struct SystemBuffer{V} + active_particle :: BitVector + eachparticle :: V # Vector{Int} + buffer_size :: Int + + function SystemBuffer(active_size, buffer_size::Integer) + active_particle = vcat(trues(active_size), falses(buffer_size)) + eachparticle = collect(1:active_size) + + return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size) + end +end + +allocate_buffer(initial_condition, buffer) = initial_condition + +function allocate_buffer(initial_condition, buffer::SystemBuffer) + (; buffer_size) = buffer + + # Initialize particles far away from simulation domain + coordinates = fill(1e16, ndims(initial_condition), buffer_size) + + if all(rho -> isapprox(rho, first(initial_condition.density), atol=eps(), rtol=eps()), + initial_condition.density) + density = first(initial_condition.density) + else + throw(ArgumentError("`initial_condition.density` needs to be constant when using `SystemBuffer`")) + end + + particle_spacing = initial_condition.particle_spacing + + buffer_ic = InitialCondition(; coordinates, density, particle_spacing) + + return union(initial_condition, buffer_ic) +end + +@inline update_system_buffer!(buffer::Nothing) = buffer + +# TODO `resize` allocates. Find a non-allocating version +@inline function update_system_buffer!(buffer::SystemBuffer) + (; active_particle) = buffer + + resize!(buffer.eachparticle, count(active_particle)) + + i = 1 + for j in eachindex(active_particle) + if active_particle[j] + buffer.eachparticle[i] = j + i += 1 + end + end + + return buffer +end + +@inline each_moving_particle(system, buffer) = buffer.eachparticle + +@inline active_coordinates(u, system, buffer) = view(u, :, buffer.active_particle) + +@inline active_particles(system, buffer) = buffer.eachparticle + +@inline function activate_next_particle(system) + (; active_particle) = system.buffer + + next_particle = findfirst(x -> !x, active_particle) + + if isnothing(next_particle) + error("0 out of $(system.buffer.buffer_size) buffer particles available") + end + + active_particle[next_particle] = true + + return next_particle +end + +@inline function deactivate_particle!(system, particle, u) + (; active_particle) = system.buffer + + active_particle[particle] = false + + # Set particle far away from simulation domain + for dim in 1:ndims(system) + # Inf or NaN causes instability outcome. + u[dim, particle] = 1e16 + end + + return system +end diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index ee89c92a6..7244253b0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -35,6 +35,19 @@ end return v[end, particle] end +# WARNING! +# These functions are intended to be used internally to set the density +# of newly activated particles in a callback. +# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` +# outside of callbacks. +@inline set_particle_density!(v, system, ::SummationDensity, particle, density) = v + +@inline function set_particle_density!(v, system, ::ContinuityDensity, particle, density) + v[end, particle] = density + + return v +end + function summation_density!(system, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) diff --git a/src/general/general.jl b/src/general/general.jl index 4a60cb52f..774d076c9 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -30,6 +30,7 @@ include("density_calculators.jl") include("corrections.jl") include("smoothing_kernels.jl") include("initial_condition.jl") +include("buffer.jl") include("system.jl") include("interpolation.jl") include("custom_quantities.jl") diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6a1d703b2..73c1fe56d 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -153,7 +153,19 @@ end return compact_support(smoothing_kernel, smoothing_length) end +@inline function compact_support(system::OpenBoundarySPHSystem, neighbor) + # Use the compact support of the fluid + return compact_support(neighbor, system) +end + +@inline function compact_support(system::OpenBoundarySPHSystem, + neighbor::OpenBoundarySPHSystem) + # This NHS is never used + return 0.0 +end + @inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem) + # This NHS is never used return 0.0 end @@ -281,6 +293,9 @@ function semidiscretize(semi, tspan; reset_threads=true, data_type=nothing) # Initialize this system initialize!(system, neighborhood_search) + + # Only for systems requiring a mandatory callback + reset_callback_flag!(system) end end @@ -343,6 +358,9 @@ function restart_with!(semi, sol; reset_threads=true) u = wrap_u(sol.u[end].x[2], system, semi) restart_with!(system, v, u) + + # Only for systems requiring a mandatory callback + reset_callback_flag!(system) end return semi @@ -429,7 +447,7 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) end # Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces -function update_systems_and_nhs(v_ode, u_ode, semi, t) +function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false) # First update step before updating the NHS # (for example for writing the current coordinates in the solid system) foreach_system(semi) do system @@ -466,7 +484,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_final!(system, v, u, v_ode, u_ode, semi, t) + update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) end end @@ -631,6 +649,32 @@ function update_nhs!(neighborhood_search, particles_moving=(true, neighbor.ismoving[])) end +function update_nhs!(neighborhood_search, + system::FluidSystem, neighbor::OpenBoundarySPHSystem, + u_system, u_neighbor) + # The current coordinates of fluids and open boundaries change over time. + + # TODO: Update only `active_coordinates` of open boundaries. + # Problem: Removing inactive particles from neighboring lists is necessary. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + +function update_nhs!(neighborhood_search, + system::OpenBoundarySPHSystem, neighbor::FluidSystem, + u_system, u_neighbor) + # The current coordinates of both open boundaries and fluids change over time. + + # TODO: Update only `active_coordinates` of open boundaries. + # Problem: Removing inactive particles from neighboring lists is necessary. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::FluidSystem, u_system, u_neighbor) @@ -725,6 +769,14 @@ function update_nhs!(neighborhood_search, return neighborhood_search end +function update_nhs!(neighborhood_search, + system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, + neighbor::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + function check_configuration(systems) foreach_system(systems) do system check_configuration(system, systems) diff --git a/src/general/system.jl b/src/general/system.jl index a4c7754c4..c3f10caad 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -17,7 +17,16 @@ initialize!(system, neighborhood_search) = system @inline n_moving_particles(system) = nparticles(system) @inline eachparticle(system) = Base.OneTo(nparticles(system)) -@inline each_moving_particle(system) = Base.OneTo(n_moving_particles(system)) + +# Wrapper for systems with `SystemBuffer` +@inline each_moving_particle(system) = each_moving_particle(system, system.buffer) +@inline each_moving_particle(system, ::Nothing) = Base.OneTo(n_moving_particles(system)) + +@inline active_coordinates(u, system) = active_coordinates(u, system, system.buffer) +@inline active_coordinates(u, system, ::Nothing) = current_coordinates(u, system) + +@inline active_particles(system) = active_particles(system, system.buffer) +@inline active_particles(system, ::Nothing) = eachparticle(system) # This should not be dispatched by system type. We always expect to get a column of `A`. @inline function extract_svector(A, system, i) @@ -64,6 +73,10 @@ end return zero(SVector{ndims(system), eltype(system)}) end +@inline set_particle_density!(v, system, particle, density) = v + +@inline set_particle_pressure!(v, system, particle, pressure) = v + @inline function smoothing_kernel(system, distance) (; smoothing_kernel, smoothing_length) = system return kernel(smoothing_kernel, distance, smoothing_length) @@ -103,6 +116,9 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t) return system end -function update_final!(system, v, u, v_ode, u_ode, semi, t) +function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system end + +# Only for systems requiring a mandatory callback +reset_callback_flag!(system) = system diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl new file mode 100644 index 000000000..586d565b1 --- /dev/null +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -0,0 +1,340 @@ +@doc raw""" + InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + +Inflow boundary zone for [`OpenBoundarySPHSystem`](@ref). + +The specified plane (line in 2D or rectangle in 3D) will be extruded in upstream +direction (the direction opposite to `flow_direction`) to create a box for the boundary zone. +There are three ways to specify the actual shape of the inflow: +1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then + be filled with inflow particles (default). +2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, + which is then extruded in upstream direction to create the inflow particles. + - In 2D, the shape must be either an initial condition with 2D coordinates, which lies + on the line specified by `plane`, or an initial condition with 1D coordinates, which lies + on the line specified by `plane` when a y-coordinate of `0` is added. + - In 3D, the shape must be either an initial condition with 3D coordinates, which lies + in the rectangle specified by `plane`, or an initial condition with 2D coordinates, + which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. +3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, + which will be used for the inflow particles. + +!!! note "Note" + Particles outside the boundary zone box will be removed. + +# Keywords +- `plane`: Tuple of points defining a part of the surface of the domain. + The points must either span a line in 2D or a rectangle in 3D. + This line or rectangle is then extruded in upstream direction to obtain + the boundary zone. + In 2D, pass two points ``(A, B)``, so that the interval ``[A, B]`` is + the inflow surface. + In 3D, pass three points ``(A, B, C)``, so that the rectangular inflow surface + is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. + These two vectors must be orthogonal. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in upstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for the inflow particles. + Particles outside the boundary zone will be removed. + Do not use together with `extrude_geometry`. +- `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane + and is extruded upstream to obtain the inflow particles. + See point 2 above for more details. + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction=[1.0, 0.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction=[0.0, 0.0, 1.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct InFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the inflow plane. + # The normal vector must point in upstream direction for an inflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to inflow plane")) + else + # Flip the normal vector to point in the opposite direction of `flow_direction` + spanning_set[:, 1] .*= -sign(dot_) + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove particles outside the boundary zone. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end +end + +@doc raw""" + OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + +Outflow boundary zone for [`OpenBoundarySPHSystem`](@ref). + +The specified plane (line in 2D or rectangle in 3D) will be extruded in downstream +direction (the direction in `flow_direction`) to create a box for the boundary zone. +There are three ways to specify the actual shape of the outflow: +1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then + be filled with outflow particles (default). +2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, + which is then extruded in downstream direction to create the outflow particles. + - In 2D, the shape must be either an initial condition with 2D coordinates, which lies + on the line specified by `plane`, or an initial condition with 1D coordinates, which lies + on the line specified by `plane` when a y-coordinate of `0` is added. + - In 3D, the shape must be either an initial condition with 3D coordinates, which lies + in the rectangle specified by `plane`, or an initial condition with 2D coordinates, + which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. +3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, + which will be used for the outflow particles. + +!!! note "Note" + Particles outside the boundary zone box will be removed. + +# Keywords +- `plane`: Tuple of points defining a part of the surface of the domain. + The points must either span a line in 2D or a rectangle in 3D. + This line or rectangle is then extruded in downstream direction to obtain + the boundary zone. + In 2D, pass two points ``(A, B)``, so that the interval ``[A, B]`` is + the outflow surface. + In 3D, pass three points ``(A, B, C)``, so that the rectangular outflow surface + is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. + These two vectors must be orthogonal. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in downstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for the outflow particles. + Particles outside the boundary zone will be removed. + Do not use together with `extrude_geometry`. +- `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane + and is extruded downstream to obtain the outflow particles. + See point 2 above for more details. + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction = [1.0, 0.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction = [0.0, 0.0, 1.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the outflow plane. + # The normal vector must point in downstream direction for an outflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to outflow plane")) + else + # Flip the normal vector to point in `flow_direction` + spanning_set[:, 1] .*= sign(dot_) + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove particles outside the boundary zone. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end +end + +@inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS + +function calculate_spanning_vectors(plane, zone_width) + return spanning_vectors(Tuple(plane), zone_width), SVector(plane[1]...) +end + +function spanning_vectors(plane_points::NTuple{2}, zone_width) + plane_size = plane_points[2] - plane_points[1] + + # Calculate normal vector of plane + b = normalize([-plane_size[2], plane_size[1]]) * zone_width + + return hcat(b, plane_size) +end + +function spanning_vectors(plane_points::NTuple{3}, zone_width) + # Vectors spanning the plane + edge1 = plane_points[2] - plane_points[1] + edge2 = plane_points[3] - plane_points[1] + + if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) + throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal")) + end + + # Calculate normal vector of plane + c = Vector(normalize(cross(edge2, edge1)) * zone_width) + + return hcat(c, edge1, edge2) +end + +@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords) + (; zone_origin, spanning_set) = boundary_zone + particle_position = particle_coords - zone_origin + + return is_in_boundary_zone(spanning_set, particle_position) +end + +@inline function is_in_boundary_zone(spanning_set::AbstractArray, + particle_position::SVector{NDIMS}) where {NDIMS} + for dim in 1:NDIMS + 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 boundary zone + return false + end + end + + # Particle is in boundary zone + return true +end + +function remove_outside_particles(initial_condition, spanning_set, zone_origin) + (; coordinates, density, particle_spacing) = initial_condition + + in_zone = trues(nparticles(initial_condition)) + + for particle in eachparticle(initial_condition) + current_position = current_coords(coordinates, initial_condition, particle) + particle_position = current_position - zone_origin + + in_zone[particle] = is_in_boundary_zone(spanning_set, particle_position) + end + + return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), + particle_spacing) +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl new file mode 100644 index 000000000..43c0a6ff1 --- /dev/null +++ b/src/schemes/boundary/open_boundary/system.jl @@ -0,0 +1,483 @@ +@doc raw""" + OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed, + fluid_system::FluidSystem, buffer_size::Integer, + reference_velocity=zeros(ndims(boundary_zone)), + reference_pressure=0.0, + reference_density=first(boundary_zone.initial_condition.density)) + +Open boundary system for in- and outflow particles. +These open boundaries use the characteristic variables to propagate the appropriate values +to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information +about the method see [description below](@ref method_of_characteristics). + +# Arguments +- `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary. + +# Keywords +- `sound_speed`: Speed of sound. +- `fluid_system`: The corresponding fluid system +- `buffer_size`: Number of buffer particles. +- `reference_velocity`: Reference velocity is either a function mapping each particle's coordinates + and time to its velocity, an array where the ``i``-th column holds + the velocity of particle ``i`` or, for a constant fluid velocity, + a vector holding this velocity. Velocity is constant zero by default. +- `reference_pressure`: Reference pressure is either a function mapping each particle's coordinates + and time to its pressure, a vector holding the pressure of each particle, + or a scalar for a constant pressure over all particles. + Pressure is constant zero by default. +- `reference_density`: Reference density is either a function mapping each particle's coordinates + and time to its density, a vector holding the density of each particle, + or a scalar for a constant density over all particles. + Density is the density of the first particle in the initial condition by default. + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D, RV, RP, + RD, B} <: System{NDIMS, IC} + initial_condition :: IC + fluid_system :: FS + mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] + density :: ARRAY1D # Array{ELTYPE, 1}: [particle] + volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] + pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] + characteristics :: ARRAY2D # Array{ELTYPE, 2}: [characteristic, particle] + previous_characteristics :: ARRAY2D # Array{ELTYPE, 2}: [characteristic, particle] + sound_speed :: ELTYPE + boundary_zone :: BZ + flow_direction :: SVector{NDIMS, ELTYPE} + reference_velocity :: RV + reference_pressure :: RP + reference_density :: RD + buffer :: B + update_callback_used :: Ref{Bool} + + function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed, + fluid_system::FluidSystem, buffer_size::Integer, + reference_velocity=zeros(ndims(boundary_zone)), + reference_pressure=0.0, + reference_density=first(boundary_zone.initial_condition.density)) + (; initial_condition) = boundary_zone + + buffer = SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + if !(reference_velocity isa Function || + (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) + throw(ArgumentError("`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "an array where the ``i``-th column holds the velocity of particle ``i`` " * + "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) + else + reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) + end + + if !(reference_pressure isa Function || reference_pressure isa Real) + throw(ArgumentError("`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "a vector holding the pressure of each particle, or a scalar")) + else + reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) + end + + if !(reference_density isa Function || reference_density isa Real) + throw(ArgumentError("`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "a vector holding the density of each particle, or a scalar")) + else + reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) + end + + mass = copy(initial_condition.mass) + pressure = [reference_pressure_(initial_condition.coordinates[:, i], 0.0) + for i in eachparticle(initial_condition)] + density = copy(initial_condition.density) + volume = similar(initial_condition.density) + + characteristics = zeros(ELTYPE, 3, length(mass)) + previous_characteristics = zeros(ELTYPE, 3, length(mass)) + + flow_direction_ = boundary_zone.flow_direction + + return new{typeof(boundary_zone), NDIMS, ELTYPE, typeof(initial_condition), + typeof(fluid_system), typeof(mass), typeof(characteristics), + typeof(reference_velocity_), typeof(reference_pressure_), + typeof(reference_density_), + typeof(buffer)}(initial_condition, fluid_system, mass, density, volume, + pressure, characteristics, previous_characteristics, + sound_speed, boundary_zone, flow_direction_, + reference_velocity_, reference_pressure_, + reference_density_, buffer, false) + end +end + +timer_name(::OpenBoundarySPHSystem) = "open_boundary" +vtkname(system::OpenBoundarySPHSystem) = "open_boundary" + +function Base.show(io::IO, system::OpenBoundarySPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "OpenBoundarySPHSystem{", ndims(system), "}(") + print(io, type2string(system.boundary_zone)) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "OpenBoundarySPHSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + summary_line(io, "fluid system", type2string(system.fluid_system)) + summary_line(io, "boundary", type2string(system.boundary_zone)) + summary_line(io, "flow direction", system.flow_direction) + summary_line(io, "prescribed velocity", string(nameof(system.reference_velocity))) + summary_line(io, "prescribed pressure", string(nameof(system.reference_pressure))) + summary_line(io, "prescribed density", string(nameof(system.reference_density))) + summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3)) + summary_footer(io) + end +end + +function reset_callback_flag(system::OpenBoundarySPHSystem) + system.update_callback_used[] = false + + return system +end + +update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true + +@inline source_terms(system::OpenBoundarySPHSystem) = nothing + +@inline hydrodynamic_mass(system::OpenBoundarySPHSystem, particle) = system.mass[particle] + +@inline function particle_density(v, system::OpenBoundarySPHSystem, particle) + return system.density[particle] +end + +@inline function particle_pressure(v, system::OpenBoundarySPHSystem, particle) + return system.pressure[particle] +end + +@inline function update_quantities!(system::OpenBoundarySPHSystem, v, u, t) + (; density, pressure, characteristics, flow_direction, sound_speed, + reference_velocity, reference_pressure, reference_density) = system + + # Update quantities based on the characteristic variables + @threaded for particle in each_moving_particle(system) + particle_position = current_coords(u, system, particle) + + J1 = characteristics[1, particle] + J2 = characteristics[2, particle] + J3 = characteristics[3, particle] + + rho_ref = reference_density(particle_position, t) + density[particle] = rho_ref + ((-J1 + 0.5 * (J2 + J3)) / sound_speed^2) + + p_ref = reference_pressure(particle_position, t) + pressure[particle] = p_ref + 0.5 * (J2 + J3) + + v_ref = reference_velocity(particle_position, t) + rho = density[particle] + v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction + + for dim in 1:ndims(system) + v[dim, particle] = v_[dim] + end + end + + return system +end + +function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + if !update_from_callback && !(system.update_callback_used[]) + throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`")) + end + + @trixi_timeit timer() "evaluate characteristics" evaluate_characteristics!(system, v, u, + v_ode, u_ode, + semi, t) +end + +# ==== Characteristics +# J1: Associated with convection and entropy and propagates at flow velocity. +# J2: Propagates downstream to the local flow +# J3: Propagates upstream to the local flow +function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) + (; volume, characteristics, previous_characteristics, boundary_zone) = system + + for particle in eachparticle(system) + previous_characteristics[1, particle] = characteristics[1, particle] + previous_characteristics[2, particle] = characteristics[2, particle] + previous_characteristics[3, particle] = characteristics[3, particle] + end + + set_zero!(characteristics) + set_zero!(volume) + + # Evaluate the characteristic variables with the fluid system + evaluate_characteristics!(system, system.fluid_system, v, u, v_ode, u_ode, semi, t) + + # Only some of the in-/outlet particles are in the influence of the fluid particles. + # Thus, we compute the characteristics for the particles that are outside the influence + # of fluid particles by using the average of the values of the previous time step. + # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119 + @threaded for particle in each_moving_particle(system) + + # Particle is outside of the influence of fluid particles + if isapprox(volume[particle], 0.0) + + # Using the average of the values at the previous time step for particles which + # are outside of the influence of fluid particles. + avg_J1 = 0.0 + avg_J2 = 0.0 + avg_J3 = 0.0 + counter = 0 + + for neighbor in each_moving_particle(system) + # Make sure that only neighbors in the influence of + # the fluid particles are used. + if volume[neighbor] > sqrt(eps()) + avg_J1 += previous_characteristics[1, neighbor] + avg_J2 += previous_characteristics[2, neighbor] + avg_J3 += previous_characteristics[3, neighbor] + counter += 1 + end + end + + characteristics[1, particle] = avg_J1 / counter + characteristics[2, particle] = avg_J2 / counter + characteristics[3, particle] = avg_J3 / counter + else + characteristics[1, particle] /= volume[particle] + characteristics[2, particle] /= volume[particle] + characteristics[3, particle] /= volume[particle] + end + prescribe_conditions!(characteristics, particle, boundary_zone) + end + + return system +end + +function evaluate_characteristics!(system, neighbor_system::FluidSystem, + v, u, v_ode, u_ode, semi, t) + (; volume, sound_speed, characteristics, flow_direction, + reference_velocity, reference_pressure, reference_density) = system + + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + + nhs = get_neighborhood_search(system, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all fluid neighbors within the kernel cutoff + for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, + nhs) do particle, neighbor, pos_diff, distance + neighbor_position = current_coords(u_neighbor_system, neighbor_system, neighbor) + + # Determine current and prescribed quantities + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_ref = reference_density(neighbor_position, t) + + p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) + p_ref = reference_pressure(neighbor_position, t) + + v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + v_neighbor_ref = reference_velocity(neighbor_position, t) + + # Determine characteristic variables + density_term = -sound_speed^2 * (rho_b - rho_ref) + pressure_term = p_b - p_ref + velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) + + kernel_ = smoothing_kernel(neighbor_system, distance) + + characteristics[1, particle] += (density_term + pressure_term) * kernel_ + characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ + characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ + + volume[particle] += kernel_ + end + + return system +end + +@inline function prescribe_conditions!(characteristics, particle, ::OutFlow) + # J3 is prescribed (i.e. determined from the exterior of the domain). + # J1 and J2 is transimtted from the domain interior. + characteristics[3, particle] = zero(eltype(characteristics)) + + return characteristics +end + +@inline function prescribe_conditions!(characteristics, particle, ::InFlow) + # Allow only J3 to propagate upstream to the boundary + characteristics[1, particle] = zero(eltype(characteristics)) + characteristics[2, particle] = zero(eltype(characteristics)) + + return characteristics +end + +# This function is called by the `UpdateCallback`, as the integrator array might be modified +function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ode, + semi, t) + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, semi) + + # Update density, pressure and velocity based on the characteristic variables. + # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 + @trixi_timeit timer() "update quantities" update_quantities!(system, v, u, t) + + @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) + + # Update buffers + update_system_buffer!(system.buffer) + update_system_buffer!(system.fluid_system.buffer) +end + +update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) = system + +function check_domain!(system, v, u, v_ode, u_ode, semi) + (; boundary_zone, fluid_system) = system + + u_fluid = wrap_u(u_ode, fluid_system, semi) + v_fluid = wrap_v(v_ode, fluid_system, semi) + + neighborhood_search = get_neighborhood_search(system, fluid_system, semi) + + for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + # Check if boundary particle is outside the boundary zone + if !is_in_boundary_zone(boundary_zone, particle_coords) + convert_particle!(system, fluid_system, boundary_zone, particle, + v, u, v_fluid, u_fluid) + end + + # Check the neighboring fluid particles whether they're entering the boundary zone + for neighbor in PointNeighbors.eachneighbor(particle_coords, neighborhood_search) + fluid_coords = current_coords(u_fluid, fluid_system, neighbor) + + # Check if neighboring fluid particle is in boundary zone + if is_in_boundary_zone(boundary_zone, fluid_coords) + convert_particle!(fluid_system, system, boundary_zone, neighbor, + v, u, v_fluid, u_fluid) + end + end + end + + return system +end + +# Outflow particle is outside the boundary zone +@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, + boundary_zone::OutFlow, particle, v, u, + v_fluid, u_fluid) + deactivate_particle!(system, particle, u) + + return system +end + +# Inflow particle is outside the boundary zone +@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, + boundary_zone::InFlow, particle, v, u, + v_fluid, u_fluid) + (; spanning_set) = boundary_zone + + # Activate a new particle in simulation domain + transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u) + + # Reset position of boundary particle + for dim in 1:ndims(system) + u[dim, particle] += spanning_set[1][dim] + end + + return system +end + +# Fluid particle is in boundary zone +@inline function convert_particle!(fluid_system::FluidSystem, system, + boundary_zone, particle, v, u, v_fluid, u_fluid) + # Activate particle in boundary zone + transfer_particle!(system, fluid_system, particle, v, u, v_fluid, u_fluid) + + # Deactivate particle in interior domain + deactivate_particle!(fluid_system, particle, u_fluid) + + return fluid_system +end + +@inline function transfer_particle!(system_new, system_old, particle_old, + v_new, u_new, v_old, u_old) + particle_new = activate_next_particle(system_new) + + # Transfer densities + density = particle_density(v_old, system_old, particle_old) + set_particle_density!(v_new, system_new, particle_new, density) + + # Transfer pressure + pressure = particle_pressure(v_old, system_old, particle_old) + set_particle_pressure!(v_new, system_new, particle_new, pressure) + + # Exchange position and velocity + for dim in 1:ndims(system_new) + u_new[dim, particle_new] = u_old[dim, particle_old] + v_new[dim, particle_new] = v_old[dim, particle_old] + end + + # TODO: Only when using TVF: set tvf + + return system_new +end + +function write_v0!(v0, system::OpenBoundarySPHSystem) + (; initial_condition) = system + + for particle in eachparticle(system) + # Write particle velocities + for dim in 1:ndims(system) + v0[dim, particle] = initial_condition.velocity[dim, particle] + end + end + + return v0 +end + +function write_u0!(u0, system::OpenBoundarySPHSystem) + (; initial_condition) = system + + for particle in eachparticle(system) + # Write particle velocities + for dim in 1:ndims(system) + u0[dim, particle] = initial_condition.coordinates[dim, particle] + end + end + + return u0 +end + +function wrap_reference_function(function_::Function, ::Val) + # Already a function + return function_ +end + +# Name the function so that the summary box does know which kind of function this is +function wrap_reference_function(constant_scalar_::Number, ::Val) + return constant_scalar(coords, t) = constant_scalar_ +end + +# For vectors and tuples +# Name the function so that the summary box does know which kind of function this is +function wrap_reference_function(constant_vector_, ::Val{NDIMS}) where {NDIMS} + return constant_vector(coords, t) = SVector{NDIMS}(constant_vector_) +end diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 201a3bbb9..d9cfb1779 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -1,7 +1,7 @@ -# Interaction of boundary with other systems +# Interaction of boundary with other systems function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::BoundarySystem, + particle_system::Union{BoundarySystem, OpenBoundarySPHSystem}, neighbor_system) # TODO Solids and moving boundaries should be considered in the continuity equation return dv diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index cade5fbe9..73af6a843 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -13,7 +13,8 @@ The interaction between fluid and boundary particles is specified by the boundar - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface. Note: currently it is assumed that all fluids have the same adhesion coefficient. """ -struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS, IC} +struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, + CA} <: BoundarySystem{NDIMS, IC} initial_condition :: IC coordinates :: CO # Array{ELTYPE, 2} boundary_model :: BM @@ -21,6 +22,7 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) adhesion_coefficient :: ELTYPE cache :: CA + buffer :: Nothing # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. @@ -28,12 +30,10 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem ismoving, adhesion_coefficient, cache) ELTYPE = eltype(coordinates) - new{typeof(boundary_model), size(coordinates, 1), ELTYPE, - typeof(initial_condition), typeof(coordinates), - typeof(movement), typeof(ismoving), typeof(cache)}(initial_condition, - coordinates, boundary_model, - movement, ismoving, - adhesion_coefficient, cache) + new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), + typeof(coordinates), typeof(movement), typeof(ismoving), + typeof(cache)}(initial_condition, coordinates, boundary_model, movement, + ismoving, adhesion_coefficient, cache, nothing) end end @@ -73,6 +73,7 @@ struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC, coordinates :: ARRAY2D # [dimension, particle] radius :: ARRAY1D # [particle] normal_stiffness :: ELTYPE + buffer :: Nothing function BoundaryDEMSystem(initial_condition, normal_stiffness) coordinates = initial_condition.coordinates @@ -80,9 +81,9 @@ struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC, ones(length(initial_condition.mass)) NDIMS = size(coordinates, 1) - return new{NDIMS, eltype(coordinates), typeof(initial_condition), - typeof(radius), typeof(coordinates)}(initial_condition, coordinates, - radius, normal_stiffness) + return new{NDIMS, eltype(coordinates), typeof(initial_condition), typeof(radius), + typeof(coordinates)}(initial_condition, coordinates, radius, + normal_stiffness, nothing) end end @@ -335,7 +336,8 @@ end # This update depends on the computed quantities of the fluid system and therefore # has to be in `update_final!` after `update_quantities!`. -function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t) +function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) (; boundary_model) = system update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f59bca7d3..4093e8ccb 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -4,7 +4,7 @@ pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, NDIMS), + acceleration=ntuple(_ -> 0.0, NDIMS), buffer_size=nothing, source_terms=nothing) System for particles of a fluid. @@ -28,6 +28,8 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) +- `buffer_size`: Number of buffer particles. + This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure)` (which are the quantities of a single particle), returning a `Tuple` @@ -40,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, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC @@ -53,6 +55,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, correction :: Nothing pressure_acceleration_formulation :: PF source_terms :: ST + buffer :: B cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, @@ -62,7 +65,12 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - source_terms=nothing) + source_terms=nothing, buffer_size=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -89,9 +97,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(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, source_terms, cache) + smoothing_length, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, source_terms, + buffer, cache) end end @@ -113,7 +123,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst show(io, system) else summary_header(io, "EntropicallyDampedSPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) + if system.buffer isa SystemBuffer + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + else + summary_line(io, "#particles", nparticles(system)) + end summary_line(io, "density calculator", system.density_calculator |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) @@ -145,6 +160,18 @@ end return v[end, particle] end +# WARNING! +# These functions are intended to be used internally to set the pressure +# of newly activated particles in a callback. +# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` +# outside of callbacks. +@inline function set_particle_pressure!(v, system::EntropicallyDampedSPHSystem, particle, + pressure) + v[end, particle] = pressure + + return v +end + @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 056662bbc..68391b285 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -1,3 +1,7 @@ +@inline function set_particle_density!(v, system::FluidSystem, particle, density) + set_particle_density!(v, system, system.density_calculator, particle, density) +end + function create_cache_density(initial_condition, ::SummationDensity) density = similar(initial_condition.density) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 554f41c14..a678b40d4 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -11,6 +11,19 @@ function dv_viscosity(particle_system, neighbor_system, sound_speed, m_a, m_b, rho_mean) end +function dv_viscosity(particle_system, neighbor_system::OpenBoundarySPHSystem, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) + # No viscosity in the open boundary system. Use viscosity of the fluid system. + viscosity = viscosity_model(particle_system) + + return dv_viscosity(viscosity, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) +end + function dv_viscosity(viscosity, particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index cb9a1da87..a8e28670b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -4,6 +4,7 @@ smoothing_kernel, smoothing_length; viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), + buffer_size=nothing, correction=nothing, source_terms=nothing) System for particles of a fluid. @@ -26,6 +27,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. See [`ArtificialViscosityMonaghan`](@ref) or [`ViscosityAdami`](@ref). - `density_diffusion`: Density diffusion terms for this system. See [`DensityDiffusion`](@ref). - `acceleration`: Acceleration vector for the system. (default: zero vector) +- `buffer_size`: Number of buffer particles. + This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure)` @@ -42,7 +45,7 @@ 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, SRFT, C} <: FluidSystem{NDIMS, IC} + V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -57,6 +60,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, pressure_acceleration_formulation :: PF source_terms :: ST surface_tension :: SRFT + buffer :: B cache :: C end @@ -66,11 +70,17 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, smoothing_length; pressure_acceleration=nothing, + buffer_size=nothing, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, surface_tension=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -108,10 +118,11 @@ function WeaklyCompressibleSPHSystem(initial_condition, return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, acceleration_, - viscosity, density_diffusion, correction, - pressure_acceleration, source_terms, surface_tension, - cache) + smoothing_kernel, smoothing_length, + acceleration_, viscosity, + density_diffusion, correction, + pressure_acceleration, + source_terms, surface_tension, buffer, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -166,7 +177,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst show(io, system) else summary_header(io, "WeaklyCompressibleSPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) + if system.buffer isa SystemBuffer + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + else + summary_line(io, "#particles", nparticles(system)) + end summary_line(io, "density calculator", system.density_calculator |> typeof |> nameof) summary_line(io, "correction method", diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index b81a7768b..fadcb240e 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -1,5 +1,8 @@ # Include all schemes without rhs first. The rhs depends on the systems to define # interactions between the different system types. +# Viscosity requires the open boundary system. +include("boundary/open_boundary/boundary_zones.jl") +include("boundary/open_boundary/system.jl") include("fluid/fluid.jl") include("boundary/boundary.jl") include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index 5aadfcfcc..eb30dfb34 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -36,6 +36,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, I damping_coefficient :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} source_terms :: ST + buffer :: Nothing function DEMSystem(initial_condition, normal_stiffness, elastic_modulus, poissons_ratio; damping_coefficient=0.0001, @@ -56,7 +57,8 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, I return new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(source_terms)}(initial_condition, mass, radius, elastic_modulus, poissons_ratio, normal_stiffness, - damping_coefficient, acceleration_, source_terms) + damping_coefficient, acceleration_, source_terms, + nothing) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 5b3021b1c..9203425fb 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -69,6 +69,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, boundary_model :: BM penalty_force :: PF source_terms :: ST + buffer :: Nothing function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, @@ -116,7 +117,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, n_moving_particles, young_modulus, poisson_ratio, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, - penalty_force, source_terms) + penalty_force, source_terms, nothing) end end @@ -246,7 +247,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode return system end -function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t) +function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) (; boundary_model) = system # Only update boundary model diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index 9b7a62401..82df9aaa3 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -8,8 +8,10 @@ Extrude either a line, a plane or a shape along a specific direction. # Arguments - `geometry`: Either particle coordinates or an [`InitialCondition`](@ref) defining a 2D shape to extrude to a 3D volume, or two 2D points - defining a line to extrude to a plane in 2D, or three 3D points defining - a parallelogram to extrude to a parallelepiped. + ``(A, B)`` defining the interval ``[A, B]`` to extrude to a plane + in 2D, or three 3D points ``(A, B, C)`` defining the parallelogram + spanned by the vectors ``\widehat{AB}`` and ``\widehat {AC}`` to extrude + to a parallelepiped. # Keywords - `particle_spacing`: Spacing between the particles. Can be omitted when `geometry` is an diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index a3fb7de7e..f7e4d5b23 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -12,10 +12,8 @@ end RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; size=(600, 400)) # Default size systems_data = map(semi.systems) do system - (; initial_condition) = system - u = wrap_u(u_ode, system, semi) - coordinates = current_coordinates(u, system) + coordinates = active_coordinates(u, system) x = collect(coordinates[1, :]) y = collect(coordinates[2, :]) @@ -24,8 +22,8 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; particle_spacing = 0.0 end - x_min, y_min = minimum(initial_condition.coordinates, dims=2) .- 0.5particle_spacing - x_max, y_max = maximum(initial_condition.coordinates, dims=2) .+ 0.5particle_spacing + x_min, y_min = minimum(coordinates, dims=2) .- 0.5particle_spacing + x_max, y_max = maximum(coordinates, dims=2) .+ 0.5particle_spacing return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, label=timer_name(system)) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 2a4c9f207..0372e7167 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -56,7 +56,8 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - @trixi_timeit timer() "update systems" update_systems_and_nhs(v_ode, u_ode, semi, t) + @trixi_timeit timer() "update systems" update_systems_and_nhs(v_ode, u_ode, semi, t; + update_from_callback=true) filenames = system_names(systems) @@ -94,7 +95,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix # Reset the collection when the iteration is 0 pvd = paraview_collection(collection_file; append=iter > 0) - points = PointNeighbors.periodic_coords(current_coordinates(u, system), + points = PointNeighbors.periodic_coords(active_coordinates(u, system), periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] @@ -112,7 +113,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) # Store particle index - vtk["index"] = eachparticle(system) + vtk["index"] = active_particles(system) vtk["time"] = t if write_meta_data @@ -216,11 +217,12 @@ function write2vtk!(vtk, v, u, t, system; write_meta_data=true) end function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) - vtk["velocity"] = view(v, 1:ndims(system), :) + vtk["velocity"] = [current_velocity(v, system, particle) + for particle in active_particles(system)] vtk["density"] = [particle_density(v, system, particle) - for particle in eachparticle(system)] + for particle in active_particles(system)] vtk["pressure"] = [particle_pressure(v, system, particle) - for particle in eachparticle(system)] + for particle in active_particles(system)] if write_meta_data vtk["acceleration"] = system.acceleration @@ -295,6 +297,37 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end +function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true) + (; reference_velocity, reference_pressure, reference_density) = system + + vtk["velocity"] = [current_velocity(v, system, particle) + for particle in active_particles(system)] + vtk["density"] = [particle_density(v, system, particle) + for particle in active_particles(system)] + vtk["pressure"] = [particle_pressure(v, system, particle) + for particle in active_particles(system)] + + NDIMS = ndims(system) + ELTYPE = eltype(system) + coords = reinterpret(reshape, SVector{NDIMS, ELTYPE}, active_coordinates(u, system)) + + vtk["prescribed_velocity"] = stack(reference_velocity.(coords, t)) + vtk["prescribed_density"] = reference_density.(coords, t) + vtk["prescribed_pressure"] = reference_pressure.(coords, t) + + if write_meta_data + vtk["boundary_zone"] = type2string(system.boundary_zone) + vtk["sound_speed"] = system.sound_speed + vtk["width"] = round(system.boundary_zone.zone_width, digits=3) + vtk["flow_direction"] = system.flow_direction + vtk["velocity_function"] = string(nameof(system.reference_velocity)) + vtk["pressure_function"] = string(nameof(system.reference_pressure)) + vtk["density_function"] = string(nameof(system.reference_density)) + end + + return vtk +end + function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true) write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 0f0032607..dbe03a130 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -114,6 +114,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/pipe_flow_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "pipe_flow_2d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/dam_break_2d_surface_tension.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/general/buffer.jl b/test/general/buffer.jl new file mode 100644 index 000000000..7a8371ea7 --- /dev/null +++ b/test/general/buffer.jl @@ -0,0 +1,68 @@ +@testset verbose=true "`SystemBuffer`" begin + # Mock fluid system + struct FluidSystemMock3 <: TrixiParticles.FluidSystem{2, Nothing} end + + inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2, + open_boundary_layers=2, density=1.0, flow_direction=[1.0, 0.0]) + system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, fluid_system=FluidSystemMock3(), + buffer_size=0) + system_buffer = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=5, + fluid_system=FluidSystemMock3()) + + n_particles = nparticles(system) + + @testset "Iterators" begin + @test TrixiParticles.each_moving_particle(system) == 1:n_particles + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:n_particles + + particle_id = TrixiParticles.activate_next_particle(system_buffer) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:(n_particles + 1) + + TrixiParticles.deactivate_particle!(system_buffer, particle_id, + ones(2, particle_id)) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:n_particles + + particle_id = 5 + TrixiParticles.deactivate_particle!(system_buffer, particle_id, + ones(2, particle_id)) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == + setdiff(1:n_particles, particle_id) + end + + @testset "Allocate Buffer" begin + initial_condition = rectangular_patch(0.1, (3, 3), perturbation_factor=0.0) + buffer = TrixiParticles.SystemBuffer(nparticles(initial_condition), 7) + + ic_with_buffer = TrixiParticles.allocate_buffer(initial_condition, buffer) + + @test nparticles(ic_with_buffer) == nparticles(initial_condition) + 7 + + masses = initial_condition.mass[1] .* ones(nparticles(ic_with_buffer)) + @test masses == ic_with_buffer.mass + + densities = initial_condition.density[1] .* ones(nparticles(ic_with_buffer)) + @test densities == ic_with_buffer.density + + pressures = initial_condition.pressure[1] .* ones(nparticles(ic_with_buffer)) + @test pressures == ic_with_buffer.pressure + + @testset "Illegal Input" begin + # The rectangular patch has a perturbed, non-constant density + ic = rectangular_patch(0.1, (3, 3)) + buffer = TrixiParticles.SystemBuffer(9, 7) + + error_str = "`initial_condition.density` needs to be constant when using `SystemBuffer`" + @test_throws ArgumentError(error_str) TrixiParticles.allocate_buffer(ic, buffer) + end + end +end diff --git a/test/general/general.jl b/test/general/general.jl index b59cff789..27f60ebba 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -3,3 +3,4 @@ include("smoothing_kernels.jl") include("density_calculator.jl") include("semidiscretization.jl") include("interpolation.jl") +include("buffer.jl") diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 633a0e240..4aff1f4fd 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -131,6 +131,9 @@ v2 = zeros(4 * 3) v_ode = vcat(vec(v1), v2) + # Avoid `SystemBuffer` barrier + TrixiParticles.each_moving_particle(system::Union{System1, System2}) = TrixiParticles.eachparticle(system) + TrixiParticles.add_source_terms!(dv_ode, v_ode, u_ode, semi) dv1 = TrixiParticles.wrap_v(dv_ode, system1, semi) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl new file mode 100644 index 000000000..7a34a4bca --- /dev/null +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -0,0 +1,221 @@ +@testset verbose=true "Boundary Zone" begin + @testset verbose=true "Boundary Zone 2D" begin + particle_spacing = 0.2 + open_boundary_layers = 4 + + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset verbose=true "Points $i" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + + plane_size = point_2 - point_1 + + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test isapprox(sign_ * flow_directions[j], + normalize(boundary_zone.spanning_set[1]), atol=1e-14) + @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]), + atol=1e-14) + end + end + end + end + + @testset verbose=true "Boundary Zone 3D" begin + particle_spacing = 0.05 + open_boundary_layers = 4 + + plane_points_1 = [ + [0.0, 0.0, 0.0], + [0.3113730847835541, 0.19079485535621643, -0.440864622592926], + ] + plane_points_2 = [ + [1.0, 0.0, 0.0], + [-0.10468611121177673, 0.252103328704834, -0.44965094327926636], + ] + plane_points_3 = [ + [0.0, 1.0, 0.0], + [0.3113730847835541, 0.25057315826416016, -0.02374829351902008], + ] + + @testset verbose=true "Points $i" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + point_3 = plane_points_3[i] + + edge1 = point_2 - point_1 + edge2 = point_3 - point_1 + + flow_directions = [ + normalize(cross(edge1, edge2)), + -normalize(cross(edge1, edge2)), + ] + + @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test plane_points_3[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[3] + @test isapprox(sign_ * flow_directions[j], + normalize(boundary_zone.spanning_set[1]), atol=1e-14) + @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]), + atol=1e-14) + end + end + end + end + + @testset verbose=true "Particle In Boundary Zone 2D" begin + plane_points = [[-0.2, -0.5], [0.3, 0.6]] + plane_size = plane_points[2] - plane_points[1] + + flow_direction = normalize([-plane_size[2], plane_size[1]]) + + inflow = InFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? sqrt(eps()) : -sqrt(eps()) + + point1 = plane_points[1] + point2 = plane_points[2] + point3 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin + + query_points = Dict( + "Behind" => ([-1.0, -1.0], false), + "Before" => ([2.0, 2.0], false), + "Closely On Point 1" => (point1 + perturb_ * flow_direction, false), + "Closely On Point 2" => (point2 + perturb_ * flow_direction, false), + "Closely On Point 3" => (point3 - perturb_ * flow_direction, false)) + + @testset verbose=true "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == + TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position) + end + end + end + + @testset verbose=true "Particle In Boundary Zone 3D" begin + point1 = [-0.2, -0.5, 0.0] + point2 = [0.3, 0.5, 0.0] + point3 = [0.13111173850909402, -0.665555869254547, 0.0] + + flow_direction = normalize(cross(point2 - point1, point3 - point1)) + + inflow = InFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? eps() : -eps() + point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin + + query_points = Dict( + "Behind" => ([-1.0, -1.0, 1.2], false), + "Before" => ([2.0, 2.0, -1.2], false), + "Closely On Point 1" => (point1 + perturb_ * flow_direction, false), + "Closely On Point 2" => (point2 + perturb_ * flow_direction, false), + "Closely On Point 3" => (point3 + perturb_ * flow_direction, false), + "Closely On Point 4" => (point4 - perturb_ * flow_direction, false)) + + @testset verbose=true "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == + TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position) + end + end + end + + @testset verbose=true "Illegal Inputs" begin + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + flow_direction = [0.0, 0.0, 1.0] + + error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal" + + @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + @test_throws ArgumentError(error_str) OutFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + rectangular_plane = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + flow_direction = [0.0, 1.0, 0.0] + + error_str = "`flow_direction` is not normal to inflow plane" + + @test_throws ArgumentError(error_str) InFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + error_str = "`flow_direction` is not normal to outflow plane" + + @test_throws ArgumentError(error_str) OutFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + end +end diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl new file mode 100644 index 000000000..1fe20722a --- /dev/null +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -0,0 +1,134 @@ +@testset verbose=true "Characteristic Variables" begin + particle_spacing = 0.1 + + # Number of boundary particles in the influence of fluid particles + influenced_particles = [20, 52, 26] + + open_boundary_layers = 8 + sound_speed = 20.0 + density = 1000.0 + pressure = 5.0 + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.2particle_spacing + + # Prescribed quantities + reference_velocity = (pos, t) -> SVector(t, 0.0) + reference_pressure = (pos, t) -> 50_000.0 * t + reference_density = (pos, t) -> 1000.0 * t + + # Plane points of open boundary + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset "Points $i" for i in eachindex(plane_points_1) + n_influenced = influenced_particles[i] + + plane_points = [plane_points_1[i], plane_points_2[i]] + + plane_size = plane_points[2] - plane_points[1] + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset "Flow Direction $j" for j in eachindex(flow_directions) + flow_direction = flow_directions[j] + inflow = InFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + outflow = OutFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "`$(nameof(typeof(boundary_zone)))`" for boundary_zone in boundary_zones + sign_ = (boundary_zone isa InFlow) ? 1 : -1 + fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, + density, pressure, + direction=(sign_ * flow_direction)) + + fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, + buffer_size=0, + density_calculator=ContinuityDensity(), + smoothing_length, sound_speed) + + boundary_system = OpenBoundarySPHSystem(boundary_zone; sound_speed, + fluid_system, buffer_size=0, + reference_velocity, + reference_pressure, + reference_density) + + semi = Semidiscretization(fluid_system, boundary_system) + + ode = semidiscretize(semi, (0.0, 5.0)) + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, boundary_system, semi) + u = TrixiParticles.wrap_u(u0_ode, boundary_system, semi) + + # ==== Characteristic Variables + # `J1 = -sound_speed^2 * (rho - rho_ref) + (p - p_ref)` + # `J2 = rho * sound_speed * (v - v_ref) + (p - p_ref)` + # `J3 = - rho * sound_speed * (v - v_ref) + (p - p_ref)` + function J1(t) + return -sound_speed^2 * (density - reference_density(0, t)) + + (pressure - reference_pressure(0, t)) + end + function J2(t) + return density * sound_speed * + dot(-reference_velocity(0, t), flow_direction) + + (pressure - reference_pressure(0, t)) + end + function J3(t) + return -density * sound_speed * + dot(-reference_velocity(0, t), flow_direction) + + (pressure - reference_pressure(0, t)) + end + + # First evaluation. + # Particles not influenced by the fluid have zero values. + t1 = 2.0 + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, t1) + evaluated_vars1 = boundary_system.characteristics + + if boundary_zone isa InFlow + @test all(isapprox.(evaluated_vars1[1, :], 0.0)) + @test all(isapprox.(evaluated_vars1[2, :], 0.0)) + @test all(isapprox.(evaluated_vars1[3, 1:n_influenced], J3(t1))) + @test all(isapprox.(evaluated_vars1[3, (n_influenced + 1):end], 0.0)) + + elseif boundary_zone isa OutFlow + @test all(isapprox.(evaluated_vars1[1, 1:n_influenced], J1(t1))) + @test all(isapprox.(evaluated_vars1[2, 1:n_influenced], J2(t1))) + @test all(isapprox.(evaluated_vars1[1:2, (n_influenced + 1):end], 0.0)) + @test all(isapprox.(evaluated_vars1[3, :], 0.0)) + end + + # Second evaluation. + # Particles not influenced by the fluid have previous values. + t2 = 3.0 + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, t2) + evaluated_vars2 = boundary_system.characteristics + + if boundary_zone isa InFlow + @test all(isapprox.(evaluated_vars2[1, :], 0.0)) + @test all(isapprox.(evaluated_vars2[2, :], 0.0)) + @test all(isapprox.(evaluated_vars2[3, 1:n_influenced], J3(t2))) + @test all(isapprox.(evaluated_vars2[3, (n_influenced + 1):end], J3(t1))) + + elseif boundary_zone isa OutFlow + @test all(isapprox.(evaluated_vars2[1, 1:n_influenced], J1(t2))) + @test all(isapprox.(evaluated_vars2[2, 1:n_influenced], J2(t2))) + @test all(isapprox.(evaluated_vars2[1, (n_influenced + 1):end], J1(t1))) + @test all(isapprox.(evaluated_vars2[2, (n_influenced + 1):end], J2(t1))) + @test all(isapprox.(evaluated_vars2[3, :], 0.0)) + end + end + end + end +end diff --git a/test/schemes/boundary/open_boundary/open_boundary.jl b/test/schemes/boundary/open_boundary/open_boundary.jl new file mode 100644 index 000000000..5fe09c703 --- /dev/null +++ b/test/schemes/boundary/open_boundary/open_boundary.jl @@ -0,0 +1,2 @@ +include("characteristic_variables.jl") +include("boundary_zone.jl") diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index 2e5e609e5..5cf7cd5bd 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -1,4 +1,5 @@ include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") include("boundary/dummy_particles/dummy_particles.jl") include("boundary/monaghan_kajtar/monaghan_kajtar.jl") +include("boundary/open_boundary/open_boundary.jl") include("fluid/fluid.jl") diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl new file mode 100644 index 000000000..1158aaff0 --- /dev/null +++ b/test/systems/open_boundary_system.jl @@ -0,0 +1,90 @@ +@testset verbose=true "`OpenBoundarySPHSystem`" begin + @testset verbose=true "Illegal Inputs" begin + plane = ([0.0, 0.0], [0.0, 1.0]) + flow_direction = (1.0, 0.0) + + # Mock fluid system + struct FluidSystemMock2 <: TrixiParticles.FluidSystem{2, Nothing} end + + inflow = InFlow(; plane, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=2) + + error_str = "`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "an array where the ``i``-th column holds the velocity of particle ``i`` " * + "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" + + reference_velocity = 1.0 + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_velocity) + + error_str = "`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "a vector holding the pressure of each particle, or a scalar" + + reference_pressure = [1.0, 1.0] + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_pressure) + + error_str = "`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "a vector holding the density of each particle, or a scalar" + + reference_density = [1.0, 1.0] + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_density) + end + @testset "`show`" begin + inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=0, + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}(InFlow) with 80 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 80 │ + │ #buffer_particles: ……………………………… 0 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary: ……………………………………………………… InFlow │ + │ flow direction: ……………………………………… [1.0, 0.0] │ + │ prescribed velocity: ………………………… constant_vector │ + │ prescribed pressure: ………………………… constant_scalar │ + │ prescribed density: …………………………… constant_scalar │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + + outflow = OutFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(outflow; sound_speed=1.0, buffer_size=0, + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}(OutFlow) with 80 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 80 │ + │ #buffer_particles: ……………………………… 0 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary: ……………………………………………………… OutFlow │ + │ flow direction: ……………………………………… [1.0, 0.0] │ + │ prescribed velocity: ………………………… constant_vector │ + │ prescribed pressure: ………………………… constant_scalar │ + │ prescribed density: …………………………… constant_scalar │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + end +end diff --git a/test/systems/systems.jl b/test/systems/systems.jl index 6aff417b3..697f575fa 100644 --- a/test/systems/systems.jl +++ b/test/systems/systems.jl @@ -2,4 +2,5 @@ include("wcsph_system.jl") include("edac_system.jl") include("solid_system.jl") include("boundary_system.jl") +include("open_boundary_system.jl") include("dem_system.jl")