From 87da03c8b49d7ea74ddb525bc5977f19f16c8344 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 2 Jul 2025 14:02:52 +0200 Subject: [PATCH 01/27] swap extrapolation --- src/schemes/boundary/open_boundary/mirroring.jl | 12 ++---------- src/schemes/boundary/open_boundary/system.jl | 3 +++ 2 files changed, 5 insertions(+), 10 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 2eb7b286de..dbe98e6a39 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -92,8 +92,7 @@ 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) + if !prescribed_pressure extrapolated_pressure_correction[] += pressure_b * R end @@ -141,19 +140,12 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, 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] - - density[particle] = f_d[1] + dot(pos_diff, df_d) + inverse_state_equation!(density, state_equation, pressure, particle) 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] diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index f299fd9286..bf07d5f252 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -450,6 +450,9 @@ end function check_reference_values!(boundary_model::BoundaryModelLastiwka, reference_density, reference_pressure, reference_velocity) + + # TODO + return boundary_model boundary_model.extrapolate_reference_values && return boundary_model if any(isnothing.([reference_density, reference_pressure, reference_velocity])) From f30f01934d1cb8eef91654e85ce21c22d4c00410 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 3 Jul 2025 12:54:52 +0200 Subject: [PATCH 02/27] swap density pressure --- src/schemes/boundary/open_boundary/mirroring.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index dbe98e6a39..5c64fdea83 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -136,13 +136,6 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, end end - if prescribed_density - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) - else - inverse_state_equation!(density, state_equation, pressure, particle) - end - if prescribed_pressure pressure[particle] = reference_value(reference_pressure, pressure[particle], particle_coords, t) @@ -152,6 +145,13 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, pressure[particle] = f_d[1] + dot(pos_diff, df_d) end + + if prescribed_density + density[particle] = reference_value(reference_density, density[particle], + particle_coords, t) + else + inverse_state_equation!(density, state_equation, pressure, particle) + end end if !(prescribed_velocity) && boundary_zone.average_inflow_velocity From 61ea4568dec9786b44033b70d37571ee0e1d5b9f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 3 Jul 2025 17:07:15 +0200 Subject: [PATCH 03/27] adapt docs --- examples/fluid/play_ground.jl | 233 ++++++++++++++++++ .../method_of_characteristics.jl | 3 +- .../boundary/open_boundary/mirroring.jl | 5 +- src/schemes/boundary/open_boundary/system.jl | 34 +-- 4 files changed, 244 insertions(+), 31 deletions(-) create mode 100644 examples/fluid/play_ground.jl diff --git a/examples/fluid/play_ground.jl b/examples/fluid/play_ground.jl new file mode 100644 index 0000000000..5bc88fb7cb --- /dev/null +++ b/examples/fluid/play_ground.jl @@ -0,0 +1,233 @@ +# Flow past a circular cylinder (vortex street), Tafuni et al. (2018). +# Other literature using this validation: +# Vacandio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998) + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +factor_d = 0.08 # Resolution in the paper is `0.01` (5M particles) + +const cylinder_diameter = 0.1 +particle_spacing = factor_d * cylinder_diameter + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 4 + +# 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 = 8 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 1.5) + +# Boundary geometry and initial fluid particle positions +domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter) + +flow_direction = [1.0, 0.0] +reynolds_number = 200 +const prescribed_velocity = 1.0 +const fluid_density = 1000.0 +sound_speed = 10 * prescribed_velocity + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2]) + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + n_layers=boundary_layers, velocity=[prescribed_velocity, 0.0], + 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 = 10 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1) + +cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2) +cylinder = SphereShape(particle_spacing, cylinder_diameter / 2, + cylinder_center, fluid_density, sphere_type=RoundSphere()) + +fluid = setdiff(pipe.fluid, cylinder) + +# ========================================================================================== +# ==== Fluid +wcsph = true + +h_factor = 1.5 +smoothing_length = h_factor * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * cylinder_diameter / reynolds_number + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +if wcsph + density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + + fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, + state_equation, smoothing_kernel, + density_diffusion=density_diffusion, + smoothing_length, viscosity=viscosity, + pressure_acceleration=tensile_instability_control, + buffer_size=n_buffer_particles) +else + state_equation = nothing + fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) +end + +# ========================================================================================== +# ==== Open Boundary +function velocity_function2d(pos, t) + return SVector(prescribed_velocity, 0.0) +end + +open_boundary_model = BoundaryModelTafuni() + +boundary_type_in = InFlow() +plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers, + density=fluid_density, particle_spacing, + boundary_type=boundary_type_in) + +reference_velocity_in = velocity_function2d +# At the inlet, neither pressure nor density are prescribed; instead, +# these values are extrapolated from the fluid domain +reference_pressure_in = nothing +reference_density_in = nothing +open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles, + reference_density=reference_density_in, + reference_pressure=reference_pressure_in, + reference_velocity=reference_velocity_in) + +boundary_type_out = OutFlow() + +plane_out = ([pipe.fluid_size[1], 0.0], [pipe.fluid_size[1], domain_size[2]]) +outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction), + open_boundary_layers, density=fluid_density, particle_spacing, + boundary_type=boundary_type_out) + +# At the outlet, we allow the flow to exit freely without imposing any boundary conditions +reference_velocity_out = nothing +reference_pressure_out = nothing +reference_density_out = nothing +open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, + boundary_model=open_boundary_model, + buffer_size=n_buffer_particles, + reference_density=reference_density_out, + reference_pressure=reference_pressure_out, + reference_velocity=reference_velocity_out) +# ========================================================================================== +# ==== Boundary +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + smoothing_kernel, smoothing_length) + +boundary_system_wall = BoundarySPHSystem(pipe.boundary, boundary_model) + +boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + viscosity=viscosity, + smoothing_kernel, smoothing_length) + +boundary_system_cylinder = BoundarySPHSystem(cylinder, boundary_model_cylinder) + +# ========================================================================================== +# ==== Postprocessing +circle = SphereShape(particle_spacing, (cylinder_diameter + particle_spacing) / 2, + cylinder_center, fluid_density, n_layers=1, + sphere_type=RoundSphere()) + +# Points for pressure interpolation, located at the wall interface +const data_points = copy(circle.coordinates) +const center = SVector(cylinder_center) + +calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[2] / (fluid_density * prescribed_velocity^2 * cylinder_diameter) +end + +calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing +function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) + force = zero(SVector{ndims(system), eltype(system)}) + + values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, + clip_negative_pressure=false) + pressure = Array(values.pressure) + + for i in axes(data_points, 2) + point = TrixiParticles.current_coords(data_points, system, i) + + # F = ∑ -p_i * A_i * n_i + force -= pressure[i] * particle_spacing .* + TrixiParticles.normalize(point - center) + end + + return 2 * force[1] / (fluid_density * prescribed_velocity^2 * cylinder_diameter) +end + +# ========================================================================================== +# ==== Simulation +min_corner = minimum(pipe.boundary.coordinates .- 5 * particle_spacing, dims=2) +max_corner = maximum(pipe.boundary.coordinates .+ 5 * particle_spacing, dims=2) +cell_list = FullGridCellList(; min_corner, max_corner) + +neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, + update_strategy=ParallelUpdate()) + +semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, + boundary_system_wall, boundary_system_cylinder; + neighborhood_search=neighborhood_search, + parallelization_backend=PolyesterBackend()) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) + +output_directory = "out" + +saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory) + +pp_callback = PostprocessCallback(; dt=0.02, + f_l=calculate_lift_force, f_d=calculate_drag_force, + output_directory, filename="resulting_force", + write_csv=true, write_file_interval=10) + +extra_callback = nothing + +callbacks = CallbackSet(info_callback, UpdateCallback(), saving_callback, + ParticleShiftingCallback(), # To obtain a near-uniform particle distribution in the wake + pp_callback, extra_callback) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-4, # 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/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 8dacc446e7..757a9328af 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -9,8 +9,7 @@ For more information about the method see [description below](@ref method_of_cha # 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. + from the fluid domain to the boundary particles. **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. """ diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 5c64fdea83..869a03bd0a 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -150,7 +150,10 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, density[particle] = reference_value(reference_density, density[particle], particle_coords, t) else - inverse_state_equation!(density, state_equation, pressure, particle) + f_d = L_inv * extrapolated_density_correction[] + df_d = f_d[two_to_end] + + density[particle] = f_d[1] + dot(pos_diff, df_d) end end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index bf07d5f252..af5a5aa4a0 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. @@ -83,9 +83,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) @@ -261,6 +258,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, @@ -270,8 +269,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 @@ -443,25 +440,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) - - # TODO - return boundary_model - 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, From 350c79172ca28e1dc799b0cf2a82407db06bf14a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 3 Jul 2025 17:08:28 +0200 Subject: [PATCH 04/27] rm example --- examples/fluid/play_ground.jl | 233 ---------------------------------- 1 file changed, 233 deletions(-) delete mode 100644 examples/fluid/play_ground.jl diff --git a/examples/fluid/play_ground.jl b/examples/fluid/play_ground.jl deleted file mode 100644 index 5bc88fb7cb..0000000000 --- a/examples/fluid/play_ground.jl +++ /dev/null @@ -1,233 +0,0 @@ -# Flow past a circular cylinder (vortex street), Tafuni et al. (2018). -# Other literature using this validation: -# Vacandio et al. (2013), Marrone et al. (2013), Calhoun (2002), Liu et al. (1998) - -using TrixiParticles -using OrdinaryDiffEq - -# ========================================================================================== -# ==== Resolution -factor_d = 0.08 # Resolution in the paper is `0.01` (5M particles) - -const cylinder_diameter = 0.1 -particle_spacing = factor_d * cylinder_diameter - -# Make sure that the kernel support of fluid particles at a boundary is always fully sampled -boundary_layers = 4 - -# 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 = 8 - -# ========================================================================================== -# ==== Experiment Setup -tspan = (0.0, 1.5) - -# Boundary geometry and initial fluid particle positions -domain_size = (25 * cylinder_diameter, 20 * cylinder_diameter) - -flow_direction = [1.0, 0.0] -reynolds_number = 200 -const prescribed_velocity = 1.0 -const fluid_density = 1000.0 -sound_speed = 10 * prescribed_velocity - -boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, - domain_size[2]) - -pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, - n_layers=boundary_layers, velocity=[prescribed_velocity, 0.0], - 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 = 10 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1) - -cylinder_center = (5 * cylinder_diameter, domain_size[2] / 2) -cylinder = SphereShape(particle_spacing, cylinder_diameter / 2, - cylinder_center, fluid_density, sphere_type=RoundSphere()) - -fluid = setdiff(pipe.fluid, cylinder) - -# ========================================================================================== -# ==== Fluid -wcsph = true - -h_factor = 1.5 -smoothing_length = h_factor * particle_spacing -smoothing_kernel = WendlandC2Kernel{2}() - -fluid_density_calculator = ContinuityDensity() - -kinematic_viscosity = prescribed_velocity * cylinder_diameter / reynolds_number - -state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, - exponent=7) -viscosity = ViscosityAdami(nu=kinematic_viscosity) - -if wcsph - density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) - - fluid_system = WeaklyCompressibleSPHSystem(fluid, fluid_density_calculator, - state_equation, smoothing_kernel, - density_diffusion=density_diffusion, - smoothing_length, viscosity=viscosity, - pressure_acceleration=tensile_instability_control, - buffer_size=n_buffer_particles) -else - state_equation = nothing - fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, smoothing_length, - sound_speed, viscosity=viscosity, - density_calculator=fluid_density_calculator, - buffer_size=n_buffer_particles) -end - -# ========================================================================================== -# ==== Open Boundary -function velocity_function2d(pos, t) - return SVector(prescribed_velocity, 0.0) -end - -open_boundary_model = BoundaryModelTafuni() - -boundary_type_in = InFlow() -plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) -inflow = BoundaryZone(; plane=plane_in, plane_normal=flow_direction, open_boundary_layers, - density=fluid_density, particle_spacing, - boundary_type=boundary_type_in) - -reference_velocity_in = velocity_function2d -# At the inlet, neither pressure nor density are prescribed; instead, -# these values are extrapolated from the fluid domain -reference_pressure_in = nothing -reference_density_in = nothing -open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_in, - reference_pressure=reference_pressure_in, - reference_velocity=reference_velocity_in) - -boundary_type_out = OutFlow() - -plane_out = ([pipe.fluid_size[1], 0.0], [pipe.fluid_size[1], domain_size[2]]) -outflow = BoundaryZone(; plane=plane_out, plane_normal=(-flow_direction), - open_boundary_layers, density=fluid_density, particle_spacing, - boundary_type=boundary_type_out) - -# At the outlet, we allow the flow to exit freely without imposing any boundary conditions -reference_velocity_out = nothing -reference_pressure_out = nothing -reference_density_out = nothing -open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles, - reference_density=reference_density_out, - reference_pressure=reference_pressure_out, - reference_velocity=reference_velocity_out) -# ========================================================================================== -# ==== Boundary -boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - smoothing_kernel, smoothing_length) - -boundary_system_wall = BoundarySPHSystem(pipe.boundary, boundary_model) - -boundary_model_cylinder = BoundaryModelDummyParticles(cylinder.density, cylinder.mass, - AdamiPressureExtrapolation(), - state_equation=state_equation, - viscosity=viscosity, - smoothing_kernel, smoothing_length) - -boundary_system_cylinder = BoundarySPHSystem(cylinder, boundary_model_cylinder) - -# ========================================================================================== -# ==== Postprocessing -circle = SphereShape(particle_spacing, (cylinder_diameter + particle_spacing) / 2, - cylinder_center, fluid_density, n_layers=1, - sphere_type=RoundSphere()) - -# Points for pressure interpolation, located at the wall interface -const data_points = copy(circle.coordinates) -const center = SVector(cylinder_center) - -calculate_lift_force(system, v_ode, u_ode, semi, t) = nothing -function calculate_lift_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) - force = zero(SVector{ndims(system), eltype(system)}) - - values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, - clip_negative_pressure=false) - pressure = Array(values.pressure) - - for i in axes(data_points, 2) - point = TrixiParticles.current_coords(data_points, system, i) - - # F = ∑ -p_i * A_i * n_i - force -= pressure[i] * particle_spacing .* - TrixiParticles.normalize(point - center) - end - - return 2 * force[2] / (fluid_density * prescribed_velocity^2 * cylinder_diameter) -end - -calculate_drag_force(system, v_ode, u_ode, semi, t) = nothing -function calculate_drag_force(system::TrixiParticles.FluidSystem, v_ode, u_ode, semi, t) - force = zero(SVector{ndims(system), eltype(system)}) - - values = interpolate_points(data_points, semi, system, v_ode, u_ode; cut_off_bnd=false, - clip_negative_pressure=false) - pressure = Array(values.pressure) - - for i in axes(data_points, 2) - point = TrixiParticles.current_coords(data_points, system, i) - - # F = ∑ -p_i * A_i * n_i - force -= pressure[i] * particle_spacing .* - TrixiParticles.normalize(point - center) - end - - return 2 * force[1] / (fluid_density * prescribed_velocity^2 * cylinder_diameter) -end - -# ========================================================================================== -# ==== Simulation -min_corner = minimum(pipe.boundary.coordinates .- 5 * particle_spacing, dims=2) -max_corner = maximum(pipe.boundary.coordinates .+ 5 * particle_spacing, dims=2) -cell_list = FullGridCellList(; min_corner, max_corner) - -neighborhood_search = GridNeighborhoodSearch{2}(; cell_list, - update_strategy=ParallelUpdate()) - -semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, - boundary_system_wall, boundary_system_cylinder; - neighborhood_search=neighborhood_search, - parallelization_backend=PolyesterBackend()) - -ode = semidiscretize(semi, tspan) - -info_callback = InfoCallback(interval=100) - -output_directory = "out" - -saving_callback = SolutionSavingCallback(; dt=0.02, prefix="", output_directory) - -pp_callback = PostprocessCallback(; dt=0.02, - f_l=calculate_lift_force, f_d=calculate_drag_force, - output_directory, filename="resulting_force", - write_csv=true, write_file_interval=10) - -extra_callback = nothing - -callbacks = CallbackSet(info_callback, UpdateCallback(), saving_callback, - ParticleShiftingCallback(), # To obtain a near-uniform particle distribution in the wake - pp_callback, extra_callback) - -sol = solve(ode, RDPK3SpFSAL35(), - abstol=1e-6, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) - reltol=1e-4, # 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); From a8b6555eac254c3af6e9001e3af5e6c0fde47aa2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 9 Jul 2025 14:05:41 +0200 Subject: [PATCH 05/27] fix missing nhs update --- src/schemes/boundary/open_boundary/system.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index af5a5aa4a0..cf747dd1fa 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -309,6 +309,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 From ac9108c632637dbdebaba6879640964ad90b9af2 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 9 Jul 2025 22:38:01 +0200 Subject: [PATCH 06/27] implement different mirror methods --- src/TrixiParticles.jl | 1 + .../boundary/open_boundary/mirroring.jl | 282 +++++++++++++++--- 2 files changed, 237 insertions(+), 46 deletions(-) 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/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 869a03bd0a..9de93a244b 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -7,10 +7,33 @@ 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. """ -struct BoundaryModelTafuni end +struct BoundaryModelTafuni{MM} + mirror_method::MM +end + +struct FirstOrderMirroring{ELTYPE} + firstorder_tolerance::ELTYPE + function FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) where {ELTYPE} + return new{typeof(firstorder_tolerance)}(firstorder_tolerance) + end +end + +struct SimpleMirroring{ELTYPE} + firstorder_tolerance::ELTYPE + function SimpleMirroring(; firstorder_tolerance::Real=1e-3) + return new{typeof(firstorder_tolerance)}(firstorder_tolerance) + end +end -function update_boundary_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, u_ode, - semi, t) +struct ZerothOrderMirroring end + +function BoundaryModelTafuni(; + mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) + return BoundaryModelTafuni(mirror_method) +end + +function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, + v, u, v_ode, u_ode, semi, t) @trixi_timeit timer() "extrapolate and correct values" begin fluid_system = corresponding_fluid_system(system, semi) @@ -19,14 +42,16 @@ function update_boundary_quantities!(system, ::BoundaryModelTafuni, v, u, v_ode, 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_open_boundary, v_fluid, + u_open_boundary, u_fluid, semi, t; system.cache...) 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, @@ -34,8 +59,6 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, 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)) @@ -74,13 +97,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, @@ -105,22 +122,189 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, 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[]) + + # pressure + if prescribed_pressure + pressure[particle] = reference_value(reference_pressure, pressure[particle], + particle_coords, t) + else + f_p = L_inv * extrapolated_pressure_correction[] + df_p = f_p[two_to_end] # f_p[2:end] as SVector + + gradient_part = mirror_method isa SimpleMirroring ? 0 : dot(pos_diff, df_p) + + pressure[particle] = f_p[1] + gradient_part + end + + # density + 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] # f_d[2:end] as SVector + + gradient_part = mirror_method isa SimpleMirroring ? 0 : dot(pos_diff, df_d) + + density[particle] = f_d[1] + gradient_part + end + + # velocity + 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] + end + else + @inbounds for dim in eachindex(pos_diff) + f_v = L_inv * extrapolated_velocity_correction[][dim, :] + df_v = f_v[two_to_end] # f_v[2:end] as SVector + + gradient_part = mirror_method isa SimpleMirroring ? 0 : + dot(pos_diff, df_v) + + v_open_boundary[dim, particle] = f_v[1] + gradient_part + end + + # 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 + + 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 + pressure[particle] = reference_value(reference_pressure, pressure[particle], + particle_coords, t) + else + pressure[particle] = first(extrapolated_pressure_correction[]) / + shepard_coefficient + end + + # density + if prescribed_density + density[particle] = reference_value(reference_density, density[particle], + particle_coords, t) + else + density[particle] = first(extrapolated_density_correction[]) / + shepard_coefficient + end + + # velocity + 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] + end + else + velocity_interpolated = extrapolated_velocity_correction[][:, 1] / + shepard_coefficient + + @inbounds for dim in eachindex(velocity_interpolated) + v_open_boundary[dim, particle] = velocity_interpolated[dim] + end + + # 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) + end + + return system +end + +function extrapolate_values!(system, ::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, reference_density, + reference_velocity, reference_pressure) = system + + fluid_system = corresponding_fluid_system(system, semi) + + # 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) + vel_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_velocity + interpolated_velocity[] += vel_b * volume_b * W_ab + end + + if !prescribed_density + interpolated_density[] += rho_b * volume_b * W_ab + end + end + + if shepard_coefficient[] > sqrt(eps()) + interpolated_density[] /= shepard_coefficient[] + interpolated_pressure[] /= shepard_coefficient[] + interpolated_velocity[] /= shepard_coefficient[] + else + interpolated_density[] = current_density(v_open_boundary, system, particle) + interpolated_pressure[] = current_pressure(v_open_boundary, system, particle) + interpolated_velocity[] = current_velocity(v_open_boundary, system, particle) end pos_diff = particle_coords - ghost_node_position - # 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) @@ -129,31 +313,29 @@ function extrapolate_values!(system, v_open_boundary, v_fluid, u_open_boundary, 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) + v_open_boundary[dim, particle] = interpolated_velocity[][dim] end - end - if prescribed_pressure - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) - 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) + # 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 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] + density[particle] = interpolated_density[] + end - density[particle] = f_d[1] + dot(pos_diff, df_d) + if prescribed_pressure + pressure[particle] = reference_value(reference_pressure, pressure[particle], + particle_coords, t) + else + pressure[particle] = interpolated_pressure[] end end @@ -252,12 +434,20 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se 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 + vel = current_velocity(v, system, particle) + vel_ = dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal + + @inbounds for dim in eachindex(vel) + v[dim, particle] = vel_[dim] + end + + return v end From d69f03d158fc0d2ac3bf2c20fbbfd173822c309a Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 10:48:32 +0200 Subject: [PATCH 07/27] add docs --- docs/src/refs.bib | 26 +++++++++++++++++ .../method_of_characteristics.jl | 23 +++++++++------ .../boundary/open_boundary/mirroring.jl | 28 ++++++++++++++++++- 3 files changed, 67 insertions(+), 10 deletions(-) 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/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 757a9328af..21ae82f426 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -1,5 +1,5 @@ @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 @@ -8,15 +8,19 @@ 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. +- `extrapolate_reference_values=nothing`: If one of the follwoing mirror methods is choosen, + 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. """ -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 @@ -31,13 +35,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 diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 9de93a244b..f820d69879 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -1,16 +1,28 @@ @doc raw""" - BoundaryModelTafuni() + BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) 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{MM} mirror_method::MM end +""" + FirstOrderMirroring + +Fluid properties are interpolated onto ghost nodes using the method propsed by [Liu and Liu (2006)](@cite Liu2006), +to retrieve first oder kernel and particle consistency. +""" struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE function FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) where {ELTYPE} @@ -18,6 +30,12 @@ struct FirstOrderMirroring{ELTYPE} end end +""" + SimpleMirroring + +This method is similar to [`FirstOrderMirroring`](@ref), but does not use +the corrected gradient as proposed by [Negi et al. (2022)](@cite Negi2022). +""" struct SimpleMirroring{ELTYPE} firstorder_tolerance::ELTYPE function SimpleMirroring(; firstorder_tolerance::Real=1e-3) @@ -25,6 +43,14 @@ struct SimpleMirroring{ELTYPE} end end +""" + ZerothOrderMirroring + +Fluid properties are interpolated onto ghost nodes using Shepard interpolation. +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 function BoundaryModelTafuni(; From 8863884afdd8afdccb8d2f2e90330d05cfd10b54 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 10:56:17 +0200 Subject: [PATCH 08/27] fix tests --- examples/fluid/pipe_flow_2d.jl | 2 +- test/schemes/boundary/open_boundary/mirroring.jl | 9 ++++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index d21a684979..430161c8b7 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -49,7 +49,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/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index cc4c970cb9..72613c4100 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, @@ -240,7 +242,8 @@ 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) From 8beb4d994ec6ba970d03276ca2cece93558cd9d3 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 11:05:49 +0200 Subject: [PATCH 09/27] fix typos --- .../boundary/open_boundary/method_of_characteristics.jl | 4 ++-- src/schemes/boundary/open_boundary/mirroring.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 21ae82f426..a2c9e347e1 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -8,8 +8,8 @@ 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=nothing`: If one of the follwoing mirror methods is choosen, - the reference values are extrapolated from the fluid domain to the boundary particles. +- `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) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index f820d69879..f3f8adf40f 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -20,8 +20,8 @@ end """ FirstOrderMirroring -Fluid properties are interpolated onto ghost nodes using the method propsed by [Liu and Liu (2006)](@cite Liu2006), -to retrieve first oder kernel and particle consistency. +Fluid properties are interpolated onto ghost nodes using the method proposed by [Liu and Liu (2006)](@cite Liu2006), +to retrieve first order kernel and particle consistency. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE From f06031842a763b41d6babddafa86ef145f4a24b4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 11:16:36 +0200 Subject: [PATCH 10/27] add docs --- src/schemes/boundary/open_boundary/mirroring.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index f3f8adf40f..95dbfa1e17 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -18,10 +18,14 @@ struct BoundaryModelTafuni{MM} end """ - FirstOrderMirroring + FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) Fluid properties are interpolated onto ghost nodes using the method proposed by [Liu and Liu (2006)](@cite Liu2006), to retrieve first order kernel and particle consistency. + +# Keywords +- `firstorder_tolerance`: If the determinant of the correction matrix is smaller than this value, + the method falls back to [`ZerothOrderMirroring`](@ref). Default is `1e-3`. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE @@ -31,10 +35,14 @@ struct FirstOrderMirroring{ELTYPE} end """ - SimpleMirroring + SimpleMirroring(; firstorder_tolerance::ELTYPE=1e-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`: If the determinant of the correction matrix is smaller than this value, + the method falls back to [`ZerothOrderMirroring`](@ref). Default is `1e-3`. """ struct SimpleMirroring{ELTYPE} firstorder_tolerance::ELTYPE From bf610b8cc54055a928233648f7cf53f17bc5d8a8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 12:10:51 +0200 Subject: [PATCH 11/27] add tests --- .../boundary/open_boundary/mirroring.jl | 274 ++++++++++++++++++ 1 file changed, 274 insertions(+) diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 72613c4100..7c93df5617 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -253,4 +253,278 @@ @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)) + 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) + + 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) + + 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) + 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 fluid system 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]) + + 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 From 8ca3825e6369cfb6be1aec3b5c868f890a4b8202 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 10 Jul 2025 12:14:18 +0200 Subject: [PATCH 12/27] apply formatter --- test/schemes/boundary/open_boundary/mirroring.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 7c93df5617..76e3829d45 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -522,8 +522,8 @@ ] @testset verbose=true "$method" for (i, method) in enumerate(("simple mirroring", - "first order mirroring", - "zeroth order mirroring")) + "first order mirroring", + "zeroth order mirroring")) @test isapprox(pressures[i], pressures_expected[i]) end end From 15c737c4423889aa88e98e83862452a5b896ffc8 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 15 Jul 2025 08:43:21 +0200 Subject: [PATCH 13/27] add `average_velocity!` for `BoundaryModelLastiwka` --- .../method_of_characteristics.jl | 29 +++++++++++++++++++ .../boundary/open_boundary/mirroring.jl | 8 +++-- 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index a2c9e347e1..37384e16f2 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -75,6 +75,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 non-uniform buffer particles 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 @@ -222,3 +229,25 @@ end return characteristics end + +function average_velocity!(v, u, system, ::BoundaryModelLastiwka, boundary_zone, semi) + return v +end + +function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone{InFlow}, + semi) + avg_velocity = sum(each_moving_particle(system)) do particle + return current_velocity(v, system, particle) + end + + avg_velocity /= system.buffer.active_particle_count[] + + @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] + end + end + + return v +end diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 95dbfa1e17..ac37a339ca 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -265,7 +265,8 @@ function extrapolate_values!(system, # 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) + average_velocity!(v_open_boundary, u_open_boundary, system, boundary_model, + boundary_zone, semi) end return system @@ -438,9 +439,10 @@ function mirror_position(particle_coords, boundary_zone) return particle_coords - 2 * dist * boundary_zone.plane_normal end -average_velocity!(v, u, system, boundary_zone, semi) = v +average_velocity!(v, u, system, ::BoundaryModelTafuni, boundary_zone, semi) = v -function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, semi) +function average_velocity!(v, u, system, ::BoundaryModelTafuni, + boundary_zone::BoundaryZone{InFlow}, semi) (; plane_normal, zone_origin, initial_condition) = boundary_zone # We only use the extrapolated velocity in the vicinity of the transition region. From f285ccd29f0670225fbc31191606ff81437ea73e Mon Sep 17 00:00:00 2001 From: LasNikas Date: Tue, 15 Jul 2025 09:00:34 +0200 Subject: [PATCH 14/27] fix --- src/schemes/boundary/open_boundary/mirroring.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index ac37a339ca..95dbfa1e17 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -265,8 +265,7 @@ function extrapolate_values!(system, # 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_model, - boundary_zone, semi) + average_velocity!(v_open_boundary, u_open_boundary, system, boundary_zone, semi) end return system @@ -439,10 +438,9 @@ function mirror_position(particle_coords, boundary_zone) return particle_coords - 2 * dist * boundary_zone.plane_normal end -average_velocity!(v, u, system, ::BoundaryModelTafuni, boundary_zone, semi) = v +average_velocity!(v, u, system, boundary_zone, semi) = v -function average_velocity!(v, u, system, ::BoundaryModelTafuni, - boundary_zone::BoundaryZone{InFlow}, semi) +function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, semi) (; plane_normal, zone_origin, initial_condition) = boundary_zone # We only use the extrapolated velocity in the vicinity of the transition region. From 5eb6442235d876233ed802cab949085e9d0b67b9 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Jul 2025 12:19:19 +0200 Subject: [PATCH 15/27] add interpolation functions --- .../boundary/open_boundary/mirroring.jl | 331 ++++++++++-------- 1 file changed, 187 insertions(+), 144 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 95dbfa1e17..6a35c748ed 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -1,22 +1,3 @@ -@doc raw""" - BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) - -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{MM} - mirror_method::MM -end - """ FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) @@ -61,6 +42,25 @@ The interpolated values at the ghost nodes are then assigned to the correspondin """ struct ZerothOrderMirroring end +@doc raw""" + BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) + +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{MM} + mirror_method::MM +end + function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) return BoundaryModelTafuni(mirror_method) @@ -68,16 +68,49 @@ end function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) + (; reference_pressure, reference_density, reference_velocity, cache) = system + (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache + @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, boundary_model.mirror_method, 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_pressure + @threaded semi for particle in each_moving_particle(system) + particle_coords = current_coords(u_open_boundary, 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_open_boundary, 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_open_boundary, system, particle) + 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[dim, particle] = v_ref[dim] + end + end end end @@ -88,8 +121,7 @@ function extrapolate_values!(system, 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) @@ -114,13 +146,13 @@ function extrapolate_values!(system, 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)})) @@ -143,16 +175,16 @@ function extrapolate_values!(system, correction_matrix[] += L - if !prescribed_pressure - 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 @@ -163,48 +195,24 @@ function extrapolate_values!(system, L_inv = inv(correction_matrix[]) # pressure - if prescribed_pressure - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) - else - f_p = L_inv * extrapolated_pressure_correction[] - df_p = f_p[two_to_end] # f_p[2:end] as SVector - - gradient_part = mirror_method isa SimpleMirroring ? 0 : dot(pos_diff, df_p) - - pressure[particle] = f_p[1] + gradient_part + if !(prescribed_pressure) + first_order_scalar_interpolation!(pressure, particle, L_inv, + interpolated_pressure_correction, + two_to_end, pos_diff, mirror_method) end # density - 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] # f_d[2:end] as SVector - - gradient_part = mirror_method isa SimpleMirroring ? 0 : dot(pos_diff, df_d) - - density[particle] = f_d[1] + gradient_part + if !(prescribed_density) + first_order_scalar_interpolation!(density, particle, L_inv, + interpolated_density_correction, + two_to_end, pos_diff, mirror_method) end # velocity - 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] - end - else - @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * extrapolated_velocity_correction[][dim, :] - df_v = f_v[two_to_end] # f_v[2:end] as SVector - - gradient_part = mirror_method isa SimpleMirroring ? 0 : - dot(pos_diff, df_v) - - v_open_boundary[dim, particle] = f_v[1] + gradient_part - end + if !(prescribed_velocity) + first_order_velocity_interpolation!(v_open_boundary, system, 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.: @@ -219,37 +227,25 @@ function extrapolate_values!(system, shepard_coefficient = correction_matrix[][1, 1] # pressure - if prescribed_pressure - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) - else - pressure[particle] = first(extrapolated_pressure_correction[]) / - shepard_coefficient + if !(prescribed_pressure) + interpolated_pressure = first(interpolated_pressure_correction[]) + zeroth_order_scalar_interpolation!(pressure, particle, shepard_coefficient, + interpolated_pressure) end # density - if prescribed_density - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) - else - density[particle] = first(extrapolated_density_correction[]) / - shepard_coefficient + if !(prescribed_density) + interpolated_density = first(interpolated_density_correction[]) + zeroth_order_scalar_interpolation!(density, particle, shepard_coefficient, + interpolated_density) end # velocity - 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] - end - else - velocity_interpolated = extrapolated_velocity_correction[][:, 1] / - shepard_coefficient - - @inbounds for dim in eachindex(velocity_interpolated) - v_open_boundary[dim, particle] = velocity_interpolated[dim] - end + if !(prescribed_velocity) + interpolated_velocity = interpolated_velocity_correction[][:, 1] + zeroth_order_velocity_interpolation!(v_open_boundary, system, 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.: @@ -271,12 +267,11 @@ function extrapolate_values!(system, return system end -function extrapolate_values!(system, ::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, reference_density, - reference_velocity, reference_pressure) = system +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) @@ -314,62 +309,46 @@ function extrapolate_values!(system, ::ZerothOrderMirroring, shepard_coefficient[] += volume_b * W_ab - if !prescribed_pressure + if !(prescribed_pressure) interpolated_pressure[] += pressure_b * volume_b * W_ab end - if !prescribed_velocity - interpolated_velocity[] += vel_b * volume_b * W_ab - end - - if !prescribed_density + if !(prescribed_density) interpolated_density[] += rho_b * volume_b * W_ab end - end - if shepard_coefficient[] > sqrt(eps()) - interpolated_density[] /= shepard_coefficient[] - interpolated_pressure[] /= shepard_coefficient[] - interpolated_velocity[] /= shepard_coefficient[] - else - interpolated_density[] = current_density(v_open_boundary, system, particle) - interpolated_pressure[] = current_pressure(v_open_boundary, system, particle) - interpolated_velocity[] = current_velocity(v_open_boundary, system, particle) + if !(prescribed_velocity) + interpolated_velocity[] += vel_b * volume_b * W_ab + end end - pos_diff = particle_coords - ghost_node_position + if shepard_coefficient[] > eps() + pos_diff = particle_coords - ghost_node_position - 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] - end - else - @inbounds for dim in eachindex(pos_diff) - v_open_boundary[dim, particle] = interpolated_velocity[][dim] + if !(prescribed_pressure) + zeroth_order_scalar_interpolation!(pressure, particle, + shepard_coefficient[], + interpolated_pressure[]) end - # 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 + if !(prescribed_density) + zeroth_order_scalar_interpolation!(density, particle, + shepard_coefficient[], + interpolated_density[]) + end - if prescribed_density - density[particle] = reference_value(reference_density, density[particle], - particle_coords, t) - else - density[particle] = interpolated_density[] - end + if !(prescribed_velocity) + zeroth_order_velocity_interpolation!(v_open_boundary, system, particle, + shepard_coefficient[], + interpolated_velocity[]) - if prescribed_pressure - pressure[particle] = reference_value(reference_pressure, pressure[particle], - particle_coords, t) - else - pressure[particle] = interpolated_pressure[] + # 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 @@ -383,6 +362,70 @@ function extrapolate_values!(system, ::ZerothOrderMirroring, return system end +function zeroth_order_scalar_interpolation!(values, particle, shepard_coefficient, + extrapolated_value) + values[particle] = extrapolated_value / shepard_coefficient + + return values +end + +function zeroth_order_velocity_interpolation!(v, system, particle, shepard_coefficient, + extrapolated_velocity) + velocity_interpolated = extrapolated_velocity / shepard_coefficient + + @inbounds for dim in eachindex(velocity_interpolated) + v[dim, particle] = velocity_interpolated[dim] + end + + return v +end + +function first_order_scalar_interpolation!(values, particle, L_inv, + extrapolated_values_correction, + two_to_end, pos_diff, ::FirstOrderMirroring) + f_s = L_inv * extrapolated_values_correction[] + df_p = f_s[two_to_end] # f_s[2:end] as SVector + + values[particle] = f_s[1] + dot(pos_diff, df_p) + + return values +end + +function first_order_scalar_interpolation!(values, particle, L_inv, + extrapolated_values_correction, + two_to_end, pos_diff, ::SimpleMirroring) + f_s = L_inv * extrapolated_values_correction[] + + values[particle] = f_s[1] + + return values +end + +function first_order_velocity_interpolation!(v, system, particle, L_inv, + interpolated_velocity_correction, + two_to_end, pos_diff, ::FirstOrderMirroring) + @inbounds for dim in eachindex(pos_diff) + f_v = L_inv * interpolated_velocity_correction[][dim, :] + df_v = f_v[two_to_end] # f_v[2:end] as SVector + + v[dim, particle] = f_v[1] + dot(pos_diff, df_v) + end + + return v +end + +function first_order_velocity_interpolation!(v, system, particle, L_inv, + interpolated_velocity_correction, + two_to_end, pos_diff, ::SimpleMirroring) + @inbounds for dim in eachindex(pos_diff) + f_v = L_inv * interpolated_velocity_correction[][dim, :] + + 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) # `pos_diff` corresponds to `x_{kl} = x_k - x_l` in the paper (Tafuni et al., 2018), # where `x_k` is the position of the ghost node and `x_l` is the position of the neighbor particle From 888ef53df8b3b7eeff0673287c0b0ade43202068 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Jul 2025 12:51:15 +0200 Subject: [PATCH 16/27] restructure again --- .../method_of_characteristics.jl | 4 +- .../boundary/open_boundary/mirroring.jl | 38 +++++++++---------- .../boundary/open_boundary/mirroring.jl | 4 ++ 3 files changed, 25 insertions(+), 21 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 37384e16f2..04f34799b0 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -3,7 +3,7 @@ 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). @@ -77,7 +77,7 @@ 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 non-uniform buffer particles distribution, + # 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 diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 6a35c748ed..3252159e8a 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -68,7 +68,8 @@ end function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni, v, u, v_ode, u_ode, semi, t) - (; reference_pressure, reference_density, reference_velocity, cache) = system + (; reference_pressure, reference_density, reference_velocity, boundary_zone, + cache) = system (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache @trixi_timeit timer() "extrapolate and correct values" begin @@ -82,9 +83,16 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni 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_open_boundary, system, particle) + particle_coords = current_coords(u, system, particle) pressure[particle] = reference_value(reference_pressure, pressure[particle], particle_coords, t) @@ -93,7 +101,7 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni if prescribed_density @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u_open_boundary, system, particle) + particle_coords = current_coords(u, system, particle) density[particle] = reference_value(reference_density, density[particle], particle_coords, t) @@ -102,8 +110,8 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni if prescribed_velocity @threaded semi for particle in each_moving_particle(system) - particle_coords = current_coords(u_open_boundary, system, particle) - v_particle = current_velocity(v_open_boundary, system, particle) + 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) @@ -228,6 +236,8 @@ function extrapolate_values!(system, # 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_interpolation!(pressure, particle, shepard_coefficient, interpolated_pressure) @@ -235,6 +245,8 @@ function extrapolate_values!(system, # 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_interpolation!(density, particle, shepard_coefficient, interpolated_density) @@ -242,6 +254,8 @@ function extrapolate_values!(system, # 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, system, particle, shepard_coefficient, @@ -257,13 +271,6 @@ function extrapolate_values!(system, 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) - end - return system end @@ -352,13 +359,6 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, 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) - end - return system end diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 76e3829d45..4efe2fe2e7 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -238,6 +238,7 @@ 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) @@ -247,6 +248,9 @@ 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] From 4fd1b32458dfeb4f10fdfb7811721b78759c31b4 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Jul 2025 16:02:15 +0200 Subject: [PATCH 17/27] fix gpu bugs --- .../method_of_characteristics.jl | 4 +-- .../boundary/open_boundary/mirroring.jl | 33 ++++++++++--------- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 04f34799b0..20aeb305b1 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -237,11 +237,9 @@ end function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone{InFlow}, semi) avg_velocity = sum(each_moving_particle(system)) do particle - return current_velocity(v, system, particle) + return current_velocity(v, system, particle) / system.buffer.active_particle_count[] end - avg_velocity /= system.buffer.active_particle_count[] - @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) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 3252159e8a..06094f9fde 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -10,7 +10,7 @@ to retrieve first order kernel and particle consistency. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) where {ELTYPE} + function FirstOrderMirroring(; firstorder_tolerance::Real=1 / 1000.0f0) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -27,7 +27,7 @@ the corrected gradient as proposed by [Negi et al. (2022)](@cite Negi2022). """ struct SimpleMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function SimpleMirroring(; firstorder_tolerance::Real=1e-3) + function SimpleMirroring(; firstorder_tolerance::Real=1 / 1000.0f0) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -62,14 +62,16 @@ struct BoundaryModelTafuni{MM} end function BoundaryModelTafuni(; - mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) + mirror_method=FirstOrderMirroring(; + firstorder_tolerance=1 / + 1000.0f0)) 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, - cache) = system + pressure, density, cache) = system (; prescribed_pressure, prescribed_density, prescribed_velocity) = cache @trixi_timeit timer() "extrapolate and correct values" begin @@ -205,21 +207,22 @@ function extrapolate_values!(system, # pressure if !(prescribed_pressure) first_order_scalar_interpolation!(pressure, particle, L_inv, - interpolated_pressure_correction, + interpolated_pressure_correction[], two_to_end, pos_diff, mirror_method) end # density if !(prescribed_density) first_order_scalar_interpolation!(density, particle, L_inv, - interpolated_density_correction, + interpolated_density_correction[], two_to_end, pos_diff, mirror_method) end # velocity if !(prescribed_velocity) first_order_velocity_interpolation!(v_open_boundary, system, particle, - L_inv, interpolated_velocity_correction, + 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). @@ -381,9 +384,9 @@ function zeroth_order_velocity_interpolation!(v, system, particle, shepard_coeff end function first_order_scalar_interpolation!(values, particle, L_inv, - extrapolated_values_correction, + interpolated_values_correction, two_to_end, pos_diff, ::FirstOrderMirroring) - f_s = L_inv * extrapolated_values_correction[] + f_s = L_inv * interpolated_values_correction df_p = f_s[two_to_end] # f_s[2:end] as SVector values[particle] = f_s[1] + dot(pos_diff, df_p) @@ -392,9 +395,9 @@ function first_order_scalar_interpolation!(values, particle, L_inv, end function first_order_scalar_interpolation!(values, particle, L_inv, - extrapolated_values_correction, + interpolated_values_correction, two_to_end, pos_diff, ::SimpleMirroring) - f_s = L_inv * extrapolated_values_correction[] + f_s = L_inv * interpolated_values_correction values[particle] = f_s[1] @@ -405,7 +408,7 @@ function first_order_velocity_interpolation!(v, system, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::FirstOrderMirroring) @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * interpolated_velocity_correction[][dim, :] + f_v = L_inv * interpolated_velocity_correction[dim, :] df_v = f_v[two_to_end] # f_v[2:end] as SVector v[dim, particle] = f_v[1] + dot(pos_diff, df_v) @@ -418,7 +421,7 @@ function first_order_velocity_interpolation!(v, system, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::SimpleMirroring) @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * interpolated_velocity_correction[][dim, :] + f_v = L_inv * interpolated_velocity_correction[dim, :] v[dim, particle] = f_v[1] end @@ -496,11 +499,9 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se active_coordinates(u, system))) 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) From fb9f33353293c5ec7435e721255619051b907604 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Jul 2025 16:14:23 +0200 Subject: [PATCH 18/27] remove unnecessary argument --- .../boundary/open_boundary/mirroring.jl | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 06094f9fde..070d6e3eb6 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -10,7 +10,7 @@ to retrieve first order kernel and particle consistency. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function FirstOrderMirroring(; firstorder_tolerance::Real=1 / 1000.0f0) + function FirstOrderMirroring(; firstorder_tolerance::Real=1/1000) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -27,7 +27,7 @@ the corrected gradient as proposed by [Negi et al. (2022)](@cite Negi2022). """ struct SimpleMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function SimpleMirroring(; firstorder_tolerance::Real=1 / 1000.0f0) + function SimpleMirroring(; firstorder_tolerance::Real=1/1000) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -63,8 +63,7 @@ end function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; - firstorder_tolerance=1 / - 1000.0f0)) + firstorder_tolerance=1/1000)) return BoundaryModelTafuni(mirror_method) end @@ -220,8 +219,7 @@ function extrapolate_values!(system, # velocity if !(prescribed_velocity) - first_order_velocity_interpolation!(v_open_boundary, system, particle, - L_inv, + first_order_velocity_interpolation!(v_open_boundary, particle, L_inv, interpolated_velocity_correction[], two_to_end, pos_diff, mirror_method) @@ -260,7 +258,7 @@ function extrapolate_values!(system, # 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, system, particle, + zeroth_order_velocity_interpolation!(v_open_boundary, particle, shepard_coefficient, interpolated_velocity) @@ -348,7 +346,7 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, end if !(prescribed_velocity) - zeroth_order_velocity_interpolation!(v_open_boundary, system, particle, + zeroth_order_velocity_interpolation!(v_open_boundary, particle, shepard_coefficient[], interpolated_velocity[]) @@ -372,7 +370,7 @@ function zeroth_order_scalar_interpolation!(values, particle, shepard_coefficien return values end -function zeroth_order_velocity_interpolation!(v, system, particle, shepard_coefficient, +function zeroth_order_velocity_interpolation!(v, particle, shepard_coefficient, extrapolated_velocity) velocity_interpolated = extrapolated_velocity / shepard_coefficient @@ -404,7 +402,7 @@ function first_order_scalar_interpolation!(values, particle, L_inv, return values end -function first_order_velocity_interpolation!(v, system, particle, L_inv, +function first_order_velocity_interpolation!(v, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::FirstOrderMirroring) @inbounds for dim in eachindex(pos_diff) @@ -417,7 +415,7 @@ function first_order_velocity_interpolation!(v, system, particle, L_inv, return v end -function first_order_velocity_interpolation!(v, system, particle, L_inv, +function first_order_velocity_interpolation!(v, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::SimpleMirroring) @inbounds for dim in eachindex(pos_diff) From c5fa956eb83d44d11da938ba35f77239524f7f31 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Wed, 23 Jul 2025 16:58:43 +0200 Subject: [PATCH 19/27] fix gpu again --- src/schemes/boundary/open_boundary/mirroring.jl | 6 +++--- test/examples/gpu.jl | 8 ++++++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 070d6e3eb6..34bc47dda1 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -10,7 +10,7 @@ to retrieve first order kernel and particle consistency. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function FirstOrderMirroring(; firstorder_tolerance::Real=1/1000) + function FirstOrderMirroring(; firstorder_tolerance::Real=1e-3) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -27,7 +27,7 @@ the corrected gradient as proposed by [Negi et al. (2022)](@cite Negi2022). """ struct SimpleMirroring{ELTYPE} firstorder_tolerance::ELTYPE - function SimpleMirroring(; firstorder_tolerance::Real=1/1000) + function SimpleMirroring(; firstorder_tolerance::Real=1e-3) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end @@ -63,7 +63,7 @@ end function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; - firstorder_tolerance=1/1000)) + firstorder_tolerance=1e-3)) return BoundaryModelTafuni(mirror_method) end diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index a9a2e791dd..df0cef63e4 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -319,7 +319,9 @@ end joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(), + open_boundary_model=BoundaryModelTafuni(; + mirror_method=FirstOrderMirroring(; + firstorder_tolerance=0.001f0)), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, @@ -340,7 +342,9 @@ end "pipe_flow_2d.jl"), wcsph=true, sound_speed=20.0f0, pressure=0.0f0, - open_boundary_model=BoundaryModelTafuni(), + open_boundary_model=BoundaryModelTafuni(; + mirror_method=FirstOrderMirroring(; + firstorder_tolerance=0.001f0)), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, From 0925255de2c6494e8c05274ab79422062af7296f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 24 Jul 2025 08:39:57 +0200 Subject: [PATCH 20/27] fix gpu --- test/examples/gpu.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index df0cef63e4..f6490f4a06 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -343,8 +343,7 @@ end wcsph=true, sound_speed=20.0f0, pressure=0.0f0, open_boundary_model=BoundaryModelTafuni(; - mirror_method=FirstOrderMirroring(; - firstorder_tolerance=0.001f0)), + mirror_method=ZerothOrderMirroring()), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, From 3a001f308b2498facb69a45cf6e4a43863c86502 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 24 Jul 2025 16:59:52 +0200 Subject: [PATCH 21/27] implement suggestions --- docs/src/refs.bib | 9 +- .../method_of_characteristics.jl | 6 +- .../boundary/open_boundary/mirroring.jl | 143 +++++++++--------- test/examples/gpu.jl | 4 +- 4 files changed, 85 insertions(+), 77 deletions(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 2f2a636dd1..b7711bc313 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -761,7 +761,14 @@ @book{Poling2001 publisher = {McGraw-Hill}, address = {New York} } - +@article{Shepard1968, + author = {Shepard, Donald}, + booktitle = {Proceedings of the 1968 23rd ACM national conference on -}, + title = {A two-dimensional interpolation function for irregularly-spaced data}, + year = {1968}, + publisher = {ACM Press}, + doi = {10.1145/800186.810616}, +} @article{Smagorinsky1963, author = {Smagorinsky, Joseph}, title = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment}, diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 20aeb305b1..061db7cd42 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -19,6 +19,7 @@ For more information about the method see [description below](@ref method_of_cha """ struct BoundaryModelLastiwka{T} extrapolate_reference_values::T + function BoundaryModelLastiwka(; extrapolate_reference_values=nothing) return new{typeof(extrapolate_reference_values)}(extrapolate_reference_values) end @@ -231,6 +232,7 @@ end end function average_velocity!(v, u, system, ::BoundaryModelLastiwka, boundary_zone, semi) + # Only apply averaging at the inflow return v end @@ -242,8 +244,8 @@ function average_velocity!(v, u, system, ::BoundaryModelLastiwka, ::BoundaryZone @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 diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 34bc47dda1..9ffa11ba91 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -1,41 +1,42 @@ """ - FirstOrderMirroring(; firstorder_tolerance::ELTYPE=1e-3) + FirstOrderMirroring(; firstorder_tolerance::=1f-3) -Fluid properties are interpolated onto ghost nodes using the method proposed by [Liu and Liu (2006)](@cite Liu2006), -to retrieve first order kernel and particle consistency. +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`: If the determinant of the correction matrix is smaller than this value, - the method falls back to [`ZerothOrderMirroring`](@ref). Default is `1e-3`. +- `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=1e-3) + function FirstOrderMirroring(; firstorder_tolerance::Real=1.0f-3) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end """ - SimpleMirroring(; firstorder_tolerance::ELTYPE=1e-3)) + 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`: If the determinant of the correction matrix is smaller than this value, - the method falls back to [`ZerothOrderMirroring`](@ref). Default is `1e-3`. +- `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=1e-3) + + function SimpleMirroring(; firstorder_tolerance::Real=1.0f-3) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end end """ - ZerothOrderMirroring + ZerothOrderMirroring() -Fluid properties are interpolated onto ghost nodes using Shepard interpolation. +Fluid properties are interpolated onto ghost nodes using Shepard interpolation [Shepard1968](@cite). 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. @@ -43,7 +44,7 @@ The interpolated values at the ghost nodes are then assigned to the correspondin struct ZerothOrderMirroring end @doc raw""" - BoundaryModelTafuni(; mirror_method=FirstOrderMirroring(; firstorder_tolerance=1e-3)) + 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 @@ -61,9 +62,7 @@ struct BoundaryModelTafuni{MM} mirror_method::MM end -function BoundaryModelTafuni(; - mirror_method=FirstOrderMirroring(; - firstorder_tolerance=1e-3)) +function BoundaryModelTafuni(; mirror_method=FirstOrderMirroring()) return BoundaryModelTafuni(mirror_method) end @@ -84,7 +83,7 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni prescribed_density, prescribed_velocity) end - if !(prescribed_velocity) && boundary_zone.average_inflow_velocity + 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. @@ -116,8 +115,8 @@ function update_boundary_quantities!(system, boundary_model::BoundaryModelTafuni v_ref = reference_value(reference_velocity, v_particle, particle_coords, t) - @inbounds for dim in eachindex(v_ref) - v[dim, particle] = v_ref[dim] + for dim in eachindex(v_ref) + @inbounds v[dim, particle] = v_ref[dim] end end end @@ -184,15 +183,15 @@ function extrapolate_values!(system, correction_matrix[] += L - if !(prescribed_pressure) + if !prescribed_pressure interpolated_pressure_correction[] += pressure_b * R end - if !(prescribed_density) + if !prescribed_density interpolated_density_correction[] += rho_b * R end - if !(prescribed_velocity) + if !prescribed_velocity interpolated_velocity_correction[] += v_b * R' end end @@ -204,22 +203,22 @@ function extrapolate_values!(system, L_inv = inv(correction_matrix[]) # pressure - if !(prescribed_pressure) - first_order_scalar_interpolation!(pressure, particle, L_inv, + if !prescribed_pressure + first_order_scalar_extrapolation!(pressure, particle, L_inv, interpolated_pressure_correction[], two_to_end, pos_diff, mirror_method) end # density - if !(prescribed_density) - first_order_scalar_interpolation!(density, particle, L_inv, + 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_interpolation!(v_open_boundary, particle, L_inv, + if !prescribed_velocity + first_order_velocity_extrapolation!(v_open_boundary, particle, L_inv, interpolated_velocity_correction[], two_to_end, pos_diff, mirror_method) @@ -231,32 +230,34 @@ function extrapolate_values!(system, boundary_zone) end + # 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) + # 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 + # components that are not required for zeroth-order interpolation. interpolated_pressure = first(interpolated_pressure_correction[]) - zeroth_order_scalar_interpolation!(pressure, particle, shepard_coefficient, + zeroth_order_scalar_extrapolation!(pressure, particle, shepard_coefficient, interpolated_pressure) end - # density - if !(prescribed_density) + # 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 + # components that are not required for zeroth-order interpolation. interpolated_density = first(interpolated_density_correction[]) - zeroth_order_scalar_interpolation!(density, particle, shepard_coefficient, + zeroth_order_scalar_extrapolation!(density, particle, shepard_coefficient, interpolated_density) end - # velocity - if !(prescribed_velocity) + # 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 + # 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, @@ -311,41 +312,41 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, rho_b = current_density(v_fluid, fluid_system, neighbor) volume_b = m_b / rho_b pressure_b = current_pressure(v_fluid, fluid_system, neighbor) - vel_b = current_velocity(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) + if !prescribed_pressure interpolated_pressure[] += pressure_b * volume_b * W_ab end - if !(prescribed_density) + if !prescribed_density interpolated_density[] += rho_b * volume_b * W_ab end - if !(prescribed_velocity) - interpolated_velocity[] += vel_b * volume_b * W_ab + if !prescribed_velocity + interpolated_velocity[] += v_b * volume_b * W_ab end end if shepard_coefficient[] > eps() pos_diff = particle_coords - ghost_node_position - if !(prescribed_pressure) - zeroth_order_scalar_interpolation!(pressure, particle, + if !prescribed_pressure + zeroth_order_scalar_extrapolation!(pressure, particle, shepard_coefficient[], interpolated_pressure[]) end - if !(prescribed_density) - zeroth_order_scalar_interpolation!(density, particle, + if !prescribed_density + zeroth_order_scalar_extrapolation!(density, particle, shepard_coefficient[], interpolated_density[]) end - if !(prescribed_velocity) + if !prescribed_velocity zeroth_order_velocity_interpolation!(v_open_boundary, particle, shepard_coefficient[], interpolated_velocity[]) @@ -363,7 +364,7 @@ function extrapolate_values!(system, mirror_method::ZerothOrderMirroring, return system end -function zeroth_order_scalar_interpolation!(values, particle, shepard_coefficient, +function zeroth_order_scalar_extrapolation!(values, particle, shepard_coefficient, extrapolated_value) values[particle] = extrapolated_value / shepard_coefficient @@ -374,25 +375,25 @@ function zeroth_order_velocity_interpolation!(v, particle, shepard_coefficient, extrapolated_velocity) velocity_interpolated = extrapolated_velocity / shepard_coefficient - @inbounds for dim in eachindex(velocity_interpolated) - v[dim, particle] = velocity_interpolated[dim] + for dim in eachindex(velocity_interpolated) + @inbounds v[dim, particle] = velocity_interpolated[dim] end return v end -function first_order_scalar_interpolation!(values, particle, L_inv, +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_p = f_s[two_to_end] # f_s[2:end] as SVector + df_s = f_s[two_to_end] # `f_s[2:end]` as SVector - values[particle] = f_s[1] + dot(pos_diff, df_p) + values[particle] = f_s[1] + dot(pos_diff, df_s) return values end -function first_order_scalar_interpolation!(values, particle, L_inv, +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 @@ -402,26 +403,26 @@ function first_order_scalar_interpolation!(values, particle, L_inv, return values end -function first_order_velocity_interpolation!(v, particle, L_inv, +function first_order_velocity_extrapolation!(v, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::FirstOrderMirroring) - @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * interpolated_velocity_correction[dim, :] - df_v = f_v[two_to_end] # f_v[2:end] as SVector + 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 - v[dim, particle] = f_v[1] + dot(pos_diff, df_v) + @inbounds v[dim, particle] = f_v[1] + dot(pos_diff, df_v) end return v end -function first_order_velocity_interpolation!(v, particle, L_inv, +function first_order_velocity_extrapolation!(v, particle, L_inv, interpolated_velocity_correction, two_to_end, pos_diff, ::SimpleMirroring) - @inbounds for dim in eachindex(pos_diff) - f_v = L_inv * interpolated_velocity_correction[dim, :] + for dim in eachindex(pos_diff) + @inbounds f_v = L_inv * interpolated_velocity_correction[dim, :] - v[dim, particle] = f_v[1] + @inbounds v[dim, particle] = f_v[1] end return v @@ -502,8 +503,8 @@ function average_velocity!(v, u, system, boundary_zone::BoundaryZone{InFlow}, se @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 @@ -518,11 +519,11 @@ function project_velocity_on_plane_normal!(v, system, particle, # 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 [...]" - vel = current_velocity(v, system, particle) - vel_ = dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal + v_particle = current_velocity(v, system, particle) + v_particle_projected = dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal - @inbounds for dim in eachindex(vel) - v[dim, particle] = vel_[dim] + for dim in eachindex(vel) + @inbounds v[dim, particle] = v_particle_projected[dim] end return v diff --git a/test/examples/gpu.jl b/test/examples/gpu.jl index f6490f4a06..0acdcb8524 100644 --- a/test/examples/gpu.jl +++ b/test/examples/gpu.jl @@ -319,9 +319,7 @@ end joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), - open_boundary_model=BoundaryModelTafuni(; - mirror_method=FirstOrderMirroring(; - firstorder_tolerance=0.001f0)), + open_boundary_model=BoundaryModelTafuni(), boundary_type_in=BidirectionalFlow(), boundary_type_out=BidirectionalFlow(), reference_density_in=nothing, From d611748d8764a39813830c1c59d24246ec84ed22 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 24 Jul 2025 17:29:13 +0200 Subject: [PATCH 22/27] fix --- src/schemes/boundary/open_boundary/mirroring.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 9ffa11ba91..829670af92 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -520,9 +520,10 @@ function project_velocity_on_plane_normal!(v, system, particle, # "Because flow from the inlet interface occurs perpendicular to the boundary, # only this component of interpolated velocity is kept [...]" v_particle = current_velocity(v, system, particle) - v_particle_projected = dot(vel, boundary_zone.plane_normal) * boundary_zone.plane_normal + v_particle_projected = dot(v_particle, boundary_zone.plane_normal) * + boundary_zone.plane_normal - for dim in eachindex(vel) + for dim in eachindex(v_particle) @inbounds v[dim, particle] = v_particle_projected[dim] end From 70105bbd88a34958673052069b994e452d90a092 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 24 Jul 2025 17:35:32 +0200 Subject: [PATCH 23/27] add comments --- .../boundary/open_boundary/method_of_characteristics.jl | 4 ++++ src/schemes/boundary/open_boundary/mirroring.jl | 1 + 2 files changed, 5 insertions(+) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 061db7cd42..addc29210a 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -16,6 +16,8 @@ For more information about the method see [description below](@ref method_of_cha **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{T} extrapolate_reference_values::T @@ -238,6 +240,8 @@ 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 diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 829670af92..6412c9150e 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -497,6 +497,7 @@ 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) / length(candidates) end From 142563750e005fdbe64d93d4882015381c2d7169 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Thu, 24 Jul 2025 18:07:52 +0200 Subject: [PATCH 24/27] implement suggestions --- src/schemes/boundary/open_boundary/mirroring.jl | 11 ++++++----- test/schemes/boundary/open_boundary/mirroring.jl | 11 +++++++++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index 6412c9150e..c5f6082bb3 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -1,5 +1,5 @@ """ - FirstOrderMirroring(; firstorder_tolerance::=1f-3) + 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. @@ -10,6 +10,7 @@ to extend the gradient into the boundary zone. """ struct FirstOrderMirroring{ELTYPE} firstorder_tolerance::ELTYPE + function FirstOrderMirroring(; firstorder_tolerance::Real=1.0f-3) return new{typeof(firstorder_tolerance)}(firstorder_tolerance) end @@ -202,21 +203,21 @@ function extrapolate_values!(system, if abs(det(correction_matrix[])) >= mirror_method.firstorder_tolerance L_inv = inv(correction_matrix[]) - # pressure + # Pressure if !prescribed_pressure first_order_scalar_extrapolation!(pressure, particle, L_inv, interpolated_pressure_correction[], two_to_end, pos_diff, mirror_method) end - # density + # Density if !prescribed_density first_order_scalar_extrapolation!(density, particle, L_inv, interpolated_density_correction[], two_to_end, pos_diff, mirror_method) end - # velocity + # Velocity if !prescribed_velocity first_order_velocity_extrapolation!(v_open_boundary, particle, L_inv, interpolated_velocity_correction[], @@ -231,7 +232,7 @@ function extrapolate_values!(system, end # No else: `correction_matrix[][1, 1] <= eps()` means no fluid neighbors - # and thus no reliable interpolation, so boundary values remain at their current state + # 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] diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 4efe2fe2e7..8f658a04a2 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -261,6 +261,8 @@ @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, @@ -283,8 +285,8 @@ 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) @@ -307,6 +309,7 @@ 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) @@ -324,6 +327,10 @@ 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) @@ -341,7 +348,7 @@ smoothing_length = 1.2 * particle_spacing smoothing_kernel = WendlandC2Kernel{2}() - # Use a fluid system to interpolate the pressure + # Use a temporary fluid system just to interpolate the pressure interpolation_system = WeaklyCompressibleSPHSystem(entire_domain, ContinuityDensity(), nothing, smoothing_kernel, From 70422b98ff3f907b4c4ca74913efa81e6c46570f Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 28 Jul 2025 12:52:12 +0200 Subject: [PATCH 25/27] implement suggestions --- src/schemes/boundary/open_boundary/mirroring.jl | 4 +++- test/schemes/boundary/open_boundary/mirroring.jl | 4 ++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/schemes/boundary/open_boundary/mirroring.jl b/src/schemes/boundary/open_boundary/mirroring.jl index c5f6082bb3..3d9fde0af7 100644 --- a/src/schemes/boundary/open_boundary/mirroring.jl +++ b/src/schemes/boundary/open_boundary/mirroring.jl @@ -37,7 +37,9 @@ end """ ZerothOrderMirroring() -Fluid properties are interpolated onto ghost nodes using Shepard interpolation [Shepard1968](@cite). +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. diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index 8f658a04a2..b0592541bb 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -367,6 +367,10 @@ 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(), From 7e0d6185ce30abf538a70227affde9c31c944aaf Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 28 Jul 2025 12:53:36 +0200 Subject: [PATCH 26/27] rm ref --- docs/src/refs.bib | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index b7711bc313..2f2a636dd1 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -761,14 +761,7 @@ @book{Poling2001 publisher = {McGraw-Hill}, address = {New York} } -@article{Shepard1968, - author = {Shepard, Donald}, - booktitle = {Proceedings of the 1968 23rd ACM national conference on -}, - title = {A two-dimensional interpolation function for irregularly-spaced data}, - year = {1968}, - publisher = {ACM Press}, - doi = {10.1145/800186.810616}, -} + @article{Smagorinsky1963, author = {Smagorinsky, Joseph}, title = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment}, From fcfd47777be2ef16280af3c49acc6b7c3a00e476 Mon Sep 17 00:00:00 2001 From: LasNikas Date: Mon, 28 Jul 2025 12:59:57 +0200 Subject: [PATCH 27/27] format --- test/schemes/boundary/open_boundary/mirroring.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index b0592541bb..6e3fb9cd6f 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -367,7 +367,6 @@ 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.