From c5dfd077efa620685032788467cf9c759b6fd2a8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Jun 2024 09:43:23 +0200 Subject: [PATCH 1/8] Precompute git hash to avoid overhead (#542) * Precompute git hash * Store git hash in callbacks * Fix --------- Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/callbacks/post_process.jl | 15 ++++++++++----- src/callbacks/solution_saving.jl | 16 ++++++++-------- src/general/file_system.jl | 31 ------------------------------- src/general/general.jl | 4 +++- src/util.jl | 18 ++++++++++++++++++ src/visualization/write2vtk.jl | 27 ++++++++++++++++++--------- 6 files changed, 57 insertions(+), 54 deletions(-) delete mode 100644 src/general/file_system.jl diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index e2e98a2a0..d68168258 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -73,6 +73,7 @@ struct PostprocessCallback{I, F} append_timestamp :: Bool write_csv :: Bool write_json :: Bool + git_hash :: Ref{String} end function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=true, @@ -94,7 +95,8 @@ function PostprocessCallback(; interval::Integer=0, dt=0.0, exclude_boundary=tru post_callback = PostprocessCallback(interval, write_file_interval, Dict{String, Vector{Any}}(), Float64[], exclude_boundary, funcs, filename, output_directory, - append_timestamp, write_csv, write_json) + append_timestamp, write_csv, write_json, + Ref("UnknownVersion")) if dt > 0 # Add a `tstop` every `dt`, and save the final solution return PeriodicCallback(post_callback, dt, @@ -209,9 +211,12 @@ function initialize_postprocess_callback!(cb, u, t, integrator) end function initialize_postprocess_callback!(cb::PostprocessCallback, u, t, integrator) + cb.git_hash[] = compute_git_hash() + # Apply the callback cb(integrator) - return nothing + + return cb end # `condition` with interval @@ -281,7 +286,7 @@ function write_postprocess_callback(pp::PostprocessCallback) mkpath(pp.output_directory) data = Dict{String, Any}() - write_meta_data!(data) + write_meta_data!(data, pp.git_hash[]) prepare_series_data!(data, pp) time_stamp = "" @@ -331,9 +336,9 @@ function create_series_dict(values, times, system_name="") "time" => times) end -function write_meta_data!(data) +function write_meta_data!(data, git_hash) meta_data = Dict("solver_name" => "TrixiParticles.jl", - "solver_version" => get_git_hash(), + "solver_version" => git_hash, "julia_version" => string(VERSION)) data["meta"] = meta_data diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 5163199d6..04552d57e 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -77,6 +77,7 @@ mutable struct SolutionSavingCallback{I, CQ} max_coordinates :: Float64 custom_quantities :: CQ latest_saved_iter :: Int + git_hash :: Ref{String} end function SolutionSavingCallback(; interval::Integer=0, dt=0.0, @@ -102,7 +103,7 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_initial_solution, save_final_solution, write_meta_data, verbose, output_directory, prefix, max_coordinates, custom_quantities, - -1) + -1, Ref("UnknownVersion")) if length(save_times) > 0 # See the large comment below for an explanation why we use `finalize` here @@ -131,6 +132,7 @@ end function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, integrator) # Reset `latest_saved_iter` solution_callback.latest_saved_iter = -1 + solution_callback.git_hash[] = compute_git_hash() # Save initial solution if solution_callback.save_initial_solution @@ -156,7 +158,7 @@ end # `affect!` function (solution_callback::SolutionSavingCallback)(integrator) - (; interval, output_directory, custom_quantities, write_meta_data, + (; interval, output_directory, custom_quantities, write_meta_data, git_hash, verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback vu_ode = integrator.u @@ -178,12 +180,10 @@ function (solution_callback::SolutionSavingCallback)(integrator) println("Writing solution to $output_directory at t = $(integrator.t)") end - @trixi_timeit timer() "save solution" trixi2vtk(vu_ode, semi, integrator.t; iter=iter, - output_directory=output_directory, - prefix=prefix, - write_meta_data=write_meta_data, - max_coordinates=max_coordinates, - custom_quantities...) + @trixi_timeit timer() "save solution" trixi2vtk(vu_ode, semi, integrator.t; + iter, output_directory, prefix, + write_meta_data, git_hash=git_hash[], + max_coordinates, custom_quantities...) # Tell OrdinaryDiffEq that `u` has not been modified u_modified!(integrator, false) diff --git a/src/general/file_system.jl b/src/general/file_system.jl deleted file mode 100644 index 0cc4cc56d..000000000 --- a/src/general/file_system.jl +++ /dev/null @@ -1,31 +0,0 @@ -vtkname(system::FluidSystem) = "fluid" -vtkname(system::SolidSystem) = "solid" -vtkname(system::BoundarySystem) = "boundary" - -function system_names(systems) - # Add `_i` to each system name, where `i` is the index of the corresponding - # system type. - # `["fluid", "boundary", "boundary"]` becomes `["fluid_1", "boundary_1", "boundary_2"]`. - cnames = systems .|> vtkname - filenames = [string(cnames[i], "_", count(==(cnames[i]), cnames[1:i])) - for i in eachindex(cnames)] - return filenames -end - -function get_git_hash() - pkg_directory = pkgdir(@__MODULE__) - git_directory = joinpath(pkg_directory, ".git") - - # Check if the .git directory exists - if !isdir(git_directory) - return "UnknownVersion" - end - - try - git_cmd = Cmd(`git describe --tags --always --first-parent --dirty`, - dir=pkg_directory) - return string(readchomp(git_cmd)) - catch e - return "UnknownVersion" - end -end diff --git a/src/general/general.jl b/src/general/general.jl index dfe52e06a..4a60cb52f 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -7,12 +7,15 @@ GPUSystem = System{NDIMS, Nothing} where {NDIMS} abstract type FluidSystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::FluidSystem) = "fluid" +vtkname(system::FluidSystem) = "fluid" abstract type SolidSystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::SolidSystem) = "solid" +vtkname(system::SolidSystem) = "solid" abstract type BoundarySystem{NDIMS, IC} <: System{NDIMS, IC} end timer_name(::BoundarySystem) = "boundary" +vtkname(system::BoundarySystem) = "boundary" @inline function set_zero!(du) du .= zero(eltype(du)) @@ -29,6 +32,5 @@ include("smoothing_kernels.jl") include("initial_condition.jl") include("system.jl") include("interpolation.jl") -include("file_system.jl") include("custom_quantities.jl") include("neighborhood_search.jl") diff --git a/src/util.jl b/src/util.jl index 2dc86a267..b3cbfb764 100644 --- a/src/util.jl +++ b/src/util.jl @@ -206,3 +206,21 @@ end function type2string(type) return string(nameof(typeof(type))) end + +function compute_git_hash() + pkg_directory = pkgdir(@__MODULE__) + git_directory = joinpath(pkg_directory, ".git") + + # Check if the .git directory exists + if !isdir(git_directory) + return "UnknownVersion" + end + + try + git_cmd = Cmd(`git describe --tags --always --first-parent --dirty`, + dir=pkg_directory) + return string(readchomp(git_cmd)) + catch e + return "UnknownVersion" + end +end diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index bd1658e32..a05e6727e 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -1,3 +1,13 @@ +function system_names(systems) + # Add `_i` to each system name, where `i` is the index of the corresponding + # system type. + # `["fluid", "boundary", "boundary"]` becomes `["fluid_1", "boundary_1", "boundary_2"]`. + cnames = vtkname.(systems) + filenames = [string(cnames[i], "_", count(==(cnames[i]), cnames[1:i])) + for i in eachindex(cnames)] + return filenames +end + """ trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", write_meta_data=true, max_coordinates=Inf, custom_quantities...) @@ -39,13 +49,14 @@ trixi2vtk(sol.u[end], semi, 0.0, iter=1, my_custom_quantity=kinetic_energy) ``` """ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="", - write_meta_data=true, max_coordinates=Inf, custom_quantities...) + write_meta_data=true, git_hash=compute_git_hash(), + max_coordinates=Inf, custom_quantities...) (; systems) = semi v_ode, u_ode = vu_ode.x # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t) + @trixi_timeit timer() "update systems" update_systems_and_nhs(v_ode, u_ode, semi, t) filenames = system_names(systems) @@ -57,17 +68,15 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix periodic_box = get_neighborhood_search(system, semi).periodic_box trixi2vtk(v, u, t, system, periodic_box; - output_directory=output_directory, - system_name=filenames[system_index], iter=iter, prefix=prefix, - write_meta_data=write_meta_data, max_coordinates=max_coordinates, - custom_quantities...) + system_name=filenames[system_index], output_directory, iter, prefix, + write_meta_data, git_hash, max_coordinates, custom_quantities...) end end # Convert data for a single TrixiParticle system to VTK format function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix="", iter=nothing, system_name=vtkname(system), write_meta_data=true, - max_coordinates=Inf, + max_coordinates=Inf, git_hash=compute_git_hash(), custom_quantities...) mkpath(output_directory) @@ -98,7 +107,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix end end - vtk_grid(file, points, cells) do vtk + @trixi_timeit timer() "write to vtk" vtk_grid(file, points, cells) do vtk # dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) @@ -107,7 +116,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix vtk["time"] = t if write_meta_data - vtk["solver_version"] = get_git_hash() + vtk["solver_version"] = git_hash vtk["julia_version"] = string(VERSION) end From 5deebdd94bd01e785d3e67d742689f5bd606d061 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 17 Jun 2024 10:17:01 +0200 Subject: [PATCH 2/8] Move timers and `@trixi_timeit` to TrixiBase.jl (#543) * Move timers and `@trixi_timeit` to TrixiBase.jl * Require latest version of TrixiBase.jl * Fix tests --------- Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- Project.toml | 2 +- examples/n_body/n_body_benchmark_trixi.jl | 4 +-- examples/n_body/n_body_solar_system.jl | 4 +-- src/TrixiParticles.jl | 3 +- src/util.jl | 42 ----------------------- test/count_allocations.jl | 4 +-- 6 files changed, 9 insertions(+), 50 deletions(-) diff --git a/Project.toml b/Project.toml index ea1bc9aa5..bff9874c1 100644 --- a/Project.toml +++ b/Project.toml @@ -42,7 +42,7 @@ SciMLBase = "1, 2" StaticArrays = "1" StrideArrays = "0.1" TimerOutputs = "0.5" -TrixiBase = "0.1" +TrixiBase = "0.1.3" PointNeighbors = "0.2" WriteVTK = "1" julia = "1.9" diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index 04c1f1af6..b47f689e7 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -108,7 +108,7 @@ end filename = tempname() open(filename, "w") do f redirect_stderr(f) do - TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + TrixiParticles.disable_debug_timings() end end @@ -124,6 +124,6 @@ end # Enable timers again open(filename, "w") do f redirect_stderr(f) do - TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + TrixiParticles.enable_debug_timings() end end diff --git a/examples/n_body/n_body_solar_system.jl b/examples/n_body/n_body_solar_system.jl index 6fa860089..b4da85aeb 100644 --- a/examples/n_body/n_body_solar_system.jl +++ b/examples/n_body/n_body_solar_system.jl @@ -49,7 +49,7 @@ callbacks = CallbackSet(info_callback, saving_callback) filename = tempname() open(filename, "w") do f redirect_stderr(f) do - TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + TrixiParticles.disable_debug_timings() end end @@ -63,6 +63,6 @@ sol = solve(ode, SymplecticEuler(), # Enable timers again open(filename, "w") do f redirect_stderr(f) do - TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + TrixiParticles.enable_debug_timings() end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index f170e08d8..2fe865d10 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -21,7 +21,8 @@ using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified! using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer! -using TrixiBase: trixi_include +using TrixiBase: trixi_include, @trixi_timeit, timer, timeit_debug_enabled, + disable_debug_timings, enable_debug_timings @reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch using PointNeighbors: PointNeighbors, for_particle_neighbor using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save diff --git a/src/util.jl b/src/util.jl index b3cbfb764..338a94415 100644 --- a/src/util.jl +++ b/src/util.jl @@ -26,48 +26,6 @@ function print_startup_message() println(s) end -# Enable debug timings `@trixi_timeit timer() "name" stuff...`. -# This allows us to disable timings completely by executing -# `TimerOutputs.disable_debug_timings(TrixiParticles)` -# and to enable them again by executing -# `TimerOutputs.enable_debug_timings(TrixiParticles)` -timeit_debug_enabled() = true - -# Store main timer for global timing of functions -const main_timer = TimerOutput() - -# Always call timer() to hide implementation details -timer() = main_timer - -# @trixi_timeit timer() "some label" expression -# -# Basically the same as a special case of `@timeit_debug` from -# [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl), -# but without `try ... finally ... end` block. Thus, it's not exception-safe, -# but it also avoids some related performance problems. Since we do not use -# exception handling in TrixiParticles, that's not really an issue. -# -# Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl). -macro trixi_timeit(timer_output, label, expr) - timeit_block = quote - if timeit_debug_enabled() - local to = $(esc(timer_output)) - local enabled = to.enabled - if enabled - local accumulated_data = $(TimerOutputs.push!)(to, $(esc(label))) - end - local b0 = $(TimerOutputs.gc_bytes)() - local t0 = $(TimerOutputs.time_ns)() - end - local val = $(esc(expr)) - if timeit_debug_enabled() && enabled - $(TimerOutputs.do_accumulate!)(accumulated_data, t0, b0) - $(TimerOutputs.pop!)(to) - end - val - end -end - """ @threaded for ... end diff --git a/test/count_allocations.jl b/test/count_allocations.jl index 3b5110ba1..1c597f9a0 100644 --- a/test/count_allocations.jl +++ b/test/count_allocations.jl @@ -45,7 +45,7 @@ function count_rhs_allocations(sol, semi) try # Disable timers, which cause extra allocations - TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles) + TrixiParticles.disable_debug_timings() # Disable multithreading, which causes extra allocations return disable_polyester_threads() do @@ -57,7 +57,7 @@ function count_rhs_allocations(sol, semi) end finally # Enable timers again - @invokelatest TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles) + @invokelatest TrixiParticles.enable_debug_timings() end end From 4af461da201fa6f817ce11d92437fb0c66beecb2 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Wed, 19 Jun 2024 09:43:16 +0200 Subject: [PATCH 3/8] Add `trixi2vtk` wrapper for `InitialCondition` (#546) * add wrapper * modify kwarg * implement suggestions --- src/visualization/write2vtk.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index a05e6727e..2a4c9f207 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -183,6 +183,32 @@ function trixi2vtk(coordinates; output_directory="out", prefix="", filename="coo return file end +""" + trixi2vtk(initial_condition::InitialCondition; output_directory="out", + prefix="", filename="initial_condition", custom_quantities...) + +Convert [`InitialCondition`](@ref) data to VTK format. + +# Arguments +- `initial_condition`: [`InitialCondition`](@ref) to be saved. + +# Keywords +- `output_directory="out"`: Output directory path. +- `prefix=""`: Prefix for the output file. +- `filename="coordinates"`: Name of the output file. + +# Returns +- `file::AbstractString`: Path to the generated VTK file. +""" +function trixi2vtk(initial_condition::InitialCondition; output_directory="out", + prefix="", filename="initial_condition", custom_quantities...) + (; coordinates, velocity, density, mass, pressure) = initial_condition + + return trixi2vtk(coordinates; output_directory, prefix, filename, + density=density, initial_velocity=velocity, mass=mass, + pressure=pressure) +end + function write2vtk!(vtk, v, u, t, system; write_meta_data=true) vtk["velocity"] = view(v, 1:ndims(system), :) From 39bc7985b8df42beae8b58984b59a9a251544bfa Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 19 Jun 2024 10:44:19 +0200 Subject: [PATCH 4/8] Example on how to save interpolated planes using postprocessing callback (#462) * add example * replace MSE with MRE * Revert "replace MSE with MRE" This reverts commit 1e0090794f7c13ea02b5c3ae298514caab93f0a4. * revert * update * format * fix and reduce runtime * Fix comment * implement suggestions * format --------- Co-authored-by: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> --- .../postprocessing/interpolation_plane.jl | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/examples/postprocessing/interpolation_plane.jl b/examples/postprocessing/interpolation_plane.jl index 5c94df996..ea273d0fc 100644 --- a/examples/postprocessing/interpolation_plane.jl +++ b/examples/postprocessing/interpolation_plane.jl @@ -96,3 +96,25 @@ plot_3d = Plots.plot(scatter_3d, xlabel="X", ylabel="Y", zlabel="Z", combined_plot = Plots.plot(plot1, plot2, plot3, plot_3d, layout=(2, 2), size=(1000, 1500), margin=3mm) + +# If we want to save planes at regular intervals, we can use the postprocessing callback. +function save_interpolated_plane(v, u, t, system) + # Size of the patch to be interpolated + interpolation_start = [0.0, 0.0] + interpolation_end = [tank_size[1], tank_size[2]] + + # The resolution the plane is interpolated to. In this case twice the original resolution. + resolution = 0.5 * fluid_particle_spacing + + file_id = ceil(Int, t * 10000) + interpolate_plane_2d_vtk(interpolation_start, interpolation_end, resolution, + semi, system, v, u, filename="plane_$file_id.vti") + return nothing +end + +save_interpolation_cb = PostprocessCallback(; dt=0.1, write_file_interval=0, + save_interpolated_plane) + +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), tspan=(0.0, 0.2), + extra_callback=save_interpolation_cb, fluid_particle_spacing=0.01) From ebd5cd98758a5f6d1afa474204621359a689bc64 Mon Sep 17 00:00:00 2001 From: Niklas Neher <73897120+LasNikas@users.noreply.github.com> Date: Wed, 19 Jun 2024 15:08:41 +0200 Subject: [PATCH 5/8] Add open boundaries (#442) * implement `ExtrudeFace` * use `sample_face` in `interpolate_plane_3d` * add docs * add tests * fix bug * create `OpenBoundarySPHSystem` * export OpenBoundarySPHSystem * fix bug * add `tlsph=true` * fix tuple bug * fix typo * fix again * fix for tests * fix layers * calculate particle spacing differently * add show-tests * add `within_boundary_zone` * change `floor` to `ceil` * boundary zone tests * add `evaluate_characteristics` * add tests * add reference functions * add `update_quantities!` * add buffer * add `SystemBuffer` tests * cosmetic change in tests * change `round` to `isapprox` * add update callback * add `UpdateCallback` * add `update_open_boundaries` * fix bugs * add timers * write prescribed quantities * improve dispatching * add example `pipe_flow_2d.jl` * add docs `UpdateCallback` * add `UpdateCallback` * fix typo * apply formatter * modify system buffer * add docstring * add docs * fix tests * fix typos * apply formatter * fix tests * add comment * fix tests * remove update bool * adapt docstring * implement suggestions * fix test * add check if callback is used * add comments in example file * generic types * merge main * undo combining compact support with DEM * remove density and pressure setter * implement suggestions * fix tests * Revert "remove density and pressure setter" This reverts commit 0203a3a21a3700a4646a327c1d7eb137fa9bb43f. * implement suggestions * rework in- and outflow * rename setter functions * rename function * adapt tests * add docs * rename function * apply formatter * fix tests * fix bug in plot recipes * apply formatter * implement suggestions for `boundary_zone.jl` * rename kwarg * implement suggestions * implement suggestions * modify `check_domain!` * add comments and check for fluid system with buffer * modify error message * fix tests * fix `callback_used` * using Measurements * fix tests * sqrt(eps()) * fix tests again * add doc to `boundary.md` * link single fluid system * fix tests * implement suggestions * fix tests * change the order of functions * modify example * fix tests * fix test again * implement doc suggestions * modify dos * fix typos * implement suggestions * implement suggestions * add EDAC to the example * implement suggestions * modify tests * adapt docs * add `initial_callback_flag` * add to `NEWS.md` and `README.md` * implement suggestions * fix typo --- NEWS.md | 10 +- README.md | 2 +- docs/src/systems/boundary.md | 83 +++ examples/fluid/pipe_flow_2d.jl | 131 +++++ examples/n_body/n_body_system.jl | 3 +- src/TrixiParticles.jl | 3 +- src/callbacks/density_reinit.jl | 2 +- src/callbacks/post_process.jl | 2 +- src/callbacks/solution_saving.jl | 2 +- src/callbacks/update.jl | 25 +- src/general/buffer.jl | 87 ++++ src/general/density_calculators.jl | 13 + src/general/general.jl | 1 + src/general/semidiscretization.jl | 56 +- src/general/system.jl | 20 +- .../boundary/open_boundary/boundary_zones.jl | 340 ++++++++++++ src/schemes/boundary/open_boundary/system.jl | 483 ++++++++++++++++++ src/schemes/boundary/rhs.jl | 4 +- src/schemes/boundary/system.jl | 24 +- .../fluid/entropically_damped_sph/system.jl | 39 +- src/schemes/fluid/fluid.jl | 4 + src/schemes/fluid/viscosity.jl | 13 + .../fluid/weakly_compressible_sph/system.jl | 28 +- src/schemes/schemes.jl | 3 + .../solid/discrete_element_method/system.jl | 4 +- .../solid/total_lagrangian_sph/system.jl | 6 +- src/setups/extrude_geometry.jl | 6 +- src/visualization/recipes_plots.jl | 8 +- src/visualization/write2vtk.jl | 45 +- test/examples/examples.jl | 8 + test/general/buffer.jl | 68 +++ test/general/general.jl | 1 + test/general/semidiscretization.jl | 3 + .../boundary/open_boundary/boundary_zone.jl | 221 ++++++++ .../open_boundary/characteristic_variables.jl | 134 +++++ .../boundary/open_boundary/open_boundary.jl | 2 + test/schemes/schemes.jl | 1 + test/systems/open_boundary_system.jl | 90 ++++ test/systems/systems.jl | 1 + 39 files changed, 1913 insertions(+), 63 deletions(-) create mode 100644 examples/fluid/pipe_flow_2d.jl create mode 100644 src/general/buffer.jl create mode 100644 src/schemes/boundary/open_boundary/boundary_zones.jl create mode 100644 src/schemes/boundary/open_boundary/system.jl create mode 100644 test/general/buffer.jl create mode 100644 test/schemes/boundary/open_boundary/boundary_zone.jl create mode 100644 test/schemes/boundary/open_boundary/characteristic_variables.jl create mode 100644 test/schemes/boundary/open_boundary/open_boundary.jl create mode 100644 test/systems/open_boundary_system.jl diff --git a/NEWS.md b/NEWS.md index ffe272e8a..5d7a295d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,8 @@ # Changelog TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) -used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. +used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. ## Version 0.2.x @@ -14,11 +14,15 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between ### Deprecated +## Version 0.1.3 + +### Added +Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC. ## Version 0.1.2 ### Added -A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH +A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH. ## Version 0.1.1 diff --git a/README.md b/README.md index b40e68fbc..76b288ed2 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ Its features include: ## Features - Incompressible Navier-Stokes - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC) - - Models: Surface Tension + - Models: Surface Tension, Open Boundaries - Solid-body mechanics - Methods: Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM) - Fluid-Structure Interaction diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 8b374a71e..ae4b9b676 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -234,3 +234,86 @@ Pages = [joinpath("schemes", "boundary", "monaghan_kajtar", "monaghan_kajtar.jl" - Alireza Valizadeh, Joseph J. Monaghan. "A study of solid wall models for weakly compressible SPH." In: Journal of Computational Physics 300 (2015), pages 5–19. [doi: 10.1016/J.JCP.2015.07.033](https://doi.org/10.1016/J.JCP.2015.07.033) + +# [Open Boundaries](@id open_boundary) + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "boundary", "open_boundary", "system.jl")] +``` + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "boundary", "open_boundary", "boundary_zones.jl")] +``` + +### [Method of characteristics](@id method_of_characteristics) + +The difficulty in non-reflecting boundary conditions, also called open boundaries, is to determine +the appropriate boundary values of the exact characteristics of the Euler equations. +Assuming the flow near the boundaries is normal to the boundary +and free of shock waves and significant viscous effects, it can be shown that three characteristic variables exist: + +- ``J_1``, associated with convection of entropy and propagates at flow velocity, +- ``J_2``, downstream-running characteristics, +- ``J_3``, upstream-running characteristics. + +Giles (1990) derived those variables based on a linearized set of governing equations: +```math +J_1 = -c_s^2 (\rho - \rho_{\text{ref}}) + (p - p_{\text{ref}}) +``` +```math +J_2 = \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}}) +``` +```math +J_3 = - \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}}) +``` +where the subscript "ref" denotes the reference flow near the boundaries, which can be prescribed. + +Specifying the reference variables is **not** equivalent to prescription of ``\rho``, ``v`` and ``p`` +directly, since the perturbation from the reference flow is allowed. + +Lastiwka et al. (2009) applied the method of characteristic to SPH and determined the number of variables that should be +**prescribed** at the boundary and the number which should be **propagated** from the fluid domain to the boundary: + +- For an **inflow** boundary: + - Prescribe *downstream*-running characteristics ``J_1`` and ``J_2`` + - Transmit ``J_3`` from the fluid domain (allow ``J_3`` to propagate upstream to the boundary). + +- For an **outflow** boundary: + - Prescribe *upstream*-running characteristic ``J_3`` + - Transmit ``J_1`` and ``J_2`` from the fluid domain. + +Prescribing is done by simply setting the characteristics to zero. To transmit the characteristics from the fluid +domain, or in other words, to carry the information of the fluid to the boundaries, Negi et al. (2020) use a Shepard Interpolation +```math +f_i = \frac{\sum_j^N f_j W_{ij}}{\sum_j^N W_{ij}}, +``` +where the ``i``-th particle is a boundary particle, ``f`` is either ``J_1``, ``J_2`` or ``J_3`` and ``N`` is the set of +neighboring fluid particles. + +To express pressure ``p``, density ``\rho`` and velocity ``v`` as functions of the characteristic variables, the system of equations +from the characteristic variables is inverted and gives +```math + \rho - \rho_{\text{ref}} = \frac{1}{c_s^2} \left( -J_1 + \frac{1}{2} J_2 + \frac{1}{2} J_3 \right), +``` +```math +u - u_{\text{ref}}= \frac{1}{2\rho c_s} \left( J_2 - J_3 \right), +``` +```math +p - p_{\text{ref}} = \frac{1}{2} \left( J_2 + J_3 \right). +``` +With ``J_1``, ``J_2`` and ``J_3`` determined, we can easily solve for the actual variables for each particle. + +### References +- M. B. Giles. "Nonreflecting boundary conditions for Euler equation calculations". + In: AIAA Journal, 28.12 pages 2050--2058. + [doi: 10.2514/3.10521](https://doi.org/10.2514/3.10521) +- M. Lastiwka, M. Basa, N. J. Quinlan. + "Permeable and non-reflecting boundary conditions in SPH". + In: International Journal for Numerical Methods in Fluids 61, (2009), pages 709--724. + [doi: 10.1002/fld.1971](https://doi.org/10.1002/fld.1971) +- P. Negi, P. Ramachandran, A. Haftu. + "An improved non-reflecting outlet boundary condition for weakly-compressible SPH". + In: Computer Methods in Applied Mechanics and Engineering 367, (2020), pages 113--119. + [doi: 10.1016/j.cma.2020.113119](https://doi.org/10.1016/j.cma.2020.113119) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl new file mode 100644 index 000000000..3d642696e --- /dev/null +++ b/examples/fluid/pipe_flow_2d.jl @@ -0,0 +1,131 @@ +# 2D channel flow simulation with open boundaries. + +using TrixiParticles +using OrdinaryDiffEq + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# 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 = 6 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +# Boundary geometry and initial fluid particle positions +domain_size = (1.0, 0.4) + +flow_direction = [1.0, 0.0] +reynolds_number = 100 +const prescribed_velocity = 2.0 + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2]) + +fluid_density = 1000.0 + +# For this particular example, it is necessary to have a background pressure. +# Otherwise the suction at the outflow is to big and the simulation becomes unstable. +pressure = 1000.0 + +sound_speed = 10 * prescribed_velocity + +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7, background_pressure=pressure) + +pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density, + pressure=pressure, n_layers=boundary_layers, + 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 = 4 * pipe.n_particles_per_dimension[2] + +# ========================================================================================== +# ==== Fluid +smoothing_length = 3.0 * particle_spacing +smoothing_kernel = WendlandC2Kernel{2}() + +fluid_density_calculator = ContinuityDensity() + +kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number + +viscosity = ViscosityAdami(nu=kinematic_viscosity) + +fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length, + sound_speed, viscosity=viscosity, + density_calculator=fluid_density_calculator, + buffer_size=n_buffer_particles) + +# Alternatively the WCSPH scheme can be used +# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed) +# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0) + +# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator, +# state_equation, smoothing_kernel, +# smoothing_length, viscosity=viscosity, +# buffer_size=n_buffer_particles) + +# ========================================================================================== +# ==== Open Boundary +function velocity_function(pos, t) + # Use this for a time-dependent inflow velocity + # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) + + return SVector(prescribed_velocity, 0.0) +end + +inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, + open_boundary_layers, density=fluid_density, particle_spacing) + +open_boundary_in = OpenBoundarySPHSystem(inflow; sound_speed, fluid_system, + buffer_size=n_buffer_particles, + reference_pressure=pressure, + reference_velocity=velocity_function) + +outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), + flow_direction, open_boundary_layers, density=fluid_density, + particle_spacing) + +open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system, + buffer_size=n_buffer_particles, + reference_pressure=pressure, + reference_velocity=velocity_function) + +# ========================================================================================== +# ==== Boundary + +boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass, + AdamiPressureExtrapolation(), + state_equation=state_equation, + #viscosity=ViscosityAdami(nu=1e-4), + smoothing_kernel, smoothing_length) + +boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model) + +# ========================================================================================== +# ==== Simulation +semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out, + boundary_system) + +ode = semidiscretize(semi, tspan) + +info_callback = InfoCallback(interval=100) +saving_callback = SolutionSavingCallback(dt=0.02, prefix="") + +callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback()) + +sol = solve(ode, RDPK3SpFSAL35(), + abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration) + reltol=1e-3, # 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/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index caf40a62e..c7226357a 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -5,12 +5,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothin initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] G :: ELTYPE + buffer :: Nothing function NBodySystem(initial_condition, G) mass = copy(initial_condition.mass) new{size(initial_condition.coordinates, 1), - eltype(mass)}(initial_condition, mass, G) + eltype(mass)}(initial_condition, mass, G, nothing) end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 2fe865d10..aa97ab720 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -44,7 +44,8 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - BoundarySPHSystem, DEMSystem, BoundaryDEMSystem + BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow, + OutFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback export ContinuityDensity, SummationDensity diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 5c8a7b5e1..74c6fe5c9 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -74,7 +74,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr # Update systems to compute quantities like density and pressure. semi = integrator.p v_ode, u_ode = u.x - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Apply the callback. cb(integrator) diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl index d68168258..9c39f3eb4 100644 --- a/src/callbacks/post_process.jl +++ b/src/callbacks/post_process.jl @@ -236,7 +236,7 @@ function (pp::PostprocessCallback)(integrator) new_data = false # Update systems to compute quantities like density and pressure - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) foreach_system(semi) do system if system isa BoundarySystem && pp.exclude_boundary diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl index 04552d57e..468a72a1c 100644 --- a/src/callbacks/solution_saving.jl +++ b/src/callbacks/solution_saving.jl @@ -139,7 +139,7 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in # Update systems to compute quantities like density and pressure semi = integrator.p v_ode, u_ode = u.x - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Apply the callback solution_callback(integrator) diff --git a/src/callbacks/update.jl b/src/callbacks/update.jl index ca3b5f114..3595d5c36 100644 --- a/src/callbacks/update.jl +++ b/src/callbacks/update.jl @@ -52,7 +52,16 @@ function initial_update!(cb, u, t, integrator) initial_update!(cb.affect!, u, t, integrator) end -initial_update!(cb::UpdateCallback, u, t, integrator) = cb(integrator) +function initial_update!(cb::UpdateCallback, u, t, integrator) + semi = integrator.p + + # Tell systems that `UpdateCallback` is used + foreach_system(semi) do system + update_callback_used!(system) + end + + return cb(integrator) +end # `condition` function (update_callback!::UpdateCallback)(u, t, integrator) @@ -69,16 +78,12 @@ function (update_callback!::UpdateCallback)(integrator) # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - update_systems_and_nhs(v_ode, u_ode, semi, t) + update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true) # Other updates might be added here later (e.g. Transport Velocity Formulation). - # @trixi_timeit timer() "update open boundary" foreach_system(semi) do system - # update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) - # end - # - # @trixi_timeit timer() "update TVF" foreach_system(semi) do system - # update_transport_velocity_eachstep!(system, v_ode, u_ode, semi, t) - # end + @trixi_timeit timer() "update open boundary" foreach_system(semi) do system + update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) + end # Tell OrdinaryDiffEq that `u` has been modified u_modified!(integrator, true) @@ -128,3 +133,5 @@ function Base.show(io::IO, ::MIME"text/plain", summary_box(io, "UpdateCallback", setup) end end + +update_callback_used!(system) = system diff --git a/src/general/buffer.jl b/src/general/buffer.jl new file mode 100644 index 000000000..b12cf05f3 --- /dev/null +++ b/src/general/buffer.jl @@ -0,0 +1,87 @@ +struct SystemBuffer{V} + active_particle :: BitVector + eachparticle :: V # Vector{Int} + buffer_size :: Int + + function SystemBuffer(active_size, buffer_size::Integer) + active_particle = vcat(trues(active_size), falses(buffer_size)) + eachparticle = collect(1:active_size) + + return new{typeof(eachparticle)}(active_particle, eachparticle, buffer_size) + end +end + +allocate_buffer(initial_condition, buffer) = initial_condition + +function allocate_buffer(initial_condition, buffer::SystemBuffer) + (; buffer_size) = buffer + + # Initialize particles far away from simulation domain + coordinates = fill(1e16, ndims(initial_condition), buffer_size) + + if all(rho -> isapprox(rho, first(initial_condition.density), atol=eps(), rtol=eps()), + initial_condition.density) + density = first(initial_condition.density) + else + throw(ArgumentError("`initial_condition.density` needs to be constant when using `SystemBuffer`")) + end + + particle_spacing = initial_condition.particle_spacing + + buffer_ic = InitialCondition(; coordinates, density, particle_spacing) + + return union(initial_condition, buffer_ic) +end + +@inline update_system_buffer!(buffer::Nothing) = buffer + +# TODO `resize` allocates. Find a non-allocating version +@inline function update_system_buffer!(buffer::SystemBuffer) + (; active_particle) = buffer + + resize!(buffer.eachparticle, count(active_particle)) + + i = 1 + for j in eachindex(active_particle) + if active_particle[j] + buffer.eachparticle[i] = j + i += 1 + end + end + + return buffer +end + +@inline each_moving_particle(system, buffer) = buffer.eachparticle + +@inline active_coordinates(u, system, buffer) = view(u, :, buffer.active_particle) + +@inline active_particles(system, buffer) = buffer.eachparticle + +@inline function activate_next_particle(system) + (; active_particle) = system.buffer + + next_particle = findfirst(x -> !x, active_particle) + + if isnothing(next_particle) + error("0 out of $(system.buffer.buffer_size) buffer particles available") + end + + active_particle[next_particle] = true + + return next_particle +end + +@inline function deactivate_particle!(system, particle, u) + (; active_particle) = system.buffer + + active_particle[particle] = false + + # Set particle far away from simulation domain + for dim in 1:ndims(system) + # Inf or NaN causes instability outcome. + u[dim, particle] = 1e16 + end + + return system +end diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index ee89c92a6..7244253b0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -35,6 +35,19 @@ end return v[end, particle] end +# WARNING! +# These functions are intended to be used internally to set the density +# of newly activated particles in a callback. +# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` +# outside of callbacks. +@inline set_particle_density!(v, system, ::SummationDensity, particle, density) = v + +@inline function set_particle_density!(v, system, ::ContinuityDensity, particle, density) + v[end, particle] = density + + return v +end + function summation_density!(system, semi, u, u_ode, density; particles=each_moving_particle(system)) set_zero!(density) diff --git a/src/general/general.jl b/src/general/general.jl index 4a60cb52f..774d076c9 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -30,6 +30,7 @@ include("density_calculators.jl") include("corrections.jl") include("smoothing_kernels.jl") include("initial_condition.jl") +include("buffer.jl") include("system.jl") include("interpolation.jl") include("custom_quantities.jl") diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 6a1d703b2..73c1fe56d 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -153,7 +153,19 @@ end return compact_support(smoothing_kernel, smoothing_length) end +@inline function compact_support(system::OpenBoundarySPHSystem, neighbor) + # Use the compact support of the fluid + return compact_support(neighbor, system) +end + +@inline function compact_support(system::OpenBoundarySPHSystem, + neighbor::OpenBoundarySPHSystem) + # This NHS is never used + return 0.0 +end + @inline function compact_support(system::BoundaryDEMSystem, neighbor::BoundaryDEMSystem) + # This NHS is never used return 0.0 end @@ -281,6 +293,9 @@ function semidiscretize(semi, tspan; reset_threads=true, data_type=nothing) # Initialize this system initialize!(system, neighborhood_search) + + # Only for systems requiring a mandatory callback + reset_callback_flag!(system) end end @@ -343,6 +358,9 @@ function restart_with!(semi, sol; reset_threads=true) u = wrap_u(sol.u[end].x[2], system, semi) restart_with!(system, v, u) + + # Only for systems requiring a mandatory callback + reset_callback_flag!(system) end return semi @@ -429,7 +447,7 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) end # Update the systems and neighborhood searches (NHS) for a simulation before calling `interact!` to compute forces -function update_systems_and_nhs(v_ode, u_ode, semi, t) +function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=false) # First update step before updating the NHS # (for example for writing the current coordinates in the solid system) foreach_system(semi) do system @@ -466,7 +484,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - update_final!(system, v, u, v_ode, u_ode, semi, t) + update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback) end end @@ -631,6 +649,32 @@ function update_nhs!(neighborhood_search, particles_moving=(true, neighbor.ismoving[])) end +function update_nhs!(neighborhood_search, + system::FluidSystem, neighbor::OpenBoundarySPHSystem, + u_system, u_neighbor) + # The current coordinates of fluids and open boundaries change over time. + + # TODO: Update only `active_coordinates` of open boundaries. + # Problem: Removing inactive particles from neighboring lists is necessary. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + +function update_nhs!(neighborhood_search, + system::OpenBoundarySPHSystem, neighbor::FluidSystem, + u_system, u_neighbor) + # The current coordinates of both open boundaries and fluids change over time. + + # TODO: Update only `active_coordinates` of open boundaries. + # Problem: Removing inactive particles from neighboring lists is necessary. + PointNeighbors.update!(neighborhood_search, + current_coordinates(u_system, system), + current_coordinates(u_neighbor, neighbor), + particles_moving=(true, true)) +end + function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::FluidSystem, u_system, u_neighbor) @@ -725,6 +769,14 @@ function update_nhs!(neighborhood_search, return neighborhood_search end +function update_nhs!(neighborhood_search, + system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, + neighbor::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, + u_system, u_neighbor) + # Don't update. This NHS is never used. + return neighborhood_search +end + function check_configuration(systems) foreach_system(systems) do system check_configuration(system, systems) diff --git a/src/general/system.jl b/src/general/system.jl index a4c7754c4..c3f10caad 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -17,7 +17,16 @@ initialize!(system, neighborhood_search) = system @inline n_moving_particles(system) = nparticles(system) @inline eachparticle(system) = Base.OneTo(nparticles(system)) -@inline each_moving_particle(system) = Base.OneTo(n_moving_particles(system)) + +# Wrapper for systems with `SystemBuffer` +@inline each_moving_particle(system) = each_moving_particle(system, system.buffer) +@inline each_moving_particle(system, ::Nothing) = Base.OneTo(n_moving_particles(system)) + +@inline active_coordinates(u, system) = active_coordinates(u, system, system.buffer) +@inline active_coordinates(u, system, ::Nothing) = current_coordinates(u, system) + +@inline active_particles(system) = active_particles(system, system.buffer) +@inline active_particles(system, ::Nothing) = eachparticle(system) # This should not be dispatched by system type. We always expect to get a column of `A`. @inline function extract_svector(A, system, i) @@ -64,6 +73,10 @@ end return zero(SVector{ndims(system), eltype(system)}) end +@inline set_particle_density!(v, system, particle, density) = v + +@inline set_particle_pressure!(v, system, particle, pressure) = v + @inline function smoothing_kernel(system, distance) (; smoothing_kernel, smoothing_length) = system return kernel(smoothing_kernel, distance, smoothing_length) @@ -103,6 +116,9 @@ function update_pressure!(system, v, u, v_ode, u_ode, semi, t) return system end -function update_final!(system, v, u, v_ode, u_ode, semi, t) +function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false) return system end + +# Only for systems requiring a mandatory callback +reset_callback_flag!(system) = system diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl new file mode 100644 index 000000000..586d565b1 --- /dev/null +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -0,0 +1,340 @@ +@doc raw""" + InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + +Inflow boundary zone for [`OpenBoundarySPHSystem`](@ref). + +The specified plane (line in 2D or rectangle in 3D) will be extruded in upstream +direction (the direction opposite to `flow_direction`) to create a box for the boundary zone. +There are three ways to specify the actual shape of the inflow: +1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then + be filled with inflow particles (default). +2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, + which is then extruded in upstream direction to create the inflow particles. + - In 2D, the shape must be either an initial condition with 2D coordinates, which lies + on the line specified by `plane`, or an initial condition with 1D coordinates, which lies + on the line specified by `plane` when a y-coordinate of `0` is added. + - In 3D, the shape must be either an initial condition with 3D coordinates, which lies + in the rectangle specified by `plane`, or an initial condition with 2D coordinates, + which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. +3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, + which will be used for the inflow particles. + +!!! note "Note" + Particles outside the boundary zone box will be removed. + +# Keywords +- `plane`: Tuple of points defining a part of the surface of the domain. + The points must either span a line in 2D or a rectangle in 3D. + This line or rectangle is then extruded in upstream direction to obtain + the boundary zone. + In 2D, pass two points ``(A, B)``, so that the interval ``[A, B]`` is + the inflow surface. + In 3D, pass three points ``(A, B, C)``, so that the rectangular inflow surface + is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. + These two vectors must be orthogonal. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in upstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for the inflow particles. + Particles outside the boundary zone will be removed. + Do not use together with `extrude_geometry`. +- `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane + and is extruded upstream to obtain the inflow particles. + See point 2 above for more details. + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction=[1.0, 0.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction=[0.0, 0.0, 1.0] + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +inflow = InFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct InFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function InFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, + density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the inflow plane. + # The normal vector must point in upstream direction for an inflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to inflow plane")) + else + # Flip the normal vector to point in the opposite direction of `flow_direction` + spanning_set[:, 1] .*= -sign(dot_) + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove particles outside the boundary zone. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end +end + +@doc raw""" + OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + +Outflow boundary zone for [`OpenBoundarySPHSystem`](@ref). + +The specified plane (line in 2D or rectangle in 3D) will be extruded in downstream +direction (the direction in `flow_direction`) to create a box for the boundary zone. +There are three ways to specify the actual shape of the outflow: +1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then + be filled with outflow particles (default). +2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D, + which is then extruded in downstream direction to create the outflow particles. + - In 2D, the shape must be either an initial condition with 2D coordinates, which lies + on the line specified by `plane`, or an initial condition with 1D coordinates, which lies + on the line specified by `plane` when a y-coordinate of `0` is added. + - In 3D, the shape must be either an initial condition with 3D coordinates, which lies + in the rectangle specified by `plane`, or an initial condition with 2D coordinates, + which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added. +3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D, + which will be used for the outflow particles. + +!!! note "Note" + Particles outside the boundary zone box will be removed. + +# Keywords +- `plane`: Tuple of points defining a part of the surface of the domain. + The points must either span a line in 2D or a rectangle in 3D. + This line or rectangle is then extruded in downstream direction to obtain + the boundary zone. + In 2D, pass two points ``(A, B)``, so that the interval ``[A, B]`` is + the outflow surface. + In 3D, pass three points ``(A, B, C)``, so that the rectangular outflow surface + is spanned by the vectors ``\widehat{AB}`` and ``\widehat{AC}``. + These two vectors must be orthogonal. +- `flow_direction`: Vector defining the flow direction. +- `open_boundary_layers`: Number of particle layers in downstream direction. +- `particle_spacing`: The spacing between the particles (see [`InitialCondition`](@ref)). +- `density`: Particle density (see [`InitialCondition`](@ref)). +- `initial_condition=nothing`: `InitialCondition` for the outflow particles. + Particles outside the boundary zone will be removed. + Do not use together with `extrude_geometry`. +- `extrude_geometry=nothing`: 1D shape in 2D or 2D shape in 3D, which lies on the plane + and is extruded downstream to obtain the outflow particles. + See point 2 above for more details. + +# Examples +```julia +# 2D +plane_points = ([0.0, 0.0], [0.0, 1.0]) +flow_direction = [1.0, 0.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D +plane_points = ([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]) +flow_direction = [0.0, 0.0, 1.0] + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + open_boundary_layers=4) + +# 3D particles sampled as cylinder +circle = SphereShape(0.1, 0.5, (0.5, 0.5), 1.0, sphere_type=RoundSphere()) + +outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, flow_direction, density=1.0, + extrude_geometry=circle, open_boundary_layers=4) +``` + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} + initial_condition :: IC + spanning_set :: S + zone_origin :: ZO + zone_width :: ZW + flow_direction :: FD + + function OutFlow(; plane, flow_direction, density, particle_spacing, + initial_condition=nothing, extrude_geometry=nothing, + open_boundary_layers::Integer) + if open_boundary_layers <= 0 + throw(ArgumentError("`open_boundary_layers` must be positive and greater than zero")) + end + + # Unit vector pointing in downstream direction + flow_direction_ = normalize(SVector(flow_direction...)) + + # Sample particles in boundary zone + if isnothing(initial_condition) && isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(plane; particle_spacing, + density, + direction=flow_direction_, + n_extrude=open_boundary_layers) + elseif !isnothing(extrude_geometry) + initial_condition = TrixiParticles.extrude_geometry(extrude_geometry; + particle_spacing, density, + direction=-flow_direction_, + n_extrude=open_boundary_layers) + end + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + zone_width = open_boundary_layers * initial_condition.particle_spacing + + # Vectors spanning the boundary zone/box + spanning_set, zone_origin = calculate_spanning_vectors(plane, zone_width) + + # First vector of `spanning_vectors` is normal to the outflow plane. + # The normal vector must point in downstream direction for an outflow boundary. + dot_ = dot(normalize(spanning_set[:, 1]), flow_direction_) + + if !isapprox(abs(dot_), 1.0, atol=1e-7) + throw(ArgumentError("`flow_direction` is not normal to outflow plane")) + else + # Flip the normal vector to point in `flow_direction` + spanning_set[:, 1] .*= sign(dot_) + end + + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) + + # Remove particles outside the boundary zone. + # This check is only necessary when `initial_condition` or `extrude_geometry` are passed. + ic = remove_outside_particles(initial_condition, spanning_set_, zone_origin) + + return new{NDIMS, typeof(ic), typeof(spanning_set_), typeof(zone_origin), + typeof(zone_width), + typeof(flow_direction_)}(ic, spanning_set_, zone_origin, zone_width, + flow_direction_) + end +end + +@inline Base.ndims(::Union{InFlow{NDIMS}, OutFlow{NDIMS}}) where {NDIMS} = NDIMS + +function calculate_spanning_vectors(plane, zone_width) + return spanning_vectors(Tuple(plane), zone_width), SVector(plane[1]...) +end + +function spanning_vectors(plane_points::NTuple{2}, zone_width) + plane_size = plane_points[2] - plane_points[1] + + # Calculate normal vector of plane + b = normalize([-plane_size[2], plane_size[1]]) * zone_width + + return hcat(b, plane_size) +end + +function spanning_vectors(plane_points::NTuple{3}, zone_width) + # Vectors spanning the plane + edge1 = plane_points[2] - plane_points[1] + edge2 = plane_points[3] - plane_points[1] + + if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) + throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal")) + end + + # Calculate normal vector of plane + c = Vector(normalize(cross(edge2, edge1)) * zone_width) + + return hcat(c, edge1, edge2) +end + +@inline function is_in_boundary_zone(boundary_zone::Union{InFlow, OutFlow}, particle_coords) + (; zone_origin, spanning_set) = boundary_zone + particle_position = particle_coords - zone_origin + + return is_in_boundary_zone(spanning_set, particle_position) +end + +@inline function is_in_boundary_zone(spanning_set::AbstractArray, + particle_position::SVector{NDIMS}) where {NDIMS} + for dim in 1:NDIMS + span_dim = spanning_set[dim] + # Checks whether the projection of the particle position + # falls within the range of the zone + if !(0 <= dot(particle_position, span_dim) <= dot(span_dim, span_dim)) + + # Particle is not in boundary zone + return false + end + end + + # Particle is in boundary zone + return true +end + +function remove_outside_particles(initial_condition, spanning_set, zone_origin) + (; coordinates, density, particle_spacing) = initial_condition + + in_zone = trues(nparticles(initial_condition)) + + for particle in eachparticle(initial_condition) + current_position = current_coords(coordinates, initial_condition, particle) + particle_position = current_position - zone_origin + + in_zone[particle] = is_in_boundary_zone(spanning_set, particle_position) + end + + return InitialCondition(; coordinates=coordinates[:, in_zone], density=first(density), + particle_spacing) +end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl new file mode 100644 index 000000000..43c0a6ff1 --- /dev/null +++ b/src/schemes/boundary/open_boundary/system.jl @@ -0,0 +1,483 @@ +@doc raw""" + OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed, + fluid_system::FluidSystem, buffer_size::Integer, + reference_velocity=zeros(ndims(boundary_zone)), + reference_pressure=0.0, + reference_density=first(boundary_zone.initial_condition.density)) + +Open boundary system for in- and outflow particles. +These open boundaries use the characteristic variables to propagate the appropriate values +to the outlet or inlet and have been proposed by Lastiwka et al. (2009). For more information +about the method see [description below](@ref method_of_characteristics). + +# Arguments +- `boundary_zone`: Use [`InFlow`](@ref) for an inflow and [`OutFlow`](@ref) for an outflow boundary. + +# Keywords +- `sound_speed`: Speed of sound. +- `fluid_system`: The corresponding fluid system +- `buffer_size`: Number of buffer particles. +- `reference_velocity`: Reference velocity is either a function mapping each particle's coordinates + and time to its velocity, an array where the ``i``-th column holds + the velocity of particle ``i`` or, for a constant fluid velocity, + a vector holding this velocity. Velocity is constant zero by default. +- `reference_pressure`: Reference pressure is either a function mapping each particle's coordinates + and time to its pressure, a vector holding the pressure of each particle, + or a scalar for a constant pressure over all particles. + Pressure is constant zero by default. +- `reference_density`: Reference density is either a function mapping each particle's coordinates + and time to its density, a vector holding the density of each particle, + or a scalar for a constant density over all particles. + Density is the density of the first particle in the initial condition by default. + +!!! warning "Experimental Implementation" + This is an experimental feature and may change in any future releases. +""" +struct OpenBoundarySPHSystem{BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, ARRAY2D, RV, RP, + RD, B} <: System{NDIMS, IC} + initial_condition :: IC + fluid_system :: FS + mass :: ARRAY1D # Array{ELTYPE, 1}: [particle] + density :: ARRAY1D # Array{ELTYPE, 1}: [particle] + volume :: ARRAY1D # Array{ELTYPE, 1}: [particle] + pressure :: ARRAY1D # Array{ELTYPE, 1}: [particle] + characteristics :: ARRAY2D # Array{ELTYPE, 2}: [characteristic, particle] + previous_characteristics :: ARRAY2D # Array{ELTYPE, 2}: [characteristic, particle] + sound_speed :: ELTYPE + boundary_zone :: BZ + flow_direction :: SVector{NDIMS, ELTYPE} + reference_velocity :: RV + reference_pressure :: RP + reference_density :: RD + buffer :: B + update_callback_used :: Ref{Bool} + + function OpenBoundarySPHSystem(boundary_zone::Union{InFlow, OutFlow}; sound_speed, + fluid_system::FluidSystem, buffer_size::Integer, + reference_velocity=zeros(ndims(boundary_zone)), + reference_pressure=0.0, + reference_density=first(boundary_zone.initial_condition.density)) + (; initial_condition) = boundary_zone + + buffer = SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + + if !(reference_velocity isa Function || + (reference_velocity isa Vector && length(reference_velocity) == NDIMS)) + throw(ArgumentError("`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "an array where the ``i``-th column holds the velocity of particle ``i`` " * + "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) + else + reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) + end + + if !(reference_pressure isa Function || reference_pressure isa Real) + throw(ArgumentError("`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "a vector holding the pressure of each particle, or a scalar")) + else + reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) + end + + if !(reference_density isa Function || reference_density isa Real) + throw(ArgumentError("`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "a vector holding the density of each particle, or a scalar")) + else + reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) + end + + mass = copy(initial_condition.mass) + pressure = [reference_pressure_(initial_condition.coordinates[:, i], 0.0) + for i in eachparticle(initial_condition)] + density = copy(initial_condition.density) + volume = similar(initial_condition.density) + + characteristics = zeros(ELTYPE, 3, length(mass)) + previous_characteristics = zeros(ELTYPE, 3, length(mass)) + + flow_direction_ = boundary_zone.flow_direction + + return new{typeof(boundary_zone), NDIMS, ELTYPE, typeof(initial_condition), + typeof(fluid_system), typeof(mass), typeof(characteristics), + typeof(reference_velocity_), typeof(reference_pressure_), + typeof(reference_density_), + typeof(buffer)}(initial_condition, fluid_system, mass, density, volume, + pressure, characteristics, previous_characteristics, + sound_speed, boundary_zone, flow_direction_, + reference_velocity_, reference_pressure_, + reference_density_, buffer, false) + end +end + +timer_name(::OpenBoundarySPHSystem) = "open_boundary" +vtkname(system::OpenBoundarySPHSystem) = "open_boundary" + +function Base.show(io::IO, system::OpenBoundarySPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "OpenBoundarySPHSystem{", ndims(system), "}(") + print(io, type2string(system.boundary_zone)) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "OpenBoundarySPHSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + summary_line(io, "fluid system", type2string(system.fluid_system)) + summary_line(io, "boundary", type2string(system.boundary_zone)) + summary_line(io, "flow direction", system.flow_direction) + summary_line(io, "prescribed velocity", string(nameof(system.reference_velocity))) + summary_line(io, "prescribed pressure", string(nameof(system.reference_pressure))) + summary_line(io, "prescribed density", string(nameof(system.reference_density))) + summary_line(io, "width", round(system.boundary_zone.zone_width, digits=3)) + summary_footer(io) + end +end + +function reset_callback_flag(system::OpenBoundarySPHSystem) + system.update_callback_used[] = false + + return system +end + +update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true + +@inline source_terms(system::OpenBoundarySPHSystem) = nothing + +@inline hydrodynamic_mass(system::OpenBoundarySPHSystem, particle) = system.mass[particle] + +@inline function particle_density(v, system::OpenBoundarySPHSystem, particle) + return system.density[particle] +end + +@inline function particle_pressure(v, system::OpenBoundarySPHSystem, particle) + return system.pressure[particle] +end + +@inline function update_quantities!(system::OpenBoundarySPHSystem, v, u, t) + (; density, pressure, characteristics, flow_direction, sound_speed, + reference_velocity, reference_pressure, reference_density) = system + + # Update quantities based on the characteristic variables + @threaded for particle in each_moving_particle(system) + particle_position = current_coords(u, system, particle) + + J1 = characteristics[1, particle] + J2 = characteristics[2, particle] + J3 = characteristics[3, particle] + + rho_ref = reference_density(particle_position, t) + density[particle] = rho_ref + ((-J1 + 0.5 * (J2 + J3)) / sound_speed^2) + + p_ref = reference_pressure(particle_position, t) + pressure[particle] = p_ref + 0.5 * (J2 + J3) + + v_ref = reference_velocity(particle_position, t) + rho = density[particle] + v_ = v_ref + ((J2 - J3) / (2 * sound_speed * rho)) * flow_direction + + for dim in 1:ndims(system) + v[dim, particle] = v_[dim] + end + end + + return system +end + +function update_final!(system::OpenBoundarySPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) + if !update_from_callback && !(system.update_callback_used[]) + throw(ArgumentError("`UpdateCallback` is required when using `OpenBoundarySPHSystem`")) + end + + @trixi_timeit timer() "evaluate characteristics" evaluate_characteristics!(system, v, u, + v_ode, u_ode, + semi, t) +end + +# ==== Characteristics +# J1: Associated with convection and entropy and propagates at flow velocity. +# J2: Propagates downstream to the local flow +# J3: Propagates upstream to the local flow +function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) + (; volume, characteristics, previous_characteristics, boundary_zone) = system + + for particle in eachparticle(system) + previous_characteristics[1, particle] = characteristics[1, particle] + previous_characteristics[2, particle] = characteristics[2, particle] + previous_characteristics[3, particle] = characteristics[3, particle] + end + + set_zero!(characteristics) + set_zero!(volume) + + # Evaluate the characteristic variables with the fluid system + evaluate_characteristics!(system, system.fluid_system, v, u, v_ode, u_ode, semi, t) + + # Only some of the in-/outlet particles are in the influence of the fluid particles. + # Thus, we compute the characteristics for the particles that are outside the influence + # of fluid particles by using the average of the values of the previous time step. + # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119 + @threaded for particle in each_moving_particle(system) + + # Particle is outside of the influence of fluid particles + if isapprox(volume[particle], 0.0) + + # Using the average of the values at the previous time step for particles which + # are outside of the influence of fluid particles. + avg_J1 = 0.0 + avg_J2 = 0.0 + avg_J3 = 0.0 + counter = 0 + + for neighbor in each_moving_particle(system) + # Make sure that only neighbors in the influence of + # the fluid particles are used. + if volume[neighbor] > sqrt(eps()) + avg_J1 += previous_characteristics[1, neighbor] + avg_J2 += previous_characteristics[2, neighbor] + avg_J3 += previous_characteristics[3, neighbor] + counter += 1 + end + end + + characteristics[1, particle] = avg_J1 / counter + characteristics[2, particle] = avg_J2 / counter + characteristics[3, particle] = avg_J3 / counter + else + characteristics[1, particle] /= volume[particle] + characteristics[2, particle] /= volume[particle] + characteristics[3, particle] /= volume[particle] + end + prescribe_conditions!(characteristics, particle, boundary_zone) + end + + return system +end + +function evaluate_characteristics!(system, neighbor_system::FluidSystem, + v, u, v_ode, u_ode, semi, t) + (; volume, sound_speed, characteristics, flow_direction, + reference_velocity, reference_pressure, reference_density) = system + + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + + nhs = get_neighborhood_search(system, neighbor_system, semi) + + system_coords = current_coordinates(u, system) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all fluid neighbors within the kernel cutoff + for_particle_neighbor(system, neighbor_system, system_coords, neighbor_coords, + nhs) do particle, neighbor, pos_diff, distance + neighbor_position = current_coords(u_neighbor_system, neighbor_system, neighbor) + + # Determine current and prescribed quantities + rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_ref = reference_density(neighbor_position, t) + + p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) + p_ref = reference_pressure(neighbor_position, t) + + v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + v_neighbor_ref = reference_velocity(neighbor_position, t) + + # Determine characteristic variables + density_term = -sound_speed^2 * (rho_b - rho_ref) + pressure_term = p_b - p_ref + velocity_term = rho_b * sound_speed * (dot(v_b - v_neighbor_ref, flow_direction)) + + kernel_ = smoothing_kernel(neighbor_system, distance) + + characteristics[1, particle] += (density_term + pressure_term) * kernel_ + characteristics[2, particle] += (velocity_term + pressure_term) * kernel_ + characteristics[3, particle] += (-velocity_term + pressure_term) * kernel_ + + volume[particle] += kernel_ + end + + return system +end + +@inline function prescribe_conditions!(characteristics, particle, ::OutFlow) + # J3 is prescribed (i.e. determined from the exterior of the domain). + # J1 and J2 is transimtted from the domain interior. + characteristics[3, particle] = zero(eltype(characteristics)) + + return characteristics +end + +@inline function prescribe_conditions!(characteristics, particle, ::InFlow) + # Allow only J3 to propagate upstream to the boundary + characteristics[1, particle] = zero(eltype(characteristics)) + characteristics[2, particle] = zero(eltype(characteristics)) + + return characteristics +end + +# This function is called by the `UpdateCallback`, as the integrator array might be modified +function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ode, + semi, t) + u = wrap_u(u_ode, system, semi) + v = wrap_v(v_ode, system, 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 quantities" update_quantities!(system, v, u, t) + + @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) + + # Update buffers + update_system_buffer!(system.buffer) + update_system_buffer!(system.fluid_system.buffer) +end + +update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t) = system + +function check_domain!(system, v, u, v_ode, u_ode, semi) + (; boundary_zone, fluid_system) = system + + u_fluid = wrap_u(u_ode, fluid_system, semi) + v_fluid = wrap_v(v_ode, fluid_system, semi) + + neighborhood_search = get_neighborhood_search(system, fluid_system, semi) + + for particle in each_moving_particle(system) + particle_coords = current_coords(u, system, particle) + + # Check if boundary particle is outside the boundary zone + if !is_in_boundary_zone(boundary_zone, particle_coords) + convert_particle!(system, fluid_system, boundary_zone, particle, + v, u, v_fluid, u_fluid) + end + + # Check the neighboring fluid particles whether they're entering the boundary zone + for neighbor in PointNeighbors.eachneighbor(particle_coords, neighborhood_search) + fluid_coords = current_coords(u_fluid, fluid_system, neighbor) + + # Check if neighboring fluid particle is in boundary zone + if is_in_boundary_zone(boundary_zone, fluid_coords) + convert_particle!(fluid_system, system, boundary_zone, neighbor, + v, u, v_fluid, u_fluid) + end + end + end + + return system +end + +# Outflow particle is outside the boundary zone +@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, + boundary_zone::OutFlow, particle, v, u, + v_fluid, u_fluid) + deactivate_particle!(system, particle, u) + + return system +end + +# Inflow particle is outside the boundary zone +@inline function convert_particle!(system::OpenBoundarySPHSystem, fluid_system, + boundary_zone::InFlow, particle, v, u, + v_fluid, u_fluid) + (; spanning_set) = boundary_zone + + # Activate a new particle in simulation domain + transfer_particle!(fluid_system, system, particle, v_fluid, u_fluid, v, u) + + # Reset position of boundary particle + for dim in 1:ndims(system) + u[dim, particle] += spanning_set[1][dim] + end + + return system +end + +# Fluid particle is in boundary zone +@inline function convert_particle!(fluid_system::FluidSystem, system, + boundary_zone, particle, v, u, v_fluid, u_fluid) + # Activate particle in boundary zone + transfer_particle!(system, fluid_system, particle, v, u, v_fluid, u_fluid) + + # Deactivate particle in interior domain + deactivate_particle!(fluid_system, particle, u_fluid) + + return fluid_system +end + +@inline function transfer_particle!(system_new, system_old, particle_old, + v_new, u_new, v_old, u_old) + particle_new = activate_next_particle(system_new) + + # Transfer densities + density = particle_density(v_old, system_old, particle_old) + set_particle_density!(v_new, system_new, particle_new, density) + + # Transfer pressure + pressure = particle_pressure(v_old, system_old, particle_old) + set_particle_pressure!(v_new, system_new, particle_new, pressure) + + # Exchange position and velocity + for dim in 1:ndims(system_new) + u_new[dim, particle_new] = u_old[dim, particle_old] + v_new[dim, particle_new] = v_old[dim, particle_old] + end + + # TODO: Only when using TVF: set tvf + + return system_new +end + +function write_v0!(v0, system::OpenBoundarySPHSystem) + (; initial_condition) = system + + for particle in eachparticle(system) + # Write particle velocities + for dim in 1:ndims(system) + v0[dim, particle] = initial_condition.velocity[dim, particle] + end + end + + return v0 +end + +function write_u0!(u0, system::OpenBoundarySPHSystem) + (; initial_condition) = system + + for particle in eachparticle(system) + # Write particle velocities + for dim in 1:ndims(system) + u0[dim, particle] = initial_condition.coordinates[dim, particle] + end + end + + return u0 +end + +function wrap_reference_function(function_::Function, ::Val) + # Already a function + return function_ +end + +# Name the function so that the summary box does know which kind of function this is +function wrap_reference_function(constant_scalar_::Number, ::Val) + return constant_scalar(coords, t) = constant_scalar_ +end + +# For vectors and tuples +# Name the function so that the summary box does know which kind of function this is +function wrap_reference_function(constant_vector_, ::Val{NDIMS}) where {NDIMS} + return constant_vector(coords, t) = SVector{NDIMS}(constant_vector_) +end diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 201a3bbb9..d9cfb1779 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -1,7 +1,7 @@ -# Interaction of boundary with other systems +# Interaction of boundary with other systems function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, neighborhood_search, - particle_system::BoundarySystem, + particle_system::Union{BoundarySystem, OpenBoundarySPHSystem}, neighbor_system) # TODO Solids and moving boundaries should be considered in the continuity equation return dv diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index cade5fbe9..73af6a843 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -13,7 +13,8 @@ The interaction between fluid and boundary particles is specified by the boundar - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface. Note: currently it is assumed that all fluids have the same adhesion coefficient. """ -struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem{NDIMS, IC} +struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, + CA} <: BoundarySystem{NDIMS, IC} initial_condition :: IC coordinates :: CO # Array{ELTYPE, 2} boundary_model :: BM @@ -21,6 +22,7 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem ismoving :: IM # Ref{Bool} (to make a mutable field compatible with GPUs) adhesion_coefficient :: ELTYPE cache :: CA + buffer :: Nothing # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. @@ -28,12 +30,10 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE, IC, CO, M, IM, CA} <: BoundarySystem ismoving, adhesion_coefficient, cache) ELTYPE = eltype(coordinates) - new{typeof(boundary_model), size(coordinates, 1), ELTYPE, - typeof(initial_condition), typeof(coordinates), - typeof(movement), typeof(ismoving), typeof(cache)}(initial_condition, - coordinates, boundary_model, - movement, ismoving, - adhesion_coefficient, cache) + new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), + typeof(coordinates), typeof(movement), typeof(ismoving), + typeof(cache)}(initial_condition, coordinates, boundary_model, movement, + ismoving, adhesion_coefficient, cache, nothing) end end @@ -73,6 +73,7 @@ struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC, coordinates :: ARRAY2D # [dimension, particle] radius :: ARRAY1D # [particle] normal_stiffness :: ELTYPE + buffer :: Nothing function BoundaryDEMSystem(initial_condition, normal_stiffness) coordinates = initial_condition.coordinates @@ -80,9 +81,9 @@ struct BoundaryDEMSystem{NDIMS, ELTYPE <: Real, IC, ones(length(initial_condition.mass)) NDIMS = size(coordinates, 1) - return new{NDIMS, eltype(coordinates), typeof(initial_condition), - typeof(radius), typeof(coordinates)}(initial_condition, coordinates, - radius, normal_stiffness) + return new{NDIMS, eltype(coordinates), typeof(initial_condition), typeof(radius), + typeof(coordinates)}(initial_condition, coordinates, radius, + normal_stiffness, nothing) end end @@ -335,7 +336,8 @@ end # This update depends on the computed quantities of the fluid system and therefore # has to be in `update_final!` after `update_quantities!`. -function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t) +function update_final!(system::BoundarySPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) (; boundary_model) = system update_pressure!(boundary_model, system, v, u, v_ode, u_ode, semi) diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index f59bca7d3..4093e8ccb 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -4,7 +4,7 @@ pressure_acceleration=inter_particle_averaged_pressure, density_calculator=SummationDensity(), alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, NDIMS), + acceleration=ntuple(_ -> 0.0, NDIMS), buffer_size=nothing, source_terms=nothing) System for particles of a fluid. @@ -28,6 +28,8 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more When set to `nothing`, the pressure acceleration formulation for the corresponding [density calculator](@ref density_calculator) is chosen. - `density_calculator`: [Density calculator](@ref density_calculator) (default: [`SummationDensity`](@ref)) +- `buffer_size`: Number of buffer particles. + This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure)` (which are the quantities of a single particle), returning a `Tuple` @@ -40,7 +42,7 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, - PF, ST, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC @@ -53,6 +55,7 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, correction :: Nothing pressure_acceleration_formulation :: PF source_terms :: ST + buffer :: B cache :: C function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, @@ -62,7 +65,12 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - source_terms=nothing) + source_terms=nothing, buffer_size=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -89,9 +97,11 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(pressure_acceleration), typeof(source_terms), + typeof(buffer), typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, source_terms, cache) + smoothing_length, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, source_terms, + buffer, cache) end end @@ -113,7 +123,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst show(io, system) else summary_header(io, "EntropicallyDampedSPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) + if system.buffer isa SystemBuffer + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + else + summary_line(io, "#particles", nparticles(system)) + end summary_line(io, "density calculator", system.density_calculator |> typeof |> nameof) summary_line(io, "viscosity", system.viscosity |> typeof |> nameof) @@ -145,6 +160,18 @@ end return v[end, particle] end +# WARNING! +# These functions are intended to be used internally to set the pressure +# of newly activated particles in a callback. +# DO NOT use outside a callback. OrdinaryDiffEq does not allow changing `v` and `u` +# outside of callbacks. +@inline function set_particle_pressure!(v, system::EntropicallyDampedSPHSystem, particle, + pressure) + v[end, particle] = pressure + + return v +end + @inline system_sound_speed(system::EntropicallyDampedSPHSystem) = system.sound_speed function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index 056662bbc..68391b285 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -1,3 +1,7 @@ +@inline function set_particle_density!(v, system::FluidSystem, particle, density) + set_particle_density!(v, system, system.density_calculator, particle, density) +end + function create_cache_density(initial_condition, ::SummationDensity) density = similar(initial_condition.density) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 554f41c14..a678b40d4 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -11,6 +11,19 @@ function dv_viscosity(particle_system, neighbor_system, sound_speed, m_a, m_b, rho_mean) end +function dv_viscosity(particle_system, neighbor_system::OpenBoundarySPHSystem, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) + # No viscosity in the open boundary system. Use viscosity of the fluid system. + viscosity = viscosity_model(particle_system) + + return dv_viscosity(viscosity, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_mean) +end + function dv_viscosity(viscosity, particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index cb9a1da87..a8e28670b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -4,6 +4,7 @@ smoothing_kernel, smoothing_length; viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), + buffer_size=nothing, correction=nothing, source_terms=nothing) System for particles of a fluid. @@ -26,6 +27,8 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. See [`ArtificialViscosityMonaghan`](@ref) or [`ViscosityAdami`](@ref). - `density_diffusion`: Density diffusion terms for this system. See [`DensityDiffusion`](@ref). - `acceleration`: Acceleration vector for the system. (default: zero vector) +- `buffer_size`: Number of buffer particles. + This is needed when simulating with [`OpenBoundarySPHSystem`](@ref). - `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure)` @@ -42,7 +45,7 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, SRFT, C} <: FluidSystem{NDIMS, IC} + V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} @@ -57,6 +60,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, pressure_acceleration_formulation :: PF source_terms :: ST surface_tension :: SRFT + buffer :: B cache :: C end @@ -66,11 +70,17 @@ function WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, smoothing_length; pressure_acceleration=nothing, + buffer_size=nothing, viscosity=nothing, density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, surface_tension=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) + + initial_condition = allocate_buffer(initial_condition, buffer) + NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -108,10 +118,11 @@ function WeaklyCompressibleSPHSystem(initial_condition, return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, acceleration_, - viscosity, density_diffusion, correction, - pressure_acceleration, source_terms, surface_tension, - cache) + smoothing_kernel, smoothing_length, + acceleration_, viscosity, + density_diffusion, correction, + pressure_acceleration, + source_terms, surface_tension, buffer, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -166,7 +177,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst show(io, system) else summary_header(io, "WeaklyCompressibleSPHSystem{$(ndims(system))}") - summary_line(io, "#particles", nparticles(system)) + if system.buffer isa SystemBuffer + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "#buffer_particles", system.buffer.buffer_size) + else + summary_line(io, "#particles", nparticles(system)) + end summary_line(io, "density calculator", system.density_calculator |> typeof |> nameof) summary_line(io, "correction method", diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index b81a7768b..fadcb240e 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -1,5 +1,8 @@ # Include all schemes without rhs first. The rhs depends on the systems to define # interactions between the different system types. +# Viscosity requires the open boundary system. +include("boundary/open_boundary/boundary_zones.jl") +include("boundary/open_boundary/system.jl") include("fluid/fluid.jl") include("boundary/boundary.jl") include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") diff --git a/src/schemes/solid/discrete_element_method/system.jl b/src/schemes/solid/discrete_element_method/system.jl index 5aadfcfcc..eb30dfb34 100644 --- a/src/schemes/solid/discrete_element_method/system.jl +++ b/src/schemes/solid/discrete_element_method/system.jl @@ -36,6 +36,7 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, I damping_coefficient :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} source_terms :: ST + buffer :: Nothing function DEMSystem(initial_condition, normal_stiffness, elastic_modulus, poissons_ratio; damping_coefficient=0.0001, @@ -56,7 +57,8 @@ struct DEMSystem{NDIMS, ELTYPE <: Real, IC, ARRAY1D, ST} <: SolidSystem{NDIMS, I return new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(source_terms)}(initial_condition, mass, radius, elastic_modulus, poissons_ratio, normal_stiffness, - damping_coefficient, acceleration_, source_terms) + damping_coefficient, acceleration_, source_terms, + nothing) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 5b3021b1c..9203425fb 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -69,6 +69,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, boundary_model :: BM penalty_force :: PF source_terms :: ST + buffer :: Nothing function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, @@ -116,7 +117,7 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, n_moving_particles, young_modulus, poisson_ratio, lame_lambda, lame_mu, smoothing_kernel, smoothing_length, acceleration_, boundary_model, - penalty_force, source_terms) + penalty_force, source_terms, nothing) end end @@ -246,7 +247,8 @@ function update_quantities!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode return system end -function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t) +function update_final!(system::TotalLagrangianSPHSystem, v, u, v_ode, u_ode, semi, t; + update_from_callback=false) (; boundary_model) = system # Only update boundary model diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index 9b7a62401..82df9aaa3 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -8,8 +8,10 @@ Extrude either a line, a plane or a shape along a specific direction. # Arguments - `geometry`: Either particle coordinates or an [`InitialCondition`](@ref) defining a 2D shape to extrude to a 3D volume, or two 2D points - defining a line to extrude to a plane in 2D, or three 3D points defining - a parallelogram to extrude to a parallelepiped. + ``(A, B)`` defining the interval ``[A, B]`` to extrude to a plane + in 2D, or three 3D points ``(A, B, C)`` defining the parallelogram + spanned by the vectors ``\widehat{AB}`` and ``\widehat {AC}`` to extrude + to a parallelepiped. # Keywords - `particle_spacing`: Spacing between the particles. Can be omitted when `geometry` is an diff --git a/src/visualization/recipes_plots.jl b/src/visualization/recipes_plots.jl index a3fb7de7e..f7e4d5b23 100644 --- a/src/visualization/recipes_plots.jl +++ b/src/visualization/recipes_plots.jl @@ -12,10 +12,8 @@ end RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; size=(600, 400)) # Default size systems_data = map(semi.systems) do system - (; initial_condition) = system - u = wrap_u(u_ode, system, semi) - coordinates = current_coordinates(u, system) + coordinates = active_coordinates(u, system) x = collect(coordinates[1, :]) y = collect(coordinates[2, :]) @@ -24,8 +22,8 @@ RecipesBase.@recipe function f(v_ode, u_ode, semi::Semidiscretization; particle_spacing = 0.0 end - x_min, y_min = minimum(initial_condition.coordinates, dims=2) .- 0.5particle_spacing - x_max, y_max = maximum(initial_condition.coordinates, dims=2) .+ 0.5particle_spacing + x_min, y_min = minimum(coordinates, dims=2) .- 0.5particle_spacing + x_max, y_max = maximum(coordinates, dims=2) .+ 0.5particle_spacing return (; x, y, x_min, x_max, y_min, y_max, particle_spacing, label=timer_name(system)) diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index 2a4c9f207..0372e7167 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -56,7 +56,8 @@ function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix # Update quantities that are stored in the systems. These quantities (e.g. pressure) # still have the values from the last stage of the previous step if not updated here. - @trixi_timeit timer() "update systems" update_systems_and_nhs(v_ode, u_ode, semi, t) + @trixi_timeit timer() "update systems" update_systems_and_nhs(v_ode, u_ode, semi, t; + update_from_callback=true) filenames = system_names(systems) @@ -94,7 +95,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix # Reset the collection when the iteration is 0 pvd = paraview_collection(collection_file; append=iter > 0) - points = PointNeighbors.periodic_coords(current_coordinates(u, system), + points = PointNeighbors.periodic_coords(active_coordinates(u, system), periodic_box) cells = [MeshCell(VTKCellTypes.VTK_VERTEX, (i,)) for i in axes(points, 2)] @@ -112,7 +113,7 @@ function trixi2vtk(v, u, t, system, periodic_box; output_directory="out", prefix write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data) # Store particle index - vtk["index"] = eachparticle(system) + vtk["index"] = active_particles(system) vtk["time"] = t if write_meta_data @@ -216,11 +217,12 @@ function write2vtk!(vtk, v, u, t, system; write_meta_data=true) end function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) - vtk["velocity"] = view(v, 1:ndims(system), :) + vtk["velocity"] = [current_velocity(v, system, particle) + for particle in active_particles(system)] vtk["density"] = [particle_density(v, system, particle) - for particle in eachparticle(system)] + for particle in active_particles(system)] vtk["pressure"] = [particle_pressure(v, system, particle) - for particle in eachparticle(system)] + for particle in active_particles(system)] if write_meta_data vtk["acceleration"] = system.acceleration @@ -295,6 +297,37 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end +function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true) + (; reference_velocity, reference_pressure, reference_density) = system + + vtk["velocity"] = [current_velocity(v, system, particle) + for particle in active_particles(system)] + vtk["density"] = [particle_density(v, system, particle) + for particle in active_particles(system)] + vtk["pressure"] = [particle_pressure(v, system, particle) + for particle in active_particles(system)] + + NDIMS = ndims(system) + ELTYPE = eltype(system) + coords = reinterpret(reshape, SVector{NDIMS, ELTYPE}, active_coordinates(u, system)) + + vtk["prescribed_velocity"] = stack(reference_velocity.(coords, t)) + vtk["prescribed_density"] = reference_density.(coords, t) + vtk["prescribed_pressure"] = reference_pressure.(coords, t) + + if write_meta_data + vtk["boundary_zone"] = type2string(system.boundary_zone) + vtk["sound_speed"] = system.sound_speed + vtk["width"] = round(system.boundary_zone.zone_width, digits=3) + vtk["flow_direction"] = system.flow_direction + vtk["velocity_function"] = string(nameof(system.reference_velocity)) + vtk["pressure_function"] = string(nameof(system.reference_pressure)) + vtk["density_function"] = string(nameof(system.reference_density)) + end + + return vtk +end + function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true) write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data) end diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 0f0032607..dbe03a130 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -114,6 +114,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/pipe_flow_2d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "pipe_flow_2d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/dam_break_2d_surface_tension.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/general/buffer.jl b/test/general/buffer.jl new file mode 100644 index 000000000..7a8371ea7 --- /dev/null +++ b/test/general/buffer.jl @@ -0,0 +1,68 @@ +@testset verbose=true "`SystemBuffer`" begin + # Mock fluid system + struct FluidSystemMock3 <: TrixiParticles.FluidSystem{2, Nothing} end + + inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.2, + open_boundary_layers=2, density=1.0, flow_direction=[1.0, 0.0]) + system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, fluid_system=FluidSystemMock3(), + buffer_size=0) + system_buffer = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=5, + fluid_system=FluidSystemMock3()) + + n_particles = nparticles(system) + + @testset "Iterators" begin + @test TrixiParticles.each_moving_particle(system) == 1:n_particles + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:n_particles + + particle_id = TrixiParticles.activate_next_particle(system_buffer) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:(n_particles + 1) + + TrixiParticles.deactivate_particle!(system_buffer, particle_id, + ones(2, particle_id)) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == 1:n_particles + + particle_id = 5 + TrixiParticles.deactivate_particle!(system_buffer, particle_id, + ones(2, particle_id)) + + TrixiParticles.update_system_buffer!(system_buffer.buffer) + + @test TrixiParticles.each_moving_particle(system_buffer) == + setdiff(1:n_particles, particle_id) + end + + @testset "Allocate Buffer" begin + initial_condition = rectangular_patch(0.1, (3, 3), perturbation_factor=0.0) + buffer = TrixiParticles.SystemBuffer(nparticles(initial_condition), 7) + + ic_with_buffer = TrixiParticles.allocate_buffer(initial_condition, buffer) + + @test nparticles(ic_with_buffer) == nparticles(initial_condition) + 7 + + masses = initial_condition.mass[1] .* ones(nparticles(ic_with_buffer)) + @test masses == ic_with_buffer.mass + + densities = initial_condition.density[1] .* ones(nparticles(ic_with_buffer)) + @test densities == ic_with_buffer.density + + pressures = initial_condition.pressure[1] .* ones(nparticles(ic_with_buffer)) + @test pressures == ic_with_buffer.pressure + + @testset "Illegal Input" begin + # The rectangular patch has a perturbed, non-constant density + ic = rectangular_patch(0.1, (3, 3)) + buffer = TrixiParticles.SystemBuffer(9, 7) + + error_str = "`initial_condition.density` needs to be constant when using `SystemBuffer`" + @test_throws ArgumentError(error_str) TrixiParticles.allocate_buffer(ic, buffer) + end + end +end diff --git a/test/general/general.jl b/test/general/general.jl index b59cff789..27f60ebba 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -3,3 +3,4 @@ include("smoothing_kernels.jl") include("density_calculator.jl") include("semidiscretization.jl") include("interpolation.jl") +include("buffer.jl") diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 633a0e240..4aff1f4fd 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -131,6 +131,9 @@ v2 = zeros(4 * 3) v_ode = vcat(vec(v1), v2) + # Avoid `SystemBuffer` barrier + TrixiParticles.each_moving_particle(system::Union{System1, System2}) = TrixiParticles.eachparticle(system) + TrixiParticles.add_source_terms!(dv_ode, v_ode, u_ode, semi) dv1 = TrixiParticles.wrap_v(dv_ode, system1, semi) diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl new file mode 100644 index 000000000..7a34a4bca --- /dev/null +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -0,0 +1,221 @@ +@testset verbose=true "Boundary Zone" begin + @testset verbose=true "Boundary Zone 2D" begin + particle_spacing = 0.2 + open_boundary_layers = 4 + + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset verbose=true "Points $i" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + + plane_size = point_2 - point_1 + + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test isapprox(sign_ * flow_directions[j], + normalize(boundary_zone.spanning_set[1]), atol=1e-14) + @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]), + atol=1e-14) + end + end + end + end + + @testset verbose=true "Boundary Zone 3D" begin + particle_spacing = 0.05 + open_boundary_layers = 4 + + plane_points_1 = [ + [0.0, 0.0, 0.0], + [0.3113730847835541, 0.19079485535621643, -0.440864622592926], + ] + plane_points_2 = [ + [1.0, 0.0, 0.0], + [-0.10468611121177673, 0.252103328704834, -0.44965094327926636], + ] + plane_points_3 = [ + [0.0, 1.0, 0.0], + [0.3113730847835541, 0.25057315826416016, -0.02374829351902008], + ] + + @testset verbose=true "Points $i" for i in eachindex(plane_points_1) + point_1 = plane_points_1[i] + point_2 = plane_points_2[i] + point_3 = plane_points_3[i] + + edge1 = point_2 - point_1 + edge2 = point_3 - point_1 + + flow_directions = [ + normalize(cross(edge1, edge2)), + -normalize(cross(edge1, edge2)), + ] + + @testset verbose=true "Flow Direction $j" for j in eachindex(flow_directions) + inflow = InFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + outflow = OutFlow(; plane=(point_1, point_2, point_3), particle_spacing, + flow_direction=flow_directions[j], density=1.0, + open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + zone_width = open_boundary_layers * + boundary_zone.initial_condition.particle_spacing + sign_ = (boundary_zone isa InFlow) ? -1 : 1 + + @test plane_points_1[i] == boundary_zone.zone_origin + @test plane_points_2[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[2] + @test plane_points_3[i] - boundary_zone.zone_origin == + boundary_zone.spanning_set[3] + @test isapprox(sign_ * flow_directions[j], + normalize(boundary_zone.spanning_set[1]), atol=1e-14) + @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]), + atol=1e-14) + end + end + end + end + + @testset verbose=true "Particle In Boundary Zone 2D" begin + plane_points = [[-0.2, -0.5], [0.3, 0.6]] + plane_size = plane_points[2] - plane_points[1] + + flow_direction = normalize([-plane_size[2], plane_size[1]]) + + inflow = InFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=plane_points, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? sqrt(eps()) : -sqrt(eps()) + + point1 = plane_points[1] + point2 = plane_points[2] + point3 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin + + query_points = Dict( + "Behind" => ([-1.0, -1.0], false), + "Before" => ([2.0, 2.0], false), + "Closely On Point 1" => (point1 + perturb_ * flow_direction, false), + "Closely On Point 2" => (point2 + perturb_ * flow_direction, false), + "Closely On Point 3" => (point3 - perturb_ * flow_direction, false)) + + @testset verbose=true "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == + TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position) + end + end + end + + @testset verbose=true "Particle In Boundary Zone 3D" begin + point1 = [-0.2, -0.5, 0.0] + point2 = [0.3, 0.5, 0.0] + point3 = [0.13111173850909402, -0.665555869254547, 0.0] + + flow_direction = normalize(cross(point2 - point1, point3 - point1)) + + inflow = InFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + outflow = OutFlow(; plane=[point1, point2, point3], particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=4) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones + perturb_ = boundary_zone isa InFlow ? eps() : -eps() + point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin + + query_points = Dict( + "Behind" => ([-1.0, -1.0, 1.2], false), + "Before" => ([2.0, 2.0, -1.2], false), + "Closely On Point 1" => (point1 + perturb_ * flow_direction, false), + "Closely On Point 2" => (point2 + perturb_ * flow_direction, false), + "Closely On Point 3" => (point3 + perturb_ * flow_direction, false), + "Closely On Point 4" => (point4 - perturb_ * flow_direction, false)) + + @testset verbose=true "$k" for k in keys(query_points) + (particle_position, evaluation) = query_points[k] + + @test evaluation == + TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position) + end + end + end + + @testset verbose=true "Illegal Inputs" begin + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + flow_direction = [0.0, 0.0, 1.0] + + error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal" + + @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + @test_throws ArgumentError(error_str) OutFlow(; plane=no_rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + rectangular_plane = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] + flow_direction = [0.0, 1.0, 0.0] + + error_str = "`flow_direction` is not normal to inflow plane" + + @test_throws ArgumentError(error_str) InFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + + error_str = "`flow_direction` is not normal to outflow plane" + + @test_throws ArgumentError(error_str) OutFlow(; plane=rectangular_plane, + particle_spacing=0.1, + flow_direction, density=1.0, + open_boundary_layers=2) + end +end diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl new file mode 100644 index 000000000..1fe20722a --- /dev/null +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -0,0 +1,134 @@ +@testset verbose=true "Characteristic Variables" begin + particle_spacing = 0.1 + + # Number of boundary particles in the influence of fluid particles + influenced_particles = [20, 52, 26] + + open_boundary_layers = 8 + sound_speed = 20.0 + density = 1000.0 + pressure = 5.0 + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 1.2particle_spacing + + # Prescribed quantities + reference_velocity = (pos, t) -> SVector(t, 0.0) + reference_pressure = (pos, t) -> 50_000.0 * t + reference_density = (pos, t) -> 1000.0 * t + + # Plane points of open boundary + plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] + plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] + + @testset "Points $i" for i in eachindex(plane_points_1) + n_influenced = influenced_particles[i] + + plane_points = [plane_points_1[i], plane_points_2[i]] + + plane_size = plane_points[2] - plane_points[1] + flow_directions = [ + normalize([-plane_size[2], plane_size[1]]), + -normalize([-plane_size[2], plane_size[1]]), + ] + + @testset "Flow Direction $j" for j in eachindex(flow_directions) + flow_direction = flow_directions[j] + inflow = InFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + outflow = OutFlow(; plane=plane_points, particle_spacing, density, + flow_direction, open_boundary_layers) + + boundary_zones = [ + inflow, + outflow, + ] + + @testset "`$(nameof(typeof(boundary_zone)))`" for boundary_zone in boundary_zones + sign_ = (boundary_zone isa InFlow) ? 1 : -1 + fluid = extrude_geometry(plane_points; particle_spacing, n_extrude=4, + density, pressure, + direction=(sign_ * flow_direction)) + + fluid_system = EntropicallyDampedSPHSystem(fluid, smoothing_kernel, + buffer_size=0, + density_calculator=ContinuityDensity(), + smoothing_length, sound_speed) + + boundary_system = OpenBoundarySPHSystem(boundary_zone; sound_speed, + fluid_system, buffer_size=0, + reference_velocity, + reference_pressure, + reference_density) + + semi = Semidiscretization(fluid_system, boundary_system) + + ode = semidiscretize(semi, (0.0, 5.0)) + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, boundary_system, semi) + u = TrixiParticles.wrap_u(u0_ode, boundary_system, semi) + + # ==== Characteristic Variables + # `J1 = -sound_speed^2 * (rho - rho_ref) + (p - p_ref)` + # `J2 = rho * sound_speed * (v - v_ref) + (p - p_ref)` + # `J3 = - rho * sound_speed * (v - v_ref) + (p - p_ref)` + function J1(t) + return -sound_speed^2 * (density - reference_density(0, t)) + + (pressure - reference_pressure(0, t)) + end + function J2(t) + return density * sound_speed * + dot(-reference_velocity(0, t), flow_direction) + + (pressure - reference_pressure(0, t)) + end + function J3(t) + return -density * sound_speed * + dot(-reference_velocity(0, t), flow_direction) + + (pressure - reference_pressure(0, t)) + end + + # First evaluation. + # Particles not influenced by the fluid have zero values. + t1 = 2.0 + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, t1) + evaluated_vars1 = boundary_system.characteristics + + if boundary_zone isa InFlow + @test all(isapprox.(evaluated_vars1[1, :], 0.0)) + @test all(isapprox.(evaluated_vars1[2, :], 0.0)) + @test all(isapprox.(evaluated_vars1[3, 1:n_influenced], J3(t1))) + @test all(isapprox.(evaluated_vars1[3, (n_influenced + 1):end], 0.0)) + + elseif boundary_zone isa OutFlow + @test all(isapprox.(evaluated_vars1[1, 1:n_influenced], J1(t1))) + @test all(isapprox.(evaluated_vars1[2, 1:n_influenced], J2(t1))) + @test all(isapprox.(evaluated_vars1[1:2, (n_influenced + 1):end], 0.0)) + @test all(isapprox.(evaluated_vars1[3, :], 0.0)) + end + + # Second evaluation. + # Particles not influenced by the fluid have previous values. + t2 = 3.0 + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, t2) + evaluated_vars2 = boundary_system.characteristics + + if boundary_zone isa InFlow + @test all(isapprox.(evaluated_vars2[1, :], 0.0)) + @test all(isapprox.(evaluated_vars2[2, :], 0.0)) + @test all(isapprox.(evaluated_vars2[3, 1:n_influenced], J3(t2))) + @test all(isapprox.(evaluated_vars2[3, (n_influenced + 1):end], J3(t1))) + + elseif boundary_zone isa OutFlow + @test all(isapprox.(evaluated_vars2[1, 1:n_influenced], J1(t2))) + @test all(isapprox.(evaluated_vars2[2, 1:n_influenced], J2(t2))) + @test all(isapprox.(evaluated_vars2[1, (n_influenced + 1):end], J1(t1))) + @test all(isapprox.(evaluated_vars2[2, (n_influenced + 1):end], J2(t1))) + @test all(isapprox.(evaluated_vars2[3, :], 0.0)) + end + end + end + end +end diff --git a/test/schemes/boundary/open_boundary/open_boundary.jl b/test/schemes/boundary/open_boundary/open_boundary.jl new file mode 100644 index 000000000..5fe09c703 --- /dev/null +++ b/test/schemes/boundary/open_boundary/open_boundary.jl @@ -0,0 +1,2 @@ +include("characteristic_variables.jl") +include("boundary_zone.jl") diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index 2e5e609e5..5cf7cd5bd 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -1,4 +1,5 @@ include("solid/total_lagrangian_sph/total_lagrangian_sph.jl") include("boundary/dummy_particles/dummy_particles.jl") include("boundary/monaghan_kajtar/monaghan_kajtar.jl") +include("boundary/open_boundary/open_boundary.jl") include("fluid/fluid.jl") diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl new file mode 100644 index 000000000..1158aaff0 --- /dev/null +++ b/test/systems/open_boundary_system.jl @@ -0,0 +1,90 @@ +@testset verbose=true "`OpenBoundarySPHSystem`" begin + @testset verbose=true "Illegal Inputs" begin + plane = ([0.0, 0.0], [0.0, 1.0]) + flow_direction = (1.0, 0.0) + + # Mock fluid system + struct FluidSystemMock2 <: TrixiParticles.FluidSystem{2, Nothing} end + + inflow = InFlow(; plane, particle_spacing=0.1, + flow_direction, density=1.0, open_boundary_layers=2) + + error_str = "`reference_velocity` must be either a function mapping " * + "each particle's coordinates and time to its velocity, " * + "an array where the ``i``-th column holds the velocity of particle ``i`` " * + "or, for a constant fluid velocity, a vector of length 2 for a 2D problem holding this velocity" + + reference_velocity = 1.0 + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_velocity) + + error_str = "`reference_pressure` must be either a function mapping " * + "each particle's coordinates and time to its pressure, " * + "a vector holding the pressure of each particle, or a scalar" + + reference_pressure = [1.0, 1.0] + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_pressure) + + error_str = "`reference_density` must be either a function mapping " * + "each particle's coordinates and time to its density, " * + "a vector holding the density of each particle, or a scalar" + + reference_density = [1.0, 1.0] + @test_throws ArgumentError(error_str) OpenBoundarySPHSystem(inflow; sound_speed=1.0, + buffer_size=0, + fluid_system=FluidSystemMock2(), + reference_density) + end + @testset "`show`" begin + inflow = InFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(inflow; sound_speed=1.0, buffer_size=0, + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}(InFlow) with 80 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 80 │ + │ #buffer_particles: ……………………………… 0 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary: ……………………………………………………… InFlow │ + │ flow direction: ……………………………………… [1.0, 0.0] │ + │ prescribed velocity: ………………………… constant_vector │ + │ prescribed pressure: ………………………… constant_scalar │ + │ prescribed density: …………………………… constant_scalar │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + + outflow = OutFlow(; plane=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, + flow_direction=(1.0, 0.0), density=1.0, open_boundary_layers=4) + system = OpenBoundarySPHSystem(outflow; sound_speed=1.0, buffer_size=0, + fluid_system=FluidSystemMock2()) + + show_compact = "OpenBoundarySPHSystem{2}(OutFlow) with 80 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ OpenBoundarySPHSystem{2} │ + │ ════════════════════════ │ + │ #particles: ………………………………………………… 80 │ + │ #buffer_particles: ……………………………… 0 │ + │ fluid system: …………………………………………… FluidSystemMock2 │ + │ boundary: ……………………………………………………… OutFlow │ + │ flow direction: ……………………………………… [1.0, 0.0] │ + │ prescribed velocity: ………………………… constant_vector │ + │ prescribed pressure: ………………………… constant_scalar │ + │ prescribed density: …………………………… constant_scalar │ + │ width: ……………………………………………………………… 0.2 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + end +end diff --git a/test/systems/systems.jl b/test/systems/systems.jl index 6aff417b3..697f575fa 100644 --- a/test/systems/systems.jl +++ b/test/systems/systems.jl @@ -2,4 +2,5 @@ include("wcsph_system.jl") include("edac_system.jl") include("solid_system.jl") include("boundary_system.jl") +include("open_boundary_system.jl") include("dem_system.jl") From ae3a2ffc5bae054a6ad88c40efa018eb77009a97 Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Wed, 19 Jun 2024 15:45:48 +0200 Subject: [PATCH 6/8] New Feature Release (#548) Open Boundaries have been added --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bff9874c1..6f6a24bb1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TrixiParticles" uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"] -version = "0.1.3-dev" +version = "0.1.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From e1e80f0f62514c65a441ebeeb7ff391a927837da Mon Sep 17 00:00:00 2001 From: Sven Berger Date: Thu, 20 Jun 2024 09:37:05 +0200 Subject: [PATCH 7/8] Set to development version 0.1.4 (#549) * Set to development version 0.1.4 * Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6f6a24bb1..8c89b8a64 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TrixiParticles" uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" authors = ["erik.faulhaber <44124897+efaulhaber@users.noreply.github.com>"] -version = "0.1.3" +version = "0.1.4-dev" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 42d854ebd82b66c33d8fc3212fbed505d17614e9 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 21 Jun 2024 09:55:48 +0200 Subject: [PATCH 8/8] Rewrite `@threaded` macro to work with GPUs (#534) * Fix wrapping for GPU arrays * Adapt `@threaded` macro to work with GPUs * Use new `@threaded` macro * Fix NHS to run GPU code with CPU arrays * Import `@index` from KernelAbstractions.jl For some reason it doesn't run on the CPU otherwise. * Add proper docs for `@threaded` * Fix `foreach_neighbor` * Reformat * Add test for GPU code on the CPU * Fix unit tests * Add unit tests for GPU interaction * Also add tests for GPU versions of `wrap_v` and `wrap_u` * Add kwarg `neighborhood_search` to dam break example * Fix tests * Reformat * Add function `wrap_array` * Fix merge * Fix open boundary code * Implement suggestions * Fix syntax error from merge * Fix `wrap_u` and `wrap_v` for GPU array types * Fix constructor of `BoundarySPHSystem` --- Project.toml | 4 +- examples/fluid/dam_break_2d.jl | 6 +- examples/n_body/n_body_system.jl | 4 +- src/TrixiParticles.jl | 6 +- src/general/corrections.jl | 2 +- src/general/general.jl | 30 +---- src/general/gpu.jl | 8 ++ src/general/neighborhood_search.jl | 11 ++ src/general/semidiscretization.jl | 42 ++++--- src/general/system.jl | 25 ++++ .../dummy_particles/dummy_particles.jl | 7 +- src/schemes/boundary/open_boundary/system.jl | 4 +- src/schemes/boundary/system.jl | 8 +- .../fluid/weakly_compressible_sph/system.jl | 2 +- .../solid/total_lagrangian_sph/system.jl | 91 +++++++-------- src/util.jl | 93 +++++++++------ test/count_allocations.jl | 17 +-- test/examples/examples.jl | 13 +++ test/general/semidiscretization.jl | 10 +- .../schemes/solid/total_lagrangian_sph/rhs.jl | 109 ++++++++++++------ 20 files changed, 302 insertions(+), 190 deletions(-) diff --git a/Project.toml b/Project.toml index 8c89b8a64..2565a0945 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,9 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c" @@ -35,6 +37,7 @@ FastPow = "0.1" ForwardDiff = "0.10" JSON = "0.21" MuladdMacro = "0.2" +PointNeighbors = "0.2.3" Polyester = "0.7.5" RecipesBase = "1" Reexport = "1" @@ -43,6 +46,5 @@ StaticArrays = "1" StrideArrays = "0.1" TimerOutputs = "0.5" TrixiBase = "0.1.3" -PointNeighbors = "0.2" WriteVTK = "1" julia = "1.9" diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index e58a7cc06..a1eeb2f38 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -73,8 +73,10 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coef # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, threaded_nhs_update=true) -ode = semidiscretize(semi, tspan) +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch, + threaded_nhs_update=true) +ode = semidiscretize(semi, tspan, data_type=nothing) info_callback = InfoCallback(interval=100) diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index c7226357a..556147ed4 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -1,7 +1,9 @@ using TrixiParticles using LinearAlgebra -struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothing} +# The second type parameter of `System` can't be `Nothing`, or TrixiParticles will launch +# GPU kernel for `for_particle_neighbor` loops. +struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, 0} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] G :: ELTYPE diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index aa97ab720..9fe23e3d4 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -9,7 +9,9 @@ using DataFrames: DataFrame using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect, PresetTimeCallback using FastPow: @fastpow using ForwardDiff: ForwardDiff +using GPUArrays: AbstractGPUArray using JSON: JSON +using KernelAbstractions: KernelAbstractions, @kernel, @index using LinearAlgebra: norm, dot, I, tr, inv, pinv, det using MuladdMacro: @muladd using Polyester: Polyester, @batch @@ -27,7 +29,9 @@ using TrixiBase: trixi_include, @trixi_timeit, timer, timeit_debug_enabled, using PointNeighbors: PointNeighbors, for_particle_neighbor using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save -# util needs to be first because of macro @trixi_timeit +# `util.jl` depends on the `GPUSystem` type defined in `system.jl` +include("general/system.jl") +# `util.jl` needs to be next because of the macros `@trixi_timeit` and `@threaded` include("util.jl") include("callbacks/callbacks.jl") include("general/general.jl") diff --git a/src/general/corrections.jl b/src/general/corrections.jl index cf64f6d16..ee894fec4 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -440,7 +440,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, end function correction_matrix_inversion_step!(corr_matrix, system) - @threaded for particle in eachparticle(system) + @threaded system for particle in eachparticle(system) L = extract_smatrix(corr_matrix, system, particle) # The matrix `L` only becomes singular when the particle and all neighbors diff --git a/src/general/general.jl b/src/general/general.jl index 774d076c9..dfd709861 100644 --- a/src/general/general.jl +++ b/src/general/general.jl @@ -1,37 +1,11 @@ -# Abstract supertype for all system types. We additionally store the type of the system's -# initial condition, which is `Nothing` when using KernelAbstractions.jl. -abstract type System{NDIMS, IC} end - -# When using KernelAbstractions.jl, the initial condition has been replaced by `nothing` -GPUSystem = System{NDIMS, Nothing} where {NDIMS} - -abstract type FluidSystem{NDIMS, IC} <: System{NDIMS, IC} end -timer_name(::FluidSystem) = "fluid" -vtkname(system::FluidSystem) = "fluid" - -abstract type SolidSystem{NDIMS, IC} <: System{NDIMS, IC} end -timer_name(::SolidSystem) = "solid" -vtkname(system::SolidSystem) = "solid" - -abstract type BoundarySystem{NDIMS, IC} <: System{NDIMS, IC} end -timer_name(::BoundarySystem) = "boundary" -vtkname(system::BoundarySystem) = "boundary" - -@inline function set_zero!(du) - du .= zero(eltype(du)) - - return du -end - -# Note that `semidiscretization.jl` depends on the system types and has to be -# included later. +# Note that `system.jl` has already been included. +# `semidiscretization.jl` depends on the system types and has to be included later. # `density_calculators.jl` needs to be included before `corrections.jl`. include("density_calculators.jl") include("corrections.jl") include("smoothing_kernels.jl") include("initial_condition.jl") include("buffer.jl") -include("system.jl") include("interpolation.jl") include("custom_quantities.jl") include("neighborhood_search.jl") diff --git a/src/general/gpu.jl b/src/general/gpu.jl index 377ea8b72..8aa481cff 100644 --- a/src/general/gpu.jl +++ b/src/general/gpu.jl @@ -13,6 +13,7 @@ Adapt.@adapt_structure DensityDiffusionAntuono Adapt.@adapt_structure BoundarySPHSystem Adapt.@adapt_structure BoundaryModelDummyParticles Adapt.@adapt_structure BoundaryModelMonaghanKajtar +Adapt.@adapt_structure TotalLagrangianSPHSystem # The initial conditions are only used for initialization, which happens before `adapt`ing # the semidiscretization, so we don't need to store `InitialCondition`s on the GPU. @@ -32,3 +33,10 @@ end function Adapt.adapt_structure(to::typeof(Array), range::UnitRange) return range end + +KernelAbstractions.get_backend(::PtrArray) = KernelAbstractions.CPU() +KernelAbstractions.get_backend(system::System) = KernelAbstractions.get_backend(system.mass) + +function KernelAbstractions.get_backend(system::BoundarySPHSystem) + KernelAbstractions.get_backend(system.coordinates) +end diff --git a/src/general/neighborhood_search.jl b/src/general/neighborhood_search.jl index 0713f153d..afcd2c9b3 100644 --- a/src/general/neighborhood_search.jl +++ b/src/general/neighborhood_search.jl @@ -9,3 +9,14 @@ function PointNeighbors.for_particle_neighbor(f, system, neighbor_system, for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particles=particles, parallel=parallel) end + +function PointNeighbors.for_particle_neighbor(f, system::GPUSystem, neighbor_system, + system_coords, neighbor_coords, + neighborhood_search; + particles=each_moving_particle(system), + parallel=true) + @threaded system for particle in particles + PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, particle) + end +end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 73c1fe56d..7f6c791db 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -368,6 +368,17 @@ end # We have to pass `system` here for type stability, # since the type of `system` determines the return type. +@inline function wrap_v(v_ode, system, semi) + (; ranges_v) = semi + + range = ranges_v[system_indices(system, semi)] + + @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) + + return wrap_array(v_ode, range, + (StaticInt(v_nvariables(system)), n_moving_particles(system))) +end + @inline function wrap_u(u_ode, system, semi) (; ranges_u) = semi @@ -375,22 +386,23 @@ end @boundscheck @assert length(range) == u_nvariables(system) * n_moving_particles(system) - # This is a non-allocating version of: - # return unsafe_wrap(Array{eltype(u_ode), 2}, pointer(view(u_ode, range)), - # (u_nvariables(system), n_moving_particles(system))) - return PtrArray(pointer(view(u_ode, range)), - (StaticInt(u_nvariables(system)), n_moving_particles(system))) + return wrap_array(u_ode, range, + (StaticInt(u_nvariables(system)), n_moving_particles(system))) end -@inline function wrap_v(v_ode, system, semi) - (; ranges_v) = semi - - range = ranges_v[system_indices(system, semi)] - - @boundscheck @assert length(range) == v_nvariables(system) * n_moving_particles(system) +@inline function wrap_array(array::Array, range, size) + # This is a non-allocating version of: + # return unsafe_wrap(Array{eltype(array), 2}, pointer(view(array, range)), size) + return PtrArray(pointer(view(array, range)), size) +end - return PtrArray(pointer(view(v_ode, range)), - (StaticInt(v_nvariables(system)), n_moving_particles(system))) +@inline function wrap_array(array, range, size) + # For non-`Array`s (typically GPU arrays), just reshape. Calling the `PtrArray` code + # above for a `CuArray` yields another `CuArray` (instead of a `PtrArray`) + # and is 8 times slower with double the allocations. + # + # Note that `size` might contain `StaticInt`s, so convert to `Int` first. + return reshape(view(array, range), Int.(size)) end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) @@ -409,7 +421,7 @@ function drift!(du_ode, v_ode, u_ode, semi, t) du = wrap_u(du_ode, system, semi) v = wrap_v(v_ode, system, semi) - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) # This can be dispatched per system add_velocity!(du, v, particle, system) end @@ -508,7 +520,7 @@ function add_source_terms!(dv_ode, v_ode, u_ode, semi) v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) # Dispatch by system type to exclude boundary systems add_acceleration!(dv, particle, system) add_source_terms_inner!(dv, v, u, particle, system, source_terms(system)) diff --git a/src/general/system.jl b/src/general/system.jl index c3f10caad..674555aea 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -1,3 +1,28 @@ +# Abstract supertype for all system types. We additionally store the type of the system's +# initial condition, which is `Nothing` when using KernelAbstractions.jl. +abstract type System{NDIMS, IC} end + +# When using KernelAbstractions.jl, the initial condition has been replaced by `nothing` +GPUSystem = System{NDIMS, Nothing} where {NDIMS} + +abstract type FluidSystem{NDIMS, IC} <: System{NDIMS, IC} end +timer_name(::FluidSystem) = "fluid" +vtkname(system::FluidSystem) = "fluid" + +abstract type SolidSystem{NDIMS, IC} <: System{NDIMS, IC} end +timer_name(::SolidSystem) = "solid" +vtkname(system::SolidSystem) = "solid" + +abstract type BoundarySystem{NDIMS, IC} <: System{NDIMS, IC} end +timer_name(::BoundarySystem) = "boundary" +vtkname(system::BoundarySystem) = "boundary" + +@inline function set_zero!(du) + du .= zero(eltype(du)) + + return du +end + initialize!(system, neighborhood_search) = system @inline Base.ndims(::System{NDIMS}) where {NDIMS} = NDIMS diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 0f7b42696..3e3abcbd9 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -291,7 +291,7 @@ function compute_pressure!(boundary_model, ::Union{SummationDensity, ContinuityD # Limit pressure to be non-negative to avoid attractive forces between fluid and # boundary particles at free surfaces (sticking artifacts). - @threaded for particle in eachparticle(system) + @threaded system for particle in eachparticle(system) apply_state_equation!(boundary_model, particle_density(v, boundary_model, particle), particle) end @@ -346,14 +346,15 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system_coords, neighbor_coords, v_neighbor_system, nhs) end - for particle in eachparticle(system) + + @threaded system for particle in eachparticle(system) # Limit pressure to be non-negative to avoid attractive forces between fluid and # boundary particles at free surfaces (sticking artifacts). pressure[particle] = max(pressure[particle], 0.0) end end - @trixi_timeit timer() "inverse state equation" @threaded for particle in eachparticle(system) + @trixi_timeit timer() "inverse state equation" @threaded system for particle in eachparticle(system) compute_adami_density!(boundary_model, system, system_coords, particle) end end diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 43c0a6ff1..b5950856e 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -171,7 +171,7 @@ end reference_velocity, reference_pressure, reference_density) = system # Update quantities based on the characteristic variables - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) particle_position = current_coords(u, system, particle) J1 = characteristics[1, particle] @@ -230,7 +230,7 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) # Thus, we compute the characteristics for the particles that are outside the influence # of fluid particles by using the average of the values of the previous time step. # See eq. 27 in Negi (2020) https://doi.org/10.1016/j.cma.2020.113119 - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) # Particle is outside of the influence of fluid particles if isapprox(volume[particle], 0.0) diff --git a/src/schemes/boundary/system.jl b/src/schemes/boundary/system.jl index 73af6a843..f36d515ac 100644 --- a/src/schemes/boundary/system.jl +++ b/src/schemes/boundary/system.jl @@ -27,13 +27,13 @@ struct BoundarySPHSystem{BM, NDIMS, ELTYPE <: Real, IC, CO, M, IM, # This constructor is necessary for Adapt.jl to work with this struct. # See the comments in general/gpu.jl for more details. function BoundarySPHSystem(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache) + ismoving, adhesion_coefficient, cache, buffer) ELTYPE = eltype(coordinates) new{typeof(boundary_model), size(coordinates, 1), ELTYPE, typeof(initial_condition), typeof(coordinates), typeof(movement), typeof(ismoving), typeof(cache)}(initial_condition, coordinates, boundary_model, movement, - ismoving, adhesion_coefficient, cache, nothing) + ismoving, adhesion_coefficient, cache, buffer) end end @@ -54,7 +54,7 @@ function BoundarySPHSystem(initial_condition, model; movement=nothing, # Because of dispatches boundary model needs to be first! return BoundarySPHSystem(initial_condition, coordinates, model, movement, - ismoving, adhesion_coefficient, cache) + ismoving, adhesion_coefficient, cache, nothing) end """ @@ -209,7 +209,7 @@ function (movement::BoundaryMovement)(system, t) is_moving(t) || return system - @threaded for particle in moving_particles + @threaded system for particle in moving_particles pos_new = initial_coords(system, particle) + movement_function(t) vel = ForwardDiff.derivative(movement_function, t) acc = ForwardDiff.derivative(t_ -> ForwardDiff.derivative(movement_function, t_), t) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index a8e28670b..44e71695b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -311,7 +311,7 @@ function reinit_density!(system, v, u, v_ode, u_ode, semi) end function compute_pressure!(system, v) - @threaded for particle in eachparticle(system) + @threaded system for particle in eachparticle(system) apply_state_equation!(system, particle_density(v, system, particle), particle) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index 9203425fb..125bb5507 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -70,55 +70,50 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, IC, ARRAY1D, ARRAY2D, penalty_force :: PF source_terms :: ST buffer :: Nothing +end - function TotalLagrangianSPHSystem(initial_condition, - smoothing_kernel, smoothing_length, - young_modulus, poisson_ratio; - n_fixed_particles=0, boundary_model=nothing, - acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel)), - penalty_force=nothing, source_terms=nothing) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) - n_particles = nparticles(initial_condition) - - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end +function TotalLagrangianSPHSystem(initial_condition, + smoothing_kernel, smoothing_length, + young_modulus, poisson_ratio; + n_fixed_particles=0, boundary_model=nothing, + acceleration=ntuple(_ -> 0.0, + ndims(smoothing_kernel)), + penalty_force=nothing, source_terms=nothing) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + n_particles = nparticles(initial_condition) - # Make acceleration an SVector - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end + if ndims(smoothing_kernel) != NDIMS + throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) + end - initial_coordinates = copy(initial_condition.coordinates) - current_coordinates = copy(initial_condition.coordinates) - mass = copy(initial_condition.mass) - material_density = copy(initial_condition.density) - correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) - pk1_corrected = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) - deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) - - n_moving_particles = n_particles - n_fixed_particles - - lame_lambda = young_modulus * poisson_ratio / - ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) - lame_mu = 0.5 * young_modulus / (1 + poisson_ratio) - - return new{typeof(boundary_model), NDIMS, ELTYPE, - typeof(initial_condition), - typeof(mass), typeof(initial_coordinates), - typeof(deformation_grad), typeof(smoothing_kernel), - typeof(penalty_force), - typeof(source_terms)}(initial_condition, initial_coordinates, - current_coordinates, mass, correction_matrix, - pk1_corrected, deformation_grad, material_density, - n_moving_particles, young_modulus, poisson_ratio, - lame_lambda, lame_mu, smoothing_kernel, - smoothing_length, acceleration_, boundary_model, - penalty_force, source_terms, nothing) + # Make acceleration an SVector + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + + initial_coordinates = copy(initial_condition.coordinates) + current_coordinates = copy(initial_condition.coordinates) + mass = copy(initial_condition.mass) + material_density = copy(initial_condition.density) + correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) + pk1_corrected = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) + deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles) + + n_moving_particles = n_particles - n_fixed_particles + + lame_lambda = young_modulus * poisson_ratio / + ((1 + poisson_ratio) * (1 - 2 * poisson_ratio)) + lame_mu = 0.5 * young_modulus / (1 + poisson_ratio) + + return TotalLagrangianSPHSystem(initial_condition, initial_coordinates, + current_coordinates, mass, correction_matrix, + pk1_corrected, deformation_grad, material_density, + n_moving_particles, young_modulus, poisson_ratio, + lame_lambda, lame_mu, smoothing_kernel, + smoothing_length, acceleration_, boundary_model, + penalty_force, source_terms, nothing) end function Base.show(io::IO, system::TotalLagrangianSPHSystem) @@ -260,7 +255,7 @@ end calc_deformation_grad!(deformation_grad, neighborhood_search, system) - @threaded for particle in eachparticle(system) + @threaded system for particle in eachparticle(system) F_particle = deformation_gradient(system, particle) pk1_particle = pk1_stress_tensor(F_particle, system) pk1_particle_corrected = pk1_particle * correction_matrix(system, particle) @@ -397,7 +392,7 @@ end function von_mises_stress(system::TotalLagrangianSPHSystem) von_mises_stress_vector = zeros(eltype(system.pk1_corrected), nparticles(system)) - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) von_mises_stress_vector[particle] = von_mises_stress(system, particle) end @@ -430,7 +425,7 @@ function cauchy_stress(system::TotalLagrangianSPHSystem) cauchy_stress_tensors = zeros(eltype(system.pk1_corrected), NDIMS, NDIMS, nparticles(system)) - @threaded for particle in each_moving_particle(system) + @threaded system for particle in each_moving_particle(system) F = deformation_gradient(system, particle) J = det(F) P = pk1_corrected(system, particle) diff --git a/src/util.jl b/src/util.jl index 338a94415..172f0a800 100644 --- a/src/util.jl +++ b/src/util.jl @@ -27,57 +27,74 @@ function print_startup_message() end """ - @threaded for ... end + @threaded system for ... end +Run either a threaded CPU loop or launch a kernel on the GPU, depending on the type of `system`. Semantically the same as `Threads.@threads` when iterating over a `AbstractUnitRange` but without guarantee that the underlying implementation uses `Threads.@threads` or works for more general `for` loops. -In particular, there may be an additional check whether only one thread is used -to reduce the overhead of serial execution or the underlying threading capabilities -might be provided by other packages such as [Polyester.jl](https://github.com/JuliaSIMD/Polyester.jl). +The first argument must either be a particle system or an array from which can be derived +if the loop has to be run threaded on the CPU or launched as a kernel on the GPU. + +In particular, the underlying threading capabilities might be provided by other packages +such as [Polyester.jl](https://github.com/JuliaSIMD/Polyester.jl). !!! warn This macro does not necessarily work for general `for` loops. For example, it does not necessarily support general iterables such as `eachline(filename)`. - -Some discussion can be found at -[https://discourse.julialang.org/t/overhead-of-threads-threads/53964](https://discourse.julialang.org/t/overhead-of-threads-threads/53964) -and -[https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435](https://discourse.julialang.org/t/threads-threads-with-one-thread-how-to-remove-the-overhead/58435). - -Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl). """ -macro threaded(expr) - # Use `esc(quote ... end)` for nested macro calls as suggested in - # https://github.com/JuliaLang/julia/issues/23221 - # - # The following code is a simple version using only `Threads.@threads` from the - # standard library with an additional check whether only a single thread is used - # to reduce some overhead (and allocations) for serial execution. - # - # return esc(quote - # let - # if Threads.nthreads() == 1 - # $(expr) - # else - # Threads.@threads $(expr) - # end - # end - # end) - # - # However, the code below using `@batch` from Polyester.jl is more efficient, - # since this packages provides threads with less overhead. Since it is written - # by Chris Elrod, the author of LoopVectorization.jl, we expect this package - # to provide the most efficient and useful implementation of threads (as we use - # them) available in Julia. - # !!! danger "Heisenbug" - # Look at the comments for `wrap_array` when considering to change this macro. - +macro threaded(system, expr) + # Reverse-engineer the for loop. + # `expr.args[1]` is the head of the for loop, like `i = eachindex(x)`. + # So, `expr.args[1].args[2]` is the iterator `eachindex(x)` + # and `expr.args[1].args[1]` is the loop variable `i`. + iterator = expr.args[1].args[2] + i = expr.args[1].args[1] + inner_loop = expr.args[2] + + # Assemble the for loop again as a call to `parallel_foreach`, using `$i` to use the + # same loop variable as used in the for loop. return esc(quote - TrixiParticles.@batch $(expr) + TrixiParticles.parallel_foreach($iterator, $system) do $i + $inner_loop + end end) end +# Use `Polyester.@batch` for low-overhead threading +@inline function parallel_foreach(f, iterator, system) + Polyester.@batch for i in iterator + @inline f(i) + end +end + +# On GPUs, execute `f` inside a GPU kernel with KernelAbstractions.jl +@inline function parallel_foreach(f, iterator, system::Union{GPUSystem, AbstractGPUArray}) + # On the GPU, we can only loop over `1:N`. Therefore, we loop over `1:length(iterator)` + # and index with `iterator[eachindex(iterator)[i]]`. + # Note that this only works with vector-like iterators that support arbitrary indexing. + indices = eachindex(iterator) + ndrange = length(indices) + + # Skip empty loops + ndrange == 0 && return + + backend = KernelAbstractions.get_backend(system) + + # Call the generic kernel that is defined below, which only calls a function with + # the global GPU index. + generic_kernel(backend)(ndrange=ndrange) do i + @inline f(iterator[indices[i]]) + end + + KernelAbstractions.synchronize(backend) +end + +@kernel function generic_kernel(f) + i = @index(Global) + @inline f(i) +end + @doc raw""" examples_dir() diff --git a/test/count_allocations.jl b/test/count_allocations.jl index 1c597f9a0..5957bc67a 100644 --- a/test/count_allocations.jl +++ b/test/count_allocations.jl @@ -16,14 +16,15 @@ function copy_semi_with_no_update_nhs(semi) neighborhood_searches) end -# Forward `for_particle_neighbor` to wrapped neighborhood search -@inline function TrixiParticles.for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search::NoUpdateNeighborhoodSearch; - particles=axes(system_coords, 2), - parallel=true) - TrixiParticles.for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search.nhs, - particles=particles, parallel=parallel) +# Forward `foreach_neighbor` to wrapped neighborhood search +@inline function TrixiParticles.PointNeighbors.foreach_neighbor(f, system_coords, + neighbor_coords, + neighborhood_search::NoUpdateNeighborhoodSearch, + particle; + search_radius=neighborhood_search.nhs.search_radius) + TrixiParticles.PointNeighbors.foreach_neighbor(f, system_coords, neighbor_coords, + neighborhood_search.nhs, particle, + search_radius=search_radius) end # No update diff --git a/test/examples/examples.jl b/test/examples/examples.jl index dbe03a130..19b6ba2b9 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -87,6 +87,19 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/dam_break_2d.jl with KernelAbstractions.jl" begin + # Emulate the GPU code on the CPU by passing `data_type = Array` + @test_nowarn_mod trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d.jl"), tspan=(0.0, 0.1), + data_type=Array) [ + r"┌ Info: The desired tank length in y-direction .*\n", + r"└ New tank length in y-direction.*\n", + ] + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/dam_break_3d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 4aff1f4fd..457b5b84d 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -1,9 +1,15 @@ # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "Semidiscretization" begin - # Mock systems - struct System1 <: TrixiParticles.System{3, Nothing} end + # Mock systems. `System1` will use the CPU backend, `System2` is a `GPUSystem`, using + # the GPU backend (emulated on the CPU). + struct System1 <: TrixiParticles.System{3, String} end struct System2 <: TrixiParticles.System{3, Nothing} end + # `System2` has no field `mass`, so we have to manually define the backend + function TrixiParticles.KernelAbstractions.get_backend(::System2) + return TrixiParticles.KernelAbstractions.CPU() + end + system1 = System1() system2 = System2() diff --git a/test/schemes/solid/total_lagrangian_sph/rhs.jl b/test/schemes/solid/total_lagrangian_sph/rhs.jl index dfc9fc234..b64fb6c97 100644 --- a/test/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/test/schemes/solid/total_lagrangian_sph/rhs.jl @@ -32,7 +32,7 @@ [0.0, 0.0], ] - @testset "Test $i" for i in 1:4 + @testset verbose=true "Test $i" for i in 1:4 #### Setup each_moving_particle = [particle[i]] # Only calculate dv for this one particle eachparticle = [particle[i], neighbor[i]] @@ -53,18 +53,27 @@ kernel_deriv = 1.0 #### Mocking - # Mock the system - system = Val{:mock_system_interact}() - TrixiParticles.ndims(::Val{:mock_system_interact}) = 2 - Base.ntuple(f, ::Symbol) = ntuple(f, 2) # Make `extract_svector` work + # Mock a CPU system to test CPU code + struct MockSystemInteractCPU <: TrixiParticles.System{2, String} end + system = MockSystemInteractCPU() - function TrixiParticles.initial_coordinates(::Val{:mock_system_interact}) + # Mock a GPU system to emulate GPU code on the CPU + struct MockSystemInteractGPU <: TrixiParticles.System{2, Nothing} end + system_gpu = MockSystemInteractGPU() + + function TrixiParticles.KernelAbstractions.get_backend(::MockSystemInteractGPU) + return TrixiParticles.KernelAbstractions.CPU() + end + + MockSystemType = Union{MockSystemInteractCPU, MockSystemInteractGPU} + + function TrixiParticles.initial_coordinates(::MockSystemType) return initial_coordinates end # Unpack calls should return predefined values or # another mock object of the type Val{:mock_property_name}. - function Base.getproperty(::Val{:mock_system_interact}, f::Symbol) + function Base.getproperty(::MockSystemType, f::Symbol) if f === :current_coordinates return current_coordinates elseif f === :material_density @@ -81,10 +90,13 @@ return Val(Symbol("mock_" * string(f))) end - function TrixiParticles.each_moving_particle(::Val{:mock_system_interact}) + function TrixiParticles.each_moving_particle(::MockSystemType) each_moving_particle end - TrixiParticles.eachparticle(::Val{:mock_system_interact}) = eachparticle + TrixiParticles.eachparticle(::MockSystemType) = eachparticle + + # Mock the neighborhood search + nhs = Val{:nhs}() TrixiParticles.PointNeighbors.eachneighbor(_, ::Val{:nhs}) = eachneighbor function Base.getproperty(::Val{:nhs}, f::Symbol) @@ -99,26 +111,25 @@ end TrixiParticles.ndims(::Val{:nhs}) = 2 - function TrixiParticles.pk1_corrected(::Val{:mock_system_dv}, particle_) - if particle_ == particle[i] - return pk1_particle_corrected[i] - end - return pk1_neighbor_corrected[i] - end - - function TrixiParticles.add_acceleration!(_, _, ::Val{:mock_system_interact}) + function TrixiParticles.add_acceleration!(_, _, ::MockSystemType) nothing end TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv #### Verification - dv = zeros(ndims(system), 10) - dv_expected = copy(dv) - dv_expected[:, particle[i]] = dv_particle_expected[i] + systems = [system, system_gpu] + names = ["CPU code", "Emulate GPU"] + @testset "$(names[j])" for j in eachindex(names) + system_ = systems[j] + + dv = zeros(ndims(system_), 10) + dv_expected = copy(dv) + dv_expected[:, particle[i]] = dv_particle_expected[i] - TrixiParticles.interact_solid_solid!(dv, Val(:nhs), system, system) + TrixiParticles.interact_solid_solid!(dv, nhs, system_, system_) - @test dv ≈ dv_expected + @test dv ≈ dv_expected + end end end @@ -140,7 +151,7 @@ 10 / 1000^2 * 1.5400218087591082 * 324.67072684047224 * 1.224, 0.0, ]) - @testset "Deformation Function: $deformation" for deformation in keys(deformations) + @testset verbose=true "Deformation Function: $deformation" for deformation in keys(deformations) J = deformations[deformation] u = zeros(2, 81) v = zeros(2, 81) @@ -176,22 +187,50 @@ semi = Semidiscretization(system) tspan = (0.0, 1.0) - semidiscretize(semi, tspan) - # Apply the deformation matrix - for particle in axes(u, 2) - # Apply deformation - u[1:2, particle] = deformations[deformation](coordinates[:, particle]) + # To make the code below work + function TrixiParticles.PtrArray{Float64}(::UndefInitializer, length) + TrixiParticles.PtrArray(zeros(length)) end - #### Verification for the particle in the middle - particle = 41 + # We can pass the data type `Array` to convert all systems to `GPUSystem`s + # and emulate the GPU kernels on the GPU. + # But this doesn't test `wrap_v` and `wrap_u` for non-`Array` types. + # In order to test this as well, we need a different data type, so we also + # pass `PtrArray`. + names = ["CPU code", "GPU code with CPU wrapping", "GPU code with GPU wrapping"] + data_types = [nothing, Array, TrixiParticles.PtrArray] + @testset "$(names[i])" for i in eachindex(names) + data_type = data_types[i] + ode = semidiscretize(semi, tspan, data_type=data_type) + + # Apply the deformation matrix + for particle in axes(u, 2) + # Apply deformation + u[1:2, particle] = deformations[deformation](coordinates[:, particle]) + end + + v_ode = ode.u0.x[1] + if isnothing(data_type) + u_ode = vec(u) + else + u_ode = data_type(vec(u)) + end + + @test typeof(v_ode) == typeof(u_ode) + @test length(v_ode) == length(u_ode) + + #### Verification for the particle in the middle + particle = 41 - dv = zeros(ndims(system), 81) - TrixiParticles.kick!(dv, v, u, semi, 0.0) + dv_ode = zero(v_ode) + TrixiParticles.kick!(dv_ode, v_ode, u_ode, ode.p, 0.0) - @test isapprox(dv[:, particle], dv_expected_41[deformation], - rtol=sqrt(eps()), atol=sqrt(eps())) + dv = TrixiParticles.wrap_v(dv_ode, system, semi) + + @test isapprox(dv[:, particle], dv_expected_41[deformation], + rtol=sqrt(eps()), atol=sqrt(eps())) + end end end -end +end;