trixi-framework · efaulhaber · Jun 19, 2024 · Feb 5, 2024 · Feb 7, 2024 · Feb 7, 2024
diff --git a/NEWS.md b/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) was 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
 

diff --git a/README.md b/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

diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
@@ -293,6 +293,9 @@ function semidiscretize(semi, tspan; reset_threads=true, data_type=nothing)
 
             # Initialize this system
             initialize!(system, neighborhood_search)
+
+            # Only for systems requiring a mandatory callback
+            initial_callback_flag!(system)
         end
     end
 
@@ -355,6 +358,9 @@ function restart_with!(semi, sol; reset_threads=true)
         u = wrap_u(sol.u[end].x[2], system, semi)
 
         restart_with!(system, v, u)
+
+        # Only for systems requiring a mandatory callback
+        initial_callback_flag!(system)
     end
 
     return semi

diff --git a/src/general/system.jl b/src/general/system.jl
@@ -119,3 +119,6 @@ end
 function update_final!(system, v, u, v_ode, u_ode, semi, t; update_from_callback=false)
     return system
 end
+
+# Only for systems requiring a mandatory callback
+initial_callback_flag!(system) = system
diff --git a/src/schemes/boundary/open_boundary/boundary_zones.jl b/src/schemes/boundary/open_boundary/boundary_zones.jl
@@ -9,17 +9,17 @@ The specified plane (line in 2D or rectangle in 3D) will be extruded in upstream
 direction (the direction opposite to `flow_direction`) to create a box for the boundary zone.
 There are three ways to specify the actual shape of the inflow:
 1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then
-   be filled with inflow particles.
-2. Pass a 1D shape in 2D or a 2D shape in 3D, which is then extruded in upstream direction
-   to create the inflow particles.
+   be filled with inflow particles (default).
+2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D,
+   which is then extruded in upstream direction to create the inflow particles.
    - In 2D, the shape must be either an initial condition with 2D coordinates, which lies
      on the line specified by `plane`, or an initial condition with 1D coordinates, which lies
      on the line specified by `plane` when a y-coordinate of `0` is added.
    - In 3D, the shape must be either an initial condition with 3D coordinates, which lies
      in the rectangle specified by `plane`, or an initial condition with 2D coordinates,
      which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added.
-3. Pass a 2D initial condition in 2D or a 3D initial condition in 3D, which will be used
-   for the inflow particles.
+3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D,
+   which will be used for the inflow particles.
 
 !!! note "Note"
     Particles outside the boundary zone box will be removed.
@@ -145,17 +145,17 @@ The specified plane (line in 2D or rectangle in 3D) will be extruded in downstre
 direction (the direction in `flow_direction`) to create a box for the boundary zone.
 There are three ways to specify the actual shape of the outflow:
 1. Don't pass `initial_condition` or `extrude_geometry`. The boundary zone box will then
-   be filled with outflow particles.
-2. Pass a 1D shape in 2D or a 2D shape in 3D, which is then extruded in downstream direction
-   to create the outflow particles.
+   be filled with outflow particles (default).
+2. Specify `extrude_geometry` by passing a 1D shape in 2D or a 2D shape in 3D,
+   which is then extruded in downstream direction to create the outflow particles.
     - In 2D, the shape must be either an initial condition with 2D coordinates, which lies
       on the line specified by `plane`, or an initial condition with 1D coordinates, which lies
       on the line specified by `plane` when a y-coordinate of `0` is added.
     -  In 3D, the shape must be either an initial condition with 3D coordinates, which lies
       in the rectangle specified by `plane`, or an initial condition with 2D coordinates,
       which lies in the rectangle specified by `plane` when a z-coordinate of `0` is added.
-3. Pass a 2D initial condition in 2D or a 3D initial condition in 3D, which will be used
-   for the outflow particles.
+3. Specify `initial_condition` by passing a 2D initial condition in 2D or a 3D initial condition in 3D,
+   which will be used for the outflow particles.
 
 !!! note "Note"
     Particles outside the boundary zone box will be removed.

