diff --git a/docs/src/refs.bib b/docs/src/refs.bib index b857c7a128..2f2a636dd1 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -406,6 +406,19 @@ @Article{Li1996 publisher = {Elsevier BV}, } +@Article{Liu2006, + author = {Liu, M.B. and Liu, G.R.}, + journal = {Applied Numerical Mathematics}, + title = {Restoring particle consistency in smoothed particle hydrodynamics}, + year = {2006}, + issn = {0168-9274}, + month = jan, + number = {1}, + pages = {19--36}, + volume = {56}, + doi = {10.1016/j.apnum.2005.02.012}, + publisher = {Elsevier BV}, +} @Article{Molteni2009, author = {Molteni, Diego and Colagrossi, Andrea}, @@ -561,6 +574,19 @@ @Article{Negi2020 publisher = {Elsevier BV}, } +@Article{Negi2022, + author = {Negi, Pawan and Ramachandran, Prabhu}, + journal = {Physics of Fluids}, + title = {How to train your solver: Verification of boundary conditions for smoothed particle hydrodynamics}, + year = {2022}, + issn = {1089-7666}, + month = nov, + number = {11}, + volume = {34}, + doi = {10.1063/5.0126234}, + publisher = {AIP Publishing}, +} + @Article{O’Connor2021, author = {O’Connor, Joseph and Rogers, Benedict D.}, journal = {Journal of Fluids and Structures}, diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 7567c8517c..3544ac7a8d 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -56,7 +56,7 @@ pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers NDIMS = ndims(pipe.fluid) -n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]^(NDIMS - 1) +n_buffer_particles = 5 * pipe.n_particles_per_dimension[2]^(NDIMS - 1) # ========================================================================================== # ==== Fluid diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 5e79798f35..bad506163b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -81,6 +81,7 @@ export tensile_instability_control export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing, BoundaryModelLastiwka, BoundaryModelTafuni, BernoulliPressureExtrapolation +export FirstOrderMirroring, ZerothOrderMirroring, SimpleMirroring export HertzContactModel, LinearContactModel export BoundaryMovement export examples_dir, validation_dir diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 8dacc446e7..addc29210a 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -1,23 +1,29 @@ @doc raw""" - BoundaryModelLastiwka(; extrapolate_reference_values::Bool=false) + BoundaryModelLastiwka(; extrapolate_reference_values=nothing) Boundary model for [`OpenBoundarySPHSystem`](@ref). This model uses the characteristic variables to propagate the appropriate values -to the outlet or inlet and have been proposed by Lastiwka et al. (2009). +to the outlet or inlet and was proposed by Lastiwka et al. (2009). It requires a specific flow direction to be passed to the [`BoundaryZone`](@ref). For more information about the method see [description below](@ref method_of_characteristics). # Keywords -- `extrapolate_reference_values=false`: If `true`, the reference values are extrapolated - from the fluid domain to the boundary particles. This is useful for open boundaries where - the reference values are not known a priori. +- `extrapolate_reference_values=nothing`: If one of the following mirroring methods is selected, + the reference values are extrapolated from the fluid domain to the boundary particles: + - [`ZerothOrderMirroring`](@ref) + - [`FirstOrderMirroring`](@ref) + - [`SimpleMirroring`](@ref) + **Note:** This feature is experimental and has not been fully validated yet. As of now, we are not aware of any published literature supporting its use. + Note that even without this extrapolation feature, + the reference values don't need to be prescribed - they're computed from the characteristics. """ -struct BoundaryModelLastiwka - extrapolate_reference_values::Bool - function BoundaryModelLastiwka(; extrapolate_reference_values::Bool=false) - return new{}(extrapolate_reference_values) +struct BoundaryModelLastiwka{T} + extrapolate_reference_values::T + + function BoundaryModelLastiwka(; extrapolate_reference_values=nothing) + return new{typeof(extrapolate_reference_values)}(extrapolate_reference_values) end end @@ -32,13 +38,14 @@ end sound_speed = system_sound_speed(fluid_system) - if boundary_model.extrapolate_reference_values + if !isnothing(boundary_model.extrapolate_reference_values) (; prescribed_pressure, prescribed_velocity, prescribed_density) = cache v_fluid = wrap_v(v_ode, fluid_system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) @trixi_timeit timer() "extrapolate and correct values" begin - extrapolate_values!(system, v, v_fluid, u, u_fluid, semi, t; + extrapolate_values!(system, boundary_model.extrapolate_reference_values, + v, v_fluid, u, u_fluid, semi, t; prescribed_pressure, prescribed_velocity, prescribed_density) end @@ -71,6 +78,13 @@ end end end + if boundary_zone.average_inflow_velocity + # Even if the velocity is prescribed, this boundary model computes the velocity for each particle individually. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_model, boundary_zone, semi) + end + return system end @@ -218,3 +232,26 @@ end return characteristics end + +function average_velocity!(v, u, system, ::BoundaryModelLastiwka, boundary_zone, semi) + # Only apply averaging at the inflow + return v +end + +function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone{InFlow}, + semi) + + # Division inside the `sum` closure to maintain GPU compatibility + avg_velocity = sum(each_moving_particle(system)) do particle + return current_velocity(v, system, particle) / system.buffer.active_particle_count[] + end + + @threaded semi for particle in each_moving_particle(system) + # Set the velocity of the ghost node to the average velocity of the fluid domain + for dim in eachindex(avg_velocity) + @inbounds v[dim, particle] = avg_velocity[dim] + end + end + + return v +end diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 2eb7b286de..3d9fde0af7 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -1,41 +1,141 @@ +""" + FirstOrderMirroring(; firstorder_tolerance=1f-3) + +Fluid properties are extrapolated onto ghost nodes using the method proposed by [Liu and Liu (2006)](@cite Liu2006), +to extend the gradient into the boundary zone. + +# Keywords +- `firstorder_tolerance=1f-3`: If the determinant of the correction matrix is smaller than this value, + the method falls back to [`ZerothOrderMirroring`](@ref). +""" +struct FirstOrderMirroring{ELTYPE} + firstorder_tolerance::ELTYPE + + function FirstOrderMirroring(; firstorder_tolerance::Real=1.0f-3) + return new{typeof(firstorder_tolerance)}(firstorder_tolerance) + end +end + +""" + SimpleMirroring(; firstorder_tolerance=1f-3)) + +This method is similar to [`FirstOrderMirroring`](@ref), but does not use +the corrected gradient as proposed by [Negi et al. (2022)](@cite Negi2022). + +# Keywords +- `firstorder_tolerance=1f-3`: If the determinant of the correction matrix is smaller than this value, + the method falls back to [`ZerothOrderMirroring`](@ref). +""" +struct SimpleMirroring{ELTYPE} + firstorder_tolerance::ELTYPE + + function SimpleMirroring(; firstorder_tolerance::Real=1.0f-3) + return new{typeof(firstorder_tolerance)}(firstorder_tolerance) + end +end + +""" + ZerothOrderMirroring() + +Fluid properties are interpolated onto ghost nodes using Shepard interpolation. +(See slide 6 from the 4th DualSPHysics Users Workshop: +[Tafuni, Lisbon 2018](https://dual.sphysics.org/4thusersworkshop/data/_uploaded/PDF_Talks_4thWorkshop/Tafuni_Lisbon2018.pdf)). +The position of the ghost nodes is obtained by mirroring the boundary particles +into the fluid along a direction that is normal to the open boundary. +The interpolated values at the ghost nodes are then assigned to the corresponding boundary particles. +""" +struct ZerothOrderMirroring end + @doc raw""" - BoundaryModelTafuni() + BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) Boundary model for the `OpenBoundarySPHSystem`. This model implements the method of [Tafuni et al. (2018)](@cite Tafuni2018) to extrapolate the properties from the fluid domain to the buffer zones (inflow and outflow) using ghost nodes. The position of the ghost nodes is obtained by mirroring the boundary particles into the fluid along a direction that is normal to the open boundary. +Fluid properties are then interpolated at these ghost node positions using surrounding fluid particles. +The values are then mirrored back to the boundary particles. +We provide three different mirroring methods: + - [`ZerothOrderMirroring`](@ref): Uses a Shepard interpolation to interpolate the values. + - [`FirstOrderMirroring`](@ref): Uses a first order correction based on the gradient of the interpolated values . + - [`SimpleMirroring`](@ref): Similar to the first order mirroring, but does not use the gradient of the interpolated values. """ -struct BoundaryModelTafuni end +struct BoundaryModelTafuni{MM} + mirror_method::MM +end + +function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) + return BoundaryModelTafuni(mirror_method) +end + +function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, + v, u, v_ode, u_ode, semi, t) + (; reference_pressure, reference_density, reference_velocity, boundary_zone, + pressure, density, cache) = system + (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache -function update_boundary_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, - semi, t) @trixi_timeit timer() "extrapolate and correct values" begin fluid_system = corresponding_fluid_system(system, semi) - v_open_boundary = wrap_v(v_ode, system, semi) v_fluid = wrap_v(v_ode, fluid_system, semi) - u_open_boundary = wrap_u(u_ode, system, semi) u_fluid = wrap_u(u_ode, fluid_system, semi) - extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, u_fluid, - semi, t; system.cache...) + extrapolate_values!(system, boundary_model.mirror_method, v, v_fluid, + u, u_fluid, semi, t; prescribed_pressure, + prescribed_density, prescribed_velocity) + end + + if !prescribed_velocity && boundary_zone.average_inflow_velocity + # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. + # Thus, turbulent flows near the inflow can lead to a non-uniform buffer particle distribution, + # resulting in a potential numerical instability. Averaging mitigates these effects. + average_velocity!(v, u, system, boundary_zone, semi) + end + + if prescribed_pressure + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + pressure[particle] = reference_value(reference_pressure, pressure[particle], + particle_coords, t) + end + end + + if prescribed_density + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + density[particle] = reference_value(reference_density, density[particle], + particle_coords, t) + end + end + + if prescribed_velocity + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + v_particle = current_velocity(v, system, particle) + + v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) + + for dim in eachindex(v_ref) + @inbounds v[dim, particle] = v_ref[dim] + end + end end end update_final!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) = system -function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, u_fluid, +function extrapolate_values!(system, + mirror_method::Union{FirstOrderMirroring, SimpleMirroring}, + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi, t; prescribed_density=false, prescribed_pressure=false, prescribed_velocity=false) - (; pressure, density, boundary_zone, reference_density, - reference_velocity, reference_pressure) = system + (; pressure, density, boundary_zone) = system fluid_system = corresponding_fluid_system(system, semi) - state_equation = system_state_equation(fluid_system) - # Static indices to avoid allocations two_to_end = SVector{ndims(system)}(2:(ndims(system) + 1)) @@ -57,13 +157,13 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, correction_matrix = Ref(zero(SMatrix{ndims(system) + 1, ndims(system) + 1, eltype(system)})) - extrapolated_density_correction = Ref(zero(SVector{ndims(system) + 1, + interpolated_density_correction = Ref(zero(SVector{ndims(system) + 1, eltype(system)})) - extrapolated_pressure_correction = Ref(zero(SVector{ndims(system) + 1, + interpolated_pressure_correction = Ref(zero(SVector{ndims(system) + 1, eltype(system)})) - extrapolated_velocity_correction = Ref(zero(SMatrix{ndims(system), + interpolated_velocity_correction = Ref(zero(SMatrix{ndims(system), ndims(system) + 1, eltype(system)})) @@ -74,13 +174,7 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, m_b = hydrodynamic_mass(fluid_system, neighbor) rho_b = current_density(v_fluid, fluid_system, neighbor) pressure_b = current_pressure(v_fluid, fluid_system, neighbor) - v_b_ = current_velocity(v_fluid, fluid_system, neighbor) - - # Project `v_b_` on the normal direction of the boundary zone (only for inflow boundaries). - # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: - # "Because flow from the inlet interface occurs perpendicular to the boundary, - # only this component of interpolated velocity is kept [...]" - v_b = project_velocity_on_plane_normal(v_b_, boundary_zone) + v_b = current_velocity(v_fluid, fluid_system, neighbor) kernel_value = smoothing_kernel(fluid_system, distance, particle) grad_kernel = smoothing_kernel_grad(fluid_system, pos_diff, distance, @@ -92,84 +186,249 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, correction_matrix[] += L - # For a WCSPH system, the pressure is determined by the state equation if it is not prescribed - if !prescribed_pressure && !(fluid_system isa WeaklyCompressibleSPHSystem) - extrapolated_pressure_correction[] += pressure_b * R + if !prescribed_pressure + interpolated_pressure_correction[] += pressure_b * R end - if !prescribed_velocity - extrapolated_velocity_correction[] += v_b * R' + if !prescribed_density + interpolated_density_correction[] += rho_b * R end - if !prescribed_density - extrapolated_density_correction[] += rho_b * R + if !prescribed_velocity + interpolated_velocity_correction[] += v_b * R' end end - # See also `correction_matrix_inversion_step!` for an explanation - if abs(det(correction_matrix[])) < 1.0f-9 - L_inv = typeof(correction_matrix[])(I) - else + pos_diff = particle_coords - ghost_node_position + + # Mirror back the ghost node values to the boundary particles + if abs(det(correction_matrix[])) >= mirror_method.firstorder_tolerance L_inv = inv(correction_matrix[]) - end - pos_diff = particle_coords - ghost_node_position + # Pressure + if !prescribed_pressure + first_order_scalar_extrapolation!(pressure, particle, L_inv, + interpolated_pressure_correction[], + two_to_end, pos_diff, mirror_method) + end - # In Negi et al. (2020) https://doi.org/10.1016/j.cma.2020.113119, - # they have modified the equation for extrapolation to - # - # f_0 = f_k - (r_0 - r_k) ⋅ ∇f_k - # - # in order to get zero gradient at the outlet interface. - # Note: This modification is mentioned here for reference only and is NOT applied in this implementation. - if prescribed_velocity - v_particle = current_velocity(v_open_boundary, system, particle) - v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) - @inbounds for dim in eachindex(v_ref) - v_open_boundary[dim, particle] = v_ref[dim] + # Density + if !prescribed_density + first_order_scalar_extrapolation!(density, particle, L_inv, + interpolated_density_correction[], + two_to_end, pos_diff, mirror_method) + end + + # Velocity + if !prescribed_velocity + first_order_velocity_extrapolation!(v_open_boundary, particle, L_inv, + interpolated_velocity_correction[], + two_to_end, pos_diff, mirror_method) + + # Project the velocity on the normal direction of the boundary zone (only for inflow boundaries). + # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: + # "Because flow from the inlet interface occurs perpendicular to the boundary, + # only this component of interpolated velocity is kept [...]" + project_velocity_on_plane_normal!(v_open_boundary, system, particle, + boundary_zone) end - else - @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * extrapolated_velocity_correction[][dim, :] - df_v = f_v[two_to_end] - v_open_boundary[dim, particle] = f_v[1] + dot(pos_diff, df_v) + # No else: `correction_matrix[][1, 1] <= eps()` means no fluid neighbors + # and thus no reliable interpolation, so boundary values remain at their current state. + elseif correction_matrix[][1, 1] > eps() + # Determinant is small, fallback to zero-th order mirroring + shepard_coefficient = correction_matrix[][1, 1] + + # Pressure + if !prescribed_pressure + # Only the first entry is used, as the subsequent entries represent gradient + # components that are not required for zeroth-order interpolation. + interpolated_pressure = first(interpolated_pressure_correction[]) + zeroth_order_scalar_extrapolation!(pressure, particle, shepard_coefficient, + interpolated_pressure) + end + + # Density + if !prescribed_density + # Only the first entry is used, as the subsequent entries represent gradient + # components that are not required for zeroth-order interpolation. + interpolated_density = first(interpolated_density_correction[]) + zeroth_order_scalar_extrapolation!(density, particle, shepard_coefficient, + interpolated_density) + end + + # Velocity + if !prescribed_velocity + # Only the first column is used, as the subsequent entries represent gradient + # components that are not required for zeroth-order interpolation. + interpolated_velocity = interpolated_velocity_correction[][:, 1] + zeroth_order_velocity_interpolation!(v_open_boundary, particle, + shepard_coefficient, + interpolated_velocity) + + # Project the velocity on the normal direction of the boundary zone (only for inflow boundaries). + # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: + # "Because flow from the inlet interface occurs perpendicular to the boundary, + # only this component of interpolated velocity is kept [...]" + project_velocity_on_plane_normal!(v_open_boundary, system, particle, + boundary_zone) end end + end - if prescribed_density - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) - else - f_d = L_inv * extrapolated_density_correction[] - df_d = f_d[two_to_end] + return system +end + +function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, + v_open_boundary, v_fluid, u_open_boundary, u_fluid, semi, t; + prescribed_density=false, prescribed_pressure=false, + prescribed_velocity=false) + (; pressure, density, boundary_zone) = system + + fluid_system = corresponding_fluid_system(system, semi) - density[particle] = f_d[1] + dot(pos_diff, df_d) + # Use the fluid-fluid nhs, since the boundary particles are mirrored into the fluid domain + nhs = get_neighborhood_search(fluid_system, fluid_system, semi) + + fluid_coords = current_coordinates(u_fluid, fluid_system) + + # We cannot use `foreach_point_neighbor` here because we are looking for neighbors + # of the ghost node positions of each particle. + # We can do this because we require the neighborhood search to support querying neighbors + # of arbitrary positions (see `PointNeighbors.requires_update`). + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u_open_boundary, system, particle) + ghost_node_position = mirror_position(particle_coords, boundary_zone) + + # Use `Ref` to ensure the variables are accessible and mutable within the closure below + # (see https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured). + shepard_coefficient = Ref(zero(eltype(system))) + interpolated_density = Ref(zero(eltype(system))) + interpolated_pressure = Ref(zero(eltype(system))) + interpolated_velocity = Ref(zero(particle_coords)) + + # TODO: Not public API + PointNeighbors.foreach_neighbor(fluid_coords, nhs, particle, ghost_node_position, + nhs.search_radius) do particle, neighbor, pos_diff, + distance + m_b = hydrodynamic_mass(fluid_system, neighbor) + rho_b = current_density(v_fluid, fluid_system, neighbor) + volume_b = m_b / rho_b + pressure_b = current_pressure(v_fluid, fluid_system, neighbor) + v_b = current_velocity(v_fluid, fluid_system, neighbor) + + W_ab = smoothing_kernel(fluid_system, distance, particle) + + shepard_coefficient[] += volume_b * W_ab + + if !prescribed_pressure + interpolated_pressure[] += pressure_b * volume_b * W_ab + end + + if !prescribed_density + interpolated_density[] += rho_b * volume_b * W_ab + end + + if !prescribed_velocity + interpolated_velocity[] += v_b * volume_b * W_ab + end end - if prescribed_pressure - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) - elseif fluid_system isa WeaklyCompressibleSPHSystem - # For a WCSPH system, the pressure is determined by the state equation - # if it is not prescribed - pressure[particle] = state_equation(density[particle]) - else - f_d = L_inv * extrapolated_pressure_correction[] - df_d = f_d[two_to_end] - - pressure[particle] = f_d[1] + dot(pos_diff, df_d) + if shepard_coefficient[] > eps() + pos_diff = particle_coords - ghost_node_position + + if !prescribed_pressure + zeroth_order_scalar_extrapolation!(pressure, particle, + shepard_coefficient[], + interpolated_pressure[]) + end + + if !prescribed_density + zeroth_order_scalar_extrapolation!(density, particle, + shepard_coefficient[], + interpolated_density[]) + end + + if !prescribed_velocity + zeroth_order_velocity_interpolation!(v_open_boundary, particle, + shepard_coefficient[], + interpolated_velocity[]) + + # Project the velocity on the normal direction of the boundary zone (only for inflow boundaries). + # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: + # "Because flow from the inlet interface occurs perpendicular to the boundary, + # only this component of interpolated velocity is kept [...]" + project_velocity_on_plane_normal!(v_open_boundary, system, particle, + boundary_zone) + end end end - if !(prescribed_velocity) && boundary_zone.average_inflow_velocity - # When no velocity is prescribed at the inflow, the velocity is extrapolated from the fluid domain. - # Thus, turbulent flows near the inflow can lead to non-uniform buffer particles distribution, - # resulting in a potential numerical instability. Averaging mitigates these effects. - average_velocity!(v_open_boundary, u_open_boundary, system, boundary_zone, semi) + return system +end + +function zeroth_order_scalar_extrapolation!(values, particle, shepard_coefficient, + extrapolated_value) + values[particle] = extrapolated_value / shepard_coefficient + + return values +end + +function zeroth_order_velocity_interpolation!(v, particle, shepard_coefficient, + extrapolated_velocity) + velocity_interpolated = extrapolated_velocity / shepard_coefficient + + for dim in eachindex(velocity_interpolated) + @inbounds v[dim, particle] = velocity_interpolated[dim] end - return system + return v +end + +function first_order_scalar_extrapolation!(values, particle, L_inv, + interpolated_values_correction, + two_to_end, pos_diff, ::FirstOrderMirroring) + f_s = L_inv * interpolated_values_correction + df_s = f_s[two_to_end] # `f_s[2:end]` as SVector + + values[particle] = f_s[1] + dot(pos_diff, df_s) + + return values +end + +function first_order_scalar_extrapolation!(values, particle, L_inv, + interpolated_values_correction, + two_to_end, pos_diff, ::SimpleMirroring) + f_s = L_inv * interpolated_values_correction + + values[particle] = f_s[1] + + return values +end + +function first_order_velocity_extrapolation!(v, particle, L_inv, + interpolated_velocity_correction, + two_to_end, pos_diff, ::FirstOrderMirroring) + for dim in eachindex(pos_diff) + @inbounds f_v = L_inv * interpolated_velocity_correction[dim, :] + df_v = f_v[two_to_end] # `f_v[2:end]` as SVector + + @inbounds v[dim, particle] = f_v[1] + dot(pos_diff, df_v) + end + + return v +end + +function first_order_velocity_extrapolation!(v, particle, L_inv, + interpolated_velocity_correction, + two_to_end, pos_diff, ::SimpleMirroring) + for dim in eachindex(pos_diff) + @inbounds f_v = L_inv * interpolated_velocity_correction[dim, :] + + @inbounds v[dim, particle] = f_v[1] + end + + return v end function correction_arrays(W_ab, grad_W_ab, pos_diff::SVector{3}, rho_b, m_b) @@ -241,28 +500,36 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se reinterpret(reshape, SVector{ndims(system), eltype(u)}, active_coordinates(u, system))) + # Division inside the `sum` closure to maintain GPU compatibility avg_velocity = sum(candidates) do particle - return current_velocity(v, system, particle) + return current_velocity(v, system, particle) / length(candidates) end - avg_velocity /= length(candidates) - @threaded semi for particle in each_moving_particle(system) # Set the velocity of the ghost node to the average velocity of the fluid domain - @inbounds for dim in eachindex(avg_velocity) - v[dim, particle] = avg_velocity[dim] + for dim in eachindex(avg_velocity) + @inbounds v[dim, particle] = avg_velocity[dim] end end return v end -project_velocity_on_plane_normal(vel, boundary_zone) = vel +project_velocity_on_plane_normal!(v, system, particle, boundary_zone) = v -function project_velocity_on_plane_normal(vel, boundary_zone::BoundaryZone{InFlow}) - # Project `v_b` on the normal direction of the boundary zone +function project_velocity_on_plane_normal!(v, system, particle, + boundary_zone::BoundaryZone{InFlow}) + # Project `vel` on the normal direction of the boundary zone # See https://doi.org/10.1016/j.jcp.2020.110029 Section 3.3.: # "Because flow from the inlet interface occurs perpendicular to the boundary, # only this component of interpolated velocity is kept [...]" - return dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal + v_particle = current_velocity(v, system, particle) + v_particle_projected = dot(v_particle, boundary_zone.plane_normal) * + boundary_zone.plane_normal + + for dim in eachindex(v_particle) + @inbounds v[dim, particle] = v_particle_projected[dim] + end + + return v end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 14034fedc2..19bf85792b 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -27,10 +27,10 @@ Open boundary system for in- and outflow particles. or a scalar for a constant density over all particles. !!! note "Note" - When using the [`BoundaryModelTafuni()`](@ref), the reference values (`reference_velocity`, - `reference_pressure`, `reference_density`) can also be set to `nothing` - since this model allows for either assigning physical quantities a priori or extrapolating them - from the fluid domaim to the buffer zones (inflow and outflow) using ghost nodes. + The reference values (`reference_velocity`, `reference_pressure`, `reference_density`) + can also be set to `nothing`. + In this case, they will either be extrapolated from the fluid domain ([BoundaryModelTafuni](@ref BoundaryModelTafuni)) + or evolved using the characteristic flow variables ([BoundaryModelLastiwka](@ref BoundaryModelLastiwka)). !!! warning "Experimental Implementation" This is an experimental feature and may change in future releases. @@ -88,9 +88,6 @@ function OpenBoundarySPHSystem(boundary_zone::BoundaryZone; reference_density=nothing) (; initial_condition) = boundary_zone - check_reference_values!(boundary_model, reference_density, reference_pressure, - reference_velocity) - buffer = SystemBuffer(nparticles(initial_condition), buffer_size) initial_condition = allocate_buffer(initial_condition, buffer) @@ -270,6 +267,8 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ u = wrap_u(u_ode, system, semi) v = wrap_v(v_ode, system, semi) + @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, 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 boundary quantities" update_boundary_quantities!(system, @@ -279,8 +278,6 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ u_ode, semi, t) - @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) - return system end @@ -351,6 +348,9 @@ function check_domain!(system, v, u, v_ode, u_ode, semi) update_system_buffer!(system.buffer, semi) update_system_buffer!(fluid_system.buffer, semi) + # Since particles have been transferred, the neighborhood searches must be updated + update_nhs!(semi, u_ode) + return system end @@ -484,22 +484,6 @@ end # instead of using the method of characteristics reference_value(value::Nothing, quantity, position, t) = quantity -function check_reference_values!(boundary_model, reference_density, reference_pressure, - reference_velocity) - return boundary_model -end - -function check_reference_values!(boundary_model::BoundaryModelLastiwka, - reference_density, reference_pressure, reference_velocity) - boundary_model.extrapolate_reference_values && return boundary_model - - if any(isnothing.([reference_density, reference_pressure, reference_velocity])) - throw(ArgumentError("for `BoundaryModelLastiwka` all reference values must be specified")) - end - - return boundary_model -end - # To account for boundary effects in the viscosity term of the RHS, use the viscosity model # of the neighboring particle systems. @inline function viscosity_model(system::OpenBoundarySPHSystem, diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index a9a2e791dd..0acdcb8524 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -340,7 +340,8 @@ end "pipe_flow_2d.jl"), wcsph=true, sound_speed=20.0f0, pressure=0.0f0, - open_boundary_model=BoundaryModelTafuni(), + open_boundary_model=BoundaryModelTafuni(; + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index cc4c970cb9..6e3fb9cd6f 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -72,7 +72,8 @@ TrixiParticles.set_zero!(open_boundary.pressure) - TrixiParticles.extrapolate_values!(open_boundary, v_open_boundary, v_fluid, + TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), + v_open_boundary, v_fluid, inflow.initial_condition.coordinates, domain_fluid.coordinates, semi, 0.0; prescribed_pressure=false, @@ -168,7 +169,8 @@ TrixiParticles.set_zero!(open_boundary.pressure) - TrixiParticles.extrapolate_values!(open_boundary, v_open_boundary, v_fluid, + TrixiParticles.extrapolate_values!(open_boundary, FirstOrderMirroring(), + v_open_boundary, v_fluid, inflow.initial_condition.coordinates, domain_fluid.coordinates, semi, 0.0; prescribed_pressure=false, @@ -236,18 +238,307 @@ TrixiParticles.initialize_neighborhood_searches!(semi) v_open_boundary = zero(inflow.initial_condition.velocity) + u_open_boundary = inflow.initial_condition.coordinates v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') TrixiParticles.set_zero!(open_boundary_in.pressure) - TrixiParticles.extrapolate_values!(open_boundary_in, v_open_boundary, v_fluid, + TrixiParticles.extrapolate_values!(open_boundary_in, FirstOrderMirroring(), + v_open_boundary, v_fluid, inflow.initial_condition.coordinates, domain_fluid.coordinates, semi, 0.0) + TrixiParticles.average_velocity!(v_open_boundary, u_open_boundary, open_boundary_in, + inflow, semi) + # Since the velocity profile increases linearly in positive x-direction, # we can use the first velocity entry as a representative value. v_x_fluid_first = v_fluid[1, 1] @test all(isapprox(-v_x_fluid_first), v_open_boundary[1, :]) end + + @testset verbose=true "Mirroring Methods" begin + function mirror(pressure_function, mirror_method; + particle_spacing=0.05, domain_size=(2.0, 1.0)) + # Initialize a fluid block with pressure according to `pressure_function` + # and a adjacent inflow and outflow open boundaries to test the pressure extrapolation. + domain_fluid = RectangularShape(particle_spacing, + round.(Int, domain_size ./ particle_spacing), + (0.0, 0.0), density=1000.0, + pressure=pressure_function) + + smoothing_length = 1.2 * particle_spacing + smoothing_kernel = WendlandC2Kernel{2}() + fluid_system = EntropicallyDampedSPHSystem(domain_fluid, smoothing_kernel, + smoothing_length, 1.0) + + fluid_system.cache.density .= domain_fluid.density + + plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) + + outflow = BoundaryZone(; plane=plane_out, boundary_type=OutFlow(), + plane_normal=[-1.0, 0.0], + open_boundary_layers=10, density=1000.0, + particle_spacing) + open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=BoundaryModelTafuni(), + buffer_size=0) + + # Temporary semidiscretization just to extrapolate the pressure into the outflow system + semi = Semidiscretization(fluid_system, open_boundary_out) + TrixiParticles.initialize_neighborhood_searches!(semi) + + v_open_boundary = zero(outflow.initial_condition.velocity) + v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure') + + TrixiParticles.set_zero!(open_boundary_out.pressure) + + TrixiParticles.extrapolate_values!(open_boundary_out, mirror_method, + v_open_boundary, v_fluid, + outflow.initial_condition.coordinates, + domain_fluid.coordinates, semi, 0.0; + prescribed_pressure=false) + + plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) + + inflow = BoundaryZone(; plane=plane_in, boundary_type=InFlow(), + plane_normal=[1.0, 0.0], + open_boundary_layers=10, density=1000.0, particle_spacing) + open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, + boundary_model=BoundaryModelTafuni(), + buffer_size=0) + + # Temporary semidiscretization just to extrapolate the pressure into the outflow system + semi = Semidiscretization(fluid_system, open_boundary_in) + TrixiParticles.initialize_neighborhood_searches!(semi) + + v_open_boundary = zero(inflow.initial_condition.velocity) + + TrixiParticles.set_zero!(open_boundary_in.pressure) + + TrixiParticles.extrapolate_values!(open_boundary_in, mirror_method, + v_open_boundary, v_fluid, + inflow.initial_condition.coordinates, + domain_fluid.coordinates, semi, 0.0; + prescribed_pressure=false) + + return fluid_system, open_boundary_in, open_boundary_out, v_fluid + end + + function interpolate_pressure(mirror_method, pressure_func; particle_spacing=0.05) + # First call the function above to initialize fluid with pressure according to the function + # and then extrapolate pressure to the inflow and outflow boundary systems. + # Then, in this function, we apply an SPH interpolation on this extrapolated pressure field + # to get a continuous representation of the extrapolated pressure field to validate. + fluid_system, open_boundary_in, open_boundary_out, + v_fluid = mirror(pressure_func, mirror_method) + + p_fluid = [TrixiParticles.current_pressure(v_fluid, fluid_system, particle) + for particle in TrixiParticles.active_particles(fluid_system)] + + fluid_system.initial_condition.pressure .= p_fluid + open_boundary_in.initial_condition.pressure .= open_boundary_in.pressure + open_boundary_out.initial_condition.pressure .= open_boundary_out.pressure + + entire_domain = union(fluid_system.initial_condition, + open_boundary_in.initial_condition, + open_boundary_out.initial_condition) + + smoothing_length = 1.2 * particle_spacing + smoothing_kernel = WendlandC2Kernel{2}() + + # Use a temporary fluid system just to interpolate the pressure + interpolation_system = WeaklyCompressibleSPHSystem(entire_domain, + ContinuityDensity(), + nothing, smoothing_kernel, + smoothing_length) + interpolation_system.pressure .= entire_domain.pressure + + semi = Semidiscretization(interpolation_system) + ode = semidiscretize(semi, (0, 0)) + v_ode, u_ode = ode.u0.x + + result = interpolate_line([-0.5, 0.5], [2.5, 0.5], 50, semi, + interpolation_system, v_ode, u_ode) + + return result.pressure + end + + pressure_func(pos) = cos(2pi * pos[1]) + + # The pressures are interpolated to obtain a unified vector of length 50, + # rather than handling three separate systems with numerous particles each. + # Additionally, it facilitates plotting for test validation purposes. + pressures = interpolate_pressure.([ + SimpleMirroring(), + FirstOrderMirroring(), + ZerothOrderMirroring() + ], pressure_func) + + pressures_expected = [ + [ + -0.961368753262176, + -0.889650754605612, + -0.6905728633912162, + -0.3906982503900935, + -0.03185103264546474, + 0.333177250635236, + 0.6476757904268203, + 0.8708718636029472, + 0.9770644543539139, + 0.930127404224564, + 0.7453239657417114, + 0.45327391778794734, + 0.0957991316946515, + -0.278044550109137, + -0.6122210347366616, + -0.8549953877169134, + -0.972402370920352, + -0.9485947076098263, + -0.7852083273404225, + -0.5067072947529646, + -0.15639554738985995, + 0.216475912074678, + 0.5602009648581898, + 0.8223333154520063, + 0.9626909343468545, + 0.9626909343468544, + 0.8223333154520058, + 0.5602009648581916, + 0.2164759120746806, + -0.15639554738985995, + -0.5067072947529652, + -0.7852083273404221, + -0.9485947076098261, + -0.9724023709203522, + -0.8549953877169136, + -0.6122210347366618, + -0.2780445501091364, + 0.0957991316946491, + 0.45327391778794657, + 0.7453239657417114, + 0.9301274042245636, + 0.9770644543539139, + 0.8708718636029484, + 0.647675790426819, + 0.33317725063523634, + -0.03185103264546614, + -0.3906982503900926, + -0.6905728633912179, + -0.8896507546056125, + -0.9613687532621761 + ], + [ + 0.06915008702843263, + 1.0737102752866752, + 2.4308935813709276, + 3.074251266025277, + 3.047730973768249, + 2.5546812207882246, + 1.8827468502566729, + 1.330793458342126, + 1.0685482295301676, + 0.9337092788362713, + 0.7453239657417114, + 0.45327391778794734, + 0.0957991316946515, + -0.278044550109137, + -0.6122210347366616, + -0.8549953877169134, + -0.972402370920352, + -0.9485947076098263, + -0.7852083273404225, + -0.5067072947529646, + -0.15639554738985995, + 0.216475912074678, + 0.5602009648581898, + 0.8223333154520063, + 0.9626909343468545, + 0.9626909343468544, + 0.8223333154520058, + 0.5602009648581916, + 0.2164759120746806, + -0.15639554738985995, + -0.5067072947529652, + -0.7852083273404221, + -0.9485947076098261, + -0.9724023709203522, + -0.8549953877169136, + -0.6122210347366618, + -0.2780445501091364, + 0.0957991316946491, + 0.45327391778794657, + 0.7453239657417114, + 0.9337092788362712, + 1.0685482295301671, + 1.3307934583421233, + 1.8827468502566747, + 2.55468122078822, + 3.047730973768248, + 3.0742512660252768, + 2.4308935813709227, + 1.0737102752866743, + 0.06915008702844015 + ], + [ + -0.961368753262176, + -0.8896507546056119, + -0.690572863391216, + -0.3906982503900935, + -0.03185103264546471, + 0.333177250635236, + 0.6476223100256845, + 0.8652532179882182, + 0.9638320291788572, + 0.929273656897224, + 0.7453239657417114, + 0.45327391778794734, + 0.0957991316946515, + -0.278044550109137, + -0.6122210347366616, + -0.8549953877169134, + -0.972402370920352, + -0.9485947076098263, + -0.7852083273404225, + -0.5067072947529646, + -0.15639554738985995, + 0.216475912074678, + 0.5602009648581898, + 0.8223333154520063, + 0.9626909343468545, + 0.9626909343468544, + 0.8223333154520058, + 0.5602009648581916, + 0.2164759120746806, + -0.15639554738985995, + -0.5067072947529652, + -0.7852083273404221, + -0.9485947076098261, + -0.9724023709203522, + -0.8549953877169136, + -0.6122210347366618, + -0.2780445501091364, + 0.0957991316946491, + 0.45327391778794657, + 0.7453239657417114, + 0.9292736568972241, + 0.9638320291788572, + 0.865253217988219, + 0.6476223100256829, + 0.33317725063523645, + -0.03185103264546612, + -0.3906982503900927, + -0.6905728633912182, + -0.8896507546056123, + -0.961368753262176 + ] + ] + + @testset verbose=true "$method" for (i, method) in enumerate(("simple mirroring", + "first order mirroring", + "zeroth order mirroring")) + @test isapprox(pressures[i], pressures_expected[i]) + end + end end