Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Determine open-boundary-velocity from conservation equation #572

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Determine open-boundary-velocity from conservation equation #572

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
91ee930
determine velocity from conservation equation
LasNikas Jul 1, 2024
c054cbf
Merge branch 'main' into determine-velocity-with-rhs
LasNikas Jul 1, 2024
1a602f7
interpolate ref values
LasNikas Jul 2, 2024
e4a8c06
undo interpolation
LasNikas Jul 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 1 addition & 3 deletions examples/fluid/pipe_flow_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,6 @@ inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,

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]]),
Expand All @@ -98,8 +97,7 @@ outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2

open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system,
buffer_size=n_buffer_particles,
reference_pressure=pressure,
reference_velocity=velocity_function)
reference_pressure=pressure)

# ==========================================================================================
# ==== Boundary
Expand Down
12 changes: 6 additions & 6 deletions src/schemes/boundary/open_boundary/boundary_zones.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ struct InFlow{NDIMS, IC, S, ZO, ZW, FD}
zone_width :: ZW
flow_direction :: FD

function InFlow(; plane, flow_direction, density, particle_spacing,
function InFlow(; plane, flow_direction, density, particle_spacing, pressure=0.0,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
if open_boundary_layers <= 0
Expand All @@ -91,13 +91,13 @@ struct InFlow{NDIMS, IC, S, ZO, ZW, FD}
# Sample particles in boundary zone
if isnothing(initial_condition) && isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
density,
density, pressure,
direction=-flow_direction_,
n_extrude=open_boundary_layers)
elseif !isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
particle_spacing,
density,
density, pressure,
direction=-flow_direction_,
n_extrude=open_boundary_layers)
end
Expand Down Expand Up @@ -214,7 +214,7 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD}
zone_width :: ZW
flow_direction :: FD

function OutFlow(; plane, flow_direction, density, particle_spacing,
function OutFlow(; plane, flow_direction, density, particle_spacing, pressure=0.0,
initial_condition=nothing, extrude_geometry=nothing,
open_boundary_layers::Integer)
if open_boundary_layers <= 0
Expand All @@ -227,11 +227,11 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD}
# Sample particles in boundary zone
if isnothing(initial_condition) && isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing,
density,
density, pressure,
direction=flow_direction_,
n_extrude=open_boundary_layers)
elseif !isnothing(extrude_geometry)
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry;
initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; pressure,
particle_spacing, density,
direction=-flow_direction_,
n_extrude=open_boundary_layers)
Expand Down
171 changes: 148 additions & 23 deletions src/schemes/boundary/open_boundary/system.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D
RD, B} <: System{NDIMS, IC}
initial_condition :: IC
fluid_system :: FS
smoothing_length :: ELTYPE
mass :: ARRAY1D # Array{ELTYPE, 1}: [particle]
density :: ARRAY1D # Array{ELTYPE, 1}: [particle]
volume :: ARRAY1D # Array{ELTYPE, 1}: [particle]
Expand All @@ -54,9 +55,9 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D

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))
reference_velocity=nothing,
reference_pressure=nothing,
reference_density=nothing)
(; initial_condition) = boundary_zone

buffer = SystemBuffer(nparticles(initial_condition), buffer_size)
Expand All @@ -66,7 +67,7 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D
NDIMS = ndims(initial_condition)
ELTYPE = eltype(initial_condition)

if !(reference_velocity isa Function ||
if !(reference_velocity isa Function || isnothing(reference_velocity) ||
(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, " *
Expand All @@ -76,25 +77,31 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D
reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS))
end

if !(reference_pressure isa Function || reference_pressure isa Real)
if !(reference_pressure isa Function || reference_pressure isa Real ||
isnothing(reference_pressure))
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)
if !(reference_density isa Function || reference_density isa Real ||
isnothing(reference_density))
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

if isnothing(reference_pressure)
pressure = copy(initial_condition.pressure)
else
pressure = [reference_pressure_(initial_condition.coordinates[:, i], 0.0)
for i in eachparticle(initial_condition)]
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)

Expand All @@ -107,7 +114,8 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D
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,
typeof(buffer)}(initial_condition, fluid_system,
fluid_system.smoothing_length, mass, density, volume,
pressure, characteristics, previous_characteristics,
sound_speed, boundary_zone, flow_direction_,
reference_velocity_, reference_pressure_,
Expand All @@ -128,6 +136,10 @@ end

function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem)
@nospecialize system # reduce precompilation time
(; reference_velocity, reference_pressure, reference_density) = system
v_ref = isnothing(reference_velocity(0, 0)) ? "-" : string(nameof(reference_velocity))
p_ref = isnothing(reference_pressure(0, 0)) ? "-" : string(nameof(reference_pressure))
rho_ref = isnothing(reference_density(0, 0)) ? "-" : string(nameof(reference_density))

if get(io, :compact, false)
show(io, system)
Expand All @@ -138,14 +150,18 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem)
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, "prescribed velocity", v_ref)
summary_line(io, "prescribed pressure", p_ref)
summary_line(io, "prescribed density", rho_ref)
summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3))
summary_footer(io)
end
end

@inline function smoothing_kernel(system::OpenBoundarySPHSystem, distance)
smoothing_kernel(system.fluid_system, distance)
end

function reset_callback_flag!(system::OpenBoundarySPHSystem)
system.update_callback_used[] = false

Expand All @@ -167,8 +183,7 @@ end
end

@inline function update_quantities!(system::OpenBoundarySPHSystem, v, u, t)
(; density, pressure, characteristics, flow_direction, sound_speed,
reference_velocity, reference_pressure, reference_density) = system
(; density, pressure, characteristics, flow_direction, sound_speed) = system

