From 91ee930cfe0d13d24a976577a1eaddb451b7f481 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 1 Jul 2024 19:29:06 +0200 Subject: [PATCH] determine velocity from conservation equation --- examples/fluid/pipe_flow_2d.jl | 4 +- .../boundary/open_boundary/boundary_zones.jl | 12 +-- src/schemes/boundary/open_boundary/system.jl | 83 ++++++++++++++----- src/schemes/boundary/rhs.jl | 45 ++++++++++ src/setups/rectangular_tank.jl | 2 +- src/visualization/write2vtk.jl | 15 ++-- 6 files changed, 120 insertions(+), 41 deletions(-) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 3d642696e..802964cd3 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -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]]), @@ -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 diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index 586d565b1..66f9087bc 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 8c580cf44..301ed46cd 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -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] @@ -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) @@ -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, " * @@ -76,7 +77,8 @@ 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")) @@ -84,7 +86,8 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D 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")) @@ -92,9 +95,13 @@ struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D 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) @@ -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_, @@ -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) @@ -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 @@ -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) @@ -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 @@ -269,8 +283,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) @@ -287,18 +300,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) @@ -466,6 +478,31 @@ 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 wrap_reference_function(::Nothing, ::Val) + return (coords, t) -> nothing +end + function wrap_reference_function(function_::Function, ::Val) # Already a function return function_ diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index d9cfb1779..f3cfc06c6 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -7,6 +7,51 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end +# Interaction of outlet open boundary with other systems +function interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, neighborhood_search, + particle_system::OpenBoundarySPHSystem{<:OutFlow}, + neighbor_system::FluidSystem) + (; fluid_system, sound_speed) = particle_system + + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + for_particle_neighbor(particle_system, neighbor_system, + system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + rho_a = particle_density(v_particle_system, particle_system, particle) + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + + p_a = particle_pressure(v_particle_system, particle_system, particle) + p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) + + m_a = hydrodynamic_mass(particle_system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance) + + dv_pressure = pressure_acceleration(fluid_system, fluid_system, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, nothing) + + dv_viscosity_ = dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + + for i in 1:ndims(particle_system) + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + end + end + + return dv +end + # For dummy particles with `ContinuityDensity`, solve the continuity equation function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index 70b179e3b..5874ee86b 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -150,7 +150,7 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real} n_particles_per_dim) boundary = InitialCondition(coordinates=boundary_coordinates, - velocity=boundary_velocities, + velocity=boundary_velocities, pressure=first(pressure), mass=boundary_masses, density=boundary_densities, particle_spacing=boundary_spacing) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index c316e8c9c..8e4a5e141 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -301,8 +301,6 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d 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) @@ -310,13 +308,14 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data 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)) + for particle in eachparticle(system) + particle_coords = current_coords(u, system, particle) + rho_ref, p_ref, v_ref = reference_values(v, system, particle, particle_coords, t) - vtk["prescribed_velocity"] = stack(reference_velocity.(coords, t)) - vtk["prescribed_density"] = reference_density.(coords, t) - vtk["prescribed_pressure"] = reference_pressure.(coords, t) + vtk["prescribed_density"] = rho_ref + vtk["prescribed_pressure"] = p_ref + vtk["prescribed_velocity"] = v_ref + end if write_meta_data vtk["boundary_zone"] = type2string(system.boundary_zone)