diff --git a/src/schemes/boundary/open_boundary/system.jl b/src/schemes/boundary/open_boundary/system.jl
@@ -146,6 +146,12 @@ function Base.show(io::IO, ::MIME"text/plain", system::OpenBoundarySPHSystem)
     end
 end
 
+function initial_callback_flag!(system::OpenBoundarySPHSystem)
+    system.update_callback_used[] = false
+
+    return system
+end
+
 update_callback_used!(system::OpenBoundarySPHSystem) = system.update_callback_used[] = true
 
 @inline source_terms(system::OpenBoundarySPHSystem) = nothing

diff --git a/test/Project.toml b/test/Project.toml
@@ -6,7 +6,6 @@ GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"

diff --git a/test/schemes/boundary/open_boundary/boundary_zone.jl b/test/schemes/boundary/open_boundary/boundary_zone.jl
@@ -6,7 +6,7 @@
         plane_points_1 = [[0.0, 0.0], [0.5, -0.5], [1.0, 0.5]]
         plane_points_2 = [[0.0, 1.0], [0.2, 2.0], [2.3, 0.5]]
 
-        @testset "Points $(i)" for i in eachindex(plane_points_1)
+        @testset verbose=true "Points $(i)" for i in eachindex(plane_points_1)
             point_1 = plane_points_1[i]
             point_2 = plane_points_2[i]
 
@@ -17,7 +17,7 @@
                 -normalize([-plane_size[2], plane_size[1]]),
             ]
 
-            @testset "Flow Direction $(j)" for j in eachindex(flow_directions)
+            @testset verbose=true "Flow Direction $(j)" for j in eachindex(flow_directions)
                 inflow = InFlow(; plane=(point_1, point_2), particle_spacing,
                                 flow_direction=flow_directions[j], density=1.0,
                                 open_boundary_layers)
@@ -30,17 +30,18 @@
                     outflow,
                 ]
 
-                @testset "$boundary_zone" for boundary_zone in boundary_zones
+                @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (boundary_zone isa InFlow) ? -1 : 1
 
                     @test plane_points_1[i] == boundary_zone.zone_origin
                     @test plane_points_2[i] - boundary_zone.zone_origin ==
                           boundary_zone.spanning_set[2]
-                    @test sign_ * flow_directions[j] ≈
-                          normalize(boundary_zone.spanning_set[1])
-                    @test zone_width ≈ norm(boundary_zone.spanning_set[1])
+                    @test isapprox(sign_ * flow_directions[j],
+                                   normalize(boundary_zone.spanning_set[1]), atol=1e-14)
+                    @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]),
+                                   atol=1e-14)
                 end
             end
         end
@@ -63,7 +64,7 @@
             [0.3113730847835541, 0.25057315826416016, -0.02374829351902008],
         ]
 
-        @testset "Points $(i)" for i in eachindex(plane_points_1)
+        @testset verbose=true "Points $(i)" for i in eachindex(plane_points_1)
             point_1 = plane_points_1[i]
             point_2 = plane_points_2[i]
             point_3 = plane_points_3[i]
@@ -76,7 +77,7 @@
                 -normalize(cross(edge1, edge2)),
             ]
 
-            @testset "Flow Direction $(j)" for j in eachindex(flow_directions)
+            @testset verbose=true "Flow Direction $(j)" for j in eachindex(flow_directions)
                 inflow = InFlow(; plane=(point_1, point_2, point_3), particle_spacing,
                                 flow_direction=flow_directions[j], density=1.0,
                                 open_boundary_layers)
@@ -89,7 +90,7 @@
                     outflow,
                 ]
 
-                @testset "$boundary_zone" for boundary_zone in boundary_zones
+                @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
                     zone_width = open_boundary_layers *
                                  boundary_zone.initial_condition.particle_spacing
                     sign_ = (boundary_zone isa InFlow) ? -1 : 1
@@ -99,9 +100,10 @@
                           boundary_zone.spanning_set[2]
                     @test plane_points_3[i] - boundary_zone.zone_origin ==
                           boundary_zone.spanning_set[3]
