diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index a748be851..5e093c7b7 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -25,7 +25,10 @@ concurrency: jobs: build-docs: name: Build and Deploy Documentation - runs-on: ubuntu-latest + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, windows-latest] steps: - name: Check out project uses: actions/checkout@v4 @@ -40,7 +43,11 @@ jobs: - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy + if: matrix.os == 'ubuntu-latest' env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key run: julia --project=docs --color=yes docs/make.jl + - name: Just Build + if: matrix.os != 'ubuntu-latest' + run: julia --project=docs --color=yes docs/make.jl diff --git a/README.md b/README.md index 0fe4741b9..c90ed5b31 100644 --- a/README.md +++ b/README.md @@ -10,8 +10,14 @@ [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10797541.svg)](https://zenodo.org/doi/10.5281/zenodo.10797541) -**TrixiParticles.jl** is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in [Julia](https://julialang.org). -A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing. +**TrixiParticles.jl** is a high-performance numerical simulation framework for particle-based methods, focused on the simulation of complex multiphysics problems, and written in [Julia](https://julialang.org). + +TrixiParticles.jl focuses on the following use cases: +- Accurate and efficient physics-based modelling of complex multiphysics problems. +- Development of new particle-based methods and models. +- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work. + +It offers intuitive configuration, robust pre- and post-processing, and vendor-agnostic GPU-support based on the Julia package [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl). [![YouTube](https://github.com/user-attachments/assets/dc2be627-a799-4bfd-9226-2077f737c4b0)](https://www.youtube.com/watch?v=V7FWl4YumcA&t=4667s) diff --git a/docs/make.jl b/docs/make.jl index 3e3411f14..ac541b0a1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -19,9 +19,12 @@ function copy_file(filename, replaces...; content = read(source_path, String) content = replace(content, replaces...) + # Use `replace` to make sure the path uses forward slashes for URLs + filename_url = replace(filename, "\\" => "/") + header = """ ```@meta - EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename" + EditURL = "https://github.com/trixi-framework/TrixiParticles.jl/blob/main/$filename_url" ``` """ content = header * content @@ -117,7 +120,8 @@ makedocs(sitename="TrixiParticles.jl", "Preprocessing" => [ "Sampling of Geometries" => joinpath("preprocessing", "preprocessing.md") ], - "Components" => [ + "GPU Support" => "gpu.md", + "API Reference" => [ "Overview" => "overview.md", "General" => [ "Semidiscretization" => joinpath("general", "semidiscretization.md"), diff --git a/docs/src/gpu.md b/docs/src/gpu.md new file mode 100644 index 000000000..8e8a12708 --- /dev/null +++ b/docs/src/gpu.md @@ -0,0 +1,61 @@ +# GPU Support + +GPU support is still an experimental feature that is actively being worked on. +As of now, the [`WeaklyCompressibleSPHSystem`](@ref) and the [`BoundarySPHSystem`](@ref) +are supported on GPUs. +We have tested this on GPUs by Nvidia and AMD. + +To run a simulation on a GPU, we need to use the [`FullGridCellList`](@ref) +as cell list for the [`GridNeighborhoodSearch`](@ref). +This cell list requires a bounding box for the domain, unlike the default cell list, which +uses an unbounded domain. +For simulations that are bounded by a closed tank, we can use the boundary of the tank +to obtain the bounding box as follows. +```jldoctest gpu; output=false, setup=:(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing)) +search_radius = TrixiParticles.compact_support(smoothing_kernel, smoothing_length) +min_corner = minimum(tank.boundary.coordinates, dims=2) .- search_radius +max_corner = maximum(tank.boundary.coordinates, dims=2) .+ search_radius +cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, max_corner) + +# output +PointNeighbors.FullGridCellList{PointNeighbors.DynamicVectorOfVectors{Int32, Matrix{Int32}, Vector{Int32}, Base.RefValue{Int32}}, Nothing, SVector{2, Float64}, SVector{2, Float64}}(Vector{Int32}[], nothing, [-0.24500000000000002, -0.24500000000000002], [1.245, 1.245]) +``` + +We then need to pass this cell list to the neighborhood search and the neighborhood search +to the [`Semidiscretization`](@ref). +```jldoctest gpu; output=false +semi = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(; cell_list)) + +# output +┌──────────────────────────────────────────────────────────────────────────────────────────────────┐ +│ Semidiscretization │ +│ ══════════════════ │ +│ #spatial dimensions: ………………………… 2 │ +│ #systems: ……………………………………………………… 2 │ +│ neighborhood search: ………………………… GridNeighborhoodSearch │ +│ total #particles: ………………………………… 636 │ +└──────────────────────────────────────────────────────────────────────────────────────────────────┘ +``` + +At this point, we should run the simulation and make sure that it still works and that +the bounding box is large enough. +For some simulations where particles move outside the initial tank coordinates, +for example when the tank is not closed or when the tank is moving, an appropriate +bounding box has to be specified. + +Then, we only need to specify the data type that is used for the simulation. +On an Nvidia GPU, we specify: +```julia +using CUDA +ode = semidiscretize(semi, tspan, data_type=CuArray) +``` +On an AMD GPU, we use: +```julia +using AMDGPU +ode = semidiscretize(semi, tspan, data_type=ROCArray) +``` +Then, we can run the simulation as usual. +All data is transferred to the GPU during initialization and all loops over particles +and their neighbors will be executed on the GPU as kernels generated by KernelAbstractions.jl. +Data is only copied to the CPU for saving VTK files via the [`SolutionSavingCallback`](@ref). diff --git a/docs/src/index.md b/docs/src/index.md index 42c2c42c6..27c1e9d47 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,6 +1,14 @@ # TrixiParticles.jl -TrixiParticles.jl is a numerical simulation framework designed for particle-based numerical methods, with an emphasis on multiphysics applications, written in Julia. A primary goal of the framework is to be user-friendly for engineering, science, and educational purposes. In addition to its extensible design and optimized implementation, we prioritize the user experience, including installation, pre- and postprocessing. Its features include: +**TrixiParticles.jl** is a high-performance particle simulation framework designed to overcome challenges of particle-based numerical methods in multiphysics applications. Existing frameworks often lack user-friendliness, involve complex configuration, and are not easily extensible for development of new methods. In the future we also want to provide seamless scalability from CPU to Exascale-level computing with GPU support. **TrixiParticles.jl** addresses these limitations with an intuitive interface, straightforward configuration, and an extensible design, facilitating efficient simulation setup and execution. + +TrixiParticles.jl focuses on the following use cases: + +- Development of new particle-based methods and models. By providing an extensible architecture to incorporate additional particle methods easily and not focusing on a single model or numerical method. +- Accurate, reliable and efficient physics-based modelling of complex multiphysics problems by providing a flexible configuration system, tools, high performance and a wide range of validation and test cases. +- Easy setup of accessible simulations for educational purposes, including student projects, coursework, and thesis work through extensive documentation, community engagement and readable configuration files. + +Its features include: ## Features - Incompressible Navier-Stokes diff --git a/docs/src/overview.md b/docs/src/overview.md index a113d1f9f..5d93c3732 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -1,10 +1,16 @@ # Overview +The actual API reference is not listed on a single page, like in most Julia packages, +but instead is split into multiple sections that follow a similar structure +as the code files themselves. +In these sections, API docs are combined with explanations of the theoretical background +of these methods. + The following page gives a rough overview of important parts of the code. ## Program flow To initiate a simulation, the goal is to solve an ordinary differential equation, for example, -by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then +by employing the time integration schemes provided by OrdinaryDiffEq.jl. These schemes are then utilized to integrate ``\mathrm{d}u/\mathrm{d}t`` and ``\mathrm{d}v/\mathrm{d}t``, where ``u`` represents the particles' positions and ``v`` their properties such as velocity and density. During a single time step or an intermediate step of the time integration scheme, the functions diff --git a/docs/src/reference-pointneighbors.md b/docs/src/reference-pointneighbors.md index f8f8c5ac7..ae184bb2a 100644 --- a/docs/src/reference-pointneighbors.md +++ b/docs/src/reference-pointneighbors.md @@ -1,4 +1,4 @@ -# PointNeighbors.jl API +# [PointNeighbors.jl API](@id pointneighbors) ```@meta CurrentModule = PointNeighbors diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index 29384b6d1..7beb184ed 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -55,13 +55,17 @@ of the boundary particle ``b``. ### Hydrodynamic density of dummy particles -We provide five options to compute the boundary density and pressure, determined by the `density_calculator`: +We provide six options to compute the boundary density and pressure, determined by the `density_calculator`: 1. (Recommended) With [`AdamiPressureExtrapolation`](@ref), the pressure is extrapolated from the pressure of the fluid according to [Adami et al., 2012](@cite Adami2012), and the density is obtained by applying the inverse of the state equation. This option usually yields the best results of the options listed here. -2. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, +2. (Only relevant for FSI) With [`BernoulliPressureExtrapolation`](@ref), the pressure is extrapolated from the + pressure similar to the [`AdamiPressureExtrapolation`](@ref), but a relative velocity-dependent pressure part + is calculated between moving solids and fluids, which increases the boundary pressure in areas prone to + penetrations. +3. With [`SummationDensity`](@ref), the density is calculated by summation over the neighboring particles, and the pressure is computed from the density with the state equation. -3. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, +4. With [`ContinuityDensity`](@ref), the density is integrated from the continuity equation, and the pressure is computed from the density with the state equation. Note that this causes a gap between fluid and boundary where the boundary is initialized without any contact to the fluid. This is due to overestimation of the boundary density @@ -69,10 +73,10 @@ We provide five options to compute the boundary density and pressure, determined contact to the fluid. Therefore, in dam break simulations, there is a visible "step", even though the boundary is supposed to be flat. See also [dual.sphysics.org/faq/#Q_13](https://dual.sphysics.org/faq/#Q_13). -4. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure +5. With [`PressureZeroing`](@ref), the density is set to the reference density and the pressure is computed from the density with the state equation. This option is not recommended. The other options yield significantly better results. -5. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure +6. With [`PressureMirroring`](@ref), the density is set to the reference density. The pressure is not used. Instead, the fluid pressure is mirrored as boundary pressure in the momentum equation. This option is not recommended due to stability issues. See [`PressureMirroring`](@ref) @@ -93,7 +97,20 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den AdamiPressureExtrapolation ``` -#### 4. [`PressureZeroing`](@ref) +#### 2. [`BernoulliPressureExtrapolation`](@ref) +Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation: +```math +p_b = \frac{\sum_f (p_f + \frac{1}{2} \, \rho_{\text{neighbor}} \left( \frac{ (\mathbf{v}_f - \mathbf{v}_{\text{body}}) \cdot (\mathbf{x}_f - \mathbf{x}_{\text{neighbor}}) }{ \left\| \mathbf{x}_f - \mathbf{x}_{\text{neighbor}} \right\| } \right)^2 \times \text{factor} +\rho_f (\bm{g} - \bm{a}_b) \cdot \bm{r}_{bf}) W(\Vert r_{bf} \Vert, h)}{\sum_f W(\Vert r_{bf} \Vert, h)} +``` +where ``\mathbf{v}_f`` is the velocity of the fluid and ``\mathbf{v}_{\text{body}}`` is the velocity of the body. +This adjustment provides a higher boundary pressure for solid bodies moving with a relative velocity to the fluid to prevent penetration. +This modification is original and not derived from any literature source. + +```@docs + BernoulliPressureExtrapolation +``` + +#### 5. [`PressureZeroing`](@ref) This is the simplest way to implement dummy boundary particles. The density of each particle is set to the reference density and the pressure to the @@ -102,7 +119,7 @@ reference pressure (the corresponding pressure to the reference density by the s PressureZeroing ``` -#### 5. [`PressureMirroring`](@ref) +#### 6. [`PressureMirroring`](@ref) Instead of calculating density and pressure for each boundary particle, we modify the momentum equation, diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index af558133f..c8778f485 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -48,7 +48,7 @@ pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_densi # Shift pipe walls in negative x-direction for the inflow pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers -n_buffer_particles = 4 * pipe.n_particles_per_dimension[2] +n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]^(ndims(pipe.fluid) - 1) # ========================================================================================== # ==== Fluid @@ -77,14 +77,16 @@ fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothi # ========================================================================================== # ==== Open Boundary -function velocity_function(pos, t) + +function velocity_function2d(pos, t) # Use this for a time-dependent inflow velocity # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) return SVector(prescribed_velocity, 0.0) end -inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction, +plane_in = ([0.0, 0.0], [0.0, domain_size[2]]) +inflow = InFlow(; plane=plane_in, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, @@ -92,9 +94,10 @@ open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system, buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) -outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]), +plane_out = ([domain_size[1], 0.0], [domain_size[1], domain_size[2]]) +outflow = OutFlow(; plane=plane_out, flow_direction, open_boundary_layers, density=fluid_density, particle_spacing) @@ -103,7 +106,7 @@ open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system, buffer_size=n_buffer_particles, reference_density=fluid_density, reference_pressure=pressure, - reference_velocity=velocity_function) + reference_velocity=velocity_function2d) # ========================================================================================== # ==== Boundary diff --git a/examples/fluid/pipe_flow_3d.jl b/examples/fluid/pipe_flow_3d.jl new file mode 100644 index 000000000..fb2cc7f0a --- /dev/null +++ b/examples/fluid/pipe_flow_3d.jl @@ -0,0 +1,45 @@ +# 3D channel flow simulation with open boundaries. + +using TrixiParticles + +# ========================================================================================== +# ==== Resolution +particle_spacing = 0.05 + +# Make sure that the kernel support of fluid particles at a boundary is always fully sampled +boundary_layers = 3 + +# Make sure that the kernel support of fluid particles at an open boundary is always +# fully sampled. +# Note: Due to the dynamics at the inlets and outlets of open boundaries, +# it is recommended to use `open_boundary_layers > boundary_layers` +open_boundary_layers = 6 + +# ========================================================================================== +# ==== Experiment Setup +tspan = (0.0, 2.0) + +function velocity_function3d(pos, t) + # Use this for a time-dependent inflow velocity + # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0) + + return SVector(prescribed_velocity, 0.0, 0.0) +end + +domain_size = (1.0, 0.4, 0.4) + +boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers, + domain_size[2], domain_size[3]) + +flow_direction = [1.0, 0.0, 0.0] + +# setup simulation +trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "pipe_flow_2d.jl"), + domain_size=domain_size, boundary_size=boundary_size, + flow_direction=flow_direction, faces=(false, false, true, true, true, true), + tspan=tspan, smoothing_kernel=WendlandC2Kernel{3}(), + reference_velocity=velocity_function3d, + plane_in=([0.0, 0.0, 0.0], [0.0, domain_size[2], 0.0], + [0.0, 0.0, domain_size[3]]), + plane_out=([domain_size[1], 0.0, 0.0], [domain_size[1], domain_size[2], 0.0], + [domain_size[1], 0.0, domain_size[3]])) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 6a2397a4b..a06056f8f 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -63,7 +63,7 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid, fluid_density_calculator, # ========================================================================================== # ==== Boundary -boundary_density_calculator = AdamiPressureExtrapolation() +boundary_density_calculator = BernoulliPressureExtrapolation() boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, state_equation=state_equation, boundary_density_calculator, diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index bcdce98ac..54bfb7957 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -3,6 +3,7 @@ module TrixiParticles using Reexport: @reexport using Adapt: Adapt +using Base: @propagate_inbounds using CSV: CSV using Dates using DataFrames: DataFrame @@ -21,7 +22,8 @@ using Printf: @printf, @sprintf using RecipesBase: RecipesBase, @series using Random: seed! using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!, - get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem + get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, + RecursiveArrayTools @reexport using StaticArrays: SVector using StaticArrays: @SMatrix, SMatrix, setindex using StrideArrays: PtrArray, StaticInt @@ -44,6 +46,7 @@ include("callbacks/callbacks.jl") include("general/general.jl") include("setups/setups.jl") include("schemes/schemes.jl") +include("refinement/refinement.jl") # Note that `semidiscretization.jl` depends on the system types and has to be # included separately. `gpu.jl` in turn depends on the semidiscretization type. @@ -51,6 +54,7 @@ include("general/semidiscretization.jl") include("general/gpu.jl") include("visualization/write2vtk.jl") include("visualization/recipes_plots.jl") +include("refinement/resize.jl") export Semidiscretization, semidiscretize, restart_with! export InitialCondition @@ -69,7 +73,9 @@ export ArtificialViscosityMonaghan, ViscosityAdami, ViscosityMorris export DensityDiffusion, DensityDiffusionMolteniColagrossi, DensityDiffusionFerrari, DensityDiffusionAntuono export BoundaryModelMonaghanKajtar, BoundaryModelDummyParticles, AdamiPressureExtrapolation, - PressureMirroring, PressureZeroing, BoundaryModelLastiwka + PressureMirroring, PressureZeroing, BoundaryModelLastiwka, + BernoulliPressureExtrapolation +export ParticleRefinement, SpatialRefinementCriterion export BoundaryMovement export examples_dir, validation_dir, trixi_include export trixi2vtk diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 1aefdff72..3d8f5701a 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -31,3 +31,4 @@ include("density_reinit.jl") include("post_process.jl") include("stepsize.jl") include("update.jl") +include("refinement.jl") diff --git a/src/callbacks/refinement.jl b/src/callbacks/refinement.jl new file mode 100644 index 000000000..130de8756 --- /dev/null +++ b/src/callbacks/refinement.jl @@ -0,0 +1,129 @@ +struct ParticleRefinementCallback{I} + interval::I +end + +function ParticleRefinementCallback(; interval::Integer=-1, dt=0.0) + if dt > 0 && interval !== -1 + throw(ArgumentError("setting both `interval` and `dt` is not supported")) + end + + # Update in intervals in terms of simulation time + if dt > 0 + interval = Float64(dt) + + # Update every time step (default) + elseif interval == -1 + interval = 1 + end + + refinement_callback = ParticleRefinementCallback(interval) + + if dt > 0 + # Add a `tstop` every `dt`, and save the final solution. + return PeriodicCallback(refinement_callback, dt, + initialize=initial_update!, + save_positions=(false, false)) + else + # The first one is the `condition`, the second the `affect!` + return DiscreteCallback(refinement_callback, refinement_callback, + initialize=initial_update!, + save_positions=(false, false)) + end +end + +# initialize +function initial_refinement!(cb, u, t, integrator) + # The `ParticleRefinementCallback` is either `cb.affect!` (with `DiscreteCallback`) + # or `cb.affect!.affect!` (with `PeriodicCallback`). + # Let recursive dispatch handle this. + + initial_refinement!(cb.affect!, u, t, integrator) +end + +function initial_refinement!(cb::ParticleRefinementCallback, u, t, integrator) + semi = integrator.p + + foreach_system(semi) do system + resize_refinement!(system) + end + + cb(integrator) +end + +# condition +function (refinement_callback::ParticleRefinementCallback)(u, t, integrator) + (; interval) = refinement_callback + + return condition_integrator_interval(integrator, interval) +end + +# affect +function (refinement_callback::ParticleRefinementCallback)(integrator) + t = integrator.t + semi = integrator.p + v_ode, u_ode = integrator.u.x + + # Update NHS + @trixi_timeit timer() "update nhs" update_nhs(u_ode, semi) + + # Basically `get_tmp_cache(integrator)` to write into in order to be non-allocating + # https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Caches + v_tmp, u_tmp = integrator.cache.tmp.x + + v_tmp .= v_ode + u_tmp .= u_ode + + refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + + resize!(integrator, (length(v_ode), length(u_ode))) + + # Tell OrdinaryDiffEq that u has been modified + u_modified!(integrator, true) + + return integrator +end + +Base.resize!(a::RecursiveArrayTools.ArrayPartition, sizes::Tuple) = resize!.(a.x, sizes) + +function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(interval=", (cb.affect!.interval), ")") +end + +function Base.show(io::IO, + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + print(io, "ParticleRefinementCallback(dt=", cb.affect!.affect!.interval, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, <:ParticleRefinementCallback}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect! + setup = [ + "interval" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end + +function Base.show(io::IO, ::MIME"text/plain", + cb::DiscreteCallback{<:Any, + <:PeriodicCallbackAffect{<:ParticleRefinementCallback}}) + @nospecialize cb # reduce precompilation time + + if get(io, :compact, false) + show(io, cb) + else + refinement_cb = cb.affect!.affect! + setup = [ + "dt" => refinement_cb.interval + ] + summary_box(io, "ParticleRefinementCallback", setup) + end +end diff --git a/src/general/corrections.jl b/src/general/corrections.jl index cf6c8bed6..c63c25a27 100644 --- a/src/general/corrections.jl +++ b/src/general/corrections.jl @@ -40,11 +40,11 @@ end k = correction.rho0 / rho_mean # Viscosity, pressure, surface_tension - return k, 1.0, k + return k, 1, k end @inline function free_surface_correction(correction, particle_system, rho_mean) - return 1.0, 1.0, 1.0 + return 1, 1, 1 end @doc raw""" @@ -220,7 +220,7 @@ function compute_correction_values!(system, kernel_correction_coefficient[particle] += volume * smoothing_kernel(system, distance) if distance > sqrt(eps()) - tmp = volume * smoothing_kernel_grad(system, pos_diff, distance) + tmp = volume * smoothing_kernel_grad(system, pos_diff, distance, particle) for i in axes(dw_gamma, 1) dw_gamma[i, particle] += tmp[i] end @@ -311,7 +311,7 @@ function compute_gradient_correction_matrix!(corr_matrix, neighborhood_search, pos_diff, distance volume = mass[neighbor] / density_fun(neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) iszero(grad_kernel) && return @@ -349,7 +349,7 @@ function compute_gradient_correction_matrix!(corr_matrix::AbstractArray, system, function compute_grad_kernel(correction, smoothing_kernel, pos_diff, distance, smoothing_length, system, particle) - return smoothing_kernel_grad(system, pos_diff, distance) + return smoothing_kernel_grad(system, pos_diff, distance, particle) end # Compute gradient of corrected kernel diff --git a/src/general/density_calculators.jl b/src/general/density_calculators.jl index 74a5f616d..e15eb33e0 100644 --- a/src/general/density_calculators.jl +++ b/src/general/density_calculators.jl @@ -23,15 +23,15 @@ difference of the coordinates, ``v_{ab} = v_a - v_b`` of the velocities of parti """ struct ContinuityDensity end -@inline function particle_density(v, system, particle) +@propagate_inbounds function particle_density(v, system, particle) particle_density(v, system.density_calculator, system, particle) end -@inline function particle_density(v, ::SummationDensity, system, particle) +@propagate_inbounds function particle_density(v, ::SummationDensity, system, particle) return system.cache.density[particle] end -@inline function particle_density(v, ::ContinuityDensity, system, particle) +@propagate_inbounds function particle_density(v, ::ContinuityDensity, system, particle) return v[end, particle] end @@ -66,7 +66,7 @@ function summation_density!(system, semi, u, u_ode, density; points=particles) do particle, neighbor, pos_diff, distance mass = hydrodynamic_mass(neighbor_system, neighbor) - density[particle] += mass * smoothing_kernel(system, distance) + density[particle] += mass * smoothing_kernel(system, distance, particle) end end end diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 600d573d5..f734a6191 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -73,14 +73,7 @@ function Semidiscretization(systems...; # Other checks might be added here later. check_configuration(systems) - sizes_u = [u_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_u = Tuple((sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) - for i in eachindex(sizes_u)) - sizes_v = [v_nvariables(system) * n_moving_particles(system) - for system in systems] - ranges_v = Tuple((sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) - for i in eachindex(sizes_v)) + ranges_v, ranges_u = ranges_vu(systems) # Create a tuple of n neighborhood searches for each of the n systems. # We will need one neighborhood search for each pair of systems. @@ -92,6 +85,16 @@ function Semidiscretization(systems...; return Semidiscretization(systems, ranges_u, ranges_v, searches) end +function ranges_vu(systems) + sizes_u = [u_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_u = [(sum(sizes_u[1:(i - 1)]) + 1):sum(sizes_u[1:i]) for i in eachindex(sizes_u)] + + sizes_v = [v_nvariables(system) * n_moving_particles(system) for system in systems] + ranges_v = [(sum(sizes_v[1:(i - 1)]) + 1):sum(sizes_v[1:i]) for i in eachindex(sizes_v)] + + return ranges_v, ranges_u +end + # Inline show function e.g. Semidiscretization(neighborhood_search=...) function Base.show(io::IO, semi::Semidiscretization) @nospecialize semi # reduce precompilation time @@ -134,8 +137,8 @@ function create_neighborhood_search(neighborhood_search, system, neighbor) end @inline function compact_support(system, neighbor) - (; smoothing_kernel, smoothing_length) = system - return compact_support(smoothing_kernel, smoothing_length) + (; smoothing_kernel) = system + return compact_support(smoothing_kernel, smoothing_length(system, 1)) # TODO end @inline function compact_support(system::OpenBoundarySPHSystem, neighbor) @@ -407,7 +410,9 @@ end function calculate_dt(v_ode, u_ode, cfl_number, semi::Semidiscretization) (; systems) = semi - return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) + return 1.0 + # TODO + # return minimum(system -> calculate_dt(v_ode, u_ode, cfl_number, system), systems) end function drift!(du_ode, v_ode, u_ode, semi, t) diff --git a/src/general/smoothing_kernels.jl b/src/general/smoothing_kernels.jl index 20268ebad..862e58ba1 100644 --- a/src/general/smoothing_kernels.jl +++ b/src/general/smoothing_kernels.jl @@ -3,7 +3,8 @@ abstract type SmoothingKernel{NDIMS} end @inline Base.ndims(::SmoothingKernel{NDIMS}) where {NDIMS} = NDIMS @inline function kernel_grad(kernel, pos_diff, distance, h) - distance < sqrt(eps()) && return zero(pos_diff) + # TODO Use `eps` relative to `h` to allow scaling of simulations + distance < sqrt(eps(typeof(h))) && return zero(pos_diff) return kernel_deriv(kernel, distance, h) / distance * pos_diff end @@ -98,7 +99,8 @@ end @inline compact_support(::GaussianKernel, h) = 3 * h @inline normalization_factor(::GaussianKernel{2}, h) = 1 / (pi * h^2) -@inline normalization_factor(::GaussianKernel{3}, h) = 1 / (pi^(3 / 2) * h^3) +# First convert `pi` to the type of `h` to preserve the type of `h` +@inline normalization_factor(::GaussianKernel{3}, h) = 1 / (oftype(h, pi)^(3 // 2) * h^3) @doc raw""" SchoenbergCubicSplineKernel{NDIMS}() @@ -136,8 +138,9 @@ struct SchoenbergCubicSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end @muladd @inline function kernel(kernel::SchoenbergCubicSplineKernel, r::Real, h) q = r / h - # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl - result = 1 / 4 * (2 - q)^3 + # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl. + # Use `//` to preserve the type of `q`. + result = 1 // 4 * (2 - q)^3 result = result - (q < 1) * (1 - q)^3 # Zero out result if q >= 2 @@ -151,7 +154,8 @@ end q = r * inner_deriv # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl - result = -3 / 4 * (2 - q)^2 + # Use `//` to preserve the type of `q`. + result = -3 // 4 * (2 - q)^2 result = result + 3 * (q < 1) * (1 - q)^2 # Zero out result if q >= 2 @@ -164,7 +168,8 @@ end @inline compact_support(::SchoenbergCubicSplineKernel, h) = 2 * h @inline normalization_factor(::SchoenbergCubicSplineKernel{1}, h) = 2 / 3h -@inline normalization_factor(::SchoenbergCubicSplineKernel{2}, h) = 10 / (7 * pi * h^2) +# `7 * pi` is always `Float64`. `pi * h^2 * 7` preserves the type of `h`. +@inline normalization_factor(::SchoenbergCubicSplineKernel{2}, h) = 10 / (pi * h^2 * 7) @inline normalization_factor(::SchoenbergCubicSplineKernel{3}, h) = 1 / (pi * h^3) @doc raw""" @@ -213,17 +218,19 @@ struct SchoenbergQuarticSplineKernel{NDIMS} <: SmoothingKernel{NDIMS} end # is a higher priority than exact precision. @fastpow @muladd @inline function kernel(kernel::SchoenbergQuarticSplineKernel, r::Real, h) q = r / h - q5_2 = (5 / 2 - q) - q3_2 = (3 / 2 - q) - q1_2 = (1 / 2 - q) + + # Preserve the type of `q` + q5_2 = (5 // 2 - q) + q3_2 = (3 // 2 - q) + q1_2 = (1 // 2 - q) # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl result = q5_2^4 - result = result - 5 * (q < 3 / 2) * q3_2^4 - result = result + 10 * (q < 1 / 2) * q1_2^4 + result = result - 5 * (q < 3 // 2) * q3_2^4 + result = result + 10 * (q < 1 // 2) * q1_2^4 # Zero out result if q >= 5/2 - result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result, zero(result)) + result = ifelse(q < 5 // 2, normalization_factor(kernel, h) * result, zero(result)) return result end @@ -232,17 +239,19 @@ end r::Real, h) inner_deriv = 1 / h q = r * inner_deriv - q5_2 = 5 / 2 - q - q3_2 = 3 / 2 - q - q1_2 = 1 / 2 - q + + # Preserve the type of `q` + q5_2 = 5 // 2 - q + q3_2 = 3 // 2 - q + q1_2 = 1 // 2 - q # We do not use `+=` or `-=` since these are not recognized by MuladdMacro.jl result = -4 * q5_2^3 - result = result + 20 * (q < 3 / 2) * q3_2^3 - result = result - 40 * (q < 1 / 2) * q1_2^3 + result = result + 20 * (q < 3 // 2) * q3_2^3 + result = result - 40 * (q < 1 // 2) * q1_2^3 # Zero out result if q >= 5/2 - result = ifelse(q < 5 / 2, normalization_factor(kernel, h) * result * inner_deriv, + result = ifelse(q < 5 // 2, normalization_factor(kernel, h) * result * inner_deriv, zero(result)) return result @@ -251,8 +260,9 @@ end @inline compact_support(::SchoenbergQuarticSplineKernel, h) = 2.5 * h @inline normalization_factor(::SchoenbergQuarticSplineKernel{1}, h) = 1 / 24h -@inline normalization_factor(::SchoenbergQuarticSplineKernel{2}, h) = 96 / (1199 * pi * h^2) -@inline normalization_factor(::SchoenbergQuarticSplineKernel{3}, h) = 1 / (20 * pi * h^3) +# `1199 * pi` is always `Float64`. `pi * h^2 * 1199` preserves the type of `h`. +@inline normalization_factor(::SchoenbergQuarticSplineKernel{2}, h) = 96 / (pi * h^2 * 1199) +@inline normalization_factor(::SchoenbergQuarticSplineKernel{3}, h) = 1 / (pi * h^3 * 20) @doc raw""" SchoenbergQuinticSplineKernel{NDIMS}() @@ -329,8 +339,9 @@ end @inline compact_support(::SchoenbergQuinticSplineKernel, h) = 3 * h @inline normalization_factor(::SchoenbergQuinticSplineKernel{1}, h) = 1 / 120h -@inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (478 * pi * h^2) -@inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (120 * pi * h^3) +# `478 * pi` is always `Float64`. `pi * h^2 * 478` preserves the type of `h`. +@inline normalization_factor(::SchoenbergQuinticSplineKernel{2}, h) = 7 / (pi * h^2 * 478) +@inline normalization_factor(::SchoenbergQuinticSplineKernel{3}, h) = 1 / (pi * h^3 * 120) abstract type WendlandKernel{NDIMS} <: SmoothingKernel{NDIMS} end @@ -400,7 +411,8 @@ end end @inline normalization_factor(::WendlandC2Kernel{2}, h) = 7 / (pi * h^2) -@inline normalization_factor(::WendlandC2Kernel{3}, h) = 21 / (2pi * h^3) +# `2 * pi` is always `Float64`. `pi * h^3 * 2` preserves the type of `h`. +@inline normalization_factor(::WendlandC2Kernel{3}, h) = 21 / (pi * h^3 * 2) @doc raw""" WendlandC4Kernel{NDIMS}() @@ -449,8 +461,10 @@ end @fastpow @muladd @inline function kernel_deriv(kernel::WendlandC4Kernel, r::Real, h) q = r / h - term1 = (1 - q)^6 * (6 + 70 / 3 * q) - term2 = 6 * (1 - q)^5 * (1 + 6q + 35 / 3 * q^2) + + # Use `//` to preserve the type of `q` + term1 = (1 - q)^6 * (6 + 70 // 3 * q) + term2 = 6 * (1 - q)^5 * (1 + 6q + 35 // 3 * q^2) derivative = term1 - term2 # Zero out result if q >= 1 @@ -461,7 +475,8 @@ end end @inline normalization_factor(::WendlandC4Kernel{2}, h) = 9 / (pi * h^2) -@inline normalization_factor(::WendlandC4Kernel{3}, h) = 495 / (32pi * h^3) +# `32 * pi` is always `Float64`. `pi * h^2 * 32` preserves the type of `h`. +@inline normalization_factor(::WendlandC4Kernel{3}, h) = 495 / (pi * h^3 * 32) @doc raw""" WendlandC6Kernel{NDIMS}() @@ -521,8 +536,9 @@ end return result end -@inline normalization_factor(::WendlandC6Kernel{2}, h) = 78 / (7pi * h^2) -@inline normalization_factor(::WendlandC6Kernel{3}, h) = 1365 / (64pi * h^3) +# `7 * pi` is always `Float64`. `pi * h^2 * 7` preserves the type of `h`. +@inline normalization_factor(::WendlandC6Kernel{2}, h) = 78 / (pi * h^2 * 7) +@inline normalization_factor(::WendlandC6Kernel{3}, h) = 1365 / (pi * h^3 * 64) @doc raw""" Poly6Kernel{NDIMS}() @@ -588,7 +604,8 @@ end @inline compact_support(::Poly6Kernel, h) = h @inline normalization_factor(::Poly6Kernel{2}, h) = 4 / (pi * h^2) -@inline normalization_factor(::Poly6Kernel{3}, h) = 315 / (64pi * h^3) +# `64 * pi` is always `Float64`. `pi * h^3 * 64` preserves the type of `h`. +@inline normalization_factor(::Poly6Kernel{3}, h) = 315 / (pi * h^3 * 64) @doc raw""" SpikyKernel{NDIMS}() diff --git a/src/general/system.jl b/src/general/system.jl index 674555aea..88b9e7bf2 100644 --- a/src/general/system.jl +++ b/src/general/system.jl @@ -54,28 +54,31 @@ initialize!(system, neighborhood_search) = system @inline active_particles(system, ::Nothing) = eachparticle(system) # This should not be dispatched by system type. We always expect to get a column of `A`. -@inline function extract_svector(A, system, i) +@propagate_inbounds function extract_svector(A, system, i) extract_svector(A, Val(ndims(system)), i) end # Return the `i`-th column of the array `A` as an `SVector`. @inline function extract_svector(A, ::Val{NDIMS}, i) where {NDIMS} - return SVector(ntuple(@inline(dim->A[dim, i]), NDIMS)) + # Explicit bounds check, which can be removed by calling this function with `@inbounds` + @boundscheck checkbounds(A, NDIMS, i) + + # Assume inbounds access now + return SVector(ntuple(@inline(dim->@inbounds A[dim, i]), NDIMS)) end # Return `A[:, :, i]` as an `SMatrix`. @inline function extract_smatrix(A, system, particle) # Extract the matrix elements for this particle as a tuple to pass to SMatrix - return SMatrix{ndims(system), ndims(system)}( - # Convert linear index to Cartesian index - ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, - div(i - 1, ndims(system)) + 1, - particle]), - Val(ndims(system)^2))) + return SMatrix{ndims(system), + ndims(system)}(ntuple(@inline(i->A[mod(i - 1, ndims(system)) + 1, + div(i - 1, ndims(system)) + 1, + particle]), + Val(ndims(system)^2))) end # Specifically get the current coordinates of a particle for all system types. -@inline function current_coords(u, system, particle) +@propagate_inbounds function current_coords(u, system, particle) return extract_svector(current_coordinates(u, system), system, particle) end @@ -91,7 +94,8 @@ end # This can be dispatched by system type. @inline initial_coordinates(system) = system.initial_condition.coordinates -@inline current_velocity(v, system, particle) = extract_svector(v, system, particle) +@propagate_inbounds current_velocity(v, system, particle) = extract_svector(v, system, + particle) @inline function current_acceleration(system, particle) # TODO: Return `dv` of solid particles @@ -102,30 +106,23 @@ end @inline set_particle_pressure!(v, system, particle, pressure) = v -@inline function smoothing_kernel(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel(smoothing_kernel, distance, smoothing_length) -end - -@inline function smoothing_kernel_deriv(system, distance) - (; smoothing_kernel, smoothing_length) = system - return kernel_deriv(smoothing_kernel, distance, smoothing_length) +@inline function smoothing_kernel(system, distance, particle) + (; smoothing_kernel) = system + return kernel(smoothing_kernel, distance, smoothing_length(system, particle)) end -@inline function smoothing_kernel_grad(system, pos_diff, distance) - return kernel_grad(system.smoothing_kernel, pos_diff, distance, system.smoothing_length) -end - -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance) +@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) (; smoothing_kernel, smoothing_length) = system.boundary_model return kernel_grad(smoothing_kernel, pos_diff, distance, smoothing_length) end @inline function smoothing_kernel_grad(system, pos_diff, distance, particle) - return corrected_kernel_grad(system.smoothing_kernel, pos_diff, distance, - system.smoothing_length, system.correction, system, - particle) + (; smoothing_kernel) = system + + return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, + smoothing_length(system, particle), + system.correction, system, particle) end # System update orders. This can be dispatched if needed. diff --git a/src/refinement/refinement.jl b/src/refinement/refinement.jl new file mode 100644 index 000000000..b76021d21 --- /dev/null +++ b/src/refinement/refinement.jl @@ -0,0 +1,142 @@ +include("refinement_criteria.jl") +include("refinement_pattern.jl") +include("split.jl") + +struct ParticleRefinement{SP, RC, ELTYPE} + refinement_pattern :: SP + refinement_criteria :: RC + max_spacing_ratio :: ELTYPE + mass_ref :: Vector{ELTYPE} # length(mass_ref) == nparticles + split_candidates :: Vector{Int} + n_particles_before_resize :: Ref{Int} + n_new_particles :: Ref{Int} +end + +function ParticleRefinement(; refinement_pattern, max_spacing_ratio, + refinement_criteria=SpatialRefinementCriterion()) + mass_ref = Vector{eltype(max_spacing_ratio)}() + + if refinement_criteria isa Tuple + refinement_criteria = (refinement_criteria,) + end + + return ParticleRefinement(refinement_pattern, refinement_criteria, max_spacing_ratio, + mass_ref, Int[], Ref(0), Ref(0)) +end + +resize_refinement!(system) = system + +function resize_refinement!(system::FluidSystem) + resize_refinement!(system, system.particle_refinement) +end + +resize_refinement!(system, ::Nothing) = system + +function resize_refinement!(system, particle_refinement) + resize!(particle_refinement.mass_ref, nparticles(system)) + + return system +end + +function refinement!(semi, v_ode, u_ode, v_tmp, u_tmp, t) + check_refinement_criteria!(semi, v_ode, u_ode) + + # Update the spacing of particles (Algorthm 1) + update_particle_spacing(semi, u_ode) + + # Split the particles (Algorithm 2) + split_particles!(semi, v_ode, u_ode, v_tmp, u_tmp) + + # Merge the particles (Algorithm 3) + + # Shift the particles + + # Correct the particles + + # Update smoothing lengths + + # Resize neighborhood search + + return semi +end + +function update_particle_spacing(semi, u_ode) + foreach_system(semi) do system + u = wrap_u(u_ode, system, semi) + update_particle_spacing(system, u, semi) + end +end + +# The methods for the `FluidSystem` are defined in `src/schemes/fluid/fluid.jl` +@inline update_particle_spacing(system, u, semi) = systemerror + +@inline function update_particle_spacing(system::FluidSystem, u, semi) + update_particle_spacing(system, system.particle_refinement, u, semi) +end + +@inline update_particle_spacing(system, ::Nothing, u, semi) = system + +@inline function update_particle_spacing(system::FluidSystem, particle_refinement, u, semi) + (; smoothing_length, smoothing_length_factor) = system.cache + (; mass_ref) = particle_refinement + + system_coords = current_coordinates(u, system) + + for particle in eachparticle(system) + dp_min, dp_max, dp_avg = min_max_avg_spacing(system, semi, system_coords, particle) + + if dp_max / dp_min < max_spacing_ratio^3 + new_spacing = min(dp_max, max_spacing_ratio * dp_min) + else + new_spacing = dp_avg + end + + smoothing_length[particle] = smoothing_length_factor * new_spacing + mass_ref[particle] = system.density[particle] * new_spacing^(ndims(system)) + end + + return system +end + +@inline function min_max_avg_spacing(system, semi, system_coords, particle) + dp_min = Inf + dp_max = zero(eltype(system)) + dp_avg = zero(eltype(system)) + counter_neighbors = 0 + + foreach_system(semi) do neighbor_system + neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, + semi) + + u_neighbor_system = wrap_u(u_ode, neighbor, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + PointNeighbors.foreach_neighbor(system_coords, neighbor_coords, neighborhood_search, + particle) do particle, neighbor, pos_diff, distance + dp_neighbor = particle_spacing(neighbor_system, neighbor) + + dp_min = min(dp_min, dp_neighbor) + dp_max = max(dp_max, dp_neighbor) + dp_avg += dp_neighbor + + counter_neighbors += 1 + end + end + + dp_avg / counter_neighbors + + return dp_min, dp_max, dp_avg +end + +@inline particle_spacing(system, particle) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system::FluidSystem, particle) + return particle_spacing(system, system.particle_refinement, particle) +end + +@inline particle_spacing(system, ::Nothing, _) = system.initial_condition.particle_spacing + +@inline function particle_spacing(system, refinement, particle) + (; smoothing_length_factor) = system.cache + return smoothing_length(system, particle) / smoothing_length_factor +end diff --git a/src/refinement/refinement_criteria.jl b/src/refinement/refinement_criteria.jl new file mode 100644 index 000000000..de4bab43a --- /dev/null +++ b/src/refinement/refinement_criteria.jl @@ -0,0 +1,64 @@ +abstract type RefinementCriteria end + +struct SpatialRefinementCriterion <: RefinementCriteria end + +struct SolutionRefinementCriterion <: RefinementCriteria end + +function check_refinement_criteria!(semi, v_ode, u_ode) + foreach_system(semi) do system + check_refinement_criteria!(system, v_ode, u_ode, semi) + end +end + +@inline check_refinement_criteria!(system, v_ode, u_ode, semi) = system + +@inline function check_refinement_criteria!(system::FluidSystem, v_ode, u_ode, semi) + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + + check_refinement_criteria!(system, system.particle_refinement, v, u, v_ode, u_ode, semi) +end + +@inline function check_refinement_criteria!(system::FluidSystem, refinement, + v, u, v_ode, u_ode, semi) + (; refinement_criteria) = system.particle_refinement + for criterion in refinement_criteria + criterion(system, v, u, v_ode, u_ode, semi) + end +end + +@inline function (criterion::SpatialRefinementCriterion)(system, v, u, v_ode, u_ode, semi) + system_coords = current_coordinates(u, system) + + foreach_system(semi) do neighbor_system + set_particle_spacing!(system, neighbor_system, system_coords, v_ode, u_ode, semi) + end + return system +end + +@inline set_particle_spacing!(system, _, _, _, _) = system + +@inline function set_particle_spacing!(particle_system, + neighbor_system::Union{BoundarySystem, SolidSystem}, + system_coords, v_ode, u_ode, semi) + (; smoothing_length, smoothing_length_factor) = particle_system.cache + + neighborhood_search = get_neighborhood_search(particle_system, neighbor_system, semi) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) + + # Loop over all pairs of particles and neighbors within the kernel cutoff. + foreach_point_neighbor(particle_system, neighbor_system, + system_coords, neighbor_coords, + neighborhood_search) do particle, neighbor, pos_diff, distance + # Only consider particles with a distance > 0. + distance < sqrt(eps()) && return + + dp_neighbor = particle_spacing(neighbor_system, neighbor) + dp_particle = smoothing_length[particle] / smoothing_length_factor + + smoothing_length[particle] = smoothing_length_factor * min(dp_neighbor, dp_particle) + end + + return particle_system +end diff --git a/src/refinement/refinement_pattern.jl b/src/refinement/refinement_pattern.jl new file mode 100644 index 000000000..759fe8262 --- /dev/null +++ b/src/refinement/refinement_pattern.jl @@ -0,0 +1,93 @@ +struct CubicSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function CubicSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1 / sqrt(2), 1 / sqrt(2)] + direction_2 = [1 / sqrt(2), -1 / sqrt(2)] + direction_3 = -direction_1 + direction_4 = -direction_2 + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system, refinement_pattern::CubicSplitting) + return 4 + refinement_pattern.center_particle +end + +struct TriangularSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function TriangularSplitting(; epsilon=0.5, alpha=0.5, center_particle=true) + ELTYPE = typeof(epsilon) + direction_1 = [0.0, 1.0] + direction_2 = [-sqrt(3) / 2, -0.5] + direction_3 = [sqrt(3) / 2, -0.5] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::TriangularSplitting) + return 3 + refinement_pattern.center_particle +end + +struct HexagonalSplitting{ELTYPE} + epsilon :: ELTYPE + alpha :: ELTYPE + center_particle :: Bool + relative_position :: Vector{SVector{2, ELTYPE}} + + function HexagonalSplitting(; epsilon=0.4, alpha=0.9, center_particle=true) + ELTYPE = typeof(epsilon) + + direction_1 = [1.0, 0.0] + direction_2 = [-1.0, 0.0] + direction_3 = [0.5, sqrt(3) / 2] + direction_4 = [0.5, -sqrt(3) / 2] + direction_5 = [-0.5, sqrt(3) / 2] + direction_6 = [-0.5, -sqrt(3) / 2] + + relative_position = [SVector{2, ELTYPE}(epsilon * direction_1), + SVector{2, ELTYPE}(epsilon * direction_2), + SVector{2, ELTYPE}(epsilon * direction_3), + SVector{2, ELTYPE}(epsilon * direction_4), + SVector{2, ELTYPE}(epsilon * direction_5), + SVector{2, ELTYPE}(epsilon * direction_6)] + + if center_particle + push!(relative_position, SVector(zero(ELTYPE), zero(ELTYPE))) + end + + new{ELTYPE}(epsilon, alpha, center_particle, relative_position) + end +end + +@inline function nchilds(system::System{2}, refinement_pattern::HexagonalSplitting) + return 6 + refinement_pattern.center_particle +end diff --git a/src/refinement/resize.jl b/src/refinement/resize.jl new file mode 100644 index 000000000..a7e70fa47 --- /dev/null +++ b/src/refinement/resize.jl @@ -0,0 +1,175 @@ +function Base.resize!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Resize all systems + foreach_system(semi) do system + resize!(system, capacity(system)) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function Base.deleteat!(semi::Semidiscretization, v_ode, u_ode, _v_ode, _u_ode) + # Delete at specific indices + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + deleteat!(system, v, u) + end + + resize!(v_ode, u_ode, _v_ode, _u_ode, semi) + + return semi +end + +function Base.resize!(v_ode, u_ode, _v_ode, _u_ode, semi::Semidiscretization) + copyto!(_v_ode, v_ode) + copyto!(_u_ode, u_ode) + + # Get ranges after resizing the systems + ranges_v_new, ranges_u_new = ranges_vu(semi.systems) + + ranges_v_old = copy(semi.ranges_v) + ranges_u_old = copy(semi.ranges_u) + + # Set ranges after resizing the systems + for i in 1:length(semi.systems) + semi.ranges_v[i] = ranges_v_new[i] + semi.ranges_u[i] = ranges_u_new[i] + end + + for i in eachindex(ranges_u_old) + length_u = min(length(ranges_u_old[i]), length(ranges_u_new[i])) + for j in 0:(length_u - 1) + u_ode[ranges_u_new[i][1] + j] = _u_ode[ranges_u_old[i][1] + j] + end + + length_v = min(length(ranges_v_old[i]), length(ranges_v_new[i])) + for j in 0:(length_v - 1) + v_ode[ranges_v_new[i][1] + j] = _v_ode[ranges_v_old[i][1] + j] + end + end + + capacity_global = sum(system -> nparticles(system), semi.systems) + + resize!(v_ode, capacity_global) + resize!(u_ode, capacity_global) + + resize!(_v_ode, capacity_global) + resize!(_u_ode, capacity_global) + + # TODO: Do the following in the callback + # resize!(integrator, (length(v_ode), length(u_ode))) + + # # Tell OrdinaryDiffEq that u has been modified + # u_modified!(integrator, true) + return v_ode +end + +Base.resize!(system::System, capacity_system) = system + +function Base.resize!(system::FluidSystem, capacity_system) + return resize!(system, system.particle_refinement, capacity_system) +end + +Base.resize!(system, ::Nothing, capacity_system) = system + +function Base.resize!(system::WeaklyCompressibleSPHSystem, refinement, capacity_system::Int) + (; mass, pressure, cache, density_calculator) = system + + refinement.n_particles_before_resize[] = nparticles(system) + + resize!(mass, capacity_system) + resize!(pressure, capacity_system) + resize_density!(system, capacity_system, density_calculator) + # TODO + # resize_cache!(system, cache, n) +end + +function Base.resize!(system::EntropicallyDampedSPHSystem, refinement, capacity_system::Int) + (; mass, cache, density_calculator) = system + + refinement.n_particles_before_resize[] = nparticles(system) + + resize!(mass, capacity_system) + resize_density!(system, capacity_system, density_calculator) + # TODO + # resize_cache!(system, capacity_system) + + return system +end + +resize_density!(system, n::Int, ::SummationDensity) = resize!(system.cache.density, n) +resize_density!(system, n::Int, ::ContinuityDensity) = system + +function resize_cache!(system, n::Int) + resize!(system.cache.smoothing_length, n) + + return system +end + +function resize_cache!(system::EntropicallyDampedSPHSystem, n) + resize!(system.cache.smoothing_length, n) + resize!(system.cache.pressure_average, n) + resize!(system.cache.neighbor_counter, n) + + return system +end + +Base.deleteat!(system::System, v, u) = system + +function Base.deleteat!(system::FluidSystem, v, u) + return deleteat!(system, system.particle_refinement, v, u) +end + +Base.deleteat!(system::FluidSystem, ::Nothing, v, u) = system + +function Base.deleteat!(system::FluidSystem, refinement, v, u) + (; delete_candidates) = refinement + + delete_counter = 0 + + for particle in eachparticle(system) + if !iszero(delete_candidates[particle]) + # swap particles (keep -> delete) + dump_id = nparticles(system) - delete_counter + + vel_keep = current_velocity(v, system, dump_id) + pos_keep = current_coords(u, system, dump_id) + + mass_keep = hydrodynamic_mass(system, dump_id) + density_keep = particle_density(v, system, dump_id) + pressure_keep = particle_pressure(v, system, dump_id) + #TODO + # smoothing_length_keep = smoothing_length(system, dump_id) + # system.cache.smoothing_length[particle] = smoothing_length_keep + + system.mass[particle] = mass_keep + + set_particle_pressure!(v, system, particle, pressure_keep) + + set_particle_density!(v, system, particle, density_keep) + + for dim in 1:ndims(system) + v[dim, particle] = vel_keep[dim] + u[dim, particle] = pos_keep[dim] + end + + delete_counter += 1 + end + end + + resize!(system, nparticles(system) - delete_counter) + + return system +end + +@inline capacity(system) = nparticles(system) + +@inline capacity(system::FluidSystem) = capacity(system, system.particle_refinement) + +@inline capacity(system, ::Nothing) = nparticles(system) + +@inline function capacity(system, particle_refinement) + return particle_refinement.n_new_particles[] + nparticles(system) +end diff --git a/src/refinement/split.jl b/src/refinement/split.jl new file mode 100644 index 000000000..47c83fef1 --- /dev/null +++ b/src/refinement/split.jl @@ -0,0 +1,98 @@ +function split_particles!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + collect_split_candidates!(system, v, u) + end + + resize!(semi, v_ode, u_ode, v_ode_tmp, u_ode_tmp) + + foreach_system(semi) do system + v = wrap_v(v_ode, system, semi) + u = wrap_u(u_ode, system, semi) + split_particles!(system, v, u) + end + + return semi +end + +@inline collect_split_candidates!(system, v, u) = System + +@inline function collect_split_candidates!(system::FluidSystem, v, u) + return collect_split_candidates!(system, system.particle_refinement, v, u) +end + +@inline collect_split_candidates!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function collect_split_candidates!(system::FluidSystem, particle_refinement, v, u) + (; mass_ref, max_spacing_ratio, split_candidates, refinement_pattern) = particle_refinement + + resize!(split_candidates, 0) + + for particle in eachparticle(system) + m_a = hydrodynamic_mass(system, particle) + m_max = max_spacing_ratio * mass_ref[particle] + + if m_a > m_max + push!(split_candidates, particle) + end + end + + n_childs_exclude_center = nchilds(system, refinement_pattern) - 1 + + particle_refinement.n_new_particles[] = length(split_candidates) * + n_childs_exclude_center + + return system +end + +@inline split_particles!(system, v, u) = System + +@inline function split_particles!(system::FluidSystem, v, u) + return split_particles!(system, system.particle_refinement, v, u) +end + +@inline split_particles!(system::FluidSystem, ::Nothing, v, u) = system + +@inline function split_particles!(system::FluidSystem, particle_refinement, v, u) + (; smoothing_length) = system.cache + (; split_candidates, refinement_pattern, n_particles_before_resize) = particle_refinement + (; alpha, relative_position) = refinement_pattern + + for particle in split_candidates + smoothing_length_old = smoothing_length[particle] + mass_old = system.mass[particle] + + system.mass[particle] = mass_old / nchilds(system, refinement_pattern) + + p_a = particle_pressure(v, system, particle) + rho_a = particle_density(v, system, particle) + + smoothing_length[particle] = alpha * smoothing_length_old + + pos_center = current_coords(u, system, particle) + vel_center = current_velocity(v, system, particle) + + for child_id_local in 1:(nchilds(system, refinement_pattern) - 1) + child = n_particles_before_resize[] + child_id_local + + system.mass[child] = mass_old / nchilds(system, refinement_pattern) + + set_particle_pressure!(v, system, child, p_a) + + set_particle_density!(v, system, child, rho_a) + + smoothing_length[child] = alpha * smoothing_length_old + + rel_pos = smoothing_length_old * relative_position[child_id_local] + new_pos = pos_center + rel_pos + + for dim in 1:ndims(system) + u[dim, child] = new_pos[dim] + v[dim, child] = vel_center[dim] + end + end + end + + return system +end diff --git a/src/schemes/boundary/dummy_particles/dummy_particles.jl b/src/schemes/boundary/dummy_particles/dummy_particles.jl index 1cd8e0dfe..ece912eeb 100644 --- a/src/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/src/schemes/boundary/dummy_particles/dummy_particles.jl @@ -69,12 +69,48 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, smoothing_length, viscosity, correction, cache) end +function smoothing_length(boundary_model::BoundaryModelDummyParticles, particle) + return boundary_model.smoothing_length +end + @doc raw""" - AdamiPressureExtrapolation() + AdamiPressureExtrapolation(; pressure_offset=0.0) `density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. + """ -struct AdamiPressureExtrapolation end +struct AdamiPressureExtrapolation{ELTYPE} + pressure_offset::ELTYPE + + function AdamiPressureExtrapolation(; pressure_offset=0.0) + return new{eltype(pressure_offset)}(pressure_offset) + end +end + +@doc raw""" + BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=1.0) + +`density_calculator` for `BoundaryModelDummyParticles`. + +# Keywords +- `pressure_offset=0.0`: Sometimes it is necessary to artificially increase the boundary pressure + to prevent penetration, which is possible by increasing this value. +- `factor=1.0` : Setting `factor` allows to just increase the strength of the dynamic + pressure part. + +""" +struct BernoulliPressureExtrapolation{ELTYPE} + pressure_offset :: ELTYPE + factor :: ELTYPE + + function BernoulliPressureExtrapolation(; pressure_offset=0.0, factor=0.0) + return new{eltype(pressure_offset)}(pressure_offset, factor) + end +end @doc raw""" PressureMirroring() @@ -134,7 +170,9 @@ end @inline create_cache_model(initial_density, ::ContinuityDensity) = (; initial_density) -function create_cache_model(initial_density, ::AdamiPressureExtrapolation) +function create_cache_model(initial_density, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}) density = copy(initial_density) volume = similar(initial_density) @@ -194,7 +232,7 @@ end # Note that the other density calculators are dispatched in `density_calculators.jl` @inline function particle_density(v, ::Union{AdamiPressureExtrapolation, PressureMirroring, - PressureZeroing}, + PressureZeroing, BernoulliPressureExtrapolation}, boundary_model, particle) (; cache) = boundary_model @@ -216,10 +254,11 @@ end function compute_density!(boundary_model, ::Union{ContinuityDensity, AdamiPressureExtrapolation, + BernoulliPressureExtrapolation, PressureMirroring, PressureZeroing}, system, v, u, v_ode, u_ode, semi) # No density update for `ContinuityDensity`, `PressureMirroring` and `PressureZeroing`. - # For `AdamiPressureExtrapolation`, the density is updated in `compute_pressure!`. + # For `AdamiPressureExtrapolation` and `BernoulliPressureExtrapolation`, the density is updated in `compute_pressure!`. return boundary_model end @@ -299,8 +338,10 @@ end boundary_model.pressure[particle] = max(boundary_model.state_equation(density), 0.0) end -function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, v, u, - v_ode, u_ode, semi) +function compute_pressure!(boundary_model, + ::Union{AdamiPressureExtrapolation, + BernoulliPressureExtrapolation}, + system, v, u, v_ode, u_ode, semi) (; pressure, cache, viscosity) = boundary_model set_zero!(pressure) @@ -327,16 +368,18 @@ function compute_pressure!(boundary_model, ::AdamiPressureExtrapolation, system, nhs = get_neighborhood_search(neighbor_system, system, semi) # Loop over fluid particles and then the neighboring boundary particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation_neighbor!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) else nhs = get_neighborhood_search(system, neighbor_system, semi) # Loop over boundary particles and then the neighboring fluid particles to extrapolate fluid pressure to the boundaries - adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, nhs) + boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, nhs) end @threaded system for particle in eachparticle(system) @@ -378,67 +421,81 @@ function compute_pressure!(boundary_model, ::Union{PressureMirroring, PressureZe return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system, system_coords, - neighbor_coords, v_neighbor_system, - neighborhood_search) +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system, system_coords, + neighbor_coords, v, + v_neighbor_system, + neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation_neighbor!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, - neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation_neighbor!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, + v, v_neighbor_system, + neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator foreach_point_neighbor(neighbor_system, system, neighbor_coords, system_coords, neighborhood_search; points=eachparticle(neighbor_system), parallel=false) do neighbor, particle, pos_diff, distance # Since neighbor and particle are switched pos_diff = -pos_diff - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_extrapolation!(boundary_model, system, neighbor_system, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) +@inline function boundary_pressure_extrapolation!(boundary_model, system, neighbor_system, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) return boundary_model end -@inline function adami_pressure_extrapolation!(boundary_model, system, - neighbor_system::FluidSystem, - system_coords, neighbor_coords, - v_neighbor_system, neighborhood_search) - (; pressure, cache, viscosity) = boundary_model +@inline function boundary_pressure_extrapolation!(boundary_model, system, + neighbor_system::FluidSystem, + system_coords, neighbor_coords, v, + v_neighbor_system, neighborhood_search) + (; pressure, cache, viscosity, density_calculator) = boundary_model + (; pressure_offset) = density_calculator # Loop over all pairs of particles and neighbors within the kernel cutoff. foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, neighborhood_search; points=eachparticle(system)) do particle, neighbor, pos_diff, distance - adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) + boundary_pressure_inner!(boundary_model, density_calculator, system, + neighbor_system, v, v_neighbor_system, particle, neighbor, + pos_diff, distance, viscosity, cache, pressure, + pressure_offset) end end -@inline function adami_pressure_inner!(boundary_model, system, neighbor_system, - v_neighbor_system, particle, neighbor, pos_diff, - distance, viscosity, cache, pressure) +@inline function boundary_pressure_inner!(boundary_model, boundary_density_calculator, + system, neighbor_system::FluidSystem, v, + v_neighbor_system, particle, neighbor, pos_diff, + distance, viscosity, cache, pressure, + pressure_offset) density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) - resulting_acc = neighbor_system.acceleration - - current_acceleration(system, particle) - - kernel_weight = smoothing_kernel(boundary_model, distance) - - pressure[particle] += (particle_pressure(v_neighbor_system, neighbor_system, - neighbor) + - dot(resulting_acc, density_neighbor * pos_diff)) * + resulting_acceleration = neighbor_system.acceleration - + current_acceleration(system, particle) + + kernel_weight = smoothing_kernel(boundary_model, distance, particle) + + pressure[particle] += (pressure_offset + + + particle_pressure(v_neighbor_system, neighbor_system, + neighbor) + + + dynamic_pressure(boundary_density_calculator, density_neighbor, + v, v_neighbor_system, pos_diff, distance, + particle, neighbor, system, neighbor_system) + + + dot(resulting_acceleration, density_neighbor * pos_diff)) * kernel_weight cache.volume[particle] += kernel_weight @@ -447,14 +504,46 @@ end kernel_weight, particle, neighbor) end +@inline function dynamic_pressure(boundary_density_calculator, density_neighbor, v, + v_neighbor_system, pos_diff, distance, particle, neighbor, + system, neighbor_system) + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::BoundarySystem, neighbor_system) + if system.ismoving[] + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) + + return 0.5 * boundary_density_calculator.factor * density_neighbor * + normal_velocity^2 / distance + end + return zero(density_neighbor) +end + +@inline function dynamic_pressure(boundary_density_calculator::BernoulliPressureExtrapolation, + density_neighbor, v, v_neighbor_system, pos_diff, + distance, particle, neighbor, + system::SolidSystem, neighbor_system) + relative_velocity = current_velocity(v, system, particle) .- + current_velocity(v_neighbor_system, neighbor_system, neighbor) + normal_velocity = dot(relative_velocity, pos_diff) / distance + + return boundary_density_calculator.factor * density_neighbor * + dot(normal_velocity, normal_velocity) / 2 +end + function compute_smoothed_velocity!(cache, viscosity, neighbor_system, v_neighbor_system, kernel_weight, particle, neighbor) return cache end -function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, - neighbor_system, v_neighbor_system, kernel_weight, - particle, neighbor) +function compute_smoothed_velocity!(cache, viscosity::ViscosityAdami, neighbor_system, + v_neighbor_system, kernel_weight, particle, neighbor) v_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) for dim in 1:ndims(neighbor_system) @@ -496,12 +585,13 @@ end return density end -@inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) - (; smoothing_kernel, smoothing_length, correction) = system.boundary_model +# TODO +# @inline function smoothing_kernel_grad(system::BoundarySystem, pos_diff, distance, particle) +# (; smoothing_kernel, smoothing_length, correction) = system.boundary_model - return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, - smoothing_length, correction, system, particle) -end +# return corrected_kernel_grad(smoothing_kernel, pos_diff, distance, +# smoothing_length, correction, system, particle) +# end @inline function correction_matrix(system::BoundarySystem, particle) extract_smatrix(system.boundary_model.cache.correction_matrix, system, particle) diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index cc9a77157..ad72045d6 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -116,11 +116,11 @@ struct InFlow{NDIMS, IC, S, ZO, ZW, FD} if !isapprox(abs(dot_), 1.0, atol=1e-7) throw(ArgumentError("`flow_direction` is not normal to inflow plane")) - else - # Flip the normal vector to point in the opposite direction of `flow_direction` - spanning_set[:, 1] .*= -sign(dot_) end + # Flip the normal vector to point in the opposite direction of `flow_direction` + spanning_set[:, 1] .*= -sign(dot_) + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) # Remove particles outside the boundary zone. @@ -251,11 +251,11 @@ struct OutFlow{NDIMS, IC, S, ZO, ZW, FD} if !isapprox(abs(dot_), 1.0, atol=1e-7) throw(ArgumentError("`flow_direction` is not normal to outflow plane")) - else - # Flip the normal vector to point in `flow_direction` - spanning_set[:, 1] .*= sign(dot_) end + # Flip the normal vector to point in `flow_direction` + spanning_set[:, 1] .*= sign(dot_) + spanning_set_ = reinterpret(reshape, SVector{NDIMS, ELTYPE}, spanning_set) # Remove particles outside the boundary zone. @@ -289,8 +289,9 @@ function spanning_vectors(plane_points::NTuple{3}, zone_width) edge1 = plane_points[2] - plane_points[1] edge2 = plane_points[3] - plane_points[1] - if !isapprox(dot(edge1, edge2), 0.0, atol=1e-7) - throw(ArgumentError("the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal")) + # Check if the edges are linearly dependent (to avoid degenerate planes) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("the vectors `AB` and `AC` must not be collinear")) end # Calculate normal vector of plane diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index 38fdd477a..f803bf9ca 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -8,8 +8,9 @@ about the method see [description below](@ref method_of_characteristics). """ struct BoundaryModelLastiwka end -@inline function update_quantities!(system, boundary_model::BoundaryModelLastiwka, - v, u, v_ode, u_ode, semi, t) +# Called from update callback via `update_open_boundary_eachstep!` +@inline function update_boundary_quantities!(system, boundary_model::BoundaryModelLastiwka, + v, u, v_ode, u_ode, semi, t) (; density, pressure, cache, flow_direction, reference_velocity, reference_pressure, reference_density) = system @@ -45,10 +46,13 @@ struct BoundaryModelLastiwka end return system end +# Called from semidiscretization function update_final!(system, ::BoundaryModelLastiwka, v, u, v_ode, u_ode, semi, t) @trixi_timeit timer() "evaluate characteristics" begin evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end + + return system end # ==== Characteristics @@ -98,9 +102,16 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) end end - characteristics[1, particle] = avg_J1 / counter - characteristics[2, particle] = avg_J2 / counter - characteristics[3, particle] = avg_J3 / counter + # To prevent NANs here if the boundary has not been in contact before. + if counter > 0 + characteristics[1, particle] = avg_J1 / counter + characteristics[2, particle] = avg_J2 / counter + characteristics[3, particle] = avg_J3 / counter + else + characteristics[1, particle] = 0 + characteristics[2, particle] = 0 + characteristics[3, particle] = 0 + end else characteristics[1, particle] /= volume[particle] characteristics[2, particle] /= volume[particle] @@ -164,7 +175,7 @@ end @inline function prescribe_conditions!(characteristics, particle, ::OutFlow) # J3 is prescribed (i.e. determined from the exterior of the domain). - # J1 and J2 is transimtted from the domain interior. + # J1 and J2 is transmitted from the domain interior. characteristics[3, particle] = zero(eltype(characteristics)) return characteristics diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9d2c79b00..c8d80ca8f 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -77,6 +77,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "an array where the ``i``-th column holds the velocity of particle ``i`` " * "or, for a constant fluid velocity, a vector of length $NDIMS for a $(NDIMS)D problem holding this velocity")) else + if reference_velocity isa Function + test_result = reference_velocity(zeros(NDIMS), 0.0) + if length(test_result) != NDIMS + throw(ArgumentError("`reference_velocity` function must be of dimension $NDIMS")) + end + end reference_velocity_ = wrap_reference_function(reference_velocity, Val(NDIMS)) end @@ -86,6 +92,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its pressure, " * "a vector holding the pressure of each particle, or a scalar")) else + if reference_pressure isa Function + test_result = reference_pressure(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_pressure` function must be a scalar function")) + end + end reference_pressure_ = wrap_reference_function(reference_pressure, Val(NDIMS)) end @@ -95,6 +107,12 @@ struct OpenBoundarySPHSystem{BM, BZ, NDIMS, ELTYPE <: Real, IC, FS, ARRAY1D, RV, "each particle's coordinates and time to its density, " * "a vector holding the density of each particle, or a scalar")) else + if reference_density isa Function + test_result = reference_density(zeros(NDIMS), 0.0) + if length(test_result) != 1 + throw(ArgumentError("`reference_density` function must be a scalar function")) + end + end reference_density_ = wrap_reference_function(reference_density, Val(NDIMS)) end @@ -190,10 +208,12 @@ function update_open_boundary_eachstep!(system::OpenBoundarySPHSystem, v_ode, u_ # Update density, pressure and velocity based on the characteristic variables. # See eq. 13-15 in Lastiwka (2009) https://doi.org/10.1002/fld.1971 - @trixi_timeit timer() "update quantities" update_quantities!(system, - system.boundary_model, - v, u, v_ode, u_ode, - semi, t) + @trixi_timeit timer() "update boundary quantities" update_boundary_quantities!(system, + system.boundary_model, + v, u, + v_ode, + u_ode, + semi, t) @trixi_timeit timer() "check domain" check_domain!(system, v, u, v_ode, u_ode, semi) diff --git a/src/schemes/boundary/rhs.jl b/src/schemes/boundary/rhs.jl index 59f634092..d90d2a291 100644 --- a/src/schemes/boundary/rhs.jl +++ b/src/schemes/boundary/rhs.jl @@ -29,7 +29,7 @@ function interact!(dv, v_particle_system, u_particle_system, v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(boundary_model, pos_diff, distance, particle) continuity_equation!(dv, fluid_density_calculator, m_b, rho_a, rho_b, v_diff, grad_kernel, particle) diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 3644c3630..b075813f2 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -30,9 +30,10 @@ function interact!(dv, v_particle_system, u_particle_system, m_a = hydrodynamic_mass(particle_system, particle) m_b = hydrodynamic_mass(neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - dv_pressure = pressure_acceleration(particle_system, neighbor_system, neighbor, + dv_pressure = pressure_acceleration(particle_system, neighbor_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) @@ -44,20 +45,29 @@ function interact!(dv, v_particle_system, u_particle_system, # Add convection term when using `TransportVelocityAdami` dv_convection = momentum_convection(particle_system, neighbor_system, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) + dv_velocity_correction = velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_convection[i] + + dv_velocity_correction[i] end v_diff = current_velocity(v_particle_system, particle_system, particle) - current_velocity(v_neighbor_system, neighbor_system, neighbor) - pressure_evolution!(dv, particle_system, v_diff, grad_kernel, - particle, pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) + pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, v_diff, grad_kernel, + particle, neighbor, pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) transport_velocity!(dv, particle_system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) @@ -69,11 +79,47 @@ function interact!(dv, v_particle_system, u_particle_system, return dv end -@inline function pressure_evolution!(dv, particle_system, v_diff, grad_kernel, particle, - pos_diff, distance, sound_speed, m_a, m_b, - p_a, p_b, rho_a, rho_b) - (; smoothing_length) = particle_system +@inline function pressure_evolution!(dv, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + v_diff, grad_kernel, particle, neighbor, + pos_diff, distance, sound_speed, + m_a, m_b, p_a, p_b, rho_a, rho_b) + (; particle_refinement) = particle_system + + # This is basically the continuity equation times `sound_speed^2` + artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) + + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + dv[end, particle] += beta_inv_a * artificial_eos + + pressure_damping_term(particle_system, neighbor_system, + particle_refinement, particle, neighbor, + pos_diff, distance, beta_inv_a, m_a, m_b, + p_a, p_b, rho_b, rho_b, grad_kernel) + + pressure_reduction(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_inv_a, beta_inv_b) + + return dv +end +@inline beta_correction(system, particle) = one(eltype(system)) + +@inline function beta_correction(system::FluidSystem, particle) + beta_correction(system, system.particle_refinement, particle) +end + +@inline beta_correction(particle_system, ::Nothing, particle) = one(eltype(particle_system)) + +@inline function beta_correction(particle_system, refinement, particle) + return inv(particle_system.cache.beta[particle]) +end +function pressure_damping_term(particle_system, neighbor_system, ::Nothing, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel) volume_a = m_a / rho_a volume_b = m_b / rho_b volume_term = (volume_a^2 + volume_b^2) / m_a @@ -81,15 +127,13 @@ end # EDAC pressure evolution pressure_diff = p_a - p_b - # This is basically the continuity equation times `sound_speed^2` - artificial_eos = m_b * rho_a / rho_b * sound_speed^2 * dot(v_diff, grad_kernel) - eta_a = rho_a * particle_system.nu_edac eta_b = rho_b * particle_system.nu_edac eta_tilde = 2 * eta_a * eta_b / (eta_a + eta_b) - # TODO For variable smoothing length use average smoothing length - tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length(particle_system, particle) + + smoothing_length(particle_system, particle)) + tmp = eta_tilde / (distance^2 + 0.01 * smoothing_length_average^2) # This formulation was introduced by Hu and Adams (2006). https://doi.org/10.1016/j.jcp.2005.09.001 # They argued that the formulation is more flexible because of the possibility to formulate @@ -102,11 +146,99 @@ end # See issue: https://github.com/trixi-framework/TrixiParticles.jl/issues/394 # # This is similar to density diffusion in WCSPH - damping_term = volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) + return volume_term * tmp * pressure_diff * dot(grad_kernel, pos_diff) +end - dv[end, particle] += artificial_eos + damping_term +function pressure_damping_term(particle_system, neighbor_system, refinement, + particle, neighbor, pos_diff, distance, beta_inv_a, + m_a, m_b, p_a, p_b, rho_a, rho_b, grad_kernel_a) + # EDAC pressure evolution + pressure_diff = p_a - p_b - return dv + # Haftu et al. (2022) uses `alpha_edac = 1.5` in all their simulations + alpha_edac = 1.5 + + # TODO: Haftu et al. (2022) use `8` but I think it depeneds on the dimension (see Monaghan, 2005) + tmp = 2 * ndims(particle_system) + 4 + + nu_edac_a = alpha_edac * sound_speed * smoothing_length(particle_system, particle) / tmp + nu_edac_a = alpha_edac * sound_speed * smoothing_length(neighbor_system, particle) / tmp + + nu_edac_ab = 4 * (nu_edac_a * nu_edac_b) / (nu_edac_a + nu_edac_b) + + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + return beta_inv_a * nu_edac_ab * pressure_diff * dot(pos_diff, grad_W_avg) * m_b / rho_b +end + +function pressure_reduction(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + return zero(eltype(particle_system)) +end + +function pressure_reduction(particle_system, neighbor_system, refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, m_b, + p_a, p_b, rho_a, rho_b, beta_a, beta_b) + beta_inv_a = beta_correction(particle_system, particle) + beta_inv_b = beta_correction(neighbor_system, neighbor) + + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = (p_a - p_a_avg) / (rho_a^2 * beta_inv_a) + P_b = (p_b - p_b_avg) / (rho_b^2 * beta_inv_b) + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + v_diff = advection_velocity(v_particle_system, particle_system, particle) - + current_velocity(v_particle_system, particle_system, particle) + + return m_b * (dot(v_diff, P_a * grad_kernel_a + P_b * grad_kernel_b)) +end + +@inline function velocity_correction(particle_system, neighbor_system, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + velocity_correction(particle_system, neighbor_system, + particle_system.particle_refinement, + pos_diff, distance, v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) +end + +@inline function velocity_correction(particle_system, neighbor_system, ::Nothing, + pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, + particle, neighbor, grad_kernel) + return zero(pos_diff) +end + +@inline function velocity_correction(particle_system, neighbor_system, + particle_refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, + rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, system, particle) + advection_velocity_a = advection_velocity(v_particle_system, system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + v_diff = momentum_velocity_a - momentum_velocity_b + v_diff_tilde = advection_velocity_a - advection_velocity_b + + beta_inv_a = beta_correction(particle_system, particle) + + return -m_b * beta_inv_a * + (dot(v_diff_tilde - v_diff, grad_kernel) * momentum_velocity_a) / rho_b end # We need a separate method for EDAC since the density is stored in `v[end-1,:]`. @@ -125,3 +257,39 @@ end grad_kernel) return dv end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + return pressure_acceleration(particle_system, particle_system.particle_refinement, + neighbor_system, particle, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, + correction) +end + +# Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + ::Nothing, particle, neighbor_system, neighbor, + m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, + distance, W_a, correction) + (; pressure_acceleration_formulation) = particle_system + return pressure_acceleration_formulation(m_a, m_b, rho_a, rho_b, p_a, p_b, W_a) +end + +@inline function pressure_acceleration(particle_system::EntropicallyDampedSPHSystem, + particle_refinement, neighbor_system, particle, + neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, + pos_diff, distance, W_a, correction) + p_a_avg = average_pressure(particle_system, particle) + p_b_avg = average_pressure(neighbor_system, neighbor) + + P_a = beta_correction(particle_system, particle) * (p_a - p_a_avg) / rho_a^2 + P_b = beta_correction(neighbor_system, neighbor) * (p_b - p_b_avg) / rho_b^2 + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (P_a * grad_kernel_a + P_b * grad_kernel_b) +end diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index eb2aa1c1d..485078c9e 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -44,12 +44,11 @@ See [Entropically Damped Artificial Compressibility for SPH](@ref edac) for more gravity-like source terms. """ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, - PF, ST, B, C} <: FluidSystem{NDIMS, IC} + PF, ST, B, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: M # Vector{ELTYPE}: [particle] density_calculator :: DC smoothing_kernel :: K - smoothing_length :: ELTYPE sound_speed :: ELTYPE viscosity :: V nu_edac :: ELTYPE @@ -59,55 +58,58 @@ struct EntropicallyDampedSPHSystem{NDIMS, ELTYPE <: Real, IC, M, DC, K, V, TV, transport_velocity :: TV source_terms :: ST buffer :: B + particle_refinement :: PR cache :: C +end - function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, - smoothing_length, sound_speed; - pressure_acceleration=inter_particle_averaged_pressure, - density_calculator=SummationDensity(), - transport_velocity=nothing, - alpha=0.5, viscosity=nothing, - acceleration=ntuple(_ -> 0.0, - ndims(smoothing_kernel)), - source_terms=nothing, buffer_size=nothing) - buffer = isnothing(buffer_size) ? nothing : - SystemBuffer(nparticles(initial_condition), buffer_size) - - initial_condition = allocate_buffer(initial_condition, buffer) +function EntropicallyDampedSPHSystem(initial_condition, smoothing_kernel, + smoothing_length, sound_speed; + pressure_acceleration=inter_particle_averaged_pressure, + density_calculator=SummationDensity(), + transport_velocity=nothing, + alpha=0.5, viscosity=nothing, + acceleration=ntuple(_ -> 0.0, + ndims(smoothing_kernel)), + particle_refinement=nothing, + source_terms=nothing, buffer_size=nothing) + buffer = isnothing(buffer_size) ? nothing : + SystemBuffer(nparticles(initial_condition), buffer_size) - NDIMS = ndims(initial_condition) - ELTYPE = eltype(initial_condition) + initial_condition = allocate_buffer(initial_condition, buffer) - mass = copy(initial_condition.mass) + NDIMS = ndims(initial_condition) + ELTYPE = eltype(initial_condition) - if ndims(smoothing_kernel) != NDIMS - throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) - end - - acceleration_ = SVector(acceleration...) - if length(acceleration_) != NDIMS - throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) - end + mass = copy(initial_condition.mass) - pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, - density_calculator, - NDIMS, ELTYPE, - nothing) - - nu_edac = (alpha * smoothing_length * sound_speed) / 8 - - cache = create_cache_density(initial_condition, density_calculator) - cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + if ndims(smoothing_kernel) != NDIMS + throw(ArgumentError("smoothing kernel dimensionality must be $NDIMS for a $(NDIMS)D problem")) + end - new{NDIMS, ELTYPE, typeof(initial_condition), typeof(mass), - typeof(density_calculator), typeof(smoothing_kernel), typeof(viscosity), - typeof(transport_velocity), typeof(pressure_acceleration), typeof(source_terms), - typeof(buffer), - typeof(cache)}(initial_condition, mass, density_calculator, smoothing_kernel, - smoothing_length, sound_speed, viscosity, nu_edac, acceleration_, - nothing, pressure_acceleration, transport_velocity, source_terms, - buffer, cache) + acceleration_ = SVector(acceleration...) + if length(acceleration_) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) end + + pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, + density_calculator, + NDIMS, ELTYPE, + nothing) + + nu_edac = (alpha * smoothing_length * sound_speed) / + (2 * ndims(initial_condition) + 4) + + cache = create_cache_density(initial_condition, density_calculator) + cache = (; create_cache_edac(initial_condition, transport_velocity)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) + + return EntropicallyDampedSPHSystem(initial_condition, mass, density_calculator, + smoothing_kernel, sound_speed, viscosity, nu_edac, + acceleration_, nothing, pressure_acceleration, + transport_velocity, source_terms, buffer, + particle_refinement, cache) end function Base.show(io::IO, system::EntropicallyDampedSPHSystem) diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index e0f07796c..f287cabba 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -13,7 +13,29 @@ function create_cache_density(ic, ::ContinuityDensity) return (;) end -@inline hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] +function create_cache_refinement(initial_condition, ::Nothing, smoothing_length) + return (; smoothing_length) +end + +function create_cache_refinement(initial_condition, refinement, smoothing_length) + smoothng_length_factor = smoothing_length / initial_condition.particle_spacing + return (; smoothing_length=smoothing_length * ones(length(initial_condition.density)), + smoothing_length_factor=smoothng_length_factor) +end + +@propagate_inbounds hydrodynamic_mass(system::FluidSystem, particle) = system.mass[particle] + +function smoothing_length(system::FluidSystem, particle) + return smoothing_length(system, system.particle_refinement, particle) +end + +function smoothing_length(system::FluidSystem, ::Nothing, particle) + return system.cache.smoothing_length +end + +function smoothing_length(system::FluidSystem, refinement, particle) + return system.cache.smoothing_length[particle] +end function write_u0!(u0, system::FluidSystem) (; initial_condition) = system @@ -89,7 +111,7 @@ include("entropically_damped_sph/entropically_damped_sph.jl") add_velocity!(du, v, particle, system, system.transport_velocity) end -@inline function momentum_convection(system, neighbor_system, +@inline function momentum_convection(system, neighbor_system, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return zero(grad_kernel) @@ -98,9 +120,11 @@ end @inline function momentum_convection(system, neighbor_system::Union{EntropicallyDampedSPHSystem, WeaklyCompressibleSPHSystem}, + pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) momentum_convection(system, neighbor_system, system.transport_velocity, - v_particle_system, v_neighbor_system, rho_a, rho_b, - m_a, m_b, particle, neighbor, grad_kernel) + system.particle_refinement, pos_diff, distance, v_particle_system, + v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, + grad_kernel) end diff --git a/src/schemes/fluid/pressure_acceleration.jl b/src/schemes/fluid/pressure_acceleration.jl index 5b0b82e11..62105cf6c 100644 --- a/src/schemes/fluid/pressure_acceleration.jl +++ b/src/schemes/fluid/pressure_acceleration.jl @@ -107,7 +107,7 @@ function choose_pressure_acceleration_formulation(pressure_acceleration::Nothing end # Formulation using symmetric gradient formulation for corrections not depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction) (; pressure_acceleration_formulation) = particle_system @@ -118,7 +118,7 @@ end end # Formulation using asymmetric gradient formulation for corrections depending on local neighborhood. -@inline function pressure_acceleration(particle_system, neighbor_system, neighbor, +@inline function pressure_acceleration(particle_system, neighbor_system, particle, neighbor, m_a, m_b, p_a, p_b, rho_a, rho_b, pos_diff, distance, W_a, correction::Union{KernelCorrection, diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 66140775e..9b736795e 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -105,10 +105,8 @@ function calc_normal_akinci!(system, neighbor_system::FluidSystem, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance m_b = hydrodynamic_mass(neighbor_system, neighbor) - density_neighbor = particle_density(v_neighbor_system, - neighbor_system, neighbor) - grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, - particle) + density_neighbor = particle_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] * smoothing_length diff --git a/src/schemes/fluid/transport_velocity.jl b/src/schemes/fluid/transport_velocity.jl index bc8fdf005..5270d1210 100644 --- a/src/schemes/fluid/transport_velocity.jl +++ b/src/schemes/fluid/transport_velocity.jl @@ -76,13 +76,14 @@ end return SVector(ntuple(@inline(dim->v[ndims(system) + dim, particle]), ndims(system))) end -@inline function momentum_convection(system, neighbor_system, ::Nothing, +@inline function momentum_convection(system, neighbor_system, ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) return SVector(ntuple(_ -> 0.0, Val(ndims(system)))) end @inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + ::Nothing, pos_diff, distance, v_particle_system, v_neighbor_system, rho_a, rho_b, m_a, m_b, particle, neighbor, grad_kernel) volume_a = m_a / rho_a @@ -101,6 +102,28 @@ end return volume_term * (0.5 * (A_a + A_b)) * grad_kernel end +@inline function momentum_convection(system, neighbor_system, ::TransportVelocityAdami, + refinement, pos_diff, distance, + v_particle_system, v_neighbor_system, rho_a, rho_b, + m_a, m_b, particle, neighbor, grad_kernel) + momentum_velocity_a = current_velocity(v_particle_system, system, particle) + advection_velocity_a = advection_velocity(v_particle_system, system, particle) + + momentum_velocity_b = current_velocity(v_neighbor_system, neighbor_system, neighbor) + advection_velocity_b = advection_velocity(v_neighbor_system, neighbor_system, neighbor) + + factor_a = beta_correction(particle_system, particle) / rho_a + factor_b = beta_correction(neighbor_system, neighbor) / rho_b + + A_a = factor_a * momentum_velocity_a * (advection_velocity_a - momentum_velocity_a)' + A_b = factor_b * momentum_velocity_b * (advection_velocity_b - momentum_velocity_b)' + + grad_kernel_a = smoothing_kernel_grad(system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + return -m_b * (A_a * grad_kernel_a + A_b * grad_kernel_b) +end + @inline transport_velocity!(dv, system, rho_a, rho_b, m_a, m_b, grad_kernel, particle) = dv @inline function transport_velocity!(dv, system::FluidSystem, diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index 09bde77f9..222fa6fd8 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -1,9 +1,9 @@ # Unpack the neighboring systems viscosity to dispatch on the viscosity type -@inline function dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) viscosity = viscosity_model(particle_system, neighbor_system) return dv_viscosity(viscosity, particle_system, neighbor_system, @@ -12,10 +12,10 @@ sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) end -@inline function dv_viscosity(viscosity, particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +@propagate_inbounds function dv_viscosity(viscosity, particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) return viscosity(particle_system, neighbor_system, v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, @@ -105,27 +105,35 @@ function kinematic_viscosity(system, viscosity::ViscosityMorris) return viscosity.nu end -@inline function (viscosity::Union{ArtificialViscosityMonaghan, - ViscosityMorris})(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, - distance, sound_speed, m_a, m_b, - rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system - +@propagate_inbounds function (viscosity::Union{ArtificialViscosityMonaghan, + ViscosityMorris})(particle_system, + neighbor_system, + v_particle_system, + v_neighbor_system, + particle, neighbor, + pos_diff, distance, + sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) rho_mean = 0.5 * (rho_a + rho_b) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) v_diff = v_a - v_b + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? nu_a = kinematic_viscosity(particle_system, - viscosity_model(neighbor_system, particle_system)) + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) + # TODO do we need the average smoothing length here? pi_ab = viscosity(sound_speed, v_diff, pos_diff, distance, rho_mean, rho_a, rho_b, - smoothing_length, grad_kernel, nu_a, nu_b) + smoothing_length_particle, grad_kernel, nu_a, nu_b) return m_b * pi_ab end @@ -166,12 +174,12 @@ end # Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". # In: Reports on Progress in Physics (2005), pages 1703-1759. # [doi: 10.1088/0034-4885/68/8/r01](http://dx.doi.org/10.1088/0034-4885/68/8/R01) -function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan) - (; smoothing_length) = system +function kinematic_viscosity(system, viscosity::ArtificialViscosityMonaghan, + smoothing_length_particle) (; alpha) = viscosity sound_speed = system_sound_speed(system) - return alpha * smoothing_length * sound_speed / (2 * ndims(system) + 4) + return alpha * smoothing_length_particle * sound_speed / (2 * ndims(system) + 4) end @doc raw""" @@ -209,13 +217,50 @@ end particle, neighbor, pos_diff, distance, sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) - (; smoothing_length) = particle_system + viscosity(particle_system, neighbor_system, particle_system.particle_refinement, + v_particle_system, v_neighbor_system, particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) +end + +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, + particle_refinement, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + beta_inv_a = beta_correction(particle_system, particle_refinement, particle) - epsilon = viscosity.epsilon nu_a = kinematic_viscosity(particle_system, viscosity_model(neighbor_system, particle_system)) + + grad_kernel_a = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) + grad_kernel_b = smoothing_kernel_grad(neighbor_system, pos_diff, distance, neighbor) + + grad_W_avg = 0.5 * (grad_kernel_a + grad_kernel_b) + + tmp = (distance^2 + 0.001 * smoothing_length(particle_system, particle)^2) + + return m_b * beta_inv_a * 4 * nu_a * dot(pos_diff, grad_W_avg) * v_diff / + (tmp * (rho_a + rho_b)) +end + +@inline function (viscosity::ViscosityAdami)(particle_system, neighbor_system, ::Nothing, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, + distance, sound_speed, m_a, m_b, + rho_a, rho_b, grad_kernel) + epsilon = viscosity.epsilon + + smoothing_length_particle = smoothing_length(particle_system, particle) + smoothing_length_neighbor = smoothing_length(particle_system, neighbor) + + # TODO is that correct with the smoothing lengths? + nu_a = kinematic_viscosity(particle_system, + viscosity_model(neighbor_system, particle_system), + smoothing_length_particle) nu_b = kinematic_viscosity(neighbor_system, - viscosity_model(particle_system, neighbor_system)) + viscosity_model(particle_system, neighbor_system), + smoothing_length_neighbor) v_a = viscous_velocity(v_particle_system, particle_system, particle) v_b = viscous_velocity(v_neighbor_system, neighbor_system, neighbor) @@ -226,8 +271,8 @@ end eta_tilde = 2 * (eta_a * eta_b) / (eta_a + eta_b) - # TODO For variable smoothing_length use average smoothing length - tmp = eta_tilde / (distance^2 + epsilon * smoothing_length^2) + smoothing_length_average = 0.5 * (smoothing_length_particle + smoothing_length_neighbor) + tmp = eta_tilde / (distance^2 + epsilon * smoothing_length_average^2) volume_a = m_a / rho_a volume_b = m_b / rho_b @@ -246,8 +291,10 @@ end return visc .* v_diff end -function kinematic_viscosity(system, viscosity::ViscosityAdami) +function kinematic_viscosity(system, viscosity::ViscosityAdami, particle) return viscosity.nu end -@inline viscous_velocity(v, system, particle) = current_velocity(v, system, particle) +@propagate_inbounds function viscous_velocity(v, system, particle) + return current_velocity(v, system, particle) +end diff --git a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl index c4ddc8b81..33a1b7787 100644 --- a/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl +++ b/src/schemes/fluid/weakly_compressible_sph/density_diffusion.jl @@ -75,9 +75,9 @@ end @inline function density_diffusion_psi(::DensityDiffusionFerrari, rho_a, rho_b, pos_diff, distance, system, particle, neighbor) - (; smoothing_length) = system - - return ((rho_a - rho_b) / 2smoothing_length) * pos_diff / distance + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(system, particle) + return ((rho_a - rho_b) / (2 * smoothing_length_particle)) * pos_diff / distance end @doc raw""" @@ -178,7 +178,7 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search foreach_point_neighbor(system, system, system_coords, system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance # Only consider particles with a distance > 0 - distance < sqrt(eps()) && return + distance < sqrt(eps(typeof(distance))) && return rho_a = particle_density(v, system, particle) rho_b = particle_density(v, system, neighbor) @@ -199,15 +199,15 @@ function update!(density_diffusion::DensityDiffusionAntuono, neighborhood_search return density_diffusion end -@inline function density_diffusion!(dv, density_diffusion::DensityDiffusion, - v_particle_system, particle, neighbor, pos_diff, - distance, m_b, rho_a, rho_b, - particle_system::FluidSystem, grad_kernel) +@propagate_inbounds function density_diffusion!(dv, density_diffusion::DensityDiffusion, + v_particle_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system::FluidSystem, grad_kernel) # Density diffusion terms are all zero for distance zero - distance < sqrt(eps()) && return + distance < sqrt(eps(typeof(distance))) && return (; delta) = density_diffusion - (; smoothing_length, state_equation) = particle_system + (; state_equation) = particle_system (; sound_speed) = state_equation volume_b = m_b / rho_b @@ -216,7 +216,10 @@ end particle_system, particle, neighbor) density_diffusion_term = dot(psi, grad_kernel) * volume_b - dv[end, particle] += delta * smoothing_length * sound_speed * density_diffusion_term + # TODO should we use the average smoothing length here? + smoothing_length_particle = smoothing_length(particle_system, particle) + dv[end, particle] += delta * smoothing_length_particle * sound_speed * + density_diffusion_term end # Density diffusion `nothing` or interaction other than fluid-fluid diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 9e40eacfa..b874737bf 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -23,19 +23,23 @@ function interact!(dv, v_particle_system, u_particle_system, foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_system_coords, neighborhood_search) do particle, neighbor, pos_diff, distance - rho_a = particle_density(v_particle_system, particle_system, particle) - rho_b = particle_density(v_neighbor_system, neighbor_system, neighbor) - rho_mean = 0.5 * (rho_a + rho_b) - - # Determine correction values - viscosity_correction, pressure_correction, surface_tension_correction = free_surface_correction(correction, - particle_system, - rho_mean) + # `foreach_point_neighbor` makes sure that `particle` and `neighbor` are + # in bounds of the respective system. For performance reasons, we use `@inbounds` + # in this hot loop to avoid bounds checking when extracting particle quantities. + rho_a = @inbounds particle_density(v_particle_system, particle_system, particle) + rho_b = @inbounds particle_density(v_neighbor_system, neighbor_system, neighbor) + rho_mean = (rho_a + rho_b) / 2 + + # Determine correction factors. + # This can be ignored, as these are all 1 when no correction is used. + (viscosity_correction, pressure_correction, + surface_tension_correction) = free_surface_correction(correction, particle_system, + rho_mean) grad_kernel = smoothing_kernel_grad(particle_system, pos_diff, distance, particle) - m_a = hydrodynamic_mass(particle_system, particle) - m_b = hydrodynamic_mass(neighbor_system, neighbor) + m_a = @inbounds hydrodynamic_mass(particle_system, particle) + m_b = @inbounds hydrodynamic_mass(neighbor_system, neighbor) # The following call is equivalent to # `p_a = particle_pressure(v_particle_system, particle_system, particle)` @@ -43,20 +47,24 @@ function interact!(dv, v_particle_system, u_particle_system, # Only when the neighbor system is a `BoundarySPHSystem` or a `TotalLagrangianSPHSystem` # with the boundary model `PressureMirroring`, this will return `p_b = p_a`, which is # the pressure of the fluid particle. - p_a, p_b = particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) + p_a, p_b = @inbounds particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) dv_pressure = pressure_correction * - pressure_acceleration(particle_system, neighbor_system, neighbor, + 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) + # Propagate `@inbounds` to the viscosity function, which accesses particle data dv_viscosity_ = viscosity_correction * - dv_viscosity(particle_system, neighbor_system, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - sound_speed, m_a, m_b, rho_a, rho_b, grad_kernel) + @inbounds dv_viscosity(particle_system, neighbor_system, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + sound_speed, m_a, m_b, rho_a, rho_b, + grad_kernel) dv_surface_tension = surface_tension_correction * surface_tension_force(surface_tension_a, surface_tension_b, @@ -66,17 +74,19 @@ function interact!(dv, v_particle_system, u_particle_system, dv_adhesion = adhesion_force(surface_tension, particle_system, neighbor_system, particle, neighbor, pos_diff, distance) - @inbounds for i in 1:ndims(particle_system) - dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + dv_surface_tension[i] + - dv_adhesion[i] + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += dv_pressure[i] + dv_viscosity_[i] + + dv_surface_tension[i] + dv_adhesion[i] # Debug example # debug_array[i, particle] += dv_pressure[i] end # TODO If variable smoothing_length is used, this should use the neighbor smoothing length - continuity_equation!(dv, density_calculator, v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, m_b, rho_a, rho_b, - particle_system, neighbor_system, grad_kernel) + # Propagate `@inbounds` to the continuity equation, which accesses particle data + @inbounds continuity_equation!(dv, density_calculator, v_particle_system, + v_neighbor_system, particle, neighbor, + pos_diff, distance, m_b, rho_a, rho_b, + particle_system, neighbor_system, grad_kernel) end # Debug example # periodic_box = neighborhood_search.periodic_box @@ -98,12 +108,12 @@ end end # This formulation was chosen to be consistent with the used pressure_acceleration formulations. -@inline function continuity_equation!(dv, density_calculator::ContinuityDensity, - v_particle_system, v_neighbor_system, - particle, neighbor, pos_diff, distance, - m_b, rho_a, rho_b, - particle_system::WeaklyCompressibleSPHSystem, - neighbor_system, grad_kernel) +@propagate_inbounds function continuity_equation!(dv, density_calculator::ContinuityDensity, + v_particle_system, v_neighbor_system, + particle, neighbor, pos_diff, distance, + m_b, rho_a, rho_b, + particle_system::WeaklyCompressibleSPHSystem, + neighbor_system, grad_kernel) (; density_diffusion) = particle_system vdiff = current_velocity(v_particle_system, particle_system, particle) - @@ -121,9 +131,10 @@ end end end -@inline function particle_neighbor_pressure(v_particle_system, v_neighbor_system, - particle_system, neighbor_system, - particle, neighbor) +@propagate_inbounds function particle_neighbor_pressure(v_particle_system, + v_neighbor_system, + particle_system, neighbor_system, + particle, neighbor) p_a = particle_pressure(v_particle_system, particle_system, particle) p_b = particle_pressure(v_neighbor_system, neighbor_system, neighbor) diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eb6c06c84..3ea8352f7 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -44,15 +44,14 @@ See [Weakly Compressible SPH](@ref wcsph) for more details on the method. """ -struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, - V, DD, COR, PF, ST, B, SRFT, C} <: FluidSystem{NDIMS, IC} +struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, V, DD, COR, + PF, ST, B, SRFT, PR, C} <: FluidSystem{NDIMS, IC} initial_condition :: IC mass :: MA # Array{ELTYPE, 1} pressure :: P # Array{ELTYPE, 1} density_calculator :: DC state_equation :: SE smoothing_kernel :: K - smoothing_length :: ELTYPE acceleration :: SVector{NDIMS, ELTYPE} viscosity :: V density_diffusion :: DD @@ -62,6 +61,7 @@ struct WeaklyCompressibleSPHSystem{NDIMS, ELTYPE <: Real, IC, MA, P, DC, SE, K, source_terms :: ST surface_tension :: SRFT buffer :: B + particle_refinement :: PR # TODO cache :: C end @@ -76,6 +76,7 @@ function WeaklyCompressibleSPHSystem(initial_condition, acceleration=ntuple(_ -> 0.0, ndims(smoothing_kernel)), correction=nothing, source_terms=nothing, + particle_refinement=nothing, surface_tension=nothing) buffer = isnothing(buffer_size) ? nothing : SystemBuffer(nparticles(initial_condition), buffer_size) @@ -116,14 +117,18 @@ function WeaklyCompressibleSPHSystem(initial_condition, cache = (; create_cache_wcsph(surface_tension, ELTYPE, NDIMS, n_particles)..., cache...) + cache = (; + create_cache_refinement(initial_condition, particle_refinement, + smoothing_length)..., cache...) return WeaklyCompressibleSPHSystem(initial_condition, mass, pressure, density_calculator, state_equation, - smoothing_kernel, smoothing_length, + smoothing_kernel, acceleration_, viscosity, density_diffusion, correction, pressure_acceleration, nothing, - source_terms, surface_tension, buffer, cache) + source_terms, surface_tension, buffer, + particle_refinement, cache) end create_cache_wcsph(correction, density, NDIMS, nparticles) = (;) @@ -211,7 +216,8 @@ end return ndims(system) + 1 end -@inline function particle_pressure(v, system::WeaklyCompressibleSPHSystem, particle) +@propagate_inbounds function particle_pressure(v, system::WeaklyCompressibleSPHSystem, + particle) return system.pressure[particle] end diff --git a/src/schemes/solid/total_lagrangian_sph/rhs.jl b/src/schemes/solid/total_lagrangian_sph/rhs.jl index 4da0f7783..114f72aed 100644 --- a/src/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/src/schemes/solid/total_lagrangian_sph/rhs.jl @@ -32,7 +32,7 @@ end rho_b = neighbor_system.material_density[neighbor] grad_kernel = smoothing_kernel_grad(particle_system, initial_pos_diff, - initial_distance) + initial_distance, particle) m_a = particle_system.mass[particle] m_b = neighbor_system.mass[neighbor] @@ -88,7 +88,7 @@ function interact!(dv, v_particle_system, u_particle_system, # Use kernel from the fluid system in order to get the same force here in # solid-fluid interaction as for fluid-solid interaction. - grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance) + grad_kernel = smoothing_kernel_grad(neighbor_system, pos_diff, distance, particle) # In fluid-solid interaction, use the "hydrodynamic pressure" of the solid particles # corresponding to the chosen boundary model. @@ -99,7 +99,8 @@ function interact!(dv, v_particle_system, u_particle_system, # are switched in the following two calls. # This way, we obtain the exact same force as for the fluid-solid interaction, # but with a flipped sign (because `pos_diff` is flipped compared to fluid-solid). - dv_boundary = pressure_acceleration(neighbor_system, particle_system, particle, + dv_boundary = pressure_acceleration(neighbor_system, particle_system, + neighbor, particle, m_b, m_a, p_b, p_a, rho_b, rho_a, pos_diff, distance, grad_kernel, neighbor_system.correction) diff --git a/src/schemes/solid/total_lagrangian_sph/system.jl b/src/schemes/solid/total_lagrangian_sph/system.jl index dfd71a2d6..8275da22b 100644 --- a/src/schemes/solid/total_lagrangian_sph/system.jl +++ b/src/schemes/solid/total_lagrangian_sph/system.jl @@ -286,7 +286,7 @@ end pos_diff = current_coords(system, particle) - current_coords(system, neighbor) grad_kernel = smoothing_kernel_grad(system, initial_pos_diff, - initial_distance) + initial_distance, particle) result = volume * pos_diff * grad_kernel' diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index c666ec3d7..9b31c463c 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -177,8 +177,8 @@ function sample_plane(plane_points::NTuple{3}, particle_spacing; tlsph=nothing) edge2 = point3_ - point1_ # Check if the points are collinear - if norm(cross(edge1, edge2)) == 0 - throw(ArgumentError("the points must not be collinear")) + if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + throw(ArgumentError("the vectors `AB` and `AC` must not be collinear")) end # Determine the number of points along each edge diff --git a/src/visualization/write2vtk.jl b/src/visualization/write2vtk.jl index ef5f68cc6..90c5fa435 100644 --- a/src/visualization/write2vtk.jl +++ b/src/visualization/write2vtk.jl @@ -252,7 +252,7 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true) vtk["viscosity"] = type2string(system.viscosity) write2vtk!(vtk, system.viscosity) vtk["smoothing_kernel"] = type2string(system.smoothing_kernel) - vtk["smoothing_length"] = system.smoothing_length + # vtk["smoothing_length"] = system.smoothing_length TODO vtk["density_calculator"] = type2string(system.density_calculator) if system isa WeaklyCompressibleSPHSystem diff --git a/test/examples/examples.jl b/test/examples/examples.jl index 5f42aa230..3673e5bd5 100644 --- a/test/examples/examples.jl +++ b/test/examples/examples.jl @@ -9,8 +9,20 @@ "hydrostatic_water_column_2d.jl"), fluid_system=nothing, sol=nothing, semi=nothing, ode=nothing) + # Neighborhood search for `FullGridCellList` test below + search_radius = TrixiParticles.compact_support(smoothing_kernel, + smoothing_length) + min_corner = minimum(tank.boundary.coordinates, dims=2) .- search_radius + max_corner = maximum(tank.boundary.coordinates, dims=2) .+ search_radius + cell_list = TrixiParticles.PointNeighbors.FullGridCellList(; min_corner, + max_corner) + semi_fullgrid = Semidiscretization(fluid_system, boundary_system, + neighborhood_search=GridNeighborhoodSearch{2}(; + cell_list)) + hydrostatic_water_column_tests = Dict( "WCSPH default" => (), + "WCSPH with FullGridCellList" => (semi=semi_fullgrid,), "WCSPH with source term damping" => (source_terms=SourceTermDamping(damping_coefficient=1e-4),), "WCSPH with SummationDensity" => (fluid_density_calculator=SummationDensity(), clip_negative_pressure=true), @@ -226,6 +238,14 @@ @test count_rhs_allocations(sol, semi) == 0 end + @trixi_testset "fluid/pipe_flow_3d.jl" begin + @test_nowarn_mod trixi_include(@__MODULE__, tspan=(0.0, 0.5), + joinpath(examples_dir(), "fluid", + "pipe_flow_3d.jl")) + @test sol.retcode == ReturnCode.Success + @test count_rhs_allocations(sol, semi) == 0 + end + @trixi_testset "fluid/lid_driven_cavity_2d.jl" begin @test_nowarn_mod trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", diff --git a/test/general/smoothing_kernels.jl b/test/general/smoothing_kernels.jl index 5d6f15c26..cb9c98f4b 100644 --- a/test/general/smoothing_kernels.jl +++ b/test/general/smoothing_kernels.jl @@ -90,4 +90,37 @@ end end end -end + + @testset verbose=false "Return Type" begin + # Test that the return type of the kernel and kernel derivative preserve + # the input type. We don't want to return `Float64` when working with `Float32`. + kernels = [ + GaussianKernel, + SchoenbergCubicSplineKernel, + SchoenbergQuarticSplineKernel, + SchoenbergQuinticSplineKernel, + WendlandC2Kernel, + WendlandC4Kernel, + WendlandC6Kernel, + SpikyKernel, + Poly6Kernel + ] + + # Test different smoothing length types + smoothing_lengths = (0.5, 0.5f0) + + @testset "$kernel_type" for kernel_type in kernels + for ndims in 2:3 + kernel_ = kernel_type{ndims}() + + for h in smoothing_lengths + result = TrixiParticles.kernel(kernel_, h / 2, h) + @test typeof(result) == typeof(h) + + result = TrixiParticles.kernel_deriv(kernel_, h / 2, h) + @test typeof(result) == typeof(h) + end + end + end + end +end; diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 27caa0293..a72ac409c 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -54,12 +54,14 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, - v_fluid .* - ones(size(fluid.coordinates)), - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, + boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid .* + ones(size(fluid.coordinates)), + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps() @@ -92,10 +94,12 @@ TrixiParticles.reset_cache!(boundary_system.boundary_model.cache, boundary_system.boundary_model.viscosity) - TrixiParticles.adami_pressure_extrapolation!(boundary_model, boundary_system, - fluid_system, boundary.coordinates, - fluid.coordinates, v_fluid, - neighborhood_search) + TrixiParticles.boundary_pressure_extrapolation!(boundary_model, boundary_system, + fluid_system, + boundary.coordinates, + fluid.coordinates, v_fluid, + v_fluid, + neighborhood_search) for particle in TrixiParticles.eachparticle(boundary_system) if volume[particle] > eps() diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index ad1b69c4b..2ed6e2810 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -187,10 +187,10 @@ end @testset verbose=true "Illegal Inputs" begin - no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.2, 2.0, -0.5]] + no_rectangular_plane = [[0.2, 0.3, -0.5], [-1.0, 1.5, 0.2], [-0.4, 0.9, -0.15]] flow_direction = [0.0, 0.0, 1.0] - error_str = "the vectors `AB` and `AC` for the provided points `A`, `B`, `C` must be orthogonal" + error_str = "the vectors `AB` and `AC` must not be collinear" @test_throws ArgumentError(error_str) InFlow(; plane=no_rectangular_plane, particle_spacing=0.1, diff --git a/test/schemes/solid/total_lagrangian_sph/rhs.jl b/test/schemes/solid/total_lagrangian_sph/rhs.jl index ea2474d89..b69f835ca 100644 --- a/test/schemes/solid/total_lagrangian_sph/rhs.jl +++ b/test/schemes/solid/total_lagrangian_sph/rhs.jl @@ -114,6 +114,7 @@ nothing end TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) = kernel_deriv + Base.eps(::Type{Val{:mock_smoothing_length}}) = eps() #### Verification systems = [system, system_gpu] diff --git a/test/systems/solid_system.jl b/test/systems/solid_system.jl index e5f0c0e05..18b887b49 100644 --- a/test/systems/solid_system.jl +++ b/test/systems/solid_system.jl @@ -131,7 +131,7 @@ Base.ntuple(f, ::Symbol) = ntuple(f, 2) # Make `extract_svector` work function TrixiParticles.current_coords(system::Val{:mock_system_tensor}, particle) - return TrixiParticles.extract_svector(current_coordinates[i], system, + return TrixiParticles.extract_svector(current_coordinates[i], Val(2), particle) end @@ -171,6 +171,7 @@ function TrixiParticles.kernel_deriv(::Val{:mock_smoothing_kernel}, _, _) return kernel_derivative end + Base.eps(::Type{Val{:mock_smoothing_length}}) = eps() # Compute deformation gradient deformation_grad = ones(2, 2, 2) @@ -371,4 +372,4 @@ @test isapprox(von_mises_stress[1], 1.4257267477533202, atol=1e-14) @test isapprox(reference_stress_tensor, cauchy_stress, atol=1e-6) end -end +end;