diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ff23f9490..4761de950 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -22,12 +22,13 @@ using Printf: @printf, @sprintf using ReadVTK: ReadVTK using RecipesBase: RecipesBase, @series using Random: seed! -using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, - get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate! +using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, + get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!, + add_tstop! @reexport using StaticArrays: SVector using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt -using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer! +using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!, @notimeit using TrixiBase: @trixi_timeit, timer, timeit_debug_enabled, disable_debug_timings, enable_debug_timings, TrixiBase @reexport using TrixiBase: trixi_include, trixi_include_changeprecision @@ -66,7 +67,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, - ParticleShiftingCallback + ParticleShiftingCallback, SplitIntegrationCallback export ContinuityDensity, SummationDensity export PenaltyForceGanzenmueller, TransportVelocityAdami export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel, diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 062a16e3a..5d058794f 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -33,3 +33,4 @@ include("stepsize.jl") include("update.jl") include("steady_state_reached.jl") include("particle_shifting.jl") +include("split_integration.jl") diff --git a/src/callbacks/split_integration.jl b/src/callbacks/split_integration.jl new file mode 100644 index 000000000..441fcf029 --- /dev/null +++ b/src/callbacks/split_integration.jl @@ -0,0 +1,254 @@ +mutable struct SplitIntegrationCallback + integrator :: Any + alg :: Any + kwargs :: Any +end + +@doc raw""" + SplitIntegrationCallback() + +Callback to... +""" +function SplitIntegrationCallback(alg; kwargs...) + split_integration_callback = SplitIntegrationCallback(nothing, alg, kwargs) + + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(split_integration_callback, split_integration_callback, + initialize=(initialize_split_integration!), + save_positions=(false, false)) +end + +function initialize_split_integration!(cb, u, t, integrator) + semi = integrator.p + split_integration_callback = cb.affect! + (; alg, kwargs) = split_integration_callback + + # Disable TLSPH integration in the original integrator + semi.integrate_tlsph[] = false + + # Create split integrator with TLSPH systems only + systems = filter(i -> i isa TotalLagrangianSPHSystem, semi.systems) + + # These neighborhood searches are never used + semi_split = Semidiscretization(systems..., + neighborhood_search=TrivialNeighborhoodSearch{ndims(first(systems))}(), + parallelization_backend=semi.parallelization_backend) + + sizes_u = (u_nvariables(system) * n_moving_particles(system) for system in systems) + sizes_v = (v_nvariables(system) * n_moving_particles(system) for system in systems) + + v_ode, u_ode = integrator.u.x + v0_ode_split = similar(v_ode, sum(sizes_v)) + u0_ode_split = similar(u_ode, sum(sizes_u)) + + # Copy the initial conditions to the split integrator + copy_to_split!(v_ode, u_ode, v0_ode_split, u0_ode_split, semi, semi_split) + + # A zero `tspan` sets `tdir` to zero, which breaks adding tstops + tspan = (integrator.t, integrator.t + 1) + p = (; v_ode, u_ode, semi, semi_split) + ode_split = DynamicalODEProblem(kick_split!, drift_split!, v0_ode_split, u0_ode_split, + tspan, p) + + # Create the split integrator. + # We need the timer here to keep the output clean because this will call `kick!` once. + @trixi_timeit timer() "split integration" begin + @trixi_timeit timer() "init" begin + TimerOutputs.@notimeit timer() split_integrator = SciMLBase.init(ode_split, alg; save_everystep=false, kwargs...) + end + end + + # Remove the `tstop` (equivalent to zero `tspan`) + SciMLBase.pop_tstop!(split_integrator) + + # Store the integrator in the callback + split_integration_callback.integrator = split_integrator + + return cb +end + +# `condition` +function (split_integration_callback::SplitIntegrationCallback)(u, t, integrator) + # Integrate TLSPH after every time step + return true +end + +# `affect!` +function (split_integration_callback::SplitIntegrationCallback)(integrator) + # Function barrier for type stability + affect_inner!(integrator, split_integration_callback.integrator) +end + +function affect_inner!(integrator, split_integrator) + semi_split = split_integrator.p.semi_split + + semi = integrator.p + new_t = integrator.t + + v_ode, u_ode = integrator.u.x + + @assert semi == split_integrator.p.semi + split_integrator.p = (; v_ode, u_ode, semi, semi_split) + + @trixi_timeit timer() "split integration" begin + # 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 and nhs" update_systems_and_nhs(v_ode, u_ode, + semi, new_t; + update_from_callback=true) + + # Integrate the split integrator to the new time + add_tstop!(split_integrator, new_t) + @trixi_timeit timer() "solve" SciMLBase.solve!(split_integrator) + + v_ode_split, u_ode_split = split_integrator.u.x + + # Copy the solutions back to the original integrator + @trixi_timeit timer() "copy back" copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split) + end + + # Tell OrdinaryDiffEq that `u` has been modified + u_modified!(integrator, true) + + return integrator +end + +function kick_split!(dv_ode_split, v_ode_split, u_ode_split, p, t) + (; v_ode, u_ode, semi, semi_split) = p + + @trixi_timeit timer() "reset ∂v/∂t" set_zero!(dv_ode_split) + + # Only loop over TLSPH systems + @trixi_timeit timer() "copy to large v,u" copy_from_split!(v_ode, u_ode, + v_ode_split, u_ode_split, + semi, semi_split) + + # Update the TLSPH systems with the other systems as neighbors + update_nhs_fun = (semi, + u_ode) -> update_nhs_split!(semi, u_ode, u_ode_split, + semi_split) + @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, + semi, t; + systems=semi_split.systems, + update_nhs_fun) + + @trixi_timeit timer() "system interaction" begin + # Only loop over systems in the split integrator + foreach_system(semi_split) do system + # Loop over all neighbors in the big integrator + foreach_system(semi) do neighbor + # Construct string for the interactions timer. + # Avoid allocations from string construction when no timers are used. + if timeit_debug_enabled() + system_index = system_indices(system, semi) + neighbor_index = system_indices(neighbor, semi) + timer_str = "$(timer_name(system))$system_index-$(timer_name(neighbor))$neighbor_index" + else + timer_str = "" + end + + dv = wrap_v(dv_ode_split, system, semi_split) + v_system = wrap_v(v_ode_split, system, semi_split) + u_system = wrap_u(u_ode_split, system, semi_split) + + v_neighbor = wrap_v(v_ode, neighbor, semi) + u_neighbor = wrap_u(u_ode, neighbor, semi) + + @trixi_timeit timer() timer_str begin + interact!(dv, v_system, u_system, v_neighbor, u_neighbor, + system, neighbor, semi; integrate_tlsph=true) + end + end + end + end + + @trixi_timeit timer() "source terms" add_source_terms!(dv_ode_split, v_ode_split, + u_ode_split, semi, t; + semi_wrap=semi_split) +end + +function drift_split!(du_ode, v_ode, u_ode, p, t) + drift!(du_ode, v_ode, u_ode, p.semi_split, t) +end + +function update_nhs_split!(semi, u_ode, u_ode_split, semi_split) + # Only loop over systems in the split integrator + foreach_system(semi_split) do system + u_system = wrap_u(u_ode_split, system, semi_split) + + # Loop over neighbors in the big integrator + foreach_system(semi) do neighbor + # Static NHS for solid-solid (same system) and no interaction for two distinct + # solid systems (TODO). + if !(neighbor isa TotalLagrangianSPHSystem) + u_neighbor = wrap_u(u_ode, neighbor, semi) + neighborhood_search = get_neighborhood_search(system, neighbor, semi) + + # Only the TLSPH particles are moving. All other systems are frozen. + # Note that this does nothing when a grid NHS is used. + update_nhs!(neighborhood_search, system, neighbor, u_system, u_neighbor, + semi; points_moving=(true, false)) + end + end + end +end + +# Copy the solution from the large integrator to the split integrator +@inline function copy_to_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split) + foreach_system(semi_split) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + v_split = wrap_v(v_ode_split, system, semi_split) + u_split = wrap_u(u_ode_split, system, semi_split) + + @threaded semi for particle in each_moving_particle(system) + for i in axes(v, 1) + v_split[i, particle] = v[i, particle] + end + + for i in axes(u, 1) + u_split[i, particle] = u[i, particle] + end + end + end +end + +# Copy the solution from the split integrator to the large integrator +@inline function copy_from_split!(v_ode, u_ode, v_ode_split, u_ode_split, semi, semi_split) + foreach_system(semi_split) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + v_split = wrap_v(v_ode_split, system, semi_split) + u_split = wrap_u(u_ode_split, system, semi_split) + + @threaded semi for particle in each_moving_particle(system) + for i in axes(v, 1) + v[i, particle] = v_split[i, particle] + end + + for i in axes(u, 1) + u[i, particle] = u_split[i, particle] + end + end + end +end + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback}) + @nospecialize cb # reduce precompilation time + print(io, "SplitIntegrationCallback(alg=", cb.affect!.alg, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:SplitIntegrationCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + update_cb = cb.affect! + setup = [ + "alg" => update_cb.alg + ] + summary_box(io, "SplitIntegrationCallback", setup) + end +end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index ae956e88d..b5ffe0c63 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -47,23 +47,25 @@ semi = Semidiscretization(fluid_system, boundary_system, └──────────────────────────────────────────────────────────────────────────────────────────────────┘ ``` """ -struct Semidiscretization{BACKEND, S, RU, RV, NS} +struct Semidiscretization{BACKEND, S, RU, RV, NS, IT} systems :: S ranges_u :: RU ranges_v :: RV neighborhood_searches :: NS parallelization_backend :: BACKEND + integrate_tlsph :: IT # `false` if TLSPH integration is decoupled # Dispatch at `systems` to distinguish this constructor from the one below when # 4 systems are passed. # This is an internal constructor only used in `test/count_allocations.jl` # and by Adapt.jl. function Semidiscretization(systems::Tuple, ranges_u, ranges_v, neighborhood_searches, - parallelization_backend::PointNeighbors.ParallelizationBackend) + parallelization_backend::PointNeighbors.ParallelizationBackend, + integrate_tlsph) new{typeof(parallelization_backend), typeof(systems), typeof(ranges_u), - typeof(ranges_v), typeof(neighborhood_searches)}(systems, ranges_u, ranges_v, - neighborhood_searches, - parallelization_backend) + typeof(ranges_v), typeof(neighborhood_searches), + typeof(integrate_tlsph)}(systems, ranges_u, ranges_v, neighborhood_searches, + parallelization_backend, integrate_tlsph) end end @@ -93,7 +95,7 @@ function Semidiscretization(systems::Union{System, Nothing}...; for system in systems) return Semidiscretization(systems, ranges_u, ranges_v, searches, - parallelization_backend) + parallelization_backend, Ref(true)) end # Inline show function e.g. Semidiscretization(neighborhood_search=...) @@ -328,7 +330,7 @@ function semidiscretize(semi, tspan; reset_threads=true) semi_new = Semidiscretization(set_system_links.(semi_.systems, Ref(semi_)), semi_.ranges_u, semi_.ranges_v, semi_.neighborhood_searches, - semi_.parallelization_backend) + semi_.parallelization_backend, semi_.integrate_tlsph) else semi_new = semi end @@ -454,7 +456,7 @@ function drift!(du_ode, v_ode, u_ode, semi, t) @threaded semi for particle in each_moving_particle(system) # This can be dispatched per system - add_velocity!(du, v, particle, system) + add_velocity!(du, v, particle, system, semi) end end end @@ -463,6 +465,18 @@ function drift!(du_ode, v_ode, u_ode, semi, t) return du_ode end +@inline function add_velocity!(du, v, particle, system, semi::Semidiscretization) + add_velocity!(du, v, particle, system) +end + +@inline function add_velocity!(du, v, particle, system::TotalLagrangianSPHSystem, + semi::Semidiscretization) + # Only add velocity for TLSPH systems if they are integrated + if semi.integrate_tlsph[] + add_velocity!(du, v, particle, system) + end +end + @inline function add_velocity!(du, v, particle, system) for i in 1:ndims(system) du[i, particle] = v[i, particle] @@ -490,11 +504,15 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) return dv_ode 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; update_from_callback=false) +# Update the systems and neighborhood searches (NHS) for a simulation +# before calling `interact!` to compute forces. +# `semi` and `update_nhs_fun` are overwritten by the `SplitIntegrationCallback`. +function update_systems_and_nhs(v_ode, u_ode, semi, t; + systems=semi, update_nhs_fun=(update_nhs!), + 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 + foreach_system(systems) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -502,13 +520,13 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals end # Update NHS - @trixi_timeit timer() "update nhs" update_nhs!(semi, u_ode) + @trixi_timeit timer() "update nhs" update_nhs_fun(semi, u_ode) # Second update step. # This is used to calculate density and pressure of the fluid systems # before updating the boundary systems, # since the fluid pressure is needed by the Adami interpolation. - foreach_system(semi) do system + foreach_system(systems) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -516,7 +534,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals end # Perform correction and pressure calculation - foreach_system(semi) do system + foreach_system(systems) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -524,7 +542,7 @@ function update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=fals end # Final update step for all remaining systems - foreach_system(semi) do system + foreach_system(systems) do system v = wrap_v(v_ode, system, semi) u = wrap_u(u_ode, system, semi) @@ -546,16 +564,22 @@ function update_nhs!(semi, u_ode) end end -function add_source_terms!(dv_ode, v_ode, u_ode, semi, t) - foreach_system(semi) do system - dv = wrap_v(dv_ode, system, semi) - v = wrap_v(v_ode, system, semi) - u = wrap_u(u_ode, system, semi) +# The `SplitIntegrationCallback` overwrites `semi_wrap` to use a different semi for wrapping +# TODO `semi` is not used yet, but will be used when the source terms API is modified +# to match the custom quantities API. +function add_source_terms!(dv_ode, v_ode, u_ode, semi, t; semi_wrap=semi) + foreach_system(semi_wrap) do system + dv = wrap_v(dv_ode, system, semi_wrap) + v = wrap_v(v_ode, system, semi_wrap) + u = wrap_u(u_ode, system, semi_wrap) @threaded semi 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), t) + # Dispatch by system type to exclude boundary systems. + # `integrate_tlsph` is extracted from the `semi_wrap`, so that this function + # can be used in the `SplitIntegrationCallback` as well. + add_acceleration!(dv, particle, system, semi_wrap.integrate_tlsph[]) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system), t, + semi_wrap.integrate_tlsph[]) end end @@ -565,6 +589,15 @@ end @inline source_terms(system) = nothing @inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms +@inline function add_acceleration!(dv, particle, system, integrate_tlsph) + add_acceleration!(dv, particle, system) +end + +@inline function add_acceleration!(dv, particle, system::TotalLagrangianSPHSystem, + integrate_tlsph) + integrate_tlsph && add_acceleration!(dv, particle, system) +end + @inline add_acceleration!(dv, particle, system) = dv @inline function add_acceleration!(dv, particle, system::Union{FluidSystem, SolidSystem}) @@ -577,6 +610,17 @@ end return dv end +@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t, + integrate_tlsph) + add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) +end + +@inline function add_source_terms_inner!(dv, v, u, particle, + system::TotalLagrangianSPHSystem, + source_terms_, t, integrate_tlsph) + integrate_tlsph && add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) +end + @inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_, t) coords = current_coords(u, system, particle) velocity = current_velocity(v, system, particle) @@ -674,163 +718,176 @@ end function update_nhs!(neighborhood_search, system::FluidSystem, neighbor::Union{FluidSystem, TotalLagrangianSPHSystem}, - u_system, u_neighbor, semi) - # The current coordinates of fluids and solids change over time + u_system, u_neighbor, semi; + # The current coordinates of fluids and solids change over time + points_moving=(true, true)) update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) + semi; points_moving, eachindex_y=active_particles(neighbor)) end function update_nhs!(neighborhood_search, system::FluidSystem, neighbor::BoundarySPHSystem, - u_system, u_neighbor, semi) - # Boundary coordinates only change over time when `neighbor.ismoving[]` + u_system, u_neighbor, semi; + # Boundary coordinates only change over time when `neighbor.ismoving[]` + points_moving=(true, neighbor.ismoving[])) update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, neighbor.ismoving[])) + semi; points_moving) end function update_nhs!(neighborhood_search, system::FluidSystem, neighbor::OpenBoundarySPHSystem, - u_system, u_neighbor, semi) - # The current coordinates of fluids and open boundaries change over time. - + u_system, u_neighbor, semi; + # The current coordinates of fluids and open boundaries change over time + points_moving=(true, true)) # TODO: Update only `active_coordinates` of open boundaries. # Problem: Removing inactive particles from neighboring lists is necessary. update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) + semi; points_moving, eachindex_y=active_particles(neighbor)) end function update_nhs!(neighborhood_search, system::OpenBoundarySPHSystem, neighbor::FluidSystem, - u_system, u_neighbor, semi) - # The current coordinates of both open boundaries and fluids change over time. - + u_system, u_neighbor, semi; + # The current coordinates of both open boundaries and fluids change over time + points_moving=(true, true)) # TODO: Update only `active_coordinates` of open boundaries. # Problem: Removing inactive particles from neighboring lists is necessary. update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) + semi; points_moving, eachindex_y=active_particles(neighbor)) end function update_nhs!(neighborhood_search, system::OpenBoundarySPHSystem, neighbor::TotalLagrangianSPHSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. This NHS is never used. return neighborhood_search end function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::OpenBoundarySPHSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. This NHS is never used. return neighborhood_search end function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::FluidSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + # The current coordinates of fluids and solids change over time + points_moving=(true, true)) # The current coordinates of fluids and solids change over time update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true), eachindex_y=active_particles(neighbor)) + semi; points_moving, eachindex_y=active_particles(neighbor)) end function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::TotalLagrangianSPHSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. Neighborhood search works on the initial coordinates, which don't change. return neighborhood_search end function update_nhs!(neighborhood_search, system::TotalLagrangianSPHSystem, neighbor::BoundarySPHSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + # The current coordinates of solids change over time + points_moving=(true, neighbor.ismoving[])) # The current coordinates of solids change over time. # Boundary coordinates only change over time when `neighbor.ismoving[]`. update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, neighbor.ismoving[])) + semi; points_moving) end # This function is the same as the one below to avoid ambiguous dispatch when using `Union` function update_nhs!(neighborhood_search, system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, - neighbor::FluidSystem, u_system, u_neighbor, semi) + neighbor::FluidSystem, u_system, u_neighbor, semi; + # Boundary coordinates only change over time when `neighbor.ismoving[]`. + # The current coordinates of fluids and solids change over time. + points_moving=(system.ismoving[], true)) # Depending on the density calculator of the boundary model, this NHS is used for # - kernel summation (`SummationDensity`) # - continuity equation (`ContinuityDensity`) # - pressure extrapolation (`AdamiPressureExtrapolation`) - # - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - # The current coordinates of fluids and solids change over time. update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], true), + semi; points_moving, eachindex_y=active_particles(neighbor)) end # This function is the same as the one above to avoid ambiguous dispatch when using `Union` function update_nhs!(neighborhood_search, system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, - neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi) + neighbor::TotalLagrangianSPHSystem, u_system, u_neighbor, semi; + # Boundary coordinates only change over time when `neighbor.ismoving[]`. + # The current coordinates of solids change over time. + points_moving=(system.ismoving[], true)) # Depending on the density calculator of the boundary model, this NHS is used for # - kernel summation (`SummationDensity`) # - continuity equation (`ContinuityDensity`) # - pressure extrapolation (`AdamiPressureExtrapolation`) - # - # Boundary coordinates only change over time when `neighbor.ismoving[]`. - # The current coordinates of fluids and solids change over time. update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], true)) + semi; points_moving) end function update_nhs!(neighborhood_search, system::BoundarySPHSystem{<:BoundaryModelDummyParticles}, neighbor::BoundarySPHSystem, - u_system, u_neighbor, semi) - # `system` coordinates only change over time when `system.ismoving[]`. - # `neighbor` coordinates only change over time when `neighbor.ismoving[]`. + u_system, u_neighbor, semi; + # `system` coordinates only change over time when `system.ismoving[]`. + # `neighbor` coordinates only change over time when `neighbor.ismoving[]`. + points_moving=(system.ismoving[], neighbor.ismoving[])) update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(system.ismoving[], neighbor.ismoving[])) + semi; points_moving) end function update_nhs!(neighborhood_search, system::DEMSystem, neighbor::DEMSystem, - u_system, u_neighbor, semi) - # Both coordinates change over time + u_system, u_neighbor, semi; + # Both coordinates change over time + points_moving=(true, true)) update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, true)) + semi; points_moving) end function update_nhs!(neighborhood_search, system::DEMSystem, neighbor::BoundaryDEMSystem, - u_system, u_neighbor, semi) - # DEM coordinates change over time, the boundary coordinates don't + u_system, u_neighbor, semi; + # DEM coordinates change over time, the boundary coordinates don't + points_moving=(true, false)) update!(neighborhood_search, current_coordinates(u_system, system), current_coordinates(u_neighbor, neighbor), - semi, points_moving=(true, false)) + semi; points_moving) end function update_nhs!(neighborhood_search, system::BoundarySPHSystem, neighbor::FluidSystem, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. This NHS is never used. return neighborhood_search end @@ -838,7 +895,8 @@ end function update_nhs!(neighborhood_search, system::BoundaryDEMSystem, neighbor::Union{DEMSystem, BoundaryDEMSystem}, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. This NHS is never used. return neighborhood_search end @@ -846,13 +904,14 @@ end function update_nhs!(neighborhood_search, system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, neighbor::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, - u_system, u_neighbor, semi) + u_system, u_neighbor, semi; + points_moving=(false, false)) # Don't update. This NHS is never used. return neighborhood_search end # Forward to PointNeighbors.jl -function update!(neighborhood_search, x, y, semi; points_moving=(true, false), +function update!(neighborhood_search, x, y, semi; points_moving=(true, true), eachindex_y=axes(y, 2)) PointNeighbors.update!(neighborhood_search, x, y; points_moving, eachindex_y, parallelization_backend=semi.parallelization_backend) diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 571517799..9f651aadd 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -2,7 +2,11 @@ function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, - neighbor_system::TotalLagrangianSPHSystem, semi) + neighbor_system::TotalLagrangianSPHSystem, semi; + integrate_tlsph=semi.integrate_tlsph[]) + # Skip interaction if TLSPH systems are integrated separately + integrate_tlsph || return dv + interact_solid_solid!(dv, particle_system, neighbor_system, semi) end @@ -62,7 +66,11 @@ end function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, - neighbor_system::FluidSystem, semi) + neighbor_system::FluidSystem, semi; + integrate_tlsph=semi.integrate_tlsph[]) + # Skip interaction if TLSPH systems are integrated separately + integrate_tlsph || return dv + sound_speed = system_sound_speed(neighbor_system) system_coords = current_coordinates(u_particle_system, particle_system) @@ -164,7 +172,8 @@ end function interact!(dv, v_particle_system, u_particle_system, v_neighbor_system, u_neighbor_system, particle_system::TotalLagrangianSPHSystem, - neighbor_system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, semi) + neighbor_system::Union{BoundarySPHSystem, OpenBoundarySPHSystem}, semi; + integrate_tlsph=semi.integrate_tlsph[]) # TODO continuity equation? return dv end diff --git a/test/count_allocations.jl b/test/count_allocations.jl index eebcaa24b..0c3bd2953 100644 --- a/test/count_allocations.jl +++ b/test/count_allocations.jl @@ -15,7 +15,7 @@ function copy_semi_with_no_update_nhs(semi) for searches in semi.neighborhood_searches) return Semidiscretization(semi.systems, semi.ranges_u, semi.ranges_v, - neighborhood_searches, SerialBackend()) + neighborhood_searches, SerialBackend(), Ref(true)) end # Forward `foreach_neighbor` to wrapped neighborhood search diff --git a/test/schemes/fluid/rhs.jl b/test/schemes/fluid/rhs.jl index 21305b73c..d013d7da3 100644 --- a/test/schemes/fluid/rhs.jl +++ b/test/schemes/fluid/rhs.jl @@ -2,6 +2,7 @@ @testset verbose=true "`pressure_acceleration`" begin # Use `@trixi_testset` to isolate the mock functions in a separate namespace @trixi_testset "Symmetry" begin + TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 masses = [[0.01, 0.01], [0.73, 0.31]] densities = [ [1000.0, 1000.0], @@ -47,7 +48,6 @@ density = ones(3) state_equation = Val(:state_equation) smoothing_kernel = Val(:smoothing_kernel) - TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 smoothing_length = -1.0 fluid = InitialCondition(; coordinates, velocity, mass, density) diff --git a/test/test_util.jl b/test/test_util.jl index a487a2d49..dab317bad 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -45,10 +45,11 @@ macro trixi_testset(name, expr) end struct DummySemidiscretization - parallelization_backend::Any + parallelization_backend :: Any + integrate_tlsph :: Any function DummySemidiscretization(; parallelization_backend=SerialBackend()) - new(parallelization_backend) + new(parallelization_backend, Ref(true)) end end