Merge branch 'main' into tafuni-open-boundaries

.gitignore

@@ -9,6 +9,7 @@ docs/src/contributing.md
 docs/src/license.md
 docs/src/code_of_conduct.md
 docs/src/news.md
+docs/src/tutorials
 *_replaced.md
 run
 *.json

NEWS.md

@@ -3,6 +3,11 @@
 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.2.3
+
+### Highlights
+Transport Velocity Formulation (TVF) based on the work of Ramachandran et al. "Entropically damped artificial compressibility for SPH" (2019) was added.
+
 ## Version 0.2.2
 
 ### Highlights

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.2.3-dev"
+version = "0.2.4-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -21,6 +21,7 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -35,19 +36,19 @@ Adapt = "3, 4"
 CSV = "0.10"
 DataFrames = "1.6"
 DelimitedFiles = "1"
-DiffEqCallbacks = "2.25, 3"
+DiffEqCallbacks = "2, 3, 4"
 FastPow = "0.1"
 FileIO = "1"
 ForwardDiff = "0.10"
-GPUArraysCore = "0.1"
+GPUArraysCore = "0.1, 0.2"
 JSON = "0.21"
 KernelAbstractions = "0.9"
 MuladdMacro = "0.2"
 PointNeighbors = "0.4.2"
 Polyester = "0.7.10"
 RecipesBase = "1"
 Reexport = "1"
-SciMLBase = "1, 2"
+SciMLBase = "2"
 StaticArrays = "1"
 StrideArrays = "0.1"
 TimerOutputs = "0.5.25"

docs/Project.toml

@@ -2,12 +2,14 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 Documenter = "1"
 DocumenterCitations = "1"
 OrdinaryDiffEq = "6"
+Plots = "1"
 PointNeighbors = "0.4"
 TrixiBase = "0.1"

docs/make.jl

@@ -61,16 +61,23 @@ function replace_with_code(filename)
         end
     end
 
+    # Check if tutorials directory exists, else create one
+    path_tutorials = joinpath("docs", "src", "tutorials")
+    if !isdir(path_tutorials)
+        mkdir(path_tutorials)
+    end
+
+    file_basename = basename(filename)
+
     # Replace all occurrences in the markdown content
