diff --git a/NEWS.md b/NEWS.md index 8abe050384..4a2cb0d941 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,11 +4,25 @@ TrixiParticles.jl follows the interpretation of [semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1) used in the Julia ecosystem. Notable changes will be documented in this file for human readability. -## Version 0.5.3 - -### Features - -- Added the computation of boundary normals for `RectangularTank`s and `SphereShape`s. +## Version 0.5.3 + +### API Changes + +- Corrected `SurfaceTensionMorris` to apply its local CSF acceleration once per particle and + retain the required one-phase surface delta. Previous coefficients compensated implicitly + for a dimensionally incomplete force repeated once per fluid neighbor and must be recalibrated. +- For `SurfaceTensionMorris` with `ColorfieldSurfaceNormal`, `ideal_density_threshold` now denotes + a fraction of the continuous complete-support kernel moment instead of an integer neighbor-count + fraction. The default zero still disables interior filtering. + +### Features + +- Added C1 interface activation for Morris CSF with `ColorfieldSurfaceNormal`. Color-gradient and + continuous support-moment indicators taper the physical surface delta without another neighbor + pass. +- Added `CorrectedCSFSurfaceNormal`, an explicit free-surface implementation of the C-CSF + interface geometry from Vergnaud et al. (2022) for `SurfaceTensionMorris`. +- Added the computation of boundary normals for `RectangularTank`s and `SphereShape`s. - Added `flush` keyword argument to `InfoCallback` to flush `stdout` after each output, useful for monitoring progress in real-time on clusters or batch systems (#1246). - Added the number of split integration time steps to the `InfoCallback` output diff --git a/README.md b/README.md index bf3c1c9cbc..415b72c23f 100644 --- a/README.md +++ b/README.md @@ -45,24 +45,24 @@ It offers intuitive configuration, robust pre- and post-processing, and vendor-a We provide several example simulation setups in the `examples` folder (which can be accessed from Julia via `examples_dir()`). - - - -
-
2D Dam Break
+
+
2D Dam Break
 -
Moving Wall
+
+
Moving Wall
-
Oscillating Beam
+
+
Oscillating Beam
 -
Dam Break with Elastic Plate
+
+
Dam Break with Elastic Plate
diff --git a/docs/literate/src/tut_2d_geometry.jl b/docs/literate/src/tut_2d_geometry.jl new file mode 100644 index 0000000000..3028e53527 --- /dev/null +++ b/docs/literate/src/tut_2d_geometry.jl @@ -0,0 +1,155 @@ +# # [Setting up a 2D simulation from geometry files](@id tut_2d_geometry) + +# In this tutorial, we build two genuine 2D setups from geometry files: +# 1. a curved pipe, where one geometry file defines the outer wall envelope and a second +# one defines the empty channel cut out of it, +# 2. a dam-break basin with a coastline profile, where one geometry file defines the +# filled coastline wall together with the seawall on the right. +# +# For a real 2D setup, we use 2D geometry formats such as `.asc` or `.dxf`. +# STL files are surface meshes and therefore naturally lead to thin 3D setups instead. + +# First, we import TrixiParticles.jl together with +# `OrdinaryDiffEqLowStorageRK` of +# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) +# and [Plots.jl](https://docs.juliaplots.org/stable/). +using TrixiParticles +using OrdinaryDiffEqLowStorageRK +using Plots + +# ## Resolution + +# We use the same particle spacing for the fluid and for the wall geometries. +particle_spacing = 0.03 +fluid_density = 1000.0 +gravity = 9.81 +sound_speed = 10.0 +state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density, + exponent=7) +nothing # hide + +# ## Loading 2D geometry files + +# The following helper loads a closed 2D geometry file and samples particles in its interior: +# 1. load the polygon with [`load_geometry`](@ref), +# 2. fill the polygon with [`ComplexShape`](@ref). +# +# This creates a true 2D solid region instead of a hollow shell around the polygon edges. +function solid_from_geometry_file(file; particle_spacing, density) + geometry = load_geometry(file) + solid = ComplexShape(geometry; particle_spacing, density, + grid_offset=0.5particle_spacing) + + return (; geometry, solid) +end + +# ## A curved pipe from two filled geometries + +# The pipe wall is a solid L-shaped region with a channel cut out of it: +# 1. one geometry file describes the outer pipe envelope, +# 2. one geometry file describes the empty channel, +# 3. the `setdiff` operation subtracts the channel from the solid envelope. +pipe_outer_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", + "curved_pipe_outer_2d.asc") +pipe_channel_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", + "curved_pipe_channel_2d.asc") + +pipe_outer = solid_from_geometry_file(pipe_outer_file; particle_spacing, + density=fluid_density) +pipe_channel = load_geometry(pipe_channel_file) + +pipe_setup = (; wall=setdiff(pipe_outer.solid, pipe_channel), + outer_geometry=pipe_outer.geometry, + channel_geometry=pipe_channel) + +# ## A dam-break basin with a coastline profile + +# In the second setup, a single 2D geometry file defines a filled coastline wall: +# the beach profile on top, a finite wall thickness below it, and the seawall on the right. +coast_file = pkgdir(TrixiParticles, "examples", "preprocessing", "data", + "coastline_profile_2d.asc") +coast = solid_from_geometry_file(coast_file; particle_spacing, density=fluid_density) + +# The geometry file gives the coastline bed and the right wall as a solid region. +# We add the left wall explicitly as a rectangular particle block and place a +# 1.5x taller rectangular dam-break water column next to it. +left_wall = RectangularShape(particle_spacing, (5, 50), (0.0, -0.12), + density=fluid_density) +reservoir = RectangularShape(particle_spacing, (28, 42), (0.15, 0.03), + acceleration=(0.0, -gravity), + state_equation=state_equation) +coast_setup = (; geometry=coast.geometry, + wall=union(coast.solid, left_wall), + fluid=setdiff(reservoir, coast.geometry)) + +p_pipe = plot(pipe_setup.wall, label="wall", title="Curved pipe", + markerstrokewidth=0, markersize=4) +plot!(p_pipe, showaxis=false, aspect_ratio=:equal, + xlims=(-0.03, 1.23), ylims=(-0.03, 1.23)) + +p_coast = plot(coast_setup.fluid, coast_setup.wall, + labels=["fluid" "wall"], title="Coastline dam break", + markerstrokewidth=0, markersize=3) +plot!(p_coast, showaxis=false, aspect_ratio=:equal, + xlims=(0.0, 2.75), ylims=(-0.15, 1.35)) + +plot(p_pipe, p_coast, layout=(1, 2), size=(900, 360)) +savefig("tut_2d_geometry_plot.png"); # hide +# ![2D geometry based initial conditions](tut_2d_geometry_plot.png) + +# ## Building the simulation systems + +# To keep the example focused, we continue with the coastline setup. +# From this point on, the simulation setup is the same as in other 2D simulation files. +setup = coast_setup +tspan = (0.0, 0.03) +nothing # hide + +# We define the state equation, smoothing kernel, and viscosity for a +# weakly compressible SPH simulation. +smoothing_length = 1.2 * particle_spacing +smoothing_kernel = SchoenbergCubicSplineKernel{2}() +viscosity = ArtificialViscosityMonaghan(alpha=0.02, beta=0.0) + +fluid_density_calculator = ContinuityDensity() +density_diffusion = DensityDiffusionMolteniColagrossi(delta=0.1) + +fluid_system = WeaklyCompressibleSPHSystem(setup.fluid; + density_calculator=fluid_density_calculator, + state_equation, smoothing_kernel, + smoothing_length, viscosity=viscosity, + density_diffusion=density_diffusion, + acceleration=(0.0, -gravity)) +nothing # hide + +# For the wall, we reuse the combined solid wall particles created above. The high-level +# constructor obtains the smoothing kernel, smoothing length, and state equation from the fluid. +boundary_model = BoundaryModelDummyParticles(setup.wall; fluid_system) +boundary_system = WallBoundarySystem(setup.wall, boundary_model) +nothing # hide + +# ## Semidiscretization + +# With fluid and wall particles defined, we can build the +# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization) exactly as in other tutorials. +semi = Semidiscretization(fluid_system, boundary_system) +ode = semidiscretize(semi, tspan) +nothing # hide + +# ## Time integration + +# The setup is now complete. +# To start the simulation, run for example +# ```julia +# callbacks = CallbackSet(InfoCallback(interval=10)) +# sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) +# ``` +# This is the same final step as in [the basic setup tutorial](@ref tut_setup). +callbacks = CallbackSet(InfoCallback(interval=10)) +nothing # hide + +sol = solve(ode, RDPK3SpFSAL35(), save_everystep=false, callback=callbacks) #!md + +# For more accurate body-fitted particles around sharper features, you can also +# apply the [particle packing workflow](@ref tut_packing) to the 2D geometry files +# before starting the simulation. diff --git a/docs/literate/src/tut_packing.jl b/docs/literate/src/tut_packing.jl index d5b9610283..9ae98729d2 100644 --- a/docs/literate/src/tut_packing.jl +++ b/docs/literate/src/tut_packing.jl @@ -75,8 +75,8 @@ plot!(right_margin=5Plots.mm) #hide # ## Creating an initial configuration of boundary particles # To create the initial configuration of the boundary particles, -# we use the sampled points of the SDF whose signed distance lies between 0 -# and `boundary_thickness`. +# we use the sampled points of the SDF whose signed distance lies between the +# geometry offset implied by `place_on_shell` and `boundary_thickness`. # Here, we need to specify the `density` of the boundary particles. # As an example, we choose `1.0` for all particles. # This gives us an [`InitialCondition`](@ref InitialCondition) for the boundary particles. @@ -125,7 +125,8 @@ plot!(geometry, linestyle=:dash, label=nothing, showaxis=false, color=:black, # ## Particle packing # In the following, we will essentially follow the same steps described in the fluid tutorials. -# That means we will generate systems that are then passed to the [`Semidiscretization`](@ref). +# That means we will generate systems that are then passed to the +# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization). # The difference from a typical physical simulation is that we use [`ParticlePackingSystem`](@ref), # which does not represent any physical law. Instead, we only use the simulation framework to time-integrate # the packing process. @@ -211,7 +212,7 @@ plot!(geometry, seriestype=:path, color=:black, label=nothing, linewidth=2) boundary_system = ParticlePackingSystem(boundary_sampled; is_boundary=true, smoothing_kernel, smoothing_length, boundary_compress_factor=0.7, signed_distance_field, - background_pressure) + boundary_thickness, background_pressure) # We can now couple the boundary system with the interior system: semi = Semidiscretization(packing_system, boundary_system) @@ -251,10 +252,12 @@ fixed_system = ParticlePackingSystem(packed_ic; smoothing_kernel, smoothing_leng # Now we define a rectangular domain that we want to pack. # In practice, you could create any `InitialCondition` that encloses your complex geometry. -tank_domain = RectangularTank(particle_spacing, (4, 4), (0, 0), min_coordinates=(-1, -2), - density) +domain_size = (4, 4) +n_particles_per_dimension = round.(Int, domain_size ./ particle_spacing) +tank_domain = RectangularShape(particle_spacing, n_particles_per_dimension, (-1, -2); + density) -sampled_outer_domain = setdiff(tank_domain.fluid, packed_ic) +sampled_outer_domain = setdiff(tank_domain, packed_ic) # If we plot these two `InitialCondition`s, we can see # that the geometry interface is not properly represented yet. diff --git a/docs/literate/src/tut_rigid_body_fsi.jl b/docs/literate/src/tut_rigid_body_fsi.jl index e87a17e88e..ea09268ed5 100644 --- a/docs/literate/src/tut_rigid_body_fsi.jl +++ b/docs/literate/src/tut_rigid_body_fsi.jl @@ -165,12 +165,8 @@ nothing # hide # See [the docs on dummy particles](@ref boundary_models) for a definition for these terms. boundary_density_calculator = AdamiPressureExtrapolation() -tank_boundary_model = BoundaryModelDummyParticles(tank.boundary.density, - tank.boundary.mass, - boundary_density_calculator, - fluid_smoothing_kernel, - fluid_smoothing_length; - state_equation) +tank_boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, + boundary_density_calculator) boundary_system = WallBoundarySystem(tank.boundary, tank_boundary_model) nothing # hide diff --git a/docs/literate/src/tut_setup.jl b/docs/literate/src/tut_setup.jl index b546c5eb89..2372e45b04 100644 --- a/docs/literate/src/tut_setup.jl +++ b/docs/literate/src/tut_setup.jl @@ -131,20 +131,18 @@ nothing # hide # To model the boundary, we use particle-based boundary conditions, in which particles # are sampled in the boundary that interact with the fluid particles to avoid penetration. -# In order to define a boundary system, we first have to choose a boundary model, -# which defines how the fluid interacts with boundary particles. -# We will use the [`BoundaryModelDummyParticles`](@ref) with [`AdamiPressureExtrapolation`](@ref). -# See [here](@ref boundary_models) for a comprehensive overview over boundary models. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length; - state_equation) +# Here, we explicitly choose the dummy-particle boundary model and use its high-level +# builder to infer kernel and equation-of-state-related settings from the adjacent +# fluid system. See [here](@ref boundary_models) for a comprehensive overview over +# boundary models. +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) nothing # hide # ## [Semidiscretization](@id tut_setup_semi) -# The key component of every simulation is the [`Semidiscretization`](@ref), +# The key component of every simulation is the +# [`Semidiscretization`](@ref TrixiParticles.Semidiscretization), # which couples all systems of the simulation. # All simulation methods in TrixiParticles.jl are semidiscretizations, which discretize # the equations in space to provide an ordinary differential equation that still diff --git a/docs/make.jl b/docs/make.jl index ee7d9f9c8e..237e4978ba 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -42,6 +42,8 @@ Literate.markdown(joinpath("docs", "literate", "src", "tut_rigid_body_fsi.jl"), joinpath("docs", "src", "tutorials")) Literate.markdown(joinpath("docs", "literate", "src", "tut_packing.jl"), joinpath("docs", "src", "tutorials")) +Literate.markdown(joinpath("docs", "literate", "src", "tut_2d_geometry.jl"), + joinpath("docs", "src", "tutorials")) copy_file("AUTHORS.md", "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)") @@ -94,6 +96,8 @@ makedocs(sitename="TrixiParticles.jl", "tut_rigid_body_fsi.md") ], "Preprocessing" => [ + "Setting up a 2D simulation from geometry files" => joinpath("tutorials", + "tut_2d_geometry.md"), "Particle packing tutorial" => joinpath("tutorials", "tut_packing.md") ] diff --git a/docs/src/development.md b/docs/src/development.md index 2e1e81fefe..ba237d4170 100644 --- a/docs/src/development.md +++ b/docs/src/development.md @@ -1,7 +1,6 @@ # [Development](@id development) - -## Preview of the documentation +## Preview of the Documentation To build the documentation, first instantiate the `docs` environment by running the following command from the TrixiParticles.jl root directory: @@ -23,30 +22,30 @@ python3 -m http.server -d docs/build ``` and open `localhost:8000` in your web browser. -## Release management +## Release Management To create a new release for TrixiParticles.jl, perform the following steps: -1) Make sure that all PRs and changes that you want to go into the release are merged to +1. Make sure that all PRs and changes that you want to go into the release are merged to `main` and that the latest commit on `main` has passed all CI tests. -2) Determine the currently released version of TrixiParticles.jl, e.g., on the +2. Determine the currently released version of TrixiParticles.jl, e.g., on the [release page](https://github.com/trixi-framework/TrixiParticles.jl/releases). For this manual, we will assume that the latest release was `v0.2.3`. -3) Decide on the next version number. We follow [semantic versioning](https://semver.org/), +3. Decide on the next version number. We follow [semantic versioning](https://semver.org/), thus each version is of the form `vX.Y.Z` where `X` is the major version, `Y` the minor version, and `Z` the patch version. In this manual, we assume that the major version is always `0`, thus the decision process on the new version is as follows: - * If the new release contains *breaking changes* (i.e., user code might not work as + - If the new release contains *breaking changes* (i.e., user code might not work as before without modifications), increase the *minor* version by one and set the *patch* version to zero. In our example, the new version should thus be `v0.3.0`. - * If the new release only contains minor modifications and/or bug fixes, the *minor* + - If the new release only contains minor modifications and/or bug fixes, the *minor* version is kept as-is and the *patch* version is increased by one. In our example, the new version should thus be `v0.2.4`. -4) Review and update the `NEWS.md` file to ensure all relevant changes, features, and bugfixes +4. Review and update the `NEWS.md` file to ensure all relevant changes, features, and bugfixes are documented for this release under the appropriate version header. -5) Edit the `version` string in the +5. Edit the `version` string in the [`Project.toml`](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/Project.toml) and set it to the new version. Push/merge this change to `main`. -6) Go to GitHub and add a comment to the commit that you would like to become the new +6. Go to GitHub and add a comment to the commit that you would like to become the new release (typically this will be the commit where you just updated the version). You can comment on a commit by going to the [commit overview](https://github.com/trixi-framework/TrixiParticles.jl/commits/main/) and clicking @@ -54,14 +53,14 @@ To create a new release for TrixiParticles.jl, perform the following steps: ``` @JuliaRegistrator register ``` -7) Wait for the magic to happen! Specifically, JuliaRegistrator will create a new PR to the +7. Wait for the magic to happen. Specifically, JuliaRegistrator will create a new PR to the Julia registry with the new release information. After a grace period of ~15 minutes, this PR will be merged automatically. A short while after, [TagBot](https://github.com/trixi-framework/TrixiParticles.jl/blob/main/.github/workflows/TagBot.yml) will create a new release of TrixiParticles.jl in our GitHub repository. -8) Once the new release has been created, the new version can be obtained through the Julia +8. Once the new release has been created, the new version can be obtained through the Julia package manager as usual. -9) To make sure people do not mistake the latest state of `main` as the latest release, we +9. To make sure people do not mistake the latest state of `main` as the latest release, we set the version in the `Project.toml` to a *development* version. The development version should be the latest released version, with the patch version incremented by one, and the `-dev` suffix added. For example, if you just released `v0.3.0`, the new development diff --git a/docs/src/general/neighborhood_search.md b/docs/src/general/neighborhood_search.md index 93019eaffc..f2e0d47d2d 100644 --- a/docs/src/general/neighborhood_search.md +++ b/docs/src/general/neighborhood_search.md @@ -8,7 +8,8 @@ different implementations. !!! note "Usage" To run a simulation with a neighborhood search implementation, pass a neighborhood - search template to the constructor of the [`Semidiscretization`](@ref). + search template to the constructor of the + [`Semidiscretization`](@ref TrixiParticles.Semidiscretization). A template is just an empty neighborhood search with search radius `0.0`. See [`copy_neighborhood_search`](@ref) and the examples below for more details. ```jldoctest semi_example; output=false, setup = :(using TrixiParticles; trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"), sol=nothing); system1 = fluid_system; system2 = boundary_system) diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index ce41637552..fae83fd48b 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -1,4 +1,5 @@ -# [Getting started](@id getting_started) +# [Getting Started](@id getting_started) + If you have not installed TrixiParticles.jl yet, please follow the instructions in [Installation](@ref installation). This page provides a short introduction. For a broader introduction, take a look at our [Tutorials](tutorial.md). @@ -28,14 +29,14 @@ This will open a new window with a 2D visualization of the final solution: For more information about visualization, see [Visualization](visualization.md). -## Running other Examples +## Running Other Examples You can find more predefined examples under [Examples](examples.md). Run them from the Julia REPL by replacing `subfolder` and `example_name`: ```julia julia> trixi_include(joinpath(examples_dir(), "subfolder", "example_name.jl")) ``` -## Modifying an example +## Modifying an Example You can pass keyword arguments to the function `trixi_include` to overwrite assignments in the file. With `trixi_include`, we can overwrite variables defined in the example file to run a different simulation without modifying the file itself. diff --git a/docs/src/gpu.md b/docs/src/gpu.md index a7bdd6c3af..4d23e5c206 100644 --- a/docs/src/gpu.md +++ b/docs/src/gpu.md @@ -23,7 +23,7 @@ FullGridCellList{PointNeighbors.DynamicVectorOfVectors{...}(...) ``` We then need to pass this cell list to the neighborhood search and the neighborhood search -to the [`Semidiscretization`](@ref). +to the [`Semidiscretization`](@ref TrixiParticles.Semidiscretization). ```jldoctest gpu; output=false semi = Semidiscretization(fluid_system, boundary_system, neighborhood_search=GridNeighborhoodSearch{2}(; cell_list)) @@ -67,7 +67,7 @@ All data is transferred to the GPU during initialization and all loops over part and their neighbors are 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). -## Run an existing example file on the GPU +## Run an Existing Example File on the GPU The example file `examples/fluid/dam_break_2d_gpu.jl` demonstrates how to run an existing example file on a GPU. @@ -102,7 +102,7 @@ trixi_include_changeprecision(Float32, coordinates_eltype=Float32) ``` -## [Single precision simulations](@id single_precision) +## [Single Precision Simulations](@id single_precision) All GPU-supported features can also be used with single precision, which is significantly faster on most GPUs and required for many Apple GPUs. diff --git a/docs/src/index.md b/docs/src/index.md index 1b16a8606f..8b0e9a699f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -28,18 +28,22 @@ Its main features include:
-
2D Dam Break
+
+
2D Dam Break
 -
Moving Wall
+
+
Moving Wall
-
Oscillating Beam
+
+
Oscillating Beam
 -
