Merge branch 'main' into edac-tvf

NEWS.md

@@ -1,8 +1,8 @@
 # Changelog
 
 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. 
-We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions. 
+used in the Julia ecosystem. Notable changes will be documented in this file for human readability.
+We aim at 3 to 4 month between major release versions and about 2 weeks between minor versions.
 
 ## Version 0.2.x
 
@@ -14,11 +14,15 @@ We aim at 3 to 4 month between major release versions and about 2 weeks between
 
 ### Deprecated
 
+## Version 0.1.3
+
+### Added
+Open boundaries using the method of characteristics based on the work of Lastiwka et al., "Permeable and non-reflecting boundary conditions in SPH" (2009) were added for WCSPH and EDAC.
 
 ## Version 0.1.2
 
 ### Added
-A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH
+A surface tension and adhesion model based on the work by Akinci et al., "Versatile Surface Tension and Adhesion for SPH Fluids" (2013) was added to WCSPH.
 
 ## Version 0.1.1
 

Project.toml

@@ -1,7 +1,7 @@
 name = "TrixiParticles"
 uuid = "66699cd8-9c01-4e9d-a059-b96c86d16b3a"
 authors = ["erik.faulhaber <[email protected]>"]
-version = "0.1.3-dev"
+version = "0.1.4-dev"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -11,7 +11,9 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastPow = "c0e83750-1142-43a8-81cf-6c956b72b4d1"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PointNeighbors = "1c4d5385-0a27-49de-8e2c-43b175c8985c"
@@ -36,14 +38,14 @@ FastPow = "0.1"
 ForwardDiff = "0.10"
 JSON = "0.21"
 MuladdMacro = "0.2"
+PointNeighbors = "0.2.3"
 Polyester = "0.7.5"
 RecipesBase = "1"
 Reexport = "1"
 SciMLBase = "1, 2"
 StaticArrays = "1"
 StrideArrays = "0.1"
 TimerOutputs = "0.5"
-TrixiBase = "0.1"
-PointNeighbors = "0.2"
+TrixiBase = "0.1.3"
 WriteVTK = "1"
 julia = "1.9"

README.md

@@ -17,7 +17,7 @@ Its features include:
 ## Features
 - Incompressible Navier-Stokes
   - Methods: Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH), Entropically Damped Artificial Compressibility (EDAC)
-  - Models: Surface Tension
+  - Models: Surface Tension, Open Boundaries
 - Solid-body mechanics
   - Methods:  Total Lagrangian SPH (TLSPH), Discrete Element Method (DEM)
 - Fluid-Structure Interaction

docs/src/systems/boundary.md

