diff --git a/NEWS.md b/NEWS.md index 601e39d313..7add08c4ba 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ 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. +## Version 0.4.1 + +### Features + +- Added an incompressible SPH solver `ImplicitIncompressibleSPHSystem` (#751). + + + ## Version 0.4 ### API Changes diff --git a/docs/make.jl b/docs/make.jl index bfb906c06c..7bce8acc3c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -149,7 +149,9 @@ makedocs(sitename="TrixiParticles.jl", "Weakly Compressible SPH (Fluid)" => joinpath("systems", "weakly_compressible_sph.md"), "Entropically Damped Artificial Compressibility for SPH (Fluid)" => joinpath("systems", - "entropically_damped_sph.md") + "entropically_damped_sph.md"), + "Implicit Incompressible SPH (Fluid)" => joinpath("systems", + "implicit_incompressible_sph.md") ], "Discrete Element Method (Solid)" => joinpath("systems", "dem.md"), "Total Lagrangian SPH (Elastic Structure)" => joinpath("systems", diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 3b853a5250..5fd53aad01 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -361,6 +361,16 @@ @Article{Hormann2001 publisher = {Elsevier BV}, } +@article{Ihmsen2013, + author = {Ihmsen, Markus and Cornelis, Jens and Solenthaler, Barbara and Horvath, Christopher and Teschner, Matthias}, + year = {2013}, + month = {07}, + pages = {1--11}, + title = {Implicit incompressible {SPH}}, + volume = {20}, + journal = {IEEE transactions on visualization and computer graphics}, + doi = {10.1109/TVCG.2013.105} +} @Article{Jacobson2013, author = {Jacobson, Alec and Kavan, Ladislav and Sorkine-Hornung, Olga}, diff --git a/docs/src/systems/implicit_incompressible_sph.md b/docs/src/systems/implicit_incompressible_sph.md new file mode 100644 index 0000000000..082f42de9c --- /dev/null +++ b/docs/src/systems/implicit_incompressible_sph.md @@ -0,0 +1,295 @@ +# [Implicit Incompressible SPH](@id iisph) + +Implicit Incompressible SPH (IISPH) as introduced by [Ihmsen et al. (2013)](@cite Ihmsen2013) +is a method that achieves incompressibility by solving the pressure Poisson equation. +The resulting linear system is iteratively solved with the relaxed Jacobi method. +Unlike the [weakly compressible SPH method](@ref wcsph), incompressible methods determine +pressure by enforcing the incompressibility constraint rather than using an equation of +state. + +```@autodocs +Modules = [TrixiParticles] +Pages = [joinpath("schemes", "fluid", "implicit_incompressible_sph", "system.jl")] +``` +## Derivation +To derive the linear system of the pressure Poisson equation, we start by discretizing the +continuity equation +```math +\frac{D\rho}{Dt} = - \rho \nabla \cdot \bm{v}. +``` + +For a particle ``i``, discretizing the left-hand side of the equation with a forward +difference yields + +```math +\frac{\rho_i(t+ \Delta t) - \rho_i(t)}{\Delta t}. +``` + +The divergence in the right-hand side is discretized with the SPH discretization for +particle ``i`` as +```math +-\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \nabla W_{ij}, +``` +where ``\bm{v}_{ij} = \bm{v}_i - \bm{v}_j``. + +Together, the following discretized version of the continuity equation for a particle ``i`` +is achieved: +```math +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij}. +``` + +The right-hand side contains the unknown velocities ``\bm{v}_{ij}(t + \Delta t)`` at the +next time step, making it an implicit formulation. + +Using the semi-implicit Euler method, we can obtain the velocity in the next time step as +```math +\bm{v}_i(t + \Delta t) = \bm{v}_i(t) + \Delta t \frac{\bm{F}_i^\text{adv}(t) + \bm{F}_i^p(t)}{m_i}, +``` + +where ``\bm{F}_i^{\text{adv}}`` denotes all non-pressure forces such as gravity, viscosity, surface +tension and more, while ``\bm{F}_i^p``denotes the unknown pressure forces, which we +want to solve for. + +Note that the IISPH is an incompressible method, which means that the density of the +fluid remain constant over time. By assuming a fixed reference density ``\rho_0`` for all +fluid particle over the whole time of the simulation, the density value at the next time +step ``\rho_i(t + \Delta t)`` also has to be this rest density. So ``\rho_0`` can be plugged in +for ``\rho_i(t + \Delta t)`` in the formula above. + +The goal is to compute the pressure values required to obtain the pressure acceleration that +ensures each particle reaches the rest density in the next time step. At the moment these +pressure values are unknown in ``t``, but all the non-pressure forces are already known. +Therefore a predicted velocity can be calculated, which depends only on the non-pressure +forces ``\bm{F}_i^{\text{adv}}``: + +```math +\bm{v}_i^{\text{adv}}(t+\Delta t)= \bm{v}_i(t) + \Delta t \frac{\bm{F}_i^{\text{adv}}(t)}{m_i}, +``` +Using this predicted velocity and the continuity equation, a predicted density can be defined +in a similar way as + +```math +\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t). +``` + +To achieve the rest density, the unknown pressure forces must counteract the compression +caused by the non-pressure forces. In other words, they must resolve the deviation between +the predicted density and the reference density. + +Therefore, the following equation needs to be fulfilled: +```math +\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. +``` + +This expression is derived by substituting the reference density ``\rho_0`` for +``\rho_i(t+\Delta t)`` in the discretized continuity equation and inserting the definitions of +the predicted velocity ``\bm{v}_{ij}(t+\Delta t)`` and predicted density ``\rho_i^{\text{adv}}``. + +The pressure force is defined as: + +```math +\bm{F}_i^p(t) = -\frac{m_i}{\rho_i} \nabla p_i = -m_i \sum_j m_j \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_j(t)}{\rho_j^2(t)} \right) \nabla W_{ij}. +``` + +Substituting this definition into the equation yields a linear system ``\bm{A}(t) \bm{p}(t) = \bm{b}(t)`` +with one equation and one unknown pressure value for each particle. +For ``n`` particles, this is a linear system with ``n`` equations and ``n`` unknowns: + +```math +\sum_j a_{ij} p_j = b_i = \rho_0 - \rho_i^{\text{adv}}. +``` + +To solve for the pressure values, a relaxed Jacobi scheme is used. +This is an iterative numerical method to solve a linear system ``Ap=b``. In each iteration, the +new values of ``p`` are obtained as computed as + +```math +p_i^{(k+1)} = (1-\omega) p_i^{(k)} + \omega \left( \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} p_j^{(k)} \right)\right). +``` + +Substituting the right-hand side, we obtain + +```math +p_i^{l+1} = (1-\omega) p_i^l + \omega \frac{\rho_0 - \rho_i^{\text{adv}} - \sum_{j \neq i} a_{ij}p_j^l}{a_{ii}}. +``` + +Therefore the diagonal elements ``a_{ii}`` and the sum ``\sum_{j \neq i} a_{ij}p_j^l`` need to +be determined. +This can be done efficiently by separating the pressure force acceleration into two +components: One that describes the influence of particle ``i``'s own pressure on its +displacement, and one that describes the displacement due to the pressure values of the +neighboring particles ``j``. + +The pressure acceleration is given by: +```math +\Delta t^2 \frac{\bm{F}_i^p}{m_i} = -\Delta t^2 \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right)\nabla W_{ij} = \underbrace{\left(- \Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij} \right) }_{d_{ii}} p_i + \sum_j \underbrace{-\Delta t^2 \frac{m_j}{\rho_j^2} \nabla W_{ij}}_{d_{ij}}p_j. +``` + +The ``d_{ii}p_i`` value describes the displacement of particle ``i`` because of the particle ``i`` +and ``d_{ij}p_j`` describes the influence from the neighboring particles ``j``. +Using this new values the linear system can be rewritten as + +```math +\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \nabla W_{ij}, +``` + +where the first sum over ``k`` loops over all neighbor particles of ``i`` and +the second sum over ``k`` loops over the neighbor particles of ``j`` +(which is a neighbor of ``i``). +So the sum over the neighboring pressure values ``p_j`` also includes the pressure values +``p_i``, since ``i`` is a neighbor of ``j``. +To separate this sum, it can be written as + +```math +\sum_k d_{jk} p_k = \sum_{k \neq i} d_{jk} p_k + d_{ji} p_i. +``` + +With this separation, the equation for the linear system can again be rewritten as + +```math +\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \nabla W_{ij}. +``` + +In this formulation all coefficients that are getting multiplied with the pressure value ``p_i`` +are separated from the other. The diagonal elements ``a_{ii}`` can therefore be defined as: + +```math +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij}. +``` + +The remaining part of the equation represents the influence of the other pressure values ``p_j``. +​Hence, the final relaxed Jacobi iteration takes the form: + +```math +p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right). +``` + +Because interactions are local, limited to particles within the kernel support defined by +the smoothing length, each particle's pressure ``p_i`` depends only on its own value and +nearby particles. +Consequently, the matrix ``A`` is sparse (most entries are zero). + +The diagonal elements ``a_{ii}`` get computed and stored at the beginning of the simulation step +and remain unchanged throughout the relaxed Jacobi iterations. The same applies for the +``d_{ii}``. The coefficients ``d_{ij}`` are computed, but not stored, as part of the calculation +of ``a_{ii}``. +During the Jacobi iterations, two loops over all particles are performed: +The first updates the values of ``\sum_j d_{ij}`` and the second computes the updated +pressure values ``p_i^{l+1}``. + +The final pressure update follows the equation of the relaxed Jacobi scheme as stated above, with two +exceptions: +#### 1. Pressure clamping +To avoid negative pressure values, one can enable pressure clamping. In this case, the +pressure update is given by + +```math +p_i^{l+1} = \max \left(0, (1-\omega) p_i^l + \omega \frac{\rho_0 - \rho_i^{\text{adv}} - \sum_{j \neq i} a_{ij}p_j^l}{a_{ii}}\right). +``` + +#### 2. Small diagonal elements +If the diagonal element ``a_{ii}`` becomes too small or even zero, the simulation may become +unstable. This can occur, for example, if a particle is isolated and receives little or no +influence from neighboring particles, or due to numerical cancellations. In such cases, the +updated pressure value is set to zero. + +There are also other options, like setting ``a_{ii}`` to the threshold value if it is beneath +and then update with the usual equation, or just don't update the pressure value at all, and +continue with the old value. The numerical error introduced by this technique remains small, +as only isolated or almost isolated particles are affected. + + +## Boundary Handling + +The previously introduced formulation only considered interactions between fluid particles +and neglected interactions between fluid and boundary particles. To account for boundary +interactions, a few modifications to the previous equations are required. + +First, the discretized form of the continuity equation must be adapted for the case in which +a neighboring particle is a boundary particle. From now on, we distinguish between +neighboring fluid particles (indexed by ``j``) and neighboring boundary particles (indexed by +``b``). + +The updated discretized continuity equation becomes: + +```math +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \nabla W_{ib}. +``` + +Since boundary particles have zero velocity, the difference between the fluid +particle's velocity and the boundary particle's velocity simplifies to just the fluid +particle's velocity ``\bm{v}_{ib}(t+\Delta t) = \bm{v}_{i}(t+\Delta t)``. +Accordingly, the predicted density ``\rho^{\text{adv}}`` becomes: + +```math +\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \nabla W_{ib}(t). +``` + +This leads to the following updated formulation of the linear system: + +```math +\Delta t^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. +``` + +Note that, since boundary particles are fixed, the force ``F_b^p`` is zero and does not appear in this equation. + +The pressure force acting on a fluid particle is computed as: + +```math +\bm{F}_i^p(t) = -\sum_j m_j \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_j(t)}{\rho_j^2(t)} \right) \nabla W_{ij}(t) - \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). +``` + +This also leads to an updated version of the equation for the diagonal elements: +```math +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_b m_b (-d_{bi}) \nabla W_{ib}. +``` + +From this point forward, the computation of the coefficients required for the Jacobi scheme +(such as ``d_{ii}``, ``d_{ij}`` etc.) depends on the specific boundary density evaluation method +used in the chosen boundary model. + + + +### Pressure Mirroring +When using pressure mirroring, the pressure value ``p_b`` of a boundary particle in the equation +above is defined to be equal to the pressure of the corresponding fluid particle ``p_i``. +In other words, the boundary particle "mirrors" the pressure of the fluid particle interacting +with it. As a result, the coefficient that describes the influence of a particle's own +pressure value ``p_i`` ​must also include contributions from boundary particles. Therefore, +the equation for calculating the coefficient ``d_{ii}`` must be adjusted as follows: + +```math +d_{ii} = -\Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij} - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. +``` + +The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: + +```math +\begin{align*} +p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} + - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj}p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right. \\ +& \quad - \left. \sum_b m_b \sum_j d_{ij} p_j^l \nabla W_{ib} \right). +\end{align*} +``` + +### Pressure Zeroing +If pressure zeroing is used instead, the pressure value of a boundary particle ``p_b`` +​is assumed to be zero. Consequently, boundary particles do not contribute to the pressure +forces acting on fluid particles. +In this case, the computation of the coefficient ``d_{ii}`` remains unchanged and is given by: + +```math +d_{ii} = -\Delta t^2 \sum_j \frac{m_j}{\rho_i^2} \nabla W_{ij}. +``` + +The equation for the relaxed Jacobi iteration remains the same as in the pressure mirroring +approach. However, the contribution from boundary particles vanishes due to their zero +pressure: + +```math +\begin{align*} +p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} + - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj}p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right. \\ +& \quad - \left. \sum_b m_b \sum_j d_{ij} p_j^l \nabla W_{ib} \right). +\end{align*} +``` \ No newline at end of file diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index a25084eca5..9394d9a33f 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -133,6 +133,7 @@ stepsize_callback = StepsizeCallback(cfl=0.9) callbacks = CallbackSet(info_callback, saving_callback, stepsize_callback, extra_callback, density_reinit_cb, saving_paper) -sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +time_integration_scheme = CarpenterKennedy2N54(williamson_condition=false) +sol = solve(ode, time_integration_scheme, dt=1.0, # This is overwritten by the stepsize callback save_everystep=false, callback=callbacks); diff --git a/examples/fluid/dam_break_2d_iisph.jl b/examples/fluid/dam_break_2d_iisph.jl new file mode 100644 index 0000000000..7a55618249 --- /dev/null +++ b/examples/fluid/dam_break_2d_iisph.jl @@ -0,0 +1,48 @@ +# 2D dam break simulation using implicit incompressible SPH (IISPH) +using TrixiParticles + +fluid_particle_spacing = 0.015 + +# Load setup from dam break example +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + fluid_particle_spacing=fluid_particle_spacing, + sol=nothing, ode=nothing) + +# IISPH doesn't require a large compact support like WCSPH and performs worse with a typical +# smoothing length used for WCSPH. +smoothing_length = 1.0 * fluid_particle_spacing +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +# This kernel slightly overestimates the density, so we reduce the mass slightly +# to obtain a density slightly below the reference density. +# Otherwise, the fluid will jump slightly at the beginning of the simulation. +tank.fluid.mass .*= 0.995 + +# Calculate kinematic viscosity for the viscosity model. +# Only ViscosityAdami and ViscosityMorris can be used for IISPH simulations since they don't +# require a speed of sound. +nu = 0.02 * smoothing_length * sound_speed / 8 +viscosity = ViscosityAdami(; nu) + +# Use IISPH as fluid system +time_step = 1e-3 +fluid_system = ImplicitIncompressibleSPHSystem(tank.fluid, smoothing_kernel, + smoothing_length, fluid_density, + viscosity=ViscosityAdami(nu=nu), + acceleration=(0.0, -gravity), + min_iterations=2, + max_iterations=30, + time_step=time_step) + +# Run the dam break simulation with these changes +trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), + viscosity_fluid=ViscosityAdami(nu=nu), + smoothing_kernel=smoothing_kernel, + smoothing_length=smoothing_length, + fluid_system=fluid_system, + boundary_density_calculator=PressureZeroing(), + tspan=tspan, + state_equation=nothing, + callbacks=CallbackSet(info_callback, saving_callback), + time_integration_scheme=SymplecticEuler(), dt=time_step) diff --git a/ext/TrixiParticlesOrdinaryDiffEqExt.jl b/ext/TrixiParticlesOrdinaryDiffEqExt.jl index 7178978156..fa9eb5fe56 100644 --- a/ext/TrixiParticlesOrdinaryDiffEqExt.jl +++ b/ext/TrixiParticlesOrdinaryDiffEqExt.jl @@ -140,7 +140,8 @@ end @threaded semi for particle in each_integrated_particle(system) for i in 1:ndims(system) du_system[i, - particle] = duprev_system[i, particle] + dt * kdu_system[i, particle] + particle] = duprev_system[i, particle] + + dt * kdu_system[i, particle] end end end diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index f25a3f0019..3517478c74 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -65,7 +65,8 @@ include("visualization/recipes_plots.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem, - WallBoundarySystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySystem + WallBoundarySystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySystem, + ImplicitIncompressibleSPHSystem export BoundaryZone, InFlow, OutFlow, BidirectionalFlow export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback, PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateReachedCallback, diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 933b034fc0..8f0ec8d21e 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -1024,6 +1024,15 @@ function check_configuration(system::TotalLagrangianSPHSystem, systems, nhs) end end +function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, nhs) + foreach_system(systems) do neighbor + if neighbor isa WeaklyCompressibleSPHSystem + throw(ArgumentError("`ImplicitIncompressibleSPHSystem` cannot be used together with + `WeaklyCompressibleSPHSystem`")) + end + end +end + function check_configuration(system::OpenBoundarySystem, systems, neighborhood_search::PointNeighbors.AbstractNeighborhoodSearch) (; boundary_model, boundary_zones) = system diff --git a/src/io/io.jl b/src/io/io.jl index ed245a4c7e..97dd2243ea 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -98,6 +98,18 @@ function add_system_data!(system_data, system::AbstractFluidSystem) end end +function add_system_data!(system_data, system::ImplicitIncompressibleSPHSystem) + system_data["system_type"] = type2string(system) + system_data["particle_spacing"] = particle_spacing(system, 1) + system_data["density_calculator"] = "SummationDensity" + system_data["smoothing_kernel"] = type2string(system.smoothing_kernel) + system_data["smoothing_length"] = system.cache.smoothing_length + system_data["acceleration"] = system.acceleration + system_data["pressure_acceleration_formulation"] = nameof(system.pressure_acceleration_formulation) + add_system_data!(system_data, shifting_technique(system)) + add_system_data!(system_data, system.viscosity) +end + function add_system_data!(system_data, system::TotalLagrangianSPHSystem) system_data["system_type"] = type2string(system) system_data["particle_spacing"] = particle_spacing(system, 1) diff --git a/src/schemes/fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl b/src/schemes/fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl new file mode 100644 index 0000000000..c67930ab8f --- /dev/null +++ b/src/schemes/fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl @@ -0,0 +1 @@ +include("system.jl") diff --git a/src/schemes/fluid/implicit_incompressible_sph/rhs.jl b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl new file mode 100644 index 0000000000..ca7c4b2c35 --- /dev/null +++ b/src/schemes/fluid/implicit_incompressible_sph/rhs.jl @@ -0,0 +1,60 @@ +# Computes the forces that particles in `particle_system` experience from particles +# in `neighbor_system` and updates `dv` accordingly. +# It takes into account pressure forces, viscosity, and for `ContinuityDensity` updates the density +# using the continuity equation. +function interact!(dv, v_particle_system, u_particle_system, + v_neighbor_system, u_neighbor_system, + particle_system::ImplicitIncompressibleSPHSystem, + neighbor_system, semi) + sound_speed = system_sound_speed(particle_system) #TODO + system_coords = current_coordinates(u_particle_system, particle_system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_neighbor(particle_system, neighbor_system, + system_coords, neighbor_system_coords, semi; + points=each_integrated_particle(particle_system)) do particle, + neighbor, + pos_diff, + distance + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds current_density(v_particle_system, particle_system, particle) + rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor) + + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + # The following call is equivalent to + # `p_a = particle_pressure(v_particle_system, particle_system, particle)` + # `p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor)` + # Only when the neighbor system is a `WallBoundarySystem` or a `TotalLagrangianSPHSystem` + # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is + # the pressure of the fluid particle. + p_a, + p_b = @inbounds particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) + + dv_pressure = pressure_acceleration(particle_system, neighbor_system, + particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, grad_kernel, nothing) + + # Propagate `@inbounds` to the viscosity function, which accesses particle data + dv_viscosity_ = @inbounds dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) + + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + end + end + return dv +end diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl new file mode 100644 index 0000000000..0e15e408f5 --- /dev/null +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -0,0 +1,565 @@ +""" + ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, smoothing_length, + reference_density; + viscosity=nothing, + acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), + omega=0.5, max_error=0.1, min_iterations=2, + max_iterations=20, time_step) + +System for particles of a fluid. +The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system +for the pressure so that density remains within a specified tolerance of the rest value. +See [Implicit Incompressible SPH](@ref iisph) for more details on the method. +!!! note "Time Integration" + This system only supports time integration with `SymplecticEuler()`. No other schemes are currently supported. + +# Arguments +- `initial_condition`: [`InitialCondition`](@ref) representing the system's particles. +- `smoothing_kernel`: Smoothing kernel to be used for this system. + See [Smoothing Kernels](@ref smoothing_kernel). +- `smoothing_length`: Smoothing length to be used for this system. + See [Smoothing Kernels](@ref smoothing_kernel). +- `reference_density`: Reference density used for the fluid particles + +# Keyword Arguments +- `viscosity`: Currently, only [`ViscosityMorris`](@ref) + and [`ViscosityAdami`](@ref) are supported. +- `acceleration`: Acceleration vector for the system. (default: zero vector) +- `omega = 0.5`: Relaxation parameter for the relaxed jacobi scheme +- `max_error = 0.1`: Maximum error (in %) for the termination condition in the relaxed jacobi scheme +- `min_iterations = 2`: Minimum number of iterations in the relaxed jacobi scheme, independent from the termination condition +- `max_iterations = 20`: Maximum number of iterations in the relaxed jacobi scheme, independent from the termination condition +- `time_step`: Time step size used for the simulation +""" +struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, + IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS} + initial_condition :: IC + mass :: ARRAY1D # Array{ELTYPE, 1} + pressure :: ARRAY1D + smoothing_kernel :: K + smoothing_length :: ELTYPE + reference_density :: ELTYPE + acceleration :: SVector{NDIMS, ELTYPE} + viscosity :: V + pressure_acceleration_formulation :: PF + surface_normal_method :: Nothing # TODO + surface_tension :: Nothing # TODO + particle_refinement :: Nothing # TODO + density :: ARRAY1D + predicted_density :: ARRAY1D + advection_velocity :: ARRAY2D # Array{ELTYPE, 2} + d_ii :: ARRAY2D # Eq. 9 + a_ii :: ARRAY1D # Diagonal elements of the implicit pressure equation (Eq. 6) + sum_d_ij_pj :: ARRAY2D # \sum_j d_{ij} p_j (Eq. 10) + sum_term :: ARRAY1D # Sum term of Eq. 13 + density_error :: ARRAY1D # Temporary storage for parallel reduction + omega :: ELTYPE # Relaxed Jacobi parameter + max_error :: ELTYPE # maximal error of the average density deviation (in %) + min_iterations :: Int # minimum number of iterations in the pressure solver + max_iterations :: Int # maximum number of iterations in the pressure solver + time_step :: ELTYPE + artificial_sound_speed :: ELTYPE # TODO + cache :: C +end + +# The default constructor needs to be accessible for Adapt.jl to work with this struct. +# See the comments in general/gpu.jl for more details. +function ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, smoothing_length, + reference_density; + viscosity=nothing, + acceleration=ntuple(_ -> 0.0, + ndims(smoothing_kernel)), + omega=0.5, max_error=0.1, min_iterations=2, + max_iterations=20, time_step, + artificial_sound_speed=1000.0) + particle_refinement = nothing # TODO + surface_tension = nothing # TODO + + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) + n_particles = nparticles(initial_condition) + + mass = copy(initial_condition.mass) + pressure = copy(initial_condition.pressure) + + if ndims(smoothing_kernel) != NDIMS + 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")) + end + + if reference_density <= 0 + throw(ArgumentError("`reference_density` must be a positive number")) + end + + if !(0 < max_error <= 100) + throw(ArgumentError("`max_error` is given in percentage, so it must be a number between 0 and 100")) + end + + if min_iterations < 1 + throw(ArgumentError("`min_iterations` must be a positive number")) + end + + if max_iterations < min_iterations + throw(ArgumentError("`min_iterations` must be smaller or equal to `max_iterations`")) + end + + if time_step <= 0 + throw(ArgumentError("`time_step` must be a positive number")) + end + + pressure_acceleration = pressure_acceleration_summation_density + + density = copy(initial_condition.density) + predicted_density = zeros(ELTYPE, n_particles) + a_ii = zeros(ELTYPE, n_particles) + d_ii = zeros(ELTYPE, NDIMS, n_particles) + advection_velocity = zeros(ELTYPE, NDIMS, n_particles) + sum_d_ij_pj = zeros(ELTYPE, NDIMS, n_particles) + sum_term = zeros(ELTYPE, n_particles) + density_error = zeros(ELTYPE, n_particles) + + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)...,) + + return ImplicitIncompressibleSPHSystem(initial_condition, mass, pressure, + smoothing_kernel, smoothing_length, + reference_density, acceleration_, viscosity, + pressure_acceleration, nothing, surface_tension, + particle_refinement, density, predicted_density, + advection_velocity, d_ii, a_ii, sum_d_ij_pj, + sum_term, density_error, omega, max_error, + min_iterations, max_iterations, time_step, + artificial_sound_speed, cache) +end + +function write_v0!(v0, system::ImplicitIncompressibleSPHSystem) + # This is as fast as a loop with `@inbounds`, but it's GPU-compatible + indices = CartesianIndices(system.initial_condition.velocity) + copyto!(v0, indices, system.initial_condition.velocity, indices) + return v0 +end + +function Base.show(io::IO, system::ImplicitIncompressibleSPHSystem) + @nospecialize system # reduce precompilation time + + print(io, "ImplicitIncompressibleSPHSystem{", ndims(system), "}(") + print(io, system.reference_density) + print(io, ", ", system.smoothing_kernel) + print(io, ", ", system.viscosity) + print(io, ", ", system.acceleration) + print(io, ", ", system.omega) + print(io, ", ", system.max_error) + print(io, ", ", system.min_iterations) + print(io, ", ", system.max_iterations) + print(io, ") with ", nparticles(system), " particles") +end + +function Base.show(io::IO, ::MIME"text/plain", system::ImplicitIncompressibleSPHSystem) + @nospecialize system # reduce precompilation time + + if get(io, :compact, false) + show(io, system) + else + summary_header(io, "ImplicitIncompressibleSPHSystem{$(ndims(system))}") + summary_line(io, "#particles", nparticles(system)) + summary_line(io, "reference density", system.reference_density) + summary_line(io, "density calculator", "SummationDensity") + summary_line(io, "smoothing kernel", system.smoothing_kernel |> typeof |> nameof) + summary_line(io, "viscosity", system.viscosity) + summary_line(io, "acceleration", system.acceleration) + summary_line(io, "omega", system.omega) + summary_line(io, "max_error", system.max_error) + summary_line(io, "min_iterations", system.min_iterations) + summary_line(io, "max_iterations", system.max_iterations) + summary_footer(io) + end +end + +@inline source_terms(system::ImplicitIncompressibleSPHSystem) = nothing + +@inline function Base.eltype(::ImplicitIncompressibleSPHSystem{<:Any, ELTYPE}) where {ELTYPE} + return ELTYPE +end + +@inline function surface_tension_model(system::ImplicitIncompressibleSPHSystem) + return nothing +end + +@propagate_inbounds function current_pressure(v, system::ImplicitIncompressibleSPHSystem) + return system.pressure +end + +@propagate_inbounds function current_density(v, system::ImplicitIncompressibleSPHSystem) + return system.density +end + +# TODO: What do we do with the sound speed? This is needed for the viscosity. +@inline system_sound_speed(system::ImplicitIncompressibleSPHSystem) = system.artificial_sound_speed + +# Calculates the pressure values by solving a linear system with a relaxed Jacobi scheme +function update_quantities!(system::ImplicitIncompressibleSPHSystem, v, u, + v_ode, u_ode, semi, t) + @trixi_timeit timer() "predict advection" predict_advection(system, v, u, v_ode, u_ode, + semi, t) + + @trixi_timeit timer() "pressure solver" pressure_solve(system, v, u, v_ode, u_ode, semi, + t) + + return system +end + +function predict_advection(system, v, u, v_ode, u_ode, semi, t) + (; density, predicted_density, advection_velocity, pressure, time_step) = system + d_ii_array = system.d_ii + + sound_speed = system_sound_speed(system) # TODO + + # Compute density by kernel summation + summation_density!(system, semi, u, u_ode, density) + + # Initialize arrays + v_particle_system = wrap_v(v_ode, system, semi) + predicted_density .= density + set_zero!(d_ii_array) + + @threaded semi for particle in each_integrated_particle(system) + # Initialize the advection velocity with the current velocity plus the system acceleration + v_particle = current_velocity(v_particle_system, system, particle) + for i in 1:ndims(system) + advection_velocity[i, + particle] = v_particle[i] + + time_step * system.acceleration[i] + end + end + + # Compute predicted velocity + foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, semi; + points=each_integrated_particle(system)) do particle, + neighbor, + pos_diff, + distance + m_a = @inbounds hydrodynamic_mass(system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) + + rho_a = @inbounds current_density(v_particle_system, system, particle) + rho_b = @inbounds current_density(v_neighbor_system, neighbor_system, neighbor) + + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + dv_viscosity_ = @inbounds dv_viscosity(system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) + # Add all other non-pressure forces + for i in 1:ndims(system) + @inbounds advection_velocity[i, particle] += time_step * dv_viscosity_[i] + end + # Calculate d_ii with eq. 9 in Ihmsen et al. (2013) + for i in 1:ndims(system) + d_ii_array[i, + particle] += calculate_d_ii(neighbor_system, m_b, rho_a, + grad_kernel[i], time_step) + end + end + end + + # Set initial pressure (p_0) to a half of the current pressure value + @threaded semi for particle in each_integrated_particle(system) + pressure[particle] = pressure[particle] / 2 + end + + # Calculation the diagonal elements (a_ii-values) + calculate_diagonal_elements!(system, v, u, v_ode, u_ode, semi) + + # Calculate the predicted density (with the continuity equation and predicted velocities) + foreach_system(semi) do neighbor_system + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, + neighbor_system_coords, semi, + points=each_integrated_particle(system)) do particle, + neighbor, + pos_diff, + distance + # Calculate the predicted velocity differences + advection_velocity_diff = predicted_velocity(system, particle) - + predicted_velocity(neighbor_system, neighbor) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + # Compute \rho_adv in eq. 4 in Ihmsen et al. (2013) + predicted_density[particle] += time_step * m_b * + dot(advection_velocity_diff, grad_kernel) + end + end +end + +function calculate_diagonal_elements!(system, v, u, v_ode, u_ode, semi) + (; a_ii, time_step) = system + + set_zero!(a_ii) + + foreach_system(semi) do neighbor_system + calculate_diagonal_elements!(a_ii, system, neighbor_system, v, u, v_ode, u_ode, + semi, time_step) + end +end + +# Calculation of the contribution of the fluid particles to the diagonal elements (a_ii-values) +# according to eq. 12 in Ihmsen et al. (2013). +function calculate_diagonal_elements!(a_ii, system, neighbor_system, v, u, v_ode, u_ode, + semi, time_step) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, semi; + points=each_integrated_particle(system)) do particle, neighbor, + pos_diff, distance + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + # Compute d_ji according to eq. 9 in Ihmsen et al. (2013). + # Note that we compute d_ji and not d_ij. We can use the antisymmetry + # of the kernel gradient and just flip the sign of W_ij to obtain W_ji. + d_ji_ = -time_step^2 * hydrodynamic_mass(system, particle) / + system.density[particle]^2 * (-grad_kernel) + + d_ii_ = d_ii(system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + # According to eq. 12 in Ihmsen et al. (2013) + a_ii[particle] += m_b * dot((d_ii_ - d_ji_), grad_kernel) + end +end + +# Calculation of the contribution of the Abstractboundary particles the diagonal elements (a_ii-values) +# according to Ihmsen et al. (2013) +function calculate_diagonal_elements!(a_ii, system, neighbor_system::AbstractBoundarySystem, + v, u, + v_ode, u_ode, semi, time_step) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, + system_coords, neighbor_system_coords, semi; + points=each_integrated_particle(system)) do particle, neighbor, + pos_diff, distance + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + + d_ii_ = d_ii(system, particle) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + + # Contribution to the diagonal elements without d_ji value (see eq. 16 in Ihmsen et al. (2013)) + a_ii[particle] += m_b * dot(d_ii_, grad_kernel) + end +end + +# Calculate pressure values with iterative pressure solver (relaxed Jacobi scheme) +function pressure_solve(system, v, u, v_ode, u_ode, semi, t) + (; reference_density, max_error, min_iterations, max_iterations, time_step) = system + + l = 1 + terminate = false + # Convert relative error in percent to absolute error + eta = max_error * reference_density / 100 + while (!terminate) + @trixi_timeit timer() "pressure solver iteration" begin + avg_density_error = pressure_solve_iteration(system, u, + u_ode, semi, time_step) + # Update termination condition + terminate = (avg_density_error <= eta && l >= min_iterations) || + l >= max_iterations + l += 1 + end + end +end + +function pressure_solve_iteration(system, u, u_ode, semi, time_step) + (; reference_density, sum_d_ij_pj, sum_term, pressure, predicted_density, a_ii, + omega, density_error) = system + + set_zero!(sum_d_ij_pj) + + system_coords = current_coordinates(u, system) + + foreach_point_neighbor(system, system, system_coords, system_coords, semi; + points=each_integrated_particle(system)) do particle, neighbor, + pos_diff, distance + # Calculate the sum d_ij * p_j over all neighbors j for each particle i + # (Ihmsen et al. 2013, eq. 13) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + p_b = pressure[neighbor] + d_ab = calculate_d_ij(system, neighbor, grad_kernel, time_step) + sum_dij_pj_ = d_ab * p_b + + for i in 1:ndims(system) + sum_d_ij_pj[i, particle] += sum_dij_pj_[i] + end + end + + # Calculate the large sum in eq. 13 of Ihmsen et al. (2013) for each particle (as `sum_term`) + set_zero!(sum_term) + + foreach_system(semi) do neighbor_system + calculate_sum_term!(sum_term, system, neighbor_system, + pressure, u, u_ode, semi, time_step) + end + + # Update the pressure values + avg_density_error = zero(eltype(system)) + @threaded semi for particle in eachparticle(system) + # Removing instabilities by avoiding to divide by very low values of `a_ii`. + # This is not mentioned in the paper but done in SPlisHSPlasH as well. + if abs(a_ii[particle]) > 1.0e-9 + pressure[particle] = max((1-omega) * pressure[particle] + + omega / a_ii[particle] * + (reference_density - predicted_density[particle] - + sum_term[particle]), 0) + else + pressure[particle] = zero(pressure[particle]) + end + # Calculate the average density error for the termination condition + if (pressure[particle] != 0.0) + new_density = a_ii[particle]*pressure[particle] + sum_term[particle] - + (reference_density - predicted_density[particle]) + + reference_density + density_error[particle] = (new_density - reference_density) + end + end + avg_density_error = sum(density_error) / nparticles(system) + + return avg_density_error +end + +# Function barrier for type stability +function calculate_sum_term!(sum_term, system, neighbor_system, + pressure, u, u_ode, semi, time_step) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + system_coords = current_coordinates(u, system) + neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) + + foreach_point_neighbor(system, neighbor_system, system_coords, + neighbor_system_coords, semi; + points=each_integrated_particle(system)) do particle, neighbor, + pos_diff, distance + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + sum_term[particle] += calculate_sum_term(system, neighbor_system, particle, + neighbor, pressure, grad_kernel, time_step) + end + + return sum_term +end + +@propagate_inbounds function predicted_velocity(system::ImplicitIncompressibleSPHSystem, + particle) + return extract_svector(system.advection_velocity, system, particle) +end + +@propagate_inbounds function predicted_velocity(system::AbstractBoundarySystem, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +@propagate_inbounds function d_ii(system::ImplicitIncompressibleSPHSystem, particle) + return extract_svector(system.d_ii, system, particle) +end + +@propagate_inbounds function d_ii(system::AbstractBoundarySystem, particle) + return zero(SVector{ndims(system), eltype(system)}) +end + +@propagate_inbounds function sum_dij_pj(system::ImplicitIncompressibleSPHSystem, particle) + return extract_svector(system.sum_d_ij_pj, system, particle) +end + +# Calculates a summand for the calculation of the d_ii values +function calculate_d_ii(system::ImplicitIncompressibleSPHSystem, m_b, rho_a, grad_kernel, + time_step) + return -time_step^2 * m_b / rho_a^2 * grad_kernel +end + +# Calculates a summand for the calculation of the d_ii values +function calculate_d_ii(system::AbstractBoundarySystem, m_b, rho_a, grad_kernel, time_step) + return calculate_d_ii(system, system.boundary_model, m_b, rho_a, + grad_kernel, time_step) +end + +# Calculates a summand for the calculation of the d_ii values +function calculate_d_ii(system, boundary_model::BoundaryModelDummyParticles, m_b, rho_a, + grad_kernel, time_step) + return calculate_d_ii(system, boundary_model, boundary_model.density_calculator, m_b, + rho_a, grad_kernel, time_step) +end + +# Calculates a summand for the calculation of the d_ii values (pressure mirroring) +function calculate_d_ii(system, boundary_model, density_calculator::PressureMirroring, m_b, + rho_a, grad_kernel, time_step) + # The linear system to solve originates from the pressure acceleration: + # ∑ m_j (p_i/ρ_i + p_b/ρ_b) ∇W_ij. + # With pressure mirroring, this becomes + # ∑ m_j (p_i/ρ_i + p_i/ρ_i) ∇W_ij, + # which simplifies to + # ∑ m_j (2 * p_i / ρ_i) ∇W_ij. + # Therefore, the diagonal element in the system now appears with a factor 2, + # whereas the entry `ij` disappears from the system. + return -time_step^2 * 2 * m_b / rho_a^2 * grad_kernel +end + +# Calculates a summand for the calculation of the d_ii values (pressure zeroing) +function calculate_d_ii(system, boundary_model, density_calculator::PressureZeroing, m_b, + rho_a, grad_kernel, time_step) + return -time_step^2 * m_b / rho_a^2 * grad_kernel +end + +# Calculates the d_ij value for a particle `i` and its neighbor `j` from eq. 9 in Ihmsen et al. (2013). +# Note that `i` is only implicitly included in the kernel gradient. +function calculate_d_ij(system, particle_j, grad_kernel, time_step) + # (Δt)^2 * m_j / ρ_j ^2 * ∇W_ij + return -time_step^2 * hydrodynamic_mass(system, particle_j) / + system.density[particle_j]^2 * grad_kernel +end + +# Calculate the large sum in eq. 13 of Ihmsen et al. (2013) for each particle (as `sum_term`) +function calculate_sum_term(system, neighbor_system::ImplicitIncompressibleSPHSystem, + particle, neighbor, pressure, grad_kernel, time_step) + m_j = hydrodynamic_mass(neighbor_system, neighbor) + sum_dik_pk = sum_dij_pj(system, particle) + d_jj = d_ii(neighbor_system, neighbor) + p_i = pressure[particle] + p_j = pressure[neighbor] + sum_djk_pk = sum_dij_pj(neighbor_system, neighbor) + d_ji = calculate_d_ij(system, particle, -grad_kernel, time_step) + + # Equation 13 of Ihmsen et al. (2013): + # m_j * (\sum_k d_ik * p_k - d_jj * p_j - \sum_{k != i} d_jk * p_k) * ∇W_ij + return m_j * dot(sum_dik_pk - d_jj * p_j - (sum_djk_pk - d_ji * p_i), grad_kernel) +end + +function calculate_sum_term(system, neighbor_system::AbstractBoundarySystem, particle, + neighbor, + pressure, grad_kernel, time_step) + sum_dik_pk = sum_dij_pj(system, particle) + m_j = hydrodynamic_mass(neighbor_system, neighbor) + + # Equation 16 of Ihmsen et al. (2013): + # m_j * sum_k d_ik * p_k * ∇W_ij + return m_j * dot(sum_dik_pk, grad_kernel) +end diff --git a/src/schemes/schemes.jl b/src/schemes/schemes.jl index c181d0ed24..026f41acbd 100644 --- a/src/schemes/schemes.jl +++ b/src/schemes/schemes.jl @@ -7,10 +7,13 @@ include("structure/discrete_element_method/discrete_element_method.jl") # Monaghan-Kajtar repulsive boundary particles require the `WallBoundarySystem` # and the `TotalLagrangianSPHSystem`. include("boundary/wall_boundary/monaghan_kajtar.jl") +# Implicit incompressible SPH requires the `BoundarySPHSystem` +include("fluid/implicit_incompressible_sph/implicit_incompressible_sph.jl") # Include rhs for all schemes include("fluid/weakly_compressible_sph/rhs.jl") include("fluid/entropically_damped_sph/rhs.jl") +include("fluid/implicit_incompressible_sph/rhs.jl") include("boundary/wall_boundary/rhs.jl") include("structure/total_lagrangian_sph/rhs.jl") include("structure/discrete_element_method/rhs.jl") diff --git a/test/examples/examples_fluid.jl b/test/examples/examples_fluid.jl index ddb9f64208..2c6d379069 100644 --- a/test/examples/examples_fluid.jl +++ b/test/examples/examples_fluid.jl @@ -204,6 +204,31 @@ end end + @trixi_testset "fluid/dam_break_2d_iisph.jl" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d_iisph.jl"), + tspan=(0.0, 0.1)) [ + 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_2d_iisph.jl with PressureMirroring" begin + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d_iisph.jl"), + tspan=(0.0, 0.1), + boundary_density_calculator=PressureMirroring()) [ + 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_2d_gpu.jl" begin @trixi_test_nowarn trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/schemes/fluid/rhs.jl b/test/schemes/fluid/rhs.jl index 291bc6fee4..d22b811479 100644 --- a/test/schemes/fluid/rhs.jl +++ b/test/schemes/fluid/rhs.jl @@ -54,7 +54,8 @@ (p_a, p_b) in pressures, grad_kernel in grad_kernels @testset verbose=true "$system_name" for system_name in [ "WCSPH", - "EDAC" + "EDAC", + "IISPH" ] if system_name == "WCSPH" system = WeaklyCompressibleSPHSystem(fluid, @@ -69,6 +70,12 @@ smoothing_length, 0.0; density_calculator, pressure_acceleration) + elseif system_name == "IISPH" + system = ImplicitIncompressibleSPHSystem(fluid, + smoothing_kernel, + smoothing_length, + 1000.0, + time_step=0.001) end # Compute accelerations a -> b and b -> a @@ -132,21 +139,35 @@ pressure_acceleration=nothing, density_calculator=density_calculator, smoothing_length, 0.0) + + system_iisph = ImplicitIncompressibleSPHSystem(fluid, smoothing_kernel, + smoothing_length, 1000.0, + time_step=0.001) + n_particles = TrixiParticles.nparticles(system_edac) # Overwrite `system.pressure` because we skip the update step system_wcsph.pressure .= fluid.pressure - @testset "`$(nameof(typeof(system)))`" for system in (system_wcsph, - system_edac) + system_iisph.pressure .= fluid.pressure + if density_calculator isa TrixiParticles.SummationDensity + systems = (system_wcsph, system_edac, system_iisph) + else + # IISPH is always using `SummationDensity`` + systems = (system_wcsph, system_edac) + end + @testset "`$(nameof(typeof(system)))`" for system in systems u = fluid.coordinates if density_calculator isa SummationDensity # Density is stored in the cache v = fluid.velocity - system.cache.density .= fluid.density - - if system isa EntropicallyDampedSPHSystem - # pressure is integrated + if system isa WeaklyCompressibleSPHSystem + system.cache.density .= fluid.density + elseif system isa EntropicallyDampedSPHSystem + # Pressure is integrated + system.cache.density .= fluid.density v = vcat(fluid.velocity, fluid.pressure') + else + system.density .= fluid.density end else # Density is integrated with `ContinuityDensity` diff --git a/test/systems/iisph_system.jl b/test/systems/iisph_system.jl new file mode 100644 index 0000000000..b8c45f7737 --- /dev/null +++ b/test/systems/iisph_system.jl @@ -0,0 +1,287 @@ +@testset verbose=true "ImplicitIncompressibleSPHSystem" begin + # Use `@trixi_testset` to isolate the mock functions in a separate namespace + @trixi_testset "Constructors" begin + coordinates_ = [ + [1.0 2.0 + 1.0 2.0], + [1.0 2.0 + 1.0 2.0 + 1.0 2.0] + ] + omegas = [ + 0.5, + 0.6 + ] + max_errors = [ + 0.1, + 0.05 + ] + min_iterations_ = [ + 2, + 10 + ] + max_iterations_ = [ + 20, + 30 + ] + time_steps_ = [ + 0.001, + 0.0001 + ] + @testset "$(i+1)D" for i in 1:2 + NDIMS = i + 1 + coordinates = coordinates_[i] + mass = [1.25, 1.5] + density = [990.0, 1000.0] + reference_density = 1000.0 + smoothing_kernel = Val(:smoothing_kernel) + omega = omegas[i] + max_error = max_errors[i] + min_iterations = min_iterations_[i] + max_iterations = max_iterations_[i] + time_step = time_steps_[i] + TrixiParticles.ndims(::Val{:smoothing_kernel}) = i + 1 + smoothing_kernel2 = Val(:smoothing_kernel2) + # The wrong dimension. 2 -> 3, 3 -> 2. + TrixiParticles.ndims(::Val{:smoothing_kernel2}) = i % 2 + 2 + smoothing_length = 0.362 + + initial_condition = InitialCondition(; coordinates, mass, density) + system = ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + omega=omega, + max_error=max_error, + min_iterations=min_iterations, + max_iterations=max_iterations, + time_step=time_step) + + @test system isa ImplicitIncompressibleSPHSystem{NDIMS} + @test system.initial_condition == initial_condition + @test system.mass == mass + @test system.reference_density == reference_density + @test system.smoothing_kernel == smoothing_kernel + @test system.omega == omega + @test TrixiParticles.initial_smoothing_length(system) == smoothing_length + @test system.viscosity === nothing + @test system.acceleration == [0.0 for _ in 1:NDIMS] + @test system.max_error == max_error + @test system.min_iterations == min_iterations + @test system.max_iterations == max_iterations + @test system.time_step == time_step + @test length(system.density) == size(coordinates, 2) + + error_str1 = "`acceleration` must be of length $NDIMS for a $(NDIMS)D problem" + @test_throws ArgumentError(error_str1) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + density, + acceleration=(0.0), + time_step=0.001) + + error_str2 = "smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem" + @test_throws ArgumentError(error_str2) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel2, + smoothing_length, + reference_density, + time_step=0.001) + + error_str3 = "`reference_density` must be a positive number" + @test_throws ArgumentError(error_str3) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + 0.0, + time_step=0.001) + + error_str4 = "`max_error` is given in percentage, so it must be a number between 0 and 100" + @test_throws ArgumentError(error_str4) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + max_error=0.0, + time_step=0.001) + + error_str5 = "`min_iterations` must be a positive number" + @test_throws ArgumentError(error_str5) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + min_iterations=0, + time_step=0.001) + + error_str6 = "`min_iterations` must be smaller or equal to `max_iterations`" + @test_throws ArgumentError(error_str6) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + min_iterations=10, + max_iterations=5, + time_step=0.001) + + error_str7 = "`time_step must be a positive number" + @test_throws ArgumentError(error_str6) ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + min_iterations=10, + max_iterations=5, + time_step=0) + end + end + + # Use `@trixi_testset` to isolate the mock functions in a separate namespace + @trixi_testset "Constructors with Setups" begin + setups = [ + RectangularShape(0.123, (2, 3), (-1.0, 0.1), density=1.0), + RectangularShape(0.123, (2, 3, 2), (-1.0, 0.1, 2.1), density=1.0), + RectangularTank(0.123, (0.369, 0.246), (0.369, 0.369), 1020.0).fluid, + RectangularTank(0.123, (0.369, 0.246, 0.246), (0.369, 0.492, 0.492), + 1020.0).fluid, + SphereShape(0.52, 0.1, (-0.2, 0.123), 1.0) + ] + setup_names = [ + "RectangularShape 2D", + "RectangularShape 3D", + "RectangularTank 2D", + "RectangularTank 3D", + "SphereShape 2D" + ] + NDIMS_ = [2, 3, 2, 3, 2] + reference_densities = [1.0, 1.0, 1020.0, 1020.0, 1.0] + @testset "$(setup_names[i])" for i in eachindex(setups) + setup = setups[i] + NDIMS = NDIMS_[i] + smoothing_kernel = Val(:smoothing_kernel) + reference_density = reference_densities[i] + TrixiParticles.ndims(::Val{:smoothing_kernel}) = NDIMS + smoothing_length = 0.362 + + density_calculator = SummationDensity() + + system = ImplicitIncompressibleSPHSystem(setup, + smoothing_kernel, + smoothing_length, + reference_density, + time_step=0.001) + + @test system isa ImplicitIncompressibleSPHSystem{NDIMS} + @test system.initial_condition == setup + @test system.mass == setup.mass + @test system.smoothing_kernel == smoothing_kernel + @test system.reference_density == reference_density + @test TrixiParticles.initial_smoothing_length(system) == smoothing_length + @test system.viscosity === nothing + @test system.acceleration == [0.0 for _ in 1:NDIMS] + @test length(system.mass) == size(setup.coordinates, 2) + @test length(system.density) == size(setup.coordinates, 2) + end + + # wrong dimension of acceleration + NDIMS_ = [2, 3] + @testset "Wrong acceleration dimension" for i in eachindex(NDIMS_) + setup = setups[i] + NDIMS = NDIMS_[i] + reference_density = 1.0 + smoothing_kernel = Val(:smoothing_kernel) + TrixiParticles.ndims(::Val{:smoothing_kernel}) = NDIMS + smoothing_length = 0.362 + + density_calculator = SummationDensity() + + error_str = "`acceleration` must be of length $NDIMS for a $(NDIMS)D problem" + @test_throws ArgumentError(error_str) ImplicitIncompressibleSPHSystem(setup, + smoothing_kernel, + smoothing_length, + reference_density, + acceleration=(0.0), + time_step=0.001) + end + end + + # Use `@trixi_testset` to isolate the mock functions in a separate namespace + @trixi_testset "show" begin + coordinates = [1.0 2.0 + 1.0 2.0] + mass = [1.25, 1.5] + density = [990.0, 1000.0] + reference_density = 1000.0 + smoothing_kernel = Val(:smoothing_kernel) + TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 + smoothing_length = 0.362 + + initial_condition = InitialCondition(; coordinates, mass, density) + system = ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + time_step=0.001) + + show_compact = "ImplicitIncompressibleSPHSystem{2}(1000.0, Val{:smoothing_kernel}(), nothing, [0.0, 0.0], 0.5, 0.1, 2, 20) with 2 particles" + @test repr(system) == show_compact + show_box = """ + ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ + │ ImplicitIncompressibleSPHSystem{2} │ + │ ══════════════════════════════════ │ + │ #particles: ………………………………………………… 2 │ + │ reference density: ……………………………… 1000.0 │ + │ density calculator: …………………………… SummationDensity │ + │ smoothing kernel: ………………………………… Val │ + │ viscosity: …………………………………………………… nothing │ + │ acceleration: …………………………………………… [0.0, 0.0] │ + │ omega: ……………………………………………………………… 0.5 │ + │ max_error: …………………………………………………… 0.1 │ + │ min_iterations: ……………………………………… 2 │ + │ max_iterations: ……………………………………… 20 │ + └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" + @test repr("text/plain", system) == show_box + end + + @testset verbose=true "write_u0!" begin + coordinates = [1.0 2.0 + 1.0 2.0] + mass = [1.25, 1.5] + density = [990.0, 1000.0] + reference_density = 1000.0 + smoothing_kernel = Val(:smoothing_kernel) + smoothing_length = 0.362 + + initial_condition = InitialCondition(; coordinates, mass, density) + system = ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, reference_density, + time_step=0.001) + + u0 = zeros(TrixiParticles.u_nvariables(system), + TrixiParticles.n_moving_particles(system)) + TrixiParticles.write_u0!(u0, system) + + @test u0 == coordinates + end + + @testset verbose=true "write_v0!" begin + coordinates = [1.0 2.0 + 1.0 2.0] + velocity = 2 * coordinates + mass = [1.25, 1.5] + density = [990.0, 1000.0] + reference_density = 1000.0 + smoothing_kernel = Val(:smoothing_kernel) + smoothing_length = 0.362 + + initial_condition = InitialCondition(; coordinates, velocity, mass, density) + + # SummationDensity (is always in use) + system = ImplicitIncompressibleSPHSystem(initial_condition, + smoothing_kernel, + smoothing_length, + reference_density, + time_step=0.001) + + v0 = zeros(TrixiParticles.v_nvariables(system), + TrixiParticles.n_moving_particles(system)) + TrixiParticles.write_v0!(v0, system) + + @test v0 == velocity + end +end