Dam Break with Elastic Plate
+
+
Dam Break with Elastic Plate
diff --git a/docs/src/install.md b/docs/src/install.md index 533fd4fc5f..be5945afe4 100644 --- a/docs/src/install.md +++ b/docs/src/install.md @@ -1,11 +1,11 @@ # [Installation](@id installation) -## Setting up Julia +## Setting Up Julia If you have not installed Julia yet, please [follow the instructions on the official website](https://julialang.org/downloads/). TrixiParticles.jl works with Julia v1.10 and newer. We recommend using the latest stable release of Julia. -## For users +## For Users TrixiParticles.jl is a registered Julia package. You can install TrixiParticles.jl, time integration sub-packages of [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) @@ -18,7 +18,7 @@ julia> using Pkg julia> Pkg.add(["TrixiParticles", "OrdinaryDiffEqLowStorageRK", "OrdinaryDiffEqSymplecticRK", "Plots"]) ``` -## [For developers](@id for-developers) +## [For Developers](@id for-developers) If you plan on editing TrixiParticles.jl itself, you can download TrixiParticles.jl to a local folder and use the code from the cloned directory: ```bash @@ -41,19 +41,19 @@ related packages (e.g., sub-packages of OrdinaryDiffEq.jl, see above) to the pro in the `run` folder and always have a reproducible environment at hand to share with others. -## Optional software/packages +## Optional Software/Packages - [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) -- A Julia package of ordinary differential equation solvers; examples in TrixiParticles.jl use sub-packages such as `OrdinaryDiffEqLowStorageRK` and `OrdinaryDiffEqSymplecticRK` -- [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia Plotting library that is used in some examples +- [Plots.jl](https://github.com/JuliaPlots/Plots.jl) -- Julia plotting library used in some examples - [PythonPlot.jl](https://github.com/JuliaPy/PythonPlot.jl) -- Plotting library that can be used instead of Plots.jl - [ParaView](https://www.paraview.org/) -- Visualization software for simulation results -## [Common issues](@id installation-issues) +## [Common Issues](@id installation-issues) If you followed the [installation instructions for developers](@ref for-developers) and run into package issues after pulling the latest version of TrixiParticles.jl, start Julia with the project in the `run` folder, ```bash - julia --project=run +julia --project=run ``` then update packages, resolve dependency conflicts, and install new dependencies: ```julia diff --git a/docs/src/preprocessing/preprocessing.md b/docs/src/preprocessing/preprocessing.md index d65fa462b2..beb56b3211 100644 --- a/docs/src/preprocessing/preprocessing.md +++ b/docs/src/preprocessing/preprocessing.md @@ -31,9 +31,9 @@ triangle = [125.0 375.0 250.0 125.0; 175.0 175.0 350.0 175.0] # Delete all edges but one -edge1 = deleteat!(TrixiParticles.Polygon(triangle), [2, 3]) -edge2 = deleteat!(TrixiParticles.Polygon(triangle), [1, 3]) -edge3 = deleteat!(TrixiParticles.Polygon(triangle), [1, 2]) +edge1 = delete_faces(TrixiParticles.Polygon(triangle), [2, 3]) +edge2 = delete_faces(TrixiParticles.Polygon(triangle), [1, 3]) +edge3 = delete_faces(TrixiParticles.Polygon(triangle), [1, 2]) algorithm = WindingNumberJacobson() @@ -269,6 +269,15 @@ For example: 0.0 1.0 ``` It is the user’s responsibility to ensure the points are ordered correctly. +For 2D `.asc` and `.dxf` files, `load_geometry` appends the first point by default +when it is not already repeated. This is only a convenience for complete, ordered +boundaries that omit the final duplicate point; it does not repair missing +segments, gaps, self-intersections, or incorrectly ordered points. Use +`load_geometry(file; close_curve=false)` for intentional open curves. Operations +that sample or classify a region, such as [`ComplexShape`](@ref), `intersect`, +and `setdiff`, require closed geometries. Boundary packing with +[`SignedDistanceField`](@ref) also requires a closed geometry, since it needs a +well-defined outside region. This format is easy to generate and inspect manually. ## DXF Format (.dxf) – recommended @@ -308,6 +317,11 @@ Modules = [TrixiParticles] Pages = [joinpath("preprocessing", "geometries", "triangle_mesh.jl")] ``` +```@docs +TrixiParticles.is_closed_geometry +delete_faces +``` + # [Particle Packing](@id particle_packing) To obtain a body-fitted and isotropic particle distribution, an initial configuration (see [Sampling of Geometries](@ref sampling_of_geometries)) is first generated. This configuration is then packed using a [`ParticlePackingSystem`](@ref) following the steps introduced in [Neher2026](@cite). @@ -325,7 +339,7 @@ The second step involves generating the SDF (see [`SignedDistanceField`](@ref)), The SDF is illustrated in Fig. 2, where the distances to the surface of the geometry are visualized as a color map. As shown, the SDF is computed only within a narrow band around the geometry’s surface, enabling a face-based neighborhood search (NHS) to be used exclusively during this step. In the third step, the initial configuration of the boundary particles is generated (orange particles in Fig. 3). -Boundary particles are created by copying the positions of SDF points located outside the geometry but within a predefined boundary thickness (see [`sample_boundary`](@ref)). +Boundary particles are created by copying the positions of SDF points located outside the geometry, starting at the offset implied by `place_on_shell` and ending at a predefined boundary thickness (see [`sample_boundary`](@ref)). In the fourth step, the initial configuration of the interior particles (green particles in Fig. 4) is generated using the hierarchical winding number approach (see [Hierarchical Winding](@ref hierarchical_winding)). After steps **1** through **4**, the initial configuration of both interior and boundary particles is obtained, as illustrated in Fig. 5. The interface of the geometry surface is not well resolved with the initial particle configuration. diff --git a/docs/src/refs.bib b/docs/src/refs.bib index c479e59210..11b15055e0 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -847,6 +847,17 @@ @Article{Valizadeh2015 publisher = {Elsevier BV}, } +@Article{Vergnaud2022, + author = {Vergnaud, A. and Oger, G. and Le Touz{\'e}, D. and DeLeffe, M. and Chiron, L.}, + title = {{C-CSF}: Accurate, robust and efficient surface tension and contact angle models for single-phase flows using {SPH}}, + journal = {Computer Methods in Applied Mechanics and Engineering}, + year = {2022}, + volume = {389}, + pages = {114292}, + doi = {10.1016/j.cma.2021.114292}, + publisher = {Elsevier BV}, +} + @Article{Wang2024, author = {Zhentong Wang and Bo Zhang and Oskar J. Haidn and Xiangyu Hu}, title = {A fourth-order kernel for improving numerical accuracy and stability in Eulerian SPH for fluids and total Lagrangian SPH for solids}, diff --git a/docs/src/systems/boundary.md b/docs/src/systems/boundary.md index cfab770268..414d069bef 100644 --- a/docs/src/systems/boundary.md +++ b/docs/src/systems/boundary.md @@ -25,7 +25,7 @@ dummy particles need to have a mass corresponding to the fluid's rest density, w "hydrodynamic mass", as opposed to mass corresponding to the material density of a [`TotalLagrangianSPHSystem`](@ref). -Here, `initial_density` and `hydrodynamic_mass` are vectors that contains the initial density +Here, `initial_density` and `hydrodynamic_mass` are vectors that contain the initial density and the hydrodynamic mass respectively for each boundary particle. Note that when used with [`SummationDensity`](@ref) (see below), this is only used to determine the element type and the number of boundary particles. @@ -37,14 +37,18 @@ This should be the same as for the adjacent fluid system with the largest smooth In the literature, this kind of boundary particles is referred to as "dummy particles" ([Adami et al., 2012](@cite Adami2012) and [Valizadeh & Monaghan, 2015](@cite Valizadeh2015)), -"frozen fluid particles" ([Akinci et al., 2012](@cite Akinci2012)) or "dynamic boundaries [Crespo et al., 2007](@cite Crespo2007). +"frozen fluid particles" ([Akinci et al., 2012](@cite Akinci2012)) or "dynamic boundaries" ([Crespo et al., 2007](@cite Crespo2007)). The key detail of this boundary condition and the only difference between the boundary models in these references is the way the density and pressure of boundary particles is computed. -Since boundary particles are treated like fluid particles, the force -on fluid particle ``a`` due to boundary particle ``b`` is given by +Since boundary particles are treated like fluid particles, their density and pressure enter +the pressure-acceleration operator selected by the interacting fluid system. For the +summation-density pressure operator, the pressure force on fluid particle ``a`` due to +boundary particle ``b`` is ```math -f_{ab} = m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_{r_a} W(\Vert r_a - r_b \Vert, h). +\bm{f}_{ab}^{p} += -m_a m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h). ``` The quantities to be defined here are the density ``\rho_b`` and pressure ``p_b`` of the boundary particle ``b``. @@ -99,13 +103,21 @@ where the sum is over all fluid particles, ``\rho_f`` and ``p_f`` denote the den ``` #### 2. [`BernoulliPressureExtrapolation`](@ref) -Identical to the pressure ``p_b `` calculated via [`AdamiPressureExtrapolation`](@ref), but it adds the dynamic pressure component of the Bernoulli equation: +Identical to the pressure ``p_b`` calculated via [`AdamiPressureExtrapolation`](@ref), +but with an added dynamic pressure term: ```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)} +p_b = \frac{\sum_f (p_f + p_{f,\mathrm{dyn}} + \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 bodies moving with a relative velocity to the fluid to prevent penetration. -This modification is original and not derived from any literature source. +with +```math +p_{f,\mathrm{dyn}} = \frac{1}{2} \, \text{factor} \, \rho_f +\left( +\frac{(\bm{v}_b - \bm{v}_f) \cdot \bm{r}_{bf}}{\Vert \bm{r}_{bf} \Vert} +\right)^2, +``` +where ``\bm{v}_f`` is the fluid velocity and ``\bm{v}_b`` is the boundary velocity. +This implementation-specific term raises the boundary pressure based on the normal +component of the relative boundary-fluid velocity and is not taken from a literature formula. ```@docs BernoulliPressureExtrapolation @@ -123,17 +135,17 @@ reference pressure (the corresponding pressure to the reference density by the s #### 6. [`PressureMirroring`](@ref) Instead of calculating density and pressure for each boundary particle, we modify the -momentum equation, +boundary pressure used in the pressure-acceleration operator. For the summation-density +pressure operator, this corresponds to modifying the pressure force ```math -\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab} +\bm{F}_a^{p} = -m_a \sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_b}{\rho_b^2} \right) \nabla_a W_{ab}, ``` -to replace the unknown density $\rho_b$ if $b$ is a boundary particle by the reference density -and the unknown pressure $p_b$ if $b$ is a boundary particle by the pressure $p_a$ of the -interacting fluid particle. -The momentum equation therefore becomes +to replace the unknown density ``\rho_b`` if ``b`` is a boundary particle by the reference density +and the unknown pressure ``p_b`` if ``b`` is a boundary particle by the pressure ``p_a`` of the +interacting fluid particle. The force therefore becomes ```math -\frac{\mathrm{d}v_a}{\mathrm{d}t} = -\sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af} --\sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab}, +\bm{F}_a^{p} = -m_a \sum_f m_f \left( \frac{p_a}{\rho_a^2} + \frac{p_f}{\rho_f^2} \right) \nabla_a W_{af} +-m_a \sum_b m_b \left( \frac{p_a}{\rho_a^2} + \frac{p_a}{\rho_0^2} \right) \nabla_a W_{ab}, ``` where the first sum is over all fluid particles and the second over all boundary particles. @@ -169,21 +181,25 @@ condition is applied. ## Repulsive Particles Boundaries modeled as boundary particles which exert forces on the fluid particles ([Monaghan, Kajtar, 2009](@cite Monaghan2009)). -The force on fluid particle ``a`` due to boundary particle ``b`` is given by +The force on fluid particle ``a`` due to boundary particle ``b`` is ```math -f_{ab} = m_a \left(\tilde{f}_{ab} - m_b \Pi_{ab} \nabla_{r_a} W(\Vert r_a - r_b \Vert, h)\right) +\bm{f}_{ab} = m_a \left(\tilde{\bm{f}}_{ab} - m_b \Pi_{ab} +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h)\right) ``` with ```math -\tilde{f}_{ab} = \frac{K}{\beta^{n-1}} \frac{r_{ab}}{\Vert r_{ab} \Vert (\Vert r_{ab} \Vert - d)} \Phi(\Vert r_{ab} \Vert, h) -\frac{2 m_b}{m_a + m_b}, +\tilde{\bm{f}}_{ab} = +\frac{K}{\beta^{n-1}} \frac{\bm{r}_{ab}} +{\Vert \bm{r}_{ab} \Vert (\Vert \bm{r}_{ab} \Vert - d)} +\Phi(\Vert \bm{r}_{ab} \Vert, h), ``` where ``m_a`` and ``m_b`` are the masses of fluid particle ``a`` and boundary particle ``b`` -respectively, ``r_{ab} = r_a - r_b`` is the difference of the coordinates of particles +respectively, ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` is the difference of the coordinates of particles ``a`` and ``b``, ``d`` denotes the boundary particle spacing and ``n`` denotes the number of dimensions (see [Monaghan & Kajtar, 2009](@cite Monaghan2009), Equation (3.1) and [Valizadeh & Monaghan, 2015](@cite Valizadeh2015)). -Note that the repulsive acceleration $\tilde{f}_{ab}$ does not depend on the masses of -the boundary particles. +The implemented repulsive acceleration ``\tilde{\bm{f}}_{ab}`` does not depend on the particle masses. +The denominator ``\Vert \bm{r}_{ab} \Vert - d`` is clipped from below by ``d/100`` in the +implementation to avoid the singularity at ``\Vert \bm{r}_{ab} \Vert = d``. Here, ``\Phi`` denotes the 1D Wendland C4 kernel, normalized to ``1.77`` for ``q=0`` ([Monaghan & Kajtar, 2009](@cite Monaghan2009), Section 4), with ``\Phi(r, h) = w(r/h)`` and ```math @@ -206,8 +222,8 @@ In [Monaghan & Kajtar (2009)](@cite Monaghan2009), a value of ``gD`` is used for where ``g`` is the gravitational acceleration and ``D`` is the depth of the fluid. The viscosity ``\Pi_{ab}`` is calculated according to the viscosity used in the -simulation, where the density of the boundary particle if needed is assumed to be -identical to the density of the fluid particle. +simulation. When a boundary density is needed, it is computed from the boundary +hydrodynamic mass and boundary particle spacing as ``m_b / d^n``. ### No-slip condition diff --git a/docs/src/systems/entropically_damped_sph.md b/docs/src/systems/entropically_damped_sph.md index 96acbad352..dd9dc8ccac 100644 --- a/docs/src/systems/entropically_damped_sph.md +++ b/docs/src/systems/entropically_damped_sph.md @@ -3,44 +3,57 @@ As opposed to the [weakly compressible SPH scheme](weakly_compressible_sph.md), which uses an equation of state, this scheme uses a pressure evolution equation to calculate the pressure ```math -\frac{\mathrm{d} p_a}{\mathrm{d}t} = - \rho c_s^2 \nabla \cdot v + \nu \nabla^2 p, +\frac{\mathrm{d} p_a}{\mathrm{d}t} = - \rho_a c_s^2 (\nabla \cdot v)_a + \nu_{\mathrm{EDAC}} (\nabla^2 p)_a, ``` which is derived by [Clausen (2013)](@cite Clausen2013). This equation is similar to the continuity equation (first term, see -[`ContinuityDensity`](@ref)), but also contains a pressure damping term (second term, similar to density diffusion +[`ContinuityDensity`](@ref)), but also contains a pressure damping term (second term, similar to density diffusion, see [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion)), which reduces acoustic pressure waves through an entropy-generation mechanism. -The pressure evolution is discretized with the SPH method by [Ramachandran (2019)](@cite Ramachandran2019) as following: +The pressure evolution is discretized with the SPH method by [Ramachandran (2019)](@cite Ramachandran2019) as follows: The first term is equivalent to the classical artificial compressible methods, which are commonly motivated by assuming the artificial equation of state ([`StateEquationCole`](@ref) with `exponent=1`) and is discretized as ```math -- \rho c_s^2 \nabla \cdot v = \sum_{b} m_b \frac{\rho_a}{\rho_b} c_s^2 v_{ab} \cdot \nabla_{r_a} W(\Vert r_a - r_b \Vert, h), +\left.- \rho c_s^2 \nabla \cdot v \right|_a += \sum_{b} m_b \frac{\rho_a}{\rho_b} c_s^2 v_{ab} \cdot \nabla_{r_a} W(\Vert r_a - r_b \Vert, h), ``` where ``\rho_a``, ``\rho_b``, ``r_a``, ``r_b``, denote the density and coordinates of particles ``a`` and ``b`` respectively, ``c_s`` is the speed of sound and ``v_{ab} = v_a - v_b`` is the difference in the velocity. The second term smooths the pressure through the introduction of entropy and is discretized as ```math -\nu \nabla^2 p = \frac{V_a^2 + V_b^2}{m_a} \tilde{\eta}_{ab} \frac{p_{ab}}{\Vert r_{ab}^2 \Vert + \eta h_{ab}^2} \nabla_{r_a} -W(\Vert r_a - r_b \Vert, h) \cdot r_{ab}, +\left.\nu_{\mathrm{EDAC}} \nabla^2 p \right|_a += \sum_b \frac{V_a^2 + V_b^2}{m_a}\, +\tilde{\eta}_{ab}\, +\frac{p_{ab}}{\Vert r_{ab} \Vert^2 + 0.01 h_{ab}^2}\, +\nabla_{r_a} W(\Vert r_a - r_b \Vert, h) \cdot r_{ab}, ``` -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively and ``p_{ab}= p_a -p_b`` is the difference in the pressure. +where ``V_a``, ``V_b`` denote the particle volumes, ``p_{ab}= p_a - p_b``, +``r_{ab} = r_a - r_b``, and ``h_{ab} = \frac{1}{2}(h_a + h_b)``. -The viscosity parameter ``\eta_a`` for a particle ``a`` is given as +The dynamic EDAC viscosity for particle ``a`` is ```math -\eta_a = \rho_a \frac{\alpha h c_s}{8}, +\eta_a = \rho_a \nu_{\mathrm{EDAC}}, ``` -where it is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5`` +with +```math +\nu_{\mathrm{EDAC}} = \frac{\alpha h c_s}{8}, +``` +and the harmonic mean +```math +\tilde{\eta}_{ab} = \frac{2 \eta_a \eta_b}{\eta_a + \eta_b}. +``` +It is found in the numerical experiments of [Ramachandran (2019)](@cite Ramachandran2019) that ``\alpha = 0.5`` is a good choice for a wide range of Reynolds numbers (0.0125 to 10000). !!! note - > The EDAC formulation keeps the density constant and this eliminates the need for the continuity equation - > or the use of a summation density to find the pressure. However, in SPH discretizations, ``m/\rho`` - > is typically used as a proxy for the particle volume. The density of the fluids can - > therefore be computed using the summation density approach. [Ramachandran2019](@cite) - + The EDAC formulation keeps the density constant and therefore eliminates the need for + the continuity equation or the use of a summation density to find the pressure. + However, in SPH discretizations, ``m/\rho`` is typically used as a proxy for the + particle volume. The density of the fluids can therefore still be computed using the + summation-density approach [Ramachandran2019](@cite). ```@autodocs Modules = [TrixiParticles] diff --git a/docs/src/systems/fluid.md b/docs/src/systems/fluid.md index c9b9860ab9..49d4349cf6 100644 --- a/docs/src/systems/fluid.md +++ b/docs/src/systems/fluid.md @@ -55,38 +55,36 @@ by Balsara ([Balsara1995](@cite)) or Morris ([Morris1997](@cite)). ##### Mathematical Formulation -The force exerted by particle ``b`` on particle ``a`` due to artificial viscosity is given by: +The acceleration contribution from particle ``b`` to particle ``a`` is ```math -F_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab} +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d}t}\right|_{ab}^{\text{AV}} = +\begin{cases} + m_b \frac{\alpha c \mu_{ab} + \beta \mu_{ab}^2}{\bar{\rho}_{ab}} + \nabla_a W_{ab}, & \text{if } \bm{v}_{ab} \cdot \bm{r}_{ab} < 0, \\ + 0, & \text{otherwise}. +\end{cases} ``` where: -- ``\Pi_{ab}`` is the artificial viscosity term defined as: - ```math - \Pi_{ab} = - \begin{cases} - -\frac{\alpha c \mu_{ab} + \beta \mu_{ab}^2}{\bar{\rho}_{ab}} & \text{if } \mathbf{v}_{ab} \cdot \mathbf{r}_{ab} < 0, \\ - 0 & \text{otherwise} - \end{cases} - ``` - ``\alpha`` and ``\beta`` are viscosity parameters, - ``c`` is the local speed of sound, - ``\bar{\rho}_{ab}`` is the arithmetic mean of the densities of particles ``a`` and ``b``. -The term ``\mu_{ab}`` is defined as: +The term ``\mu_{ab}`` is defined as ```math -\mu_{ab} = \frac{h \, v_{ab} \cdot r_{ab}}{\Vert r_{ab} \Vert^2 + \epsilon h^2}, +\mu_{ab} = \frac{h \, \bm{v}_{ab} \cdot \bm{r}_{ab}} + {\Vert \bm{r}_{ab} \Vert^2 + \epsilon h^2}, ``` with: - ``h`` being the smoothing length, - ``\epsilon`` a small parameter to prevent singularities, -- ``r_{ab} = r_a - r_b`` representing the difference of the coordinate vectors, -- ``v_{ab} = v_a - v_b`` representing the relative velocity between particles. +- ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` representing the difference of the coordinate vectors, +- ``\bm{v}_{ab} = \bm{v}_a - \bm{v}_b`` representing the relative velocity between particles. ##### Resolution Dependency and Effective Viscosity @@ -109,19 +107,21 @@ This results in a more realistic representation of flow dynamics in weakly compr ##### Mathematical Formulation -An additional force term ``\tilde{f}_{ab}`` is introduced to the pressure gradient force ``f_{ab}`` between particles ``a`` and ``b``: +An additional force term ``\tilde{\bm{F}}_{ab}`` is introduced in the momentum equation: ```math -\tilde{f}_{ab} = m_a m_b \frac{(\mu_a + \mu_b)\, r_{ab} \cdot \nabla W_{ab}}{\rho_a \rho_b (\Vert r_{ab} \Vert^2 + \epsilon h^2)}\, v_{ab}, +\tilde{\bm{F}}_{ab} = +m_a m_b \frac{(\mu_a + \mu_b)\, \bm{r}_{ab} \cdot \nabla_a W_{ab}} +{\rho_a \rho_b (\Vert \bm{r}_{ab} \Vert^2 + \epsilon h^2)}\, \bm{v}_{ab}, ``` where: -- ``\mu_a = \rho_a \nu`` and ``\mu_b = \rho_b \nu`` represent the dynamic viscosities of particles ``a``and ``b`` (with ``\nu`` being the kinematic viscosity), -- ``r_{ab} = r_a - r_b`` represents the difference of the coordinate vectors, -- ``v_{ab} = v_a - v_b`` represents the relative velocity between particles. +- ``\mu_a = \rho_a \nu`` and ``\mu_b = \rho_b \nu`` represent the dynamic viscosities of particles ``a`` and ``b`` (with ``\nu`` being the kinematic viscosity), +- ``\bm{r}_{ab} = \bm{r}_a - \bm{r}_b`` represents the difference of the coordinate vectors, +- ``\bm{v}_{ab} = \bm{v}_a - \bm{v}_b`` represents the relative velocity between particles, - `` h `` is the smoothing length, -- `` \nabla W_{ab} `` is the gradient of the smoothing kernel, +- `` \nabla_a W_{ab} `` is the gradient of the smoothing kernel, - `` \epsilon `` is a small parameter to prevent singularities. #### ViscosityAdami @@ -132,19 +132,24 @@ while minimizing compressibility effects. This results in accurate laminar flow ##### Mathematical Formulation -The viscous interaction is modeled through a shear force for incompressible flows: +The viscous interaction is modeled through the following pairwise force: ```math -f_{ab} = \sum_w \bar{\eta}_{ab} \left( V_a^2 + V_b^2 \right) \frac{v_{ab}}{||r_{ab}||^2 + \epsilon h_{ab}^2} \, (\nabla W_{ab} \cdot r_{ab}), +\bm{F}_{ab}^{\nu} = +\left( V_a^2 + V_b^2 \right)\, +\bar{\eta}_{ab}\, +\frac{\nabla_a W_{ab} \cdot \bm{r}_{ab}} +{\Vert \bm{r}_{ab} \Vert^2 + \epsilon h_{ab}^2}\, +\bm{v}_{ab}. ``` where: -- `` r_{ab} = r_a - r_b `` is the difference of the coordinate vectors, -- `` v_{ab} = v_a - v_b `` is their relative velocity, +- `` \bm{r}_{ab} = \bm{r}_a - \bm{r}_b `` is the difference of the coordinate vectors, +- `` \bm{v}_{ab} = \bm{v}_a - \bm{v}_b `` is their relative velocity, - `` V_a = m_a / \rho_a`` and `` V_b = m_b / \rho_b`` are the particle volumes, -- `` h_{ab} `` is the smoothing length, -- `` \nabla W_{ab} `` is the gradient of the smoothing kernel, +- `` h_{ab} = \frac{1}{2}(h_a + h_b) `` is the arithmetic mean of the smoothing lengths, +- `` \nabla_a W_{ab} `` is the gradient of the smoothing kernel, - `` \epsilon `` is a small parameter that prevents singularities (see [Ramachandran (2019)](@cite Ramachandran2019)). The inter-particle-averaged shear stress is defined as: @@ -225,10 +230,10 @@ The surface normal at a particle is derived from the color field, a scalar field to distinguish between different fluid phases or between fluid and air. The color field gradients point towards the interface, and the normalized gradient defines the surface normal direction. -The simplest SPH formulation for a surface normal, ``n_a`` is given as +In the literature, the unnormalized surface normal ``\bm{n}_a`` is commonly written as ```math -n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, +\bm{n}_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab}, ``` where: @@ -238,12 +243,14 @@ where: - ``\rho_b`` is the density of particle ``b``, - ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel ``W_{ab}`` with respect to particle ``a``. +For single-fluid surface-normal calculations, ``c_b = 1`` for neighboring fluid particles. + #### Normalization of surface normals The calculated normals are normalized to unit vectors: ```math -\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}. +\hat{\bm{n}}_a = \frac{\bm{n}_a}{\Vert \bm{n}_a \Vert}. ``` Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces. @@ -291,11 +298,12 @@ In the following table some values are shown for reference. The values marked wi ### [Akinci-based intra-particle force surface tension and wall adhesion model](@id akinci_ipf) -The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components: +The [Akinci](@cite Akinci2013) model divides surface tension into distinct force components, +which TrixiParticles.jl applies as acceleration contributions. -#### Cohesion force +#### Cohesion contribution -The cohesion force captures the attraction between particles at the fluid interface, creating the effect of surface tension. +The cohesion contribution captures the attraction between particles at the fluid interface, creating the effect of surface tension. It is defined by the distance between particles and the support radius ``h_c``, using a kernel-based formulation. **Key features:** @@ -303,10 +311,11 @@ It is defined by the distance between particles and the support radius ``h_c``, - Particles within half the support radius experience a repulsive force to prevent clustering. - Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion. -Mathematically: +In the acceleration form used by TrixiParticles.jl, the pairwise cohesion contribution is ```math -F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert}, +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d} t}\right|_{ab}^{\text{cohesion}} += -\sigma m_b C(r) \frac{\bm{r}}{\Vert \bm{r} \Vert}, ``` where ``C(r)``, the cohesion kernel, is defined as: @@ -315,29 +324,29 @@ where ``C(r)``, the cohesion kernel, is defined as: C(r)=\frac{32}{\pi h_c^9} \begin{cases} (h_c-r)^3 r^3, & \text{if } 2r > h_c, \\ -2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ +2(h_c-r)^3 r^3 - \frac{h_c^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ 0, & \text{otherwise.} \end{cases} ``` -#### Surface area minimization force - -The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. -It acts based on the difference in surface normals: +#### Surface area minimization contribution +The surface area minimization contribution models curvature reduction and acts on the +difference in surface normals. TrixiParticles.jl uses the local smoothing length: ```math -F_{\text{curvature}} = -\sigma (n_a - n_b), +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d} t}\right|_{ab}^{\text{curvature}} += -\sigma h (\bm{n}_a - \bm{n}_b), ``` +where ``\bm{n}_a`` and ``\bm{n}_b`` are the surface normals of the interacting particles. -where ``n_a`` and ``n_b`` are the surface normals of the interacting particles. +#### Wall adhesion contribution -#### Wall adhesion force - -This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. +This contribution models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. It uses a custom kernel with a peak at 0.75 times the support radius: ```math -F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert}, +\left.\frac{\mathrm{d}\bm{v}_a}{\mathrm{d} t}\right|_{ab}^{\text{adhesion}} += -\beta m_b A(r) \frac{\bm{r}}{\Vert \bm{r} \Vert}, ``` where ``A(r)`` is the adhesion kernel: @@ -358,15 +367,73 @@ The method described by [Morris](@cite Morris2000) estimates curvature by combin The computed curvature is then used to determine forces acting perpendicular to the interface. While this method provides accurate surface tension forces, it does not explicitly conserve momentum. -In the Morris model, surface tension is computed based on local interface curvature ``\kappa`` and the unit surface normal ``\hat{n}.`` -By estimating ``\hat{n}`` and ``\kappa`` at each particle near the interface, the surface tension force for particle a can be written as: +In the Morris model, surface tension is computed from local interface curvature ``\kappa``, the +unit surface normal ``\hat{\bm{n}}``, and the surface delta ``\delta_s``. The acceleration is a +particle-local source evaluated once per right-hand side evaluation: ```math -F_{\text{surface tension}} = - \sigma \frac{\kappa_a}{\rho_a}\hat{n}_a +\frac{\mathrm d\bm v_a}{\mathrm dt}\bigg|_\sigma += -\frac{\sigma}{\rho_a}\kappa_a\delta_{s,a}\hat{\bm n}_a. ``` -This formulation focuses directly on geometric properties of the interface, making it relatively straightforward to implement when a reliable interface detection -(e.g., a color function) is available. However, accurately estimating ``\kappa`` and ``n`` may require fine resolutions. +The factors have dimensions ``[\sigma]=kg/s^2``, ``[\kappa]=1/m``, +``[\delta_s]=1/m``, and ``[\rho]=kg/m^3``, giving acceleration in ``m/s^2``. This formulation +does not explicitly conserve momentum, and accurately estimating curvature still requires +adequate resolution. + +#### Smooth colorfield interface activity + +With `ColorfieldSurfaceNormal`, Morris CSF uses a C1 interface activity ``\lambda_a``. Let +``h_c`` be the compact-support radius, ``\gamma_a=h_c\Vert\bm g_a\Vert``, ``\epsilon_n`` be +`interface_threshold`, and ``\alpha`` be `interface_taper_start` (default `0.8`). With + +```math +S(x)=\begin{cases} +0,&x\le0,\\ +3x^2-2x^3,&0 0`, the +support activity is ``\lambda_{q,a}=1-S((q_a-\tau)/\Delta q)``, where ``\Delta q`` is +`support_taper_width` (default `0.025`). The final activity and one-phase surface delta are + +```math +\lambda_a=\lambda_{g,a}\lambda_{q,a},\qquad +\delta_{s,a}=2\Vert\bm g_a\Vert\lambda_a. +``` + +Setting `ideal_density_threshold=0` disables support filtering. For Morris CSF with +`ColorfieldSurfaceNormal`, this keyword represents a continuous fraction of complete kernel +support instead of an integer neighbor-count fraction. Dummy boundary particles complete +``q_a`` near walls without carrying capillary stress. These controls do not alter the separate +C-CSF geometry described below. + +[`CorrectedCSFSurfaceNormal`](@ref) selects the corrected continuous-surface-force (C-CSF) +interface geometry of [Vergnaud et al.](@cite Vergnaud2022) for [`SurfaceTensionMorris`](@ref). It +computes the outward normal from the renormalized gradient of the minimum eigenvalue of the +first-order kernel moment. Curvature uses a renormalized divergence and the published thin-jet +angular filter; the surface delta uses the published Shepard correction. + +```julia +surface_tension = SurfaceTensionMorris(surface_tension_coefficient=0.072) +surface_normal_method = CorrectedCSFSurfaceNormal() +``` + +This explicit opt-in supports one fluid system and free-surface geometry. Boundary-integral and +contact-angle terms are not included. --- @@ -379,28 +446,35 @@ where accumulated numerical error can be significant. #### Stress tensor formulation -The surface tension force can be seen as a divergence of a stress tensor ``S`` +The surface tension force can be written as the divergence of a stress tensor ``\bm{S}``: ```math -F_{\text{surface tension}} = \nabla \cdot S, +\bm{F}_{a}^{\sigma} = m_a \nabla \cdot \bm{S}, ``` -with ``S`` defined as +with ```math -S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}), +\bm{S} = \sigma \delta_s (I - \hat{\bm{n}} \otimes \hat{\bm{n}}). ``` with: - ``\delta_s``: Surface delta function, -- ``\hat{n}``: Unit normal vector, +- ``\hat{\bm{n}}``: Unit normal vector, - ``I``: Identity matrix. This divergence can be computed numerically in the SPH framework as ```math -\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab} +\bm{F}_{a}^{\sigma} += m_a \sum_b \frac{m_b}{\rho_a \rho_b} (\bm{S}_a + \bm{S}_b) \nabla_a W_{ab}. +``` + +TrixiParticles.jl stores ``\sigma`` outside the tensor and uses the stabilized tensor +```math +\bm{S}_a^{\text{impl}} += \delta_{s,a} (I - \hat{\bm{n}}_a \otimes \hat{\bm{n}}_a) - \delta_{s,\max} I, ``` #### Advantages and limitations diff --git a/docs/src/systems/implicit_incompressible_sph.md b/docs/src/systems/implicit_incompressible_sph.md index 0927aa8821..deccbfe4d4 100644 --- a/docs/src/systems/implicit_incompressible_sph.md +++ b/docs/src/systems/implicit_incompressible_sph.md @@ -28,14 +28,15 @@ difference yields The divergence in the right-hand side is discretized with the SPH discretization for particle ``i`` as ```math --\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \nabla W_{ij}, +-\frac{1}{\rho_i} \sum_j m_j \bm{v}_{ij} \cdot \nabla W_{ij}, ``` where ``\bm{v}_{ij} = \bm{v}_i - \bm{v}_j``. Together, the following discretized version of the continuity equation for a particle ``i`` is achieved: ```math -\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_j m_j \bm{v}_{ij}(t+\Delta t) \nabla W_{ij}. +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} += \sum_j m_j \bm{v}_{ij}(t+\Delta t) \cdot \nabla W_{ij}. ``` Note that the linear system is only solved for fluid particles, so ``i`` always represents @@ -50,12 +51,12 @@ Using the semi-implicit Euler method, we can obtain the velocity in the next tim ``` where ``\bm{F}_i^{\text{adv}}`` denotes all non-pressure forces such as gravity, viscosity, surface -tension and more, while ``\bm{F}_i^p``denotes the unknown pressure forces, which we +tension and more, while ``\bm{F}_i^p`` denotes the unknown pressure forces, which we want to solve for. Note that the IISPH is an incompressible method, which means that the density of the -fluid remain constant over time. By assuming a fixed reference density ``\rho_0`` for all -fluid particle over the whole time of the simulation, the density value at the next time +fluid remains constant over time. By assuming a fixed reference density ``\rho_0`` for all +fluid particles over the whole simulation, the density value at the next time step ``\rho_i(t + \Delta t)`` also has to be this rest density. So ``\rho_0`` can be plugged in for ``\rho_i(t + \Delta t)`` in the equation above. @@ -72,7 +73,7 @@ Using this predicted velocity and the continuity equation, a predicted density c in a similar way as ```math -\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}} \nabla W_{ij}(t). +\rho_i^{\text{adv}}(t + \Delta t)= \rho_i(t) + \Delta t \sum_j m_j \bm{v}_{ij}^{\text{adv}}(t+\Delta t) \cdot \nabla W_{ij}(t). ``` To achieve the rest density, the unknown pressure forces must counteract the compression @@ -81,7 +82,7 @@ the predicted density and the reference density. Therefore, the following equation needs to be fulfilled: ```math -\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. +\Delta t ^2 \sum_j m_j \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_j^p(t)}{m_j} \right) \cdot \nabla W_{ij}(t) = \rho_0 - \rho_i^{\text{adv}}. ``` This expression is derived by substituting the reference density ``\rho_0`` for @@ -130,10 +131,10 @@ The pressure acceleration is given by: The ``d_{ii}p_i`` value describes the displacement of particle ``i`` because of the particle ``i`` and ``d_{ij}p_j`` describes the influence from the neighboring particles ``j``. -Using this new values the linear system can be rewritten as +Using these values, the linear system can be rewritten as ```math -\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \nabla W_{ij}, +\rho_0 - \rho_i^{\text{adv}} = \sum_j m_j \left( d_{ii}p_i + \sum_k d_{ik}p_k - d_{jj}p_j - \sum_k d_{jk}p_k \right) \cdot \nabla W_{ij}, ``` where the first sum over ``k`` loops over all neighbor particles of ``i`` and @@ -150,21 +151,21 @@ To separate this sum, it can be written as With this separation, the equation for the linear system can again be rewritten as ```math -\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \nabla W_{ij}. +\rho_0 - \rho_i^{\text{adv}} = p_i \sum_j m_j ( d_{ii} - d_{ji}) \cdot \nabla W_{ij} + \sum_j m_j \left ( \sum_k d_{ik} p_k - d_{jj} p_j - \sum_{k \neq i} d_{jk}p_k \right) \cdot \nabla W_{ij}. ``` In this formulation all coefficients that are getting multiplied with the pressure value ``p_i`` are separated from the other. The diagonal elements ``a_{ii}`` can therefore be defined as: ```math -a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij}. +a_{ii} = \sum_j m_j ( d_{ii} - d_{ji}) \cdot \nabla W_{ij}. ``` The remaining part of the equation represents the influence of the other pressure values ``p_j``. -​Hence, the final relaxed Jacobi iteration takes the form: +Hence, the final relaxed Jacobi iteration takes the form: ```math -p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \nabla W_{ij} \right). +p_i^{l+1} = (1 - \omega) p_i^{l} + \omega \frac{1}{a_{ii}} \left( \rho_0 -\rho_i^{\text{adv}} - \sum_j m_j \left( \sum_k d_{ik} p_k^l - d_{jj} p_j^l - \sum_{k \neq i} d_{jk} p_k^l \right) \cdot \nabla W_{ij} \right). ``` Because interactions are local, limited to particles within the kernel support defined by @@ -209,7 +210,7 @@ as only isolated or almost isolated particles are affected. ## Boundary Handling The previously introduced formulation did not distinguish between fluid and boundary -particles. To account boundary interactions correctly, a few modifications to the previous +particles. To account for boundary interactions correctly, a few modifications to the previous equations are required. First, the discretized form of the continuity equation must be adapted for the case in which @@ -220,7 +221,7 @@ neighboring fluid particles (indexed by ``f``) and neighboring boundary particle The updated discretized continuity equation becomes: ```math -\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_f m_f \bm{v}_{if}(t+\Delta t) \nabla W_{if} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \nabla W_{ib}. +\frac{\rho_i(t + \Delta t) - \rho_i(t)}{\Delta t} = \sum_f m_f \bm{v}_{if}(t+\Delta t) \cdot \nabla W_{if} + \sum_b m_b \bm{v}_{ib}(t+\Delta t) \cdot \nabla W_{ib}. ``` Since boundary particles have zero velocity, the difference between the fluid @@ -229,13 +230,13 @@ particle's velocity ``\bm{v}_{ib}(t+\Delta t) = \bm{v}_{i}(t+\Delta t)``. Accordingly, the predicted density ``\rho^{\text{adv}}`` becomes: ```math -\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_f m_f \bm{v}_{if}^{\text{adv}} \nabla W_{if}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \nabla W_{ib}(t). +\rho_i^{\text{adv}} = \rho_i (t) + \Delta t \sum_f m_f \bm{v}_{if}^{\text{adv}} \cdot \nabla W_{if}(t) + \Delta t \sum_b m_b \bm{v}_{i}^{\text{adv}} \cdot \nabla W_{ib}(t). ``` This leads to the following updated formulation of the linear system: ```math -\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \nabla W_{if} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \cdot \nabla W_{if} + \Delta t^2 \sum_b m_b \frac{\bm{F}_i^p(t)}{m_i} \cdot \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. ``` Note that, since boundary particles are fixed, the force ``F_b^p`` is zero and does not appear @@ -244,12 +245,12 @@ in this equation. The pressure force acting on a fluid particle is computed as: ```math -\bm{F}_i^p(t) = -\sum_f m_f \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_f(t)}{\rho_f^2(t)} \right) \nabla W_{if}(t) - \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). +\bm{F}_i^p(t) = -m_i \sum_f m_f \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_f(t)}{\rho_f^2(t)} \right) \nabla W_{if}(t) - m_i \sum_b m_b \left( \frac{p_i(t)}{\rho_i^2(t)} + \frac{p_b(t)}{\rho_b^2(t)} \right) \nabla W_{ib}(t). ``` This also leads to an updated version of the equation for the diagonal elements: ```math -a_{ii} = \sum_j m_j ( d_{ii} - d_{ji})\nabla W_{ij} + \sum_b m_b (-d_{bi}) \nabla W_{ib}. +a_{ii} = \sum_f m_f ( d_{ii} - d_{fi}) \cdot \nabla W_{if} + \sum_b m_b d_{ii} \cdot \nabla W_{ib}. ``` From this point forward, the computation of the coefficients required for the Jacobi scheme @@ -263,11 +264,12 @@ When using pressure mirroring, the pressure value ``p_b`` of a boundary particle above is defined to be equal to the pressure of the corresponding fluid particle ``p_i``. In other words, the boundary particle "mirrors" the pressure of the fluid particle interacting with it. As a result, the coefficient that describes the influence of a particle's own -pressure value ``p_i`` ​must also include contributions from boundary particles. Therefore, -the equation for calculating the coefficient ``d_{ii}`` must be adjusted as follows: +pressure value ``p_i`` must include a doubled contribution from each boundary particle. +Therefore, ``d_{ii}`` becomes ```math -d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. +d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} + - 2\Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. ``` The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: @@ -275,30 +277,29 @@ The corresponding relaxed Jacobi iteration for pressure mirroring then becomes: ```math \begin{align*} p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} - - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \nabla W_{if} \right. \\ -& \quad - \left. \sum_b m_b \sum_f d_{if} p_f^l \nabla W_{ib} \right). + - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \cdot \nabla W_{if} \right. \\ +& \quad - \left. \sum_b m_b \left( \sum_f d_{if} p_f^l \right) \cdot \nabla W_{ib} \right). \end{align*} ``` ### Pressure Zeroing If pressure zeroing is used instead, the pressure value of a boundary particle ``p_b`` -​is assumed to be zero. Consequently, boundary particles do not contribute to the pressure -forces acting on fluid particles. -In this case, the computation of the coefficient ``d_{ii}`` remains unchanged and is given by: +is assumed to be zero. In the linear system, this removes the boundary pressure unknowns, +but boundary particles still contribute through the ``p_i/\rho_i^2`` part of the pressure +acceleration. Therefore ``d_{ii}`` is ```math -d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if}. +d_{ii} = -\Delta t^2 \sum_f \frac{m_f}{\rho_i^2} \nabla W_{if} + - \Delta t^2 \sum_b \frac{m_b}{\rho_i^2} \nabla W_{ib}. ``` -The equation for the relaxed Jacobi iteration remains the same as in the pressure mirroring -approach. However, the contribution from boundary particles vanishes due to their zero -pressure: +The corresponding relaxed Jacobi iteration reads ```math \begin{align*} p_i^{l+1} = (1 - \omega) p_i^l + \omega \frac{1}{a_{ii}} &\left( \rho_0 - \rho_i^{\text{adv}} - - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \nabla W_{if} \right. \\ -& \quad - \left. \sum_b m_b \sum_j d_{if} p_f^l \nabla W_{ib} \right). + - \sum_f m_f \left( \sum_k d_{ik} p_k^l - d_{ff}p_f^l - \sum_{k \neq i} d_{fk} p_k^l \right) \cdot \nabla W_{if} \right. \\ +& \quad - \left. \sum_b m_b \left( \sum_f d_{if} p_f^l \right) \cdot \nabla W_{ib} \right). \end{align*} ``` @@ -307,8 +308,8 @@ The density calculators [`AdamiPressureExtrapolation`](@ref) and [`BernoulliPres can also be used with IISPH. When using one of these pressure extrapolation methods the calculation of the PPE is exactly the same as when using pressure zeroing. -So within the linear systems the pressure values are equal to zero (``p_b=0``) and therefore -are not considered in the calculations. Only in the pressure acceleration, the extrapolated +So within the linear system the boundary pressures are treated as zero (``p_b=0``), exactly +as in pressure zeroing. Only in the pressure acceleration, the extrapolated pressure values are used for the boundary particles. For more information on these two methods, refer to the docs for the [boundary models](@ref boundary_models). @@ -318,7 +319,7 @@ The [`PressureBoundaries`](@ref) density calculator was introduced by only be used with IISPH. In the standard IISPH method the PPE is solved only for fluid particles. The pressure values for the boundary particles are then approximated, for example by using -pressure mirroing. +pressure mirroring. With `PressureBoundaries`, however, the linear system is extended to include the boundary particles as well. This means that the pressure values of both the fluid and the boundary particles are computed directly by solving the PPE. @@ -333,22 +334,22 @@ also solved as part of the linear system. This leads to the following condition for the boundary particles ``b``: ```math -\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \nabla W_{if} + \Delta t^2 \sum_b m_b \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_b^p(t)}{m_b} \right) \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_f^p(t)}{m_f} \right) \cdot \nabla W_{if} + \Delta t^2 \sum_b m_b \left( \frac{\bm{F}_i^p(t)}{m_i} - \frac{\bm{F}_b^p(t)}{m_b} \right) \cdot \nabla W_{ib} = \rho_0 - \rho_i^{\text{adv}}. ``` -Note that in this case ``i`` is a boundary particle,``f`` are its fluid neighbors, and ``b`` +Note that in this case ``i`` is a boundary particle, ``f`` are its fluid neighbors, and ``b`` its boundary neighbors. Since the pressure force for boundary particles is zero (as mentioned before), and because in this case ``i`` and ``b`` are both boundary particles, the PPE simplifies to ```math -\Delta t^2 \sum_f m_f - \frac{\bm{F}_f^p(t)}{m_f} \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. +-\Delta t^2 \sum_f m_f \frac{\bm{F}_f^p(t)}{m_f} \cdot \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. ``` If we substitute the definition of the pressure force from above, we obtain ```math -\Delta t^2 \sum_f m_f \left( \sum_k m_k \left( \frac{p_f(t)}{\rho_j^2(t)} + \frac{p_k(t)}{\rho_k^2(t)} \right) \nabla W_{fk}\right) \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. +\Delta t^2 \sum_f m_f \left( \sum_k m_k \left( \frac{p_f(t)}{\rho_f^2(t)} + \frac{p_k(t)}{\rho_k^2(t)} \right) \nabla W_{fk}\right) \cdot \nabla W_{if} = \rho_0 - \rho_i^{\text{adv}}. ``` where ``k`` represents all neighboring particles (fluid and boundary) of fluid particle ``f``. @@ -358,23 +359,23 @@ indirectly as neighbors of fluid particles), their ``d_{ii}`` values are zero. T diagonal elements simplify to ```math -a_{ii} = \sum_f \left( -d_{fi}\right) \nabla W_{if}. +a_{ii} = - \sum_f m_f d_{fi} \cdot \nabla W_{if}. ``` The off-diagonal term ``\sum_{j \neq i} a_{ij} p_j`` in the relaxed Jacobi iteration for boundary particles takes the following form ```math -\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( d_{ff} - \sum_{f_j} d_{ff_j}p_{f_j}\right) \nabla W_{if}. +\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( - d_{ff} p_f - \sum_{k \neq i} d_{fk}p_k \right) \cdot \nabla W_{if}. ``` But not only the addition of the boundary particles to the linear system changes when using pressure boundaries, also the PPE for the fluid particles is changing slightly. Since the boundary particles have now their own pressure values, ``p_b`` is no longer -eliminated (as in pressure zeroing or pressure extrapolation with ``p_b=0``) nor simply +eliminated (as in pressure zeroing or pressure extrapolation with ``p_b=0``) nor simply replaced by the fluid particle's pressure (as in pressure mirroring with ``p_b=p_i``). -Instead, ``p_b``remains part of the PPE. +Instead, ``p_b`` remains part of the PPE. This has no effect on the diagonal elements ``a_{ii}``, which are identical to those of the other density calculators. The ``d_{ii}`` values are also computed in the same way as for @@ -383,11 +384,11 @@ However, the off-diagonal term ``\sum_{j \neq i} a_{ij} p_j`` in the relaxed Ja changes slightly. For pressure boundaries it takes the form ```math - \sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( \sum_k d_{ik} p_k - d_{ff} p_f - \sum_{k \neq f} d_{fk} p_k \right) \nabla W_{if} \\ - + \sum_b m_b \left( \sum_k d_{ik} p_k \right) \nabla W_{ib}, +\sum_{j \neq i} a_{ij} p_j = \sum_f m_f \left( \sum_k d_{ik} p_k - d_{ff} p_f - \sum_{k \neq i} d_{fk} p_k \right) \cdot \nabla W_{if} ++ \sum_b m_b \left( \sum_k d_{ik} p_k \right) \cdot \nabla W_{ib}, ``` -where ``k``represents all neighboring particles of ``i`` (both fluid and boundary). +where ``k`` represents all neighboring particles of ``i`` (both fluid and boundary). ```@docs PressureBoundaries -``` \ No newline at end of file +``` diff --git a/docs/src/systems/total_lagrangian_sph.md b/docs/src/systems/total_lagrangian_sph.md index b6dc4e6f04..2965aa9713 100644 --- a/docs/src/systems/total_lagrangian_sph.md +++ b/docs/src/systems/total_lagrangian_sph.md @@ -9,7 +9,7 @@ The governing equations with respect to the initial configuration are given by: \frac{\mathrm{D}\bm{v}}{\mathrm{D}t} = \frac{1}{\rho_0} \nabla_0 \cdot \bm{P} + \bm{g}, ``` where the zero subscript denotes a derivative with respect to the initial configuration -and $\bm{P}$ is the first Piola-Kirchhoff (PK1) stress tensor. +and ``\bm{P}`` is the first Piola-Kirchhoff (PK1) stress tensor. The discretized version of this equation is given by [O’Connor & Rogers (2021)](@cite OConnor2021): ```math @@ -21,17 +21,17 @@ with the correction matrix (see also [`GradientCorrection`](@ref)) ```math \bm{L}_{0a} := \left( -\sum_{b} \frac{m_{0b}}{\rho_{0b}} \nabla_{0a} W(\bm{X}_{ab}) \bm{X}_{ab}^T \right)^{-1} \in \R^{d \times d}. ``` -The subscripts $a$ and $b$ denote quantities of particle $a$ and $b$, respectively. +The subscripts ``a`` and ``b`` denote quantities of particles ``a`` and ``b``, respectively. The zero subscript on quantities denotes that the quantity is to be measured in the initial configuration. -The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$, -the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$. +The difference in the initial coordinates is denoted by ``\bm{X}_{ab} = \bm{X}_a - \bm{X}_b``, +the difference in the current coordinates is denoted by ``\bm{x}_{ab} = \bm{x}_a - \bm{x}_b``. -For the computation of the PK1 stress tensor, the deformation gradient $\bm{F}$ is computed per particle as +For the computation of the PK1 stress tensor, the deformation gradient ``\bm{F}`` is computed per particle as ```math \bm{F}_a = \sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ba} (\bm{L}_{0a}\nabla_{0a} W(\bm{X}_{ab}))^T \\ \qquad = -\left(\sum_b \frac{m_{0b}}{\rho_{0b}} \bm{x}_{ab} (\nabla_{0a} W(\bm{X}_{ab}))^T \right) \bm{L}_{0a}^T ``` -with $1 \leq i,j \leq d$. +with ``1 \leq i,j \leq d``. From the deformation gradient, the Green-Lagrange strain ```math \bm{E} = \frac{1}{2}(\bm{F}^T\bm{F} - \bm{I}) @@ -53,9 +53,9 @@ and ```math \lambda = \frac{E\nu}{(1 + \nu)(1 - 2\nu)} ``` -are the Lamé coefficients, where $E$ is the Young's modulus and $\nu$ is the Poisson ratio. +are the Lamé coefficients, where ``E`` is the Young's modulus and ``\nu`` is the Poisson ratio. -The term $\bm{f}_a^{PF}$ is an optional penalty force. See e.g. [`PenaltyForceGanzenmueller`](@ref). +The term ``\bm{f}_a^{PF}`` is an optional penalty force. See e.g. [`PenaltyForceGanzenmueller`](@ref). ```@autodocs Modules = [TrixiParticles] @@ -69,26 +69,26 @@ This is caused by the stiffness matrix having zero eigenvalues (so-called hourgl The name "hourglass modes" comes from the fact that elements can deform into an hourglass shape. Similar effects can occur in SPH as well. -Particles can change positions without changing the SPH approximation of the deformation gradient $\bm{F}$, +Particles can change positions without changing the SPH approximation of the deformation gradient ``\bm{F}``, thus, without causing an increase of energy. To ensure regular particle positions, we can apply similar correction forces as are used in FEM. -[Ganzenmüller (2015)](@cite Ganzenmueller2015) introduced a so-called hourglass correction force or penalty force $f^{PF}$, +[Ganzenmüller (2015)](@cite Ganzenmueller2015) introduced a so-called hourglass correction force or penalty force ``f^{PF}``, which is given by ```math \bm{f}_a^{PF} = \frac{1}{2} \alpha \sum_b \frac{m_{0a} m_{0b} W_{0ab}}{\rho_{0a}\rho_{0b} |\bm{X}_{ab}|^2} \left( E \delta_{ab}^a + E \delta_{ba}^b \right) \frac{\bm{x}_{ab}}{|\bm{x}_{ab}|} ``` -The subscripts $a$ and $b$ denote quantities of particle $a$ and $b$, respectively. +The subscripts ``a`` and ``b`` denote quantities of particles ``a`` and ``b``, respectively. The zero subscript on quantities denotes that the quantity is to be measured in the initial configuration. -The difference in the initial coordinates is denoted by $\bm{X}_{ab} = \bm{X}_a - \bm{X}_b$, -the difference in the current coordinates is denoted by $\bm{x}_{ab} = \bm{x}_a - \bm{x}_b$. -Note that [Ganzenmüller (2015)](@cite Ganzenmueller2015) has a flipped sign here because they define $\bm{x}_{ab}$ the other way around. +The difference in the initial coordinates is denoted by ``\bm{X}_{ab} = \bm{X}_a - \bm{X}_b``, +the difference in the current coordinates is denoted by ``\bm{x}_{ab} = \bm{x}_a - \bm{x}_b``. +Note that [Ganzenmüller (2015)](@cite Ganzenmueller2015) has a flipped sign here because they define ``\bm{x}_{ab}`` the other way around. This correction force is based on the potential energy density of a Hookean material. -Thus, $E$ is the Young's modulus and $\alpha$ is a dimensionless coefficient that controls +Thus, ``E`` is the Young's modulus and ``\alpha`` is a dimensionless coefficient that controls the amplitude of hourglass correction. -The separation vector $\delta_{ab}^a$ indicates the change of distance which the particle separation should attain +The separation vector ``\delta_{ab}^a`` indicates the change of distance which the particle separation should attain in order to minimize the error and is given by ```math \delta_{ab}^a = \frac{\bm{\epsilon}_{ab}^a \cdot \bm{x_{ab}}}{|\bm{x}_{ab}|}, @@ -106,16 +106,16 @@ Pages = [joinpath("schemes", "structure", "total_lagrangian_sph", "penalty_force ## Viscosity Another technique that is used to correct the hourglass instability is artificial viscosity. -Hereby, a viscosity term designed for fluids (see [Viscosity](@ref viscosity_sph)) is applied. -First, the force ``f_{ab}^{\text{fluid}}`` exerted by particle ``b`` on particle ``a`` +Here, a viscosity term designed for fluids (see [Viscosity](@ref viscosity_sph)) is applied. +First, the force ``\bm{F}_{ab}^{\text{fluid}}`` exerted by particle ``b`` on particle ``a`` due to artificial viscosity is computed as if both particles were fluid particles (see [Viscosity](@ref viscosity_sph) for the relevant equations). Then, according to [Lin et al. (2015)](@cite Lin2015), this force can be applied to TLSPH with the following conversion: ```math -f_{ab}^{\text{AV}} = \det(F_a) F_a^{-1} f_{ab}^{\text{fluid}}, +\bm{F}_{ab}^{\text{AV}} = \det(\bm{F}_a) \bm{F}_a^{-T} \bm{F}_{ab}^{\text{fluid}}, ``` -where ``F_a`` is the deformation gradient at particle ``a``. +where ``\bm{F}_a`` is the deformation gradient at particle ``a``. We found that artificial viscosity is not effective at correcting the incorrect particle positions due to hourglass modes. diff --git a/docs/src/systems/weakly_compressible_sph.md b/docs/src/systems/weakly_compressible_sph.md index 2c69f4b3c9..86edf2590c 100644 --- a/docs/src/systems/weakly_compressible_sph.md +++ b/docs/src/systems/weakly_compressible_sph.md @@ -52,22 +52,24 @@ pressure field. It is highly recommended to use density diffusion when using WCS ### Formulation All density diffusion terms extend the continuity equation (see [`ContinuityDensity`](@ref)) -by an additional term +by an additional term: ```math -\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla W_{ab} - + \delta h c \sum_{b} V_b \psi_{ab} \cdot \nabla W_{ab}, +\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = + \sum_{b} m_b \frac{\rho_a}{\rho_b} v_{ab} \cdot \nabla W_{ab} + + \delta c \sum_{b} \bar{h}_{ab} V_b \psi_{ab} \cdot \nabla W_{ab}, ``` -where ``V_b = m_b / \rho_b`` is the volume of particle ``b`` and ``\psi_{ab}`` depends on +where ``\bar{h}_{ab} = \frac{1}{2}(h_a + h_b)`` is the averaged smoothing length, +``V_b = m_b / \rho_b`` is the volume of particle ``b`` and ``\psi_{ab}`` depends on the density diffusion method (see [`AbstractDensityDiffusion`](@ref TrixiParticles.AbstractDensityDiffusion) for available terms). Also, ``\rho_a`` denotes the density of particle ``a`` and ``r_{ab} = r_a - r_b`` is the difference of the coordinates, ``v_{ab} = v_a - v_b`` of the velocities of particles -``a`` and ``b``. +``a`` and ``b``. For fixed smoothing length, ``\bar{h}_{ab} = h``. ### Numerical Results All density diffusion terms remove numerical noise in the pressure field and produce more -accurate results than weakly commpressible SPH without density diffusion. +accurate results than weakly compressible SPH without density diffusion. This can be demonstrated with dam break examples in 2D and 3D. Here, ``δ = 0.1`` has been used for all terms. Note that, due to added stability, the adaptive time integration method that was used here @@ -129,11 +131,14 @@ in such simulations. ### Mathematical formulation We use the following formulation by [Sun et al. (2018)](@cite Sun2018). -After each time step, a correction term ``\delta r_a`` is added to the position ``r_a`` +The relation ``\text{CFL} \cdot \text{Ma} = \Delta t \, v_\text{max} / h`` +is stated there on page 29, immediately above Equation 9, and gives the +dimensional form below. +After each time step, a correction term ``\delta \bm{r}_a`` is added to the position ``\bm{r}_a`` of particle ``a``, which is given by ```math -\delta r_a = -4 \Delta t \, v_\text{max} h - \sum_b \left( 1 + R \left( \frac{W_{ab}}{W(\Delta x_a)} \right)^n \right) \nabla W_{ab} +\delta \bm{r}_a = -4 \Delta t \, v_\text{max} h + \sum_b \left( 1 + R \left( \frac{W_{ab}}{W(\Delta x_a)} \right)^n \right) \nabla_a W_{ab} \frac{m_b}{\rho_a + \rho_b}, ``` where: @@ -141,14 +146,26 @@ where: - ``v_\text{max}`` is the maximum velocity over all particles, - ``h`` is the smoothing length, - ``R`` and ``n`` are constants, which are set to ``0.2`` and ``4`` respectively, -- ``W(\Delta x_a)`` is the smoothing kernel of the particle size of particle ``a``, - which can be interpreted as the target particle spacing that we want to achieve. -- ``\nabla W_{ab}`` is the gradient of the smoothing kernel, +- ``\Delta x_a`` is the target particle spacing associated with particle ``a``, +- ``W(\Delta x_a)`` is the smoothing kernel evaluated at that target particle spacing, +- ``\nabla_a W_{ab}`` is the gradient of the smoothing kernel with respect to particle ``a``, - ``m_b`` is the mass of particle ``b``, - ``\rho_a, \rho_b`` is the density of particles ``a`` and ``b``, respectively. -Note that we replaced ``\text{CFL} \cdot \text{Ma}`` by ``\Delta t \cdot v_\text{max} / h``, -as explained in [Sun2018](@cite Sun2018) on page 29, right above Equation 9. +TrixiParticles.jl applies this correction through a shifting velocity +```math +\delta \bm{r}_a = \Delta t \, \delta \bm{v}_a, +``` +with +```math +\delta \bm{v}_a = - v_* \frac{(2h)^2}{2\Delta x} + \sum_b \left( 1 + \frac{2}{10} \left( \frac{W_{ab}}{W(\Delta x)} \right)^4 \right) + \frac{m_b}{\rho_a + \rho_b} \nabla_a W_{ab}. +``` +Here, ``v_*`` is the velocity scale configured by the shifting technique. It is either +``v_\text{factor}\max_a \Vert \bm{v}_a \Vert`` when `v_max_factor` is used, or +``v_\text{factor} c`` when `sound_speed_factor` is used. +The constants are fixed to ``R = 0.2`` and ``n = 4``. The ``\delta``-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as ``\delta^+``-SPH. @@ -178,18 +195,29 @@ is a constant background pressure field. The tilde in the second term of the right-hand side indicates that the material derivative has an advection part. -The discretized form of the last term is +In the literature, the discretized form of the last term is ```math -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab}, ``` -where ``V_a``, ``V_b`` denote the volume of particles ``a`` and ``b`` respectively. +where ``V_a`` and ``V_b`` denote the particle volumes of particles ``a`` and ``b`` respectively. Note that although in the continuous case ``\nabla p_{\text{background}} = 0``, -the discretization is not 0th-order consistent for **non**-uniform particle distribution, -which means that there is a non-vanishing contribution only when particles are disordered. -That also means that ``p_{\text{background}}`` occurs as pre-factor to correct -the trajectory of a particle resulting in uniform pressure distributions. -Suggested is a background pressure which is in the order of the reference pressure, -but it can be chosen arbitrarily large when the time-step criterion is adjusted. +the discretization is not 0th-order consistent for non-uniform particle distributions. +This means that a non-vanishing contribution appears only when the particles are disordered, +so ``p_{\text{background}}`` acts as a prefactor that regularizes the trajectories and promotes +more uniform particle distributions. + +TrixiParticles.jl evaluates this term with the selected pressure-acceleration operator. +For the default [`ContinuityDensity`](@ref) pressure acceleration and the CFL estimate +used by [Adami et al. (2013)](@cite Adami2013), +```math +\Delta t \leq \frac{1}{4} \frac{h}{c_s}, +``` +used as an equality, this gives +```math +\delta \bm{v}_a = - \frac{p_{\text{background}}}{8} \frac{h}{c_s} +\sum_b \frac{2m_b}{\rho_a \rho_b} \nabla_a W_{ab}, +``` +where ``h`` is the smoothing length and ``c_s`` is the speed of sound. The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is @@ -200,18 +228,31 @@ where the tensor ``\bm{A} = \rho v\left(\tilde{v}-v\right)^T`` is a consequence of the modified advection velocity and can be interpreted as the convection of momentum with the relative velocity ``\tilde{v}-v``. -The discretized form of the momentum equation for a particle ``a`` reads as +The discretized form of the momentum equation for a particle ``a`` reads ```math -\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. +\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} += \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) +\left[ -\tilde{p}_{ab} \nabla_a W_{ab} ++ \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right]. ``` Here, ``\tilde{p}_{ab}`` is the density-weighted pressure ```math \tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b}, ``` -with the density ``\rho_a``, ``\rho_b`` and the pressure ``p_a``, ``p_b`` of particles ``a`` -and ``b``, respectively. ``\bm{A}_a`` and ``\bm{A}_b`` are the convection tensors -for particle ``a`` and ``b``, respectively, and are given, e.g., for particle ``a``, -as ``\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T``. +with ``\rho_a``, ``\rho_b`` and ``p_a``, ``p_b`` denoting the densities and pressures +of particles ``a`` and ``b``, respectively. + +TrixiParticles.jl evaluates this additional term with the selected pressure-acceleration +operator. For the default [`ContinuityDensity`](@ref) pressure acceleration, this gives +```math +\left.\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t}\right|_{\bm{A}} += - \sum_b \frac{m_b}{\rho_a \rho_b} +\left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab}. +``` +Here, for example, +```math +\bm{A}_a = \rho_a \bm{v}_a \left(\tilde{\bm{v}}_a - \bm{v}_a\right)^T. +``` To apply the TVF, use the keyword argument `shifting_technique` in the constructor of a system that supports it. @@ -247,13 +288,13 @@ Only the combination of PST and TIC is able to produce physical results. The force that particle ``a`` experiences from particle ``b`` due to pressure is given by ```math -f_{ab} = -m_a m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla W_{ab} +\bm{f}_{ab} = -m_a m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla_a W_{ab} ``` for the WCSPH method with [`ContinuityDensity`](@ref). -The TIC formulation changes this force to +The TIC formulation changes this term to ```math -f_{ab} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla W_{ab}. +\bm{f}_{ab}^{\mathrm{TIC}} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla_a W_{ab}. ``` Note that this formulation is asymmetric and sacrifices conservation of linear and angular momentum. diff --git a/docs/src/time_integration.md b/docs/src/time_integration.md index f8a2842c13..58d9a47908 100644 --- a/docs/src/time_integration.md +++ b/docs/src/time_integration.md @@ -107,7 +107,7 @@ half step for ``v``, yielding u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, \operatorname{drift}(v^0, u^0, t^0), \\ v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\ v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\ -u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t). +u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1/2}, t^0 + \Delta t). \end{align*} ``` This scheme is implemented in `OrdinaryDiffEqSymplecticRK` as `LeapfrogDriftKickDrift` and yields @@ -133,7 +133,7 @@ v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, \operatorname{kick}(v^0, u^0, t^0), \\ \rho^{1/2} &= \rho^0 + \frac{1}{2} \Delta t\, R(v^0, u^0, t^0), \\ v^1 &= v^0 + \Delta t\, \operatorname{kick} \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t \right), \\ \rho^1 &= \rho^0 \frac{2 - \varepsilon^{1/2}}{2 + \varepsilon^{1/2}}, \\ -u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1}, t^0 + \Delta t), +u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, \operatorname{drift}(v^{1}, u^{1/2}, t^0 + \Delta t), \end{align*} ``` where diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 4225faf486..8a27d363b2 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -1,12 +1,16 @@ # Tutorials +Choose a tutorial based on the task in front of you. + > New to TrixiParticles.jl? Start with [Setting up your simulation from scratch](tutorials/tut_setup.md). ## Recommended Path 1. [Setting up your simulation from scratch](tutorials/tut_setup.md): learn the structure of a simulation file and run a complete WCSPH example. 2. [Modifying or extending components of TrixiParticles.jl within a simulation file](tutorials/tut_custom_kernel.md): replace selected parts of an existing setup without cloning the package. -3. [Particle packing tutorial](tutorials/tut_packing.md): build a body-fitted particle configuration for complex geometries. +3. [Setting up a 2D simulation from geometry files](tutorials/tut_2d_geometry.md): load 2D geometry files, turn them into filled wall regions, and combine them with standard 2D fluid blocks. +4. [Particle packing tutorial](tutorials/tut_packing.md): build a body-fitted particle configuration for complex geometries. +5. [Fluid-structure interaction with rigid bodies](tutorials/tut_rigid_body_fsi.md): simulate objects moving in a fluid. ## Tutorials @@ -38,6 +42,19 @@ directly in the file you run. - Focus: `trixi_include`, custom kernels, rapid iteration - Choose this if: you want to prototype changes without cloning and modifying the package +### [Setting up a 2D simulation from geometry files](tutorials/tut_2d_geometry.md) + +```@raw html +2D pipe and coastline geometries converted to wall and fluid particles +``` + +Load 2D geometry files, fill them with particles using `ComplexShape`, and build genuine 2D setups such as a curved pipe and a coastline dam break. + +- Focus: `load_geometry`, `ComplexShape`, `setdiff`, 2D `Polygon`s +- Choose this if: you want a true 2D setup from line-based geometry data + ### [Particle packing tutorial](tutorials/tut_packing.md) ```@raw html diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 9dda94bfba..ac58596792 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -1,6 +1,6 @@ # Visualization -## Export VTK files +## Export VTK Files You can export particle data as VTK files by using the [`SolutionSavingCallback`](@ref). All [predefined examples](examples.md) already use this callback to export VTK files to the `out` directory relative to the current working directory. @@ -18,6 +18,7 @@ Follow these steps to view the exported VTK files in ParaView: 5. Hold the left mouse button to move the solution around. You will now see the following: + ![image](https://github.com/user-attachments/assets/383d323a-3020-4232-9dc3-682b0afe8653) It is useful to make the dot size dependent on the actual particle size. @@ -28,9 +29,10 @@ Then, in the Properties panel (bottom left), adjust the following settings: 3. Activate "Scale by Array" and select "`particle_spacing`" in "Gaussian Scale Array". 4. Deactivate "Use Scale Function". 5. Set the "Gaussian Radius" to "`0.5`". + ![image](https://github.com/user-attachments/assets/194d9a09-5937-4ee4-b229-07078afe3ff0) -#### Visualization with Macro +#### Visualization with a Macro To simplify visualization of particle data in ParaView, you can use a macro. It reduces the manual steps from the previous section to a single click. Install the macro as follows. @@ -44,7 +46,6 @@ Install the macro as follows. 6. Click on the macro name in the **Macros** menu (or toolbar, if pinned) to run it. 7. The Point Gaussian representation with `particle_spacing` scaling will be applied automatically. - --- #### Macro Code @@ -78,11 +79,12 @@ sourceDisplay.UseScaleFunction = 0 sourceDisplay.GaussianRadius = 0.5 ``` -#### Show results +#### Show Results To view the result variables, first make sure that "fluid_1.pvd" is highlighted in the "Pipeline Browser", then select a variable in the variable-selection combo box (see the image below). For example, choose "density". To view the time evolution, press the play button (also shown below). + ![image](https://github.com/user-attachments/assets/10dcf7eb-5808-4d4d-9db8-4beb25b5e51a) ## API diff --git a/examples/fluid/dam_break_2d.jl b/examples/fluid/dam_break_2d.jl index 11ee987599..22969a427c 100644 --- a/examples/fluid/dam_break_2d.jl +++ b/examples/fluid/dam_break_2d.jl @@ -80,12 +80,8 @@ viscosity_wall = nothing # viscosity_wall = viscosity_fluid # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, - correction=nothing, - reference_particle_spacing=0, viscosity=viscosity_wall, clip_negative_pressure=true) diff --git a/examples/fluid/dam_break_3d.jl b/examples/fluid/dam_break_3d.jl index ae9b580960..0e5438ac53 100644 --- a/examples/fluid/dam_break_3d.jl +++ b/examples/fluid/dam_break_3d.jl @@ -57,10 +57,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothi boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fluid/falling_water_column_2d.jl b/examples/fluid/falling_water_column_2d.jl index 694c13ab76..c9c4b00fa5 100644 --- a/examples/fluid/falling_water_column_2d.jl +++ b/examples/fluid/falling_water_column_2d.jl @@ -55,10 +55,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothi boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl index 1f0016ac6d..657f496a48 100644 --- a/examples/fluid/falling_water_spheres_2d.jl +++ b/examples/fluid/falling_water_spheres_2d.jl @@ -82,12 +82,11 @@ boundary_density_calculator = AdamiPressureExtrapolation() wall_viscosity = nu # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; + fluid_system=sphere_surface_tension, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length; state_equation, viscosity=ViscosityAdami(nu=wall_viscosity), - reference_particle_spacing=fluid_particle_spacing, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model; diff --git a/examples/fluid/hydrostatic_water_column_2d.jl b/examples/fluid/hydrostatic_water_column_2d.jl index 58b381c0c0..93e6665c8f 100644 --- a/examples/fluid/hydrostatic_water_column_2d.jl +++ b/examples/fluid/hydrostatic_water_column_2d.jl @@ -60,10 +60,8 @@ boundary_density_calculator = AdamiPressureExtrapolation() # This is to set wall viscosity with `trixi_include` viscosity_wall = nothing -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, viscosity=viscosity_wall) boundary_system = WallBoundarySystem(tank.boundary, boundary_model, prescribed_motion=nothing) diff --git a/examples/fluid/lid_driven_cavity_2d.jl b/examples/fluid/lid_driven_cavity_2d.jl index c1aa4ae2a8..747caf2763 100644 --- a/examples/fluid/lid_driven_cavity_2d.jl +++ b/examples/fluid/lid_driven_cavity_2d.jl @@ -84,20 +84,14 @@ is_moving(t) = true lid_movement = PrescribedMotion(lid_movement_function, is_moving) -boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density, - cavity.boundary.mass, - AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length; - viscosity, state_equation) - -boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass, - AdamiPressureExtrapolation(), - smoothing_kernel, smoothing_length; - viscosity, state_equation) - -boundary_system_cavity = WallBoundarySystem(cavity.boundary, boundary_model_cavity) - -boundary_system_lid = WallBoundarySystem(lid, boundary_model_lid, +cavity_boundary_model = BoundaryModelDummyParticles(cavity.boundary; + fluid_system=fluid_system, + viscosity=viscosity) +boundary_system_cavity = WallBoundarySystem(cavity.boundary, cavity_boundary_model) + +lid_boundary_model = BoundaryModelDummyParticles(lid; fluid_system=fluid_system, + viscosity=viscosity) +boundary_system_lid = WallBoundarySystem(lid, lid_boundary_model, prescribed_motion=lid_movement) # ========================================================================================== diff --git a/examples/fluid/moving_wall_2d.jl b/examples/fluid/moving_wall_2d.jl index 6e0b0f9dfe..253386981c 100644 --- a/examples/fluid/moving_wall_2d.jl +++ b/examples/fluid/moving_wall_2d.jl @@ -59,10 +59,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothi # ========================================================================================== # ==== Boundary boundary_density_calculator = AdamiPressureExtrapolation() -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, - boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation) +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, + boundary_density_calculator) boundary_system = WallBoundarySystem(tank.boundary, boundary_model, prescribed_motion=boundary_movement) diff --git a/examples/fluid/periodic_array_of_cylinders_2d.jl b/examples/fluid/periodic_array_of_cylinders_2d.jl index f2899b6bfd..ec1344eafc 100644 --- a/examples/fluid/periodic_array_of_cylinders_2d.jl +++ b/examples/fluid/periodic_array_of_cylinders_2d.jl @@ -71,10 +71,8 @@ fluid_system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_le # ========================================================================================== # ==== Boundary -boundary_model = BoundaryModelDummyParticles(boundary.density, boundary.mass, - AdamiPressureExtrapolation(), smoothing_kernel, - smoothing_length; - viscosity=ViscosityAdami(; nu), state_equation) +boundary_model = BoundaryModelDummyParticles(boundary; fluid_system=fluid_system, + viscosity=ViscosityAdami(; nu)) boundary_system = WallBoundarySystem(boundary, boundary_model) diff --git a/examples/fluid/periodic_channel_2d.jl b/examples/fluid/periodic_channel_2d.jl index 9635f48a78..097e9b5f16 100644 --- a/examples/fluid/periodic_channel_2d.jl +++ b/examples/fluid/periodic_channel_2d.jl @@ -59,10 +59,8 @@ viscosity_wall = nothing # Activate to switch to no-slip walls #viscosity_wall = ViscosityAdami(nu=0.0025 * smoothing_length * sound_speed / 8) -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, viscosity=viscosity_wall) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fluid/pipe_flow_2d.jl b/examples/fluid/pipe_flow_2d.jl index 6022b980e4..1539d9c0f0 100644 --- a/examples/fluid/pipe_flow_2d.jl +++ b/examples/fluid/pipe_flow_2d.jl @@ -136,16 +136,13 @@ outflow = BoundaryZone(; boundary_face=face_out, face_normal=(-flow_direction), initial_condition=outlet.fluid, boundary_type=boundary_type_out) open_boundary = OpenBoundarySystem(inflow, outflow; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles) + boundary_model=open_boundary_model) # ========================================================================================== # ==== Boundary wall = union(pipe.boundary, inlet.boundary, outlet.boundary) viscosity_boundary = viscosity -boundary_model = BoundaryModelDummyParticles(wall.density, wall.mass, - AdamiPressureExtrapolation(), smoothing_kernel, - smoothing_length; state_equation, +boundary_model = BoundaryModelDummyParticles(wall; fluid_system=fluid_system, viscosity=viscosity_boundary) boundary_system = WallBoundarySystem(wall, boundary_model) diff --git a/examples/fluid/poiseuille_flow_2d.jl b/examples/fluid/poiseuille_flow_2d.jl index ba36a2f471..175571a8f2 100644 --- a/examples/fluid/poiseuille_flow_2d.jl +++ b/examples/fluid/poiseuille_flow_2d.jl @@ -145,16 +145,14 @@ outlet_boundary_zone = BoundaryZone(; boundary_face=outlet_face, open_boundary = OpenBoundarySystem(inlet_boundary_zone, outlet_boundary_zone; fluid_system, boundary_model=open_boundary_model, - calculate_flow_rate=true, - buffer_size=n_buffer_particles) + calculate_flow_rate=true) # ========================================================================================== # ==== Boundary wall_boundary = union(channel.boundary) -boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass, - AdamiPressureExtrapolation(), smoothing_kernel, - smoothing_length; state_equation, viscosity) +boundary_model = BoundaryModelDummyParticles(wall_boundary; fluid_system=fluid_system, + viscosity) boundary_system = WallBoundarySystem(wall_boundary, boundary_model) diff --git a/examples/fluid/poiseuille_flow_3d.jl b/examples/fluid/poiseuille_flow_3d.jl index 4459df0648..7d76cc8ffb 100644 --- a/examples/fluid/poiseuille_flow_3d.jl +++ b/examples/fluid/poiseuille_flow_3d.jl @@ -162,14 +162,12 @@ outlet_zone = BoundaryZone(; boundary_face=outlet_face, boundary_type=outlet_boundary_type) open_boundary = OpenBoundarySystem(inlet_zone, outlet_zone; fluid_system, - boundary_model=open_boundary_model, - buffer_size=n_buffer_particles) + boundary_model=open_boundary_model) # ========================================================================================== # ==== Boundary -boundary_model = BoundaryModelDummyParticles(wall_boundary.density, wall_boundary.mass, - AdamiPressureExtrapolation(), smoothing_kernel, - smoothing_length; state_equation, viscosity) +boundary_model = BoundaryModelDummyParticles(wall_boundary; fluid_system=fluid_system, + viscosity) boundary_system = WallBoundarySystem(wall_boundary, boundary_model) diff --git a/examples/fsi/dam_break_gate_2d.jl b/examples/fsi/dam_break_gate_2d.jl index 49d72f2aee..9cc9e66271 100644 --- a/examples/fsi/dam_break_gate_2d.jl +++ b/examples/fsi/dam_break_gate_2d.jl @@ -120,20 +120,15 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothi boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model_tank = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +tank_boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, clip_negative_pressure=true) - -boundary_model_gate = BoundaryModelDummyParticles(gate.density, gate.mass, +gate_boundary_model = BoundaryModelDummyParticles(gate; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, clip_negative_pressure=true) -boundary_system_tank = WallBoundarySystem(tank.boundary, boundary_model_tank) -boundary_system_gate = WallBoundarySystem(gate, boundary_model_gate, +boundary_system_tank = WallBoundarySystem(tank.boundary, tank_boundary_model) +boundary_system_gate = WallBoundarySystem(gate, gate_boundary_model, prescribed_motion=gate_movement) # ========================================================================================== diff --git a/examples/fsi/dam_break_plate_2d.jl b/examples/fsi/dam_break_plate_2d.jl index 00edc5c216..2bfd4f609e 100644 --- a/examples/fsi/dam_break_plate_2d.jl +++ b/examples/fsi/dam_break_plate_2d.jl @@ -94,10 +94,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; smoothing_kernel, smoothi boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length; - state_equation, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/falling_rigid_spheres_2d.jl b/examples/fsi/falling_rigid_spheres_2d.jl index 7eaa2b4dcb..59c05ce448 100644 --- a/examples/fsi/falling_rigid_spheres_2d.jl +++ b/examples/fsi/falling_rigid_spheres_2d.jl @@ -70,10 +70,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length; - state_equation, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/falling_rotating_rigid_squares_2d.jl b/examples/fsi/falling_rotating_rigid_squares_2d.jl index 4a14905cf0..0ab46fa9fd 100644 --- a/examples/fsi/falling_rotating_rigid_squares_2d.jl +++ b/examples/fsi/falling_rotating_rigid_squares_2d.jl @@ -86,10 +86,9 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; boundary_density_calculator = AdamiPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length; - state_equation, clip_negative_pressure=true) + clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl index 035f9ad3da..136161624f 100644 --- a/examples/fsi/falling_spheres_2d.jl +++ b/examples/fsi/falling_spheres_2d.jl @@ -74,10 +74,8 @@ fluid_system = WeaklyCompressibleSPHSystem(tank.fluid; boundary_density_calculator = BernoulliPressureExtrapolation() # Clip negative boundary pressure values to avoid sticking artifacts at the boundary. -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - fluid_smoothing_kernel, fluid_smoothing_length; - state_equation, clip_negative_pressure=true) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) diff --git a/examples/fsi/hydrostatic_water_column_2d.jl b/examples/fsi/hydrostatic_water_column_2d.jl index 9fd7e15a36..3815a28d68 100644 --- a/examples/fsi/hydrostatic_water_column_2d.jl +++ b/examples/fsi/hydrostatic_water_column_2d.jl @@ -113,9 +113,8 @@ else damping_coefficient=0.05)) end -boundary_model = BoundaryModelDummyParticles(tank.boundary.density, tank.boundary.mass, +boundary_model = BoundaryModelDummyParticles(tank.boundary; fluid_system=fluid_system, boundary_density_calculator, - smoothing_kernel, smoothing_length_fluid; state_equation) boundary_system = WallBoundarySystem(tank.boundary, boundary_model) boundary_model_structure = BoundaryModelDummyParticles(hydrodynamic_densities, diff --git a/examples/preprocessing/complex_shape_2d.jl b/examples/preprocessing/complex_shape_2d.jl index fa2b762988..62de84ad49 100644 --- a/examples/preprocessing/complex_shape_2d.jl +++ b/examples/preprocessing/complex_shape_2d.jl @@ -7,8 +7,8 @@ # 3. Utilize the Winding Number algorithm to determine if points are inside or outside. # 4. Visualize the sampled particles and the winding number field. # -# The example uses an "inverted_open_curve" geometry, where standard inside/outside -# definitions might be ambiguous without a robust point-in-polygon test like winding numbers. +# The example uses a polygonal star geometry, where standard inside/outside +# definitions benefit from a robust point-in-polygon test like winding numbers. # ========================================================================================== using TrixiParticles @@ -16,7 +16,7 @@ using Plots particle_spacing = 0.05 -filename = "inverted_open_curve" +filename = "star" file = joinpath("examples", "preprocessing", "data", filename * ".asc") geometry = load_geometry(file) @@ -24,7 +24,6 @@ geometry = load_geometry(file) trixi2vtk(geometry) point_in_geometry_algorithm = WindingNumberJacobson(; geometry, - winding_number_factor=0.4, hierarchical_winding=true) # Returns `InitialCondition` diff --git a/examples/preprocessing/data/coastline_profile_2d.asc b/examples/preprocessing/data/coastline_profile_2d.asc new file mode 100644 index 0000000000..c250c87be4 --- /dev/null +++ b/examples/preprocessing/data/coastline_profile_2d.asc @@ -0,0 +1,20 @@ +# ASCII +0.15 -0.12 0 +2.68 -0.12 0 +2.68 1.08 0 +2.62 0.66 0 +2.53 0.52 0 +2.42 0.40 0 +2.30 0.42 0 +2.18 0.33 0 +2.05 0.24 0 +1.92 0.26 0 +1.78 0.18 0 +1.62 0.11 0 +1.46 0.14 0 +1.28 0.06 0 +1.05 0.02 0 +0.82 0.05 0 +0.55 0.03 0 +0.15 0.03 0 +0.15 -0.12 0 diff --git a/examples/preprocessing/data/curved_pipe_channel_2d.asc b/examples/preprocessing/data/curved_pipe_channel_2d.asc new file mode 100644 index 0000000000..76a1d7b8fe --- /dev/null +++ b/examples/preprocessing/data/curved_pipe_channel_2d.asc @@ -0,0 +1,20 @@ +# ASCII +0.00 0.12 0 +0.60 0.12 0 +0.72423 0.13646 0 +0.84000 0.18431 0 +0.93941 0.26059 0 +1.01569 0.36000 0 +1.06354 0.47577 0 +1.08 0.60 0 +1.08 1.20 0 +0.72 1.20 0 +0.72 0.60 0 +0.71591 0.56894 0 +0.70392 0.54000 0 +0.68485 0.51515 0 +0.66000 0.49608 0 +0.63106 0.48409 0 +0.60 0.48 0 +0.00 0.48 0 +0.00 0.12 0 diff --git a/examples/preprocessing/data/curved_pipe_outer_2d.asc b/examples/preprocessing/data/curved_pipe_outer_2d.asc new file mode 100644 index 0000000000..1b26afac92 --- /dev/null +++ b/examples/preprocessing/data/curved_pipe_outer_2d.asc @@ -0,0 +1,14 @@ +# ASCII +0.00 0.00 0 +0.60 0.00 0 +0.75529 0.02044 0 +0.90000 0.08038 0 +1.02426 0.17574 0 +1.11962 0.30000 0 +1.17956 0.44471 0 +1.20 0.60 0 +1.20 1.20 0 +0.60 1.20 0 +0.60 0.60 0 +0.00 0.60 0 +0.00 0.00 0 diff --git a/examples/preprocessing/packing_2d.jl b/examples/preprocessing/packing_2d.jl index fcf0d40c64..7055a2dd98 100644 --- a/examples/preprocessing/packing_2d.jl +++ b/examples/preprocessing/packing_2d.jl @@ -71,7 +71,8 @@ packing_system = ParticlePackingSystem(shape_sampled; smoothing_length, boundary_system = ParticlePackingSystem(boundary_sampled; smoothing_length, is_boundary=true, signed_distance_field, - place_on_shell, boundary_compress_factor=0.8, + place_on_shell, boundary_thickness, + boundary_compress_factor=0.8, background_pressure) # ========================================================================================== diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl index 7101be8e49..d66c22d898 100644 --- a/src/TrixiParticles.jl +++ b/src/TrixiParticles.jl @@ -16,12 +16,13 @@ using ForwardDiff: ForwardDiff using GPUArraysCore: AbstractGPUArray using JSON: JSON using KernelAbstractions: KernelAbstractions, @kernel, @index -using LinearAlgebra: norm, normalize, cross, dot, I, tr, inv, pinv, det +using LinearAlgebra: norm, normalize, cross, dot, I, tr, inv, pinv, det, eigvals, + Symmetric using Polyester: Polyester, @batch using Printf: @printf, @sprintf using ReadVTK: ReadVTK using RecipesBase: RecipesBase, @series -using Random: seed! +using Random: MersenneTwister using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, DynamicalODEProblem, derivative_discontinuity!, get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!, add_tstop! @@ -102,7 +103,7 @@ export trixi2vtk, vtk2trixi export RectangularTank, RectangularShape, SphereShape, ComplexShape export ParticlePackingSystem, SignedDistanceField export WindingNumberHormann, WindingNumberJacobson -export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, +export VoxelSphere, RoundSphere, reset_wall!, extrude_geometry, load_geometry, delete_faces, sample_boundary, planar_geometry_to_face export SourceTermDamping export ShepardKernelCorrection, KernelCorrection, AkinciFreeSurfaceCorrection, @@ -114,7 +115,7 @@ export interpolate_line, interpolate_points, interpolate_plane_3d, interpolate_p interpolate_plane_2d_vtk export SurfaceTensionAkinci, CohesionForceAkinci, SurfaceTensionMorris, SurfaceTensionMomentumMorris -export ColorfieldSurfaceNormal +export ColorfieldSurfaceNormal, CorrectedCSFSurfaceNormal export SymplecticPositionVerlet export coordinates_eltype diff --git a/src/callbacks/density_reinit.jl b/src/callbacks/density_reinit.jl index 430f3b5644..14f850dec9 100644 --- a/src/callbacks/density_reinit.jl +++ b/src/callbacks/density_reinit.jl @@ -32,7 +32,8 @@ end Callback to reinitialize the density field when using [`ContinuityDensity`](@ref) [Panizzo2007](@cite). -Pass `system` and the [`Semidiscretization`](@ref) containing it. The callback stores +Pass `system` and the [`Semidiscretization`](@ref TrixiParticles.Semidiscretization) +containing it. The callback stores the system index and uses the corresponding system from the integrator semidiscretization at runtime, which remains valid if [`semidiscretize`](@ref) replaces systems internally. diff --git a/src/callbacks/mechanical_work_calculator.jl b/src/callbacks/mechanical_work_calculator.jl index f20672af09..07895673fd 100644 --- a/src/callbacks/mechanical_work_calculator.jl +++ b/src/callbacks/mechanical_work_calculator.jl @@ -28,7 +28,8 @@ The accumulated value can be retrieved via [`calculated_mechanical_work`](@ref). # Arguments - `system`: The [`TotalLagrangianSPHSystem`](@ref) whose particles should be monitored. -- `semi`: The [`Semidiscretization`](@ref) that contains `system`. +- `semi`: The [`Semidiscretization`](@ref TrixiParticles.Semidiscretization) + that contains `system`. # Keywords - `interval=1`: Interval (in number of time steps) at which to compute the instantaneous power. diff --git a/src/callbacks/stepsize.jl b/src/callbacks/stepsize.jl index c34c3d666b..73c0cfe164 100644 --- a/src/callbacks/stepsize.jl +++ b/src/callbacks/stepsize.jl @@ -16,11 +16,11 @@ The step size is therefore only applied once at the beginning of the simulation. The step size ``\Delta t`` is chosen as the minimum ```math - \Delta t = \min(\Delta t_\eta, \Delta t_a, \Delta t_c), + \Delta t = \min(\Delta t_\nu, \Delta t_a, \Delta t_c), ``` where ```math - \Delta t_\eta = 0.125 \, h^2 / \eta, \quad \Delta t_a = 0.25 \sqrt{h / \lVert g \rVert}, + \Delta t_\nu = 0.125 \, h^2 / \nu, \quad \Delta t_a = 0.25 \sqrt{h / \lVert g \rVert}, \quad \Delta t_c = \text{CFL} \, h / c, ``` with ``\nu = \alpha h c / (2n + 4)``, where ``\alpha`` is the parameter of the viscosity diff --git a/src/general/custom_quantities.jl b/src/general/custom_quantities.jl index 3f1f74c637..0b28b48bdb 100644 --- a/src/general/custom_quantities.jl +++ b/src/general/custom_quantities.jl @@ -20,11 +20,10 @@ function kinetic_energy(system::AbstractStructureSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) mass = system.mass - energy = zero(eltype(system)) - return sum(each_active_particle(system)) do particle + return sum(each_active_particle(system); init=zero(eltype(system))) do particle v_i = current_velocity(v, system, particle) - energy += mass[particle] * dot(v_i, v_i) / 2 + return mass[particle] * dot(v_i, v_i) / 2 end end @@ -39,7 +38,7 @@ end Returns the total mass of all particles in a system. """ function total_mass(system, dv_ode, du_ode, v_ode, u_ode, semi, t) - return sum(system.mass) + return sum(active_values(system.mass, system)) end function total_mass(system::AbstractBoundarySystem, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -63,7 +62,7 @@ Returns the maximum pressure over all particles in a system. """ function max_pressure(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - return maximum(current_pressure(v, system)) + return maximum(active_values(current_pressure(v, system), system)) end function max_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -77,7 +76,7 @@ Returns the minimum pressure over all particles in a system. """ function min_pressure(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - return minimum(current_pressure(v, system)) + return minimum(active_values(current_pressure(v, system), system)) end function min_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -91,8 +90,8 @@ Returns the average pressure over all particles in a system. """ function avg_pressure(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - sum_ = sum(current_pressure(v, system)) - return sum_ / nparticles(system) + pressure = active_values(current_pressure(v, system), system) + return sum(pressure) / length(pressure) end function avg_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -106,7 +105,7 @@ Returns the maximum density over all particles in a system. """ function max_density(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - return maximum(current_density(v, system)) + return maximum(active_values(current_density(v, system), system)) end function max_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -120,7 +119,7 @@ Returns the minimum density over all particles in a system. """ function min_density(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - return minimum(current_density(v, system)) + return minimum(active_values(current_density(v, system), system)) end function min_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) @@ -134,10 +133,12 @@ Returns the average_density over all particles in a system. """ function avg_density(system::AbstractFluidSystem, dv_ode, du_ode, v_ode, u_ode, semi, t) v = wrap_v(v_ode, system, semi) - sum_ = sum(current_density(v, system)) - return sum_ / nparticles(system) + density = active_values(current_density(v, system), system) + return sum(density) / length(density) end function avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) return NaN end + +@inline active_values(values, system) = view(values, each_active_particle(system)) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 959ed70deb..10fca415e3 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -189,14 +189,17 @@ end Create an `ODEProblem` from the semidiscretization with the specified `tspan`. # Arguments -- `semi`: A [`Semidiscretization`](@ref) holding the systems involved in the simulation. +- `semi`: A [`Semidiscretization`](@ref TrixiParticles.Semidiscretization) + holding the systems involved in the simulation. - `tspan`: The time span over which the simulation will be run. # Keywords - `restart_with`: Can be used to restart the simulation from VTK solution files (see [`SolutionSavingCallback`](@ref)). This can be either `nothing` (default, no restart) or a `Tuple` of filenames, - one for each system in the [`Semidiscretization`](@ref). - The order of the filenames must match the order of the systems in the [`Semidiscretization`](@ref). + one for each system in the + [`Semidiscretization`](@ref TrixiParticles.Semidiscretization). + The order of the filenames must match the order of the systems in the + [`Semidiscretization`](@ref TrixiParticles.Semidiscretization). Note that `semidiscretize` replaces the initial time (`tspan[1]`) with the timestamp read from the VTK files. If the user-provided `tspan[1]` does not match the restart time, it is adjusted and an info message is logged. If multiple files are provided, their @@ -356,7 +359,8 @@ end Set the initial coordinates and velocities of all systems in `semi` to the final values in the solution `sol`. [`semidiscretize`](@ref) has to be called again afterwards, or another -[`Semidiscretization`](@ref) can be created with the updated systems. +[`Semidiscretization`](@ref TrixiParticles.Semidiscretization) can be created +with the updated systems. # Arguments - `semi`: The semidiscretization diff --git a/src/io/io.jl b/src/io/io.jl index 692dfd8528..188f5f617a 100644 --- a/src/io/io.jl +++ b/src/io/io.jl @@ -319,7 +319,15 @@ function add_system_data!(system_data, surface_normal_method::ColorfieldSurfaceN system_data["surface_normal_method"] = Dict{String, Any}() system_data["surface_normal_method"]["model"] = type2string(surface_normal_method) system_data["surface_normal_method"]["boundary_contact_threshold"] = surface_normal_method.boundary_contact_threshold + system_data["surface_normal_method"]["interface_threshold"] = surface_normal_method.interface_threshold system_data["surface_normal_method"]["ideal_density_threshold"] = surface_normal_method.ideal_density_threshold + system_data["surface_normal_method"]["interface_taper_start"] = surface_normal_method.interface_taper_start + system_data["surface_normal_method"]["support_taper_width"] = surface_normal_method.support_taper_width +end + +function add_system_data!(system_data, surface_normal_method::CorrectedCSFSurfaceNormal) + system_data["surface_normal_method"] = Dict{String, Any}() + system_data["surface_normal_method"]["model"] = type2string(surface_normal_method) end function add_system_data!(system_data, boundary_zone::BoundaryZone, indice) diff --git a/src/preprocessing/geometries/geometries.jl b/src/preprocessing/geometries/geometries.jl index c1c0948522..5336c10c67 100644 --- a/src/preprocessing/geometries/geometries.jl +++ b/src/preprocessing/geometries/geometries.jl @@ -4,6 +4,70 @@ include("io.jl") @inline eachface(mesh) = Base.OneTo(nfaces(mesh)) +""" + is_closed_geometry(geometry) + +Return `true` if a polygon or triangle mesh forms a closed region or surface. +""" +function is_closed_geometry(polygon::Polygon) + vertex_degrees = polygon_vertex_degrees(polygon) + + return !isempty(vertex_degrees) && all(==(2), values(vertex_degrees)) +end + +function polygon_vertex_degrees(polygon) + VERTEX = typeof(first(first(polygon.edge_vertices))) + vertex_degrees = Dict{VERTEX, Int}() + + for edge in polygon.edge_vertices + for vertex in edge + vertex_degrees[vertex] = get(vertex_degrees, vertex, 0) + 1 + end + end + + return vertex_degrees +end + +function is_closed_geometry(mesh::TriangleMesh) + return all(==(2), edge_face_counts(mesh)) +end + +function require_closed_geometry(geometry, operation) + is_closed_geometry(geometry) && return nothing + + msg = "`$operation` requires a closed geometry. " * + closure_error_detail(geometry) + + throw(ArgumentError(msg)) +end + +function closure_error_detail(polygon::Polygon) + invalid_vertices = count(!=(2), values(polygon_vertex_degrees(polygon))) + + return "Found $invalid_vertices polygon vertices with an incident-edge count " * + "different from 2. If the vertices already trace a complete 2D boundary, " * + "construct or load the geometry with `close_curve=true`; otherwise provide " * + "a closed boundary." +end + +function closure_error_detail(mesh::TriangleMesh) + invalid_edges = count(!=(2), edge_face_counts(mesh)) + + return "Found $invalid_edges mesh edges with an incident-face count different from 2." +end + +function edge_face_counts(mesh::TriangleMesh) + edge_face_counts = zeros(Int, length(mesh.edge_vertices_ids)) + + for face_edges in mesh.face_edges_ids + edge_face_counts[face_edges[1]] += 1 + edge_face_counts[face_edges[2]] += 1 + edge_face_counts[face_edges[3]] += 1 + end + + return edge_face_counts +end + function Base.setdiff(initial_condition::InitialCondition, geometries::Union{Polygon, TriangleMesh}...) geometry = first(geometries) @@ -11,6 +75,7 @@ function Base.setdiff(initial_condition::InitialCondition, if ndims(geometry) != ndims(initial_condition) throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition")) end + require_closed_geometry(geometry, "setdiff") coords = reinterpret(reshape, SVector{ndims(geometry), eltype(initial_condition.coordinates)}, @@ -41,6 +106,7 @@ function Base.intersect(initial_condition::InitialCondition, if ndims(geometry) != ndims(initial_condition) throw(ArgumentError("all passed geometries must have the same dimensionality as the initial condition")) end + require_closed_geometry(geometry, "intersect") coords = reinterpret(reshape, SVector{ndims(geometry), eltype(initial_condition.coordinates)}, diff --git a/src/preprocessing/geometries/io.jl b/src/preprocessing/geometries/io.jl index bbccc276ec..a7c6106827 100644 --- a/src/preprocessing/geometries/io.jl +++ b/src/preprocessing/geometries/io.jl @@ -1,5 +1,5 @@ """ - load_geometry(filename; element_type=Float64) + load_geometry(filename; element_type=Float64, close_curve=true) Load file and return corresponding type for [`ComplexShape`](@ref). Supported file formats are `.stl`, `.asc` and `dxf`. @@ -18,16 +18,20 @@ For comprehensive information about the supported file formats, refer to the doc # Keywords - `element_type`: Element type (default is `Float64`) +- `close_curve`: Close 2D `.asc` and `.dxf` curves by appending the first point + when it is not already repeated. This assumes the vertices already + trace a complete, ordered boundary. Set this to `false` for intentional + open curves. Region sampling and classification reject open geometries. """ -function load_geometry(filename; element_type=Float64) +function load_geometry(filename; element_type=Float64, close_curve=true) ELTYPE = element_type file_extension = splitext(filename)[end] if file_extension == ".asc" - geometry = load_ascii(filename; ELTYPE, skipstart=1) + geometry = load_ascii(filename; ELTYPE, skipstart=1, close_curve) elseif file_extension == ".dxf" - geometry = load_dxf(filename; ELTYPE) + geometry = load_dxf(filename; ELTYPE, close_curve) elseif file_extension == ".stl" geometry = load(FileIO.query(filename); ELTYPE) else @@ -37,21 +41,21 @@ function load_geometry(filename; element_type=Float64) return geometry end -function load_ascii(filename; ELTYPE=Float64, skipstart=1) +function load_ascii(filename; ELTYPE=Float64, skipstart=1, close_curve=true) # Read the data from the ASCII file in as a matrix of coordinates. # Ignore the first `skipstart` lines of the file (e.g. headers). points = DelimitedFiles.readdlm(filename, ' ', ELTYPE, '\n'; skipstart)[:, 1:2] - return Polygon(copy(points')) + return Polygon(copy(points'); close_curve) end -function load_dxf(filename; ELTYPE=Float64) +function load_dxf(filename; ELTYPE=Float64, close_curve=true) points = Tuple{ELTYPE, ELTYPE}[] load_dxf!(points, filename) - return Polygon(stack(points)) + return Polygon(stack(points); close_curve) end function load_dxf!(points::Vector{Tuple{T, T}}, filename) where {T} diff --git a/src/preprocessing/geometries/polygon.jl b/src/preprocessing/geometries/polygon.jl index c56315b729..c2614f4571 100644 --- a/src/preprocessing/geometries/polygon.jl +++ b/src/preprocessing/geometries/polygon.jl @@ -8,21 +8,30 @@ struct Polygon{NDIMS, ELTYPE} min_corner :: SVector{NDIMS, ELTYPE} max_corner :: SVector{NDIMS, ELTYPE} - function Polygon(vertices) + function Polygon(vertices; close_curve=true) NDIMS = size(vertices, 1) - return Polygon{NDIMS}(vertices) + return Polygon{NDIMS}(vertices; close_curve) end # Function barrier to make `NDIMS` static and therefore `SVector`s type-stable - function Polygon{NDIMS}(vertices_) where {NDIMS} - n_vertices = size(vertices_, 2) + function Polygon{NDIMS}(vertices_; close_curve=true) where {NDIMS} ELTYPE = eltype(vertices_) - min_corner = SVector{NDIMS}(minimum(vertices_, dims=2)) - max_corner = SVector{NDIMS}(maximum(vertices_, dims=2)) + vertices = collect(reinterpret(reshape, SVector{NDIMS, ELTYPE}, vertices_)) - vertices = reinterpret(reshape, SVector{NDIMS, ELTYPE}, vertices_) + if length(vertices) < 3 + throw(ArgumentError("polygon requires at least three vertices")) + end + + if close_curve && !isapprox(first(vertices), last(vertices)) + push!(vertices, first(vertices)) + end + + n_vertices = length(vertices) + + min_corner = SVector([minimum(v[i] for v in vertices) for i in 1:NDIMS]...) + max_corner = SVector([maximum(v[i] for v in vertices) for i in 1:NDIMS]...) # Sum over all the edges and determine if the vertices are in clockwise order # to make sure that all normals pointing outwards. @@ -63,6 +72,10 @@ struct Polygon{NDIMS, ELTYPE} push!(edge_normals, edge_normal) end + if length(edge_vertices) < 3 + throw(ArgumentError("polygon requires at least three non-degenerate edges")) + end + vertex_normals = Vector{NTuple{2, SVector{NDIMS, ELTYPE}}}() # Calculate vertex pseudo-normals. @@ -95,6 +108,63 @@ struct Polygon{NDIMS, ELTYPE} return new{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals, edge_vertices_ids, min_corner, max_corner) end + + function Polygon{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, + edge_normals, edge_vertices_ids, + min_corner, max_corner) where {NDIMS, ELTYPE} + return new{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals, + edge_vertices_ids, min_corner, max_corner) + end +end + +function vertex_normals_from_edges(edge_vertices, edge_normals) + VERTEX = typeof(first(first(edge_vertices))) + normal_sums = Dict{VERTEX, VERTEX}() + + for (edge, edge_normal) in zip(edge_vertices, edge_normals) + for vertex in edge + normal_sums[vertex] = get(normal_sums, vertex, zero(edge_normal)) + edge_normal + end + end + + return map(edge_vertices, edge_normals) do edge, edge_normal + normals = map(edge) do vertex + normal_sum = normal_sums[vertex] + normal_norm = norm(normal_sum) + + return iszero(normal_norm) ? edge_normal : normal_sum / normal_norm + end + + return Tuple(normals) + end +end + +function rebuild_polygon_from_edges(edge_vertices, edge_normals) + NDIMS = length(first(edge_normals)) + ELTYPE = eltype(first(edge_normals)) + vertices = SVector{NDIMS, ELTYPE}[] + vertex_ids = Dict{SVector{NDIMS, ELTYPE}, Int}() + + edge_vertices_ids = map(edge_vertices) do edge + v1, v2 = edge + id1 = get!(vertex_ids, v1) do + push!(vertices, v1) + return length(vertices) + end + id2 = get!(vertex_ids, v2) do + push!(vertices, v2) + return length(vertices) + end + + return (id1, id2) + end + + min_corner = SVector([minimum(v[i] for v in vertices) for i in 1:NDIMS]...) + max_corner = SVector([maximum(v[i] for v in vertices) for i in 1:NDIMS]...) + vertex_normals = vertex_normals_from_edges(edge_vertices, edge_normals) + + return Polygon{NDIMS, ELTYPE}(vertices, edge_vertices, vertex_normals, edge_normals, + edge_vertices_ids, min_corner, max_corner) end function Base.show(io::IO, geometry::Polygon) @@ -119,14 +189,23 @@ end @inline Base.eltype(::Polygon{NDIMS, ELTYPE}) where {NDIMS, ELTYPE} = ELTYPE -@inline function Base.deleteat!(polygon::Polygon, indices) - (; edge_vertices, edge_normals, edge_vertices_ids) = polygon +""" + delete_faces(geometry, indices) + +Return a geometry with the faces at `indices` removed and derived geometry data rebuilt. +""" +@inline function delete_faces(polygon::Polygon, indices) + edge_vertices = copy(polygon.edge_vertices) + edge_normals = copy(polygon.edge_normals) deleteat!(edge_vertices, indices) - deleteat!(edge_vertices_ids, indices) deleteat!(edge_normals, indices) - return polygon + if isempty(edge_vertices) + throw(ArgumentError("cannot delete all polygon edges")) + end + + return rebuild_polygon_from_edges(edge_vertices, edge_normals) end @inline nfaces(mesh::Polygon) = length(mesh.edge_normals) diff --git a/src/preprocessing/geometries/triangle_mesh.jl b/src/preprocessing/geometries/triangle_mesh.jl index ede02d4682..7677784854 100644 --- a/src/preprocessing/geometries/triangle_mesh.jl +++ b/src/preprocessing/geometries/triangle_mesh.jl @@ -129,19 +129,16 @@ struct TriangleMesh{NDIMS, ELTYPE} min_corner = SVector([minimum(v[i] for v in vertices) for i in 1:NDIMS]...) max_corner = SVector([maximum(v[i] for v in vertices) for i in 1:NDIMS]...) - for i in eachindex(edge_normals) - # Skip zero normals, which would be normalized to `NaN` vectors. - # The edge normals are only used for the `SignedDistanceField`, which is - # essential for the packing. - # Zero normals are caused by exactly or nearly duplicated faces. - if !iszero(norm(edge_normals[i])) - edge_normals[i] = normalize(edge_normals[i]) + for normals in (edge_normals, vertex_normals) + for i in eachindex(normals) + normals_norm = norm(normals[i]) + !iszero(normals_norm) && (normals[i] = normals[i] / normals_norm) end end return new{NDIMS, ELTYPE}(vertices, face_vertices, face_vertices_ids, face_edges_ids, edge_vertices_ids, - normalize.(vertex_normals), edge_normals, + vertex_normals, edge_normals, face_normals, min_corner, max_corner) end end @@ -171,15 +168,19 @@ end @inline face_normal(triangle, geometry::TriangleMesh) = geometry.face_normals[triangle] -@inline function Base.deleteat!(mesh::TriangleMesh, indices) - (; face_vertices, face_vertices_ids, face_edges_ids, face_normals) = mesh +@inline function delete_faces(mesh::TriangleMesh, indices) + face_vertices = copy(mesh.face_vertices) + face_normals = copy(mesh.face_normals) deleteat!(face_vertices, indices) - deleteat!(face_vertices_ids, indices) - deleteat!(face_edges_ids, indices) deleteat!(face_normals, indices) - return mesh + if isempty(face_vertices) + throw(ArgumentError("cannot delete all triangle mesh faces")) + end + + vertices = collect(Iterators.flatten(face_vertices)) + return TriangleMesh(face_vertices, face_normals, vertices) end @inline nfaces(mesh::TriangleMesh) = length(mesh.face_normals) diff --git a/src/preprocessing/particle_packing/signed_distance.jl b/src/preprocessing/particle_packing/signed_distance.jl index 01e862365f..2ef2c5eae4 100644 --- a/src/preprocessing/particle_packing/signed_distance.jl +++ b/src/preprocessing/particle_packing/signed_distance.jl @@ -16,11 +16,15 @@ to this surface. distance of `abs(max_signed_distance)` to the surface of the shape will be sampled. - `points`: Points on which the signed distance is computed. + Pass a collection of static vectors or an `NDIMS`-by-`N` matrix with + one point per column. When set to `nothing` (default), the bounding box of the shape will be sampled with a uniform grid of points. - `use_for_boundary_packing`: Set to `true` if [`SignedDistanceField`] is used to pack a boundary [`ParticlePackingSystem`](@ref). Use the default of `false` when packing without a boundary. + This requires a closed geometry, since boundary packing + needs a well-defined outside region. """ struct SignedDistanceField{ELTYPE, P, N, D} positions :: P @@ -38,6 +42,11 @@ function SignedDistanceField(geometry, particle_spacing; NDIMS = ndims(geometry) ELTYPE = eltype(particle_spacing) + if use_for_boundary_packing + require_closed_geometry(geometry, + "SignedDistanceField with `use_for_boundary_packing=true`") + end + sdf_factor = use_for_boundary_packing ? 2 : 1 search_radius = sdf_factor * max_signed_distance @@ -62,9 +71,11 @@ function SignedDistanceField(geometry, particle_spacing; min_corner; place_on_shell=true) points = reinterpret(reshape, SVector{NDIMS, eltype(grid)}, grid) + else + points = wrap_points(points, Val(NDIMS)) end - positions = copy(points) + positions = collect(points) # This gives a performance boost for large geometries delete_positions_in_empty_cells!(positions, nhs) diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl index a4bb7e9aa7..45585541ee 100644 --- a/src/preprocessing/particle_packing/system.jl +++ b/src/preprocessing/particle_packing/system.jl @@ -5,6 +5,7 @@ smoothing_length=shape.particle_spacing, smoothing_length_interpolation=smoothing_length, is_boundary=false, boundary_compress_factor=1, + boundary_thickness=nothing, neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(), background_pressure, place_on_shell=false, fixed_system=false) @@ -26,10 +27,6 @@ For more information on the methods, see [particle packing](@ref particle_packin - `is_boundary`: When `shape` is inside the geometry that was used to create `signed_distance_field`, set `is_boundary=false`. Otherwise (`shape` is the sampled boundary), set `is_boundary=true`. - The thickness of the boundary is specified by creating - `signed_distance_field` with: - - `use_for_boundary_packing=true` - - `max_signed_distance=boundary_thickness` See [`SignedDistanceField`](@ref). - `fixed_system`: When set to `true`, the system remains static, meaning particles will not move and the `InitialCondition` will stay unchanged. @@ -54,6 +51,10 @@ For more information on the methods, see [particle packing](@ref particle_packin Compression can be useful for highly convex geometries, where the boundary volume increases significantly while the mass of the boundary particles remains constant. Recommended values are `0.8` or `0.9`. +- `boundary_thickness`: Thickness of the sampled boundary when `is_boundary=true`. + By default, this is `signed_distance_field.max_signed_distance`. + If [`sample_boundary`](@ref) used a smaller `boundary_thickness` + than the `SignedDistanceField`, pass the same value here. """ struct ParticlePackingSystem{S, F, NDIMS, ELTYPE <: Real, PR, C, AV, IC, M, D, K, N, SD} <: AbstractFluidSystem{NDIMS} @@ -106,6 +107,7 @@ function ParticlePackingSystem(shape::InitialCondition; smoothing_length=shape.particle_spacing, smoothing_length_interpolation=smoothing_length, is_boundary=false, boundary_compress_factor=1, + boundary_thickness=nothing, neighborhood_search=GridNeighborhoodSearch{ndims(shape)}(), background_pressure, place_on_shell=false, fixed_system=false) @@ -147,10 +149,30 @@ function ParticlePackingSystem(shape::InitialCondition; # Its value is negative if the particle is inside the geometry. # Otherwise (if outside), the value is positive. if is_boundary - offset = place_on_shell ? shape.particle_spacing : shape.particle_spacing / 2 + if isnothing(signed_distance_field) + fixed_system || + throw(ArgumentError("`signed_distance_field` is required when `is_boundary=true`")) + + shift_length = zero(ELTYPE) + else + boundary_thickness_ = isnothing(boundary_thickness) ? + signed_distance_field.max_signed_distance : + convert(ELTYPE, boundary_thickness) + + if boundary_thickness_ > signed_distance_field.max_signed_distance + throw(ArgumentError("`boundary_thickness` is greater than " * + "`max_signed_distance` of `SignedDistanceField`.")) + end + + if boundary_thickness_ < zero(boundary_thickness_) + throw(ArgumentError("`boundary_thickness` must be non-negative")) + end - shift_length = -boundary_compress_factor * - signed_distance_field.max_signed_distance - offset + offset = place_on_shell ? shape.particle_spacing : shape.particle_spacing / 2 + + shift_length = -boundary_compress_factor * + boundary_thickness_ - offset + end else shift_length = place_on_shell ? zero(ELTYPE) : shape.particle_spacing / 2 end diff --git a/src/preprocessing/point_in_poly/winding_number_hormann.jl b/src/preprocessing/point_in_poly/winding_number_hormann.jl index a93cf789c4..9f9b849acb 100644 --- a/src/preprocessing/point_in_poly/winding_number_hormann.jl +++ b/src/preprocessing/point_in_poly/winding_number_hormann.jl @@ -14,6 +14,7 @@ struct WindingNumberHormann end # https://doi.org/10.1016/S0925-7721(01)00012-8 function (point_in_poly::WindingNumberHormann)(geometry, points; store_winding_number=false) (; edge_vertices) = geometry + points = wrap_points(points, Val(ndims(geometry))) # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe inpoly = fill(false, length(points)) diff --git a/src/preprocessing/point_in_poly/winding_number_jacobson.jl b/src/preprocessing/point_in_poly/winding_number_jacobson.jl index 145a6d9a28..14a0d62cae 100644 --- a/src/preprocessing/point_in_poly/winding_number_jacobson.jl +++ b/src/preprocessing/point_in_poly/winding_number_jacobson.jl @@ -51,14 +51,16 @@ end """ WindingNumberJacobson(; geometry=nothing, winding_number_factor=sqrt(eps()), - hierarchical_winding=false) + hierarchical_winding=!isnothing(geometry)) Algorithm for inside-outside segmentation of a complex geometry proposed by [Jacobson2013](@cite). # Keywords - `geometry`: Complex geometry returned by [`load_geometry`](@ref) and is only required when using `hierarchical_winding=true`. - `hierarchical_winding`: If set to `true`, an optimized hierarchical approach will be used, - which gives a significant speedup. For further information see [Hierarchical Winding](@ref hierarchical_winding). + which gives a significant speedup. It defaults to `true` when `geometry` + is passed and `false` otherwise. For further information see + [Hierarchical Winding](@ref hierarchical_winding). - `winding_number_factor`: For leaky geometries, a factor of `0.4` will give a better inside-outside segmentation. !!! warning "Experimental Implementation" @@ -69,7 +71,7 @@ struct WindingNumberJacobson{ELTYPE, W} winding :: W function WindingNumberJacobson(; geometry=nothing, winding_number_factor=sqrt(eps()), - hierarchical_winding=true) + hierarchical_winding=(!isnothing(geometry))) if hierarchical_winding && geometry isa Nothing throw(ArgumentError("`geometry` must be of type `Polygon` (2D) or `TriangleMesh` (3D) when using hierarchical winding")) end @@ -104,6 +106,7 @@ end function (point_in_poly::WindingNumberJacobson)(geometry, points; store_winding_number=false) (; winding_number_factor, winding) = point_in_poly + points = wrap_points(points, Val(ndims(geometry))) # We cannot use a `BitVector` here, as writing to a `BitVector` is not thread-safe inpoly = fill(false, length(points)) diff --git a/src/preprocessing/preprocessing.jl b/src/preprocessing/preprocessing.jl index c538280bb4..1f29263cae 100644 --- a/src/preprocessing/preprocessing.jl +++ b/src/preprocessing/preprocessing.jl @@ -1,3 +1,19 @@ +function wrap_points(points, ::Val{NDIMS}) where {NDIMS} + if points isa AbstractMatrix + if size(points, 1) != NDIMS + throw(ArgumentError("point matrix must have $NDIMS rows")) + end + + # Interpret an `NDIMS`-by-`N` matrix as one static vector per column. Constructing + # the vectors explicitly also supports non-contiguous matrix views. + return map(eachcol(points)) do point + return SVector{NDIMS, eltype(points)}(point) + end + end + + return points +end + include("geometries/geometries.jl") include("point_in_poly/point_in_poly.jl") include("particle_packing/particle_packing.jl") diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl index c7f29b1e36..175eee289b 100644 --- a/src/schemes/boundary/open_boundary/boundary_zones.jl +++ b/src/schemes/boundary/open_boundary/boundary_zones.jl @@ -384,6 +384,16 @@ function set_up_boundary_zone(boundary_face, face_normal, density, particle_spac flow_direction = zero(face_normal) end + # Validate boundary geometry before sampling particles. + unit_spanning_set, + _ = calculate_spanning_vectors(boundary_face, + one(eltype(face_normal))) + dot_face_normal = dot(normalize(unit_spanning_set[:, 1]), face_normal) + + if !isapprox(abs(dot_face_normal), 1) + throw(ArgumentError("`face_normal` is not normal to the boundary face")) + end + # Sample particles in boundary zone if isnothing(initial_condition) && isnothing(extrude_geometry) initial_condition = TrixiParticles.extrude_geometry(boundary_face; particle_spacing, @@ -408,13 +418,6 @@ function set_up_boundary_zone(boundary_face, face_normal, density, particle_spac # Vectors spanning the boundary zone/box spanning_set, zone_origin = calculate_spanning_vectors(boundary_face, zone_width) - # First vector of `spanning_vectors` is normal to the boundary face. - dot_face_normal = dot(normalize(spanning_set[:, 1]), face_normal) - - if !isapprox(abs(dot_face_normal), 1) - throw(ArgumentError("`face_normal` is not normal to the boundary face")) - end - if boundary_type isa InFlow # First vector of `spanning_vectors` is normal to the boundary face dot_flow = dot(normalize(spanning_set[:, 1]), flow_direction) @@ -463,11 +466,22 @@ function spanning_vectors(face_vertices::NTuple{3}, zone_width) edge1 = face_vertices[2] - face_vertices[1] edge2 = face_vertices[3] - face_vertices[1] + edge1_norm = norm(edge1) + edge2_norm = norm(edge2) + edge_tolerance = sqrt(eps(typeof(edge1_norm * edge2_norm))) * edge1_norm * + edge2_norm + # Check if the edges are linearly dependent (to avoid degenerate planes) - if isapprox(norm(cross(edge1, edge2)), 0.0; atol=eps()) + cross_norm = norm(cross(edge1, edge2)) + if isapprox(cross_norm, zero(cross_norm); atol=edge_tolerance) throw(ArgumentError("the vectors `AB` and `AC` must not be collinear")) end + edge_dot = dot(edge1, edge2) + if !isapprox(edge_dot, zero(edge_dot); atol=edge_tolerance) + throw(ArgumentError("the vectors `AB` and `AC` must be orthogonal")) + end + # Calculate normal vector of `boundary_face` c = Vector(normalize(cross(edge2, edge1)) * zone_width) diff --git a/src/schemes/boundary/open_boundary/method_of_characteristics.jl b/src/schemes/boundary/open_boundary/method_of_characteristics.jl index ce04478a0c..5023b9023f 100644 --- a/src/schemes/boundary/open_boundary/method_of_characteristics.jl +++ b/src/schemes/boundary/open_boundary/method_of_characteristics.jl @@ -169,6 +169,7 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) # Particle is outside of the influence of fluid particles. # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1. if volume[particle] * smoothing_length^ndims(system) < eps(eltype(smoothing_length)) + zone_id = system.boundary_zone_indices[particle] # Using the average of the values at the previous time step for particles which # are outside of the influence of fluid particles. @@ -178,6 +179,8 @@ function evaluate_characteristics!(system, v, u, v_ode, u_ode, semi, t) counter = 0 for neighbor in each_integrated_particle(system) + system.boundary_zone_indices[neighbor] == zone_id || continue + # Make sure that only neighbors in the influence of # the fluid particles are used. # `volume` is in the order of 1 / h^d, so volume * h^d is in the order of 1. diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl index 9da559481e..92c9b4d2d6 100644 --- a/src/schemes/boundary/open_boundary/system.jl +++ b/src/schemes/boundary/open_boundary/system.jl @@ -1,7 +1,9 @@ @doc raw""" OpenBoundarySystem(boundary_zone::BoundaryZone; - fluid_system::AbstractFluidSystem, buffer_size::Integer, - boundary_model, calculate_flow_rate=false) + fluid_system::AbstractFluidSystem, + buffer_size=default_open_boundary_buffer_size(fluid_system), + boundary_model=BoundaryModelMirroringTafuni(), + calculate_flow_rate=false) Open boundary system for in- and outflow particles. @@ -10,7 +12,10 @@ Open boundary system for in- and outflow particles. # Keywords - `fluid_system`: The corresponding fluid system -- `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models)) +- `buffer_size`: Number of buffer particles for the boundary system. + Defaults to the buffer size of `fluid_system`. +- `boundary_model`: Boundary model (see [Open Boundary Models](@ref open_boundary_models)). + Defaults to [`BoundaryModelMirroringTafuni`](@ref). - `calculate_flow_rate=false`: Set to `true` to calculate the volumetric flow rate through each boundary zone. This value is automatically enabled when using [`RCRWindkesselModel`](@ref). Otherwise, it is useful only for postprocessing. @@ -49,6 +54,18 @@ struct OpenBoundarySystem{BM, ELTYPE, NDIMS, IC, FS, FSI, K, ARRAY1D, BC, FC, BZ cache :: C end +function default_open_boundary_buffer_size(fluid_system) + fluid_buffer = buffer(fluid_system) + + if fluid_buffer isa SystemBuffer + return fluid_buffer.buffer_size + end + + throw(ArgumentError("`buffer_size` could not be inferred for `OpenBoundarySystem` " * + "because `fluid_system` has no buffer. Pass `buffer_size=...` " * + "explicitly or construct `fluid_system` with `buffer_size=...`.")) +end + function OpenBoundarySystem(boundary_model, initial_condition, fluid_system, fluid_system_index, smoothing_kernel, smoothing_length, mass, volume, boundary_candidates, fluid_candidates, @@ -70,8 +87,10 @@ function OpenBoundarySystem(boundary_model, initial_condition, fluid_system, end function OpenBoundarySystem(boundary_zones::Union{BoundaryZone, Nothing}...; - fluid_system::AbstractFluidSystem, buffer_size::Integer, - boundary_model, calculate_flow_rate=false, + fluid_system::AbstractFluidSystem, + buffer_size=default_open_boundary_buffer_size(fluid_system), + boundary_model=BoundaryModelMirroringTafuni(), + calculate_flow_rate=false, pressure_acceleration=fluid_system.pressure_acceleration_formulation, shifting_technique=boundary_model isa BoundaryModelDynamicalPressureZhang ? @@ -807,9 +826,9 @@ function check_configuration(system::OpenBoundarySystem, systems, neighborhood_s system.fluid_system_index[] = fluid_system_index if boundary_model isa BoundaryModelCharacteristicsLastiwka && - any(zone -> isnothing(zone.flow_direction), boundary_zones) - throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a specific flow direction. " * - "Please specify `InFlow()` and `OutFlow()`.")) + any(zone -> zone.is_bidirectional, boundary_zones) + throw(ArgumentError("`BoundaryModelCharacteristicsLastiwka` needs a directed boundary zone. " * + "Please specify `InFlow()` or `OutFlow()` instead of `BidirectionalFlow()`.")) end if first(PointNeighbors.requires_update(neighborhood_search)) diff --git a/src/schemes/boundary/wall_boundary/dummy_particles.jl b/src/schemes/boundary/wall_boundary/dummy_particles.jl index 55aa20226b..e0ec835c1e 100644 --- a/src/schemes/boundary/wall_boundary/dummy_particles.jl +++ b/src/schemes/boundary/wall_boundary/dummy_particles.jl @@ -76,6 +76,46 @@ struct BoundaryModelDummyParticles{DC, SE, CLIP, ELTYPE <: Real, VECTOR, K, V, C end end +@doc raw""" + BoundaryModelDummyParticles(initial_condition; + fluid_system::AbstractFluidSystem, + initial_density=initial_condition.density, + hydrodynamic_mass=initial_condition.mass, + boundary_density_calculator=AdamiPressureExtrapolation(), + smoothing_kernel=system_smoothing_kernel(fluid_system), + smoothing_length=initial_smoothing_length(fluid_system), + viscosity=nothing, + state_equation=system_state_equation(fluid_system), + correction=system_correction(fluid_system), + clip_negative_pressure=false, + reference_particle_spacing=default_reference_particle_spacing(fluid_system)) + +High-level convenience constructor for dummy-particle wall models that infers the kernel, +smoothing length, correction, and equation-of-state-related settings from the adjacent +`fluid_system`. +""" +function BoundaryModelDummyParticles(initial_condition; + fluid_system::AbstractFluidSystem, + initial_density=initial_condition.density, + hydrodynamic_mass=initial_condition.mass, + boundary_density_calculator=AdamiPressureExtrapolation(), + smoothing_kernel=system_smoothing_kernel(fluid_system), + smoothing_length=initial_smoothing_length(fluid_system), + viscosity=nothing, + state_equation=system_state_equation(fluid_system), + correction=system_correction(fluid_system), + clip_negative_pressure=false, + reference_particle_spacing=default_reference_particle_spacing(fluid_system)) + return BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, + boundary_density_calculator, smoothing_kernel, + smoothing_length; + viscosity, state_equation, correction, + clip_negative_pressure, + reference_particle_spacing) +end + +# The default constructor needs to be accessible for Adapt.jl to work with this struct. +# See the comments in general/gpu.jl for more details. function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, density_calculator, smoothing_kernel, smoothing_length; viscosity=nothing, @@ -109,6 +149,15 @@ function BoundaryModelDummyParticles(initial_density, hydrodynamic_mass, clip_negative_pressure) end +@inline function default_reference_particle_spacing(fluid_system) + if hasproperty(fluid_system, :cache) && + hasproperty(fluid_system.cache, :reference_particle_spacing) + return fluid_system.cache.reference_particle_spacing + end + + return zero(eltype(fluid_system)) +end + @inline function Base.ndims(boundary_model::BoundaryModelDummyParticles) return ndims(boundary_model.smoothing_kernel) end diff --git a/src/schemes/fluid/entropically_damped_sph/rhs.jl b/src/schemes/fluid/entropically_damped_sph/rhs.jl index 2ea1e4dc4d..380d3fed68 100644 --- a/src/schemes/fluid/entropically_damped_sph/rhs.jl +++ b/src/schemes/fluid/entropically_damped_sph/rhs.jl @@ -20,6 +20,19 @@ function interact!(dv, v_particle_system, u_particle_system, h = initial_smoothing_length(particle_system) almostzero = sqrt(eps(h^2)) + if particle_system === neighbor_system + @threaded semi for particle in each_integrated_particle(particle_system) + rho_a = @inbounds current_density(v_particle_system, particle_system, + particle) + v_a = @inbounds current_velocity(v_particle_system, particle_system, particle) + acceleration = surface_tension_acceleration(surface_tension_a, particle_system, + particle, rho_a, v_a) + for i in 1:ndims(particle_system) + @inbounds dv[i, particle] += acceleration[i] + end + end + end + # Loop over all pairs of particles and neighbors within the kernel cutoff foreach_point_neighbor(particle_system, neighbor_system, system_coords, neighbor_coords, semi; diff --git a/src/schemes/fluid/entropically_damped_sph/system.jl b/src/schemes/fluid/entropically_damped_sph/system.jl index 2b085d3ebb..aee6e1ee16 100644 --- a/src/schemes/fluid/entropically_damped_sph/system.jl +++ b/src/schemes/fluid/entropically_damped_sph/system.jl @@ -122,9 +122,10 @@ function EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, smooth if surface_tension !== nothing && surface_normal_method === nothing surface_normal_method = ColorfieldSurfaceNormal() end + validate_corrected_csf(surface_normal_method, surface_tension) if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method or a surface tension model")) end if correction isa ShepardKernelCorrection && @@ -299,11 +300,54 @@ function update_quantities!(system::EntropicallyDampedSPHSystem, v, u, end function update_pressure!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t) + (; correction, density_calculator) = system + + # These are only computed when using corrections + compute_correction_values!(system, correction, u, v_ode, u_ode, semi) + compute_gradient_correction_matrix!(correction, system, u, v_ode, u_ode, semi) + # `kernel_correct_density!` only performed for `SummationDensity` + kernel_correct_density!(system, v, u, v_ode, u_ode, semi, correction, + density_calculator) + compute_surface_normal!(system, system.surface_normal_method, v, u, v_ode, u_ode, semi, t) compute_surface_delta_function!(system, system.surface_tension, semi) end +function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, + semi, correction, density_calculator) + return system +end + +function kernel_correct_density!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, + semi, corr::ShepardKernelCorrection, ::SummationDensity) + system.cache.density ./= system.cache.kernel_correction_coefficient +end + +function compute_gradient_correction_matrix!(correction, + system::EntropicallyDampedSPHSystem, u, + v_ode, u_ode, semi) + return system +end + +function compute_gradient_correction_matrix!(corr::Union{GradientCorrection, + BlendedGradientCorrection, + MixedKernelGradientCorrection}, + system::EntropicallyDampedSPHSystem, u, + v_ode, u_ode, semi) + (; cache, correction, smoothing_kernel) = system + (; correction_matrix) = cache + + system_coords = current_coordinates(u, system) + + compute_gradient_correction_matrix!(correction_matrix, system, system_coords, + v_ode, u_ode, semi, correction, smoothing_kernel) +end + +@inline function correction_matrix(system::EntropicallyDampedSPHSystem, particle) + extract_smatrix(system.cache.correction_matrix, system, particle) +end + function update_final!(system::EntropicallyDampedSPHSystem, v, u, v_ode, u_ode, semi, t; kwargs...) (; surface_tension) = system @@ -381,6 +425,6 @@ function restart_with!(system::EntropicallyDampedSPHSystem, v, u) for particle in each_integrated_particle(system) system.initial_condition.coordinates[:, particle] .= u[:, particle] system.initial_condition.velocity[:, particle] .= v[1:ndims(system), particle] - system.initial_condition.pressure[particle] = v[end, particle] + system.initial_condition.pressure[particle] = v[ndims(system) + 1, particle] end end diff --git a/src/schemes/fluid/fluid.jl b/src/schemes/fluid/fluid.jl index fb49a65359..777480c923 100644 --- a/src/schemes/fluid/fluid.jl +++ b/src/schemes/fluid/fluid.jl @@ -297,11 +297,13 @@ function restart_v(system::AbstractFluidSystem, data) end function check_configuration(fluid_system::AbstractFluidSystem, systems, nhs) - if !(fluid_system isa ParticlePackingSystem) && !isnothing(fluid_system.surface_tension) + if !(fluid_system isa ParticlePackingSystem) && + (!isnothing(fluid_system.surface_tension) || + !isnothing(fluid_system.surface_normal_method)) foreach_system(systems) do neighbor - if neighbor isa AbstractFluidSystem && - isnothing(fluid_system.surface_tension) && - isnothing(fluid_system.surface_normal_method) + if neighbor isa AbstractFluidSystem && !(neighbor isa ParticlePackingSystem) && + isnothing(neighbor.surface_tension) && + isnothing(neighbor.surface_normal_method) throw(ArgumentError("either none or all fluid systems in a simulation need " * "to use a surface tension model or a surface normal method.")) end diff --git a/src/schemes/fluid/implicit_incompressible_sph/system.jl b/src/schemes/fluid/implicit_incompressible_sph/system.jl index 9b8fedc288..c7df1ce185 100644 --- a/src/schemes/fluid/implicit_incompressible_sph/system.jl +++ b/src/schemes/fluid/implicit_incompressible_sph/system.jl @@ -199,6 +199,8 @@ end return system.density end +@inline system_state_equation(system::ImplicitIncompressibleSPHSystem) = nothing + # TODO: What do we do with the sound speed? This is needed for the viscosity. @inline system_sound_speed(system::ImplicitIncompressibleSPHSystem) = system.artificial_sound_speed @@ -487,7 +489,8 @@ function calculate_sum_d_ij_pj!(sum_d_ij_pj, system, (; time_step) = system system_coords = current_coordinates(u, system) - neighbor_coords = current_coordinates(u, neighbor_system) + u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) + neighbor_coords = current_coordinates(u_neighbor_system, neighbor_system) foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_coords, semi; points=each_integrated_particle(system)) do particle, neighbor, @@ -574,11 +577,13 @@ function pressure_update(system, pressure, reference_density, a_ii, sum_term, om pressure[particle] = zero(pressure[particle]) end # Calculate the average density error for the termination condition - if (pressure[particle] != 0.0) + if pressure[particle] != 0.0 new_density = a_ii[particle] * pressure[particle] + sum_term[particle] - iisph_source_term(system, particle) + reference_density - density_error[particle] = (new_density - reference_density) + density_error[particle] = abs(new_density - reference_density) + else + density_error[particle] = zero(eltype(density_error)) end end relative_density_error = sum(density_error) / reference_density @@ -736,9 +741,11 @@ end function check_configuration(system::ImplicitIncompressibleSPHSystem, systems, nhs) (; time_step, omega) = system foreach_system(systems) do neighbor - if neighbor isa WeaklyCompressibleSPHSystem - throw(ArgumentError("`ImplicitIncompressibleSPHSystem` cannot be used together with - `WeaklyCompressibleSPHSystem`")) + if neighbor isa WeaklyCompressibleSPHSystem || + neighbor isa EntropicallyDampedSPHSystem + neighbor_name = neighbor |> typeof |> nameof + throw(ArgumentError("`ImplicitIncompressibleSPHSystem` cannot be used " * + "together with `$neighbor_name`")) end if neighbor isa WallBoundarySystem if (neighbor.boundary_model isa BoundaryModelDummyParticles && diff --git a/src/schemes/fluid/surface_normal_sph.jl b/src/schemes/fluid/surface_normal_sph.jl index 4db94ea763..7334f58c75 100644 --- a/src/schemes/fluid/surface_normal_sph.jl +++ b/src/schemes/fluid/surface_normal_sph.jl @@ -1,26 +1,140 @@ @doc raw""" ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) + ideal_density_threshold=0.0, interface_taper_start=0.8, + support_taper_width=0.025) Color field based computation of the interface normals. # Keywords - `boundary_contact_threshold=0.1`: If this threshold is reached the fluid is assumed to be in contact with the boundary. - `interface_threshold=0.01`: Threshold for normals to be removed as being invalid. -- `ideal_density_threshold=0.0`: Assume particles are inside if they are above this threshold, which is relative to the `ideal_neighbor_count`. +- `ideal_density_threshold=0.0`: For Morris CSF, assume particles are inside when their + continuous kernel-support moment is above this fraction of + complete support. Zero disables this filter. Other models + retain their existing neighbor-count interpretation. +- `interface_taper_start=0.8`: Start Morris CSF interface activation at this fraction of + `interface_threshold`. +- `support_taper_width=0.025`: Width of the Morris CSF support-moment transition above + `ideal_density_threshold`. """ struct ColorfieldSurfaceNormal{ELTYPE} boundary_contact_threshold::ELTYPE interface_threshold::ELTYPE ideal_density_threshold::ELTYPE + interface_taper_start::ELTYPE + support_taper_width::ELTYPE end -function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, - ideal_density_threshold=0.0) - return ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, +@doc raw""" + CorrectedCSFSurfaceNormal() + +Interface geometry for the corrected continuous-surface-force (C-CSF) method of Vergnaud +et al. (2022). The outward unit normal is computed from the renormalized gradient of the +smallest eigenvalue of the first-order kernel moment. Curvature uses the corresponding +renormalized divergence with the published thin-jet angular filter, and the surface delta +uses the published Shepard correction. + +This explicit opt-in implements the single-fluid free-surface core (equations 15--25) with +[`SurfaceTensionMorris`](@ref). Boundary-integral and contact-angle terms are not included. +""" +struct CorrectedCSFSurfaceNormal end + +@inline validate_corrected_csf(surface_normal_method, surface_tension) = nothing + +function validate_corrected_csf(::CorrectedCSFSurfaceNormal, surface_tension) + surface_tension isa SurfaceTensionMorris || + throw(ArgumentError("`CorrectedCSFSurfaceNormal` requires `SurfaceTensionMorris`")) + return nothing +end + +function ColorfieldSurfaceNormal(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + return ColorfieldSurfaceNormal(; boundary_contact_threshold, interface_threshold, ideal_density_threshold) end +function ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01, + ideal_density_threshold=0.0, interface_taper_start=0.8, + support_taper_width=0.025) + if !(boundary_contact_threshold isa Real) || isnan(boundary_contact_threshold) || + boundary_contact_threshold < 0 + throw(ArgumentError("`boundary_contact_threshold` must be non-negative and not NaN")) + end + if !(interface_threshold isa Real) || !isfinite(interface_threshold) || + interface_threshold < 0 + throw(ArgumentError("`interface_threshold` must be finite and non-negative")) + end + if !(ideal_density_threshold isa Real) || !isfinite(ideal_density_threshold) || + ideal_density_threshold < 0 + throw(ArgumentError("`ideal_density_threshold` must be finite and non-negative")) + end + if !(interface_taper_start isa Real) || !isfinite(interface_taper_start) || + !(0 <= interface_taper_start < 1) + throw(ArgumentError("`interface_taper_start` must be finite and in [0, 1)")) + end + if !(support_taper_width isa Real) || !isfinite(support_taper_width) || + support_taper_width <= 0 + throw(ArgumentError("`support_taper_width` must be finite and positive")) + end + + thresholds = promote(boundary_contact_threshold, interface_threshold, + ideal_density_threshold) + ELTYPE = typeof(first(thresholds)) + if ELTYPE <: Integer + thresholds = float.(thresholds) + ELTYPE = typeof(first(thresholds)) + end + + taper_start = convert(ELTYPE, interface_taper_start) + taper_width = convert(ELTYPE, support_taper_width) + return ColorfieldSurfaceNormal(thresholds..., taper_start, taper_width) +end + +@inline function cubic_smoothstep(value) + value <= zero(value) && return zero(value) + value >= one(value) && return one(value) + return value^2 * (3 - 2value) +end + +@inline function gradient_interface_activity(normal_norm, support_radius, + surface_normal_method::ColorfieldSurfaceNormal) + threshold = surface_normal_method.interface_threshold + dimensionless_norm = support_radius * normal_norm + if iszero(threshold) + return iszero(dimensionless_norm) ? zero(dimensionless_norm) : + one(dimensionless_norm) + end + + lower_bound = surface_normal_method.interface_taper_start * threshold + transition_coordinate = (dimensionless_norm - lower_bound) / + (threshold - lower_bound) + return cubic_smoothstep(transition_coordinate) +end + +@inline function support_interface_activity(support_moment, + surface_normal_method::ColorfieldSurfaceNormal) + threshold = surface_normal_method.ideal_density_threshold + iszero(threshold) && return one(support_moment) + + transition_coordinate = (support_moment - threshold) / + surface_normal_method.support_taper_width + return one(support_moment) - cubic_smoothstep(transition_coordinate) +end + +@inline function surface_interface_activity(system, particle) + return surface_interface_activity(surface_tension_model(system), system, particle) +end + +@inline function surface_interface_activity(::SurfaceTensionMorris, system, particle) + return @inbounds system.cache.interface_activity[particle] +end + +@inline function surface_interface_activity(surface_tension, system, particle) + normal = surface_normal(system, particle) + return dot(normal, normal) > eps(eltype(normal)) ? one(eltype(normal)) : + zero(eltype(normal)) +end + function create_cache_surface_normal(surface_normal_method, ELTYPE, NDIMS, nparticles) return (;) end @@ -33,6 +147,20 @@ function create_cache_surface_normal(::ColorfieldSurfaceNormal, ELTYPE, NDIMS, n return (; surface_normal, neighbor_count, colorfield, correction_factor) end +function create_cache_surface_normal(::CorrectedCSFSurfaceNormal, ELTYPE, NDIMS, nparticles) + surface_normal = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + neighbor_count = Array{ELTYPE, 1}(undef, nparticles) + correction_factor = Array{ELTYPE, 1}(undef, nparticles) + ccsf_correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, nparticles) + ccsf_minimum_eigenvalue = Array{ELTYPE, 1}(undef, nparticles) + ccsf_lambda_gradient = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + ccsf_color_gradient = Array{ELTYPE, 2}(undef, NDIMS, nparticles) + ccsf_shepard_sum = Array{ELTYPE, 1}(undef, nparticles) + return (; surface_normal, neighbor_count, correction_factor, + ccsf_correction_matrix, ccsf_minimum_eigenvalue, + ccsf_lambda_gradient, ccsf_color_gradient, ccsf_shepard_sum) +end + @inline function surface_normal(particle_system::AbstractFluidSystem, particle) (; cache) = particle_system return extract_svector(cache.surface_normal, particle_system, particle) @@ -67,6 +195,8 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui for i in 1:ndims(system) cache.surface_normal[i, particle] += m_b / density_neighbor * grad_kernel[i] end + accumulate_surface_support_moment!(system, surface_tension_model(system), particle, + m_b / density_neighbor, pos_diff, grad_kernel) cache.neighbor_count[particle] += 1 end @@ -74,11 +204,45 @@ function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractFlui return system end +@inline function accumulate_surface_support_moment!(system, surface_tension, particle, + volume, pos_diff, grad_kernel) + return system +end + +@inline function accumulate_surface_support_moment!(system, ::SurfaceTensionMorris, + particle, volume, pos_diff, + grad_kernel) + value = -volume * dot(pos_diff, grad_kernel) / ndims(system) + @inbounds system.cache.support_moment[particle] += value + return system +end + +@inline function accumulate_boundary_surface_support_moment!(system, surface_tension, + neighbor_system, + v_neighbor_system, particle, + neighbor, pos_diff, distance) + return system +end + +@inline function accumulate_boundary_surface_support_moment!(system, + ::SurfaceTensionMorris, + neighbor_system, + v_neighbor_system, particle, + neighbor, pos_diff, distance) + m_b = hydrodynamic_mass(neighbor_system, neighbor) + density_neighbor = current_density(v_neighbor_system, neighbor_system, neighbor) + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + value = -m_b / density_neighbor * dot(pos_diff, grad_kernel) / ndims(system) + @inbounds system.cache.support_moment[particle] += value + return system +end + # Section 2.2 in Akinci et al. 2013 "Versatile Surface Tension and Adhesion for SPH Fluids" # Note: This is the simplest form of normal approximation commonly used in SPH and comes # with serious deficits in accuracy especially at corners, small neighborhoods and boundaries function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_system, v, - u_neighbor_system, semi, surface_normal_method) + v_neighbor_system, u_neighbor_system, semi, + surface_normal_method) (; cache) = system (; colorfield, initial_colorfield) = neighbor_system.boundary_model.cache (; boundary_contact_threshold) = surface_normal_method @@ -107,6 +271,10 @@ function calc_boundary_normal!(system::AbstractFluidSystem, neighbor_system, u_s foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_system_coords, semi) do particle, neighbor, pos_diff, distance + accumulate_boundary_surface_support_moment!(system, surface_tension_model(system), + neighbor_system, v_neighbor_system, + particle, neighbor, pos_diff, distance) + # We assume that we are in contact with the boundary if the color of the boundary particle # is larger than the threshold if colorfield[neighbor] / maximum_colorfield > boundary_contact_threshold @@ -126,8 +294,8 @@ end function calc_normal!(system::AbstractFluidSystem, neighbor_system::AbstractBoundarySystem, u_system, v, v_neighbor_system, u_neighbor_system, semi, surface_normal_method, neighbor_surface_normal_method) - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) + return calc_boundary_normal!(system, neighbor_system, u_system, v, v_neighbor_system, + u_neighbor_system, semi, surface_normal_method) end function remove_invalid_normals!(system::AbstractFluidSystem, surface_tension, @@ -147,8 +315,7 @@ end # See Morris 2000 "Simulating surface tension with smoothed particle hydrodynamics" function remove_invalid_normals!(system::AbstractFluidSystem, - surface_tension::Union{SurfaceTensionMorris, - SurfaceTensionMomentumMorris}, + surface_tension::SurfaceTensionMomentumMorris, surface_normal_method::ColorfieldSurfaceNormal) (; cache, smoothing_kernel) = system (; ideal_density_threshold, interface_threshold) = surface_normal_method @@ -188,6 +355,54 @@ function remove_invalid_normals!(system::AbstractFluidSystem, return system end +function remove_invalid_normals!(system::AbstractFluidSystem, + ::SurfaceTensionMorris, + surface_normal_method::ColorfieldSurfaceNormal) + (; cache, smoothing_kernel) = system + support_radius = compact_support(smoothing_kernel, initial_smoothing_length(system)) + + for particle in each_integrated_particle(system) + cache.delta_s[particle] = zero(eltype(system)) + cache.interface_activity[particle] = zero(eltype(system)) + + particle_surface_normal = surface_normal(system, particle) + norm2 = dot(particle_surface_normal, particle_surface_normal) + if !(norm2 > eps(norm2)) + cache.surface_normal[1:ndims(system), particle] .= 0 + continue + end + + normal_norm = sqrt(norm2) + gradient_activity = gradient_interface_activity(normal_norm, support_radius, + surface_normal_method) + support_moment = cache.support_moment[particle] + support_activity = support_interface_activity(support_moment, + surface_normal_method) + activity = gradient_activity * support_activity + if !(activity > zero(activity)) + cache.surface_normal[1:ndims(system), particle] .= 0 + continue + end + + cache.interface_activity[particle] = activity + # A one-phase free surface samples one half of the kernel-smoothed interface. + cache.delta_s[particle] = 2 * normal_norm * activity + cache.surface_normal[1:ndims(system), + particle] = particle_surface_normal / normal_norm + end + + return system +end + +@inline reset_surface_interface_data!(system, surface_tension) = system + +@inline function reset_surface_interface_data!(system, ::SurfaceTensionMorris) + set_zero!(system.cache.support_moment) + set_zero!(system.cache.interface_activity) + set_zero!(system.cache.delta_s) + return system +end + function compute_surface_normal!(system, surface_normal_method, v, u, v_ode, u_ode, semi, t) return system end @@ -200,6 +415,7 @@ function compute_surface_normal!(system::AbstractFluidSystem, # Reset surface normal set_zero!(cache.surface_normal) set_zero!(cache.neighbor_count) + reset_surface_interface_data!(system, surface_tension) # TODO: if color values are set only different systems need to be called @trixi_timeit timer() "compute surface normal" foreach_system(semi) do neighbor_system @@ -215,6 +431,130 @@ function compute_surface_normal!(system::AbstractFluidSystem, return system end +@inline function ccsf_store_matrix!(matrix_cache, system, particle, matrix) + for column in 1:ndims(system), row in 1:ndims(system) + @inbounds matrix_cache[row, column, particle] = matrix[row, column] + end + return matrix_cache +end + +@inline function ccsf_minimum_eigenvalue(matrix) + # `eigmin` falls back to an allocating dense eigensolver for static matrices. + symmetric_matrix = (matrix + transpose(matrix)) / 2 + return minimum(eigvals(Symmetric(symmetric_matrix))) +end + +@inline function ccsf_corrected_divergence(normal_difference, renormalization, + kernel_direction) + return dot(renormalization * normal_difference, kernel_direction) +end + +@inline function ccsf_lambda_difference(lambda_i, lambda_j) + return lambda_i >= oftype(lambda_i, 0.7) ? lambda_j - lambda_i : lambda_j +end + +function compute_surface_normal!(system::AbstractFluidSystem, + ::CorrectedCSFSurfaceNormal, + v, u, v_ode, u_ode, semi, t) + system.surface_tension isa SurfaceTensionMorris || + throw(ArgumentError("`CorrectedCSFSurfaceNormal` requires `SurfaceTensionMorris`")) + cache = system.cache + matrix_cache = cache.ccsf_correction_matrix + lambda = cache.ccsf_minimum_eigenvalue + lambda_gradient = cache.ccsf_lambda_gradient + color_gradient = cache.ccsf_color_gradient + shepard_sum = cache.ccsf_shepard_sum + coordinates = current_coordinates(u, system) + + set_zero!(cache.surface_normal) + set_zero!(cache.neighbor_count) + set_zero!(matrix_cache) + set_zero!(lambda) + set_zero!(lambda_gradient) + set_zero!(color_gradient) + set_zero!(shepard_sum) + set_zero!(cache.support_moment) + + @trixi_timeit timer() "compute C-CSF moments" begin + foreach_point_neighbor(system, system, coordinates, coordinates, semi; + points=each_integrated_particle(system)) do particle, + neighbor, + pos_diff, + distance + m_b = hydrodynamic_mass(system, neighbor) + rho_b = current_density(v, system, neighbor) + volume_b = m_b / rho_b + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + kernel = smoothing_kernel(system, distance, particle) + + moment = -volume_b * grad_kernel * permutedims(pos_diff) + for column in 1:ndims(system), row in 1:ndims(system) + @inbounds matrix_cache[row, column, particle] += moment[row, column] + end + for dimension in 1:ndims(system) + @inbounds color_gradient[dimension, + particle] += volume_b * grad_kernel[dimension] + end + @inbounds shepard_sum[particle] += volume_b * kernel + @inbounds cache.neighbor_count[particle] += 1 + end + end + + @threaded semi for particle in each_integrated_particle(system) + inverse_renormalization = extract_smatrix(matrix_cache, system, particle) + @inbounds lambda[particle] = ccsf_minimum_eigenvalue(inverse_renormalization) + renormalization = abs(det(inverse_renormalization)) < 1.0f-9 ? + one(inverse_renormalization) : inv(inverse_renormalization) + ccsf_store_matrix!(matrix_cache, system, particle, renormalization) + end + + @trixi_timeit timer() "compute C-CSF normal" begin + foreach_point_neighbor(system, system, coordinates, coordinates, semi; + points=each_integrated_particle(system)) do particle, + neighbor, + pos_diff, + distance + rho_b = current_density(v, system, neighbor) + volume_b = hydrodynamic_mass(system, neighbor) / rho_b + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + renormalization = extract_smatrix(matrix_cache, system, particle) + lambda_a = @inbounds lambda[particle] + lambda_b = @inbounds lambda[neighbor] + coefficient = ccsf_lambda_difference(lambda_a, lambda_b) + contribution = coefficient * volume_b * renormalization * grad_kernel + for dimension in 1:ndims(system) + @inbounds lambda_gradient[dimension, particle] += contribution[dimension] + end + end + end + + set_zero!(cache.interface_activity) + set_zero!(cache.delta_s) + for particle in each_integrated_particle(system) + gradient = extract_svector(lambda_gradient, system, particle) + gradient_norm = norm(gradient) + lambda_i = @inbounds lambda[particle] + threshold = oftype(lambda_i, 0.1) * lambda_i / + smoothing_length(system, particle) + if gradient_norm > threshold + normal = -gradient / gradient_norm + for dimension in 1:ndims(system) + @inbounds cache.surface_normal[dimension, particle] = normal[dimension] + end + @inbounds cache.interface_activity[particle] = one(lambda_i) + end + + raw_gradient = extract_svector(color_gradient, system, particle) + shepard = @inbounds shepard_sum[particle] + correction = shepard > eps(shepard) ? + max(one(shepard), inv(2shepard)) : one(shepard) + @inbounds cache.delta_s[particle] = 2correction * norm(raw_gradient) + @inbounds cache.support_moment[particle] = lambda_i + end + + return system +end + function calc_curvature!(system, neighbor_system, u_system, v, v_neighbor_system, u_neighbor_system, semi, surface_normal_method, neighbor_surface_normal_method) @@ -231,8 +571,6 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF system_coords = current_coordinates(u_system, system) neighbor_system_coords = current_coordinates(u_neighbor_system, neighbor_system) - set_zero!(correction_factor) - foreach_point_neighbor(system, neighbor_system, system_coords, neighbor_system_coords, semi) do particle, neighbor, pos_diff, distance @@ -241,28 +579,72 @@ function calc_curvature!(system::AbstractFluidSystem, neighbor_system::AbstractF n_a = surface_normal(system, particle) n_b = surface_normal(neighbor_system, neighbor) v_b = m_b / rho_b + activity_a = surface_interface_activity(system, particle) + activity_b = surface_interface_activity(neighbor_system, neighbor) - # Eq. 22: we can test against `eps()` here since the surface normals that are invalid have been removed - if dot(n_a, n_a) > eps() && dot(n_b, n_b) > eps() + if activity_a > zero(activity_a) && activity_b > zero(activity_b) w = smoothing_kernel(system, distance, particle) grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + weighted_volume = v_b * activity_b for i in 1:ndims(system) - curvature[particle] += v_b * (n_b[i] - n_a[i]) * grad_kernel[i] + curvature[particle] += weighted_volume * (n_b[i] - n_a[i]) * + grad_kernel[i] end - # Eq. 24 - correction_factor[particle] += v_b * w + correction_factor[particle] += weighted_volume * w end end - # Eq. 23 - for particle in each_integrated_particle(system) - curvature[particle] /= (correction_factor[particle] + eps()) - end + return system +end +function calc_curvature!(system::AbstractFluidSystem, + neighbor_system::AbstractFluidSystem, + u_system, v, v_neighbor_system, u_neighbor_system, semi, + ::CorrectedCSFSurfaceNormal, + ::CorrectedCSFSurfaceNormal) + system === neighbor_system || + throw(ArgumentError("`CorrectedCSFSurfaceNormal` currently supports one fluid system")) + cache = system.cache + coordinates = current_coordinates(u_system, system) + cosine_threshold = -inv(convert(eltype(system), ndims(system))) + + foreach_point_neighbor(system, system, coordinates, coordinates, semi; + points=each_integrated_particle(system)) do particle, neighbor, + pos_diff, distance + n_a = surface_normal(system, particle) + n_b = surface_normal(system, neighbor) + dot(n_a, n_a) > eps(eltype(n_a)) || return + dot(n_b, n_b) > eps(eltype(n_b)) || return + dot(n_a, n_b) >= cosine_threshold || return + + rho_b = current_density(v, system, neighbor) + volume_b = hydrodynamic_mass(system, neighbor) / rho_b + grad_kernel = smoothing_kernel_grad(system, pos_diff, distance, particle) + renormalization = extract_smatrix(cache.ccsf_correction_matrix, system, particle) + @inbounds cache.curvature[particle] += volume_b * + ccsf_corrected_divergence(n_b - n_a, + renormalization, + grad_kernel) + end return system end +@inline function normalized_surface_curvature(curvature_numerator, denominator) + denominator > sqrt(eps(typeof(denominator))) || return zero(curvature_numerator) + return curvature_numerator / denominator +end + +@inline function finalize_surface_curvature(curvature_numerator, denominator, + surface_normal_method) + return normalized_surface_curvature(curvature_numerator, denominator) +end + +@inline function finalize_surface_curvature(curvature_numerator, denominator, + ::CorrectedCSFSurfaceNormal) + return curvature_numerator +end + function compute_curvature!(system, surface_tension, v, u, v_ode, u_ode, semi, t) return system end @@ -271,17 +653,25 @@ function compute_curvature!(system::AbstractFluidSystem, surface_tension::SurfaceTensionMorris, v, u, v_ode, u_ode, semi, t) (; cache, surface_tension) = system + normal_method = surface_normal_method(system) - # Reset surface curvature + # Reset once so contributions from multiple fluid systems accumulate consistently. set_zero!(cache.curvature) + set_zero!(cache.correction_factor) @trixi_timeit timer() "compute surface curvature" foreach_system(semi) do neighbor_system u_neighbor_system = wrap_u(u_ode, neighbor_system, semi) v_neighbor_system = wrap_v(v_ode, neighbor_system, semi) calc_curvature!(system, neighbor_system, u, v, v_neighbor_system, - u_neighbor_system, semi, surface_normal_method(system), + u_neighbor_system, semi, normal_method, surface_normal_method(neighbor_system)) end + + for particle in each_integrated_particle(system) + denominator = cache.correction_factor[particle] + cache.curvature[particle] = finalize_surface_curvature(cache.curvature[particle], + denominator, normal_method) + end return system end diff --git a/src/schemes/fluid/surface_tension.jl b/src/schemes/fluid/surface_tension.jl index 5656e95e12..a19947e805 100644 --- a/src/schemes/fluid/surface_tension.jl +++ b/src/schemes/fluid/surface_tension.jl @@ -51,6 +51,11 @@ It calculates surface tension forces based on the curvature of the fluid interfa using particle normals and their divergence, making it suitable for simulating phenomena like droplet formation and capillary wave dynamics. +The one-phase color-gradient magnitude is retained as a surface delta. The local +continuum-surface-force acceleration is evaluated once per particle as +``-sigma * kappa * delta_s * n_hat / rho``. A smooth interface activity avoids discrete normal +and curvature-stencil switches when using [`ColorfieldSurfaceNormal`](@ref). + See [`surface_tension`](@ref) for more details. @@ -72,7 +77,10 @@ end function create_cache_surface_tension(::SurfaceTensionMorris, ELTYPE, NDIMS, nparticles) curvature = Array{ELTYPE, 1}(undef, nparticles) - return (; curvature) + delta_s = Array{ELTYPE, 1}(undef, nparticles) + interface_activity = Array{ELTYPE, 1}(undef, nparticles) + support_moment = Array{ELTYPE, 1}(undef, nparticles) + return (; curvature, delta_s, interface_activity, support_moment) end @doc raw""" @@ -229,18 +237,26 @@ end particle, neighbor, pos_diff, distance, rho_a, rho_b, grad_kernel, surface_tension_correction) - (; surface_tension_coefficient) = surface_tension_a - - # No surface tension with oneself. See `src/general/smoothing_kernels.jl` for more details. - distance^2 < eps(initial_smoothing_length(particle_system)^2) && return dv_particle + # Morris CSF is a particle-local continuum force. It is added once outside the + # neighbor loop by `surface_tension_acceleration`. + return dv_particle +end - n_a = surface_normal(particle_system, particle) - curvature_a = curvature(particle_system, particle) +@inline function surface_tension_acceleration(surface_tension, particle_system, particle, + rho_a, vector_template) + return zero(vector_template) +end - dv_particle[] -= surface_tension_correction * surface_tension_coefficient / rho_a * - curvature_a * n_a +@inline function surface_tension_acceleration(surface_tension::SurfaceTensionMorris, + particle_system, particle, rho_a, + vector_template) + delta_s = @inbounds particle_system.cache.delta_s[particle] + iszero(delta_s) && return zero(vector_template) - return dv_particle + normal = surface_normal(particle_system, particle) + curvature_a = curvature(particle_system, particle) + return -surface_tension.surface_tension_coefficient / rho_a * curvature_a * delta_s * + normal end function compute_stress_tensors!(system, surface_tension, v, u, v_ode, u_ode, semi, t) diff --git a/src/schemes/fluid/viscosity.jl b/src/schemes/fluid/viscosity.jl index b39ccb0ef4..b81d637239 100644 --- a/src/schemes/fluid/viscosity.jl +++ b/src/schemes/fluid/viscosity.jl @@ -285,7 +285,7 @@ end end @doc raw""" - ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.01) + ViscosityAdamiSGS(; nu, C_S=0.1, epsilon=0.001) Viscosity model that extends the standard [Adami formulation](@ref ViscosityAdami) by incorporating a subgrid-scale (SGS) eddy viscosity via a Smagorinsky-type [Smagorinsky (1963)](@cite Smagorinsky1963) closure. @@ -325,7 +325,7 @@ This model is appropriate for turbulent flows where unresolved scales contribute # Keywords - `nu`: Standard kinematic viscosity. - `C_S`: Smagorinsky constant. -- `epsilon=0.01`: Parameter to prevent singularities +- `epsilon=0.001`: Parameter to prevent singularities """ struct ViscosityAdamiSGS{ELTYPE} nu :: ELTYPE # Kinematic viscosity [e.g., 1e-6 m²/s] @@ -447,7 +447,7 @@ This model is appropriate for turbulent flows where unresolved scales contribute # Keywords - `nu`: Standard kinematic viscosity. - `C_S`: Smagorinsky constant. -- `epsilon=0.01`: Parameter to prevent singularities +- `epsilon=0.001`: Parameter to prevent singularities """ struct ViscosityMorrisSGS{ELTYPE} nu :: ELTYPE # Kinematic viscosity [e.g., 1e-6 m²/s] diff --git a/src/schemes/fluid/weakly_compressible_sph/rhs.jl b/src/schemes/fluid/weakly_compressible_sph/rhs.jl index 836063538b..8217711dec 100644 --- a/src/schemes/fluid/weakly_compressible_sph/rhs.jl +++ b/src/schemes/fluid/weakly_compressible_sph/rhs.jl @@ -40,6 +40,11 @@ function interact!(dv, v_particle_system, u_particle_system, # inside the closure in the `foreach_neighbor` loop. dv_particle = Ref(zero(v_a)) drho_particle = Ref(zero(rho_a)) + if particle_system === neighbor_system + dv_particle[] += surface_tension_acceleration(surface_tension_a, + particle_system, particle, + rho_a, v_a) + end # Loop over all neighbors within the kernel cutoff @inbounds foreach_neighbor(system_coords, neighbor_system_coords, diff --git a/src/schemes/fluid/weakly_compressible_sph/system.jl b/src/schemes/fluid/weakly_compressible_sph/system.jl index eea0607d7d..bc3a09f81b 100644 --- a/src/schemes/fluid/weakly_compressible_sph/system.jl +++ b/src/schemes/fluid/weakly_compressible_sph/system.jl @@ -133,9 +133,10 @@ function WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, if surface_tension !== nothing && surface_normal_method === nothing surface_normal_method = ColorfieldSurfaceNormal() end + validate_corrected_csf(surface_normal_method, surface_tension) if surface_normal_method !== nothing && reference_particle_spacing < eps() - throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using `ColorfieldSurfaceNormal` or a surface tension model")) + throw(ArgumentError("`reference_particle_spacing` must be set to a positive value when using a surface-normal method or a surface tension model")) end pressure_acceleration = choose_pressure_acceleration_formulation(pressure_acceleration, diff --git a/src/schemes/structure/rigid_body/system.jl b/src/schemes/structure/rigid_body/system.jl index 033ea9f1da..604d601b8b 100644 --- a/src/schemes/structure/rigid_body/system.jl +++ b/src/schemes/structure/rigid_body/system.jl @@ -281,8 +281,9 @@ function calc_normal!(system::AbstractFluidSystem, surface_normal_method, neighbor_surface_normal_method) haskey(neighbor_system.boundary_model.cache, :initial_colorfield) || return system - return calc_boundary_normal!(system, neighbor_system, u_system, v, u_neighbor_system, - semi, surface_normal_method) + return calc_boundary_normal!(system, neighbor_system, u_system, v, + v_neighbor_system, u_neighbor_system, semi, + surface_normal_method) end @inline function adhesion_force!(dv_particle, diff --git a/src/setups/complex_shape.jl b/src/setups/complex_shape.jl index 6a78b412e2..d0da406f7d 100644 --- a/src/setups/complex_shape.jl +++ b/src/setups/complex_shape.jl @@ -52,6 +52,8 @@ function ComplexShape(geometry; particle_spacing, density, throw(ArgumentError("`WindingNumberHormann` only supports 2D geometries")) end + require_closed_geometry(geometry, "ComplexShape") + if grid_offset < 0.0 throw(ArgumentError("only a positive `grid_offset` is supported")) end @@ -91,8 +93,10 @@ of the geometry. - `boundary_density`: Density of each boundary particle. - `place_on_shell`: When `place_on_shell=true`, boundary particles will be placed one particle spacing from the surface of the geometry. - Otherwise when `place_on_shell=true` (simulating fluid particles), + Otherwise when `place_on_shell=false` (simulating fluid particles), boundary particles will be placed half particle spacing away from the surface. + Thus, `boundary_thickness` must be at least one particle spacing + for `place_on_shell=true` and half a particle spacing otherwise. # Examples @@ -111,7 +115,7 @@ boundary_sampled = sample_boundary(signed_distance_field; boundary_density=1.0, │ InitialCondition │ │ ════════════════ │ │ #dimensions: ……………………………………………… 2 │ -│ #particles: ………………………………………………… 889 │ +│ #particles: ………………………………………………… 677 │ │ particle spacing: ………………………………… 0.03 │ │ eltype: …………………………………………………………… Float64 │ │ coordinate eltype: ……………………………… Float64 │ @@ -133,10 +137,22 @@ function sample_boundary(signed_distance_field; end # Only keep the required part of the signed distance field - distance_to_boundary = zero(particle_spacing) - keep_indices = (distance_to_boundary .< distances .<= max_signed_distance) + distance_to_boundary = place_on_shell ? particle_spacing : particle_spacing / 2 + if boundary_thickness < distance_to_boundary + throw(ArgumentError("`boundary_thickness` must be at least " * + "`particle_spacing` for `place_on_shell=true` and " * + "half `particle_spacing` for `place_on_shell=false`.")) + end + + keep_indices = (distance_to_boundary .<= distances .<= boundary_thickness) + boundary_positions = positions[keep_indices] + + if isempty(boundary_positions) + throw(ArgumentError("No boundary particles were sampled. Increase " * + "`boundary_thickness` or generate a denser `SignedDistanceField`.")) + end - boundary_coordinates = stack(positions[keep_indices]) + boundary_coordinates = stack(boundary_positions) return InitialCondition(; coordinates=boundary_coordinates, density=boundary_density, particle_spacing) end diff --git a/src/setups/extrude_geometry.jl b/src/setups/extrude_geometry.jl index 460583c38e..43455678bb 100644 --- a/src/setups/extrude_geometry.jl +++ b/src/setups/extrude_geometry.jl @@ -87,6 +87,14 @@ shape = extrude_geometry(shape; direction, particle_spacing=0.1, n_extrude=4, de function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::Integer, velocity=zeros(length(direction)), place_on_shell=false, mass=nothing, density=nothing, pressure=0.0) + if all(iszero, direction) + throw(ArgumentError("`direction` needs to be non-zero")) + end + + if n_extrude < 1 + throw(ArgumentError("`n_extrude` needs to be positive")) + end + direction_ = normalize(direction) NDIMS = length(direction_) @@ -105,6 +113,11 @@ function extrude_geometry(geometry; particle_spacing=-1, direction, n_extrude::I face_coords = sample_plane(geometry, particle_spacing; place_on_shell) + if size(face_coords, 1) != NDIMS + throw(ArgumentError("`direction` must be of length $(size(face_coords, 1)) " * + "for the sampled geometry")) + end + coords = (face_coords .+ i * particle_spacing * direction_ for i in 0:(n_extrude - 1)) # In this context, `stack` is faster than `hcat(coords...)` @@ -235,6 +248,10 @@ end function shift_plane_corners(plane_points::NTuple{2}, direction, particle_spacing, place_on_shell) + if length(direction) != 2 + throw(ArgumentError("`direction` must be 2D when extruding 2D points")) + end + # With `place_on_shell`, particles need to be AT the min coordinates and not half a particle # spacing away from it. (place_on_shell) && (return plane_points) @@ -254,6 +271,10 @@ end function shift_plane_corners(plane_points::NTuple{3}, direction, particle_spacing, place_on_shell) + if length(direction) != 3 + throw(ArgumentError("`direction` must be 3D when extruding 3D points")) + end + # With `place_on_shell`, particles need to be AT the min coordinates and not half a particle # spacing away from it. (place_on_shell) && (return plane_points) diff --git a/src/setups/rectangular_shape.jl b/src/setups/rectangular_shape.jl index c9043a35e6..a1daff091d 100644 --- a/src/setups/rectangular_shape.jl +++ b/src/setups/rectangular_shape.jl @@ -92,7 +92,7 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord throw(ArgumentError("`min_coordinates` must be of length $NDIMS for a $(NDIMS)D problem")) end - if density !== nothing && any(density .< eps()) + if density !== nothing && !(density isa Function) && any(density .< eps()) throw(ArgumentError("`density` needs to be positive and larger than $(eps())")) end @@ -105,34 +105,50 @@ function RectangularShape(particle_spacing, n_particles_per_dimension, min_coord place_on_shell, loop_order) if !isnothing(coordinates_perturbation) - seed!(1) amplitude = coordinates_perturbation * particle_spacing - coordinates .+= rand((-amplitude):(particle_spacing * 1e-3):(amplitude), + coordinates .+= rand(MersenneTwister(1), + (-amplitude):(particle_spacing * 1e-3):(amplitude), NDIMS, n_particles) end # Allow zero acceleration with state equation, but interpret `nothing` acceleration # with state equation as a likely mistake. if acceleration isa AbstractVector || acceleration isa Tuple + if length(acceleration) != NDIMS + throw(ArgumentError("`acceleration` must be of length $NDIMS for a $(NDIMS)D problem")) + end + if pressure != 0.0 throw(ArgumentError("`pressure` cannot be used together with `acceleration` " * "and `state_equation` (hydrostatic pressure gradient)")) end if state_equation === nothing - density_fun = pressure -> density + if density === nothing + throw(ArgumentError("`density` must be specified when using " * + "`acceleration` without `state_equation`")) + end else if density !== nothing throw(ArgumentError("`density` cannot be used together with `acceleration` " * "and `state_equation` (hydrostatic pressure gradient)")) end - density_fun = pressure -> inverse_state_equation(state_equation, pressure) end # Initialize hydrostatic pressure pressure = Vector{ELTYPE}(undef, n_particles) - initialize_pressure!(pressure, particle_spacing, acceleration, - density_fun, n_particles_per_dimension, loop_order) + if state_equation === nothing && density isa Function + initialize_pressure_with_coordinate_density!(pressure, particle_spacing, + acceleration, density, + coordinates, + n_particles_per_dimension, + loop_order) + else + density_fun = state_equation === nothing ? (pressure -> density) : + (pressure -> inverse_state_equation(state_equation, pressure)) + initialize_pressure!(pressure, particle_spacing, acceleration, + density_fun, n_particles_per_dimension, loop_order) + end if state_equation !== nothing # Weakly compressible case: get density from inverse state equation @@ -223,15 +239,23 @@ function rectangular_shape_coords(particle_spacing, n_particles_per_dimension, return coordinates end -function initialize_pressure!(pressure, particle_spacing, acceleration, density_fun, - n_particles_per_dimension, loop_order) +function acceleration_dimension(acceleration) if count(a -> abs(a) > eps(), acceleration) > 1 throw(ArgumentError("hydrostatic pressure calculation is not supported with " * "diagonal acceleration")) end + return findfirst(a -> abs(a) > eps(), acceleration) +end + +function initialize_pressure!(pressure, particle_spacing, acceleration, density_fun, + n_particles_per_dimension, loop_order) # Dimension in which the acceleration is acting - accel_dim = findfirst(a -> abs(a) > eps(), acceleration) + accel_dim = acceleration_dimension(acceleration) + if accel_dim === nothing + fill!(pressure, zero(eltype(pressure))) + return pressure + end # Compute 1D pressure gradient with explicit Euler method factor = particle_spacing * abs(acceleration[accel_dim]) @@ -265,3 +289,72 @@ function initialize_pressure!(pressure, particle_spacing, acceleration, density_ pressure[particle] = pressure_1d[index_in_accel_dim] end end + +function particle_indices_by_cartesian_index(n_particles_per_dimension, loop_order) + NDIMS = length(n_particles_per_dimension) + particle_indices = Array{Int}(undef, n_particles_per_dimension) + cartesian_indices = CartesianIndices(n_particles_per_dimension) + permutation = loop_permutation(loop_order, Val(NDIMS)) + permuted_indices = permutedims(cartesian_indices, permutation) + + for particle in eachindex(permuted_indices) + particle_indices[permuted_indices[particle]] = particle + end + + return particle_indices +end + +# This is needed for `density = coords -> ...`. The pressure-dependent path above can reuse +# one 1D pressure profile for every column. Coordinate-dependent density may vary between +# columns, so each gravity-aligned column needs its own explicit Euler integration. +function initialize_pressure_with_coordinate_density!(pressure, particle_spacing, + acceleration, density_fun, + coordinates, + n_particles_per_dimension, + loop_order) + # Dimension in which the acceleration is acting + accel_dim = acceleration_dimension(acceleration) + if accel_dim === nothing + fill!(pressure, zero(eltype(pressure))) + return pressure + end + + NDIMS = length(n_particles_per_dimension) + factor = particle_spacing * abs(acceleration[accel_dim]) + particle_indices = particle_indices_by_cartesian_index(n_particles_per_dimension, + loop_order) + + accel_indices = if sign(acceleration[accel_dim]) < 0 + n_particles_per_dimension[accel_dim]:-1:1 + else + 1:n_particles_per_dimension[accel_dim] + end + surface_index = first(accel_indices) + column_starts = ntuple(dim -> dim == accel_dim ? (surface_index:surface_index) : + axes(particle_indices, dim), Val(NDIMS)) + + for column_start in CartesianIndices(column_starts) + pressure_prev = zero(eltype(pressure)) + density_prev = zero(eltype(pressure)) + for (i, accel_index) in enumerate(accel_indices) + index = ntuple(dim -> dim == accel_dim ? accel_index : column_start[dim], + Val(NDIMS)) + particle = particle_indices[index...] + coords = SVector{NDIMS, eltype(coordinates)}(ntuple(dim -> coordinates[dim, + particle], + Val(NDIMS))) + density = density_fun(coords) + + if i == 1 + pressure[particle] = 0.5factor * density + else + pressure[particle] = pressure_prev + factor * density_prev + end + + pressure_prev = pressure[particle] + density_prev = density + end + end + + return pressure +end diff --git a/src/setups/rectangular_tank.jl b/src/setups/rectangular_tank.jl index 40abd3a9ca..155dece92f 100644 --- a/src/setups/rectangular_tank.jl +++ b/src/setups/rectangular_tank.jl @@ -117,10 +117,27 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real, F, B} throw(ArgumentError("`fluid_density` needs to be positive and larger than $(eps()).")) end + if any(<(0), fluid_size_) + throw(ArgumentError("`fluid_size` dimensions need to be non-negative")) + end + + if !(n_layers isa Integer) || n_layers < 1 + throw(ArgumentError("`n_layers` needs to be a positive integer")) + end + n_layers = Int(n_layers) + + if spacing_ratio < eps() + throw(ArgumentError("`spacing_ratio` needs to be positive and larger than $(eps()).")) + end + if length(tank_size) != NDIMS throw(ArgumentError("`tank_size` must be of length $NDIMS for a $(NDIMS)D problem")) end + if any(<(0), tank_size_) + throw(ArgumentError("`tank_size` dimensions need to be non-negative")) + end + # Fluid particle data n_particles_per_dim, fluid_size_ = fluid_particles_per_dimension(fluid_size_, particle_spacing) @@ -168,7 +185,7 @@ struct RectangularTank{NDIMS, NDIMSt2, ELTYPE <: Real, F, B} # Move the tank corner in the negative coordinate directions to the desired position boundary.coordinates .+= min_coordinates - if norm(fluid_size) > eps() + if all(>(0), n_particles_per_dim) if state_equation !== nothing # Use hydrostatic pressure gradient and calculate density from inverse state # equation, so don't pass fluid density. @@ -480,14 +497,20 @@ function check_tank_overlap(fluid_size::NTuple{2}, tank_size, particle_spacing, fluid_size_x, fluid_size_y = fluid_size if tank_size[1] < fluid_size[1] - 1e-5 * particle_spacing - n_particles_x -= 1 + n_particles_x = max(0, + floor(Int, + (tank_size[1] + 1e-5 * particle_spacing) / + particle_spacing)) fluid_size_x = n_particles_x * particle_spacing @info "The fluid was overlapping.\n New fluid length in x-direction is set to $fluid_size_x." end if tank_size[2] < fluid_size[2] - 1e-5 * particle_spacing - n_particles_y -= 1 + n_particles_y = max(0, + floor(Int, + (tank_size[2] + 1e-5 * particle_spacing) / + particle_spacing)) fluid_size_y = n_particles_y * particle_spacing @info "The fluid was overlapping.\n New fluid length in y-direction is set to $fluid_size_y." @@ -502,21 +525,30 @@ function check_tank_overlap(fluid_size::NTuple{3}, tank_size, particle_spacing, fluid_size_x, fluid_size_y, fluid_size_z = fluid_size if tank_size[1] < fluid_size[1] - 1e-5 * particle_spacing - n_particles_x -= 1 + n_particles_x = max(0, + floor(Int, + (tank_size[1] + 1e-5 * particle_spacing) / + particle_spacing)) fluid_size_x = n_particles_x * particle_spacing @info "The fluid was overlapping.\n New fluid length in x-direction is set to $fluid_size_x." end if tank_size[2] < fluid_size[2] - 1e-5 * particle_spacing - n_particles_y -= 1 + n_particles_y = max(0, + floor(Int, + (tank_size[2] + 1e-5 * particle_spacing) / + particle_spacing)) fluid_size_y = n_particles_y * particle_spacing @info "The fluid was overlapping.\n New fluid length in y-direction is set to $fluid_size_y." end if tank_size[3] < fluid_size[3] - 1e-5 * particle_spacing - n_particles_z -= 1 + n_particles_z = max(0, + floor(Int, + (tank_size[3] + 1e-5 * particle_spacing) / + particle_spacing)) fluid_size_z = n_particles_z * particle_spacing @info "The fluid was overlapping.\n New fluid length in z-direction is set to $fluid_size_z." diff --git a/src/setups/sphere_shape.jl b/src/setups/sphere_shape.jl index fea261e127..b547897888 100644 --- a/src/setups/sphere_shape.jl +++ b/src/setups/sphere_shape.jl @@ -36,7 +36,8 @@ coordinate directions as `cutout_min` and `cutout_max`. - `cutout_min`: Corner in negative coordinate directions of a cuboid that is to be cut out of the sphere. - `cutout_max`: Corner in positive coordinate directions of a cuboid that is to be - cut out of the sphere. + cut out of the sphere. If the cutout has zero volume, no particles + are removed. - `place_on_shell = false`: If `place_on_shell=true`, particles will be placed on the shell of the shape. For example, the [`TotalLagrangianSPHSystem`](@ref) requires particles to be placed on the shell of the shape and @@ -114,9 +115,22 @@ function SphereShape(particle_spacing, radius, center_position, density; cutout_min_ = collect(cutout_min) cutout_max_ = collect(cutout_max) + # A zero-volume cutout means no cutout. This keeps the 2D zero default valid for + # 3D shapes while still validating dimensionality once a real cutout is requested. + has_cutout = length(cutout_min_) != length(cutout_max_) || + norm(cutout_max_ - cutout_min_) > eps() + + if has_cutout && (length(cutout_min_) != NDIMS || length(cutout_max_) != NDIMS) + throw(ArgumentError("`cutout_min` and `cutout_max` must be of length $NDIMS " * + "for a $(NDIMS)D problem")) + end + + if has_cutout && any(cutout_min_ .> cutout_max_) + throw(ArgumentError("`cutout_min` must be smaller than or equal to `cutout_max`")) + end + # Remove particles in cutout # TODO This should consider the particle radius as well - has_cutout = norm(cutout_max_ - cutout_min_) > eps() function in_cutout(particle) return has_cutout && all(cutout_min_ .<= view(coordinates, :, particle) .<= cutout_max_) diff --git a/test/general/custom_quantities.jl b/test/general/custom_quantities.jl index 9181269a6b..b5d9f63830 100644 --- a/test/general/custom_quantities.jl +++ b/test/general/custom_quantities.jl @@ -121,4 +121,104 @@ @test isnan(avg_density(boundary_system, dv_ode, du_ode, v_ode, u_ode, semi, t)) end end + + @testset "Structure kinetic energy" begin + struct EnergyStructureMock{IC, M} <: TrixiParticles.AbstractStructureSystem{2} + initial_condition::IC + mass::M + end + + Base.eltype(::EnergyStructureMock) = Float64 + TrixiParticles.compact_support(::EnergyStructureMock, neighbor) = 1.0 + function TrixiParticles.write_u0!(u0, system::EnergyStructureMock) + u0 .= system.initial_condition.coordinates + return u0 + end + function TrixiParticles.write_v0!(v0, system::EnergyStructureMock) + v0 .= system.initial_condition.velocity + return v0 + end + + coordinates = [0.0 1.0 2.0 + 0.0 0.0 0.0] + velocity = [1.0 2.0 3.0 + 4.0 5.0 6.0] + mass = [1.0, 2.0, 3.0] + ic = InitialCondition(; coordinates, velocity, mass, density=ones(3)) + system = EnergyStructureMock(ic, mass) + semi = Semidiscretization(system; neighborhood_search=nothing) + ode = semidiscretize(semi, (0.0, 1.0)) + v_ode, u_ode = ode.u0.x + dv_ode, du_ode = similar(v_ode), similar(u_ode) + + expected = sum(axes(velocity, 2)) do particle + return mass[particle] * dot(velocity[:, particle], velocity[:, particle]) / 2 + end + + @test kinetic_energy(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == expected + end + + @testset "Active particle reductions" begin + struct ActiveFluidQuantityMock{IC, M, P, B} <: TrixiParticles.AbstractFluidSystem{2} + initial_condition::IC + mass::M + pressure::P + buffer::B + surface_tension::Nothing + surface_normal_method::Nothing + end + + Base.eltype(::ActiveFluidQuantityMock) = Float64 + TrixiParticles.v_nvariables(::ActiveFluidQuantityMock) = 3 + TrixiParticles.buffer(system::ActiveFluidQuantityMock) = system.buffer + TrixiParticles.compact_support(::ActiveFluidQuantityMock, neighbor) = 1.0 + function TrixiParticles.current_velocity(v, ::ActiveFluidQuantityMock) + return view(v, 1:2, :) + end + function TrixiParticles.current_density(v, ::ContinuityDensity, + ::ActiveFluidQuantityMock) + return view(v, 3, :) + end + function TrixiParticles.current_density(v, system::ActiveFluidQuantityMock) + return TrixiParticles.current_density(v, ContinuityDensity(), system) + end + function TrixiParticles.current_pressure(v, system::ActiveFluidQuantityMock) + return system.pressure + end + function TrixiParticles.write_u0!(u0, system::ActiveFluidQuantityMock) + u0 .= system.initial_condition.coordinates + return u0 + end + function TrixiParticles.write_v0!(v0, system::ActiveFluidQuantityMock) + v0[1:2, :] .= system.initial_condition.velocity + v0[3, :] .= system.initial_condition.density + return v0 + end + + coordinates = [0.0 1.0 2.0 + 0.0 0.0 0.0] + velocity = [1.0 10.0 3.0 + 2.0 20.0 4.0] + mass = [1.0, 2.0, 4.0] + density = [10.0, 50.0, 30.0] + pressure = [100.0, 500.0, 300.0] + ic = InitialCondition(; coordinates, velocity, mass, density, pressure) + buffer = TrixiParticles.SystemBuffer(nparticles(ic), 0) + buffer.active_particle[2] = false + TrixiParticles.update_system_buffer!(buffer) + system = ActiveFluidQuantityMock(ic, mass, pressure, buffer, nothing, nothing) + semi = Semidiscretization(system; neighborhood_search=nothing) + ode = semidiscretize(semi, (0.0, 1.0)) + v_ode, u_ode = ode.u0.x + dv_ode, du_ode = similar(v_ode), similar(u_ode) + + @test total_mass(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 5.0 + @test max_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 300.0 + @test min_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 100.0 + @test avg_pressure(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 200.0 + @test max_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 30.0 + @test min_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 10.0 + @test avg_density(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 20.0 + @test kinetic_energy(system, dv_ode, du_ode, v_ode, u_ode, semi, t) == 52.5 + end end diff --git a/test/general/semidiscretization.jl b/test/general/semidiscretization.jl index 4b466ad66a..d8cd080a1a 100644 --- a/test/general/semidiscretization.jl +++ b/test/general/semidiscretization.jl @@ -141,6 +141,29 @@ @test_throws ArgumentError(error_str) Semidiscretization(fluid_system, boundary_system) end + + @testset verbose=true "Fluid Surface Tension Consistency" begin + struct FluidSurfaceMock <: TrixiParticles.AbstractFluidSystem{2} + surface_tension::Any + surface_normal_method::Any + end + + system_with_surface = FluidSurfaceMock(SurfaceTensionMorris(), + ColorfieldSurfaceNormal()) + system_with_normal = FluidSurfaceMock(nothing, ColorfieldSurfaceNormal()) + system_without_surface = FluidSurfaceMock(nothing, nothing) + + error_str = "either none or all fluid systems in a simulation need " * + "to use a surface tension model or a surface normal method." + @test_throws ArgumentError(error_str) TrixiParticles.check_configuration(system_with_surface, + (system_with_surface, + system_without_surface), + nothing) + @test_nowarn TrixiParticles.check_configuration(system_with_surface, + (system_with_surface, + system_with_normal), + nothing) + end end @testset verbose=true "`show`" begin diff --git a/test/preprocessing/geometries/geometries.jl b/test/preprocessing/geometries/geometries.jl index 4f428d266c..e01099dbe0 100644 --- a/test/preprocessing/geometries/geometries.jl +++ b/test/preprocessing/geometries/geometries.jl @@ -54,12 +54,98 @@ end end + @testset verbose=true "Open Polygon Closure" begin + open_square = [1.0 2.0 2.0 1.0; + 1.0 1.0 2.0 2.0] + + geometry = TrixiParticles.Polygon(open_square) + + @test TrixiParticles.nfaces(geometry) == 4 + @test first(geometry.vertices) == last(geometry.vertices) + @test TrixiParticles.volume(geometry) ≈ 1.0 + + mktempdir() do dir + filename = joinpath(dir, "open_square.asc") + open(filename, "w") do io + println(io, "# ASCII") + for vertex in eachcol(open_square) + println(io, vertex[1], " ", vertex[2]) + end + end + + geometry_from_file = load_geometry(filename) + + @test TrixiParticles.nfaces(geometry_from_file) == 4 + @test first(geometry_from_file.vertices) == last(geometry_from_file.vertices) + @test TrixiParticles.volume(geometry_from_file) ≈ 1.0 + end + end + + @testset verbose=true "Closed Geometry Detection" begin + open_square = [1.0 2.0 2.0 1.0; + 1.0 1.0 2.0 2.0] + + closed_polygon = TrixiParticles.Polygon(open_square) + open_polygon = TrixiParticles.Polygon(open_square; close_curve=false) + partial_polygon = delete_faces(closed_polygon, 2) + rebuilt_closed_polygon = delete_faces(closed_polygon, Int[]) + + @test TrixiParticles.is_closed_geometry(closed_polygon) + @test TrixiParticles.is_closed_geometry(rebuilt_closed_polygon) + @test !TrixiParticles.is_closed_geometry(open_polygon) + @test !TrixiParticles.is_closed_geometry(partial_polygon) + + shape = RectangularShape(0.5, (2, 2), (1.0, 1.0), density=1.0) + @test_throws ArgumentError intersect(shape, open_polygon) + @test_throws ArgumentError setdiff(shape, open_polygon) + + file = pkgdir(TrixiParticles, "test", "preprocessing", "data") + planar_geometry = load_geometry(joinpath(file, "inflow_geometry.stl")) + closed_mesh = extrude_geometry(planar_geometry, 0.8) + open_mesh = extrude_geometry(planar_geometry, 0.8; omit_top_face=true) + + @test !TrixiParticles.is_closed_geometry(planar_geometry) + @test TrixiParticles.is_closed_geometry(closed_mesh) + @test !TrixiParticles.is_closed_geometry(open_mesh) + end + + @testset verbose=true "`delete_faces` Rebuilds Derived Data" begin + triangle = [0.0 1.0 0.5 0.0; + 0.0 0.0 0.7 0.0] + + edge_only = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle), [1, 2]) + + @test TrixiParticles.nfaces(edge_only) == 1 + @test length(edge_only.vertices) == 2 + @test length(edge_only.vertex_normals) == 1 + @test edge_only.min_corner == min.(edge_only.edge_vertices[1]...) + @test edge_only.max_corner == max.(edge_only.edge_vertices[1]...) + @test edge_only.vertex_normals[1] == (edge_only.edge_normals[1], + edge_only.edge_normals[1]) + + A = SVector(0.0, 0.0, 0.0) + B = SVector(1.0, 0.0, 0.0) + C = SVector(0.0, 1.0, 0.0) + D = SVector(1.0, 1.0, 0.0) + face_vertices = [(A, B, C), (B, D, C)] + face_normals = [SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 1.0)] + mesh = TrixiParticles.TriangleMesh(face_vertices, face_normals, [A, B, C, D]) + + mesh = TrixiParticles.delete_faces(mesh, 1) + + @test TrixiParticles.nfaces(mesh) == 1 + @test length(mesh.vertices) == 3 + @test length(mesh.edge_normals) == 3 + @test mesh.face_vertices == [face_vertices[2]] + end + @testset verbose=true "Real World Data" begin data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data") @testset verbose=true "2D" begin files = ["hexagon", "circle", "inverted_open_curve"] + close_curves = [true, true, false] n_edges = [6, 63, 240] volumes = [2.5980750000000006, 3.1363805763454, 2.6153740535469048] @@ -74,7 +160,8 @@ points = vcat((data.var"Points:0")', (data.var"Points:1")') - geometry = load_geometry(joinpath(data_dir, files[i] * ".asc")) + geometry = load_geometry(joinpath(data_dir, files[i] * ".asc"); + close_curve=close_curves[i]) @test TrixiParticles.nfaces(geometry) == n_edges[i] @@ -157,6 +244,19 @@ end end + @testset verbose=true "Degenerate Triangle Normals" begin + vertex = SVector(0.0, 0.0, 0.0) + normal = SVector(0.0, 0.0, 0.0) + + geometry = TrixiParticles.TriangleMesh([(vertex, vertex, vertex)], + [normal], [vertex, vertex, vertex]) + + @test all(iszero, geometry.vertex_normals) + @test all(iszero, geometry.edge_normals) + @test all(all(isfinite, normal) for normal in geometry.vertex_normals) + @test all(all(isfinite, normal) for normal in geometry.edge_normals) + end + @testset verbose=true "Union" begin # Build a single geometry by uniting multiple STL patches (cuboid.stl contains separate solids). # The union should produce a closed volume. @@ -202,47 +302,22 @@ omit_bottom_face=true) winding_number_factor = 0.2 - @testset verbose=true "Omit Top Face" begin - expected_min_corner = [-0.036399998962879196; 0.24624998748302457; -0.5233639197487431;;] - expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] - - ic_1 = ComplexShape(geometry_extruded_1; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_1, - winding_number_factor)) - - @test nparticles(ic_1) == 2994 - @test isapprox(maximum(ic_1.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_1.coordinates, dims=2), expected_min_corner) - end - @testset verbose=true "Omit Bottom Face" begin - expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.49336391974874305;;] - expected_max_corner = [0.38360000103712083; 1.0562499874830245; -0.07336391974874301;;] - - ic_2 = ComplexShape(geometry_extruded_2; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_2, - winding_number_factor)) - - @test nparticles(ic_2) == 2988 - @test isapprox(maximum(ic_2.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_2.coordinates, dims=2), expected_min_corner) - end - - @testset verbose=true "Omit Both" begin - expected_min_corner = [-0.0663999989628792; 0.1562499874830246; -0.5233639197487431;;] - expected_max_corner = [0.38360000103712083; 1.1462499874830245; -0.07336391974874301;;] - - ic_3 = ComplexShape(geometry_extruded_3; particle_spacing=0.03, density=1.0, - point_in_geometry_algorithm=WindingNumberJacobson(; - geometry=geometry_extruded_3, - winding_number_factor)) - - @test nparticles(ic_3) == 3258 - @test isapprox(maximum(ic_3.coordinates, dims=2), expected_max_corner) - @test isapprox(minimum(ic_3.coordinates, dims=2), expected_min_corner) - end + @test_throws ArgumentError ComplexShape(geometry_extruded_1; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_1, + winding_number_factor)) + @test_throws ArgumentError ComplexShape(geometry_extruded_2; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_2, + winding_number_factor)) + @test_throws ArgumentError ComplexShape(geometry_extruded_3; + particle_spacing=0.03, density=1.0, + point_in_geometry_algorithm=WindingNumberJacobson(; + geometry=geometry_extruded_3, + winding_number_factor)) end end diff --git a/test/preprocessing/packing/nhs_faces.jl b/test/preprocessing/packing/nhs_faces.jl index 292522a91b..a9b038f2ed 100644 --- a/test/preprocessing/packing/nhs_faces.jl +++ b/test/preprocessing/packing/nhs_faces.jl @@ -4,7 +4,7 @@ 0.0 0.0 0.7 0.0] # Only use the third edge of the triangle, i.e. the edge from [0.1, 0.0] to [0.0, 0.0] - edge_aligned = deleteat!(TrixiParticles.Polygon(triangle), [1, 2]) + edge_aligned = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle), [1, 2]) edge_id = 1 # Only one edge in `Polygon` cell_sizes = [1.0 + sqrt(eps()), 0.1] @@ -27,7 +27,8 @@ end # Only use the first edge of the triangle, i.e. the edge from [0.0, 0.0] to [0.5, 0.7] - edge_arbitrary = deleteat!(TrixiParticles.Polygon(triangle), [2, 3]) + edge_arbitrary = TrixiParticles.delete_faces(TrixiParticles.Polygon(triangle), + [2, 3]) edge_id = 1 # Only one edge in `Polygon` expected_ncells_bbox = [(1, 1), (6, 7)] diff --git a/test/preprocessing/packing/signed_distance.jl b/test/preprocessing/packing/signed_distance.jl index cb02753b26..d49ae57fcb 100644 --- a/test/preprocessing/packing/signed_distance.jl +++ b/test/preprocessing/packing/signed_distance.jl @@ -44,6 +44,16 @@ @test repr("text/plain", signed_distance_field) == show_box end + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + + @test SignedDistanceField(geometry, 0.1) isa SignedDistanceField + @test_throws ArgumentError SignedDistanceField(geometry, 0.1; + use_for_boundary_packing=true) + end + @testset verbose=true "Real World Data" begin data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data") @@ -135,4 +145,15 @@ end end end + + @testset verbose=true "Point Matrix Input" begin + data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") + geometry = load_geometry(joinpath(data_dir, "hexagon.asc")) + + point = first(geometry.vertices) + signed_distance_field = SignedDistanceField(geometry, 0.1; points=hcat(point)) + + @test signed_distance_field.positions == [point] + @test signed_distance_field.distances == [0.0] + end end diff --git a/test/preprocessing/point_in_poly/winding_number_jacobson.jl b/test/preprocessing/point_in_poly/winding_number_jacobson.jl index ee2119a62d..3622b48f52 100644 --- a/test/preprocessing/point_in_poly/winding_number_jacobson.jl +++ b/test/preprocessing/point_in_poly/winding_number_jacobson.jl @@ -3,6 +3,11 @@ data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") geometry = load_geometry(joinpath(data_dir, "circle.asc")) + winding = WindingNumberJacobson() + + show_compact = "WindingNumberJacobson{NaiveWinding}()" + @test repr(winding) == show_compact + winding = WindingNumberJacobson(; hierarchical_winding=false) show_compact = "WindingNumberJacobson{NaiveWinding}()" @@ -30,4 +35,33 @@ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", winding) == show_box end + + @testset verbose=true "Point Matrix Input" begin + geometry = TrixiParticles.Polygon([0.0 1.0 1.0 0.0 0.0; + 0.0 0.0 1.0 1.0 0.0]) + point_storage = [0.5 0.0 1.5; + 0.5 0.0 1.5] + points = @view point_storage[:, 1:2:3] + + expected = Bool[true, false] + + inpoly_jacobson, _ = WindingNumberJacobson()(geometry, points) + inpoly_hormann, _ = WindingNumberHormann()(geometry, points) + + @test inpoly_jacobson == expected + @test inpoly_hormann == expected + end + + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + points = [SVector(0.5, 0.5)] + + jacobson = WindingNumberJacobson(; hierarchical_winding=false) + hormann = WindingNumberHormann() + + @test jacobson(geometry, points)[1] isa Vector{Bool} + @test hormann(geometry, points)[1] isa Vector{Bool} + end end diff --git a/test/schemes/boundary/dummy_particles/dummy_particles.jl b/test/schemes/boundary/dummy_particles/dummy_particles.jl index 8317bd144a..3c046dfe69 100644 --- a/test/schemes/boundary/dummy_particles/dummy_particles.jl +++ b/test/schemes/boundary/dummy_particles/dummy_particles.jl @@ -458,8 +458,8 @@ (width_reference, height_reference), (width_reference, height_reference), density; acceleration=[0.0, -9.81], - state_equation, n_layers=0, - faces=(true, true, true, false)) + state_equation, + faces=(false, false, false, false)) # Because it is a pain to deal with the linear indices of the pressure arrays, # we convert the matrices to Cartesian indices based on the coordinates. diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl index 8542fd554d..869f565f87 100644 --- a/test/schemes/boundary/open_boundary/boundary_zone.jl +++ b/test/schemes/boundary/open_boundary/boundary_zone.jl @@ -256,6 +256,24 @@ end end + @testset verbose=true "Boundary Zone 3D Float32 Tolerance" begin + edge1 = Float32[0.6208666, 0.6295315, 0.46713477] + edge2 = Float32[-0.48528308, 0.7766439, -0.4016525] + boundary_face = (zeros(Float32, 3), edge1, edge2) + face_normal = normalize(cross(edge1, edge2)) + + # This is orthogonal to Float32 precision, but not to a Float64-based tolerance. + @test abs(dot(edge1, edge2)) > sqrt(eps()) * norm(edge1) * norm(edge2) + @test abs(dot(edge1, edge2)) <= sqrt(eps(Float32)) * norm(edge1) * norm(edge2) + + boundary_zone = BoundaryZone(; boundary_face, particle_spacing=0.5f0, + face_normal, density=1.0f0, + open_boundary_layers=1, boundary_type=InFlow(), + sample_points=nothing) + + @test size(boundary_zone.initial_condition.coordinates, 2) > 0 + end + @testset verbose=true "Particle In Boundary Zone 2D" begin face_vertices = [[-0.2, -0.5], [0.3, 0.6]] face_size = face_vertices[2] - face_vertices[1] @@ -367,6 +385,26 @@ open_boundary_layers=2, boundary_type=OutFlow()) + non_orthogonal_face = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 1.0, 0.0]] + flow_direction = [0.0, 0.0, 1.0] + + error_str = "the vectors `AB` and `AC` must be orthogonal" + + @test_throws ArgumentError(error_str) BoundaryZone(; + boundary_face=non_orthogonal_face, + particle_spacing=0.1, + face_normal=flow_direction, + density=1.0, + open_boundary_layers=2, + boundary_type=InFlow()) + @test_throws ArgumentError(error_str) BoundaryZone(; + boundary_face=non_orthogonal_face, + particle_spacing=0.1, + face_normal=(-flow_direction), + density=1.0, + open_boundary_layers=2, + boundary_type=OutFlow()) + rectangular_face = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] flow_direction = [0.0, 1.0, 0.0] diff --git a/test/schemes/boundary/open_boundary/characteristic_variables.jl b/test/schemes/boundary/open_boundary/characteristic_variables.jl index 9d2e9fa03e..7e21b0efc4 100644 --- a/test/schemes/boundary/open_boundary/characteristic_variables.jl +++ b/test/schemes/boundary/open_boundary/characteristic_variables.jl @@ -18,6 +18,72 @@ # Add small offset to avoid "ArgumentError: density must be positive and larger than `eps()`" reference_density = (pos, t) -> 1000.0 * (t + sqrt(eps())) + @testset "Reject bidirectional flow" begin + initial_condition = rectangular_patch(particle_spacing, (2, 2)) + fluid_system = WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, + density_calculator=ContinuityDensity(), + state_equation=nothing) + + bidirectional = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), + face_normal=[1.0, 0.0], open_boundary_layers, + density, particle_spacing) + boundary_system = OpenBoundarySystem(bidirectional; fluid_system, buffer_size=0, + boundary_model=BoundaryModelCharacteristicsLastiwka()) + + error_str = "`BoundaryModelCharacteristicsLastiwka` needs a directed boundary zone. " * + "Please specify `InFlow()` or `OutFlow()` instead of `BidirectionalFlow()`." + @test_throws ArgumentError(error_str) Semidiscretization(fluid_system, + boundary_system) + end + + @testset "Fallback is zone-local" begin + face_vertices = ([0.0, 0.0], [0.0, 0.5]) + face_vertices_far = ([10.0, 0.0], [10.0, 0.5]) + flow_direction = SVector(1.0, 0.0) + + inflow = BoundaryZone(; boundary_face=face_vertices, face_normal=flow_direction, + open_boundary_layers, boundary_type=InFlow(), + reference_velocity, reference_pressure, reference_density, + density, particle_spacing) + inflow_far = BoundaryZone(; boundary_face=face_vertices_far, + face_normal=flow_direction, + open_boundary_layers, boundary_type=InFlow(), + reference_velocity, reference_pressure, reference_density, + density, particle_spacing) + fluid = extrude_geometry(face_vertices; particle_spacing, n_extrude=4, + density, pressure, direction=flow_direction) + fluid_system = EntropicallyDampedSPHSystem(fluid; smoothing_kernel, + smoothing_length, + sound_speed, + buffer_size=0, + density_calculator=ContinuityDensity()) + boundary_system = OpenBoundarySystem(inflow, inflow_far; + fluid_system, buffer_size=0, + boundary_model=BoundaryModelCharacteristicsLastiwka()) + semi = Semidiscretization(fluid_system, boundary_system) + ode = semidiscretize(semi, (0.0, 5.0)) + + v0_ode, u0_ode = ode.u0.x + v = TrixiParticles.wrap_v(v0_ode, boundary_system, semi) + u = TrixiParticles.wrap_u(u0_ode, boundary_system, semi) + + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, 2.0) + TrixiParticles.evaluate_characteristics!(boundary_system, + v, u, v0_ode, u0_ode, semi, 3.0) + + zone_1_particles = findall(==(1), boundary_system.boundary_zone_indices) + zone_2_particles = findall(==(2), boundary_system.boundary_zone_indices) + + @test any(!isapprox(characteristic, 0.0) + for characteristic in boundary_system.cache.characteristics[:, + zone_1_particles]) + @test all(isapprox(characteristic, 0.0) + for characteristic in boundary_system.cache.characteristics[:, + zone_2_particles]) + end + # Face vertices of open boundary face_vertices_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]] face_vertices_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]] @@ -144,7 +210,7 @@ initial_condition = rectangular_patch(particle_spacing, ntuple(_ -> 2, n_dims)) boundary_face = n_dims == 2 ? ([0.0, 0.0], [0.0, 1.0]) : - ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 1.0, 1.0]) + ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]) face_normal = n_dims == 2 ? [1.0, 0.0] : [1.0, 0.0, 0.0] inflow = BoundaryZone(; boundary_face, boundary_type=InFlow(), face_normal, open_boundary_layers=10, density=1.0, particle_spacing) diff --git a/test/schemes/boundary/open_boundary/dynamical_pressure.jl b/test/schemes/boundary/open_boundary/dynamical_pressure.jl index bd325eb690..dbe80f7406 100644 --- a/test/schemes/boundary/open_boundary/dynamical_pressure.jl +++ b/test/schemes/boundary/open_boundary/dynamical_pressure.jl @@ -123,7 +123,7 @@ initial_condition = rectangular_patch(particle_spacing, ntuple(_ -> 2, n_dims)) boundary_face = n_dims == 2 ? ([0.0, 0.0], [0.0, 1.0]) : - ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 1.0, 1.0]) + ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]) face_normal = n_dims == 2 ? [1.0, 0.0] : [1.0, 0.0, 0.0] inflow = BoundaryZone(; boundary_face, boundary_type=InFlow(), face_normal, open_boundary_layers=10, density=1.0, particle_spacing) diff --git a/test/schemes/boundary/open_boundary/mirroring.jl b/test/schemes/boundary/open_boundary/mirroring.jl index d9e41ee0a2..5bd9a5e5d7 100644 --- a/test/schemes/boundary/open_boundary/mirroring.jl +++ b/test/schemes/boundary/open_boundary/mirroring.jl @@ -548,7 +548,7 @@ initial_condition = rectangular_patch(particle_spacing, ntuple(_ -> 2, n_dims)) boundary_face = n_dims == 2 ? ([0.0, 0.0], [0.0, 1.0]) : - ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 1.0, 1.0]) + ([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]) face_normal = n_dims == 2 ? [1.0, 0.0] : [1.0, 0.0, 0.0] inflow = BoundaryZone(; boundary_face, boundary_type=InFlow(), face_normal, open_boundary_layers=10, density=1.0, particle_spacing) diff --git a/test/schemes/fluid/surface_normal_sph.jl b/test/schemes/fluid/surface_normal_sph.jl index 5eb8a81704..ce7b2154c3 100644 --- a/test/schemes/fluid/surface_normal_sph.jl +++ b/test/schemes/fluid/surface_normal_sph.jl @@ -113,9 +113,6 @@ function compute_and_test_surface_values(system, semi, ode; NDIMS=2) TrixiParticles.compute_surface_normal!(system, system.surface_normal_method, v, u, v0_ode, u0_ode, semi, 0.0) - TrixiParticles.remove_invalid_normals!(system, system.surface_tension, - system.surface_normal_method) - # After computation, check that surface normals have been computed and are not NaN or Inf @test all(isfinite, system.cache.surface_normal) @test all(isfinite, system.cache.neighbor_count) @@ -143,6 +140,94 @@ function compute_curvature!(system, semi, ode) v, u, v0_ode, u0_ode, semi, 0.0) end +@testset "Corrected C-CSF interface geometry" begin + particle_spacing = 0.05 + radius = 0.5 + reference_density = 1000.0 + smoothing_kernel = WendlandC2Kernel{2}() + smoothing_length = 1.4particle_spacing + fluid = SphereShape(particle_spacing, radius, (0.0, 0.0), reference_density; + sphere_type=RoundSphere()) + state_equation = StateEquationCole(; sound_speed=10.0, reference_density, + exponent=7) + surface_tension = SurfaceTensionMorris(; surface_tension_coefficient=1.0) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), + state_equation, surface_tension, + surface_normal_method=CorrectedCSFSurfaceNormal(), + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + cache = system.cache + active = findall(>(0), cache.interface_activity) + @test !isempty(active) + @test all(isfinite, cache.ccsf_minimum_eigenvalue) + @test 0.4 < minimum(cache.ccsf_minimum_eigenvalue) < 0.6 + @test maximum(cache.ccsf_minimum_eigenvalue) > 0.99 + @test all(isfinite, cache.surface_normal) + @test all(isfinite, cache.curvature) + @test all(>=(0), cache.delta_s) + @test all(active) do particle + dot(TrixiParticles.surface_normal(system, particle), + fluid.coordinates[:, particle]) > 0 + end + + weighted_curvature = sum(cache.curvature[active] .* cache.delta_s[active]) / + sum(cache.delta_s[active]) + @test isapprox(weighted_curvature, inv(radius); rtol=0.15) + + system_data = Dict{String, Any}() + TrixiParticles.add_system_data!(system_data, system.surface_normal_method) + @test system_data["surface_normal_method"]["model"] == + "CorrectedCSFSurfaceNormal" + + @test_throws ArgumentError WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, + smoothing_length, + density_calculator=ContinuityDensity(), + state_equation, + surface_tension=SurfaceTensionMomentumMorris(), + surface_normal_method=CorrectedCSFSurfaceNormal(), + reference_particle_spacing=particle_spacing) +end + +@testset "Corrected C-CSF 3D curvature" begin + particle_spacing = 0.05 + radius = 0.5 + reference_density = 1000.0 + smoothing_kernel = WendlandC2Kernel{3}() + smoothing_length = 1.4particle_spacing + fluid = SphereShape(particle_spacing, radius, (0.0, 0.0, 0.0), reference_density; + sphere_type=RoundSphere()) + system = WeaklyCompressibleSPHSystem(fluid; smoothing_kernel, smoothing_length, + density_calculator=ContinuityDensity(), + state_equation=StateEquationCole(; + sound_speed=10.0, + reference_density, + exponent=7), + surface_tension=SurfaceTensionMorris(; + surface_tension_coefficient=1.0), + surface_normal_method=CorrectedCSFSurfaceNormal(), + reference_particle_spacing=particle_spacing) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + TrixiParticles.update_systems_and_nhs(ode.u0.x..., semi, 0.0) + + cache = system.cache + active = findall(>(0), cache.interface_activity) + @test !isempty(active) + @test all(active) do particle + dot(TrixiParticles.surface_normal(system, particle), + fluid.coordinates[:, particle]) > 0 + end + @test minimum(cache.curvature[active]) > 0 + + weighted_curvature = sum(cache.curvature[active] .* cache.delta_s[active]) / + sum(cache.delta_s[active]) + @test isapprox(weighted_curvature, 2 / radius; rtol=0.15) +end + @testset verbose=true "Rigid Dummy Boundary Matches Wall Boundary" begin NDIMS = 2 particle_spacing = 0.2 @@ -163,7 +248,8 @@ end NDIMS, smoothing_length, smoothing_kernel, surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, ideal_density_threshold=0.9), - wall=true, walldistance=2.0, boundary_system_type=:wall) + wall=true, walldistance=particle_spacing, + boundary_system_type=:wall) rigid_system, rigid_boundary, rigid_semi, rigid_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, @@ -171,11 +257,19 @@ end NDIMS, smoothing_length, smoothing_kernel, surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, ideal_density_threshold=0.9), - wall=true, walldistance=2.0, + wall=true, walldistance=particle_spacing, boundary_system_type=:rigid) + free_system, _, free_semi, + free_ode = create_fluid_system(coordinates, velocity, mass, density, particle_spacing, + SurfaceTensionMorris(surface_tension_coefficient=0.072); + NDIMS, smoothing_length, smoothing_kernel, + surface_normal_method=ColorfieldSurfaceNormal(interface_threshold=0.1, + ideal_density_threshold=0.9)) + compute_and_test_surface_values(wall_system, wall_semi, wall_ode; NDIMS) compute_and_test_surface_values(rigid_system, rigid_semi, rigid_ode; NDIMS) + compute_and_test_surface_values(free_system, free_semi, free_ode; NDIMS) @test isapprox(rigid_boundary.boundary_model.cache.initial_colorfield, wall_boundary.boundary_model.cache.initial_colorfield, @@ -186,6 +280,14 @@ end @test isapprox(rigid_system.cache.neighbor_count, wall_system.cache.neighbor_count, rtol=sqrt(eps()), atol=sqrt(eps())) + @test all(isfinite, wall_system.cache.support_moment) + @test maximum(abs, wall_system.cache.support_moment) > 0 + @test isapprox(rigid_system.cache.support_moment, + wall_system.cache.support_moment, + rtol=sqrt(eps()), atol=sqrt(eps())) + @test maximum(abs, + wall_system.cache.support_moment - free_system.cache.support_moment) > + sqrt(eps()) end @testset verbose=true "CSS/CSF: Sphere Surface Normals" begin diff --git a/test/schemes/fluid/surface_tension.jl b/test/schemes/fluid/surface_tension.jl index 7fe8abbd97..089340f0f6 100644 --- a/test/schemes/fluid/surface_tension.jl +++ b/test/schemes/fluid/surface_tension.jl @@ -1,5 +1,81 @@ @testset verbose=true "Surface Tension" begin + @testset "smooth interface activity" begin + method = ColorfieldSurfaceNormal(; boundary_contact_threshold=1, + interface_threshold=0.1f0, + ideal_density_threshold=0.9, + interface_taper_start=0.8, + support_taper_width=0.05) + @test method isa ColorfieldSurfaceNormal{Float64} + @test method.interface_taper_start === 0.8 + @test method.support_taper_width === 0.05 + @test ColorfieldSurfaceNormal(1, 1, 0) isa ColorfieldSurfaceNormal{Float64} + @test ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1f0, + interface_threshold=0.01f0, + ideal_density_threshold=0.0f0, + interface_taper_start=0.8f0, + support_taper_width=0.025f0) isa + ColorfieldSurfaceNormal{Float32} + + for ELTYPE in (Float32, Float64) + @test TrixiParticles.cubic_smoothstep(ELTYPE(-1)) === ELTYPE(0) + @test TrixiParticles.cubic_smoothstep(ELTYPE(0)) === ELTYPE(0) + @test TrixiParticles.cubic_smoothstep(ELTYPE(0.5)) === ELTYPE(0.5) + @test TrixiParticles.cubic_smoothstep(ELTYPE(1)) === ELTYPE(1) + @test TrixiParticles.cubic_smoothstep(ELTYPE(2)) === ELTYPE(1) + + method_ = ColorfieldSurfaceNormal(; boundary_contact_threshold=ELTYPE(0.1), + interface_threshold=ELTYPE(0.1), + ideal_density_threshold=ELTYPE(0.9), + interface_taper_start=ELTYPE(0.8), + support_taper_width=ELTYPE(0.05)) + @test TrixiParticles.gradient_interface_activity(ELTYPE(0.08), one(ELTYPE), + method_) === ELTYPE(0) + @test TrixiParticles.gradient_interface_activity(ELTYPE(0.09), one(ELTYPE), + method_) ≈ ELTYPE(0.5) + @test TrixiParticles.gradient_interface_activity(ELTYPE(0.1), one(ELTYPE), + method_) === ELTYPE(1) + @test TrixiParticles.support_interface_activity(ELTYPE(0.9), method_) === + ELTYPE(1) + @test TrixiParticles.support_interface_activity(ELTYPE(0.925), method_) ≈ + ELTYPE(0.5) + @test TrixiParticles.support_interface_activity(ELTYPE(0.95), method_) === + ELTYPE(0) + + step = sqrt(eps(ELTYPE)) + derivative_at_zero = TrixiParticles.cubic_smoothstep(step) / step + derivative_at_one = (one(ELTYPE) - + TrixiParticles.cubic_smoothstep(one(ELTYPE) - step)) / step + @test abs(derivative_at_zero) < 4step + @test abs(derivative_at_one) < 4step + end + + disabled = ColorfieldSurfaceNormal(; ideal_density_threshold=0.0) + @test TrixiParticles.support_interface_activity(10.0, disabled) == 1.0 + @test TrixiParticles.normalized_surface_curvature(1.0, 0.0) == 0.0 + @test TrixiParticles.normalized_surface_curvature(1.0, eps()) == 0.0 + @test TrixiParticles.normalized_surface_curvature(2.0, 0.5) == 4.0 + + for threshold in (-1, NaN, Inf) + @test_throws ArgumentError ColorfieldSurfaceNormal(interface_threshold=threshold) + @test_throws ArgumentError ColorfieldSurfaceNormal(ideal_density_threshold=threshold) + end + for taper_start in (-0.1, 1.0, NaN, Inf) + @test_throws ArgumentError ColorfieldSurfaceNormal(; + interface_taper_start=taper_start) + end + for taper_width in (0.0, -0.1, NaN, Inf) + @test_throws ArgumentError ColorfieldSurfaceNormal(; + support_taper_width=taper_width) + end + + system_data = Dict{String, Any}() + TrixiParticles.add_system_data!(system_data, method) + @test system_data["surface_normal_method"]["interface_threshold"] ≈ 0.1 + @test system_data["surface_normal_method"]["interface_taper_start"] === 0.8 + @test system_data["surface_normal_method"]["support_taper_width"] === 0.05 + end + @testset verbose=true "`cohesion_force_akinci`" begin surface_tension = SurfaceTensionAkinci(surface_tension_coefficient=1.0) support_radius = 1.0 @@ -90,6 +166,151 @@ @test isapprox(zero[2], 0.0, atol=6e-15) end + @testset "Morris CSF local force" begin + function build_morris_system(solver, particle_count) + coordinates = zeros(2, particle_count) + coordinates[1, :] .= range(0.0; step=0.25, length=particle_count) + initial_condition = InitialCondition(; coordinates, + velocity=zeros(2, particle_count), + mass=ones(particle_count), + density=ones(particle_count), + particle_spacing=0.25) + smoothing_kernel = WendlandC2Kernel{2}() + surface_tension = SurfaceTensionMorris(; surface_tension_coefficient=0.7) + normal_method = ColorfieldSurfaceNormal(; interface_threshold=0.1) + if solver == :wcsph + return WeaklyCompressibleSPHSystem(initial_condition; smoothing_kernel, + smoothing_length=0.5, + density_calculator=ContinuityDensity(), + state_equation=StateEquationCole(; + sound_speed=10.0, + reference_density=1.0, + exponent=1), + surface_tension, + surface_normal_method=normal_method, + reference_particle_spacing=0.25) + end + return EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length=0.5, sound_speed=10.0, + density_calculator=ContinuityDensity(), + surface_tension, + surface_normal_method=normal_method, + reference_particle_spacing=0.25) + end + + function morris_rhs_effect(system) + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + system.cache.surface_normal[1, :] .= 1.0 + system.cache.surface_normal[2, :] .= 0.0 + system.cache.curvature .= 3.0 + system.cache.delta_s .= 2.0 + system.cache.interface_activity .= 1.0 + + return GC.@preserve v_ode u_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + rho_a = TrixiParticles.current_density(v, system, 1) + expected = TrixiParticles.surface_tension_acceleration(system.surface_tension, + system, 1, rho_a, + SVector(0.0, 0.0)) + with_surface_tension = zeros(eltype(v), size(v)) + TrixiParticles.interact!(with_surface_tension, v, u, v, u, + system, system, semi) + system.cache.delta_s .= 0 + without_surface_tension = zeros(eltype(v), size(v)) + TrixiParticles.interact!(without_surface_tension, v, u, v, u, + system, system, semi) + return (with_surface_tension - without_surface_tension)[1:2, :], + expected + end + end + + effects = [] + for solver in (:wcsph, :edac), particle_count in (2, 4) + effect, + expected = morris_rhs_effect(build_morris_system(solver, particle_count)) + @test all(particle -> effect[:, particle] ≈ expected, axes(effect, 2)) + push!(effects, effect[:, 1]) + end + @test all(effect -> effect ≈ first(effects), effects) + + system = build_morris_system(:wcsph, 2) + system.cache.surface_normal .= [2.0 1.0; 0.0 1.0] + system.cache.support_moment .= 0 + TrixiParticles.remove_invalid_normals!(system, system.surface_tension, + system.surface_normal_method) + @test system.cache.delta_s ≈ [4.0, 2sqrt(2)] + @test system.cache.interface_activity == [1.0, 1.0] + @test system.cache.surface_normal[:, 1] ≈ [1.0, 0.0] + @test system.cache.surface_normal[:, 2] ≈ [1 / sqrt(2), 1 / sqrt(2)] + system.cache.surface_normal[:, 1] .= [NaN, 0.0] + TrixiParticles.remove_invalid_normals!(system, system.surface_tension, + system.surface_normal_method) + @test iszero(system.cache.surface_normal[:, 1]) + @test iszero(system.cache.delta_s[1]) + @test iszero(system.cache.interface_activity[1]) + + system.cache.surface_normal[1, :] .= 1.0 + system.cache.surface_normal[2, :] .= 0.0 + system.cache.curvature .= 3.0 + system.cache.delta_s .= 2.0 + acceleration = TrixiParticles.surface_tension_acceleration(system.surface_tension, + system, 1, 1.0, + SVector(0.0, 0.0)) + @test acceleration ≈ SVector(-4.2, 0.0) + system.cache.curvature[1] /= 2 + system.cache.delta_s[1] /= 2 + scaled_acceleration = TrixiParticles.surface_tension_acceleration(system.surface_tension, + system, 1, 1.0, + SVector(0.0, + 0.0)) + @test scaled_acceleration ≈ acceleration / 4 + + semi = Semidiscretization(system) + ode = semidiscretize(semi, (0.0, 0.01)) + v_ode, u_ode = ode.u0.x + TrixiParticles.update_systems_and_nhs(v_ode, u_ode, semi, 0.0) + system.cache.surface_normal .= [1.0 0.0; 0.0 1.0] + + function curvature_with_neighbor_activity(activity) + system.cache.interface_activity .= [1.0, activity] + fill!(system.cache.curvature, 0) + fill!(system.cache.correction_factor, 0) + GC.@preserve v_ode u_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + TrixiParticles.calc_curvature!(system, system, u, v, v, u, semi, + system.surface_normal_method, + system.surface_normal_method) + end + denominator = system.cache.correction_factor[1] + return TrixiParticles.normalized_surface_curvature(system.cache.curvature[1], + denominator) + end + + curvature_zero = curvature_with_neighbor_activity(0.0) + curvature_small = curvature_with_neighbor_activity(1.0e-6) + curvature_full = curvature_with_neighbor_activity(1.0) + @test iszero(curvature_zero) + @test abs(curvature_small) < 1.0e-4 * abs(curvature_full) + @test isfinite(curvature_full) + + curvature_numerator = copy(system.cache.curvature) + correction_factor = copy(system.cache.correction_factor) + GC.@preserve v_ode u_ode begin + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + TrixiParticles.calc_curvature!(system, system, u, v, v, u, semi, + system.surface_normal_method, + system.surface_normal_method) + end + @test system.cache.curvature ≈ 2curvature_numerator + @test system.cache.correction_factor ≈ 2correction_factor + end + @testset "compute_stress_tensors! (MomentumMorris)" begin # 1. Define Minimal Initial Condition with 2 Particles in 2D coords = [0.0 1.0; diff --git a/test/setups/complex_shape.jl b/test/setups/complex_shape.jl index 2917661225..ca6b879cb1 100644 --- a/test/setups/complex_shape.jl +++ b/test/setups/complex_shape.jl @@ -2,6 +2,44 @@ data_dir = pkgdir(TrixiParticles, "examples", "preprocessing", "data") validation_dir = pkgdir(TrixiParticles, "test", "preprocessing", "data") + @testset verbose=true "Sample Boundary" begin + particle_spacing = 0.1 + positions = [ + SVector(0.0, 0.0), + SVector(0.1, 0.0), + SVector(0.2, 0.0), + SVector(0.3, 0.0), + SVector(0.4, 0.0) + ] + distances = [0.01, 0.05, 0.1, 0.2, 0.3] + + signed_distance_field = (; positions, distances, particle_spacing, + boundary_packing=true, max_signed_distance=0.3) + + boundary = sample_boundary(signed_distance_field; boundary_density=1.0, + boundary_thickness=0.2, place_on_shell=false) + @test boundary.coordinates ≈ stack(positions[2:4]) + + boundary = sample_boundary(signed_distance_field; boundary_density=1.0, + boundary_thickness=0.2, place_on_shell=true) + @test boundary.coordinates ≈ stack(positions[3:4]) + + @test_throws ArgumentError sample_boundary(signed_distance_field; + boundary_density=1.0, + boundary_thickness=0.04, + place_on_shell=false) + + too_thin_sdf = (; positions, distances, particle_spacing, + boundary_packing=true, max_signed_distance=0.1) + @test_throws ArgumentError sample_boundary(too_thin_sdf; boundary_density=1.0, + boundary_thickness=0.2) + + not_boundary_sdf = (; positions, distances, particle_spacing, + boundary_packing=false, max_signed_distance=0.3) + @test_throws ArgumentError sample_boundary(not_boundary_sdf; boundary_density=1.0, + boundary_thickness=0.2) + end + @testset verbose=true "2D" begin @testset verbose=true "Shifted Rectangle" begin algorithms = [ @@ -41,7 +79,7 @@ end @testset verbose=true "Real World Data" begin - files = ["hexagon", "circle", "inverted_open_curve"] + files = ["hexagon", "circle"] algorithms = [ WindingNumberHormann(), WindingNumberJacobson(; hierarchical_winding=false) @@ -72,7 +110,8 @@ # See https://docs.julialang.org/en/v1/base/base/#var%22name%22 coords = vcat((data.var"Points:0")', (data.var"Points:1")') - geometry = load_geometry(joinpath(data_dir, files[j] * ".asc")) + geometry = load_geometry(joinpath(data_dir, files[j] * ".asc"); + close_curve=true) shape_sampled = ComplexShape(geometry; particle_spacing=0.05, density=1.0, point_in_geometry_algorithm) @@ -82,6 +121,15 @@ end end + @testset verbose=true "Open Geometry Validation" begin + open_square = [0.0 1.0 1.0 0.0; + 0.0 0.0 1.0 1.0] + geometry = TrixiParticles.Polygon(open_square; close_curve=false) + + @test_throws ArgumentError ComplexShape(geometry; particle_spacing=0.1, + density=1.0) + end + @testset verbose=true "Intersect of Overlapping Shapes and Geometries" begin shape = RectangularShape(0.1, (10, 10), (0.0, 0.0), density=1.0) geometry = load_geometry(joinpath(data_dir, "circle.asc")) diff --git a/test/setups/extrude_geometry.jl b/test/setups/extrude_geometry.jl index 4146ab6dd6..b61b14549b 100644 --- a/test/setups/extrude_geometry.jl +++ b/test/setups/extrude_geometry.jl @@ -47,6 +47,28 @@ @test shape.coordinates ≈ expected_coords end + + @testset verbose=true "Errors" begin + point1 = [0.0, 0.0] + point2 = [0.0, 1.0] + + @test_throws ArgumentError extrude_geometry((point1, point2); + direction=[0.0, 0.0], + particle_spacing=0.1, + n_extrude=1, density=1.0) + @test_nowarn extrude_geometry((point1, point2); + direction=[1e-20, 0.0], + particle_spacing=0.1, + n_extrude=1, density=1.0) + @test_throws ArgumentError extrude_geometry((point1, point2); + direction=[1.0, 0.0], + particle_spacing=0.1, + n_extrude=0, density=1.0) + @test_throws ArgumentError extrude_geometry((point1, point2); + direction=[0.0, 0.0, 1.0], + particle_spacing=0.1, + n_extrude=1, density=1.0) + end end # 3D diff --git a/test/setups/rectangular_shape.jl b/test/setups/rectangular_shape.jl index 4e2419487f..aae2e8a290 100644 --- a/test/setups/rectangular_shape.jl +++ b/test/setups/rectangular_shape.jl @@ -48,6 +48,33 @@ @test shape.coordinates == expected_coords[i] end end + + @testset "Function Density" begin + shape = RectangularShape(0.1, (2, 1), (0.0, 0.0), + density=coords -> 1000.0 + coords[1]) + + @test shape.density ≈ [1000.05, 1000.15] + end + + @testset "Coordinates Perturbation Does Not Reset Random State" begin + Random.seed!(42) + first_random_number = rand() + next_random_number = rand() + + Random.seed!(42) + @test rand() == first_random_number + + RectangularShape(0.1, (2, 2), (0.0, 0.0), density=1.0, + coordinates_perturbation=0.1) + + @test rand() == next_random_number + end + + @testset "Errors" begin + @test_throws ArgumentError RectangularShape(0.1, (2, 2), (0.0, 0.0), + density=1000.0, + acceleration=(0.0, -9.81, 0.0)) + end end # Only show all of these nested testsets in case of errors @@ -123,6 +150,26 @@ @test shape.pressure ≈ 4.71 * 1000.0 * vec(reverse(pressure')) end end + + @testset "Function Density" begin + density_function = coords -> 1000.0 + 100coords[1] + 10coords[2] + shape = RectangularShape(particle_spacing, (2, 3), (0.0, 0.0), + density=density_function, + acceleration=(0.0, -1.0)) + + @test shape.density ≈ [1005.5, 1015.5, 1006.5, 1016.5, 1007.5, 1017.5] + @test shape.pressure ≈ [251.775, 254.275, 151.125, 152.625, 50.375, + 50.875] + @test shape.mass ≈ particle_spacing^2 * shape.density + end + + @testset "Zero Acceleration" begin + shape = RectangularShape(particle_spacing, (2, 5), (0.0, 0.0), + density=1000.0, acceleration=(0.0, 0.0)) + + @test shape.pressure == zeros(10) + @test shape.density == 1000 * ones(10) + end end # Use `@trixi_testset` to isolate the mock functions in a separate namespace @@ -186,6 +233,16 @@ shape.pressure) @test shape.mass == particle_spacing^2 * shape.density end + + @testset "Zero Acceleration" begin + shape = RectangularShape(particle_spacing, (2, 5), (0.0, 0.0); + acceleration=(0.0, 0.0), state_equation) + + @test shape.pressure == zeros(10) + @test shape.density == + TrixiParticles.inverse_state_equation.(Ref(state_equation), + shape.pressure) + end end end diff --git a/test/setups/rectangular_tank.jl b/test/setups/rectangular_tank.jl index 1a1fc88508..33775a5873 100644 --- a/test/setups/rectangular_tank.jl +++ b/test/setups/rectangular_tank.jl @@ -119,6 +119,43 @@ (water_width, water_height, 0.5), (tank_width, tank_height), water_density, spacing_ratio=3) + + error = ArgumentError("`fluid_size` dimensions need to be non-negative") + @test_throws error RectangularTank(particle_spacing, + (-water_width, water_height), + (tank_width, tank_height), + water_density) + + error = ArgumentError("`tank_size` dimensions need to be non-negative") + @test_throws error RectangularTank(particle_spacing, + (water_width, water_height), + (-tank_width, tank_height), + water_density) + + @test_throws ArgumentError RectangularTank(particle_spacing, + (water_width, water_height), + (tank_width, tank_height), + water_density, spacing_ratio=0) + + @test_throws ArgumentError RectangularTank(particle_spacing, + (water_width, water_height), + (tank_width, tank_height), + water_density, n_layers=0) + + @test_throws ArgumentError RectangularTank(particle_spacing, + (water_width, water_height), + (tank_width, tank_height), + water_density, n_layers=1.5) + + tank = RectangularTank(0.1, (1.0, 1.0), (0.3, 0.3), water_density) + @test tank.n_particles_per_dimension == (3, 3) + @test all(tank.fluid_size .≈ (0.3, 0.3)) + + tank = RectangularTank(0.1, (1.0, 1.0), (0.05, 0.3), water_density; + acceleration=(1.0, 0.0)) + @test isempty(tank.fluid.coordinates) + @test tank.n_particles_per_dimension == (0, 3) + @test all(tank.fluid_size .≈ (0.0, 0.3)) end end diff --git a/test/setups/sphere_shape.jl b/test/setups/sphere_shape.jl index c8b87d76db..42fc3e7b70 100644 --- a/test/setups/sphere_shape.jl +++ b/test/setups/sphere_shape.jl @@ -79,6 +79,18 @@ end end + @testset verbose=true "Errors" begin + @test_throws ArgumentError SphereShape(0.1, 0.5, (0.0, 0.0), 1000.0; + cutout_min=(0.2, 0.0), + cutout_max=(0.1, 0.1)) + @test_throws ArgumentError SphereShape(0.1, 0.5, (0.0, 0.0, 0.0), + 1000.0; cutout_min=(0.0, 0.0), + cutout_max=(0.1, 0.1)) + @test_nowarn SphereShape(0.1, 0.5, (0.0, 0.0, 0.0), + 1000.0; cutout_min=(0.0, 0.0), + cutout_max=(0.0, 0.0)) + end + @testset verbose=true "SphereShape 3D" begin shape_names = [ "1-particle VoxelSphere", diff --git a/test/systems/boundary_system.jl b/test/systems/boundary_system.jl index 8c55dfa029..67512a19bb 100644 --- a/test/systems/boundary_system.jl +++ b/test/systems/boundary_system.jl @@ -28,6 +28,55 @@ end end + @testset verbose=true "High-level Dummy-Particle Builder" begin + boundary_coordinates = [1.0 2.0 + 1.0 2.0] + fluid_coordinates = [0.0 0.5 + 0.0 0.0] + + boundary_ic = InitialCondition(; coordinates=boundary_coordinates, mass, density) + fluid_ic = InitialCondition(; coordinates=fluid_coordinates, mass, density) + + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.8 + state_equation = StateEquationCole(; sound_speed=15.0, reference_density=1000.0, + exponent=1) + viscosity = ViscosityAdami(nu=1e-6) + + fluid_system = WeaklyCompressibleSPHSystem(fluid_ic; + density_calculator=ContinuityDensity(), + state_equation, smoothing_kernel, + smoothing_length, + correction=KernelCorrection(), + reference_particle_spacing=0.1) + + boundary_model = BoundaryModelDummyParticles(boundary_ic; + fluid_system=fluid_system, + viscosity=viscosity) + system = WallBoundarySystem(boundary_ic, boundary_model, + adhesion_coefficient=0.3, + color_value=2) + + @test system isa WallBoundarySystem + @test system.boundary_model isa BoundaryModelDummyParticles + @test system.boundary_model.hydrodynamic_mass == boundary_ic.mass + @test system.boundary_model.density_calculator isa AdamiPressureExtrapolation + @test system.boundary_model.smoothing_kernel === smoothing_kernel + @test system.boundary_model.smoothing_length == smoothing_length + @test system.boundary_model.viscosity == viscosity + @test system.boundary_model.state_equation == state_equation + @test system.boundary_model.correction isa KernelCorrection + @test system.boundary_model.cache.reference_particle_spacing == 0.1 + @test system.adhesion_coefficient == 0.3 + @test system.cache.color == 2 + + edac_system = EntropicallyDampedSPHSystem(fluid_ic; smoothing_kernel, + smoothing_length, sound_speed=15.0) + edac_boundary_model = BoundaryModelDummyParticles(boundary_ic; + fluid_system=edac_system) + @test edac_boundary_model.state_equation === nothing + end + @testset verbose=true "Moving Boundaries" begin @testset "$(i+1)D" for i in 1:2 NDIMS = i + 1 diff --git a/test/systems/edac_system.jl b/test/systems/edac_system.jl index ce04cd6774..33316af9ae 100644 --- a/test/systems/edac_system.jl +++ b/test/systems/edac_system.jl @@ -213,6 +213,57 @@ @test v0 == vcat(velocity, [0.8, 1.0]') end + @trixi_testset "Correction cache updates" begin + coordinates = [0.0 0.1 0.0 + 0.0 0.0 0.1] + velocity = zeros(2, 3) + mass = ones(3) + density = fill(1000.0, 3) + pressure = zeros(3) + initial_condition = InitialCondition(; coordinates, velocity, mass, density, + pressure) + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.5 + + corrections = (ShepardKernelCorrection(), KernelCorrection(), GradientCorrection(), + MixedKernelGradientCorrection()) + + @testset "$(typeof(correction))" for correction in corrections + system = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, sound_speed=10.0, + correction, pressure_acceleration=nothing) + semi = Semidiscretization(system) + + TrixiParticles.initialize_neighborhood_searches!(semi) + + u_ode = vec(coordinates) + v0 = zeros(TrixiParticles.v_nvariables(system), + TrixiParticles.n_integrated_particles(system)) + TrixiParticles.write_v0!(v0, system) + v_ode = vec(v0) + + v = TrixiParticles.wrap_v(v_ode, system, semi) + u = TrixiParticles.wrap_u(u_ode, system, semi) + + for cache_key in (:kernel_correction_coefficient, :dw_gamma, + :correction_matrix) + hasproperty(system.cache, cache_key) || continue + + fill!(getproperty(system.cache, cache_key), NaN) + end + + TrixiParticles.update_quantities!(system, v, u, v_ode, u_ode, semi, 0.0) + TrixiParticles.update_pressure!(system, v, u, v_ode, u_ode, semi, 0.0) + + for cache_key in (:kernel_correction_coefficient, :dw_gamma, + :correction_matrix) + hasproperty(system.cache, cache_key) || continue + + @test all(isfinite, getproperty(system.cache, cache_key)) + end + end + end + @trixi_testset "Average Pressure" begin particle_spacing = 0.1 smoothing_kernel = SchoenbergCubicSplineKernel{2}() @@ -244,4 +295,35 @@ nparticles(system)) end end + + @trixi_testset "restart_with! with ContinuityDensity" begin + coordinates = [0.5 2.0 + 1.0 2.0] + velocity = 2 * coordinates + mass = [1.25, 1.5] + density = [990.0, 1000.0] + pressure = [5.0, 7.8] + smoothing_kernel = Val(:smoothing_kernel) + TrixiParticles.ndims(::Val{:smoothing_kernel}) = 2 + smoothing_length = 0.362 + sound_speed = 10.0 + + initial_condition = InitialCondition(; coordinates, velocity, mass, density, + pressure) + system = EntropicallyDampedSPHSystem(initial_condition; smoothing_kernel, + smoothing_length, sound_speed, + density_calculator=ContinuityDensity()) + + u_new = coordinates .+ 1 + velocity_new = velocity .+ 2 + pressure_new = [11.0, 13.0] + density_new = [980.0, 970.0] + v_new = vcat(velocity_new, pressure_new', density_new') + + TrixiParticles.restart_with!(system, v_new, u_new) + + @test system.initial_condition.coordinates == u_new + @test system.initial_condition.velocity == velocity_new + @test system.initial_condition.pressure == pressure_new + end end diff --git a/test/systems/iisph_system.jl b/test/systems/iisph_system.jl index 2bcd740552..de38f2cd2c 100644 --- a/test/systems/iisph_system.jl +++ b/test/systems/iisph_system.jl @@ -68,6 +68,7 @@ @test system.max_iterations == max_iterations @test system.time_step == time_step @test length(system.density) == size(coordinates, 2) + @test TrixiParticles.system_state_equation(system) === nothing # A too-short acceleration vector triggers dimension validation error_str1 = "`acceleration` must be of length $NDIMS for a $(NDIMS)D problem" @@ -439,12 +440,12 @@ system_pressure.predicted_density .= [990.0, 1010.0] system_pressure.sum_term .= [5.0, -2.0] system_pressure.a_ii .= [0.5, 1.0e-10] - fill!(system_pressure.density_error, 0.0) + system_pressure.density_error .= [0.0, 99.0] semi = DummySemidiscretization() # First particle uses standard Jacobi update; second hits the safeguarded zero-a_ii path. # For particle 1: (1-omega)*0 + omega/a_ii * (source - sum_term) with omega=0.4, - # source=(1000-990)=10, a_ii=0.5, sum_term=5 gives pressure 4 and density_error -3 + # source=(1000-990)=10, a_ii=0.5, sum_term=5 gives pressure 4 and abs(density_error) 3 relative_error = TrixiParticles.pressure_update(system_pressure, system_pressure.pressure, system_pressure.reference_density, @@ -454,9 +455,50 @@ system_pressure.density_error, semi) - @test isapprox(relative_error, -0.003) + @test isapprox(relative_error, 0.003) @test isapprox(system_pressure.pressure, [4.0, 0.0]) - @test isapprox(system_pressure.density_error, [-3.0, 0.0]) + @test isapprox(system_pressure.density_error, [3.0, 0.0]) + end + + @testset "Cross-system pressure sums use neighbor coordinates" begin + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.5 + time_step = 0.5 + + coordinates_a = reshape([0.0, 0.0], 2, 1) + ic_a = InitialCondition(; coordinates=coordinates_a, + velocity=zeros(2, 1), + mass=[1.0], + density=[1000.0], + pressure=[1.0]) + system = ImplicitIncompressibleSPHSystem(ic_a; + smoothing_kernel, + smoothing_length, + reference_density=1000.0, + time_step) + + coordinates_b = [0.1 0.2 + 0.0 0.0] + ic_b = InitialCondition(; coordinates=coordinates_b, + velocity=zeros(2, 2), + mass=[1.0, 1.0], + density=[1000.0, 1000.0], + pressure=[1.0, 2.0]) + neighbor_system = ImplicitIncompressibleSPHSystem(ic_b; + smoothing_kernel, + smoothing_length, + reference_density=1000.0, + time_step) + + semi = Semidiscretization(system, neighbor_system) + TrixiParticles.initialize_neighborhood_searches!(semi) + u_ode = vcat(vec(coordinates_a), vec(coordinates_b)) + u = TrixiParticles.wrap_u(u_ode, system, semi) + + @test_nowarn TrixiParticles.calculate_sum_d_ij_pj!(system.sum_d_ij_pj, + system, neighbor_system, + u, u_ode, semi) + @test !iszero(system.sum_d_ij_pj[1, 1]) end @testset "Source term and iteration limits" begin @@ -489,4 +531,30 @@ @test TrixiParticles.maximum_iisph_iterations(system_iters) == 7 end end + + @testset "Reject incompatible fluid systems" begin + smoothing_kernel = SchoenbergCubicSplineKernel{2}() + smoothing_length = 0.5 + coordinates = [0.0 0.1 + 0.0 0.2] + velocity = zeros(2, 2) + mass = [1.0, 1.0] + density = [1000.0, 1000.0] + pressure = [0.0, 0.0] + ic = InitialCondition(; coordinates, velocity, mass, density, pressure) + + iisph_system = ImplicitIncompressibleSPHSystem(ic; + smoothing_kernel, + smoothing_length, + reference_density=1000.0, + time_step=0.5) + edac_system = EntropicallyDampedSPHSystem(ic; smoothing_kernel, + smoothing_length, + sound_speed=10.0) + + error_str = "`ImplicitIncompressibleSPHSystem` cannot be used together with " * + "`EntropicallyDampedSPHSystem`" + @test_throws ArgumentError(error_str) Semidiscretization(iisph_system, + edac_system) + end end diff --git a/test/systems/open_boundary_system.jl b/test/systems/open_boundary_system.jl index d1b2b3fc1d..eca71a9d6e 100644 --- a/test/systems/open_boundary_system.jl +++ b/test/systems/open_boundary_system.jl @@ -2,14 +2,16 @@ @testset "`show`" begin # Mock fluid system - struct FluidSystemMock2 <: TrixiParticles.AbstractFluidSystem{2} + struct FluidSystemMock2{B} <: TrixiParticles.AbstractFluidSystem{2} pressure_acceleration_formulation::Nothing density_diffusion::Nothing + buffer::B end TrixiParticles.initial_smoothing_length(system::FluidSystemMock2) = 1.0 TrixiParticles.nparticles(system::FluidSystemMock2) = 1 TrixiParticles.system_smoothing_kernel(system::FluidSystemMock2) = nothing TrixiParticles.density_calculator(system::FluidSystemMock2) = TrixiParticles.ContinuityDensity() + TrixiParticles.buffer(system::FluidSystemMock2) = system.buffer inflow = BoundaryZone(; boundary_face=([0.0, 0.0], [0.0, 1.0]), particle_spacing=0.05, @@ -17,7 +19,8 @@ open_boundary_layers=4, boundary_type=InFlow()) system = OpenBoundarySystem(inflow; buffer_size=0, boundary_model=BoundaryModelCharacteristicsLastiwka(), - fluid_system=FluidSystemMock2(nothing, nothing)) + fluid_system=FluidSystemMock2(nothing, nothing, + nothing)) show_compact = "OpenBoundarySystem{2}() with 80 particles" @test repr(system) == show_compact @@ -40,7 +43,8 @@ boundary_type=OutFlow()) system = OpenBoundarySystem(outflow; buffer_size=0, boundary_model=BoundaryModelMirroringTafuni(), - fluid_system=FluidSystemMock2(nothing, nothing)) + fluid_system=FluidSystemMock2(nothing, nothing, + nothing)) show_compact = "OpenBoundarySystem{2}() with 80 particles" @test repr(system) == show_compact @@ -59,7 +63,8 @@ system = OpenBoundarySystem(outflow, inflow; buffer_size=0, boundary_model=BoundaryModelMirroringTafuni(), - fluid_system=FluidSystemMock2(nothing, nothing)) + fluid_system=FluidSystemMock2(nothing, nothing, + nothing)) show_compact = "OpenBoundarySystem{2}() with 160 particles" @test repr(system) == show_compact @@ -78,7 +83,8 @@ system = OpenBoundarySystem(outflow, inflow; buffer_size=0, boundary_model=BoundaryModelDynamicalPressureZhang(), - fluid_system=FluidSystemMock2(nothing, nothing)) + fluid_system=FluidSystemMock2(nothing, nothing, + nothing)) show_compact = "OpenBoundarySystem{2}() with 160 particles" @test repr(system) == show_compact @@ -96,6 +102,20 @@ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box + + fluid_system_with_buffer = FluidSystemMock2(nothing, nothing, + TrixiParticles.SystemBuffer(1, 3)) + system = OpenBoundarySystem(outflow; fluid_system=fluid_system_with_buffer) + @test system.boundary_model isa BoundaryModelMirroringTafuni + @test system.buffer.buffer_size == 3 + + error_str = "`buffer_size` could not be inferred for `OpenBoundarySystem` " * + "because `fluid_system` has no buffer. Pass `buffer_size=...` " * + "explicitly or construct `fluid_system` with `buffer_size=...`." + @test_throws ArgumentError(error_str) OpenBoundarySystem(outflow; + fluid_system=FluidSystemMock2(nothing, + nothing, + nothing)) end @testset "boundary zone width" begin diff --git a/test/systems/packing_system.jl b/test/systems/packing_system.jl index 20c915fb4c..7b80a723b4 100644 --- a/test/systems/packing_system.jl +++ b/test/systems/packing_system.jl @@ -41,6 +41,22 @@ └──────────────────────────────────────────────────────────────────────────────────────────────────┘""" @test repr("text/plain", system) == show_box + signed_distance_field = SignedDistanceField(geometry, 0.1; + use_for_boundary_packing=true, + max_signed_distance=0.3) + boundary_sampled = sample_boundary(signed_distance_field; boundary_density=1.0, + boundary_thickness=0.2, + place_on_shell=false) + system = ParticlePackingSystem(boundary_sampled; signed_distance_field, + background_pressure=1.0, is_boundary=true, + boundary_thickness=0.2) + @test system.shift_length == -0.25 + @test_throws ArgumentError ParticlePackingSystem(boundary_sampled; + signed_distance_field, + background_pressure=1.0, + is_boundary=true, + boundary_thickness=0.4) + system = ParticlePackingSystem(initial_condition, signed_distance_field=nothing, background_pressure=1.0)