# Update quantities based on the characteristic variables
@threaded system for particle in each_moving_particle(system)
Expand All @@ -178,13 +193,12 @@ end
J2 = characteristics[2, particle]
J3 = characteristics[3, particle]

rho_ref = reference_density(particle_position, t)
rho_ref, p_ref, v_ref = reference_values(v, system, particle, 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

Expand All @@ -202,11 +216,92 @@ function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi,
throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`"))
end

# interpolate_reference_values!(system, system.fluid_system, v, u, v_ode, u_ode, semi, t)

@trixi_timeit timer() "evaluate characteristics" evaluate_characteristics!(system, v, u,
v_ode, u_ode,
semi, t)
end

# DualSPHysics (Tafuni et al.) interpolate the boundary values directly from the interior
# fluid domain to the boundary particles by mirroring the boundary particles into
# the fluid domain with ghost nodes.
# Thus, they have to do a taylor series approximation `f_0 = f_k + (r_0 - r_k) ⋅ ∇f_k`,
# where `∇f_k` is the corrected gradient calculated at the ghost nodes.
#
# The method of characteristic do not calculated the boundary values directly but evaluate
# the characteristic variables with reference values `rho_ref`, `v_ref` and `p_ref`.
# Thus, it should be possible to interpolate the `ref` values directly at the ghost node position
# without any corrections.
#
# Tafuni et al. (2018)
# "A versatile algorithm for the treatment of open boundary conditions in Smoothed particle hydrodynamics GPU models"
# https://doi.org/10.1016/j.cma.2018.08.004
#
function interpolate_reference_values!(system, fluid_system::FluidSystem,
v, u, v_ode, u_ode, semi, t)
(; volume, initial_condition, boundary_zone) = system
(; density, pressure, velocity) = initial_condition

set_zero!(volume)
set_zero!(density)
set_zero!(pressure)
set_zero!(velocity)

v_fluid_system = wrap_v(v_ode, fluid_system, semi)
u_fluid_system = wrap_u(u_ode, fluid_system, semi)

fluid_coords = current_coordinates(u_fluid_system, fluid_system)
nhs = get_neighborhood_search(system, fluid_system, semi)

for particle in each_moving_particle(system)
particle_coords = current_coords(u, system, particle)

ghost_node_position = mirror_position(particle_coords, boundary_zone)

# Check the neighboring fluid particles whether they're entering the boundary zone
for neighbor in PointNeighbors.eachneighbor(ghost_node_position, nhs)
fluid_coords = current_coords(u_fluid_system, fluid_system, neighbor)

pos_diff = ghost_node_position - fluid_coords
distance2 = dot(pos_diff, pos_diff)
if distance2 <= nhs.search_radius^2
distance = sqrt(distance2)
rho_b = particle_density(v_fluid_system, fluid_system, neighbor)

p_b = particle_pressure(v_fluid_system, fluid_system, neighbor)

v_fluid = current_velocity(v_fluid_system, fluid_system, neighbor)

kernel_weight = smoothing_kernel(fluid_system, distance)

pressure[particle] += p_b * kernel_weight

density[particle] += rho_b * kernel_weight

for dim in 1:ndims(system)
velocity[dim, particle] += v_fluid[dim] * kernel_weight
end

volume[particle] += kernel_weight
end
end

end

for particle in each_moving_particle(system)
if volume[particle] > eps()
pressure[particle] /= volume[particle]
density[particle] /= volume[particle]
for dim in 1:ndims(system)
velocity[dim, particle] /= volume[particle]
end
end
end

return system
end

# ==== Characteristics
# J1: Associated with convection and entropy and propagates at flow velocity.
# J2: Propagates downstream to the local flow
Expand Down Expand Up @@ -269,8 +364,7 @@ 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
(; volume, sound_speed, characteristics, flow_direction) = system

v_neighbor_system = wrap_v(v_ode, neighbor_system, semi)
u_neighbor_system = wrap_u(u_ode, neighbor_system, semi)
Expand All @@ -287,18 +381,17 @@ function evaluate_characteristics!(system, neighbor_system::FluidSystem,

# 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)

rho_ref, p_ref, v_ref = reference_values(v, system, particle, 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))
velocity_term = rho_b * sound_speed * (dot(v_b - v_ref, flow_direction))

kernel_ = smoothing_kernel(neighbor_system, distance)

Expand Down Expand Up @@ -466,6 +559,38 @@ function write_u0!(u0, system::OpenBoundarySPHSystem)
return u0
end

function reference_values(v, system, particle, position, t)
(; reference_velocity, reference_pressure, reference_density,
initial_condition) = system

rho_ref = reference_density(position, t)
p_ref = reference_pressure(position, t)
v_ref = reference_velocity(position, t)

if isnothing(rho_ref)
rho_ref = initial_condition.density[particle]
end
if isnothing(p_ref)
p_ref = initial_condition.pressure[particle]
end
if isnothing(v_ref)
v_ref = current_velocity(v, system, particle)
end

return rho_ref, p_ref, v_ref
end

function mirror_position(particle_coords, boundary_zone)
particle_position = particle_coords - boundary_zone.zone_origin
dist = dot(particle_position, boundary_zone.flow_direction)

return particle_coords - 2 * dist * boundary_zone.flow_direction
end

function wrap_reference_function(::Nothing, ::Val)
return (coords, t) -> nothing
end

function wrap_reference_function(function_::Function, ::Val)
# Already a function
return function_
Expand Down
Loading
Loading