@@ -234,3 +234,86 @@ Pages = [joinpath("schemes", "boundary", "monaghan_kajtar", "monaghan_kajtar.jl"
 - Alireza Valizadeh, Joseph J. Monaghan. "A study of solid wall models for weakly compressible SPH."
   In: Journal of Computational Physics 300 (2015), pages 5–19.
   [doi: 10.1016/J.JCP.2015.07.033](https://doi.org/10.1016/J.JCP.2015.07.033)
+
+# [Open Boundaries](@id open_boundary)
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "boundary", "open_boundary", "system.jl")]
+```
+
+```@autodocs
+Modules = [TrixiParticles]
+Pages = [joinpath("schemes", "boundary", "open_boundary", "boundary_zones.jl")]
+```
+
+### [Method of characteristics](@id method_of_characteristics)
+
+The difficulty in non-reflecting boundary conditions, also called open boundaries, is to determine
+the appropriate boundary values of the exact characteristics of the Euler equations.
+Assuming the flow near the boundaries is normal to the boundary
+and free of shock waves and significant viscous effects, it can be shown that three characteristic variables exist:
+
+- ``J_1``, associated with convection of entropy and propagates at flow velocity,
+- ``J_2``, downstream-running characteristics,
+- ``J_3``, upstream-running characteristics.
+
+Giles (1990) derived those variables based on a linearized set of governing equations:
+```math
+J_1 = -c_s^2 (\rho - \rho_{\text{ref}}) + (p - p_{\text{ref}})
+```
+```math
+J_2 = \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}})
+```
+```math
+J_3 = - \rho c_s (v - v_{\text{ref}}) + (p - p_{\text{ref}})
+```
+where the subscript "ref" denotes the reference flow near the boundaries, which can be prescribed.
+
+Specifying the reference variables is **not** equivalent to prescription of ``\rho``, ``v`` and ``p``
+directly, since the perturbation from the reference flow is allowed.
+
+Lastiwka et al. (2009) applied the method of characteristic to SPH and determined the number of variables that should be
+**prescribed** at the boundary and the number which should be **propagated** from the fluid domain to the boundary:
+
+- For an **inflow** boundary:
+    - Prescribe *downstream*-running characteristics ``J_1`` and ``J_2``
+    - Transmit ``J_3`` from the fluid domain (allow ``J_3`` to propagate upstream to the boundary).
+
+- For an **outflow** boundary:
+    - Prescribe *upstream*-running characteristic ``J_3``
+    - Transmit ``J_1`` and ``J_2`` from the fluid domain.
+
+Prescribing is done by simply setting the characteristics to zero. To transmit the characteristics from the fluid
+domain, or in other words, to carry the information of the fluid to the boundaries, Negi et al. (2020) use a Shepard Interpolation
+```math
+f_i = \frac{\sum_j^N f_j W_{ij}}{\sum_j^N W_{ij}},
+```
+where the ``i``-th particle is a boundary particle, ``f`` is either  ``J_1``, ``J_2`` or ``J_3`` and ``N`` is the set of
+neighboring fluid particles.
+
+To express pressure ``p``, density ``\rho`` and velocity ``v`` as functions of the characteristic variables, the system of equations
+from the characteristic variables is inverted and gives
+```math
+ \rho - \rho_{\text{ref}} = \frac{1}{c_s^2} \left( -J_1 + \frac{1}{2} J_2 + \frac{1}{2} J_3 \right),
+```
+```math
+u - u_{\text{ref}}= \frac{1}{2\rho c_s} \left( J_2 - J_3 \right),
+```
+```math
+p - p_{\text{ref}} = \frac{1}{2} \left( J_2 + J_3 \right).
+```
+With ``J_1``, ``J_2`` and ``J_3`` determined, we can easily solve for the actual variables for each particle.
+
+### References
+- M. B. Giles. "Nonreflecting boundary conditions for Euler equation calculations".
+  In: AIAA Journal, 28.12 pages 2050--2058.
+  [doi: 10.2514/3.10521](https://doi.org/10.2514/3.10521)
+- M. Lastiwka, M. Basa, N. J. Quinlan.
+  "Permeable and non-reflecting boundary conditions in SPH".
+  In: International Journal for Numerical Methods in Fluids 61, (2009), pages 709--724.
+  [doi: 10.1002/fld.1971](https://doi.org/10.1002/fld.1971)
+- P. Negi, P. Ramachandran, A. Haftu.
+  "An improved non-reflecting outlet boundary condition for weakly-compressible SPH".
+  In: Computer Methods in Applied Mechanics and Engineering 367, (2020), pages 113--119.
+  [doi: 10.1016/j.cma.2020.113119](https://doi.org/10.1016/j.cma.2020.113119)

examples/fluid/dam_break_2d.jl

@@ -73,8 +73,10 @@ boundary_system = BoundarySPHSystem(tank.boundary, boundary_model, adhesion_coef
 
 # ==========================================================================================
 # ==== Simulation
-semi = Semidiscretization(fluid_system, boundary_system, threaded_nhs_update=true)
-ode = semidiscretize(semi, tspan)
+semi = Semidiscretization(fluid_system, boundary_system,
+                          neighborhood_search=GridNeighborhoodSearch,
+                          threaded_nhs_update=true)
+ode = semidiscretize(semi, tspan, data_type=nothing)
 
 info_callback = InfoCallback(interval=100)
 

examples/fluid/pipe_flow_2d.jl

@@ -0,0 +1,131 @@
+# 2D channel flow simulation with open boundaries.
+
+using TrixiParticles
+using OrdinaryDiffEq
+
+# ==========================================================================================
+# ==== Resolution
+particle_spacing = 0.05
+
+# Make sure that the kernel support of fluid particles at a boundary is always fully sampled
+boundary_layers = 3
+
+# Make sure that the kernel support of fluid particles at an open boundary is always
+# fully sampled.
+# Note: Due to the dynamics at the inlets and outlets of open boundaries,
+# it is recommended to use `open_boundary_layers > boundary_layers`
+open_boundary_layers = 6
+
+# ==========================================================================================
+# ==== Experiment Setup
+tspan = (0.0, 2.0)
+
+# Boundary geometry and initial fluid particle positions
+domain_size = (1.0, 0.4)
+
+flow_direction = [1.0, 0.0]
+reynolds_number = 100
+const prescribed_velocity = 2.0
+
+boundary_size = (domain_size[1] + 2 * particle_spacing * open_boundary_layers,
+                 domain_size[2])
+
+fluid_density = 1000.0
+
+# For this particular example, it is necessary to have a background pressure.
+# Otherwise the suction at the outflow is to big and the simulation becomes unstable.
+pressure = 1000.0
+
+sound_speed = 10 * prescribed_velocity
+
+state_equation = StateEquationCole(; sound_speed, reference_density=fluid_density,
+                                   exponent=7, background_pressure=pressure)
+
+pipe = RectangularTank(particle_spacing, domain_size, boundary_size, fluid_density,
+                       pressure=pressure, n_layers=boundary_layers,
+                       faces=(false, false, true, true))
+
+# Shift pipe walls in negative x-direction for the inflow
+pipe.boundary.coordinates[1, :] .-= particle_spacing * open_boundary_layers
+
+n_buffer_particles = 4 * pipe.n_particles_per_dimension[2]
+
+# ==========================================================================================
+# ==== Fluid
+smoothing_length = 3.0 * particle_spacing
+smoothing_kernel = WendlandC2Kernel{2}()
+
+fluid_density_calculator = ContinuityDensity()
+
+kinematic_viscosity = prescribed_velocity * domain_size[2] / reynolds_number
+
+viscosity = ViscosityAdami(nu=kinematic_viscosity)
+
+fluid_system = EntropicallyDampedSPHSystem(pipe.fluid, smoothing_kernel, smoothing_length,
+                                           sound_speed, viscosity=viscosity,
+                                           density_calculator=fluid_density_calculator,
+                                           buffer_size=n_buffer_particles)
+
+# Alternatively the WCSPH scheme can be used
+# alpha = 8 * kinematic_viscosity / (smoothing_length * sound_speed)
+# viscosity = ArtificialViscosityMonaghan(; alpha, beta=0.0)
+
+# fluid_system = WeaklyCompressibleSPHSystem(pipe.fluid, fluid_density_calculator,
+#                                            state_equation, smoothing_kernel,
+#                                            smoothing_length, viscosity=viscosity,
+#                                            buffer_size=n_buffer_particles)
+
+# ==========================================================================================
+# ==== Open Boundary
+function velocity_function(pos, t)
+    # Use this for a time-dependent inflow velocity
+    # return SVector(0.5prescribed_velocity * sin(2pi * t) + prescribed_velocity, 0)
+
+    return SVector(prescribed_velocity, 0.0)
+end
+
+inflow = InFlow(; plane=([0.0, 0.0], [0.0, domain_size[2]]), flow_direction,
+                open_boundary_layers, density=fluid_density, particle_spacing)
+
+open_boundary_in = OpenBoundarySPHSystem(inflow; sound_speed, fluid_system,
+                                         buffer_size=n_buffer_particles,
+                                         reference_pressure=pressure,
+                                         reference_velocity=velocity_function)
+
+outflow = OutFlow(; plane=([domain_size[1], 0.0], [domain_size[1], domain_size[2]]),
+                  flow_direction, open_boundary_layers, density=fluid_density,
+                  particle_spacing)
+
+open_boundary_out = OpenBoundarySPHSystem(outflow; sound_speed, fluid_system,
+                                          buffer_size=n_buffer_particles,
+                                          reference_pressure=pressure,
+                                          reference_velocity=velocity_function)
+
+# ==========================================================================================
+# ==== Boundary
+
+boundary_model = BoundaryModelDummyParticles(pipe.boundary.density, pipe.boundary.mass,
+                                             AdamiPressureExtrapolation(),
+                                             state_equation=state_equation,
+                                             #viscosity=ViscosityAdami(nu=1e-4),
+                                             smoothing_kernel, smoothing_length)
+
+boundary_system = BoundarySPHSystem(pipe.boundary, boundary_model)
+
+# ==========================================================================================
+# ==== Simulation
+semi = Semidiscretization(fluid_system, open_boundary_in, open_boundary_out,
+                          boundary_system)
+
+ode = semidiscretize(semi, tspan)
+
+info_callback = InfoCallback(interval=100)
+saving_callback = SolutionSavingCallback(dt=0.02, prefix="")
+
+callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
+
+sol = solve(ode, RDPK3SpFSAL35(),
+            abstol=1e-5, # Default abstol is 1e-6 (may need to be tuned to prevent boundary penetration)
+            reltol=1e-3, # Default reltol is 1e-3 (may need to be tuned to prevent boundary penetration)
+            dtmax=1e-2, # Limit stepsize to prevent crashing
+            save_everystep=false, callback=callbacks);

examples/fluid/taylor_green_vortex_2d.jl

@@ -139,8 +139,7 @@ info_callback = InfoCallback(interval=100)
 saving_callback = SolutionSavingCallback(dt=0.02,
                                          L1v=compute_L1v_error,
                                          L1p=compute_L1p_error,
-                                         diff_p_loc_p_avg=diff_p_loc_p_avg,
-                                         t=time_vector)
+                                         diff_p_loc_p_avg=diff_p_loc_p_avg)
 
 callbacks = CallbackSet(info_callback, saving_callback, UpdateCallback())
 

examples/n_body/n_body_benchmark_trixi.jl

@@ -108,7 +108,7 @@ end
 filename = tempname()
 open(filename, "w") do f
     redirect_stderr(f) do
-        TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+        TrixiParticles.disable_debug_timings()
     end
 end
 
@@ -124,6 +124,6 @@ end
 # Enable timers again
 open(filename, "w") do f
     redirect_stderr(f) do
-        TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+        TrixiParticles.enable_debug_timings()
     end
 end

examples/n_body/n_body_solar_system.jl

@@ -49,7 +49,7 @@ callbacks = CallbackSet(info_callback, saving_callback)
 filename = tempname()
 open(filename, "w") do f
     redirect_stderr(f) do
-        TrixiParticles.TimerOutputs.disable_debug_timings(TrixiParticles)
+        TrixiParticles.disable_debug_timings()
     end
 end
 
@@ -63,6 +63,6 @@ sol = solve(ode, SymplecticEuler(),
 # Enable timers again
 open(filename, "w") do f
     redirect_stderr(f) do
-        TrixiParticles.TimerOutputs.enable_debug_timings(TrixiParticles)
+        TrixiParticles.enable_debug_timings()
     end
 end

examples/n_body/n_body_system.jl

@@ -1,16 +1,19 @@
 using TrixiParticles
 using LinearAlgebra
 
-struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothing}
+# The second type parameter of `System` can't be `Nothing`, or TrixiParticles will launch
+# GPU kernel for `for_particle_neighbor` loops.
+struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, 0}
     initial_condition :: InitialCondition{ELTYPE}
     mass              :: Array{ELTYPE, 1} # [particle]
     G                 :: ELTYPE
+    buffer            :: Nothing
 
     function NBodySystem(initial_condition, G)
         mass = copy(initial_condition.mass)
 
         new{size(initial_condition.coordinates, 1),
-            eltype(mass)}(initial_condition, mass, G)
+            eltype(mass)}(initial_condition, mass, G, nothing)
     end
 end
 

examples/postprocessing/interpolation_plane.jl

@@ -96,3 +96,25 @@ plot_3d = Plots.plot(scatter_3d, xlabel="X", ylabel="Y", zlabel="Z",
 
 combined_plot = Plots.plot(plot1, plot2, plot3, plot_3d, layout=(2, 2),
                            size=(1000, 1500), margin=3mm)
+
+# If we want to save planes at regular intervals, we can use the postprocessing callback.
+function save_interpolated_plane(v, u, t, system)
+    # Size of the patch to be interpolated
+    interpolation_start = [0.0, 0.0]
+    interpolation_end = [tank_size[1], tank_size[2]]
+
+    # The resolution the plane is interpolated to. In this case twice the original resolution.
+    resolution = 0.5 * fluid_particle_spacing
+
+    file_id = ceil(Int, t * 10000)
+    interpolate_plane_2d_vtk(interpolation_start, interpolation_end, resolution,
+                             semi, system, v, u, filename="plane_$file_id.vti")
+    return nothing
+end
+
+save_interpolation_cb = PostprocessCallback(; dt=0.1, write_file_interval=0,
+                                            save_interpolated_plane)
+
+trixi_include(@__MODULE__,
+              joinpath(examples_dir(), "fluid", "dam_break_2d.jl"), tspan=(0.0, 0.2),
+              extra_callback=save_interpolation_cb, fluid_particle_spacing=0.01)