diff --git a/NEWS.md b/NEWS.md index ad3e44c715..136dd3fd29 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,9 +4,21 @@ 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.5.3 - -### Features +## Version 0.6.0 + +### API Changes + +- Renamed the fluid-system keyword `surface_normal_method` to `surface_method` and added + detection-only surface methods. The old constructor keyword and accessor are deprecated. + +### Features + +- Added an experimental surface-activity-weighted pressure-difference model for corrected + WCSPH and EDAC free-surface simulations. + +## Version 0.5.3 + +### Features - Added normal vectors to `InitialCondition`, with automatic computation for boundary particles generated by `RectangularTank` and `SphereShape` (#1036). diff --git a/docs/src/refs.bib b/docs/src/refs.bib index c479e59210..aac001a9bd 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -727,6 +727,18 @@ @Article{Price2012 publisher = {Elsevier BV}, } +@Article{Quinlan2006, + author = {Quinlan, Nathan J. and Basa, Mihai and Lastiwka, Martin}, + title = {Truncation error in mesh-free particle methods}, + journal = {International Journal for Numerical Methods in Engineering}, + year = {2006}, + volume = {66}, + number = {13}, + pages = {2064--2085}, + doi = {10.1002/nme.1617}, + publisher = {Wiley}, +} + @Article{Ramachandran2019, author = {Ramachandran, Prabhu and Puri, Kunal}, title = {Entropically damped artificial compressibility for SPH}, @@ -764,6 +776,17 @@ @Article{Schoenberg1946 publisher = {American Mathematical Society (AMS)}, } +@Article{Sigalotti2021, + author = {Sigalotti, Leonardo Di G. and Klapp, Jaime and G{\'o}mez Gesteira, Moncho}, + title = {The Mathematics of Smoothed Particle Hydrodynamics ({SPH}) Consistency}, + journal = {Frontiers in Applied Mathematics and Statistics}, + year = {2021}, + volume = {7}, + pages = {797455}, + doi = {10.3389/fams.2021.797455}, + publisher = {Frontiers Media SA}, +} + @Article{Smagorinsky1963, author = {Smagorinsky, Joseph}, title = {General Circulation Experiments with the Primitive Equations. I. The Basic Experiment}, @@ -926,6 +949,18 @@ @Article{Zhang2025 publisher = {AIP Publishing}, } +@Article{Zhu2015, + author = {Zhu, Qirong and Hernquist, Lars and Li, Yuexing}, + title = {Numerical Convergence in Smoothed Particle Hydrodynamics}, + journal = {The Astrophysical Journal}, + year = {2015}, + volume = {800}, + number = {1}, + pages = {6}, + doi = {10.1088/0004-637X/800/1/6}, + publisher = {IOP Publishing}, +} + @article{Zhu2021, author = {Zhu, Yujie and Zhang, Chi and Yu, Yongchuan and Hu, Xiangyu}, title = {A {CAD}-compatible body-fitted particle generator for arbitrarily complex geometry and its application to wave-structure interaction}, diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..116f83d1f2 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -204,6 +204,177 @@ Pages = [joinpath("schemes", "fluid", "viscosity.jl")] ## [Corrections](@id corrections) +### Configuration + +Density and gradient corrections can be configured independently for WCSPH and EDAC: + +```julia +fluid_system = WeaklyCompressibleSPHSystem(initial_condition; + density_calculator=SummationDensity(), + state_equation, smoothing_kernel, + smoothing_length, + density_correction=ShepardKernelCorrection(), + gradient_correction=MixedKernelGradientCorrection()) +``` + +The legacy `correction` keyword remains available for selecting one correction, but it cannot be +combined with the two role-specific keywords. + +| System | Density handling | Supported correction behavior | +|:-------|:-----------------|:------------------------------| +| WCSPH with [`SummationDensity`](@ref) | Density is recomputed algebraically at every RHS evaluation | Shepard can filter density; all gradient corrections are supported | +| WCSPH with [`ContinuityDensity`](@ref) | Density is an evolved ODE variable | Gradient corrections are supported; Shepard is applied only by [`DensityReinitializationCallback`](@ref) | +| EDAC with [`SummationDensity`](@ref) | Density is algebraic and pressure is evolved independently | Same one-pass Shepard limitation as WCSPH; all gradient corrections are supported | +| EDAC with [`ContinuityDensity`](@ref) | Density and pressure are evolved ODE variables | Gradient corrections are supported; continuous Shepard density overwrite is rejected | +| [`ImplicitIncompressibleSPHSystem`](@ref) | Summation density is coupled to the pressure projection | Corrections are not supported | + +IISPH relies on antisymmetric raw kernel gradients throughout its pressure matrix. Supporting an +asymmetric corrected gradient would require rederiving every projection term, so corrections are +intentionally not exposed for IISPH. + +The default pressure acceleration is selected to match the density evolution law. When passing a +formulation explicitly, the caller is responsible for choosing the corresponding pairing: + +| Pressure acceleration | Consistent density calculator | Asymmetric gradient corrections | +|:----------------------|:-------------------|:--------------------------------| +| `pressure_acceleration_summation_density` | [`SummationDensity`](@ref) | Supported | +| `pressure_acceleration_continuity_density` | [`ContinuityDensity`](@ref) | Supported | +| `inter_particle_averaged_pressure` | Either | Supported | +| [`tensile_instability_control`](@ref) | [`ContinuityDensity`](@ref) | Not supported | + +For a correction whose gradients differ at particles ``a`` and ``b``, the conservative extensions +use both ``\widetilde{\nabla}W_{ab}^{(a)}`` and +``\widetilde{\nabla}W_{ba}^{(b)}``. They reduce algebraically to the original formulas when the +gradient is antisymmetric and give equal-and-opposite pair forces for arbitrary corrected +gradients. Tensile instability control has no such extension and is therefore rejected with +`KernelCorrection`, `GradientCorrection`, `BlendedGradientCorrection`, and +`MixedKernelGradientCorrection`. + +### Consistency validation + +The consistency degree of an SPH correction describes which polynomial fields its discrete +operator reproduces exactly. This is different from a convergence order: zeroth-order +consistency means exact constants, while first-order consistency means exact affine fields +([Bonet and Lok (1999)](@cite Bonet1999); [Sigalotti et al. (2021)](@cite Sigalotti2021)). + +The correction operators are validated on regular and perturbed particle patches by comparing +their discrete moments with the analytical identities + +```math +\sum_b V_b \widetilde{\nabla}W_{ab} = \bm{0}, \qquad +\sum_b V_b \widetilde{\nabla}W_{ab}(\bm{x}_b-\bm{x}_a)^T = \bm{I}. +``` + +The local truncation scaling follows by inserting a Taylor expansion into the discrete +interpolation ``I_h f`` and direct gradient ``G_h f``: + +```math +\bm{M}_k = \sum_b V_b(\bm{x}_b-\bm{x}_a)^{\otimes k}W_{ab}, \qquad +\bm{G}_k = \sum_b V_b\widetilde{\nabla}W_{ab} + \otimes(\bm{x}_b-\bm{x}_a)^{\otimes k}. +``` + +```math +I_h f-f_a = (M_0-1)f_a + \bm{M}_1 \cdot \nabla f_a + + \frac{1}{2}\bm{M}_2 : \nabla^2 f_a + \cdots, +``` + +```math +G_h f-\nabla f_a = f_a\bm{G}_0 + (\bm{G}_1-\bm{I})\nabla f_a + + \frac{1}{2}\bm{G}_2 : \nabla^2 f_a + \cdots. +``` + +Here ``M_k=O(h^k)`` and ``G_k=O(h^{k-1})``. Consequently, exact constants give an +``O(h)`` interpolation on a generic one-sided support, while an exact linear gradient gives an +``O(h)`` derivative there. On a symmetric interior support, odd moments cancel and both can +display ``O(h^2)`` local truncation errors. The expected behavior of the implemented operators is: + +| Correction | Enforced discrete moment | Generic/truncated support | Symmetric interior support | +|:-----------|:-------------------------|:--------------------------|:---------------------------| +| [`ShepardKernelCorrection`](@ref) | ``M_0=1`` | ``O(h)`` interpolation | ``O(h^2)`` interpolation | +| [`KernelCorrection`](@ref) | ``\bm{G}_0=\bm{0}`` | Removes the ``O(h^{-1})`` constant leakage, but leaves an ``O(1)`` first-moment error | No guaranteed rate at fixed ``\Delta x/h`` | +| [`GradientCorrection`](@ref) | ``\bm{G}_1=\bm{I}`` for the difference gradient | ``O(h)`` gradient | ``O(h^2)`` gradient | +| [`BlendedGradientCorrection`](@ref) | ``\bm{G}_1`` error scaled by ``1-\lambda`` | Fixed ``\lambda<1`` leaves an ``O(1)`` error | No guaranteed asymptotic rate at fixed ``\Delta x/h`` | +| [`MixedKernelGradientCorrection`](@ref) | ``\bm{G}_0=\bm{0}`` and ``\bm{G}_1=\bm{I}`` | ``O(h)`` gradient | ``O(h^2)`` gradient | + +The Shepard interpolation scalings in this table assume prescribed volumes ``V_b`` that are +consistent with the interpolated field. The current [`SummationDensity`](@ref) implementation +instead forms ``V_b=m_b/\rho_b`` from the uncorrected summation density and performs one +normalization pass. At a truncated free surface with fixed ``\Delta x/h``, this reduces the error +constant but does not remove the ``O(1)`` boundary error. The validation therefore reports the +ideal normalized interpolation and the production summation-density update as separate operators. +The continuity-density reinitialization uses the evolved density as an independent volume source +and therefore recovers the expected Shepard scaling. + +These are local operator scalings on self-similar regular particle patches with +``h\propto\Delta x``; they are not convergence rates of the complete SPH scheme. Classical SPH +also has a particle quadrature error depending on ``\Delta x/h`` and the particle distribution. +Formal convergence without consistency correction generally requires the joint limit +``h\to0``, ``\Delta x/h\to0``, and an increasing neighbor count +([Quinlan et al. (2006)](@cite Quinlan2006); [Zhu et al. (2015)](@cite Zhu2015)). + +The reproducible study reports boundary and interior scalings separately. It prints a Markdown +table and writes its complete data to +`out/correction_convergence.csv`: + +```bash +julia --project=. validation/corrections/convergence.jl +``` + +#### Measured operator scaling + +The following values are from the finest refinement (``N=96`` particles per coordinate +direction) of a cubic manufactured field with a [`WendlandC6Kernel`](@ref), +``h/\Delta x=2``, and prescribed particle volumes. +The boundary sample uses a one-sided kernel support away from the corners; the interior sample +has a complete symmetric support. The measured scaling is calculated between ``N=48`` and +``N=96``. The summation-density rows use the production [`SummationDensity`](@ref) update instead +of prescribed volumes; the reinitialization row uses the evolved continuity density as its volume +source. + +| Method | Operator | Boundary ``L_2`` error | Boundary scaling | Interior ``L_2`` error | Interior scaling | +|:-------|:---------|-----------------------:|-----------------:|-----------------------:|-----------------:| +| Uncorrected | Interpolation | ``2.622e-1`` | ``-0.007`` | ``2.956e-4`` | ``0.543`` | +| [`ShepardKernelCorrection`](@ref) | Normalized interpolation | ``1.567e-3`` | ``1.052`` | ``4.466e-5`` | ``2.014`` | +| Uncorrected | Direct gradient | ``6.339e1`` | ``-1.009`` | ``7.226e-4`` | ``-0.090`` | +| [`KernelCorrection`](@ref) | Direct gradient | ``3.200e-1`` | ``-0.005`` | ``9.735e-4`` | ``-0.068`` | +| [`GradientCorrection`](@ref) | Difference gradient | ``1.123e-2`` | ``1.010`` | ``2.381e-5`` | ``2.016`` | +| [`BlendedGradientCorrection`](@ref), ``\lambda=0.5`` | Difference gradient | ``1.765e-1`` | ``-0.049`` | ``3.539e-4`` | ``-0.170`` | +| [`MixedKernelGradientCorrection`](@ref) | Direct gradient | ``9.173e-3`` | ``1.012`` | ``2.381e-5`` | ``2.016`` | +| Uncorrected | Summation density | ``2.631e-1`` | ``-0.002`` | ``2.537e-4`` | ``0.040`` | +| [`ShepardKernelCorrection`](@ref) | Summation density | ``1.918e-1`` | ``-0.004`` | ``2.557e-4`` | ``0.076`` | +| [`ShepardKernelCorrection`](@ref) | Continuity-density reinitialization | ``3.810e-4`` | ``1.010`` | ``2.401e-6`` | ``2.001`` | + +#### Measured pressure acceleration + +The same study evaluates every supported pressure-acceleration pairing with a manufactured +pressure that vanishes at the left free surface and the exact acceleration +``-\nabla p/\rho``. Since a conservative asymmetric pair uses correction data from both particles, +the interior sample keeps both kernel neighborhoods complete. Each entry below is the interior +``L_2`` error at ``N=96`` followed by the scaling from ``N=48`` to ``N=96``. + +| Correction | Summation-density pressure | Inter-particle, summation density | Continuity-density pressure | Inter-particle, continuity density | +|:-----------|----------------------------:|----------------------------------:|----------------------------:|----------------------------------:| +| None | ``9.557e-4 / -0.162`` | ``9.557e-4 / -0.162`` | ``7.046e-4 / -0.224`` | ``7.046e-4 / -0.224`` | +| [`ShepardKernelCorrection`](@ref) | ``9.557e-4 / -0.162`` | ``9.557e-4 / -0.162`` | Not applicable | Not applicable | +| [`KernelCorrection`](@ref) | ``9.557e-4 / -0.162`` | ``9.557e-4 / -0.162`` | ``9.557e-4 / -0.162`` | ``9.557e-4 / -0.162`` | +| [`GradientCorrection`](@ref) | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | +| [`BlendedGradientCorrection`](@ref), ``\lambda=0.5`` | ``4.613e-4 / -0.355`` | ``4.613e-4 / -0.355`` | ``3.358e-4 / -0.509`` | ``3.358e-4 / -0.509`` | +| [`MixedKernelGradientCorrection`](@ref) | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | ``3.439e-5 / 2.021`` | +| Shepard density + mixed gradient | ``3.465e-5 / 2.004`` | ``3.465e-5 / 2.005`` | Not applicable | Not applicable | + +For positive pressure, [`tensile_instability_control`](@ref) reduces exactly to the uncorrected +continuity-density pressure law and produces the same measured errors. All formulations have a +constant-pressure null response in the complete interior to an absolute acceleration error below +``1e-7``. + +On the truncated free-surface row, none of the conservative pressure operators converges at fixed +``h/\Delta x``; the observed scaling remains approximately zero. The local correction moments +constrain one particle's gradient operator, but do not impose consistency on a conservative pair +assembled from two differently truncated neighborhoods. This limitation is reported rather than +hidden by applying a non-conservative one-sided pressure difference. The complete boundary and +interior data for every variation are written to `out/correction_convergence.csv`. + ```@autodocs Modules = [TrixiParticles] Pages = [joinpath("general", "corrections.jl")] @@ -211,7 +382,7 @@ Pages = [joinpath("general", "corrections.jl")] --- -## [Surface Normals](@id surface_normal) +## [Surface Detection And Normals](@id surface_normal) ### Overview of surface normal calculation in SPH @@ -248,6 +419,44 @@ The calculated normals are normalized to unit vectors: Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. +#### Surface methods and activity + +Fluid systems configure interface geometry with the `surface_method` keyword. Every surface +method computes a smooth `surface_activity`. Methods derived from +`AbstractSurfaceNormalMethod` additionally provide a normal, so normal calculation always +includes detection. + +`ColorfieldSurfaceDetection` computes activity only. `ColorfieldSurfaceNormal` uses the same +colorfield accumulation and additionally filters and stores the gradient as a surface normal. +Different `color_value`s detect interfaces between represented liquids. A constant nonzero +color detects a free surface because its kernel support ends at the unrepresented exterior. +Equal colors do not create an internal interface. + +`surface_activity` is available in particle VTK output and as a custom quantity for +`SolutionSavingCallback` and `PostprocessCallback`. `surf_normal` is written only for +normal-capable methods. + +Point and plane interpolation evaluate the same colorfield gradient. With `cut_off_bnd=true`, +kernel-weighted color contributions also determine whether a point belongs to the reference +phase. This prevents extrapolation through both free surfaces and interfaces with another +liquid while retaining the existing solid-boundary cutoff. + +#### Interface-aware pressure + +`SurfacePressureDifference()` can be passed as `surface_pressure` to WCSPH and EDAC systems. +It blends the configured conservative pressure acceleration with the consistent difference +operator + +```math +-\frac{1}{\rho_a}\sum_b \frac{m_b}{\rho_b} + (p_b - p_a)\nabla_a W_{ab} +``` + +using the target particle's `surface_activity`. Only interactions within the same fluid +system are blended. Boundary and cross-system interactions retain the conservative pressure +formulation. The model is experimental, requires `GradientCorrection` or +`MixedKernelGradientCorrection`, and is not supported by IISPH. + #### Handling noise and errors in normal calculation In regions distant from the interface, the calculated normals may be small or inaccurate due to the diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 11ee987599..54f3bf4342 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -65,12 +65,18 @@ viscosity_fluid = ArtificialViscosityMonaghan(; alpha, beta=0.0) density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) # density_diffusion = DensityDiffusionAntuono(delta=0.1) +correction = nothing +surface_method = nothing +surface_pressure = nothing +reference_particle_spacing = isnothing(surface_method) ? 0 : fluid_particle_spacing + fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothing_length, density_calculator=fluid_density_calculator, state_equation, viscosity=viscosity_fluid, density_diffusion, acceleration=(0.0, -gravity), - correction=nothing, surface_tension=nothing, - reference_particle_spacing=0) + correction, surface_tension=nothing, + surface_method, surface_pressure, + reference_particle_spacing) # ========================================================================================== # ==== Boundary @@ -84,8 +90,8 @@ boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundar boundary_density_calculator, smoothing_kernel, smoothing_length; state_equation, - correction=nothing, - reference_particle_spacing=0, + correction, + reference_particle_spacing, viscosity=viscosity_wall, clip_negative_pressure=true) diff --git a/examples/fluid/sphere_surface_tension_2d.jl b/examples/fluid/sphere_surface_tension_2d.jl index 1abcdcb149..b22e519e1d 100644 --- a/examples/fluid/sphere_surface_tension_2d.jl +++ b/examples/fluid/sphere_surface_tension_2d.jl @@ -53,7 +53,7 @@ fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel=fluid_smoothi density_calculator=ContinuityDensity(), reference_particle_spacing=particle_spacing, acceleration=zeros(length(fluid_size)), - surface_normal_method=ColorfieldSurfaceNormal(), + surface_method=ColorfieldSurfaceNormal(), surface_tension=SurfaceTensionMorris(surface_tension_coefficient=50 * 0.0728)) diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ee3bf532b5..efbca9400c 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -107,7 +107,8 @@ export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, sample_boundary, planar_geometry_to_face export SourceTermDamping export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, - GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection + GradientCorrection, BlendedGradientCorrection, MixedKernelGradientCorrection, + CorrectionConfiguration export nparticles, eachparticle export available_data, kinetic_energy, total_mass, max_pressure, min_pressure, avg_pressure, max_density, min_density, avg_density @@ -115,7 +116,10 @@ export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_p interpolate_plane_2d_vtk export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, SurfaceTensionMomentumMorris -export ColorfieldSurfaceNormal +export AbstractSurfaceMethod, AbstractSurfaceNormalMethod, ColorfieldSurfaceDetection, + ColorfieldSurfaceNormal, surface_method, computes_surface_normal, surface_activity, + surface_normal +export AbstractSurfacePressure, SurfacePressureDifference export SymplecticPositionVerlet export coordinates_eltype diff --git a/src/general/corrections.jl b/src/general/corrections.jl index d97eeaf341..52c32e6f13 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -52,7 +52,8 @@ end ShepardKernelCorrection() Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation -to obtain a 0-th order accurate result, which was first proposed by [Li et al. (1996)](@cite Li1996). +to obtain a zeroth-order consistent result (exact reproduction of constants), which was first +proposed by [Li et al. (1996)](@cite Li1996). The kernel correction coefficient is determined by ```math @@ -61,7 +62,10 @@ c(x) = \sum_{b=1} V_b W_b(x), where ``V_b = m_b / \rho_b`` is the volume of particle ``b``. This correction is applied with [`SummationDensity`](@ref) to correct the density and leads -to an improvement, especially at free surfaces. +to an improvement, especially at free surfaces. With summation density, the current one-pass +implementation uses the provisional density in ``V_b`` and therefore reduces the free-surface +error without guaranteeing convergence. [`DensityReinitializationCallback`](@ref) instead uses +the independently evolved continuity density and realizes the consistent Shepard operator. !!! note - It is also referred to as "0th order correction". @@ -73,7 +77,8 @@ struct ShepardKernelCorrection end KernelCorrection() Kernel correction, as explained by [Bonet (1999)](@cite Bonet1999), uses Shepard interpolation -to obtain a 0-th order accurate result, which was first proposed by Li et al. +to obtain a zeroth-order consistent kernel gradient (an exact zero gradient for constants), +which was first proposed by Li et al. This can be further extended to obtain a kernel corrected gradient as shown by [Basa et al. (2008)](@cite Basa2008). The kernel correction coefficient is determined by @@ -100,7 +105,9 @@ struct KernelCorrection end MixedKernelGradientCorrection() Combines [`GradientCorrection`](@ref) and [`KernelCorrection`](@ref), -which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1999)). +which gives a first-order consistent gradient that exactly differentiates affine fields +(see [Bonet, 1999](@cite Bonet1999)). This consistency degree does not by itself specify the +convergence order of a full SPH discretization. # Notes: - Stability issues, especially when particles separate into small clusters. @@ -108,6 +115,65 @@ which results in a 1st-order-accurate SPH method (see [Bonet, 1999](@cite Bonet1 """ struct MixedKernelGradientCorrection end +@doc raw""" + CorrectionConfiguration(; density=nothing, gradient=nothing) + +Configure density and gradient corrections independently. `density` can be `nothing` or +[`ShepardKernelCorrection`](@ref). `gradient` can be `nothing`, [`KernelCorrection`](@ref), +[`GradientCorrection`](@ref), [`BlendedGradientCorrection`](@ref), or +[`MixedKernelGradientCorrection`](@ref). +""" +struct CorrectionConfiguration{D, G} + density::D + gradient::G + + function CorrectionConfiguration(density::D, gradient::G) where {D, G} + if !(density === nothing || density isa ShepardKernelCorrection) + throw(ArgumentError("`density` must be `nothing` or `ShepardKernelCorrection()`")) + end + if !(gradient === nothing || + gradient isa Union{KernelCorrection, GradientCorrection, + BlendedGradientCorrection, MixedKernelGradientCorrection}) + throw(ArgumentError("unsupported gradient correction `$(typeof(gradient))`")) + end + + return new{D, G}(density, gradient) + end +end + +function CorrectionConfiguration(; density=nothing, gradient=nothing) + return CorrectionConfiguration(density, gradient) +end + +correction_density(::Any) = nothing +correction_density(correction::ShepardKernelCorrection) = correction +correction_density(correction::CorrectionConfiguration) = correction.density + +correction_gradient(::Nothing) = nothing +correction_gradient(::ShepardKernelCorrection) = nothing +correction_gradient(::AkinciFreeSurfaceCorrection) = nothing +correction_gradient(correction) = correction +correction_gradient(correction::CorrectionConfiguration) = correction.gradient + +correction_force(::Any) = nothing +correction_force(correction::AkinciFreeSurfaceCorrection) = correction + +function resolve_correction_configuration(correction, density_correction, + gradient_correction) + if correction !== nothing && + (density_correction !== nothing || gradient_correction !== nothing) + throw(ArgumentError("`correction` cannot be combined with `density_correction` or " * + "`gradient_correction`")) + end + + if density_correction === nothing && gradient_correction === nothing + return correction + end + + return CorrectionConfiguration(; density=density_correction, + gradient=gradient_correction) +end + function kernel_correction_coefficient(system::AbstractFluidSystem, particle) return system.cache.kernel_correction_coefficient[particle] end @@ -166,9 +232,24 @@ function compute_shepard_coeff!(system, system_coords, v_ode, u_ode, semi, end end + sanitize_kernel_correction_coefficient!(kernel_correction_coefficient, system, semi) + return kernel_correction_coefficient end +function sanitize_kernel_correction_coefficient!(coefficient, system, semi) + minimum_coefficient = sqrt(eps(eltype(coefficient))) + + @threaded semi for particle in eachparticle(system) + value = coefficient[particle] + if !isfinite(value) || value <= minimum_coefficient + coefficient[particle] = one(value) + end + end + + return coefficient +end + function dw_gamma(system::AbstractFluidSystem, particle) return extract_svector(system.cache.dw_gamma, system, particle) end @@ -255,9 +336,22 @@ function compute_correction_values!(system, end end - for particle in eachparticle(system), i in axes(dw_gamma, 1) - dw_gamma[i, particle] /= kernel_correction_coefficient[particle] + minimum_coefficient = sqrt(eps(eltype(kernel_correction_coefficient))) + @threaded semi for particle in eachparticle(system) + coefficient = kernel_correction_coefficient[particle] + if !isfinite(coefficient) || coefficient <= minimum_coefficient + kernel_correction_coefficient[particle] = one(coefficient) + for i in axes(dw_gamma, 1) + dw_gamma[i, particle] = zero(eltype(dw_gamma)) + end + else + for i in axes(dw_gamma, 1) + dw_gamma[i, particle] /= coefficient + end + end end + + return kernel_correction_coefficient end @doc raw""" @@ -284,6 +378,9 @@ The gradient correction, as commonly proposed, involves multiplying this gradien The correction matrix $\bm{L}_a$ is computed based on the provided particle configuration, aiming to make the corrected gradient more accurate, especially near domain boundaries. +It gives a first-order consistent gradient by differentiating every affine field exactly. +For smooth fields, the local truncation error is generally ``O(h)`` on asymmetric supports and +``O(h^2)`` on symmetric interior supports. To satisfy ```math @@ -313,6 +410,8 @@ This calculates the following, \tilde\nabla A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i ``` with ``0 \leq \lambda \leq 1`` being the blending factor. +For a fixed ``\lambda < 1``, the uncorrected first-moment error remains and no asymptotic order +improvement is guaranteed. # Arguments - `blending_factor`: Blending factor between corrected and regular SPH gradient. @@ -321,6 +420,10 @@ struct BlendedGradientCorrection{ELTYPE <: Real} blending_factor::ELTYPE function BlendedGradientCorrection(blending_factor) + if !(zero(blending_factor) <= blending_factor <= one(blending_factor)) + throw(ArgumentError("`blending_factor` must be between 0 and 1")) + end + return new{eltype(blending_factor)}(blending_factor) end end @@ -376,8 +479,10 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, semi) do particle, neighbor, pos_diff, distance function kernel_grad_local(correction, smoothing_kernel, pos_diff, distance, smoothing_length_, system, particle) - return smoothing_kernel_grad_unsafe(system, pos_diff, distance, - particle) + # Do not dispatch through `system`: the correction matrix being used + # by that path is the matrix currently being assembled here. + return kernel_grad_unsafe(smoothing_kernel, pos_diff, distance, + smoothing_length_) end # Compute gradient of corrected kernel @@ -426,8 +531,9 @@ function correction_matrix_inversion_step!(corr_matrix, system, semi) @threaded semi for particle in eachparticle(system) L = extract_smatrix(corr_matrix, system, particle) - # The matrix `L` only becomes singular when the particle and all neighbors - # are collinear (in 2D) or lie all in the same plane (in 3D). + # The matrix `L` becomes singular when the particle and all neighbors are collinear + # (in 2D) or lie all in the same plane (in 3D). Nearly singular matrices are also + # rejected below to avoid amplifying particle disorder. # This happens only when two (in 2D) or three (in 3D) particles are isolated, # or in cases where there is only one layer of fluid particles on a wall. # In these edge cases, we just disable the correction and set the corrected @@ -441,7 +547,12 @@ function correction_matrix_inversion_step!(corr_matrix, system, semi) # so `L` is singular if and only if the position vectors X_ab don't span the # full space, i.e., particle a and all neighbors lie on the same line (in 2D) # or plane (in 3D). - if abs(det(L)) < 1.0f-9 + scale = maximum(abs, L) + relative_determinant = abs(det(L)) / scale^ndims(system) + minimum_relative_determinant = sqrt(eps(eltype(L))) + + if !isfinite(relative_determinant) || + relative_determinant < minimum_relative_determinant L_inv = I else L_inv = inv(L) @@ -458,6 +569,14 @@ end create_cache_correction(correction, density, NDIMS, nparticles) = (;) +function create_cache_correction(correction::CorrectionConfiguration, density, NDIMS, + n_particles) + density_cache = create_cache_correction(correction.density, density, NDIMS, n_particles) + gradient_cache = create_cache_correction(correction.gradient, density, NDIMS, + n_particles) + return merge(density_cache, gradient_cache) +end + function create_cache_correction(::ShepardKernelCorrection, density, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density)) end diff --git a/src/general/custom_quantities.jl b/src/general/custom_quantities.jl index 3f1f74c637..0715678c2a 100644 --- a/src/general/custom_quantities.jl +++ b/src/general/custom_quantities.jl @@ -141,3 +141,24 @@ end function avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) return NaN end + +""" + surface_activity + +Return the per-particle smooth surface activity, or `nothing` for systems without a +configured surface method. +""" +function surface_activity(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, + semi, t) + hasproperty(system.cache, :surface_activity) || return nothing + return view(system.cache.surface_activity, each_active_particle(system)) +end + +surface_activity(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing + +function surface_normal(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) + computes_surface_normal(surface_method(system)) || return nothing + return view(system.cache.surface_normal, :, each_active_particle(system)) +end + +surface_normal(system, dv_ode, du_ode, v_ode, u_ode, semi, t) = nothing diff --git a/src/general/interpolation.jl b/src/general/interpolation.jl index 961224fab7..75c5f9e777 100644 --- a/src/general/interpolation.jl +++ b/src/general/interpolation.jl @@ -23,11 +23,11 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` outside the + interpolated fluid domain. Boundaries are detected by comparing their + kernel-weighted mass to that of the fluid. When `ref_system` uses + a colorfield surface method, free and multiphase surfaces are detected + from the same color contributions used for particle surface detection. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -41,7 +41,7 @@ See also: [`interpolate_plane_2d_vtk`](@ref), [`interpolate_plane_3d`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used, which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used, which is not as accurate as a real surface reconstruction. # Examples @@ -115,11 +115,9 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref), - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. - `output_directory="out"`: Directory to save the VTI file. - `filename="plane"`: Name of the VTI file. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -130,7 +128,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_3d`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used, which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used, which is not as accurate as a real surface reconstruction. # Examples @@ -177,6 +175,10 @@ function interpolate_plane_2d_vtk(min_corner, max_corner, resolution, semi, ref_ vtk["density"] = density vtk["velocity"] = velocity vtk["pressure"] = pressure + if hasproperty(results, :surface_activity) + vtk["surface_activity"] = reshape(results.surface_activity, + length(x_range), length(y_range)) + end end end @@ -250,11 +252,9 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -268,7 +268,7 @@ See also: [`interpolate_plane_2d`](@ref), [`interpolate_plane_2d_vtk`](@ref), !!! note - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. # Examples @@ -346,11 +346,9 @@ See also: [`interpolate_points`](@ref), [`interpolate_plane_2d`](@ref), # Keywords - `endpoint=true`: A boolean to include (`true`) or exclude (`false`) the end point in the interpolation. - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -366,7 +364,7 @@ See also: [`interpolate_points`](@ref), [`interpolate_plane_2d`](@ref), - This function is particularly useful for analyzing gradients or creating visualizations along a specified line in the SPH simulation domain. - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. # Examples @@ -431,11 +429,9 @@ See also: [`interpolate_line`](@ref), [`interpolate_plane_2d`](@ref), # Keywords - `smoothing_length=initial_smoothing_length(ref_system)`: The smoothing length used in the interpolation. -- `cut_off_bnd=true`: Boolean to indicate if quantities should be set to `NaN` when the point - is "closer" to the boundary than to the fluid in a kernel-weighted sense. - Or, in more detail, when the boundary has more influence than the fluid - on the density summation in this point, i.e., when the boundary particles - add more kernel-weighted mass than the fluid particles. +- `cut_off_bnd=true`: Set quantities outside the interpolated fluid domain to `NaN`. + Solid boundaries use kernel-weighted mass. A configured colorfield + surface method additionally uses phase-color contributions. - `clip_negative_pressure=false`: One common approach in SPH models is to clip negative pressure values, but this is unphysical. Instead we clip here during interpolation thus only impacting the local interpolated value. @@ -463,7 +459,7 @@ results = interpolate_points(points, semi, ref_system, sol) - This function is particularly useful for analyzing gradients or creating visualizations along a specified line in the SPH simulation domain. - The interpolation accuracy is subject to the density of particles and the chosen smoothing length. - - With `cut_off_bnd`, a density-based estimation of the surface is used which is not as + - With `cut_off_bnd`, a kernel-based estimation of the surface is used which is not as accurate as a real surface reconstruction. """ @inline function interpolate_points(point_coords, semi, ref_system, sol::ODESolution; @@ -555,10 +551,22 @@ end n_points = size(point_coords, 2) ELTYPE = eltype(point_coords) + surface_method_ = surface_method(ref_system) + detect_surface = is_colorfield_surface_method(surface_method_) computed_density = allocate(semi.parallelization_backend, ELTYPE, n_points) other_density = allocate(semi.parallelization_backend, ELTYPE, n_points) shepard_coefficient = allocate(semi.parallelization_backend, ELTYPE, n_points) neighbor_count = allocate(semi.parallelization_backend, Int, n_points) + surface_gradient = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, + (ndims(ref_system), n_points)) : nothing + surface_activity_ = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : nothing + reference_colorfield = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : + nothing + other_colorfield = detect_surface ? + allocate(semi.parallelization_backend, ELTYPE, n_points) : nothing # The wall velocity considers more neighbors, so we need to use # a different Shepard coefficient. shepard_coefficient_wall = allocate(semi.parallelization_backend, ELTYPE, n_points) @@ -568,15 +576,26 @@ end set_zero!(other_density) set_zero!(shepard_coefficient) set_zero!(neighbor_count) + if detect_surface + set_zero!(surface_gradient) + set_zero!(surface_activity_) + set_zero!(reference_colorfield) + set_zero!(other_colorfield) + end cache = create_cache_interpolation(ref_system, n_points, semi) ref_id = system_indices(ref_system, semi) ref_smoothing_kernel = ref_system.smoothing_kernel + interpolation_surface_threshold = detect_surface ? + surface_method_.interpolation_surface_threshold : + zero(ELTYPE) + reference_color = detect_surface ? ref_system.cache.color : 0 # If we neither cut off at the boundary nor include the boundary wall velocity, # we only need to iterate over the reference system. - systems = (cut_off_bnd || include_wall_velocity) ? semi : (ref_system,) + systems = (cut_off_bnd || include_wall_velocity || detect_surface) ? + semi : (ref_system,) foreach_system(systems) do neighbor_system system_id = system_indices(neighbor_system, semi) @@ -586,6 +605,12 @@ end u = wrap_u(u_ode, neighbor_system, semi) neighbor_coords = current_coordinates(u, neighbor_system) + contributes_surface = detect_surface && + has_system_interaction(ref_system, neighbor_system, semi) && + (contributes_to_colorfield(neighbor_system) || + contributes_boundary_colorfield(neighbor_system)) + surface_color = contributes_to_colorfield(neighbor_system) ? + neighbor_system.cache.color : reference_color foreach_point_neighbor(point_coords, neighbor_coords, nhs; parallelization_backend) do point, neighbor, pos_diff, @@ -594,6 +619,22 @@ end volume_b = m_b / current_density(v, neighbor_system, neighbor) W_ab = kernel(ref_smoothing_kernel, distance, smoothing_length) + if contributes_surface + grad_kernel = kernel_grad(ref_smoothing_kernel, pos_diff, distance, + smoothing_length) + for i in 1:ndims(ref_system) + surface_gradient[i, point] += volume_b * surface_color * grad_kernel[i] + end + + if neighbor_system isa AbstractFluidSystem + if surface_color == reference_color + reference_colorfield[point] += volume_b * W_ab + else + other_colorfield[point] += volume_b * W_ab + end + end + end + if include_wall_velocity # The wall velocity considers more neighbors, so we need to use # a different Shepard coefficient. @@ -610,7 +651,10 @@ end interpolate_system!(cache, v, neighbor_system, point, neighbor, volume_b, W_ab, clip_negative_pressure) else - other_density[point] += m_b * W_ab + if cut_off_bnd && + (!detect_surface || !(neighbor_system isa AbstractFluidSystem)) + other_density[point] += m_b * W_ab + end if include_wall_velocity velocity_neighbor_ = current_velocity(v, neighbor_system, neighbor) @@ -622,17 +666,38 @@ end end end - neighbor_count[point] += 1 + if system_id == ref_id || cut_off_bnd || include_wall_velocity + neighbor_count[point] += 1 + end end end @threaded parallelization_backend for point in axes(point_coords, 2) + if detect_surface + normal_norm = zero(ELTYPE) + for i in 1:ndims(ref_system) + normal_norm += surface_gradient[i, point]^2 + end + surface_activity_[point] = gradient_surface_activity(sqrt(normal_norm), + compact_support(ref_smoothing_kernel, + smoothing_length), + surface_method_) + end + + outside_reference_phase = detect_surface && + (reference_colorfield[point] < + interpolation_surface_threshold || + other_colorfield[point] > reference_colorfield[point]) cut_off = computed_density[point] < eps() || - (cut_off_bnd && other_density[point] > computed_density[point]) + (cut_off_bnd && (other_density[point] > computed_density[point] || + outside_reference_phase)) if cut_off # Return NaN values that can be filtered out in ParaView computed_density[point] = NaN neighbor_count[point] = 0 + if detect_surface + surface_activity_[point] = NaN + end # We need to convert the `NamedTuple` to a `Tuple` for GPU compatibility foreach(Tuple(cache)) do field @@ -658,7 +723,10 @@ end end end - return (; computed_density, point_coords, neighbor_count, cache...) + surface_detection_output = detect_surface ? + (; surface_activity=surface_activity_) : (;) + return (; computed_density, point_coords, neighbor_count, surface_detection_output..., + cache...) end @inline function create_cache_interpolation(ref_system::AbstractFluidSystem, n_points, semi) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index a9ed769dc6..4ef14ebf7c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -920,22 +920,15 @@ function check_system_color(systems) system isa AbstractFluidSystem || return false system isa ParticlePackingSystem && return false - return !isnothing(system.surface_tension) || - system.surface_normal_method isa ColorfieldSurfaceNormal + return is_colorfield_surface_method(surface_method(system)) end if requires_color_check + system_ids = findall(system -> system isa AbstractFluidSystem && + !(system isa ParticlePackingSystem), systems) - # Systems that contribute to the colorfield/contact logic. - system_ids = findall(system -> (system isa AbstractFluidSystem && - !(system isa ParticlePackingSystem)) || - system isa WallBoundarySystem || - system isa - RigidBodySystem{<:BoundaryModelDummyParticles}, - systems) - - if length(system_ids) > 1 && sum(i -> systems[i].cache.color, system_ids) == 0 - throw(ArgumentError("If `ColorfieldSurfaceNormal` or a surface tension model is used, at least one participating system must have a color different from 0.")) + if all(i -> iszero(systems[i].cache.color), system_ids) + throw(ArgumentError("If a colorfield surface method or a surface tension model is used, at least one participating system must have a color different from 0.")) end end end diff --git a/src/io/io.jl b/src/io/io.jl index 692dfd8528..b60c2220f9 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -86,7 +86,8 @@ function add_system_data!(system_data, system::AbstractFluidSystem) system_data["pressure_acceleration_formulation"] = nameof(system.pressure_acceleration_formulation) add_system_data!(system_data, shifting_technique(system)) add_system_data!(system_data, system.surface_tension) - add_system_data!(system_data, system.surface_normal_method) + add_system_data!(system_data, system.surface_method) + add_system_data!(system_data, surface_pressure_model(system)) add_system_data!(system_data, system.viscosity) add_system_data!(system_data, system.correction) add_system_data!(system_data, system_state_equation(system)) @@ -107,6 +108,7 @@ function add_system_data!(system_data, system::ImplicitIncompressibleSPHSystem) 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.surface_method) add_system_data!(system_data, system.viscosity) end @@ -306,6 +308,13 @@ function add_system_data!(system_data, system_data["correction_method"]["model"] = type2string(correction) end +function add_system_data!(system_data, correction::CorrectionConfiguration) + system_data["correction_method"] = Dict{String, Any}() + system_data["correction_method"]["model"] = type2string(correction) + system_data["correction_method"]["density"] = type2string(correction.density) + system_data["correction_method"]["gradient"] = type2string(correction.gradient) +end + function add_system_data!(system_data, surface_tension::Union{CohesionForceAkinci, SurfaceTensionAkinci, SurfaceTensionMorris, @@ -315,11 +324,19 @@ function add_system_data!(system_data, system_data["surface_tension"]["surface_tension_coefficient"] = surface_tension.surface_tension_coefficient end -function add_system_data!(system_data, surface_normal_method::ColorfieldSurfaceNormal) - system_data["surface_normal_method"] = Dict{String, Any}() - system_data["surface_normal_method"]["model"] = type2string(surface_normal_method) - system_data["surface_normal_method"]["boundary_contact_threshold"] = surface_normal_method.boundary_contact_threshold - system_data["surface_normal_method"]["ideal_density_threshold"] = surface_normal_method.ideal_density_threshold +function add_system_data!(system_data, surface_method_::ColorfieldSurfaceMethod) + system_data["surface_method"] = Dict{String, Any}() + system_data["surface_method"]["model"] = type2string(surface_method_) + system_data["surface_method"]["computes_surface_normal"] = computes_surface_normal(surface_method_) + system_data["surface_method"]["boundary_contact_threshold"] = surface_method_.boundary_contact_threshold + system_data["surface_method"]["interface_threshold"] = surface_method_.interface_threshold + system_data["surface_method"]["ideal_density_threshold"] = surface_method_.ideal_density_threshold + system_data["surface_method"]["interface_taper_start"] = surface_method_.interface_taper_start + system_data["surface_method"]["interpolation_surface_threshold"] = surface_method_.interpolation_surface_threshold +end + +function add_system_data!(system_data, surface_pressure::SurfacePressureDifference) + system_data["surface_pressure"] = Dict("model" => type2string(surface_pressure)) end function add_system_data!(system_data, boundary_zone::BoundaryZone, indice) diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl index f242fd74b2..43320b2e4e 100644 --- a/src/io/write_vtk.jl +++ b/src/io/write_vtk.jl @@ -329,11 +329,15 @@ function write2vtk!(vtk, v, u, t, system::AbstractFluidSystem) vtk["pressure"] = [current_pressure(v, system, particle) for particle in eachparticle(system)] - if system.surface_normal_method !== nothing - vtk["surf_normal"] = [surface_normal(system, particle) - for particle in eachparticle(system)] + if system.surface_method !== nothing + vtk["surface_activity"] = system.cache.surface_activity vtk["neighbor_count"] = system.cache.neighbor_count vtk["color"] = system.cache.color + + if computes_surface_normal(system.surface_method) + vtk["surf_normal"] = [surface_normal(system, particle) + for particle in eachparticle(system)] + end end if system.surface_tension isa SurfaceTensionMorris || diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index fcde467f2e..0337ada231 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -3,6 +3,7 @@ density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, state_equation=nothing, correction=nothing, + density_correction=nothing, gradient_correction=nothing, clip_negative_pressure=false, reference_particle_spacing=0.0) @@ -20,7 +21,10 @@ Boundary model for [`WallBoundarySystem`](@ref). # Keywords - `state_equation`: This should be the same as for the adjacent fluid system (see e.g. [`StateEquationCole`](@ref)). -- `correction`: Correction method of the adjacent fluid system (see [Corrections](@ref corrections)). +- `correction`: Legacy keyword for one correction method. Cannot be combined with + `density_correction` or `gradient_correction`. +- `density_correction`: Density correction of the adjacent fluid system. +- `gradient_correction`: Gradient correction of the adjacent fluid system. - `viscosity`: Slip (default) or no-slip condition. See description below for further information. - `clip_negative_pressure=false`: Clip negative boundary pressures to avoid sticking @@ -80,6 +84,8 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, state_equation=nothing, correction=nothing, + density_correction=nothing, + gradient_correction=nothing, clip_negative_pressure=false, reference_particle_spacing=0.0) pressure = initial_boundary_pressure(initial_density, density_calculator, @@ -88,6 +94,8 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, ELTYPE = eltype(smoothing_length) @assert length(initial_density) == length(hydrodynamic_mass) n_particles = length(initial_density) + correction = resolve_correction_configuration(correction, density_correction, + gradient_correction) cache = (; create_cache_model(viscosity, n_particles, NDIMS)..., create_cache_model(initial_density, density_calculator, NDIMS)..., @@ -232,24 +240,31 @@ struct PressureBoundaries{ELTYPE} end @inline create_cache_model(correction, density, NDIMS, nparticles) = (;) +function create_cache_model(correction::CorrectionConfiguration, density, NDIMS, + n_particles) + density_cache = create_cache_model(correction.density, density, NDIMS, n_particles) + gradient_cache = create_cache_model(correction.gradient, density, NDIMS, n_particles) + return merge(density_cache, gradient_cache) +end + function create_cache_model(::ShepardKernelCorrection, density, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density)) end function create_cache_model(::KernelCorrection, density, NDIMS, n_particles) - dw_gamma = Array{Float64}(undef, NDIMS, n_particles) + dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density), dw_gamma) end function create_cache_model(::Union{GradientCorrection, BlendedGradientCorrection}, density, NDIMS, n_particles) - correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles) return (; correction_matrix) end function create_cache_model(::MixedKernelGradientCorrection, density, NDIMS, n_particles) - dw_gamma = Array{Float64}(undef, NDIMS, n_particles) - correction_matrix = Array{Float64, 3}(undef, NDIMS, NDIMS, n_particles) + dw_gamma = Array{eltype(density)}(undef, NDIMS, n_particles) + correction_matrix = Array{eltype(density), 3}(undef, NDIMS, NDIMS, n_particles) return (; kernel_correction_coefficient=similar(density), dw_gamma, correction_matrix) end @@ -394,20 +409,43 @@ end @inline function update_pressure!(boundary_model::BoundaryModelDummyParticles, system, v, u, v_ode, u_ode, semi) (; correction, density_calculator) = boundary_model + density_correction = correction_density(correction) + gradient_correction = correction_gradient(correction) compute_pressure!(boundary_model, density_calculator, system, v, u, v_ode, u_ode, semi) - # These are only computed when using corrections - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) - compute_gradient_correction_matrix!(correction, boundary_model, system, u, v_ode, u_ode, - semi) - # `kernel_correct_density!` only performed for `SummationDensity` - kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction, + # Density correction must be applied before assembling gradient corrections so all + # gradient moments use the density that will be used by the RHS. + compute_correction_values!(system, density_correction, u, v_ode, u_ode, semi) + kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, + density_correction, density_calculator) + recompute_pressure_after_density_correction!(boundary_model, density_calculator, + density_correction, system, v, u, v_ode, + u_ode, + semi) + compute_correction_values!(system, gradient_correction, u, v_ode, u_ode, semi) + compute_gradient_correction_matrix!(gradient_correction, boundary_model, system, u, + v_ode, + u_ode, semi) return boundary_model end +function recompute_pressure_after_density_correction!(boundary_model, density_calculator, + correction, system, v, u, v_ode, + u_ode, + semi) + return boundary_model +end + +function recompute_pressure_after_density_correction!(boundary_model, ::SummationDensity, + ::ShepardKernelCorrection, system, v, + u, + v_ode, u_ode, semi) + compute_pressure!(boundary_model, SummationDensity(), system, v, u, v_ode, u_ode, semi) +end + function kernel_correct_density!(boundary_model, v, u, v_ode, u_ode, semi, correction, density_calculator) return boundary_model @@ -428,13 +466,13 @@ function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, MixedKernelGradientCorrection}, boundary_model, system, u, v_ode, u_ode, semi) - (; cache, correction, smoothing_kernel) = boundary_model + (; cache, smoothing_kernel) = boundary_model (; correction_matrix) = cache system_coords = current_coordinates(u, system) compute_gradient_correction_matrix!(correction_matrix, system, system_coords, - v_ode, u_ode, semi, correction, smoothing_kernel) + v_ode, u_ode, semi, corr, smoothing_kernel) end function compute_density!(boundary_model, ::SummationDensity, system, v, u, v_ode, u_ode, diff --git a/src/schemes/boundary/wall_boundary/system.jl b/src/schemes/boundary/wall_boundary/system.jl index f9864ecc05..4072e98abd 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -14,9 +14,9 @@ The interaction between fluid and boundary particles is specified by the boundar - `prescribed_motion`: For moving boundaries, a [`PrescribedMotion`](@ref) can be passed. - `adhesion_coefficient`: Coefficient specifying the adhesion of a fluid to the surface. Note: currently it is assumed that all fluids have the same adhesion coefficient. -- `color_value`: Integer label used for calculation of surface normals. +- `color_value`: Integer label used for colorfield surface calculations. Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation + colorfield surface methods: fluid-boundary surface evaluation reads the resulting boundary colorfield to detect wall contact. """ struct WallBoundarySystem{BM, ELTYPE <: Real, NDIMS, IC, CO, M, IM, @@ -329,7 +329,7 @@ function system_smoothing_kernel(system::WallBoundarySystem{<:BoundaryModelDummy end function system_correction(system::WallBoundarySystem{<:BoundaryModelDummyParticles}) - return system.boundary_model.correction + return correction_gradient(system.boundary_model.correction) end @inline function density_calculator(system::WallBoundarySystem) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 2ea1e4dc4d..25c161b8fa 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -4,6 +4,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system::EntropicallyDampedSPHSystem, neighbor_system, semi) (; sound_speed, density_calculator, correction, nu_edac) = particle_system + gradient_correction = correction_gradient(correction) system_coords = current_coordinates(u_particle_system, particle_system) neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -63,7 +64,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle, neighbor, m_a, m_b, p_a - p_avg, p_b - p_avg, rho_a, rho_b, pos_diff, distance, grad_kernel, - correction) + gradient_correction) dv_particle = Ref(dv_pressure) @inbounds dv_viscosity!(dv_particle, particle_system, neighbor_system, @@ -77,7 +78,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, m_a, m_b, rho_a, rho_b, v_a, v_b, - pos_diff, distance, grad_kernel, correction) + pos_diff, distance, grad_kernel, gradient_correction) @inbounds surface_tension_force!(dv_particle, surface_tension_a, surface_tension_b, diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 6e03cf7837..b2b4069e8a 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -6,9 +6,11 @@ shifting_technique=nothing, alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, + surface_method=nothing, buffer_size=nothing, + surface_pressure=nothing, reference_particle_spacing=0.0, color_value=1, - source_terms=nothing) + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing) System for particles of a fluid. As opposed to the [weakly compressible SPH scheme](@ref wcsph), which uses an equation of state, @@ -38,7 +40,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more from the local pressure (default: `true` when using shifting, `false` otherwise). - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySystem`](@ref). -- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) +- `correction`: Legacy keyword for configuring one correction method. Cannot be + combined with `density_correction` or `gradient_correction`. +- `density_correction`: Density correction method. Currently supports + [`ShepardKernelCorrection`](@ref) with [`SummationDensity`](@ref). +- `gradient_correction`: Gradient correction method. See [Corrections](@ref corrections). - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure, t)` (which are the quantities of a single particle), returning a `Tuple` @@ -50,17 +56,19 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more The keyword argument `acceleration` should be used instead for gravity-like source terms. - `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) -- `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: Integer label used for calculation of surface normals. - Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation - reads the resulting boundary colorfield to detect wall contact. +- `surface_method`: Surface detection or normal method used by this system. + Methods that compute normals always also compute surface + activity. The default is `nothing`, or + `ColorfieldSurfaceNormal()` when required by surface tension. +- `surface_pressure`: Optional interface-aware pressure model. The experimental + [`SurfacePressureDifference`](@ref) blends the conservative + pressure acceleration with a pressure-difference operator. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface detection. Different + values identify represented fluid-fluid interfaces. """ -struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, TV, +struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, PF, SP, TV, AVGP, ST, SRFT, SRFN, B, PR, C} <: AbstractFluidSystem{NDIMS} initial_condition :: IC @@ -73,11 +81,12 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, COR, acceleration :: SVector{NDIMS, ELTYPE} correction :: COR pressure_acceleration_formulation :: PF + surface_pressure :: SP shifting_technique :: TV average_pressure_reduction :: AVGP source_terms :: ST surface_tension :: SRFT - surface_normal_method :: SRFN + surface_method :: SRFN buffer :: B particle_refinement :: PR cache :: C @@ -93,9 +102,12 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth alpha=0.5, viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), - correction=nothing, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing, surface_tension=nothing, - surface_normal_method=nothing, buffer_size=nothing, + surface_method=nothing, surface_normal_method=nothing, + surface_pressure=nothing, + buffer_size=nothing, reference_particle_spacing=0.0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -110,6 +122,11 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth mass = copy(initial_condition.mass) n_particles = length(initial_condition.mass) + correction = resolve_correction_configuration(correction, density_correction, + gradient_correction) + density_correction_ = correction_density(correction) + gradient_correction_ = correction_gradient(correction) + if ndims(smoothing_kernel) != NDIMS throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) end @@ -119,15 +136,16 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end + surface_method = select_surface_method(surface_tension, surface_method, + surface_normal_method) + surface_pressure = validate_surface_pressure(surface_pressure, surface_method, + gradient_correction_) - if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) end - if correction isa ShepardKernelCorrection && + if density_correction_ isa ShepardKernelCorrection && density_calculator isa ContinuityDensity throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end @@ -135,7 +153,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, - correction) + gradient_correction_) avg_pressure_reduction = Val(average_pressure_reduction) @@ -145,8 +163,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth create_cache_shifting(initial_condition, shifting_technique)..., create_cache_avg_pressure_reduction(initial_condition, avg_pressure_reduction)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, @@ -166,17 +183,19 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth EntropicallyDampedSPHSystem{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), typeof(correction), - typeof(pressure_acceleration), typeof(shifting_technique), + typeof(pressure_acceleration), typeof(surface_pressure), + typeof(shifting_technique), typeof(avg_pressure_reduction), typeof(source_terms), - typeof(surface_tension), typeof(surface_normal_method), + typeof(surface_tension), typeof(surface_method), typeof(buffer), Nothing, typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, sound_speed, viscosity, nu_edac, acceleration_, correction, - pressure_acceleration, shifting_technique, + pressure_acceleration, surface_pressure, + shifting_technique, avg_pressure_reduction, source_terms, surface_tension, - surface_normal_method, buffer, + surface_method, buffer, particle_refinement, cache) end @@ -199,7 +218,8 @@ function Base.show(io::IO, system::EntropicallyDampedSPHSystem) print(io, ", ", system.smoothing_kernel) print(io, ", ", system.acceleration) print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) + print(io, ", ", system.surface_method) + print(io, ", ", system.surface_pressure) print(io, ") with ", nparticles(system), " particles") end @@ -228,7 +248,8 @@ function Base.show(io::IO, ::MIME"text/plain", system::EntropicallyDampedSPHSyst typeof(system.average_pressure_reduction).parameters[1] ? "yes" : "no") summary_line(io, "acceleration", system.acceleration) summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) + summary_line(io, "surface method", system.surface_method) + summary_line(io, "surface pressure", system.surface_pressure) summary_footer(io) end end @@ -251,7 +272,9 @@ end @inline buffer(system::EntropicallyDampedSPHSystem) = system.buffer -system_correction(system::EntropicallyDampedSPHSystem) = system.correction +function system_correction(system::EntropicallyDampedSPHSystem) + correction_gradient(system.correction) +end @inline function current_velocity(v, system::EntropicallyDampedSPHSystem) return view(v, 1:ndims(system), :) @@ -299,11 +322,56 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) - compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, - t) + (; correction, density_calculator) = system + density_correction = correction_density(correction) + gradient_correction = correction_gradient(correction) + + # Density correction must be applied before assembling gradient corrections so all + # gradient moments use the density that will be used by the RHS. + compute_correction_values!(system, density_correction, u, v_ode, u_ode, semi) + kernel_correct_density!(system, v, u, v_ode, u_ode, semi, density_correction, + density_calculator) + compute_correction_values!(system, gradient_correction, u, v_ode, u_ode, semi) + compute_gradient_correction_matrix!(gradient_correction, system, u, v_ode, u_ode, semi) + + compute_surface!(system, system.surface_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end +function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, + semi, correction, density_calculator) + return system +end + +function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, + semi, ::ShepardKernelCorrection, ::SummationDensity) + system.cache.density ./= system.cache.kernel_correction_coefficient +end + +function compute_gradient_correction_matrix!(correction, + system::EntropicallyDampedSPHSystem, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + system::EntropicallyDampedSPHSystem, u, + v_ode, u_ode, semi) + (; cache, smoothing_kernel) = system + (; correction_matrix) = cache + + system_coords = current_coordinates(u, system) + + compute_gradient_correction_matrix!(correction_matrix, system, system_coords, + v_ode, u_ode, semi, corr, smoothing_kernel) +end + +@inline function correction_matrix(system::EntropicallyDampedSPHSystem, particle) + extract_smatrix(system.cache.correction_matrix, system, particle) +end + function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fb49a65359..317aff3b9c 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -246,14 +246,27 @@ end return nothing end -@inline function surface_normal_method(system::AbstractFluidSystem) - return system.surface_normal_method +@inline function surface_method(system::AbstractFluidSystem) + hasproperty(system, :surface_method) || return nothing + return system.surface_method end -@inline function surface_normal_method(system) +@inline function surface_method(system) return nothing end +function surface_normal_method(system) + Base.depwarn("`surface_normal_method(system)` is deprecated; use `surface_method(system)`", + :surface_normal_method) + method = surface_method(system) + return computes_surface_normal(method) ? method : nothing +end + +@inline contributes_to_colorfield(system) = false +@inline function contributes_to_colorfield(system::AbstractFluidSystem) + return hasproperty(system.cache, :color) +end + function restart_u(system::AbstractFluidSystem, data) inactive_coords = convert(coordinates_eltype(system), 1e16) coords_total = fill(inactive_coords, u_nvariables(system), @@ -297,13 +310,13 @@ function restart_v(system::AbstractFluidSystem, data) end function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs) - if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension) + if requires_surface_normal(fluid_system.surface_tension) foreach_system(systems) do neighbor - if neighbor isa AbstractFluidSystem && - isnothing(fluid_system.surface_tension) && - isnothing(fluid_system.surface_normal_method) - throw(ArgumentError("either none or all fluid systems in a simulation need " * - "to use a surface tension model or a surface normal method.")) + if neighbor isa AbstractFluidSystem && !(neighbor isa ParticlePackingSystem) && + !computes_surface_normal(surface_method(neighbor)) + throw(ArgumentError("all interacting fluid systems must use a surface method " * + "that computes normals when a surface-tension model " * + "requires interface normals")) end end end diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index e4ef03584f..a4a8234316 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -4,7 +4,11 @@ viscosity=nothing, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), omega=0.5, max_error=0.1, min_iterations=2, - max_iterations=20, time_step) + max_iterations=20, time_step, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, + surface_method=nothing, + reference_particle_spacing=0.0, color_value=1) System for particles of a fluid. The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system @@ -30,9 +34,16 @@ See [Implicit Incompressible SPH](@ref iisph) for more details on the method. - `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 +- `correction`: Corrections are currently unsupported and passing a non-`nothing` + value throws an error. +- `density_correction`: Currently unsupported. +- `gradient_correction`: Currently unsupported. +- `surface_method`: Optional surface detection or normal method. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface calculations. """ struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, - IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS} + IC, K, V, PF, SM, C} <: AbstractFluidSystem{NDIMS} initial_condition :: IC mass :: ARRAY1D # Array{ELTYPE, 1} pressure :: ARRAY1D @@ -42,7 +53,7 @@ struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V pressure_acceleration_formulation :: PF - surface_normal_method :: Nothing # TODO + surface_method :: SM surface_tension :: Nothing # TODO particle_refinement :: Nothing # TODO density :: ARRAY1D @@ -71,10 +82,20 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, ndims(smoothing_kernel)), omega=0.5, max_error=0.1, min_iterations=2, max_iterations=20, time_step, - artificial_sound_speed=1000.0) + artificial_sound_speed=1000.0, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, + surface_method=nothing, + surface_normal_method=nothing, + reference_particle_spacing=0.0, color_value=1) particle_refinement = nothing # TODO surface_tension = nothing # TODO + if correction !== nothing || density_correction !== nothing || + gradient_correction !== nothing + throw(ArgumentError("corrections are not supported by `ImplicitIncompressibleSPHSystem`")) + end + NDIMS = ndims(initial_condition) ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) @@ -112,6 +133,12 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`time_step` must be a positive number")) end + surface_method = select_surface_method(surface_tension, surface_method, + surface_normal_method) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) + end + pressure_acceleration = pressure_acceleration_summation_density density = copy(initial_condition.density) @@ -124,13 +151,19 @@ function ImplicitIncompressibleSPHSystem(initial_condition; smoothing_kernel, density_error = zeros(ELTYPE, n_particles) cache = (; + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, - smoothing_length)...,) + smoothing_length)..., + color=Int(color_value)) + if reference_particle_spacing > 0 + cache = (; cache..., reference_particle_spacing) + end return ImplicitIncompressibleSPHSystem(initial_condition, mass, pressure, smoothing_kernel, smoothing_length, reference_density, acceleration_, viscosity, - pressure_acceleration, nothing, surface_tension, + pressure_acceleration, surface_method, + surface_tension, particle_refinement, density, predicted_density, advection_velocity, d_ii, a_ii, sum_d_ij_pj, sum_term, density_error, omega, max_error, @@ -216,6 +249,12 @@ function update_quantities!(system::ImplicitIncompressibleSPHSystem, v, u, semi) end +function update_pressure!(system::ImplicitIncompressibleSPHSystem, v, u, v_ode, u_ode, + semi, t) + compute_surface!(system, system.surface_method, v, u, v_ode, u_ode, semi, t) + return system +end + function update_implicit_sph!(semi, v_ode, u_ode, t) # This check is performed statically by the compiler and has no overhead if !any(system -> system isa ImplicitIncompressibleSPHSystem, semi.systems) diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index b6114c0bc1..4873143197 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -1,3 +1,67 @@ +abstract type AbstractSurfacePressure end + +""" + SurfacePressureDifference() + +Blend the configured conservative pressure acceleration with a first-order consistent +pressure-difference operator according to the target particle's [`surface_activity`](@ref). +Only interactions within the same fluid system are blended; boundaries and interactions +between different fluid systems remain conservative. + +This experimental model requires `GradientCorrection` or `MixedKernelGradientCorrection` +and a configured surface method. +""" +struct SurfacePressureDifference <: AbstractSurfacePressure end + +@inline function surface_pressure_model(system) + hasproperty(system, :surface_pressure) || return nothing + return system.surface_pressure +end + +function validate_surface_pressure(surface_pressure, surface_method_, gradient_correction) + isnothing(surface_pressure) && return surface_pressure + surface_pressure isa AbstractSurfacePressure || + throw(ArgumentError("`surface_pressure` must be an `AbstractSurfacePressure` or `nothing`")) + surface_method_ isa AbstractSurfaceMethod || + throw(ArgumentError("`SurfacePressureDifference` requires a configured surface method")) + if !(gradient_correction isa Union{GradientCorrection, MixedKernelGradientCorrection}) + throw(ArgumentError("`SurfacePressureDifference` currently supports only " * + "`GradientCorrection` and `MixedKernelGradientCorrection`")) + end + + return surface_pressure +end + +@inline function pressure_acceleration_difference(m_b, rho_a, rho_b, p_a, p_b, W_a) + return -m_b * div_fast(p_b - p_a, rho_a * rho_b) * W_a +end + +@inline function blend_surface_pressure(surface_pressure, conservative_acceleration, + activity, same_system, m_b, rho_a, rho_b, p_a, p_b, + W_a) + return conservative_acceleration +end + +@inline function blend_surface_pressure(::SurfacePressureDifference, + conservative_acceleration, activity, same_system, + m_b, rho_a, rho_b, p_a, p_b, W_a) + same_system || return conservative_acceleration + difference_acceleration = pressure_acceleration_difference(m_b, rho_a, rho_b, p_a, p_b, + W_a) + return (1 - activity) * conservative_acceleration + activity * difference_acceleration +end + +@inline function apply_surface_pressure(particle_system, neighbor_system, particle, + conservative_acceleration, m_b, rho_a, rho_b, p_a, + p_b, W_a) + model = surface_pressure_model(particle_system) + isnothing(model) && return conservative_acceleration + activity = surface_activity(particle_system, particle) + return blend_surface_pressure(model, conservative_acceleration, activity, + particle_system === neighbor_system, m_b, rho_a, rho_b, + p_a, p_b, W_a) +end + # As shown in "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic # formulations" by Bonet and Lok (1999), for a consistent formulation this form has to be # used with `SummationDensity`. @@ -93,6 +157,18 @@ end return -volume_term * pressure_tilde * W_a end +# Conservative extension for correction methods with asymmetric kernel gradients. This reduces +# to the symmetric formulation above when `W_b == -W_a`. +@inline function inter_particle_averaged_pressure(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a, + W_b) + volume_a = m_a / rho_a + volume_b = m_b / rho_b + volume_term = (volume_a^2 + volume_b^2) / m_a + pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b) + + return -0.5 * volume_term * pressure_tilde * (W_a - W_b) +end + function choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, correction) @@ -146,11 +222,19 @@ end # Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. @inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, - distance, W_a, correction) + distance, W_a, correction; + use_surface_pressure=true) # Without correction or with `AkinciFreeSurfaceCorrection`, the kernel gradient is # symmetric, so call the symmetric version of the pressure acceleration formulation. - return pressure_acceleration_formulation(particle_system)(m_a, m_b, rho_a, rho_b, - p_a, p_b, W_a) + conservative_acceleration = pressure_acceleration_formulation(particle_system)(m_a, m_b, + rho_a, + rho_b, + p_a, p_b, + W_a) + use_surface_pressure || return conservative_acceleration + return apply_surface_pressure(particle_system, neighbor_system, particle, + conservative_acceleration, m_b, rho_a, rho_b, p_a, p_b, + W_a) end # Formulation using asymmetric gradient formulation for corrections depending on local neighborhood. @@ -160,11 +244,19 @@ end correction::Union{KernelCorrection, GradientCorrection, BlendedGradientCorrection, - MixedKernelGradientCorrection}) + MixedKernelGradientCorrection}; + use_surface_pressure=true) W_b = smoothing_kernel_grad(neighbor_system, -pos_diff, distance, neighbor) # With correction, the kernel gradient is not necessarily symmetric, so call the # asymmetric version of the pressure acceleration formulation. - return pressure_acceleration_formulation(particle_system)(m_a, m_b, rho_a, rho_b, - p_a, p_b, W_a, W_b) + conservative_acceleration = pressure_acceleration_formulation(particle_system)(m_a, m_b, + rho_a, + rho_b, + p_a, p_b, + W_a, W_b) + use_surface_pressure || return conservative_acceleration + return apply_surface_pressure(particle_system, neighbor_system, particle, + conservative_acceleration, m_b, rho_a, rho_b, p_a, p_b, + W_a) end diff --git a/src/schemes/fluid/shifting_techniques.jl b/src/schemes/fluid/shifting_techniques.jl index 315f288c04..3833fd6d69 100644 --- a/src/schemes/fluid/shifting_techniques.jl +++ b/src/schemes/fluid/shifting_techniques.jl @@ -610,7 +610,8 @@ end # `pressure_acceleration` function. dv_particle[] += pressure_acceleration(system, neighbor_system, particle, neighbor, m_a, m_b, A_a, A_b, rho_a, rho_b, - pos_diff, distance, grad_kernel, correction) + pos_diff, distance, grad_kernel, correction; + use_surface_pressure=false) return dv_particle end @@ -707,7 +708,8 @@ function update_shifting!(system, shifting::TransportVelocityAdami, v, u, v_ode, pressure_acceleration(system, neighbor_system, particle, neighbor, m_a, m_b, 1, 1, rho_a, rho_b, pos_diff, distance, grad_kernel, - system_correction(system)) + system_correction(system); + use_surface_pressure=false) # Write into the buffer for i in eachindex(delta_v_) diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index adbc9d7dbe..5d394c762f 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -1,58 +1,235 @@ +abstract type AbstractSurfaceMethod end +abstract type AbstractSurfaceNormalMethod <: AbstractSurfaceMethod end + @doc raw""" - ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) + ColorfieldSurfaceDetection(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + +Detect fluid surfaces from the magnitude of a colorfield gradient. Different +`color_value`s detect represented fluid-fluid interfaces, while incomplete constant-color +support detects a free surface. This method computes [`surface_activity`](@ref), but does +not expose a surface normal. +""" +struct ColorfieldSurfaceDetection{ELTYPE} <: AbstractSurfaceMethod + boundary_contact_threshold::ELTYPE + interface_threshold::ELTYPE + ideal_density_threshold::ELTYPE + interface_taper_start::ELTYPE + interpolation_surface_threshold::ELTYPE +end -Color field based computation of the interface normals. +@doc raw""" + ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + +Compute colorfield surface normals and [`surface_activity`](@ref). The detection stage is +identical to [`ColorfieldSurfaceDetection`](@ref). The raw gradient is filtered after its +magnitude has been converted to activity and is normalized when required by the configured +surface-tension model. # Keywords -- `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary. -- `interface_threshold=0.01`: Threshold for normals to be removed as being invalid. -- `ideal_density_threshold=0.0`: Assume particles are inside if they are above this threshold, which is relative to the `ideal_neighbor_count`. +- `boundary_contact_threshold=0.1`: Finite value in `[0, 1]` used to detect contact with + dummy-particle boundaries. +- `interface_threshold=0.01`: Finite, non-negative dimensionless gradient threshold. +- `ideal_density_threshold=0.0`: Optional neighbor-count heuristic for unrepresented exterior + phases. Zero disables it; keep it disabled for fully supported multiphase interfaces. +- `interface_taper_start=0.8`: Start of the smooth activity transition as a fraction of + `interface_threshold`. +- `interpolation_surface_threshold=0.45`: Minimum normalized reference-color contribution + retained by interpolated output. """ -struct ColorfieldSurfaceNormal{ELTYPE} +struct ColorfieldSurfaceNormal{ELTYPE} <: AbstractSurfaceNormalMethod boundary_contact_threshold::ELTYPE interface_threshold::ELTYPE ideal_density_threshold::ELTYPE + interface_taper_start::ELTYPE + interpolation_surface_threshold::ELTYPE end -function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) - return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, +const ColorfieldSurfaceMethod = Union{ColorfieldSurfaceDetection, ColorfieldSurfaceNormal} + +function validate_surface_threshold(threshold, name; upper_bound=nothing, + strict_upper_bound=false) + interval = if isnothing(upper_bound) + "non-negative" + elseif strict_upper_bound + "in [0, $upper_bound)" + else + "in [0, $upper_bound]" + end + threshold isa Real || + throw(ArgumentError("`$name` must be a finite real number $interval")) + + valid_upper_bound = isnothing(upper_bound) || + (strict_upper_bound ? threshold < upper_bound : + threshold <= upper_bound) + if !isfinite(threshold) || threshold < 0 || !valid_upper_bound + throw(ArgumentError("`$name` must be a finite real number $interval")) + end + + return threshold +end + +function colorfield_surface_parameters(; boundary_contact_threshold=0.1, + interface_threshold=0.01, + ideal_density_threshold=0.0, + interface_taper_start=0.8, + interpolation_surface_threshold=0.45) + boundary_threshold = validate_surface_threshold(boundary_contact_threshold, + "boundary_contact_threshold"; + upper_bound=1) + normal_threshold = validate_surface_threshold(interface_threshold, + "interface_threshold") + density_threshold = validate_surface_threshold(ideal_density_threshold, + "ideal_density_threshold"; + upper_bound=1) + taper_start = validate_surface_threshold(interface_taper_start, + "interface_taper_start"; + upper_bound=1, + strict_upper_bound=true) + interpolation_threshold = validate_surface_threshold(interpolation_surface_threshold, + "interpolation_surface_threshold"; + upper_bound=1) + parameters = promote(boundary_threshold, normal_threshold, density_threshold, + taper_start, interpolation_threshold) + return eltype(parameters) <: Integer ? float.(parameters) : parameters +end + +function ColorfieldSurfaceDetection(; kwargs...) + return ColorfieldSurfaceDetection(colorfield_surface_parameters(; kwargs...)...) +end + +function ColorfieldSurfaceNormal(; kwargs...) + return ColorfieldSurfaceNormal(colorfield_surface_parameters(; kwargs...)...) +end + +function ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + return ColorfieldSurfaceNormal(; boundary_contact_threshold, interface_threshold, ideal_density_threshold) end -function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) +@inline computes_surface_normal(surface_method) = false +@inline computes_surface_normal(::AbstractSurfaceNormalMethod) = true + +@inline is_colorfield_surface_method(surface_method) = false +@inline is_colorfield_surface_method(::ColorfieldSurfaceMethod) = true + +@inline contributes_boundary_colorfield(system) = false +@inline contributes_boundary_colorfield(::AbstractBoundarySystem) = true + +@inline function default_surface_method(surface_tension, surface_method) + if isnothing(surface_method) && requires_surface_normal(surface_tension) + return ColorfieldSurfaceNormal() + end + + return surface_method +end + +function select_surface_method(surface_tension, surface_method, surface_normal_method) + if !isnothing(surface_method) && !isnothing(surface_normal_method) + throw(ArgumentError("`surface_method` and deprecated `surface_normal_method` cannot both be set")) + end + + if !isnothing(surface_normal_method) + Base.depwarn("`surface_normal_method` is deprecated; use `surface_method` instead", + :surface_normal_method) + surface_method = surface_normal_method + end + + surface_method = default_surface_method(surface_tension, surface_method) + if !(surface_method isa Union{Nothing, AbstractSurfaceMethod}) + throw(ArgumentError("`surface_method` must be an `AbstractSurfaceMethod` or `nothing`")) + end + if requires_surface_normal(surface_tension) && !computes_surface_normal(surface_method) + throw(ArgumentError("$(typeof(surface_tension)) requires a surface method that computes surface normals")) + end + + return surface_method +end + +@inline function cubic_smoothstep(value) + value <= zero(value) && return zero(value) + value >= one(value) && return one(value) + return value^2 * (3 - 2value) +end + +@inline function gradient_surface_activity(normal_norm, support_radius, + surface_method::ColorfieldSurfaceMethod) + threshold = surface_method.interface_threshold + dimensionless_norm = support_radius * normal_norm + if iszero(threshold) + return iszero(dimensionless_norm) ? zero(dimensionless_norm) : + one(dimensionless_norm) + end + + lower_bound = surface_method.interface_taper_start * threshold + transition_coordinate = (dimensionless_norm - lower_bound) / + (threshold - lower_bound) + return cubic_smoothstep(transition_coordinate) +end + +function create_cache_surface(surface_method, ELTYPE, NDIMS, nparticles) return (;) end -function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) +function create_cache_surface(::ColorfieldSurfaceDetection, ELTYPE, NDIMS, nparticles) + surface_gradient = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + surface_activity = Array{ELTYPE, 1}(undef, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + colorfield = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_gradient, surface_activity, neighbor_count, colorfield) +end + +function create_cache_surface(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, nparticles) surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + surface_activity = Array{ELTYPE, 1}(undef, nparticles) neighbor_count = Array{ELTYPE, 1}(undef, nparticles) colorfield = Array{ELTYPE, 1}(undef, nparticles) correction_factor = Array{ELTYPE, 1}(undef, nparticles) - return (; surface_normal, neighbor_count, colorfield, correction_factor) + return (; surface_normal, surface_activity, neighbor_count, colorfield, + correction_factor) +end + +@inline function surface_gradient(cache, ::ColorfieldSurfaceDetection) + return cache.surface_gradient +end + +@inline function surface_gradient(cache, ::ColorfieldSurfaceNormal) + return cache.surface_normal end @inline function surface_normal(particle_system::AbstractFluidSystem, particle) - (; cache) = particle_system - return extract_svector(cache.surface_normal, particle_system, particle) + return extract_svector(particle_system.cache.surface_normal, particle_system, particle) +end + +@inline function surface_activity(particle_system::AbstractFluidSystem, particle) + return @inbounds particle_system.cache.surface_activity[particle] end -function calc_normal!(system, neighbor_system, u_system, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method, - neighbor_surface_normal_method) - # Normal not needed +function calc_surface!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surface_method, neighbor_surface_method) return system end # Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" # and Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics". -function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFluidSystem, - u_system, v, - v_neighbor_system, u_neighbor_system, semi, surface_normal_method, - ::ColorfieldSurfaceNormal) - (; cache) = system +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod, + neighbor_surface_method) + contributes_to_colorfield(neighbor_system) || return system + (; cache) = system + gradient = surface_gradient(cache, surface_method) + color_b = neighbor_system.cache.color system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) @@ -60,61 +237,53 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui system_coords, neighbor_system_coords, semi; points=each_integrated_particle(system)) do particle, neighbor, pos_diff, distance - m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = current_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + volume_b = hydrodynamic_mass(neighbor_system, neighbor) / + current_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = kernel_grad(system_smoothing_kernel(system), pos_diff, distance, + smoothing_length(system, particle)) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + gradient[i, particle] += volume_b * color_b * grad_kernel[i] end - cache.neighbor_count[particle] += 1 end return system end -# Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" -# Note: This is the simplest form of normal approximation commonly used in SPH and comes -# with serious deficits in accuracy especially at corners, small neighborhoods and boundaries -function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_system, v, - u_neighbor_system, semi, surface_normal_method) +function calc_boundary_surface!(system::AbstractFluidSystem, neighbor_system, u_system, v, + u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod) (; cache) = system + gradient = surface_gradient(cache, surface_method) (; colorfield, initial_colorfield) = neighbor_system.boundary_model.cache - (; boundary_contact_threshold) = surface_normal_method + (; boundary_contact_threshold) = surface_method + color_a = system.cache.color system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - # First we need to calculate the smoothed colorfield values of the boundary - # TODO: move colorfield to extra step - # TODO: this is only correct for a single fluid - - # Reset to the constant boundary interpolated color values - colorfield .= initial_colorfield - - # Accumulate fluid neighbors + colorfield .= abs.(initial_colorfield) foreach_point_neighbor(neighbor_system, system, neighbor_system_coords, system_coords, semi) do particle, neighbor, pos_diff, distance colorfield[particle] += hydrodynamic_mass(system, neighbor) / - current_density(v, system, neighbor) * system.cache.color * + current_density(v, system, neighbor) * abs(color_a) * smoothing_kernel(system, distance, particle) end maximum_colorfield = maximum(colorfield) + iszero(maximum_colorfield) && return system foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_system_coords, semi) do particle, neighbor, pos_diff, distance - # We assume that we are in contact with the boundary if the color of the boundary particle - # is larger than the threshold if colorfield[neighbor] / maximum_colorfield > boundary_contact_threshold - m_b = hydrodynamic_mass(system, particle) - density_neighbor = current_density(v, system, particle) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + volume_a = hydrodynamic_mass(system, particle) / + current_density(v, system, particle) + grad_kernel = kernel_grad(system_smoothing_kernel(system), pos_diff, distance, + smoothing_length(system, particle)) for i in 1:ndims(system) - cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] + gradient[i, particle] += volume_a * color_a * grad_kernel[i] end cache.neighbor_count[particle] += 1 end @@ -123,114 +292,128 @@ function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_s return system end -function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractBoundarySystem, - u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method, neighbor_surface_normal_method) - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::AbstractBoundarySystem, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method::ColorfieldSurfaceMethod, + neighbor_surface_method) + return calc_boundary_surface!(system, neighbor_system, u_system, v, u_neighbor_system, + semi, surface_method) end -function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, - surface_normal_method) - (; cache) = system +function invalid_surface_particle(system, surface_method::ColorfieldSurfaceMethod, + particle, support_radius) + neighbor_count = system.cache.neighbor_count[particle] + minimum_neighbor_count = 2^ndims(system) + 1 + neighbor_count < minimum_neighbor_count && return true + + threshold = surface_method.ideal_density_threshold + return threshold > 0 && + threshold * ideal_neighbor_count(Val(ndims(system)), + system.cache.reference_particle_spacing, + support_radius) < neighbor_count +end + +function finalize_surface!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceMethod) + gradient = surface_gradient(system.cache, surface_method) + support_radius = compact_support(system_smoothing_kernel(system), + initial_smoothing_length(system)) - # We remove invalid normals (too few neighbors) to reduce the impact of underdefined normals for particle in each_integrated_particle(system) - # A corner has that many neighbors assuming a regular 2 * r distribution and a compact_support of 4r - if cache.neighbor_count[particle] < 2^ndims(system) + 1 - cache.surface_normal[1:ndims(system), particle] .= 0 + particle_gradient = extract_svector(gradient, system, particle) + gradient_norm = norm(particle_gradient) + activity = gradient_surface_activity(gradient_norm, support_radius, surface_method) + + if invalid_surface_particle(system, surface_method, particle, support_radius) + system.cache.surface_activity[particle] = zero(activity) + gradient[1:ndims(system), particle] .= 0 + else + system.cache.surface_activity[particle] = activity end end return system end -# See Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" -function remove_invalid_normals!(system::AbstractFluidSystem, - surface_tension::Union{SurfaceTensionMorris, - SurfaceTensionMomentumMorris}, - surface_normal_method::ColorfieldSurfaceNormal) - (; cache, smoothing_kernel) = system - (; ideal_density_threshold, interface_threshold) = surface_normal_method - (; neighbor_count) = cache - - smoothing_length_ = initial_smoothing_length(system) - - # We remove invalid normals i.e. they have a small norm (eq. 20) - normal_condition2 = (interface_threshold / - compact_support(smoothing_kernel, smoothing_length_))^2 +function finalize_surface!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceNormal) + gradient = surface_gradient(system.cache, surface_method) + support_radius = compact_support(system_smoothing_kernel(system), + initial_smoothing_length(system)) + normal_condition2 = (surface_method.interface_threshold / support_radius)^2 for particle in each_integrated_particle(system) - - # Heuristic condition if there is no gas phase to find the free surface. - # We remove normals for particles which have a lot of support e.g. they are in the interior. - if ideal_density_threshold > 0 && - ideal_density_threshold * - ideal_neighbor_count(Val(ndims(system)), cache.reference_particle_spacing, - compact_support(smoothing_kernel, smoothing_length_)) < - neighbor_count[particle] - cache.surface_normal[1:ndims(system), particle] .= 0 - continue - end - - particle_surface_normal = surface_normal(system, particle) - norm2 = dot(particle_surface_normal, particle_surface_normal) - - # See eq. 21 - if norm2 > normal_condition2 - cache.surface_normal[1:ndims(system), - particle] = particle_surface_normal / sqrt(norm2) + particle_gradient = extract_svector(gradient, system, particle) + norm2 = dot(particle_gradient, particle_gradient) + gradient_norm = sqrt(norm2) + activity = gradient_surface_activity(gradient_norm, support_radius, surface_method) + + if invalid_surface_particle(system, surface_method, particle, support_radius) + system.cache.surface_activity[particle] = zero(activity) + gradient[1:ndims(system), particle] .= 0 + elseif norm2 > normal_condition2 + system.cache.surface_activity[particle] = activity + if normalize_surface_normals(surface_tension) + gradient[1:ndims(system), particle] = particle_gradient / gradient_norm + end else - cache.surface_normal[1:ndims(system), particle] .= 0 + system.cache.surface_activity[particle] = activity + gradient[1:ndims(system), particle] .= 0 end end return system end -function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) +@inline normalize_surface_normals(surface_tension) = false +@inline normalize_surface_normals(::SurfaceTensionMorris) = true +@inline normalize_surface_normals(::SurfaceTensionMomentumMorris) = true + +function compute_surface!(system, surface_method, v, u, v_ode, u_ode, semi, t) return system end -function compute_surface_normal!(system::AbstractFluidSystem, - surface_normal_method_::ColorfieldSurfaceNormal, - v, u, v_ode, u_ode, semi, t) +function compute_surface!(system::AbstractFluidSystem, + surface_method_::ColorfieldSurfaceMethod, + v, u, v_ode, u_ode, semi, t) (; cache, surface_tension) = system - # Reset surface normal - set_zero!(cache.surface_normal) + set_zero!(surface_gradient(cache, surface_method_)) + set_zero!(cache.surface_activity) set_zero!(cache.neighbor_count) - # TODO: if color values are set only different systems need to be called - @trixi_timeit timer() "compute surface normal" begin + @trixi_timeit timer() "compute surface" begin foreach_system_wrapped(semi, v_ode, u_ode) do neighbor_system, v_neighbor_system, u_neighbor_system - if !has_system_interaction(system, neighbor_system, semi) - # No interaction between these systems. - return - end + has_system_interaction(system, neighbor_system, semi) || return - calc_normal!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method_, - surface_normal_method(neighbor_system)) + calc_surface!(system, neighbor_system, u, v, v_neighbor_system, + u_neighbor_system, semi, surface_method_, + surface_method(neighbor_system)) end end - remove_invalid_normals!(system, surface_tension, surface_normal_method_) + finalize_surface!(system, surface_tension, surface_method_) return system end +function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, + surface_method::ColorfieldSurfaceNormal) + return finalize_surface!(system, surface_tension, surface_method) +end + function calc_curvature!(system, neighbor_system, u_system, v, - v_neighbor_system, u_neighbor_system, semi, surface_normal_method, - neighbor_surface_normal_method) + v_neighbor_system, u_neighbor_system, semi, surface_method, + neighbor_surface_method) end # Section 5 in Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractFluidSystem, u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method::ColorfieldSurfaceNormal, - neighbor_surface_normal_method::ColorfieldSurfaceNormal) + surface_method_::ColorfieldSurfaceNormal, + neighbor_surface_method::ColorfieldSurfaceNormal) (; cache) = system (; curvature, correction_factor) = cache @@ -248,7 +431,6 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF n_b = surface_normal(neighbor_system, neighbor) v_b = m_b / rho_b - # Eq. 22: we can test against `eps()` here since the surface normals that are invalid have been removed if dot(n_a, n_a) > eps() && dot(n_b, n_b) > eps() w = smoothing_kernel(system, distance, particle) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) @@ -256,12 +438,10 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF for i in 1:ndims(system) curvature[particle] += v_b * (n_b[i] - n_a[i]) * grad_kernel[i] end - # Eq. 24 correction_factor[particle] += v_b * w end end - # Eq. 23 for particle in each_integrated_particle(system) curvature[particle] /= (correction_factor[particle] + eps()) end @@ -276,23 +456,19 @@ end function compute_curvature!(system::AbstractFluidSystem, surface_tension::SurfaceTensionMorris, v, u, v_ode, u_ode, semi, t) - (; cache, surface_tension) = system + (; cache) = system - # Reset surface curvature set_zero!(cache.curvature) @trixi_timeit timer() "compute surface curvature" begin foreach_system_wrapped(semi, v_ode, u_ode) do neighbor_system, v_neighbor_system, u_neighbor_system - if !has_system_interaction(system, neighbor_system, semi) - # No interaction between these systems. - return - end + has_system_interaction(system, neighbor_system, semi) || return calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method(system), - surface_normal_method(neighbor_system)) + u_neighbor_system, semi, surface_method(system), + surface_method(neighbor_system)) end end return system diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..3e0682328d 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -98,6 +98,11 @@ struct SurfaceTensionMomentumMorris{ELTYPE} <: AbstractSurfaceTension end end +@inline requires_surface_normal(surface_tension) = false +@inline requires_surface_normal(::SurfaceTensionAkinci) = true +@inline requires_surface_normal(::SurfaceTensionMorris) = true +@inline requires_surface_normal(::SurfaceTensionMomentumMorris) = true + function create_cache_surface_tension(::SurfaceTensionMomentumMorris, ELTYPE, NDIMS, nparticles) delta_s = Array{ELTYPE, 1}(undef, nparticles) diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 8ceb29d27f..e63fabbace 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -8,6 +8,8 @@ function interact!(dv, v_particle_system, u_particle_system, eachparticle=each_integrated_particle(particle_system), kwargs...) (; density_calculator, correction) = particle_system + gradient_correction = correction_gradient(correction) + force_correction = correction_force(correction) sound_speed = system_sound_speed(particle_system) @@ -72,7 +74,7 @@ function interact!(dv, v_particle_system, u_particle_system, # Determine correction factors. # This can usually be ignored, as these are all 1 when no correction is used. (viscosity_correction, pressure_correction, - surface_tension_correction) = free_surface_correction(correction, + surface_tension_correction) = free_surface_correction(force_correction, particle_system, rho_a, rho_b) @@ -81,7 +83,7 @@ function interact!(dv, v_particle_system, u_particle_system, 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, correction) + distance, grad_kernel, gradient_correction) dv_particle[] += dv_pressure * pressure_correction # Propagate `@inbounds` to the viscosity function, which accesses particle data @@ -96,7 +98,7 @@ function interact!(dv, v_particle_system, u_particle_system, particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, m_a, m_b, rho_a, rho_b, v_a, v_b, - pos_diff, distance, grad_kernel, correction) + pos_diff, distance, grad_kernel, gradient_correction) @inbounds surface_tension_force!(dv_particle, surface_tension_a, surface_tension_b, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eea0607d7d..c6b4c3a389 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -6,8 +6,10 @@ pressure_acceleration=nothing, shifting_technique=nothing, buffer_size=nothing, - correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing, + surface_tension=nothing, surface_method=nothing, + surface_pressure=nothing, reference_particle_spacing=0.0, color_value=1)) System for particles of a fluid. @@ -41,7 +43,11 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. with this system. Default is no shifting. - `buffer_size`: Number of buffer particles. This is needed when simulating with [`OpenBoundarySystem`](@ref). -- `correction`: Correction method used for this system. (default: no correction, see [Corrections](@ref corrections)) +- `correction`: Legacy keyword for configuring one correction method. Cannot be + combined with `density_correction` or `gradient_correction`. +- `density_correction`: Density correction method. Currently supports + [`ShepardKernelCorrection`](@ref) with [`SummationDensity`](@ref). +- `gradient_correction`: Gradient correction method. See [Corrections](@ref corrections). - `source_terms`: Additional source terms for this system. Has to be either `nothing` (by default), or a function of `(coords, velocity, density, pressure, t)` (which are the quantities of a single particle), returning a `Tuple` @@ -53,17 +59,19 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. The keyword argument `acceleration` should be used instead for gravity-like source terms. - `surface_tension`: Surface tension model used for this SPH system. (default: no surface tension) -- `surface_normal_method`: The surface normal method to be used for this SPH system. - (default: no surface normal method or `ColorfieldSurfaceNormal()` if a surface_tension model is used) -- `reference_particle_spacing`: The reference particle spacing used for weighting values at the boundary, - which currently is only needed when using surface tension. -- `color_value`: Integer label used for calculation of surface normals. - Currently this is only used together with [`BoundaryModelDummyParticles`](@ref) and - [`ColorfieldSurfaceNormal`](@ref): fluid-boundary normal evaluation - reads the resulting boundary colorfield to detect wall contact. +- `surface_method`: Surface detection or normal method used by this system. + Methods that compute normals always also compute surface + activity. The default is `nothing`, or + `ColorfieldSurfaceNormal()` when required by surface tension. +- `surface_pressure`: Optional interface-aware pressure model. The experimental + [`SurfacePressureDifference`](@ref) blends the conservative + pressure acceleration with a pressure-difference operator. +- `reference_particle_spacing`: Reference spacing required by colorfield surface methods. +- `color_value`: Scalar contributed to colorfield surface detection. Different + values identify represented fluid-fluid interfaces. """ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, - PF, SC, ST, B, SRFT, SRFN, PR, + PF, SP, SC, ST, B, SRFT, SRFN, PR, C} <: AbstractFluidSystem{NDIMS} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} @@ -76,10 +84,11 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, density_diffusion :: DD correction :: COR pressure_acceleration_formulation :: PF + surface_pressure :: SP shifting_technique :: SC source_terms :: ST surface_tension :: SRFT - surface_normal_method :: SRFN + surface_method :: SRFN buffer :: B particle_refinement :: PR # TODO cache :: C @@ -96,8 +105,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, pressure_acceleration=nothing, shifting_technique=nothing, buffer_size=nothing, - correction=nothing, source_terms=nothing, - surface_tension=nothing, surface_normal_method=nothing, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing, + surface_tension=nothing, surface_method=nothing, + surface_normal_method=nothing, + surface_pressure=nothing, reference_particle_spacing=0, color_value=1) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -112,6 +124,11 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, ELTYPE = eltype(initial_condition) n_particles = nparticles(initial_condition) + correction = resolve_correction_configuration(correction, density_correction, + gradient_correction) + density_correction_ = correction_density(correction) + gradient_correction_ = correction_gradient(correction) + mass = copy(initial_condition.mass) pressure = similar(initial_condition.pressure) @@ -125,29 +142,29 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end - if correction isa ShepardKernelCorrection && + if density_correction_ isa ShepardKernelCorrection && density_calculator isa ContinuityDensity throw(ArgumentError("`ShepardKernelCorrection` cannot be used with `ContinuityDensity`")) end - if surface_tension !== nothing && surface_normal_method === nothing - surface_normal_method = ColorfieldSurfaceNormal() - end + surface_method = select_surface_method(surface_tension, surface_method, + surface_normal_method) + surface_pressure = validate_surface_pressure(surface_pressure, surface_method, + gradient_correction_) - if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + if is_colorfield_surface_method(surface_method) && reference_particle_spacing < eps() + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a colorfield surface method")) end pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, density_calculator, NDIMS, ELTYPE, - correction) + gradient_correction_) cache = (; create_cache_density(initial_condition, density_calculator)..., create_cache_correction(correction, initial_condition.density, NDIMS, n_particles)..., - create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, - n_particles)..., + create_cache_surface(surface_method, ELTYPE, NDIMS, n_particles)..., create_cache_surface_tension(surface_tension, ELTYPE, NDIMS, n_particles)..., create_cache_refinement(initial_condition, particle_refinement, @@ -168,8 +185,9 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, density_calculator, state_equation, smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, - shifting_technique, source_terms, surface_tension, - surface_normal_method, buffer, particle_refinement, + surface_pressure, shifting_technique, source_terms, + surface_tension, + surface_method, buffer, particle_refinement, cache) end @@ -185,8 +203,9 @@ function Base.show(io::IO, system::WeaklyCompressibleSPHSystem) print(io, ", ", system.density_diffusion) print(io, ", ", system.shifting_technique) print(io, ", ", system.surface_tension) - print(io, ", ", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal + print(io, ", ", system.surface_method) + print(io, ", ", system.surface_pressure) + if is_colorfield_surface_method(system.surface_method) print(io, ", ", system.cache.color) end print(io, ", ", system.acceleration) @@ -217,8 +236,9 @@ function Base.show(io::IO, ::MIME"text/plain", system::WeaklyCompressibleSPHSyst summary_line(io, "density diffusion", system.density_diffusion) summary_line(io, "shifting technique", system.shifting_technique) summary_line(io, "surface tension", system.surface_tension) - summary_line(io, "surface normal method", system.surface_normal_method) - if system.surface_normal_method isa ColorfieldSurfaceNormal + summary_line(io, "surface method", system.surface_method) + summary_line(io, "surface pressure", system.surface_pressure) + if is_colorfield_surface_method(system.surface_method) summary_line(io, "color", system.cache.color) end summary_line(io, "acceleration", system.acceleration) @@ -243,7 +263,9 @@ end @inline buffer(system::WeaklyCompressibleSPHSystem) = system.buffer -system_correction(system::WeaklyCompressibleSPHSystem) = system.correction +function system_correction(system::WeaklyCompressibleSPHSystem) + correction_gradient(system.correction) +end @propagate_inbounds function current_velocity(v, system::WeaklyCompressibleSPHSystem) return current_velocity(v, system.density_calculator, system) @@ -321,19 +343,20 @@ end end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) - (; density_calculator, correction, surface_normal_method, surface_tension) = system - - compute_pressure!(system, v, semi) - - # These are only computed when using corrections - compute_correction_values!(system, correction, u, v_ode, u_ode, semi) - compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) - # `kernel_correct_density!` only performed for `SummationDensity` - kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, + (; density_calculator, correction, surface_method, surface_tension) = system + density_correction = correction_density(correction) + gradient_correction = correction_gradient(correction) + + # Density correction must be applied before assembling gradient corrections so all + # gradient moments use the density that will be used by the RHS. + compute_correction_values!(system, density_correction, u, v_ode, u_ode, semi) + kernel_correct_density!(system, v, u, v_ode, u_ode, semi, density_correction, density_calculator) + compute_correction_values!(system, gradient_correction, u, v_ode, u_ode, semi) + compute_gradient_correction_matrix!(gradient_correction, system, u, v_ode, u_ode, semi) + compute_pressure!(system, v, semi) - # These are only computed when using surface tension - compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) + compute_surface!(system, surface_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, surface_tension, semi) return system end @@ -369,13 +392,13 @@ function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, MixedKernelGradientCorrection}, system::WeaklyCompressibleSPHSystem, u, v_ode, u_ode, semi) - (; cache, correction, smoothing_kernel) = system + (; cache, smoothing_kernel) = system (; correction_matrix) = cache system_coords = current_coordinates(u, system) compute_gradient_correction_matrix!(correction_matrix, system, system_coords, - v_ode, u_ode, semi, correction, smoothing_kernel) + v_ode, u_ode, semi, corr, smoothing_kernel) end function reinit_density!(vu_ode, semi) @@ -397,14 +420,16 @@ end function reinit_density!(system::WeaklyCompressibleSPHSystem, ::ContinuityDensity, v, u, v_ode, u_ode, semi) + # Use the independently evolved density to determine particle volumes before replacing it + # with the reinitialized summation density. + kernel_correction_coefficient = similar(v, size(v, 2)) + compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi, + kernel_correction_coefficient) + # Compute density with `SummationDensity` and store the result in `v`, # overwriting the previous integrated density. summation_density!(system, semi, u, u_ode, v[end, :]) - # Apply `ShepardKernelCorrection` - kernel_correction_coefficient = zeros(size(v[end, :])) - compute_shepard_coeff!(system, current_coordinates(u, system), v_ode, u_ode, semi, - kernel_correction_coefficient) @threaded semi for particle in eachparticle(system) v[end, particle] /= kernel_correction_coefficient[particle] end diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index dccacc2668..c6e9fe7c22 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -32,9 +32,9 @@ torque and applied consistently to all rigid particles. only evaluated for fluid-structure interaction with surface-tension-enabled fluid systems. - `color_value`: Integer label stored as `system.cache.color`. - Currently this is used with `BoundaryModelDummyParticles` during - colorfield initialization so fluids using - [`ColorfieldSurfaceNormal`](@ref) can detect contact with rigid + Currently this is used with `BoundaryModelDummyParticles` during + colorfield initialization so fluids using a colorfield surface method + can detect contact with rigid bodies, it participates in the multi-system color sanity check for surface-tension setups, and it is written to VTK output as `"color"`. """ @@ -267,7 +267,7 @@ end end @inline function system_correction(system::RigidBodySystem{<:BoundaryModelDummyParticles}) - return system.boundary_model.correction + return correction_gradient(system.boundary_model.correction) end function initialize!(system::RigidBodySystem, semi) @@ -275,14 +275,18 @@ function initialize!(system::RigidBodySystem, semi) return system end -function calc_normal!(system::AbstractFluidSystem, - neighbor_system::RigidBodySystem{<:BoundaryModelDummyParticles}, - u_system, v, v_neighbor_system, u_neighbor_system, semi, - surface_normal_method, neighbor_surface_normal_method) +function calc_surface!(system::AbstractFluidSystem, + neighbor_system::RigidBodySystem{<:BoundaryModelDummyParticles}, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + surface_method_, neighbor_surface_method) haskey(neighbor_system.boundary_model.cache, :initial_colorfield) || return system - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) + return calc_boundary_surface!(system, neighbor_system, u_system, v, u_neighbor_system, + semi, surface_method_) +end + +@inline function contributes_boundary_colorfield(::RigidBodySystem{<:BoundaryModelDummyParticles}) + return true end @inline function adhesion_force!(dv_particle, @@ -631,10 +635,10 @@ function check_configuration(system::RigidBodySystem, systems, nhs) end if neighbor isa AbstractFluidSystem && - neighbor.surface_normal_method isa ColorfieldSurfaceNormal + is_colorfield_surface_method(surface_method(neighbor)) if !(boundary_model isa BoundaryModelDummyParticles) throw(ArgumentError("`RigidBodySystem` is only compatible with " * - "`ColorfieldSurfaceNormal` when using " * + "colorfield surface methods when using " * "`BoundaryModelDummyParticles`.")) end @@ -642,7 +646,7 @@ function check_configuration(system::RigidBodySystem, systems, nhs) throw(ArgumentError("`RigidBodySystem` with `BoundaryModelDummyParticles` " * "requires `reference_particle_spacing` to be set on " * "the boundary model when used together with " * - "`ColorfieldSurfaceNormal` or a surface tension model.")) + "a colorfield surface method or a surface tension model.")) end end end diff --git a/test/examples/dam_break_2d_corrections.jl b/test/examples/dam_break_2d_corrections.jl index 3e70458b81..ea53b08d21 100644 --- a/test/examples/dam_break_2d_corrections.jl +++ b/test/examples/dam_break_2d_corrections.jl @@ -76,6 +76,34 @@ @test count_rhs_allocations(sol) == 0 end + @testset verbose=true "surface pressure $correction_name" for (correction_name, + correction) in + (("gradient", + GradientCorrection()), + ("mixed", + MixedKernelGradientCorrection())) + @trixi_test_nowarn trixi_include(@__MODULE__, + joinpath(examples_dir(), "fluid", + "dam_break_2d.jl"); + fluid_particle_spacing=particle_spacing, + smoothing_length=2 * particle_spacing, + smoothing_kernel=WendlandC6Kernel{2}(), + boundary_density_calculator=SummationDensity(), + fluid_density_calculator=SummationDensity(), + correction, + surface_method=ColorfieldSurfaceDetection(ideal_density_threshold=0.9), + surface_pressure=SurfacePressureDifference(), + use_reinit=false, clip_negative_pressure=true, + prefix="surface_pressure_$(correction_name)", + tspan, + fluid_density, density_diffusion=nothing, + boundary_layers=5, sol=nothing) + + sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol) == 0 + end + @testset verbose=true "$correction_name" for correction_name in keys(correction_dict) local fluid_density_calculator = density_calculator_dict[correction_name] local correction = correction_dict[correction_name] diff --git a/test/general/corrections.jl b/test/general/corrections.jl new file mode 100644 index 0000000000..7ed9b44d3a --- /dev/null +++ b/test/general/corrections.jl @@ -0,0 +1,652 @@ +@trixi_testset "Correction Consistency" begin + function correction_setup(correction=nothing; n=9, perturbation=false, + density_calculator=ContinuityDensity(), edac=false, + density_correction=nothing, gradient_correction=nothing, + pressure_acceleration=:default, + velocity=(pos -> SVector(pos[1], pos[2]))) + particle_spacing = 1.0 / n + smoothing_length = 2.0 * particle_spacing + smoothing_kernel = WendlandC6Kernel{2}() + fluid = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); + density=1000.0, velocity, + coordinates_perturbation=perturbation ? 0.1 : nothing) + + if edac + if pressure_acceleration === :default + system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel, + smoothing_length, sound_speed=10.0, + density_calculator, correction, + density_correction, + gradient_correction) + else + system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel, + smoothing_length, sound_speed=10.0, + density_calculator, correction, + density_correction, + gradient_correction, + pressure_acceleration) + end + else + state_equation = StateEquationCole(; sound_speed=10.0, + reference_density=1000.0, exponent=1) + if pressure_acceleration === :default + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, + smoothing_length, density_calculator, + state_equation, correction, + density_correction, + gradient_correction) + else + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, + smoothing_length, density_calculator, + state_equation, correction, + density_correction, + gradient_correction, + pressure_acceleration) + end + end + + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + semi = ode.p.semi + system = first(semi.systems) + + return (; system, semi, v_ode, u_ode, particle_spacing) + end + + function fill_correction_cache!(system, value) + for name in (:kernel_correction_coefficient, :dw_gamma, :correction_matrix) + hasproperty(system.cache, name) || continue + fill!(getproperty(system.cache, name), value) + end + return system + end + + function update_correction!(setup) + (; system, semi, v_ode, u_ode) = setup + fill_correction_cache!(system, NaN) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + return setup + end + + function correction_moments(setup; field=(pos -> 1.0)) + (; system, semi, v_ode, u_ode) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + values = [field(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + n_particles = TrixiParticles.nparticles(system) + + zeroth_gradient_moment = zeros(2, n_particles) + first_gradient_moment = zeros(2, 2, n_particles) + direct_gradient = zeros(2, n_particles) + difference_gradient = zeros(2, n_particles) + + GC.@preserve v_ode u_ode begin + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + semi) do particle, neighbor, pos_diff, + distance + pos_diff_ = SVector(pos_diff) + volume = TrixiParticles.hydrodynamic_mass(system, neighbor) / + TrixiParticles.current_density(v, system, neighbor) + gradient = TrixiParticles.smoothing_kernel_grad(system, pos_diff_, distance, + particle) + neighbor_offset = -pos_diff_ + + for i in 1:2 + zeroth_gradient_moment[i, particle] += volume * gradient[i] + direct_gradient[i, particle] += volume * values[neighbor] * gradient[i] + difference_gradient[i, + particle] += volume * + (values[neighbor] - values[particle]) * + gradient[i] + for j in 1:2 + first_gradient_moment[i, j, + particle] += volume * gradient[i] * + neighbor_offset[j] + end + end + end + end + + return (; zeroth_gradient_moment, first_gradient_moment, direct_gradient, + difference_gradient) + end + + function corner_particle(system) + coordinates = TrixiParticles.initial_coordinates(system) + return argmin(eachindex(axes(coordinates, 2))) do particle + coordinates[1, particle] + coordinates[2, particle] + end + end + + @testset "Cache lifecycle" begin + corrections = (KernelCorrection(), GradientCorrection(), + BlendedGradientCorrection(0.4), MixedKernelGradientCorrection()) + + for edac in (false, true), correction in corrections + setup = correction_setup(correction; edac) + update_correction!(setup) + + for name in (:kernel_correction_coefficient, :dw_gamma, :correction_matrix) + hasproperty(setup.system.cache, name) || continue + @test all(isfinite, getproperty(setup.system.cache, name)) + end + end + + setup = correction_setup(ShepardKernelCorrection(); + density_calculator=SummationDensity()) + update_correction!(setup) + density = TrixiParticles.current_density(TrixiParticles.wrap_v(setup.v_ode, + setup.system, + setup.semi), + setup.system) + @test setup.system.pressure ≈ setup.system.state_equation.(density) + + setup_edac = correction_setup(ShepardKernelCorrection(); + density_calculator=SummationDensity(), edac=true) + update_correction!(setup_edac) + @test all(isfinite, setup_edac.system.cache.kernel_correction_coefficient) + @test all(isfinite, setup_edac.system.cache.density) + + for edac in (false, true) + combined = correction_setup(; density_calculator=SummationDensity(), edac, + density_correction=ShepardKernelCorrection(), + gradient_correction=MixedKernelGradientCorrection()) + update_correction!(combined) + @test combined.system.correction isa CorrectionConfiguration + @test all(isfinite, combined.system.cache.kernel_correction_coefficient) + @test all(isfinite, combined.system.cache.dw_gamma) + @test all(isfinite, combined.system.cache.correction_matrix) + end + + @test_throws ArgumentError correction_setup(GradientCorrection(); + gradient_correction=GradientCorrection()) + @test_throws ArgumentError correction_setup(; + density_calculator=ContinuityDensity(), + density_correction=ShepardKernelCorrection()) + @test_throws ArgumentError CorrectionConfiguration(; density=GradientCorrection()) + @test_throws ArgumentError CorrectionConfiguration(; + gradient=ShepardKernelCorrection()) + @test_throws ArgumentError BlendedGradientCorrection(-0.1) + @test_throws ArgumentError BlendedGradientCorrection(1.1) + + iisph_particles = RectangularShape(0.1, (2, 2), (0.0, 0.0); density=1000.0) + @test_throws ArgumentError ImplicitIncompressibleSPHSystem(iisph_particles; + smoothing_kernel=WendlandC6Kernel{2}(), + smoothing_length=0.2, + reference_density=1000.0, + time_step=0.01, + correction=GradientCorrection()) + + coefficients = ones(TrixiParticles.nparticles(setup.system)) + coefficients[1] = 0.0 + coefficients[2] = NaN + TrixiParticles.sanitize_kernel_correction_coefficient!(coefficients, setup.system, + setup.semi) + @test coefficients[1:2] == ones(2) + + n = 5 + particle_spacing = 1.0 / n + smoothing_kernel = WendlandC6Kernel{2}() + particles = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); + density=1000.0) + state_equation = StateEquationCole(; sound_speed=10.0, + reference_density=1000.0, exponent=1) + boundary_model = BoundaryModelDummyParticles(particles.density, particles.mass, + SummationDensity(), smoothing_kernel, + 2particle_spacing; state_equation, + correction=ShepardKernelCorrection()) + boundary = WallBoundarySystem(particles, boundary_model) + semi = Semidiscretization(boundary; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + boundary = first(ode.p.semi.systems) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, ode.p.semi, 0.0) + @test boundary.boundary_model.pressure ≈ + state_equation.(boundary.boundary_model.cache.density) + + combined_boundary_model = BoundaryModelDummyParticles(particles.density, + particles.mass, + SummationDensity(), + smoothing_kernel, + 2particle_spacing; + state_equation, + density_correction=ShepardKernelCorrection(), + gradient_correction=MixedKernelGradientCorrection()) + combined_boundary = WallBoundarySystem(particles, combined_boundary_model) + combined_semi = Semidiscretization(combined_boundary; + parallelization_backend=SerialBackend()) + combined_ode = semidiscretize(combined_semi, (0.0, 1.0); reset_threads=false) + combined_v = Array(combined_ode.u0.x[1]) + combined_u = Array(combined_ode.u0.x[2]) + combined_boundary = first(combined_ode.p.semi.systems) + TrixiParticles.update_systems_and_nhs(combined_v, combined_u, combined_ode.p.semi, + 0.0) + @test combined_boundary.boundary_model.correction isa CorrectionConfiguration + @test all(isfinite, + combined_boundary.boundary_model.cache.kernel_correction_coefficient) + @test all(isfinite, combined_boundary.boundary_model.cache.dw_gamma) + @test all(isfinite, combined_boundary.boundary_model.cache.correction_matrix) + + density32 = fill(1000.0f0, 4) + mass32 = fill(10.0f0, 4) + boundary32 = BoundaryModelDummyParticles(density32, mass32, SummationDensity(), + WendlandC6Kernel{2}(), 0.2f0; + state_equation, + correction=MixedKernelGradientCorrection()) + @test eltype(boundary32.cache.dw_gamma) == Float32 + @test eltype(boundary32.cache.correction_matrix) == Float32 + end + + @testset "Discrete moments and polynomial reproduction" begin + identity_matrix = Matrix{Float64}(I, 2, 2) + linear_field(pos) = 2.0 + 3.0 * pos[1] - 2.0 * pos[2] + exact_gradient = [3.0, -2.0] + + for perturbation in (false, true) + raw_setup = update_correction!(correction_setup(nothing; perturbation)) + raw_moments = correction_moments(raw_setup; field=linear_field) + + kernel_setup = update_correction!(correction_setup(KernelCorrection(); + perturbation)) + kernel_moments = correction_moments(kernel_setup; field=linear_field) + @test maximum(abs, kernel_moments.zeroth_gradient_moment) < 2e-12 + + gradient_setup = update_correction!(correction_setup(GradientCorrection(); + perturbation)) + gradient_moments = correction_moments(gradient_setup; field=linear_field) + @test maximum(particle -> norm(gradient_moments.first_gradient_moment[:, :, + particle] - + identity_matrix), + TrixiParticles.eachparticle(gradient_setup.system)) < 2e-12 + @test maximum(particle -> norm(gradient_moments.difference_gradient[:, + particle] - + exact_gradient), + TrixiParticles.eachparticle(gradient_setup.system)) < 5e-12 + + mixed_setup = update_correction!(correction_setup(MixedKernelGradientCorrection(); + perturbation)) + mixed_moments = correction_moments(mixed_setup; field=linear_field) + @test maximum(abs, mixed_moments.zeroth_gradient_moment) < 3e-12 + @test maximum(particle -> norm(mixed_moments.first_gradient_moment[:, :, + particle] - + identity_matrix), + TrixiParticles.eachparticle(mixed_setup.system)) < 3e-12 + @test maximum(particle -> norm(mixed_moments.direct_gradient[:, particle] - + exact_gradient), + TrixiParticles.eachparticle(mixed_setup.system)) < 1e-11 + + blending_factor = 0.4 + blended_setup = update_correction!(correction_setup(BlendedGradientCorrection(blending_factor); + perturbation)) + blended_moments = correction_moments(blended_setup; field=linear_field) + expected = (1 - blending_factor) * raw_moments.first_gradient_moment + for particle in TrixiParticles.eachparticle(blended_setup.system) + expected[:, :, particle] .+= blending_factor * identity_matrix + end + @test maximum(abs, blended_moments.first_gradient_moment - expected) < 2e-12 + + corner = corner_particle(raw_setup.system) + @test norm(raw_moments.first_gradient_moment[:, :, corner] - identity_matrix) > + 1e-2 + end + + combined_setup = update_correction!(correction_setup(; + density_calculator=SummationDensity(), + density_correction=ShepardKernelCorrection(), + gradient_correction=MixedKernelGradientCorrection())) + combined_moments = correction_moments(combined_setup; field=linear_field) + @test maximum(abs, combined_moments.zeroth_gradient_moment) < 3e-12 + @test maximum(particle -> norm(combined_moments.first_gradient_moment[:, :, + particle] - + identity_matrix), + TrixiParticles.eachparticle(combined_setup.system)) < 3e-12 + + particle_spacing = 0.25 + smoothing_kernel = WendlandC6Kernel{3}() + particles = RectangularShape(particle_spacing, (4, 4, 4), (0.0, 0.0, 0.0); + density=1000.0) + state_equation = StateEquationCole(; sound_speed=10.0, + reference_density=1000.0, exponent=1) + system = WeaklyCompressibleSPHSystem(particles; smoothing_kernel, + smoothing_length=2particle_spacing, + density_calculator=ContinuityDensity(), + state_equation, + correction=GradientCorrection()) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + system = first(ode.p.semi.systems) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, ode.p.semi, 0.0) + v = TrixiParticles.wrap_v(v_ode, system, ode.p.semi) + u = TrixiParticles.wrap_u(u_ode, system, ode.p.semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + first_moment = zeros(3, 3) + GC.@preserve v_ode u_ode begin + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + ode.p.semi; + points=1:1) do particle, + neighbor, + pos_diff, + distance + volume = TrixiParticles.hydrodynamic_mass(system, neighbor) / + TrixiParticles.current_density(v, system, neighbor) + gradient = TrixiParticles.smoothing_kernel_grad(system, SVector(pos_diff), + distance, particle) + for j in 1:3, i in 1:3 + first_moment[i, j] -= volume * gradient[i] * pos_diff[j] + end + end + end + @test first_moment ≈ Matrix{Float64}(I, 3, 3) atol = 3e-12 + + collinear_coordinates = [0.0 0.1 0.2; 0.0 0.0 0.0] + collinear = InitialCondition(; coordinates=collinear_coordinates, + velocity=zeros(2, 3), density=fill(1000.0, 3), + particle_spacing=0.1) + collinear_system = WeaklyCompressibleSPHSystem(collinear; + smoothing_kernel=WendlandC6Kernel{2}(), + smoothing_length=0.2, + density_calculator=ContinuityDensity(), + state_equation, + correction=GradientCorrection()) + collinear_semi = Semidiscretization(collinear_system; + parallelization_backend=SerialBackend()) + collinear_ode = semidiscretize(collinear_semi, (0.0, 1.0); + reset_threads=false) + collinear_v = Array(collinear_ode.u0.x[1]) + collinear_u = Array(collinear_ode.u0.x[2]) + collinear_system = first(collinear_ode.p.semi.systems) + TrixiParticles.update_systems_and_nhs(collinear_v, collinear_u, + collinear_ode.p.semi, 0.0) + for particle in TrixiParticles.eachparticle(collinear_system) + @test TrixiParticles.correction_matrix(collinear_system, particle) == I + end + + nearly_collinear_coordinates = [0.0 0.1 0.2; 0.0 1.0e-12 0.0] + nearly_collinear = InitialCondition(; coordinates=nearly_collinear_coordinates, + velocity=zeros(2, 3), + density=fill(1000.0, 3), + particle_spacing=0.1) + nearly_collinear_system = WeaklyCompressibleSPHSystem(nearly_collinear; + smoothing_kernel=WendlandC6Kernel{2}(), + smoothing_length=0.2, + density_calculator=ContinuityDensity(), + state_equation, + correction=GradientCorrection()) + nearly_collinear_semi = Semidiscretization(nearly_collinear_system; + parallelization_backend=SerialBackend()) + nearly_collinear_ode = semidiscretize(nearly_collinear_semi, (0.0, 1.0); + reset_threads=false) + nearly_collinear_v = Array(nearly_collinear_ode.u0.x[1]) + nearly_collinear_u = Array(nearly_collinear_ode.u0.x[2]) + nearly_collinear_system = first(nearly_collinear_ode.p.semi.systems) + TrixiParticles.update_systems_and_nhs(nearly_collinear_v, nearly_collinear_u, + nearly_collinear_ode.p.semi, 0.0) + for particle in TrixiParticles.eachparticle(nearly_collinear_system) + @test TrixiParticles.correction_matrix(nearly_collinear_system, particle) == I + end + end + + @testset "Shepard partition of unity" begin + setup = correction_setup(nothing) + (; system, semi, v_ode, u_ode) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coefficient = zeros(TrixiParticles.nparticles(system)) + numerator = zero(coefficient) + + TrixiParticles.compute_shepard_coeff!(system, + TrixiParticles.current_coordinates(u, system), + v_ode, u_ode, semi, coefficient) + coordinates = TrixiParticles.current_coordinates(u, system) + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + semi) do particle, neighbor, pos_diff, + distance + numerator[particle] += TrixiParticles.hydrodynamic_mass(system, neighbor) * + TrixiParticles.smoothing_kernel(system, distance, + particle) + end + + @test numerator ./ coefficient ≈ fill(1000.0, length(numerator)) atol = 2e-12 + @test TrixiParticles.current_density(v, system) == fill(1000.0, length(numerator)) + end + + @testset "Manufactured continuity equation" begin + analytic_density_rate = -2000.0 + errors = Dict{Any, Float64}() + corrections = (nothing, KernelCorrection(), GradientCorrection(), + BlendedGradientCorrection(0.4), MixedKernelGradientCorrection()) + + for correction in corrections + setup = correction_setup(correction) + dv_ode = zero(setup.v_ode) + TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode, + (; semi=setup.semi, split_integration_data=nothing), 0.0) + dv = TrixiParticles.wrap_v(dv_ode, setup.system, setup.semi) + error = dv[end, :] .- analytic_density_rate + errors[correction] = sqrt(sum(abs2, error) / length(error)) + end + + @test errors[GradientCorrection()] < 2e-10 + @test errors[MixedKernelGradientCorrection()] < 2e-10 + @test errors[BlendedGradientCorrection(0.4)] < errors[nothing] + @test errors[nothing] > 1.0 + + for correction in (GradientCorrection(), MixedKernelGradientCorrection()) + setup = correction_setup(correction; edac=true) + dv_ode = zero(setup.v_ode) + TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode, + (; semi=setup.semi, split_integration_data=nothing), 0.0) + dv = TrixiParticles.wrap_v(dv_ode, setup.system, setup.semi) + pressure_error = dv[3, :] .+ 200000.0 + density_error = dv[4, :] .+ 2000.0 + @test sqrt(sum(abs2, pressure_error) / length(pressure_error)) < 2e-8 + @test sqrt(sum(abs2, density_error) / length(density_error)) < 2e-10 + end + end + + @testset "Supported pressure variation matrix" begin + function set_pressure_field!(setup, edac) + pressure = range(1.0, 2.0; + length=TrixiParticles.nparticles(setup.system)) + if edac + v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi) + v[3, :] .= pressure + elseif setup.system.density_calculator isa ContinuityDensity + v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi) + v[end, :] .= 1000.0 .+ pressure + else + setup.system.pressure .= pressure + end + return setup + end + + summation_corrections = ((; correction=nothing, density_correction=nothing, + gradient_correction=nothing), + (; correction=ShepardKernelCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=KernelCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=GradientCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=BlendedGradientCorrection(0.5), + density_correction=nothing, + gradient_correction=nothing), + (; correction=MixedKernelGradientCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=nothing, + density_correction=ShepardKernelCorrection(), + gradient_correction=MixedKernelGradientCorrection())) + continuity_corrections = ((; correction=nothing, density_correction=nothing, + gradient_correction=nothing), + (; correction=KernelCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=GradientCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=BlendedGradientCorrection(0.5), + density_correction=nothing, + gradient_correction=nothing), + (; correction=MixedKernelGradientCorrection(), + density_correction=nothing, + gradient_correction=nothing), + (; correction=nothing, density_correction=nothing, + gradient_correction=MixedKernelGradientCorrection())) + summation_pressure = (nothing, + TrixiParticles.pressure_acceleration_summation_density, + TrixiParticles.inter_particle_averaged_pressure) + continuity_pressure = (nothing, + TrixiParticles.pressure_acceleration_continuity_density, + TrixiParticles.inter_particle_averaged_pressure) + + for edac in (false, true), configuration in summation_corrections, + pressure_acceleration in summation_pressure + setup = correction_setup(configuration.correction; n=4, edac, + density_calculator=SummationDensity(), + density_correction=configuration.density_correction, + gradient_correction=configuration.gradient_correction, + pressure_acceleration) + set_pressure_field!(setup, edac) + dv_ode = zero(setup.v_ode) + TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode, + (; semi=setup.semi, split_integration_data=nothing), 0.0) + @test all(isfinite, dv_ode) + @test any(!iszero, view(dv_ode, 1:2, :)) + end + + for edac in (false, true), configuration in continuity_corrections, + pressure_acceleration in continuity_pressure + setup = correction_setup(configuration.correction; n=4, edac, + density_calculator=ContinuityDensity(), + density_correction=configuration.density_correction, + gradient_correction=configuration.gradient_correction, + pressure_acceleration) + set_pressure_field!(setup, edac) + dv_ode = zero(setup.v_ode) + TrixiParticles.kick!(dv_ode, setup.v_ode, setup.u_ode, + (; semi=setup.semi, split_integration_data=nothing), 0.0) + @test all(isfinite, dv_ode) + @test any(!iszero, view(dv_ode, 1:2, :)) + end + + for edac in (false, true) + setup_tensile = correction_setup(; n=4, edac, + density_calculator=ContinuityDensity(), + pressure_acceleration=tensile_instability_control) + set_pressure_field!(setup_tensile, edac) + dv_tensile = zero(setup_tensile.v_ode) + TrixiParticles.kick!(dv_tensile, setup_tensile.v_ode, setup_tensile.u_ode, + (; semi=setup_tensile.semi, + split_integration_data=nothing), 0.0) + @test all(isfinite, dv_tensile) + @test any(!iszero, view(dv_tensile, 1:2, :)) + + for configuration in continuity_corrections[2:end] + @test_throws ArgumentError correction_setup(configuration.correction; + n=4, edac, + density_calculator=ContinuityDensity(), + density_correction=configuration.density_correction, + gradient_correction=configuration.gradient_correction, + pressure_acceleration=tensile_instability_control) + end + end + end + + @testset "Continuity density reinitialization" begin + setup = correction_setup() + v = TrixiParticles.wrap_v(setup.v_ode, setup.system, setup.semi) + u = TrixiParticles.wrap_u(setup.u_ode, setup.system, setup.semi) + TrixiParticles.reinit_density!(setup.system, v, u, setup.v_ode, setup.u_ode, + setup.semi) + + @test TrixiParticles.current_density(v, setup.system) ≈ fill(1000.0, 81) atol = 2e-12 + @test maximum(abs, setup.system.pressure) < 2e-10 + end + + @testset "Analytical operator scaling" begin + include(joinpath(validation_dir(), "corrections", "convergence.jl")) + results = CorrectionConvergence.run_convergence(; resolutions=(12, 24, 48)) + @test all(result -> isfinite(result.error), results) + + function finest(method, operator, region) + return last(filter(result -> result.method == method && + result.operator == operator && + result.region == region, + results)) + end + + raw_interpolation_boundary = finest(:none, :interpolation, :boundary) + shepard_interpolation_boundary = finest(:shepard, :interpolation, :boundary) + raw_difference_boundary = finest(:none, :difference_gradient, :boundary) + gradient_difference_boundary = finest(:gradient, :difference_gradient, :boundary) + blended_difference_boundary = finest(:blended, :difference_gradient, :boundary) + mixed_difference_boundary = finest(:mixed, :difference_gradient, :boundary) + raw_direct_boundary = finest(:none, :direct_gradient, :boundary) + kernel_direct_boundary = finest(:kernel, :direct_gradient, :boundary) + mixed_direct_boundary = finest(:mixed, :direct_gradient, :boundary) + + @test shepard_interpolation_boundary.order > + raw_interpolation_boundary.order + 0.9 + @test gradient_difference_boundary.order > raw_difference_boundary.order + 0.9 + @test mixed_difference_boundary.order > raw_difference_boundary.order + 0.9 + @test kernel_direct_boundary.order > raw_direct_boundary.order + 0.9 + @test mixed_direct_boundary.order > kernel_direct_boundary.order + 0.9 + @test blended_difference_boundary.error < raw_difference_boundary.error + + shepard_interpolation_interior = finest(:shepard, :interpolation, :interior) + gradient_difference_interior = finest(:gradient, :difference_gradient, :interior) + mixed_difference_interior = finest(:mixed, :difference_gradient, :interior) + mixed_direct_interior = finest(:mixed, :direct_gradient, :interior) + raw_density_boundary = finest(:none, :summation_density, :boundary) + shepard_density_boundary = finest(:shepard, :summation_density, :boundary) + reinitialized_density_boundary = finest(:shepard, :density_reinitialization, + :boundary) + reinitialized_density_interior = finest(:shepard, :density_reinitialization, + :interior) + + @test shepard_interpolation_interior.order > 1.8 + @test gradient_difference_interior.order > 1.8 + @test mixed_difference_interior.order > 1.8 + @test mixed_direct_interior.order > 1.8 + @test shepard_density_boundary.error < raw_density_boundary.error + @test abs(shepard_density_boundary.order) < 0.1 + @test reinitialized_density_boundary.order > 0.9 + @test reinitialized_density_interior.order > 1.8 + + pressure_operators = (:pressure_summation, + :pressure_interparticle_summation, + :pressure_continuity, + :pressure_interparticle_continuity) + for operator in pressure_operators + @test finest(:gradient, operator, :interior).order > 1.8 + @test finest(:mixed, operator, :interior).order > 1.8 + end + for operator in (:pressure_summation, :pressure_interparticle_summation) + @test finest(:shepard_mixed, operator, :interior).order > 1.8 + end + + constant_pressure_results = filter(results) do result + startswith(string(result.operator), "constant_pressure_") && + result.region == :interior && result.resolution == 48 + end + @test !isempty(constant_pressure_results) + @test maximum(result -> result.error, constant_pressure_results) < 1e-7 + + for region in (:boundary, :interior) + tensile = finest(:none, :pressure_tensile_positive, region) + continuity = finest(:none, :pressure_continuity, region) + @test tensile.error ≈ continuity.error rtol = 5e-13 + end + end +end diff --git a/test/general/general.jl b/test/general/general.jl index acb07de8b9..2ab768000e 100644 --- a/test/general/general.jl +++ b/test/general/general.jl @@ -1,6 +1,7 @@ include("initial_condition.jl") include("smoothing_kernels.jl") include("density_calculator.jl") +include("corrections.jl") include("semidiscretization.jl") include("interpolation.jl") include("buffer.jl") diff --git a/test/general/interpolation.jl b/test/general/interpolation.jl index 5469876724..bd8420ecd4 100644 --- a/test/general/interpolation.jl +++ b/test/general/interpolation.jl @@ -131,6 +131,111 @@ semi_boundary = Semidiscretization(fluid_system, boundary_system) TrixiParticles.initialize_neighborhood_searches!(semi_boundary) + surface_detection_ic = RectangularShape(particle_spacing, (nx, ny), (0.0, 0.0), + density=1000.0) + surface_detection_system = WeaklyCompressibleSPHSystem(surface_detection_ic; + smoothing_kernel, + smoothing_length=1.5 * + particle_spacing, + density_calculator=ContinuityDensity(), + state_equation, viscosity, + acceleration=(0.0, -9.81), + surface_method=ColorfieldSurfaceDetection(ideal_density_threshold=0.9), + reference_particle_spacing=particle_spacing) + surface_detection_system.pressure .= surface_detection_ic.pressure + semi_surface_detection = Semidiscretization(surface_detection_system) + TrixiParticles.initialize_neighborhood_searches!(semi_surface_detection) + + @testset verbose=true "Interpolated Free Surface Detection" begin + min_x, max_x = extrema(view(surface_detection_ic.coordinates, 1, :)) + max_y = maximum(view(surface_detection_ic.coordinates, 2, :)) + point_coords = [(min_x + max_x) / 2 (min_x + max_x) / 2; + max_y max_y + particle_spacing] + v_surface_detection = vcat(surface_detection_ic.velocity, + surface_detection_ic.density') + u_surface_detection = surface_detection_ic.coordinates + + uncut = interpolate_points(point_coords, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; + cut_off_bnd=false) + cut = interpolate_points(point_coords, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; + cut_off_bnd=true) + + @test all(isfinite, uncut.density) + @test all(isfinite, uncut.surface_activity) + @test uncut.surface_activity[2] > 0.9 + @test isfinite(cut.density[1]) + @test cut.neighbor_count[1] > 0 + @test isnan(cut.density[2]) + @test isnan(cut.surface_activity[2]) + @test cut.neighbor_count[2] == 0 + + mktempdir() do output_directory + interpolate_plane_2d_vtk([min_x, max_y - particle_spacing], + [max_x, max_y + particle_spacing], + particle_spacing, semi_surface_detection, + surface_detection_system, v_surface_detection, + u_surface_detection; output_directory, + filename="surface_detection") + vtk_file = TrixiParticles.ReadVTK.VTKFile(joinpath(output_directory, + "surface_detection.vti")) + point_data = TrixiParticles.ReadVTK.get_point_data(vtk_file) + @test "surface_activity" in keys(point_data) + end + end + + @testset verbose=true "Interpolated Multiphase Surface Detection" begin + interface_spacing = 0.1 + y_coordinates = collect(-0.5:interface_spacing:0.5) + coordinates_a = hcat(([x, y] for x in -0.5:interface_spacing:-0.1 + for y in y_coordinates)...) + coordinates_b = hcat(([x, y] for x in 0.0:interface_spacing:0.5 + for y in y_coordinates)...) + initial_condition_a = InitialCondition(; coordinates=coordinates_a, + density=fill(1000.0, + size(coordinates_a, 2)), + particle_spacing=interface_spacing) + initial_condition_b = InitialCondition(; coordinates=coordinates_b, + density=fill(1000.0, + size(coordinates_b, 2)), + particle_spacing=interface_spacing) + interface_kernel = WendlandC2Kernel{2}() + interface_state_equation = StateEquationCole(sound_speed=10.0, + reference_density=1000.0, + exponent=1) + system_a = WeaklyCompressibleSPHSystem(initial_condition_a; + smoothing_kernel=interface_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation=interface_state_equation, + surface_method=ColorfieldSurfaceDetection(interface_threshold=1.0e-6), + reference_particle_spacing=interface_spacing, + color_value=1) + system_b = WeaklyCompressibleSPHSystem(initial_condition_b; + smoothing_kernel=interface_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation=interface_state_equation, + color_value=2) + interface_semi = Semidiscretization(system_a, system_b) + interface_ode = semidiscretize(interface_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(interface_ode.u0.x..., interface_semi, + 0.0) + interface_points = [-0.1 0.1; 0.0 0.0] + uncut = interpolate_points(interface_points, interface_semi, system_a, + interface_ode.u0.x...; cut_off_bnd=false) + cut = interpolate_points(interface_points, interface_semi, system_a, + interface_ode.u0.x...; cut_off_bnd=true) + + @test all(>(0.9), uncut.surface_activity) + @test isfinite(cut.density[1]) + @test isnan(cut.density[2]) + @test isnan(cut.surface_activity[2]) + end + # Some simple results expected_zero(y) = (density=[NaN], neighbor_count=[0], point_coords=[0.0; y;;], velocity=[NaN; NaN;;], pressure=[NaN]) diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 5b92d133ba..d9ee26f079 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -132,7 +132,7 @@ # Mock fluid system struct FluidSystemMock <: TrixiParticles.AbstractFluidSystem{2} surface_tension::Nothing - surface_normal_method::Nothing + surface_method::Nothing FluidSystemMock() = new(nothing, nothing) end @@ -342,7 +342,7 @@ u = TrixiParticles.wrap_u(u_ode, system, semi) TrixiParticles.compute_correction_values!(system, - TrixiParticles.system_correction(system), + TrixiParticles.correction_density(system.correction), u, v_ode, u_ode, semi) return copy(system.cache.kernel_correction_coefficient), semi diff --git a/test/schemes/fluid/pressure_acceleration.jl b/test/schemes/fluid/pressure_acceleration.jl index 80de3d9a66..fca175cde0 100644 --- a/test/schemes/fluid/pressure_acceleration.jl +++ b/test/schemes/fluid/pressure_acceleration.jl @@ -12,6 +12,183 @@ @test f_2 == TrixiParticles.pressure_acceleration_continuity_density end + @testset "Algebraic formulations and asymmetric conservation" begin + m_a, m_b = 1.2, 0.8 + rho_a, rho_b = 1000.0, 980.0 + p_a, p_b = 2.0, 3.0 + W_a = SVector(0.2, -0.1) + W_b = -W_a + + summation = TrixiParticles.pressure_acceleration_summation_density + continuity = TrixiParticles.pressure_acceleration_continuity_density + interparticle = TrixiParticles.inter_particle_averaged_pressure + difference = TrixiParticles.pressure_acceleration_difference + + @test summation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈ + -m_b * (p_a / rho_a^2 + p_b / rho_b^2) * W_a + @test continuity(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈ + -m_b * (p_a + p_b) / (rho_a * rho_b) * W_a + + volume_term = ((m_a / rho_a)^2 + (m_b / rho_b)^2) / m_a + pressure_tilde = (rho_b * p_a + rho_a * p_b) / (rho_a + rho_b) + @test interparticle(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) ≈ + -volume_term * pressure_tilde * W_a + @test tensile_instability_control(m_a, m_b, rho_a, rho_b, -p_a, p_b, W_a) ≈ + -m_b * (p_a + p_b) / (rho_a * rho_b) * W_a + + for pressure_formulation in (summation, continuity, interparticle) + symmetric = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) + asymmetric = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, + W_a, W_b) + @test asymmetric ≈ symmetric + @test pressure_formulation(m_a, m_b, rho_a, rho_b, 0.0, 0.0, + W_a, W_b) == zero(W_a) + + acceleration_a = pressure_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, + W_a, W_b) + acceleration_b = pressure_formulation(m_b, m_a, rho_b, rho_a, p_b, p_a, + W_b, W_a) + @test m_a * acceleration_a + m_b * acceleration_b ≈ zero(W_a) atol = eps() + end + @test tensile_instability_control(m_a, m_b, rho_a, rho_b, 0.0, 0.0, + W_a) == zero(W_a) + + difference_acceleration = difference(m_b, rho_a, rho_b, p_a, p_b, W_a) + conservative_acceleration = continuity(m_a, m_b, rho_a, rho_b, p_a, p_b, + W_a, W_b) + model = SurfacePressureDifference() + @test difference_acceleration ≈ + -m_b * (p_b - p_a) / (rho_a * rho_b) * W_a + @test difference(m_b, rho_a, rho_b, p_a, p_a, W_a) == zero(W_a) + @test TrixiParticles.blend_surface_pressure(model, conservative_acceleration, + 0.0, true, m_b, rho_a, rho_b, p_a, + p_b, W_a) == conservative_acceleration + @test TrixiParticles.blend_surface_pressure(model, conservative_acceleration, + 1.0, true, m_b, rho_a, rho_b, p_a, + p_b, W_a) == difference_acceleration + @test TrixiParticles.blend_surface_pressure(model, conservative_acceleration, + 0.5, true, m_b, rho_a, rho_b, p_a, + p_b, W_a) ≈ + 0.5 * (conservative_acceleration + difference_acceleration) + @test TrixiParticles.blend_surface_pressure(model, conservative_acceleration, + 1.0, false, m_b, rho_a, rho_b, p_a, + p_b, W_a) == conservative_acceleration + end + + @testset "Surface pressure configuration" begin + particle_spacing = 0.1 + initial_condition = RectangularShape(particle_spacing, (4, 4), (0.0, 0.0), + density=1000.0) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=1000.0, + exponent=1) + surface_method_ = ColorfieldSurfaceDetection(ideal_density_threshold=0.9) + model = SurfacePressureDifference() + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation, + gradient_correction=GradientCorrection(), + surface_method=surface_method_, + surface_pressure=:invalid, + reference_particle_spacing=particle_spacing) + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation, + gradient_correction=GradientCorrection(), + surface_pressure=model) + @test_throws ArgumentError WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation, + gradient_correction=KernelCorrection(), + surface_method=surface_method_, + surface_pressure=model, + reference_particle_spacing=particle_spacing) + + for correction in (GradientCorrection(), MixedKernelGradientCorrection()) + wcsph = WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation, + gradient_correction=correction, + surface_method=surface_method_, + surface_pressure=model, + reference_particle_spacing=particle_spacing) + edac = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length=0.15, sound_speed=10.0, + density_calculator=SummationDensity(), + gradient_correction=correction, + surface_method=surface_method_, + surface_pressure=model, + reference_particle_spacing=particle_spacing) + @test wcsph.surface_pressure === model + @test edac.surface_pressure === model + end + + system = WeaklyCompressibleSPHSystem(initial_condition; + smoothing_kernel, + smoothing_length=0.15, + density_calculator=SummationDensity(), + state_equation, + gradient_correction=GradientCorrection(), + surface_method=surface_method_, + surface_pressure=model, + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 0.01); reset_threads=false) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = TrixiParticles.current_coordinates(u, system) + particle, neighbor = 1, 2 + pos_diff = SVector{2}(coordinates[:, particle] - coordinates[:, neighbor]) + distance = norm(pos_diff) + W_a = TrixiParticles.smoothing_kernel_grad(system, pos_diff, distance, particle) + m_a = TrixiParticles.hydrodynamic_mass(system, particle) + m_b = TrixiParticles.hydrodynamic_mass(system, neighbor) + rho_a = TrixiParticles.current_density(v, system, particle) + rho_b = TrixiParticles.current_density(v, system, neighbor) + p_a, p_b = 2.0, 3.0 + gradient_correction = TrixiParticles.correction_gradient(system.correction) + + system.cache.surface_activity[particle] = 1 + acceleration = TrixiParticles.pressure_acceleration(system, system, particle, + neighbor, m_a, m_b, p_a, p_b, + rho_a, rho_b, pos_diff, + distance, + W_a, gradient_correction) + @test acceleration ≈ TrixiParticles.pressure_acceleration_difference(m_b, rho_a, + rho_b, p_a, + p_b, W_a) + W_b = TrixiParticles.smoothing_kernel_grad(system, -pos_diff, distance, neighbor) + conservative_acceleration = TrixiParticles.pressure_acceleration(system, system, + particle, neighbor, + m_a, m_b, p_a, p_b, + rho_a, rho_b, + pos_diff, distance, + W_a, + gradient_correction; + use_surface_pressure=false) + @test conservative_acceleration ≈ + system.pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, + W_a, W_b) + constant_acceleration = TrixiParticles.pressure_acceleration(system, system, + particle, neighbor, + m_a, m_b, p_a, p_a, + rho_a, rho_b, pos_diff, + distance, W_a, + gradient_correction) + @test constant_acceleration == zero(W_a) + end + @testset verbose=true "Illegal Inputs" begin correction_dict_1 = Dict( "KernelCorrection" => KernelCorrection(), diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 5eb8a81704..cb0e770051 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -61,7 +61,7 @@ end function create_fluid_system(coordinates, velocity, mass, density, particle_spacing, surface_tension; - surface_normal_method=ColorfieldSurfaceNormal(), + surface_method=ColorfieldSurfaceNormal(), color_value=1, NDIMS=2, smoothing_length=1.0, wall=false, walldistance=0.0, boundary_system_type=:wall, smoothing_kernel=SchoenbergCubicSplineKernel{NDIMS}()) @@ -76,9 +76,9 @@ function create_fluid_system(coordinates, velocity, mass, density, particle_spac system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, density_calculator=SummationDensity(), state_equation, - surface_normal_method, + surface_method, reference_particle_spacing=particle_spacing, - surface_tension) + surface_tension, color_value) if wall boundary_system = if boundary_system_type == :wall @@ -109,12 +109,11 @@ function compute_and_test_surface_values(system, semi, ode; NDIMS=2) v = TrixiParticles.wrap_v(v0_ode, system, semi) u = TrixiParticles.wrap_u(u0_ode, system, semi) - # Compute the surface normals - TrixiParticles.compute_surface_normal!(system, system.surface_normal_method, v, u, - v0_ode, u0_ode, semi, 0.0) + TrixiParticles.compute_surface!(system, system.surface_method, v, u, + v0_ode, u0_ode, semi, 0.0) TrixiParticles.remove_invalid_normals!(system, system.surface_tension, - system.surface_normal_method) + system.surface_method) # After computation, check that surface normals have been computed and are not NaN or Inf @test all(isfinite, system.cache.surface_normal) @@ -143,6 +142,225 @@ function compute_curvature!(system, semi, ode) v, u, v0_ode, u0_ode, semi, 0.0) end +@testset verbose=true "Colorfield Surface Detection" begin + normal_method = ColorfieldSurfaceNormal(ideal_density_threshold=0.9) + detection_method = ColorfieldSurfaceDetection(ideal_density_threshold=0.9) + @test ColorfieldSurfaceNormal() == ColorfieldSurfaceNormal(0.1, 0.01, 0.0) + @test normal_method.interpolation_surface_threshold == 0.45 + @test detection_method.interpolation_surface_threshold == 0.45 + @test TrixiParticles.computes_surface_normal(normal_method) + @test !TrixiParticles.computes_surface_normal(detection_method) + + @test_throws ArgumentError ColorfieldSurfaceNormal(boundary_contact_threshold=-0.1) + @test_throws ArgumentError ColorfieldSurfaceNormal(interface_threshold=Inf) + @test_throws ArgumentError ColorfieldSurfaceNormal(interface_taper_start=1.0) + @test_throws ArgumentError ColorfieldSurfaceNormal(interpolation_surface_threshold=1.1) + @test_throws ArgumentError ColorfieldSurfaceDetection(interface_threshold=-0.1) + @test_throws ArgumentError ColorfieldSurfaceDetection(interface_threshold="invalid") + + particle_spacing = 0.1 + coordinates = RectangularShape(particle_spacing, (21, 11), (0.0, 0.0), + density=1000.0) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + system, _, semi, + ode = create_fluid_system(coordinates.coordinates, coordinates.velocity, + coordinates.mass, coordinates.density, particle_spacing, + nothing; smoothing_length=1.5 * particle_spacing, + smoothing_kernel, surface_method=normal_method) + + activity = system.cache.surface_activity + x = coordinates.coordinates + min_x, max_x = extrema(view(x, 1, :)) + min_y, max_y = extrema(view(x, 2, :)) + interior = [particle + for particle in eachindex(activity) + if min_x + 3particle_spacing < x[1, particle] < + max_x - 3particle_spacing && + min_y + 3particle_spacing < x[2, particle] < + max_y - 3particle_spacing] + top_surface = [particle + for particle in eachindex(activity) + if x[2, particle] == max_y && + min_x + 3particle_spacing < x[1, particle] < + max_x - 3particle_spacing] + + @test all(isfinite, activity) + @test all(iszero, activity[interior]) + @test all(==(1), activity[top_surface]) + + detection_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + detection_semi = Semidiscretization(detection_system) + detection_ode = semidiscretize(detection_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(detection_ode.u0.x..., detection_semi, 0.0) + @test !haskey(detection_system.cache, :surface_normal) + @test detection_system.cache.surface_activity == activity + @test isapprox(detection_system.cache.surface_gradient, system.cache.surface_normal; + rtol=10eps(), atol=10eps()) + + edac_system = EntropicallyDampedSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + sound_speed=10.0, + density_calculator=SummationDensity(), + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + edac_semi = Semidiscretization(edac_system) + edac_ode = semidiscretize(edac_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(edac_ode.u0.x..., edac_semi, 0.0) + @test edac_system.cache.surface_activity == activity + + iisph_system = ImplicitIncompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + reference_density=1000.0, + time_step=0.001, + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + iisph_semi = Semidiscretization(iisph_system) + iisph_ode = semidiscretize(iisph_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(iisph_ode.u0.x..., iisph_semi, 0.0) + @test iisph_system.cache.surface_activity == activity + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * + particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_tension=SurfaceTensionMorris(), + surface_method=detection_method, + reference_particle_spacing=particle_spacing) + cohesion_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + surface_tension=CohesionForceAkinci()) + @test isnothing(cohesion_system.surface_method) + @test !haskey(cohesion_system.cache, :surface_activity) + + corrected_system = WeaklyCompressibleSPHSystem(coordinates; + smoothing_kernel, + smoothing_length=1.5 * particle_spacing, + density_calculator=SummationDensity(), + state_equation=system.state_equation, + correction=GradientCorrection(), + surface_method=normal_method, + reference_particle_spacing=particle_spacing) + corrected_semi = Semidiscretization(corrected_system) + corrected_ode = semidiscretize(corrected_semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(corrected_ode.u0.x..., corrected_semi, 0.0) + @test corrected_system.cache.surface_activity == activity + @test isapprox(corrected_system.cache.surface_normal, system.cache.surface_normal; + rtol=10eps(), atol=10eps()) + + v_ode, u_ode = ode.u0.x + @test TrixiParticles.surface_activity(system, nothing, nothing, v_ode, u_ode, + semi, 0.0) == activity + @test TrixiParticles.surface_normal(detection_system, nothing, nothing, + detection_ode.u0.x..., detection_semi, 0.0) === + nothing + + metadata = Dict{String, Any}() + TrixiParticles.add_system_data!(metadata, normal_method) + @test metadata["surface_method"]["computes_surface_normal"] + @test metadata["surface_method"]["interpolation_surface_threshold"] == 0.45 + + mktempdir() do output_directory + trixi2vtk(ode.u0, semi, 0.0; output_directory, + prefix="surface_detection", overwrite=true) + vtk_data = vtk2trixi(joinpath(output_directory, + "surface_detection_fluid_1_current.vtu")) + @test vtk_data.surface_activity == activity + @test hasproperty(vtk_data, :surf_normal) + + trixi2vtk(detection_ode.u0, detection_semi, 0.0; output_directory, + prefix="detection_only", overwrite=true) + detection_vtk_data = vtk2trixi(joinpath(output_directory, + "detection_only_fluid_1_current.vtu")) + @test detection_vtk_data.surface_activity == activity + @test !hasproperty(detection_vtk_data, :surf_normal) + end +end + +@testset verbose=true "Multicolor Surface Activity" begin + particle_spacing = 0.1 + smoothing_length = 0.15 + y_coordinates = collect(-0.5:particle_spacing:0.5) + coordinates_a = hcat(([x, y] for x in -0.5:particle_spacing:-0.1 + for y in y_coordinates)...) + coordinates_b = hcat(([x, y] for x in 0.0:particle_spacing:0.5 + for y in y_coordinates)...) + smoothing_kernel = WendlandC2Kernel{2}() + state_equation = StateEquationCole(sound_speed=10.0, reference_density=1000.0, + exponent=1) + + function interface_geometry(color_a, color_b, surface_method_) + initial_condition_a = InitialCondition(; coordinates=coordinates_a, + density=fill(1000.0, + size(coordinates_a, 2)), + particle_spacing) + initial_condition_b = InitialCondition(; coordinates=coordinates_b, + density=fill(1000.0, + size(coordinates_b, 2)), + particle_spacing) + system_a = WeaklyCompressibleSPHSystem(initial_condition_a; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, + surface_method=surface_method_, + reference_particle_spacing=particle_spacing, + color_value=color_a) + system_b = WeaklyCompressibleSPHSystem(initial_condition_b; smoothing_kernel, + smoothing_length, + density_calculator=SummationDensity(), + state_equation, color_value=color_b) + semi = Semidiscretization(system_a, system_b) + ode = semidiscretize(semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + interface_particle = argmin(eachindex(eachcol(coordinates_a))) do particle + abs(coordinates_a[1, particle] + particle_spacing) + + abs(coordinates_a[2, particle]) + end + gradient = if surface_method_ isa ColorfieldSurfaceNormal + TrixiParticles.surface_normal(system_a, interface_particle) + else + TrixiParticles.extract_svector(system_a.cache.surface_gradient, system_a, + interface_particle) + end + return gradient, TrixiParticles.surface_activity(system_a, interface_particle), + system_b + end + + detection_method = ColorfieldSurfaceDetection(interface_threshold=1.0e-6) + normal_method = ColorfieldSurfaceNormal(interface_threshold=1.0e-6) + increasing_gradient, increasing_activity, + non_surface_neighbor = interface_geometry(0, 2, detection_method) + unit_gradient, _, _ = interface_geometry(0, 1, detection_method) + decreasing_gradient, decreasing_activity, _ = interface_geometry(2, 0, + detection_method) + equal_gradient, equal_activity, _ = interface_geometry(1, 1, detection_method) + normal_gradient, normal_activity, _ = interface_geometry(0, 2, normal_method) + + @test isnothing(non_surface_neighbor.surface_method) + @test increasing_gradient[1] > 0 + @test decreasing_gradient[1] < 0 + @test isapprox(norm(increasing_gradient), 2norm(unit_gradient); rtol=1.0e-12) + @test norm(equal_gradient) < 100eps() + @test increasing_activity == 1 + @test decreasing_activity == 1 + @test equal_activity == 0 + @test normal_activity == increasing_activity + @test normal_gradient == increasing_gradient +end + @testset verbose=true "Rigid Dummy Boundary Matches Wall Boundary" begin NDIMS = 2 particle_spacing = 0.2 @@ -161,16 +379,16 @@ end wall_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0, boundary_system_type=:wall) rigid_system, rigid_boundary, rigid_semi, rigid_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0, boundary_system_type=:rigid) @@ -186,6 +404,9 @@ end @test isapprox(rigid_system.cache.neighbor_count, wall_system.cache.neighbor_count, rtol=sqrt(eps()), atol=sqrt(eps())) + @test isapprox(rigid_system.cache.surface_activity, + wall_system.cache.surface_activity, + rtol=sqrt(eps()), atol=sqrt(eps())) end @testset verbose=true "CSS/CSF: Sphere Surface Normals" begin @@ -219,8 +440,8 @@ end particle_spacing, SurfaceTensionMorris(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0) compute_and_test_surface_values(system, semi, ode; NDIMS) @@ -319,8 +540,8 @@ end particle_spacing, SurfaceTensionAkinci(surface_tension_coefficient=0.072); NDIMS, smoothing_length, smoothing_kernel, - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), wall=true, walldistance=2.0) compute_and_test_surface_values(system, semi, ode; NDIMS) diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 7fe8abbd97..404d934552 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -115,8 +115,8 @@ density_calculator=density_calc, state_equation=eq_state, surface_tension=SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0), - surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, - ideal_density_threshold=0.9), + surface_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9), reference_particle_spacing=1.0,) # 4. Verify Cache Contains Necessary Fields diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index ce04cd6774..a17061da7d 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -127,7 +127,7 @@ system = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smoothing_length, sound_speed) - show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, nothing, Val{:smoothing_kernel}(), [0.0, 0.0], nothing, nothing) with 2 particles" + show_compact = "EntropicallyDampedSPHSystem{2}(SummationDensity(), nothing, nothing, Val{:smoothing_kernel}(), [0.0, 0.0], nothing, nothing, nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -143,7 +143,8 @@ │ average pressure reduction: ……… no │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ surface tension: …………………………………… nothing │ - │ surface normal method: …………………… nothing │ + │ surface method: ……………………………………… nothing │ + │ surface pressure: ………………………………… nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box end diff --git a/test/systems/wcsph_system.jl b/test/systems/wcsph_system.jl index 37b94f3f64..cee4cba3b5 100644 --- a/test/systems/wcsph_system.jl +++ b/test/systems/wcsph_system.jl @@ -198,7 +198,7 @@ smoothing_length, density_calculator, state_equation, density_diffusion) - show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, nothing, [0.0, 0.0], nothing) with 2 particles" + show_compact = "WeaklyCompressibleSPHSystem{2}(SummationDensity(), nothing, Val{:state_equation}(), Val{:smoothing_kernel}(), nothing, Val{:density_diffusion}(), nothing, nothing, nothing, nothing, [0.0, 0.0], nothing) with 2 particles" @test repr(system) == show_compact show_box = """ ┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ @@ -213,7 +213,8 @@ │ density diffusion: ……………………………… Val{:density_diffusion}() │ │ shifting technique: …………………………… nothing │ │ surface tension: …………………………………… nothing │ - │ surface normal method: …………………… nothing │ + │ surface method: ……………………………………… nothing │ + │ surface pressure: ………………………………… nothing │ │ acceleration: …………………………………………… [0.0, 0.0] │ │ source terms: …………………………………………… Nothing │ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" diff --git a/validation/corrections/convergence.jl b/validation/corrections/convergence.jl new file mode 100644 index 0000000000..d935246d89 --- /dev/null +++ b/validation/corrections/convergence.jl @@ -0,0 +1,501 @@ +module CorrectionConvergence + +using TrixiParticles +using LinearAlgebra: norm +using Printf: @printf, @sprintf + +export run_convergence, print_report, write_csv + +function field(position) + 2.0 + position[1] + 0.5 * position[2] + position[1]^2 - + 0.25 * position[1] * position[2] + 0.75 * position[2]^2 + + 0.2 * position[1]^3 - 0.1 * position[1]^2 * position[2] + + 0.15 * position[1] * position[2]^2 - 0.05 * position[2]^3 +end + +function field_gradient(position) + return SVector(1.0 + 2.0 * position[1] - 0.25 * position[2] + + 0.6 * position[1]^2 - 0.2 * position[1] * position[2] + + 0.15 * position[2]^2, + 0.5 - 0.25 * position[1] + 1.5 * position[2] - + 0.1 * position[1]^2 + 0.3 * position[1] * position[2] - + 0.15 * position[2]^2) +end + +function density_field(position) + return 1000.0 * (1.0 + 0.1 * position[1] + 0.05 * position[2] + + 0.02 * position[1]^2 - 0.01 * position[1] * position[2] + + 0.015 * position[2]^2) +end + +function free_surface_pressure(position) + return position[1] + 0.5 * position[1] * position[2] + 0.2 * position[1]^3 + + 0.1 * position[1] * position[2]^2 +end + +function free_surface_pressure_gradient(position) + return SVector(1.0 + 0.5 * position[2] + 0.6 * position[1]^2 + + 0.1 * position[2]^2, + 0.5 * position[1] + 0.2 * position[1] * position[2]) +end + +function setup_operator(n, correction; density_calculator=ContinuityDensity()) + particle_spacing = 1.0 / n + smoothing_length = 2.0 * particle_spacing + smoothing_kernel = WendlandC6Kernel{2}() + fluid = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); density=1000.0) + state_equation = StateEquationCole(; sound_speed=10.0, reference_density=1000.0, + exponent=1) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator, + state_equation, correction) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + semi = ode.p.semi + system = first(semi.systems) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + + return (; system, semi, v_ode, u_ode, particle_spacing) +end + +function pressure_operator_errors(n, correction, density_calculator, + pressure_formulation) + setup = setup_operator(n, correction; density_calculator) + (; system, semi, v_ode, u_ode, particle_spacing) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + n_particles = TrixiParticles.nparticles(system) + pressure = [free_surface_pressure(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + constant_pressure = fill(2.0, n_particles) + exact_acceleration = zeros(2, n_particles) + for particle in TrixiParticles.eachparticle(system) + exact_acceleration[:, + particle] = -free_surface_pressure_gradient(SVector{2}(view(coordinates, + :, + particle))) / + TrixiParticles.current_density(v, system, particle) + end + + acceleration = zeros(2, n_particles) + constant_acceleration = zeros(2, n_particles) + gradient_correction = TrixiParticles.correction_gradient(system.correction) + asymmetric = gradient_correction isa Union{KernelCorrection, GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection} + GC.@preserve v_ode u_ode begin + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + semi) do particle, neighbor, pos_diff, + distance + m_a = TrixiParticles.hydrodynamic_mass(system, particle) + m_b = TrixiParticles.hydrodynamic_mass(system, neighbor) + rho_a = TrixiParticles.current_density(v, system, particle) + rho_b = TrixiParticles.current_density(v, system, neighbor) + W_a = TrixiParticles.smoothing_kernel_grad(system, SVector(pos_diff), distance, + particle) + + pressure_acceleration = if asymmetric + W_b = TrixiParticles.smoothing_kernel_grad(system, SVector(-pos_diff), + distance, neighbor) + pressure_formulation(m_a, m_b, rho_a, rho_b, pressure[particle], + pressure[neighbor], W_a, W_b) + else + pressure_formulation(m_a, m_b, rho_a, rho_b, pressure[particle], + pressure[neighbor], W_a) + end + constant_pressure_acceleration = if asymmetric + W_b = TrixiParticles.smoothing_kernel_grad(system, SVector(-pos_diff), + distance, neighbor) + pressure_formulation(m_a, m_b, rho_a, rho_b, + constant_pressure[particle], + constant_pressure[neighbor], W_a, W_b) + else + pressure_formulation(m_a, m_b, rho_a, rho_b, + constant_pressure[particle], + constant_pressure[neighbor], W_a) + end + + for dimension in 1:2 + acceleration[dimension, particle] += pressure_acceleration[dimension] + constant_acceleration[dimension, + particle] += constant_pressure_acceleration[dimension] + end + end + end + + support = TrixiParticles.compact_support(system, system) + boundary = [particle + for particle in axes(coordinates, 2) + if isapprox(coordinates[1, particle], minimum(view(coordinates, 1, :)); + atol=eps()) && + 2 * support < coordinates[2, particle] < 1.0 - 2 * support] + # The conservative asymmetric formulation uses the correction matrix of both particles. + # Keep both neighborhoods away from unrelated boundaries in both samples. + interior = [particle + for particle in axes(coordinates, 2) + if 2 * support < coordinates[1, particle] < 1.0 - 2 * support && + 2 * support < coordinates[2, particle] < 1.0 - 2 * support] + isempty(boundary) && error("resolution $n has no particles in the boundary sample") + isempty(interior) && + error("resolution $n has no particles in the pressure interior sample") + + function sample_errors(particles) + manufactured = normalized_l2(acceleration[:, particles], + exact_acceleration[:, particles]) + constant = norm(constant_acceleration[:, particles]) / sqrt(length(particles)) + return (; manufactured, constant) + end + + return (; particle_spacing, boundary=sample_errors(boundary), + interior=sample_errors(interior)) +end + +function setup_summation_density(n, correction) + particle_spacing = 1.0 / n + smoothing_length = 2.0 * particle_spacing + smoothing_kernel = WendlandC6Kernel{2}() + shape = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); density=1000.0) + fluid = InitialCondition(; coordinates=shape.coordinates, density=density_field, + particle_spacing) + state_equation = StateEquationCole(; sound_speed=10.0, reference_density=1000.0, + exponent=1) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=SummationDensity(), + state_equation, correction) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + semi = ode.p.semi + system = first(semi.systems) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + + return (; system, semi, v_ode, u_ode, particle_spacing) +end + +function setup_continuity_density(n) + particle_spacing = 1.0 / n + smoothing_length = 2.0 * particle_spacing + smoothing_kernel = WendlandC6Kernel{2}() + shape = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); density=1000.0) + fluid = InitialCondition(; coordinates=shape.coordinates, density=density_field, + particle_spacing) + state_equation = StateEquationCole(; sound_speed=10.0, reference_density=1000.0, + exponent=1) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), + state_equation) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + semi = ode.p.semi + system = first(semi.systems) + + return (; system, semi, v_ode, u_ode, particle_spacing) +end + +function sample_regions(coordinates, support, n) + min_x = minimum(view(coordinates, 1, :)) + boundary = [particle + for particle in axes(coordinates, 2) + if isapprox(coordinates[1, particle], min_x; atol=eps()) && + support < coordinates[2, particle] < 1.0 - support] + isempty(boundary) && error("resolution $n has no particles in the boundary sample") + interior = [particle + for particle in axes(coordinates, 2) + if support < coordinates[1, particle] < 1.0 - support && + support < coordinates[2, particle] < 1.0 - support] + isempty(interior) && error("resolution $n has no particles in the interior sample") + + return (; boundary, interior) +end + +function operator_errors(n, correction) + setup = setup_operator(n, correction) + (; system, semi, v_ode, u_ode, particle_spacing) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + n_particles = TrixiParticles.nparticles(system) + values = [field(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + exact_gradients = zeros(2, n_particles) + for particle in TrixiParticles.eachparticle(system) + exact_gradients[:, + particle] = field_gradient(SVector{2}(view(coordinates, :, + particle))) + end + + interpolation = zeros(n_particles) + kernel_coefficient = zeros(n_particles) + direct_gradient = zeros(2, n_particles) + difference_gradient = zeros(2, n_particles) + GC.@preserve v_ode u_ode begin + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + semi) do particle, neighbor, pos_diff, + distance + pos_diff_ = SVector(pos_diff) + volume = TrixiParticles.hydrodynamic_mass(system, neighbor) / + TrixiParticles.current_density(v, system, neighbor) + kernel = TrixiParticles.smoothing_kernel(system, distance, particle) + gradient = TrixiParticles.smoothing_kernel_grad(system, pos_diff_, distance, + particle) + + interpolation[particle] += volume * values[neighbor] * kernel + kernel_coefficient[particle] += volume * kernel + for dimension in 1:2 + direct_gradient[dimension, + particle] += volume * values[neighbor] * gradient[dimension] + difference_gradient[dimension, + particle] += volume * + (values[neighbor] - values[particle]) * + gradient[dimension] + end + end + end + + support = TrixiParticles.compact_support(system, system) + (; boundary, interior) = sample_regions(coordinates, support, n) + normalized_interpolation = interpolation ./ kernel_coefficient + + function sample_errors(particles) + exact_values = values[particles] + + interpolation_error = normalized_l2(interpolation[particles], exact_values) + shepard_error = normalized_l2(normalized_interpolation[particles], exact_values) + direct_error = normalized_l2(direct_gradient[:, particles], + exact_gradients[:, particles]) + difference_error = normalized_l2(difference_gradient[:, particles], + exact_gradients[:, particles]) + + return (; interpolation_error, shepard_error, direct_error, difference_error) + end + + return (; particle_spacing, boundary=sample_errors(boundary), + interior=sample_errors(interior)) +end + +function summation_density_errors(n, correction) + setup = setup_summation_density(n, correction) + (; system, semi, v_ode, u_ode, particle_spacing) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + density = GC.@preserve v_ode u_ode begin + collect(TrixiParticles.current_density(v, system)) + end + exact_density = [density_field(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + support = TrixiParticles.compact_support(system, system) + regions = sample_regions(coordinates, support, n) + + function sample_error(particles) + return normalized_l2(density[particles], exact_density[particles]) + end + + return (; particle_spacing, boundary=sample_error(regions.boundary), + interior=sample_error(regions.interior)) +end + +function reinitialized_density_errors(n) + setup = setup_continuity_density(n) + (; system, semi, v_ode, u_ode, particle_spacing) = setup + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + TrixiParticles.reinit_density!(system, v, u, v_ode, u_ode, semi) + density = GC.@preserve v_ode u_ode begin + collect(TrixiParticles.current_density(v, system)) + end + exact_density = [density_field(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + support = TrixiParticles.compact_support(system, system) + regions = sample_regions(coordinates, support, n) + + function sample_error(particles) + return normalized_l2(density[particles], exact_density[particles]) + end + + return (; particle_spacing, boundary=sample_error(regions.boundary), + interior=sample_error(regions.interior)) +end + +function normalized_l2(approximation, exact) + return norm(approximation - exact) / norm(exact) +end + +function correction_name(correction) + isnothing(correction) && return :none + correction isa ShepardKernelCorrection && return :shepard + correction isa KernelCorrection && return :kernel + correction isa GradientCorrection && return :gradient + correction isa BlendedGradientCorrection && return :blended + correction isa MixedKernelGradientCorrection && return :mixed + correction isa CorrectionConfiguration && return :shepard_mixed + error("unsupported correction $(typeof(correction))") +end + +function append_result!(results, previous, method, operator, region, resolution, spacing, + error) + key = (method, operator, region) + order = if haskey(previous, key) + previous_spacing, previous_error = previous[key] + log(previous_error / error) / log(previous_spacing / spacing) + else + NaN + end + push!(results, (; method, operator, region, resolution, spacing, error, order)) + previous[key] = (spacing, error) + return results +end + +""" + run_convergence(; resolutions=(12, 24, 48, 96)) + +Measure local correction and pressure-acceleration scaling on self-similar regular patches with +fixed `h / Δx`. Boundary and symmetric-interior samples are reported separately. These +measurements are not convergence rates of the complete SPH discretization. +""" +function run_convergence(; resolutions=(12, 24, 48, 96)) + corrections = (nothing, KernelCorrection(), GradientCorrection(), + BlendedGradientCorrection(0.5), MixedKernelGradientCorrection()) + results = NamedTuple[] + previous = Dict{Tuple{Symbol, Symbol, Symbol}, Tuple{Float64, Float64}}() + + for resolution in resolutions + for correction in corrections + errors = operator_errors(resolution, correction) + method = correction_name(correction) + for region in (:boundary, :interior) + region_errors = getproperty(errors, region) + append_result!(results, previous, method, :difference_gradient, region, + resolution, errors.particle_spacing, + region_errors.difference_error) + + if method in (:none, :kernel, :mixed) + append_result!(results, previous, method, :direct_gradient, region, + resolution, errors.particle_spacing, + region_errors.direct_error) + end + + if method == :none + append_result!(results, previous, :none, :interpolation, region, + resolution, errors.particle_spacing, + region_errors.interpolation_error) + append_result!(results, previous, :shepard, :interpolation, region, + resolution, errors.particle_spacing, + region_errors.shepard_error) + end + end + end + + for (method, correction) in ((:none, nothing), + (:shepard, ShepardKernelCorrection())) + errors = summation_density_errors(resolution, correction) + for region in (:boundary, :interior) + append_result!(results, previous, method, :summation_density, region, + resolution, errors.particle_spacing, + getproperty(errors, region)) + end + end + + reinitialization_errors = reinitialized_density_errors(resolution) + for region in (:boundary, :interior) + append_result!(results, previous, :shepard, :density_reinitialization, region, + resolution, reinitialization_errors.particle_spacing, + getproperty(reinitialization_errors, region)) + end + + resolution < 24 && continue + + summation_corrections = (nothing, ShepardKernelCorrection(), KernelCorrection(), + GradientCorrection(), BlendedGradientCorrection(0.5), + MixedKernelGradientCorrection(), + CorrectionConfiguration(; + density=ShepardKernelCorrection(), + gradient=MixedKernelGradientCorrection())) + continuity_corrections = (nothing, KernelCorrection(), GradientCorrection(), + BlendedGradientCorrection(0.5), + MixedKernelGradientCorrection()) + pressure_cases = ((:pressure_summation, :constant_pressure_summation, + SummationDensity(), + TrixiParticles.pressure_acceleration_summation_density, + summation_corrections), + (:pressure_interparticle_summation, + :constant_pressure_interparticle_summation, SummationDensity(), + TrixiParticles.inter_particle_averaged_pressure, + summation_corrections), + (:pressure_continuity, :constant_pressure_continuity, + ContinuityDensity(), + TrixiParticles.pressure_acceleration_continuity_density, + continuity_corrections), + (:pressure_interparticle_continuity, + :constant_pressure_interparticle_continuity, ContinuityDensity(), + TrixiParticles.inter_particle_averaged_pressure, + continuity_corrections), + (:pressure_tensile_positive, :constant_pressure_tensile, + ContinuityDensity(), TrixiParticles.tensile_instability_control, + (nothing,))) + + for (operator, constant_operator, density_calculator, + pressure_formulation, corrections) in pressure_cases + + for correction in corrections + errors = pressure_operator_errors(resolution, correction, + density_calculator, + pressure_formulation) + method = correction_name(correction) + for region in (:boundary, :interior) + region_errors = getproperty(errors, region) + append_result!(results, previous, method, operator, region, resolution, + errors.particle_spacing, region_errors.manufactured) + append_result!(results, previous, method, constant_operator, region, + resolution, errors.particle_spacing, + region_errors.constant) + end + end + end + end + + return results +end + +function print_report(results; io=stdout) + println(io, "| Method | Operator | Region | N | L2 error | Observed scaling |") + println(io, "|:--|:--|:--|--:|--:|--:|") + for result in results + order = isnan(result.order) ? "-" : @sprintf("%.3f", result.order) + @printf(io, "| %s | %s | %s | %d | %.6e | %s |\n", result.method, + result.operator, result.region, result.resolution, result.error, order) + end + return results +end + +function write_csv(filename, results) + directory = dirname(filename) + isempty(directory) || mkpath(directory) + open(filename, "w") do io + println(io, + "method,operator,region,resolution,spacing,l2_error,observed_scaling") + for result in results + println(io, + join((result.method, result.operator, result.region, result.resolution, + result.spacing, result.error, result.order), ',')) + end + end + return filename +end + +end # module CorrectionConvergence + +if abspath(PROGRAM_FILE) == @__FILE__ + results = CorrectionConvergence.run_convergence() + CorrectionConvergence.print_report(results) + output_file = joinpath("out", "correction_convergence.csv") + CorrectionConvergence.write_csv(output_file, results) + println("\nWrote $output_file") +end diff --git a/validation/corrections/surface_pressure_difference.jl b/validation/corrections/surface_pressure_difference.jl new file mode 100644 index 0000000000..59c0fe6f10 --- /dev/null +++ b/validation/corrections/surface_pressure_difference.jl @@ -0,0 +1,158 @@ +module SurfacePressureDifferenceValidation + +using TrixiParticles +using LinearAlgebra: norm +using Printf: @printf + +function pressure_field(position) + x, y = position + return x * (1 - x) * y * (1 - y) * (1 + 0.2x + 0.1y) +end + +function pressure_gradient(position) + x, y = position + factor = 1 + 0.2x + 0.1y + return SVector((1 - 2x) * y * (1 - y) * factor + + 0.2x * (1 - x) * y * (1 - y), + x * (1 - x) * (1 - 2y) * factor + + 0.1x * (1 - x) * y * (1 - y)) +end + +function errors(n, correction) + particle_spacing = 1.0 / n + smoothing_length = 2.0 * particle_spacing + smoothing_kernel = WendlandC6Kernel{2}() + fluid = RectangularShape(particle_spacing, (n, n), (0.0, 0.0); density=1000.0) + state_equation = StateEquationCole(; sound_speed=10.0, reference_density=1000.0, + exponent=1) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), + state_equation, + gradient_correction=correction, + surface_method=ColorfieldSurfaceDetection(ideal_density_threshold=0.9), + surface_pressure=SurfacePressureDifference(), + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system; parallelization_backend=SerialBackend()) + ode = semidiscretize(semi, (0.0, 1.0); reset_threads=false) + v_ode = Array(ode.u0.x[1]) + u_ode = Array(ode.u0.x[2]) + semi = ode.p.semi + system = first(semi.systems) + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + coordinates = Array(TrixiParticles.current_coordinates(u, system)) + n_particles = TrixiParticles.nparticles(system) + pressure = [pressure_field(SVector{2}(view(coordinates, :, particle))) + for particle in TrixiParticles.eachparticle(system)] + acceleration = zeros(2, n_particles) + constant_acceleration = zeros(2, n_particles) + + GC.@preserve v_ode u_ode begin + TrixiParticles.foreach_point_neighbor(system, system, coordinates, coordinates, + semi) do particle, neighbor, pos_diff, + distance + m_a = TrixiParticles.hydrodynamic_mass(system, particle) + m_b = TrixiParticles.hydrodynamic_mass(system, neighbor) + rho_a = TrixiParticles.current_density(v, system, particle) + rho_b = TrixiParticles.current_density(v, system, neighbor) + W_a = TrixiParticles.smoothing_kernel_grad(system, SVector(pos_diff), distance, + particle) + gradient_correction = TrixiParticles.correction_gradient(system.correction) + acceleration[:, + particle] .+= TrixiParticles.pressure_acceleration(system, + system, + particle, + neighbor, + m_a, m_b, + pressure[particle], + pressure[neighbor], + rho_a, rho_b, + pos_diff, + distance, W_a, + gradient_correction) + constant_acceleration[:, + particle] .+= TrixiParticles.pressure_acceleration(system, + system, + particle, + neighbor, + m_a, + m_b, + 2.0, + 2.0, + rho_a, + rho_b, + pos_diff, + distance, + W_a, + gradient_correction) + end + end + + exact = zeros(2, n_particles) + for particle in TrixiParticles.eachparticle(system) + exact[:, + particle] = -pressure_gradient(SVector{2}(view(coordinates, :, particle))) / + TrixiParticles.current_density(v, system, particle) + end + + support = TrixiParticles.compact_support(system, system) + min_x = minimum(view(coordinates, 1, :)) + boundary = [particle + for particle in axes(coordinates, 2) + if isapprox(coordinates[1, particle], min_x; atol=eps()) && + support < coordinates[2, particle] < 1.0 - support] + interior = [particle + for particle in axes(coordinates, 2) + if 2support < coordinates[1, particle] < 1.0 - 2support && + 2support < coordinates[2, particle] < 1.0 - 2support] + + relative_error(particles) = norm(acceleration[:, particles] - exact[:, particles]) / + norm(exact[:, particles]) + constant_error(particles) = norm(constant_acceleration[:, particles]) / + sqrt(length(particles)) + mass = TrixiParticles.hydrodynamic_mass(system, 1) + momentum_residual = norm(mass * vec(sum(acceleration; dims=2))) + + return (; spacing=particle_spacing, boundary=relative_error(boundary), + interior=relative_error(interior), + constant_boundary=constant_error(boundary), momentum_residual) +end + +function run(; resolutions=(24, 48, 96)) + results = NamedTuple[] + for (method, correction) in ((:gradient, GradientCorrection()), + (:mixed, MixedKernelGradientCorrection())) + previous = nothing + println("\n$method") + println("| N | Boundary error | Boundary order | Interior error | Interior order | " * + "Constant boundary | Momentum residual |") + println("|--:|---------------:|---------------:|---------------:|---------------:|" * + "------------------:|------------------:|") + for n in resolutions + result = errors(n, correction) + boundary_order = isnothing(previous) ? NaN : + log(previous.boundary / result.boundary) / + log(previous.spacing / result.spacing) + interior_order = isnothing(previous) ? NaN : + log(previous.interior / result.interior) / + log(previous.spacing / result.spacing) + @printf("| %d | %.6e | %s | %.6e | %s | %.6e | %.6e |\n", n, + result.boundary, + isnan(boundary_order) ? "-" : string(boundary_order), + result.interior, + isnan(interior_order) ? "-" : string(interior_order), + result.constant_boundary, result.momentum_residual) + push!(results, (; method, n, result..., boundary_order, interior_order)) + previous = result + end + end + return results +end + +end + +if abspath(PROGRAM_FILE) == @__FILE__ + SurfacePressureDifferenceValidation.run() +end