Convective exit

src/Flow.jl

@@ -70,17 +70,18 @@ struct Flow{D, T, Sf<:AbstractArray{T}, Vf<:AbstractArray{T}, Tf<:AbstractArray{
     U :: NTuple{D, T} # domain boundary values
     Δt:: Vector{T} # time step (stored in CPU memory)
     ν :: T # kinematic viscosity
-    function Flow(N::NTuple{D}, U::NTuple{D}; f=Array, Δt=0.25, ν=0., uλ::Function=(i, x) -> 0., T=Float64) where D
+    exitBC :: Bool # Convection exit
+    function Flow(N::NTuple{D}, U::NTuple{D}; f=Array, Δt=0.25, ν=0., uλ::Function=(i, x) -> 0., T=Float64, exitBC = false) where D
         Ng = N .+ 2
         Nd = (Ng..., D)
-        u = Array{T}(undef, Nd...) |> f; apply!(uλ, u); BC!(u, U)
+        u = Array{T}(undef, Nd...) |> f; apply!(uλ, u); BC!(u, U, exitBC); exitBC!(u,u,U,0.)
         u⁰ = copy(u)
         fv, p, σ = zeros(T, Nd) |> f, zeros(T, Ng) |> f, zeros(T, Ng) |> f
         V, σᵥ = zeros(T, Nd) |> f, zeros(T, Ng) |> f
         μ₀ = ones(T, Nd) |> f
         BC!(μ₀,ntuple(zero, D))
         μ₁ = zeros(T, Ng..., D, D) |> f
-        new{D,T,typeof(p),typeof(u),typeof(μ₁)}(u,u⁰,fv,p,σ,V,σᵥ,μ₀,μ₁,U,T[Δt],ν)
+        new{D,T,typeof(p),typeof(u),typeof(μ₁)}(u,u⁰,fv,p,σ,V,σᵥ,μ₀,μ₁,U,T[Δt],ν,exitBC)
     end
 end
 
@@ -91,7 +92,7 @@ function BDIM!(a::Flow)
 end
 
 function project!(a::Flow{n},b::AbstractPoisson,w=1) where n
-    dt = a.Δt[end]/w
+    dt = w*a.Δt[end]
     @inside b.z[I] = div(I,a.u)+a.σᵥ[I]; b.x .*= dt # set source term & solution IC
     solver!(b)
     for i ∈ 1:n  # apply solution and unscale to recover pressure
@@ -107,17 +108,19 @@ Integrate the `Flow` one time step using the [Boundary Data Immersion Method](ht
 and the `AbstractPoisson` pressure solver to project the velocity onto an incompressible flow.
 """
 @fastmath function mom_step!(a::Flow,b::AbstractPoisson)
-    a.u⁰ .= a.u; a.u .= 0
+    a.u⁰ .= a.u; scale_u!(a,0)
     # predictor u → u'
     conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν)
-    BDIM!(a); BC!(a.u,a.U)
-    project!(a,b); BC!(a.u,a.U)
+    BDIM!(a); BC!(a.u,a.U,a.exitBC)
+    a.exitBC && exitBC!(a.u,a.u⁰,a.U,a.Δt[end]) # convective exit
+    project!(a,b); BC!(a.u,a.U,a.exitBC)
     # corrector u → u¹
     conv_diff!(a.f,a.u,a.σ,ν=a.ν)
-    BDIM!(a); a.u ./= 2; BC!(a.u,a.U)
-    project!(a,b,2); BC!(a.u,a.U)
+    BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC)
+    project!(a,b,0.5); BC!(a.u,a.U,a.exitBC)
     push!(a.Δt,CFL(a))
 end
+scale_u!(a,scale) = @loop a.u[Ii] *= scale over Ii ∈ inside_u(size(a.p))
 
 function CFL(a::Flow)
     @inside a.σ[I] = flux_out(I,a.u)

src/WaterLily.jl

@@ -55,9 +55,9 @@ struct Simulation
     pois :: AbstractPoisson
     function Simulation(dims::NTuple{N}, u_BC::NTuple{N}, L::Number;
                         Δt=0.25, ν=0., U=√sum(abs2,u_BC), ϵ=1,
-                        uλ::Function=(i,x)->u_BC[i],
+                        uλ::Function=(i,x)->u_BC[i],exitBC=false,
                         body::AbstractBody=NoBody(),T=Float32,mem=Array) where N
-        flow = Flow(dims,u_BC;uλ,Δt,ν,T,f=mem)
+        flow = Flow(dims,u_BC;uλ,Δt,ν,T,f=mem,exitBC)
         measure!(flow,body;ϵ)
         new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ))
     end

src/util.jl

@@ -112,27 +112,23 @@
 
 using StaticArrays
 """
-    loc(i,I)
+    loc(i,I) = loc(Ii)
 
 Location in space of the cell at CartesianIndex `I` at face `i`.
 Using `i=0` returns the cell center s.t. `loc = I`.
 """
 @inline loc(i,I::CartesianIndex{N},T=Float64) where N = SVector{N,T}(I.I .- 0.5 .* δ(i,I).I)
-
+@inline loc(Ii::CartesianIndex,T=Float64) = loc(last(Ii),Base.front(Ii),T)
+Base.last(I::CartesianIndex) = last(I.I)
+Base.front(I::CartesianIndex) = CI(Base.front(I.I))
 """
     apply!(f, c)
 
 Apply a vector function `f(i,x)` to the faces of a uniform staggered array `c`.
 """
-function apply!(f,c)
-    N,n = size_u(c)
-    for i ∈ 1:n
-        @loop c[I,i] = f(i,loc(i,I)) over I ∈ CartesianIndices(N)
-    end
-end
-
+apply!(f,c) = @loop c[Ii] = f(last(Ii),loc(Ii)) over Ii ∈ CartesianIndices(c)
 """
-    slice(dims,i,j,low=1,trim=0)
+    slice(dims,i,j,low=1)
 
 Return `CartesianIndices` range slicing through an array of size `dims` in
 dimension `j` at index `i`. `low` optionally sets the lower extent of the range
@@ -150,20 +146,29 @@
 boundary. For example `aₓ(x)=Aₓ ∀ x ∈ minmax(X)`. A zero Neumann condition
 is applied to the tangential components.
 """
-function BC!(a,A)
+function BC!(a,A,saveexit=false)
     N,n = size_u(a)
     for j ∈ 1:n, i ∈ 1:n
         if i==j # Normal direction, Dirichlet
-            for s ∈ (1,2,N[j])
+            for s ∈ (1,2)
                 @loop a[I,i] = A[i] over I ∈ slice(N,s,j)
             end
+            (!saveexit || i>1) && (@loop a[I,i] = A[i] over I ∈ slice(N,N[j],j)) # overwrite exit
         else    # Tangential directions, Neumann
             @loop a[I,i] = a[I+δ(j,I),i] over I ∈ slice(N,1,j)
             @loop a[I,i] = a[I-δ(j,I),i] over I ∈ slice(N,N[j],j)
         end
     end
 end
 
+function exitBC!(u,u⁰,U,Δt)
+    N,_ = size_u(u)
+    exitR = slice(N.-1,N[1],1,2)              # exit slice excluding ghosts
+    @loop u[I,1] = u⁰[I,1]-U[1]*Δt*(u⁰[I,1]-u⁰[I-δ(1,I),1]) over I ∈ exitR
+    ∮u = sum(u[exitR,1])/length(exitR)-U[1]   # mass flux imbalance
+    @loop u[I,1] -= ∮u over I ∈ exitR         # correct flux
+end
+
 """
     BC!(a)
 Apply zero Neumann boundary conditions to the ghost cells of a _scalar_ field.

test/runtests.jl

@@ -32,7 +32,8 @@ arrays = setup_backends()
     @test ex == :(a[I, i] = Math.add(b[I], func(I, q)))
     @test sym == [:a, :I, :i, :(p.b), :q]
 
-    for f ∈ arrays
+    # for f ∈ arrays
+    for f ∈ [Array]
         p = Float64[i+j  for i ∈ 1:4, j ∈ 1:5] |> f
         @test inside(p) == CartesianIndices((2:3,2:4))
         @test L₂(p) == 187
@@ -46,14 +47,19 @@ arrays = setup_backends()
         σ = rand(Ng...) |> f # scalar
         BC!(u, U)
         BC!(σ)
-        @allowscalar() do
-            @test all(u[1, :, 1] .== U[1]) && all(u[2, :, 1] .== U[1]) && all(u[end, :, 1] .== U[1]) &&
+        @allowscalar @test all(u[1, :, 1] .== U[1]) && all(u[2, :, 1] .== U[1]) && all(u[end, :, 1] .== U[1]) &&
                 all(u[3:end-1, 1, 1] .== u[3:end-1, 2, 1]) && all(u[3:end-1, end, 1] .== u[3:end-1, end-1, 1])
-            @test all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) &&
+        @allowscalar @test all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) &&
                 all(u[1, 3:end-1, 2] .== u[2, 3:end-1, 2]) && all(u[end, 3:end-1, 2] .== u[end-1, 3:end-1, 2])
-            @test all(σ[1, 2:end-1] .== σ[2, 2:end-1]) && all(σ[end, 2:end-1] .== σ[end-1, 2:end-1]) &&
+        @allowscalar @test all(σ[1, 2:end-1] .== σ[2, 2:end-1]) && all(σ[end, 2:end-1] .== σ[end-1, 2:end-1]) &&
                 all(σ[2:end-1, 1] .== σ[2:end-1, 2]) && all(σ[2:end-1, end] .== σ[2:end-1, end-1])
-        end
+
+        @allowscalar u[end,:,1] .= 3
+        BC!(u, U, true) # save exit values
+        @allowscalar @test all(u[end, :, 1] .== 3)
+
+        WaterLily.exitBC!(u,u,U,0) # conservative exit check
+        @allowscalar @test all(u[end,2:end-1, 1] .== U[1])
     end
 end
 
@@ -205,13 +211,13 @@ end
     end
 end
 
-function sphere_sim(radius = 8; mem=Array)
-    body = AutoBody((x,t)-> √sum(abs2,x .- 2radius) - radius)
-    return Simulation(radius.*(6,4),(1,0),radius; body, ν=radius/250, T=Float32, mem)
+function sphere_sim(radius = 8; mem=Array, exitBC=false)
+    body = AutoBody((x,t)-> √sum(abs2,x .- (2radius+1.5)) - radius)
+    return Simulation(radius.*(6,4),(1,0),radius; body, ν=radius/250, T=Float32, mem, exitBC)
 end
 @testset "WaterLily.jl" begin
-    for mem ∈ arrays
-        sim = sphere_sim(32;mem);
+    for mem ∈ arrays, exitBC ∈ (true,false)
+        sim = sphere_sim(;mem,exitBC);
         @test sim_time(sim) == 0
         sim_step!(sim,0.1,remeasure=false)
         @test length(sim.flow.Δt)-1 == length(sim.pois.n)÷2