-    filename_noext, extension = splitext(filename)
-    copy_file(filename, new_file="$(filename_noext)_replaced$extension",
+    copy_file(filename, new_file=joinpath(path_tutorials, file_basename),
               pattern => replace_include)
 end
 
-replace_with_code("docs/src/tutorials/tut_setup.md")
-replace_with_code("docs/src/tutorials/tut_dam_break.md")
-replace_with_code("docs/src/tutorials/tut_beam.md")
-replace_with_code("docs/src/tutorials/tut_falling.md")
+replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_setup.md"))
+replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_dam_break.md"))
+replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_beam.md"))
+replace_with_code(joinpath("docs", "src", "tutorials_template", "tut_falling.md"))
 
 copy_file("AUTHORS.md",
           "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")

docs/src/getting_started.md

@@ -18,6 +18,11 @@ Then start the simulation by executing
 julia> trixi_include(joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl"))
 ```
 
+The easiest way to quickly visualize the result is to use Plots.jl:
+```jldoctest getting_started; filter = r".*"s, setup=:(trixi_include(@__MODULE__, joinpath(examples_dir(), "fluid", "hydrostatic_water_column_2d.jl")))
+julia> using Plots; plot(sol)
+```
+
 This will open a new window with a 2D visualization of the final solution:
 ![plot_hydrostatic_water_column](https://github.com/trixi-framework/TrixiParticles.jl/assets/44124897/95821154-577d-4323-ba57-16ef02ea24e0)
 

docs/src/refs.bib

@@ -13,6 +13,20 @@ @Article{Adami2012
 }
 
 
+@Article{Adami2013,
+  author    = {S. Adami and X.Y. Hu and N.A. Adams},
+  journal   = {Journal of Computational Physics},
+  title     = {A transport-velocity formulation for smoothed particle hydrodynamics},
+  year      = {2013},
+  month     = {may},
+  pages     = {292--307},
+  volume    = {241},
+  doi       = {10.1016/j.jcp.2013.01.043},
+  groups    = {enhancement},
+  publisher = {Elsevier {BV}},
+}
+
+
 @Article{Akinci2012,
   author    = {Akinci, Nadir and Ihmsen, Markus and Akinci, Gizem and Solenthaler, Barbara and Teschner, Matthias},
   journal   = {ACM Transactions on Graphics},

docs/src/systems/entropically_damped_sph.md

@@ -45,3 +45,64 @@ is a good choice for a wide range of Reynolds numbers (0.0125 to 10000).
 Modules = [TrixiParticles]
 Pages = [joinpath("schemes", "fluid", "entropically_damped_sph", "system.jl")]
 ```
+
+## [Transport Velocity Formulation (TVF)](@id transport_velocity_formulation)
+Standard SPH suffers from problems like tensile instability or the creation of void regions in the flow.
+To address these problems, [Adami (2013)](@cite Adami2013) modified the advection velocity and added an extra term to the momentum equation.
+The authors introduced the so-called Transport Velocity Formulation (TVF) for WCSPH. [Ramachandran (2019)](@cite Ramachandran2019) applied the TVF
+also for the [EDAC](@ref edac) scheme.
+
+The transport velocity ``\tilde{v}_a`` of particle ``a`` is used to evolve the position of the particle ``r_a`` from one time step to the next by
+
+```math
+\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a
+```
+
+and is obtained at every time-step ``\Delta t`` from
+
+```math
+\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right),
+```
+
+where ``\rho_a`` is the density of particle ``a`` and ``p_{\text{background}}`` 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
+
+```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.
+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 prefactor 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 can be chosen arbitrarily large when the time-step criterion is adjusted.
+
+The inviscid momentum equation with an additional convection term for a particle moving with ``\tilde{v}`` is
+
+```math
+\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p +  \nabla \cdot \bm{A},
+```
+
+ 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
+
+```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].
+```
+
+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``.
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "fluid", "transport_velocity.jl")]
+```

docs/src/tutorial.md

@@ -2,12 +2,12 @@
 
 
 ## Fluid
-- [Setting up your simulation from scratch](tutorials/tut_setup_replaced.md)
-- [Setting up a dam break simulation](tutorials/tut_dam_break_replaced.md)
+- [Setting up your simulation from scratch](tutorials/tut_setup.md)
+- [Setting up a dam break simulation](tutorials/tut_dam_break.md)
 
 ## Mechanics
-- [Deforming a beam](tutorials/tut_beam_replaced.md)
+- [Deforming a beam](tutorials/tut_beam.md)
 
 
 ## Fluid-Structure Interaction
-- [Setting up a falling structure](tutorials/tut_falling_replaced.md)
+- [Setting up a falling structure](tutorials/tut_falling.md)

examples/fluid/lid_driven_cavity_2d.jl

@@ -0,0 +1,104 @@
+# Lid-driven cavity
+#
+# S. Adami et al
+# "A transport-velocity formulation for smoothed particle hydrodynamics".
+# In: Journal of Computational Physics, Volume 241 (2013), pages 292-307.
+# https://doi.org/10.1016/j.jcp.2013.01.043
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.02
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 4
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 5.0)
+reynolds_number = 100.0
+
+cavity_size = (1.0, 1.0)
+
+fluid_density = 1.0
+
+const VELOCITY_LID = 1.0
+sound_speed = 10 * VELOCITY_LID
+
+pressure = sound_speed^2 * fluid_density
+
+viscosity = ViscosityAdami(; nu=VELOCITY_LID / reynolds_number)
+
+cavity = RectangularTank(particle_spacing, cavity_size, cavity_size, fluid_density,
+                         n_layers=boundary_layers,
+                         faces=(true, true, true, false), pressure=pressure)
+
+lid_position = 0.0 - particle_spacing * boundary_layers
+lid_length = cavity.n_particles_per_dimension[1] + 2boundary_layers
+
+lid = RectangularShape(particle_spacing, (lid_length, 3),
+                       (lid_position, cavity_size[2]), density=fluid_density)
+
+# ==========================================================================================
+# ==== Fluid
+
+smoothing_length = 1.0 * particle_spacing
+smoothing_kernel = SchoenbergQuinticSplineKernel{2}()
+
+fluid_system = EntropicallyDampedSPHSystem(cavity.fluid, smoothing_kernel, smoothing_length,
+                                           density_calculator=ContinuityDensity(),
+                                           sound_speed, viscosity=viscosity,
+                                           transport_velocity=TransportVelocityAdami(pressure))
+
+# ==========================================================================================
+# ==== Boundary
+
+lid_movement_function(t) = SVector(VELOCITY_LID * t, 0.0)
+
+is_moving(t) = true
+
+lid_movement = BoundaryMovement(lid_movement_function, is_moving)
+
+boundary_model_cavity = BoundaryModelDummyParticles(cavity.boundary.density,
+                                                    cavity.boundary.mass,
+                                                    AdamiPressureExtrapolation(),
+                                                    viscosity=viscosity,
+                                                    smoothing_kernel, smoothing_length)
+
+boundary_model_lid = BoundaryModelDummyParticles(lid.density, lid.mass,
+                                                 AdamiPressureExtrapolation(),
+                                                 viscosity=viscosity,
+                                                 smoothing_kernel, smoothing_length)
+
+boundary_system_cavity = BoundarySPHSystem(cavity.boundary, boundary_model_cavity)
+
+boundary_system_lid = BoundarySPHSystem(lid, boundary_model_lid, movement=lid_movement)
+
+# ==========================================================================================
+# ==== Simulation
+bnd_thickness = boundary_layers * particle_spacing
+periodic_box = PeriodicBox(min_corner=[-bnd_thickness, -bnd_thickness],
+                           max_corner=cavity_size .+ [bnd_thickness, bnd_thickness])
+
+semi = Semidiscretization(fluid_system, boundary_system_cavity, boundary_system_lid,
+                          neighborhood_search=GridNeighborhoodSearch{2}(; periodic_box))
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+
+saving_callback = SolutionSavingCallback(dt=0.02)
+
+pp_callback = nothing
+
+callbacks = CallbackSet(info_callback, saving_callback, pp_callback, UpdateCallback())
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(),
+            abstol=1e-6, # Default abstol is 1e-6 (may needs to be tuned to prevent boundary penetration)
+            reltol=1e-4, # Default reltol is 1e-3 (may needs to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            maxiters=Int(1e7),
+            save_everystep=false, callback=callbacks);