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..63797d8df8 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")] diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index ee3bf532b5..54c923000f 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 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/io/io.jl b/src/io/io.jl index 692dfd8528..5d1c0b980b 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -306,6 +306,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, 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..f8c6d8e624 100644 --- a/src/schemes/boundary/wall_boundary/system.jl +++ b/src/schemes/boundary/wall_boundary/system.jl @@ -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..84200bd066 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -8,7 +8,8 @@ acceleration=ntuple(_ -> 0.0, NDIMS), surface_tension=nothing, surface_normal_method=nothing, buffer_size=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 +39,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` @@ -93,7 +98,8 @@ 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, reference_particle_spacing=0.0, color_value=1) @@ -110,6 +116,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 @@ -127,7 +138,7 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) 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 +146,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) @@ -251,7 +262,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 +312,57 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, 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_normal!(system, system.surface_normal_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/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index e4ef03584f..596d322fe7 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -4,7 +4,9 @@ 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) System for particles of a fluid. The system employs implicit incompressible SPH (IISPH), iteratively solving a linear system @@ -30,6 +32,10 @@ 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. """ struct ImplicitIncompressibleSPHSystem{NDIMS, ELTYPE <: Real, ARRAY1D, ARRAY2D, IC, K, V, PF, C} <: AbstractFluidSystem{NDIMS} @@ -71,10 +77,17 @@ 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) 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) diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index b6114c0bc1..253d7c2b41 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -93,6 +93,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) 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..e84d56fc20 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -6,7 +6,8 @@ pressure_acceleration=nothing, shifting_technique=nothing, buffer_size=nothing, - correction=nothing, source_terms=nothing, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, reference_particle_spacing=0.0, color_value=1)) @@ -41,7 +42,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` @@ -96,7 +101,8 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, pressure_acceleration=nothing, shifting_technique=nothing, buffer_size=nothing, - correction=nothing, source_terms=nothing, + correction=nothing, density_correction=nothing, + gradient_correction=nothing, source_terms=nothing, surface_tension=nothing, surface_normal_method=nothing, reference_particle_spacing=0, color_value=1) buffer = isnothing(buffer_size) ? nothing : @@ -112,6 +118,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,7 +136,7 @@ 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 @@ -141,7 +152,7 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, 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, @@ -243,7 +254,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) @@ -322,15 +335,17 @@ end function update_pressure!(system::WeaklyCompressibleSPHSystem, v, u, v_ode, u_ode, semi, t) (; density_calculator, correction, surface_normal_method, surface_tension) = system + density_correction = correction_density(correction) + gradient_correction = correction_gradient(correction) - 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 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) @@ -369,13 +384,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 +412,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..c5dfb2cdfd 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -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) 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/semidiscretization.jl b/test/general/semidiscretization.jl index 5b92d133ba..9b7ddf3cb9 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -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..1563dc2d15 100644 --- a/test/schemes/fluid/pressure_acceleration.jl +++ b/test/schemes/fluid/pressure_acceleration.jl @@ -12,6 +12,47 @@ @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 + + @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) + end + @testset verbose=true "Illegal Inputs" begin correction_dict_1 = Dict( "KernelCorrection" => KernelCorrection(), 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