Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Make Adami loop flipping optimization optional #581

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

Make Adami loop flipping optimization optional #581

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
ab18cfb
Make Adami loop flipping optimization optional
efaulhaber Jul 24, 2024
278e92e
Disable loop flipping in validation
efaulhaber Jul 24, 2024
9ef452c
Update vaildation files
efaulhaber Jul 24, 2024
83be60a
Reformat
efaulhaber Dec 12, 2024
10df125
Add option to bernoulli model as well
efaulhaber Dec 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 46 additions & 15 deletions src/schemes/boundary/dummy_particles/dummy_particles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,20 +70,34 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass,
end

@doc raw"""
AdamiPressureExtrapolation(; pressure_offset=0.0)
AdamiPressureExtrapolation(; pressure_offset=0.0, allow_loop_flipping=true)

`density_calculator` for `BoundaryModelDummyParticles`.

# Keywords
- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure
to prevent penetration, which is possible by increasing this value.
- `allow_loop_flipping=true`: Allow to flip the loop order for the pressure extrapolation.
Disable to prevent error variations between simulations with
different numbers of threads.
Usually, the first (multithreaded) loop is over the boundary
particles and the second loop over the fluid neighbors.
When the number of boundary particles is larger than
`ceil(0.5 * nthreads())` times the number of fluid particles,
it is usually more efficient to flip the loop order and loop
over the fluid particles first.
The factor depends on the number of threads, as the flipped
loop is not thread parallelizable.
This can cause error variations between simulations with
different numbers of threads.

"""
struct AdamiPressureExtrapolation{ELTYPE}
pressure_offset::ELTYPE
pressure_offset :: ELTYPE
allow_loop_flipping :: Bool

function AdamiPressureExtrapolation(; pressure_offset=0.0)
return new{eltype(pressure_offset)}(pressure_offset)
function AdamiPressureExtrapolation(; pressure_offset=0.0, allow_loop_flipping=true)
return new{eltype(pressure_offset)}(pressure_offset, allow_loop_flipping)
end
end

Expand All @@ -97,14 +111,29 @@ end
to prevent penetration, which is possible by increasing this value.
- `factor=1.0` : Setting `factor` allows to just increase the strength of the dynamic
pressure part.
- `allow_loop_flipping=true`: Allow to flip the loop order for the pressure extrapolation.
Disable to prevent error variations between simulations with
different numbers of threads.
Usually, the first (multithreaded) loop is over the boundary
particles and the second loop over the fluid neighbors.
When the number of boundary particles is larger than
`ceil(0.5 * nthreads())` times the number of fluid particles,
it is usually more efficient to flip the loop order and loop
over the fluid particles first.
The factor depends on the number of threads, as the flipped
loop is not thread parallelizable.
This can cause error variations between simulations with
different numbers of threads.

"""
struct BernoulliPressureExtrapolation{ELTYPE}
pressure_offset :: ELTYPE
factor :: ELTYPE
pressure_offset :: ELTYPE
factor :: ELTYPE
allow_loop_flipping :: Bool

function BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=0.0)
return new{eltype(pressure_offset)}(pressure_offset, factor)
function BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=0.0,
allow_loop_flipping=true)
return new{eltype(pressure_offset)}(pressure_offset, factor, allow_loop_flipping)
end
end

Expand Down Expand Up @@ -339,6 +368,7 @@ function compute_pressure!(boundary_model,
BernoulliPressureExtrapolation},
system, v, u, v_ode, u_ode, semi)
(; pressure, cache, viscosity) = boundary_model
(; allow_loop_flipping) = boundary_model.density_calculator

set_zero!(pressure)

Expand All @@ -355,23 +385,24 @@ function compute_pressure!(boundary_model,
neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system)

# This is an optimization for simulations with large and complex boundaries.
# Especially, in 3D simulations with large and/or complex structures outside
# Especially in 3D simulations with large and/or complex structures outside
# of areas with permanent flow.
# Note: The version iterating neighbors first is not thread parallelizable.
# The factor is based on the achievable speed-up of the thread parallelizable version.
if nparticles(system) >
ceil(Int, 0.5 * Threads.nthreads()) * nparticles(neighbor_system)
nhs = get_neighborhood_search(neighbor_system, system, semi)
if allow_loop_flipping &&
nparticles(system) > div(Threads.nthreads(), 2) * nparticles(neighbor_system)

# Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries
# Loop over fluid particles and then the neighboring boundary particles
# to extrapolate fluid pressure to the boundaries.
nhs = get_neighborhood_search(neighbor_system, system, semi)
boundary_pressure_extrapolation_neighbor!(boundary_model, system,
neighbor_system,
system_coords, neighbor_coords, v,
v_neighbor_system, nhs)
else
# Loop over boundary particles and then the neighboring fluid particles
# to extrapolate fluid pressure to the boundaries.
nhs = get_neighborhood_search(system, neighbor_system, semi)

# Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries
boundary_pressure_extrapolation!(boundary_model, system,
neighbor_system,
system_coords, neighbor_coords, v,
Expand Down
4 changes: 4 additions & 0 deletions validation/dam_break_2d/validation_dam_break_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,9 +98,12 @@ fluid_system_edac = EntropicallyDampedSPHSystem(tank_edac.fluid, smoothing_kerne
pressure_acceleration=nothing,
acceleration=(0.0, -gravity))

# Disable loop flipping to produce consistent results over different thread numbers
boundary_density_calculator = AdamiPressureExtrapolation(allow_loop_flipping=false)
trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
fluid_particle_spacing=particle_spacing,
smoothing_length=smoothing_length, smoothing_kernel=smoothing_kernel,
boundary_density_calculator=boundary_density_calculator,
boundary_layers=4, state_equation=nothing,
solution_prefix="validation_" * method * "_" * formatted_string,
extra_callback=postprocessing_cb, tspan=tspan,
Expand Down Expand Up @@ -137,6 +140,7 @@ postprocessing_cb = PostprocessCallback(; dt=0.02, output_directory="out",
trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "dam_break_2d.jl"),
fluid_particle_spacing=particle_spacing,
smoothing_length=smoothing_length, smoothing_kernel=smoothing_kernel,
boundary_density_calculator=boundary_density_calculator,
boundary_layers=4,
solution_prefix="validation_" * method * "_" * formatted_string,
extra_callback=postprocessing_cb, tspan=tspan,
Expand Down
Loading
Loading