Add open boundaries

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
 

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/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/n_body/n_body_system.jl

@@ -5,12 +5,13 @@ struct NBodySystem{NDIMS, ELTYPE <: Real} <: TrixiParticles.System{NDIMS, Nothin
     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
 

src/TrixiParticles.jl

@@ -44,7 +44,8 @@ include("visualization/recipes_plots.jl")
 export Semidiscretization, semidiscretize, restart_with!
 export InitialCondition
 export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangianSPHSystem,
-       BoundarySPHSystem, DEMSystem, BoundaryDEMSystem
+       BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow,
+       OutFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
        PostprocessCallback, StepsizeCallback, UpdateCallback
 export ContinuityDensity, SummationDensity

src/callbacks/density_reinit.jl

@@ -74,7 +74,7 @@ function initialize_reinit_cb!(cb::DensityReinitializationCallback, u, t, integr
         # Update systems to compute quantities like density and pressure.
         semi = integrator.p
         v_ode, u_ode = u.x
-        update_systems_and_nhs(v_ode, u_ode, semi, t)
+        update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
 
         # Apply the callback.
         cb(integrator)

src/callbacks/post_process.jl

@@ -236,7 +236,7 @@ function (pp::PostprocessCallback)(integrator)
     new_data = false
 
     # Update systems to compute quantities like density and pressure
-    update_systems_and_nhs(v_ode, u_ode, semi, t)
+    update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
 
     foreach_system(semi) do system
         if system isa BoundarySystem && pp.exclude_boundary

src/callbacks/solution_saving.jl

@@ -139,7 +139,7 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in
         # Update systems to compute quantities like density and pressure
         semi = integrator.p
         v_ode, u_ode = u.x
-        update_systems_and_nhs(v_ode, u_ode, semi, t)
+        update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
 
         # Apply the callback
         solution_callback(integrator)

src/callbacks/update.jl

@@ -52,7 +52,16 @@ function initial_update!(cb, u, t, integrator)
     initial_update!(cb.affect!, u, t, integrator)
 end
 
-initial_update!(cb::UpdateCallback, u, t, integrator) = cb(integrator)
+function initial_update!(cb::UpdateCallback, u, t, integrator)
+    semi = integrator.p
+
+    # Tell systems that `UpdateCallback` is used
+    foreach_system(semi) do system
+        update_callback_used!(system)
+    end
+
+    return cb(integrator)
+end
 
 # `condition`
 function (update_callback!::UpdateCallback)(u, t, integrator)
@@ -69,16 +78,12 @@ function (update_callback!::UpdateCallback)(integrator)
 
     # Update quantities that are stored in the systems. These quantities (e.g. pressure)
     # still have the values from the last stage of the previous step if not updated here.
-    update_systems_and_nhs(v_ode, u_ode, semi, t)
+    update_systems_and_nhs(v_ode, u_ode, semi, t; update_from_callback=true)
 
     # Other updates might be added here later (e.g. Transport Velocity Formulation).
-    # @trixi_timeit timer() "update open boundary" foreach_system(semi) do system
-    #     update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t)
-    # end
-    #
-    # @trixi_timeit timer() "update TVF" foreach_system(semi) do system
-    #     update_transport_velocity_eachstep!(system, v_ode, u_ode, semi, t)
-    # end
+    @trixi_timeit timer() "update open boundary" foreach_system(semi) do system
+        update_open_boundary_eachstep!(system, v_ode, u_ode, semi, t)
+    end
 
     # Tell OrdinaryDiffEq that `u` has been modified
     u_modified!(integrator, true)
@@ -128,3 +133,5 @@ function Base.show(io::IO, ::MIME"text/plain",
         summary_box(io, "UpdateCallback", setup)
     end
 end
+
+update_callback_used!(system) = system