-                    @test sign_ * flow_directions[j] ≈
-                          normalize(boundary_zone.spanning_set[1])
-                    @test zone_width ≈ norm(boundary_zone.spanning_set[1])
+                    @test isapprox(sign_ * flow_directions[j],
+                                   normalize(boundary_zone.spanning_set[1]), atol=1e-14)
+                    @test isapprox(zone_width, norm(boundary_zone.spanning_set[1]),
+                                   atol=1e-14)
                 end
             end
         end
@@ -123,26 +125,25 @@
             outflow,
         ]
 
-        @testset "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
-            perturb_ = boundary_zone isa InFlow ? eps() : -eps()
+        @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
+            perturb_ = boundary_zone isa InFlow ? sqrt(eps()) : -sqrt(eps())
 
-            point_3 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin
+            point1 = plane_points[1]
+            point2 = plane_points[2]
+            point3 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin
 
             query_points = Dict(
                 "Behind" => ([-1.0, -1.0], false),
                 "Before" => ([2.0, 2.0], false),
-                "On Point 1" => (plane_points[1], true),
-                "On Point 2" => (plane_points[2], true),
-                "On Point 3" => (point_3, true),
-                "Closely On Point 1" => (plane_points[1] .+ perturb_ * flow_direction,
-                                         false))
+                "Closely On Point 1" => (point1 + perturb_ * flow_direction, false),
+                "Closely On Point 2" => (point2 + perturb_ * flow_direction, false),
+                "Closely On Point 3" => (point3 - perturb_ * flow_direction, false))
 
-            @testset "$k" for k in keys(query_points)
+            @testset verbose=true "$k" for k in keys(query_points)
                 (particle_position, evaluation) = query_points[k]
 
                 @test evaluation ==
-                      TrixiParticles.is_in_boundary_zone(boundary_zone,
-                                                         particle_position .± eps())
+                      TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position)
             end
         end
     end
@@ -164,24 +165,23 @@
             outflow,
         ]
 
-        @testset "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
+        @testset verbose=true "$(nameof(typeof(boundary_zone)))" for boundary_zone in boundary_zones
             perturb_ = boundary_zone isa InFlow ? eps() : -eps()
+            point4 = boundary_zone.spanning_set[1] + boundary_zone.zone_origin
+
             query_points = Dict(
                 "Behind" => ([-1.0, -1.0, 1.2], false),
                 "Before" => ([2.0, 2.0, -1.2], false),
-                "On Point 1" => (point1, true),
-                "On Point 2" => (point2, true),
-                "On Point 3" => (point3, true),
-                "On Point 4" => (boundary_zone.spanning_set[1] + boundary_zone.zone_origin,
-                                 true),
-                "Closely On Point 1" => (point1 .+ perturb_ * flow_direction, false))
-
-            @testset "$k" for k in keys(query_points)
+                "Closely On Point 1" => (point1 + perturb_ * flow_direction, false),
+                "Closely On Point 2" => (point2 + perturb_ * flow_direction, false),
+                "Closely On Point 3" => (point3 + perturb_ * flow_direction, false),
+                "Closely On Point 4" => (point4 - perturb_ * flow_direction, false))
+
+            @testset verbose=true "$k" for k in keys(query_points)
                 (particle_position, evaluation) = query_points[k]
 
                 @test evaluation ==
-                      TrixiParticles.is_in_boundary_zone(boundary_zone,
-                                                         particle_position .± eps())
+                      TrixiParticles.is_in_boundary_zone(boundary_zone, particle_position)
             end
         end
     end

diff --git a/test/test_util.jl b/test/test_util.jl
@@ -7,7 +7,6 @@ using DataFrames: DataFrame
 using JSON: JSON
 using QuadGK: quadgk # For integration in smoothing kernel tests
 using Random: Random # For rectangular patch
-using Measurements: ±
 using Polyester: disable_polyester_threads # For `count_rhs_allocations`
 
 """