From e34b76c0488c26efb2e0df3e46e057edf1d4aec1 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 29 Jan 2024 11:08:03 +0100 Subject: [PATCH] Use source terms instead of damping (#344) * Use source terms instead of damping * Update docs of `Semidiscretization` * Update docs of systems * [skip ci] Add docs for source terms * Remove `neighborhood_search=GridNeighborhoodSearch` * Fix tests * Add check for positive density to `InitialCondition` * Add tests for source terms * Fix tests * Fix examples * Make n-body fast again * Implement suggestions * Add source terms after interactions * Implement suggestions * Reformat --- examples/fluid/accelerated_tank_2d.jl | 3 +- examples/fluid/dam_break_2d.jl | 3 +- examples/fluid/dam_break_3d.jl | 3 +- examples/fluid/falling_water_column_2d.jl | 3 +- examples/fluid/moving_wall_2d.jl | 3 +- examples/fluid/periodic_channel_2d.jl | 1 - examples/fluid/rectangular_tank_2d.jl | 3 +- examples/fluid/rectangular_tank_edac_2d.jl | 3 +- examples/fsi/dam_break_2d.jl | 3 +- examples/fsi/dam_break_gate_2d.jl | 3 +- examples/fsi/falling_spheres_2d.jl | 3 +- examples/fsi/falling_water_column_2d.jl | 3 +- examples/n_body/n_body_benchmark_trixi.jl | 4 +- examples/n_body/n_body_solar_system.jl | 4 +- examples/n_body/n_body_system.jl | 11 +- examples/solid/oscillating_beam_2d.jl | 2 +- src/TrixiParticles.jl | 3 +- src/general/initial_condition.jl | 4 + src/general/semidiscretization.jl | 109 +++++++++++++----- src/neighborhood_search/trivial_nhs.jl | 5 + .../fluid/entropically_damped_sph/system.jl | 26 +++-- .../fluid/weakly_compressible_sph/system.jl | 55 ++++++--- .../solid/total_lagrangian_sph/system.jl | 23 ++-- src/setups/sphere_shape.jl | 8 +- test/general/density_calculator.jl | 3 +- test/general/initial_condition.jl | 6 + test/general/semidiscretization.jl | 46 ++++++-- .../fluid/weakly_compressible_sph/rhs.jl | 2 +- 28 files changed, 221 insertions(+), 124 deletions(-) diff --git a/examples/fluid/accelerated_tank_2d.jl b/examples/fluid/accelerated_tank_2d.jl index 6d759c919..2cc797fbe 100644 --- a/examples/fluid/accelerated_tank_2d.jl +++ b/examples/fluid/accelerated_tank_2d.jl @@ -62,8 +62,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index bef693333..7f696b324 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -65,8 +65,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=250) diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl index 8f6c17b19..8ee08f4b5 100644 --- a/examples/fluid/dam_break_3d.jl +++ b/examples/fluid/dam_break_3d.jl @@ -56,8 +56,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl index 0b4e0e6e4..5346a70d5 100644 --- a/examples/fluid/falling_water_column_2d.jl +++ b/examples/fluid/falling_water_column_2d.jl @@ -56,8 +56,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/examples/fluid/moving_wall_2d.jl b/examples/fluid/moving_wall_2d.jl index 09f2382dd..553c3e929 100644 --- a/examples/fluid/moving_wall_2d.jl +++ b/examples/fluid/moving_wall_2d.jl @@ -77,8 +77,7 @@ boundary_system_wall = BoundarySPHSystem(wall, boundary_model_wall, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_wall, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system_tank, boundary_system_wall) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 0b8f8b9de..65152f7f7 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -52,7 +52,6 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch, periodic_box_min_corner=[0.0, -0.25], periodic_box_max_corner=[1.0, 0.75]) ode = semidiscretize(semi, tspan) diff --git a/examples/fluid/rectangular_tank_2d.jl b/examples/fluid/rectangular_tank_2d.jl index 0feab0a53..37d81f1bb 100644 --- a/examples/fluid/rectangular_tank_2d.jl +++ b/examples/fluid/rectangular_tank_2d.jl @@ -63,8 +63,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) diff --git a/examples/fluid/rectangular_tank_edac_2d.jl b/examples/fluid/rectangular_tank_edac_2d.jl index 186c1493c..9b1a63746 100644 --- a/examples/fluid/rectangular_tank_edac_2d.jl +++ b/examples/fluid/rectangular_tank_edac_2d.jl @@ -49,8 +49,7 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=50) diff --git a/examples/fsi/dam_break_2d.jl b/examples/fsi/dam_break_2d.jl index 9ac3e86c8..1e886aff6 100644 --- a/examples/fsi/dam_break_2d.jl +++ b/examples/fsi/dam_break_2d.jl @@ -122,8 +122,7 @@ solid_system = TotalLagrangianSPHSystem(solid, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, solid_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system, solid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 6716aff9a..0300cd529 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -151,8 +151,7 @@ solid_system = TotalLagrangianSPHSystem(solid, # ========================================================================================== # ==== Simulation semi = Semidiscretization(fluid_system, boundary_system_tank, - boundary_system_gate, solid_system, - neighborhood_search=GridNeighborhoodSearch) + boundary_system_gate, solid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index f06dcfa1a..c6fb8f6f2 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -109,8 +109,7 @@ solid_system_2 = TotalLagrangianSPHSystem(sphere2, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_system_2) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=10) diff --git a/examples/fsi/falling_water_column_2d.jl b/examples/fsi/falling_water_column_2d.jl index 8554a2e2f..15d597504 100644 --- a/examples/fsi/falling_water_column_2d.jl +++ b/examples/fsi/falling_water_column_2d.jl @@ -64,8 +64,7 @@ solid_system = TotalLagrangianSPHSystem(solid, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(fluid_system, solid_system, - neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(fluid_system, solid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/examples/n_body/n_body_benchmark_trixi.jl b/examples/n_body/n_body_benchmark_trixi.jl index b6a5ccf8c..17868d201 100644 --- a/examples/n_body/n_body_benchmark_trixi.jl +++ b/examples/n_body/n_body_benchmark_trixi.jl @@ -68,7 +68,7 @@ masses = [ # Offset sun momentum velocity[:, 1] = -velocity[:, 2:end] * masses[2:end] / SOLAR_MASS -initial_condition = InitialCondition(; coordinates, velocity, density=0.0, mass=masses) +initial_condition = InitialCondition(; coordinates, velocity, density=1.0, mass=masses) G = 1.0 particle_system = NBodySystem(initial_condition, G) @@ -76,7 +76,7 @@ particle_system = NBodySystem(initial_condition, G) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(particle_system) +semi = Semidiscretization(particle_system, neighborhood_search=nothing) # This is significantly faster than using OrdinaryDiffEq. function symplectic_euler!(velocity, coordinates, semi) diff --git a/examples/n_body/n_body_solar_system.jl b/examples/n_body/n_body_solar_system.jl index 72f4f178d..6fa860089 100644 --- a/examples/n_body/n_body_solar_system.jl +++ b/examples/n_body/n_body_solar_system.jl @@ -23,7 +23,7 @@ masses = [ 1.99e30, 3.30e23, 4.87e24, 5.97e24, 6.42e23, 1.90e27, 5.68e26, 8.68e25, 1.02e26, ] -initial_condition = InitialCondition(; coordinates, velocity, density=0.0, mass=masses) +initial_condition = InitialCondition(; coordinates, velocity, density=1.0, mass=masses) G = 6.6743e-11 particle_system = NBodySystem(initial_condition, G) @@ -31,7 +31,7 @@ particle_system = NBodySystem(initial_condition, G) # ========================================================================================== # ==== Simulation -semi = Semidiscretization(particle_system) +semi = Semidiscretization(particle_system, neighborhood_search=nothing) day = 24 * 3600.0 year = 365day diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl index 9a37e03ce..068eae953 100644 --- a/examples/n_body/n_body_system.jl +++ b/examples/n_body/n_body_system.jl @@ -18,10 +18,6 @@ TrixiParticles.timer_name(::NBodySystem) = "nbody" @inline Base.eltype(system::NBodySystem) = eltype(system.initial_condition.coordinates) -@inline function TrixiParticles.add_acceleration!(dv, particle, system::NBodySystem) - return dv -end - function TrixiParticles.write_u0!(u0, system::NBodySystem) u0 .= system.initial_condition.coordinates @@ -53,11 +49,8 @@ function TrixiParticles.interact!(dv, v_particle_system, u_particle_system, neighbor_system::NBodySystem) (; mass, G) = neighbor_system - system_coords = TrixiParticles.current_coordinates(u_particle_system, - particle_system) - - neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, - neighbor_system) + system_coords = TrixiParticles.current_coordinates(u_particle_system, particle_system) + neighbor_coords = TrixiParticles.current_coordinates(u_neighbor_system, neighbor_system) # Loop over all pairs of particles and neighbors within the kernel cutoff. TrixiParticles.for_particle_neighbor(particle_system, neighbor_system, diff --git a/examples/solid/oscillating_beam_2d.jl b/examples/solid/oscillating_beam_2d.jl index 56b09b011..3791197ff 100644 --- a/examples/solid/oscillating_beam_2d.jl +++ b/examples/solid/oscillating_beam_2d.jl @@ -58,7 +58,7 @@ solid_system = TotalLagrangianSPHSystem(solid, # ========================================================================================== # ==== Simulation -semi = Semidiscretization(solid_system, neighborhood_search=GridNeighborhoodSearch) +semi = Semidiscretization(solid_system) ode = semidiscretize(semi, tspan) info_callback = InfoCallback(interval=100) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 434e4ef15..f232a664b 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -50,11 +50,12 @@ export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerr export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, PressureMirroring, PressureZeroing export BoundaryMovement -export GridNeighborhoodSearch +export GridNeighborhoodSearch, TrivialNeighborhoodSearch export examples_dir, trixi_include export trixi2vtk export RectangularTank, RectangularShape, SphereShape export VoxelSphere, RoundSphere, reset_wall! +export SourceTermDamping export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection export nparticles diff --git a/src/general/initial_condition.jl b/src/general/initial_condition.jl index 764b82beb..53496603a 100644 --- a/src/general/initial_condition.jl +++ b/src/general/initial_condition.jl @@ -126,6 +126,10 @@ struct InitialCondition{ELTYPE} densities = density_fun.(coordinates_svector) end + if any(densities .< eps()) + throw(ArgumentError("density must be positive and larger than `eps()`")) + end + if pressure isa AbstractVector if length(pressure) != n_particles throw(ArgumentError("Expected: length(pressure) == size(coordinates, 2)\n" * diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 0c816add3..b79ea3044 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -1,27 +1,42 @@ """ - Semidiscretization(systems...; neighborhood_search=nothing, damping_coefficient=nothing) + Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch, + periodic_box_min_corner=nothing, periodic_box_max_corner=nothing) The semidiscretization couples the passed systems to one simulation. The type of neighborhood search to be used in the simulation can be specified with the keyword argument `neighborhood_search`. A value of `nothing` means no neighborhood search. +# Arguments +- `systems`: Systems to be coupled in this semidiscretization + +# Keywords +- `neighborhood_search`: The type of neighborhood search to be used in the simulation. + By default, the [`GridNeighborhoodSearch`](@ref) is used. + Use [`TrivialNeighborhoodSearch`](@ref) to loop over all particles + (no neighborhood search). +- `periodic_box_min_corner`: In order to use a periodic domain, pass the coordinates + of the domain corner in negative coordinate directions. +- `periodic_box_max_corner`: In order to use a periodic domain, pass the coordinates + of the domain corner in positive coordinate directions. + # Examples ```julia -semi = Semidiscretization(fluid_system, boundary_system; neighborhood_search=GridNeighborhoodSearch, damping_coefficient=nothing) +semi = Semidiscretization(fluid_system, boundary_system) + +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=TrivialNeighborhoodSearch) ``` """ -struct Semidiscretization{S, RU, RV, NS, DC} +struct Semidiscretization{S, RU, RV, NS} systems :: S ranges_u :: RU ranges_v :: RV neighborhood_searches :: NS - damping_coefficient :: DC - function Semidiscretization(systems...; neighborhood_search=nothing, + function Semidiscretization(systems...; neighborhood_search=GridNeighborhoodSearch, periodic_box_min_corner=nothing, - periodic_box_max_corner=nothing, - damping_coefficient=nothing) + periodic_box_max_corner=nothing) sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems] ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) @@ -44,9 +59,8 @@ struct Semidiscretization{S, RU, RV, NS, DC} for neighbor in systems) for system in systems) - new{typeof(systems), typeof(ranges_u), typeof(ranges_v), - typeof(searches), typeof(damping_coefficient)}(systems, ranges_u, ranges_v, - searches, damping_coefficient) + new{typeof(systems), typeof(ranges_u), + typeof(ranges_v), typeof(searches)}(systems, ranges_u, ranges_v, searches) end end @@ -75,7 +89,6 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) summary_line(io, "#systems", length(semi.systems)) summary_line(io, "neighborhood search", semi.neighborhood_searches |> eltype |> eltype |> nameof) - summary_line(io, "damping coefficient", semi.damping_coefficient) summary_line(io, "total #particles", sum(nparticles.(semi.systems))) summary_footer(io) end @@ -310,11 +323,10 @@ function kick!(dv_ode, v_ode, u_ode, semi, t) @trixi_timeit timer() "update systems and nhs" update_systems_and_nhs(v_ode, u_ode, semi, t) - @trixi_timeit timer() "gravity and damping" gravity_and_damping!(dv_ode, v_ode, - semi) - @trixi_timeit timer() "system interaction" system_interaction!(dv_ode, v_ode, u_ode, semi) + + @trixi_timeit timer() "source terms" add_source_terms!(dv_ode, v_ode, u_ode, semi) end return dv_ode @@ -374,25 +386,28 @@ function update_nhs(u_ode, semi) end end -function gravity_and_damping!(dv_ode, v_ode, semi) - (; damping_coefficient) = semi - - # Set velocity and add acceleration for each system +function add_source_terms!(dv_ode, v_ode, u_ode, semi) 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) @threaded for particle in each_moving_particle(system) - # This can be dispatched per system + # Dispatch by system type to exclude boundary systems add_acceleration!(dv, particle, system) - add_damping_force!(dv, damping_coefficient, v, particle, system) + add_source_terms_inner!(dv, v, u, particle, system, source_terms(system)) end end return dv_ode end -@inline function add_acceleration!(dv, particle, system) +@inline source_terms(system) = nothing +@inline source_terms(system::Union{FluidSystem, SolidSystem}) = system.source_terms + +@inline add_acceleration!(dv, particle, system) = dv + +@inline function add_acceleration!(dv, particle, system::Union{FluidSystem, SolidSystem}) (; acceleration) = system for i in 1:ndims(system) @@ -402,20 +417,56 @@ end return dv end -@inline add_acceleration!(dv, particle, system::BoundarySPHSystem) = dv +@inline function add_source_terms_inner!(dv, v, u, particle, system, source_terms_) + coords = current_coords(u, system, particle) + velocity = current_velocity(v, system, particle) + density = particle_density(v, system, particle) + pressure = particle_pressure(v, system, particle) -@inline function add_damping_force!(dv, damping_coefficient, v, particle, - system::FluidSystem) - for i in 1:ndims(system) - dv[i, particle] -= damping_coefficient * v[i, particle] + source = source_terms_(coords, velocity, density, pressure) + + # Loop over `eachindex(source)`, so that users could also pass source terms for + # the density when using `ContinuityDensity`. + for i in eachindex(source) + dv[i, particle] += source[i] end return dv end -# Currently no damping for non-fluid systems -@inline add_damping_force!(dv, damping_coefficient, v, particle, system) = dv -@inline add_damping_force!(dv, ::Nothing, v, particle, system::FluidSystem) = dv +@inline add_source_terms_inner!(dv, v, u, particle, system, source_terms_::Nothing) = dv + +@doc raw""" + SourceTermDamping(; damping_coefficient) + +A source term to be used when a damping step is required before running a full simulation. +The term ``-c \cdot v_a`` is added to the acceleration ``\frac{\mathrm{d}v_a}{\mathrm{d}t}`` +of particle ``a``, where ``c`` is the damping coefficient and ``v_a`` is the velocity of +particle ``a``. + +# Keywords +- `damping_coefficient`: The coefficient ``d`` above. A higher coefficient means more + damping. A coefficient of `1e-4` is a good starting point for + damping a fluid at rest. + +# Examples +```julia +source_terms = SourceTermDamping(; damping_coefficient=1e-4) +``` +""" +struct SourceTermDamping{ELTYPE} + damping_coefficient::ELTYPE + + function SourceTermDamping(; damping_coefficient) + return new{typeof(damping_coefficient)}(damping_coefficient) + end +end + +@inline function (source_term::SourceTermDamping)(coords, velocity, density, pressure) + (; damping_coefficient) = source_term + + return -damping_coefficient * velocity +end function system_interaction!(dv_ode, v_ode, u_ode, semi) # Call `interact!` for each pair of systems diff --git a/src/neighborhood_search/trivial_nhs.jl b/src/neighborhood_search/trivial_nhs.jl index 95d7ef0ea..c915cfd56 100644 --- a/src/neighborhood_search/trivial_nhs.jl +++ b/src/neighborhood_search/trivial_nhs.jl @@ -1,3 +1,8 @@ +@doc raw""" + TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle) + +Trivial neighborhood search that simply loops over all particles. +""" struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} search_radius :: ELTYPE eachparticle :: EP diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 56dfb7eb1..12dd2f8cd 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -1,7 +1,10 @@ @doc raw""" - EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, - sound_speed; alpha=0.5, viscosity=NoViscosity(), - acceleration=ntuple(_ -> 0.0, NDIMS)) + EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, sound_speed; + alpha=0.5, viscosity=NoViscosity(), + initial_pressure_function=nothing, + acceleration=ntuple(_ -> 0.0, NDIMS), + source_terms=nothing) Entropically damped artificial compressibility (EDAC) for SPH introduced by (Ramachandran 2019). As opposed to the weakly compressible SPH scheme, which uses an equation of state @@ -31,7 +34,8 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). In: Computers and Fluids 179 (2019), pages 579-594. [doi: 10.1016/j.compfluid.2018.11.023](https://doi.org/10.1016/j.compfluid.2018.11.023) """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V} <: FluidSystem{NDIMS} +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V, ST} <: + FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] density :: Array{ELTYPE, 1} # [particle] @@ -42,12 +46,14 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V} <: FluidSyst viscosity :: V nu_edac :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} + source_terms :: ST function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, smoothing_length, sound_speed; alpha=0.5, viscosity=NoViscosity(), acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel))) + ndims(smoothing_kernel)), + source_terms=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) @@ -58,7 +64,6 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V} <: FluidSyst throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) end - # Make acceleration an SVector acceleration_ = SVector(acceleration...) if length(acceleration_) != NDIMS throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) @@ -68,10 +73,11 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, DC, K, V} <: FluidSyst density_calculator = SummationDensity() - new{NDIMS, ELTYPE, typeof(density_calculator), typeof(smoothing_kernel), - typeof(viscosity)}(initial_condition, mass, density, density_calculator, - smoothing_kernel, smoothing_length, sound_speed, viscosity, - nu_edac, acceleration_) + new{NDIMS, ELTYPE, typeof(density_calculator), + typeof(smoothing_kernel), typeof(viscosity), + typeof(source_terms)}(initial_condition, mass, density, density_calculator, + smoothing_kernel, smoothing_length, sound_speed, + viscosity, nu_edac, acceleration_, source_terms) end end diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index 65614236f..5ceb2e206 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -2,8 +2,9 @@ WeaklyCompressibleSPHSystem(initial_condition, density_calculator, state_equation, smoothing_kernel, smoothing_length; - viscosity=NoViscosity(), - acceleration=ntuple(_ -> 0.0, NDIMS)) + viscosity=NoViscosity(), density_diffusion=nothing, + acceleration=ntuple(_ -> 0.0, NDIMS), + correction=nothing, source_terms=nothing) Weakly compressible SPH introduced by (Monaghan, 1994). This formulation relies on a stiff equation of state (see [`StateEquationCole`](@ref)) that generates large pressure changes @@ -12,21 +13,38 @@ see [`ContinuityDensity`](@ref) and [`SummationDensity`](@ref). # Arguments - `initial_condition`: Initial condition representing the system's particles. -- `density_calculator`: Density calculator for the SPH system. See [`ContinuityDensity`](@ref) and [`SummationDensity`](@ref). -- `state_equation`: Equation of state for the SPH system. See [`StateEquationCole`](@ref). +- `density_calculator`: Density calculator for the system. + See [`ContinuityDensity`](@ref) and [`SummationDensity`](@ref). +- `state_equation`: Equation of state for the system. See [`StateEquationCole`](@ref). +- `smoothing_kernel`: Smoothing kernel to be used for this system. + See [`SmoothingKernel`](@ref). +- `smoothing_length`: Smoothing length to be used for this system. + See [`SmoothingKernel`](@ref). # Keyword Arguments -- `viscosity`: Viscosity model for the SPH system (default: no viscosity). See [`ArtificialViscosityMonaghan`](@ref) or [`ViscosityAdami`](@ref). -- `acceleration`: Acceleration vector for the SPH system. (default: zero vector) -- `correction`: Correction method used for this SPH system. (default: no correction) +- `viscosity`: Viscosity model for this system (default: no viscosity). + 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) +- `correction`: Correction method used for this system. (default: no correction) +- `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` + or `SVector` that is to be added to the acceleration of that particle. + See, for example, [`SourceTermDamping`](@ref). + Note that these source terms will not be used in the calculation of the + boundary pressure when using a boundary with + [`BoundaryModelDummyParticles`](@ref) and [`AdamiPressureExtrapolation`](@ref). + The keyword argument `acceleration` should be used instead for + gravity-like source terms. ## References: - Joseph J. Monaghan. "Simulating Free Surface Flows in SPH". In: Journal of Computational Physics 110 (1994), pages 399-406. [doi: 10.1006/jcph.1994.1034](https://doi.org/10.1006/jcph.1994.1034) """ -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, DD, COR, PF, C} <: - FluidSystem{NDIMS} +struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, + V, DD, COR, PF, ST, C} <: FluidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} mass :: Array{ELTYPE, 1} # [particle] pressure :: Array{ELTYPE, 1} # [particle] @@ -39,6 +57,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, DD, COR, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF + source_terms :: ST cache :: C function WeaklyCompressibleSPHSystem(initial_condition, @@ -48,7 +67,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, DD, COR, viscosity=NoViscosity(), density_diffusion=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - correction=nothing) + correction=nothing, source_terms=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -81,13 +100,17 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, DC, SE, K, V, DD, COR, create_cache_wcsph(correction, initial_condition.density, NDIMS, n_particles)..., cache...) - return new{NDIMS, ELTYPE, typeof(density_calculator), typeof(state_equation), - typeof(smoothing_kernel), typeof(viscosity), typeof(density_diffusion), + return new{NDIMS, ELTYPE, typeof(density_calculator), + typeof(state_equation), typeof(smoothing_kernel), + typeof(viscosity), typeof(density_diffusion), typeof(correction), typeof(pressure_acceleration), - typeof(cache)}(initial_condition, mass, pressure, density_calculator, - state_equation, smoothing_kernel, smoothing_length, - acceleration_, viscosity, density_diffusion, correction, - pressure_acceleration, cache) + typeof(source_terms), typeof(cache)}(initial_condition, mass, pressure, + density_calculator, state_equation, + smoothing_kernel, smoothing_length, + acceleration_, viscosity, + density_diffusion, correction, + pressure_acceleration, + source_terms, cache) end end diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index d8f0d9fe0..85d1cc45c 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -4,7 +4,7 @@ young_modulus, poisson_ratio, boundary_model; n_fixed_particles=0, acceleration=ntuple(_ -> 0.0, NDIMS), - penalty_force=nothing) + penalty_force=nothing, source_terms=nothing) System for particles of an elastic solid. @@ -75,7 +75,7 @@ The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGa In: International Journal for Numerical Methods in Engineering 48 (2000), pages 1359–1400. [doi: 10.1002/1097-0207](https://doi.org/10.1002/1097-0207) """ -struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem{NDIMS} +struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF, ST} <: SolidSystem{NDIMS} initial_condition :: InitialCondition{ELTYPE} initial_coordinates :: Array{ELTYPE, 2} # [dimension, particle] current_coordinates :: Array{ELTYPE, 2} # [dimension, particle] @@ -94,13 +94,15 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem acceleration :: SVector{NDIMS, ELTYPE} boundary_model :: BM penalty_force :: PF + source_terms :: ST + function TotalLagrangianSPHSystem(initial_condition, smoothing_kernel, smoothing_length, young_modulus, poisson_ratio, boundary_model; n_fixed_particles=0, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - penalty_force=nothing) + penalty_force=nothing, source_terms=nothing) NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -131,13 +133,14 @@ struct TotalLagrangianSPHSystem{BM, NDIMS, ELTYPE <: Real, K, PF} <: SolidSystem return new{typeof(boundary_model), NDIMS, ELTYPE, typeof(smoothing_kernel), - typeof(penalty_force)}(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) + 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) end end diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index 573644be4..12695b423 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -86,10 +86,6 @@ function SphereShape(particle_spacing, radius, center_position, density; throw(ArgumentError("`particle_spacing` needs to be positive and larger than $(eps())")) end - if density < eps() - throw(ArgumentError("`density` needs to be positive and larger than $(eps())")) - end - NDIMS = length(center_position) coordinates = sphere_shape_coords(sphere_type, particle_spacing, radius, @@ -306,7 +302,7 @@ function round_sphere(sphere, particle_spacing, radius, center::SVector{3}) n_particles = round(Int, 4pi * radius^2 / particle_spacing^2) end - # With less than 5 particles, this doesn't work properly + # With fewer than 5 particles, this doesn't work properly if n_particles < 5 if n_particles == 4 # Return tetrahedron @@ -342,7 +338,7 @@ function round_sphere(sphere, particle_spacing, radius, center::SVector{3}) # In: Electronic Transactions on Numerical Analysis 25 (2006), pages 309-327. # [http://eudml.org/doc/129860](http://eudml.org/doc/129860). - # This is the Θ function, which is only defined by Leopardi as the inverse of V, without + # This is the Θ function, which is defined by Leopardi only as the inverse of V, without # giving a closed formula. theta(v) = acos(1 - v / 2pi) diff --git a/test/general/density_calculator.jl b/test/general/density_calculator.jl index a8e3feec0..2ded39476 100644 --- a/test/general/density_calculator.jl +++ b/test/general/density_calculator.jl @@ -23,8 +23,7 @@ using OrdinaryDiffEq (; cache) = fluid_system (; density) = cache # Density is in the cache for SummationDensity - semi = Semidiscretization(fluid_system, neighborhood_search=GridNeighborhoodSearch, - damping_coefficient=1e-5) + semi = Semidiscretization(fluid_system) tspan = (0.0, 0.01) ode = semidiscretize(semi, tspan) diff --git a/test/general/initial_condition.jl b/test/general/initial_condition.jl index de1688fc1..468ebe65b 100644 --- a/test/general/initial_condition.jl +++ b/test/general/initial_condition.jl @@ -63,6 +63,12 @@ @test_throws ArgumentError(error_str) InitialCondition(coordinates=zeros(2, 2), velocity=zeros(2, 2), density=ones(2)) + + error_str = "density must be positive and larger than `eps()`" + @test_throws ArgumentError(error_str) InitialCondition(coordinates=zeros(2, 2), + velocity=zeros(2, 2), + density=zeros(2), + mass=ones(2)) end @testset "Constant Quantities" begin diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index dbd15dda3..d828425bb 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -7,9 +7,6 @@ system1 = System1() system2 = System2() - Base.ndims(::System1) = 2 - Base.ndims(::System2) = 2 - TrixiParticles.u_nvariables(::System1) = 3 TrixiParticles.u_nvariables(::System2) = 4 TrixiParticles.v_nvariables(::System1) = 3 @@ -23,21 +20,21 @@ TrixiParticles.compact_support(::System2, neighbor) = 0.2 @testset verbose=true "Constructor" begin - semi = Semidiscretization(system1, system2) + semi = Semidiscretization(system1, system2, neighborhood_search=nothing) # Verification @test semi.ranges_u == (1:6, 7:18) @test semi.ranges_v == (1:6, 7:12) - nhs = ((TrixiParticles.TrivialNeighborhoodSearch{2}(0.2, Base.OneTo(2)), - TrixiParticles.TrivialNeighborhoodSearch{2}(0.2, Base.OneTo(3))), - (TrixiParticles.TrivialNeighborhoodSearch{2}(0.2, Base.OneTo(2)), - TrixiParticles.TrivialNeighborhoodSearch{2}(0.2, Base.OneTo(3)))) + nhs = ((TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(2)), + TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(3))), + (TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(2)), + TrixiParticles.TrivialNeighborhoodSearch{3}(0.2, Base.OneTo(3)))) @test semi.neighborhood_searches == nhs end - @testset verbose=true "show" begin - semi = Semidiscretization(system1, system2) + @testset verbose=true "`show`" begin + semi = Semidiscretization(system1, system2, neighborhood_search=nothing) show_compact = "Semidiscretization($System1(), $System2(), neighborhood_search=TrivialNeighborhoodSearch)" @test repr(semi) == show_compact @@ -46,12 +43,37 @@ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ │ Semidiscretization │ │ ══════════════════ │ - │ #spatial dimensions: ………………………… 2 │ + │ #spatial dimensions: ………………………… 3 │ │ #systems: ……………………………………………………… 2 │ │ neighborhood search: ………………………… TrivialNeighborhoodSearch │ - │ damping coefficient: ………………………… nothing │ │ total #particles: ………………………………… 5 │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", semi) == show_box end + + @testset verbose=true "Source Terms" begin + TrixiParticles.source_terms(::System1) = SourceTermDamping(damping_coefficient=0.1) + TrixiParticles.particle_density(v, system::System1, particle) = 0.0 + TrixiParticles.particle_pressure(v, system::System1, particle) = 0.0 + + semi = Semidiscretization(system1, system2, neighborhood_search=nothing) + + dv_ode = zeros(3 * 2 + 2 * 3) + du_ode = zeros(3 * 2 + 4 * 3) + u_ode = zero(du_ode) + + v1 = [1.0 2.0 + 3.0 4.0 + 5.0 6.0] + v2 = zeros(4 * 3) + v_ode = vcat(vec(v1), v2) + + TrixiParticles.add_source_terms!(dv_ode, v_ode, u_ode, semi) + + dv1 = TrixiParticles.wrap_v(dv_ode, system1, semi) + @test dv1 == -0.1 * v1 + + dv2 = TrixiParticles.wrap_v(dv_ode, system2, semi) + @test iszero(dv2) + end end diff --git a/test/schemes/fluid/weakly_compressible_sph/rhs.jl b/test/schemes/fluid/weakly_compressible_sph/rhs.jl index b3963e8c8..cd580c2a0 100644 --- a/test/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/test/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -32,7 +32,7 @@ coordinates = zeros(2, 3) velocity = zeros(2, 3) mass = zeros(3) - density = zeros(3) + density = ones(3) state_equation = Val(:state_equation) smoothing_kernel = Val(:smoothing_kernel) TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2