Working exit prototype

The main problem with the previous version were the `a.u .= 0` and `a.u ./= 2` lines in `mom_step!` which kept modifying the exit plane values. I changed these lines to only update the values inside the domain.

I added an `flow.exit` to the struct and initialization functions and use this to trigger the convection exit behavior.

The `apply!` function had the opposite problem - it WASN'T filling in the boundary value. Now it does.

It also seemed weird to me that the convection exit value in Lotus was based on the old value at the boundary but the u* value (before projection) upstream. I'm not sure it matters much, but it seems more consistent to just use the old values in both places, so that's what it does now.

Finally, I added a Lamp vortex dipole test case, which is nice since it should leave the exit without changing size or speed.

Convective-Exit.jl

@@ -1,106 +1,35 @@
-using WaterLily
-using BenchmarkTools
-include("examples/TwoD_plots.jl")
-
-function test_conv_BC!(sim)
-    # duration of the simulation
-    duration = 16
-    step = 0.4
-    t₀ = 0.0
-    plot()
-
-    @time @gif for tᵢ in range(t₀,t₀+duration;step)
-
-        # update until time tᵢ in the background
-        t = sum(sim.flow.Δt[1:end-1])
-        while t < tᵢ*sim.L/sim.U
-
-            # update flow
-            # mom_step_BC!(sim.flow,sim.pois)
-            WaterLily.mom_step!(sim.flow,sim.pois)
-
-            # compute motion and acceleration 1DOF
-            Δt = sim.flow.Δt[end]
-            t += Δt
-        end
-
-        # plot vorticity
-        # @inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
-        # flood(sim.flow.σ; shift=(-0.5,-0.5), clims=(-5,5))
-        # flood(sim.flow.p; shift=(-0.5,-0.5), clims=(-1,1))
-        # addbody(real(z.+n/4),imag(z.+im*m/2))
-        plot!(sim.flow.u[end,:,1])
-        # plot!(sim.flow.u[end,:,1])
-
-        # print time step
-        println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3))
+using WaterLily,SpecialFunctions,ForwardDiff,StaticArrays
+
+function lamb_dipole(N;D=N/3,U=1,exit=false)
+    β = 2.4394π/D
+    @show besselj1(β*D/2)
+    C = -2U/(β*besselj0(β*D/2))
+    function ψ(x,y)
+        r = √(x^2+y^2)
+        ifelse(r ≥ D/2, U*((D/2r)^2-1)*y, C*besselj1(β*r)*y/r)
     end
+    center = SA[N/2,N/2]
+    function uλ(i,xy)
+        x,y = xy-center
+        ifelse(i==1,ForwardDiff.derivative(y->ψ(x,y),y)+1+U,-ForwardDiff.derivative(x->ψ(x,y),x))
+    end
+    Simulation((N, N), (1,0), D; uλ, exit)
 end
 
-L = 4
-N = (L,L)
-Nd = N .+ 2
-Nd = (Nd...,2)
-U = (1.0,0.0)
-u = Array{Float64}(undef, Nd...)
-
-# vertical shear
-println("--Before - BC(u,U)--")
-apply!((i, x) -> (i==1 ? x[2] : 0), u)
-BC!(u,zeros(2))
-display(u[:,:,1]')
-
-println("--BCΔt!(u,U;Δt=1)--")
-BC!(u,zeros(2)) # reset BC values
-WaterLily.BCΔt!(u, (1,0); Δt=1.0) # convect and mass check
-display(u[:,:,1]')
-
-println("--BCΔt!(u,U)--") # same as BC!(u,1)
-BC!(u,zeros(2)) # reset BC values
-WaterLily.BCΔt!(u, (1,0)) # only mass check
-display(u[:,:,1]')
-
-# test on flow sim
-# some definitons
-U = 1
-Re = 250
-
-m, n = 32,48
-println("$n x $m: ", prod((m,n)))
-radius = 4
-
-# make a circle body
-body = AutoBody((x,t)->√sum(abs2, x .- [n/4,m/2]) - radius)
-z = [radius*exp(im*θ) for θ ∈ range(0,2π,length=33)]
-
-# make a simulation
-sim = Simulation((n,m), (U,0), radius; ν=U*radius/Re, body, T=Float64)
-test_conv_BC!(sim)
-
-# flood(sim.flow.u[2:end-1,2:end-1,1])
-# flood(sim.flow.p[2:end-1,2:end-1])
-# flood(sim.flow.μ₀[:,:,1])
-# plot(sim.flow.u[end,:,1])
-# plot!(sim.flow.u[end-1,:,1])
-
-
-# test all contributions to mom_step!
-# @benchmark WaterLily.mom_step!(sim.flow,sim.pois)
-
-# # step by step
-# println("conv_diff!")
-# @benchmark WaterLily.conv_diff!(sim.flow.f,sim.flow.u⁰,sim.flow.σ,ν=sim.flow.ν)
-
-# println("BDIM!")
-# @benchmark WaterLily.BDIM!(sim.flow)
-
-# println("BCΔt!")
-# @benchmark WaterLily.BCΔt!(sim.flow.u,sim.flow.U;Δt=sim.flow.Δt[end]) 
-
-# println("project!")
-# @benchmark WaterLily.project!(sim.flow,sim.pois)
-
-# println("BCΔt!")
-# @benchmark WaterLily.BC!(sim.flow.u,sim.flow.U)
-
+include("examples/TwoD_plots.jl")
 
+sim = lamb_dipole(64,U=0.25);
+@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
+flood(sim.flow.σ,clims=(-5,5))
+sim_step!(sim,1.2,remeasure=false)
+@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
+flood(sim.flow.σ,clims=(-5,5))
+
+sim = lamb_dipole(64,U=0.25,exit=true);
+sim_step!(sim,1.2,remeasure=false)
+@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
+flood(sim.flow.σ,clims=(-5,5))
+
+U=2
+sim = lamb_dipole(64,U=U,exit=true);
+sim_gif!(sim,duration=2.5/(1+U),step=0.05,clims=(-20,20))

Project.toml

@@ -11,6 +11,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

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
+    exit :: Bool # Convection exit
+    function Flow(N::NTuple{D}, U::NTuple{D}; f=Array, Δt=0.25, ν=0., uλ::Function=(i, x) -> 0., T=Float64, exit = 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, exit); 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],ν,exit)
     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Δt!(a.u,a.U;Δt=a.Δt[end]) # convective exit
-    project!(a,b); BCΔt!(a.u,a.U)
+    BDIM!(a); BC!(a.u,a.U,a.exit)
+    a.exit && exitBC!(a.u,a.u⁰,a.U,a.Δt[end]) # convective exit
+    project!(a,b); BC!(a.u,a.U,a.exit)
     # corrector u → u¹
     conv_diff!(a.f,a.u,a.σ,ν=a.ν)
-    BDIM!(a); a.u ./= 2; BCΔt!(a.u,a.U)
-    project!(a,b,2); BCΔt!(a.u,a.U)
+    BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exit)
+    project!(a,b,0.5); BC!(a.u,a.U,a.exit)
     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],exit=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,exit)
         measure!(flow,body;ϵ)
         new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ))
     end

src/util.jl

@@ -112,27 +112,23 @@ rep(ex::Expr) = ex.head == :. ? Symbol(ex.args[2].value) : ex
 
 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 @@ condition `a[I,i]=A[i]` is applied to the vector component _normal_ to the domai
 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)
+    exit = 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 ∈ exit
+    ∮u = sum(u[exit,1])/length(exit)-U[1]     # mass flux imbalance
+    @loop u[I,1] -= ∮u over I ∈ exit          # correct flux
+end
+
 """
     BC!(a)
 Apply zero Neumann boundary conditions to the ghost cells of a _scalar_ field.
@@ -174,57 +179,4 @@ function BC!(a)
         @loop a[I] = a[I+δ(j,I)] over I ∈ slice(N,1,j)
         @loop a[I] = a[I-δ(j,I)] over I ∈ slice(N,N[j],j)
     end
-end
-
-"""
-    BCΔt!(a,A;Δt=-1)
-Apply the boundary condition on a field `a` where time integration of the 1D convection
-equation is used for the exit in the i==j==1 direction. Time integration is performed 
-if Δt>0. Other Boundary conditions are applied as in `BC!(a,A)`. A global mass flux 
-check is performed in the 1-direction to ensure global mass conservation (not required 
-for the other direction).
-"""
-function BCΔt!(a,A;Δt=-1)
-    N,n = size_u(a)
-    area = ntuple(i->prod(N.-2)/(N[i]-2),n) # are of domain (without ghosts)
-    for j ∈ 1:n, i ∈ 1:n
-        if i==j # Normal direction, Dirichlet
-            for s ∈ (1,2)
-                @loop a[I,i] = A[i] over I ∈ slice(N,s,j)
-            end
-            if i==1 && Δt > 0 # convective exit
-                @loop a[I,i] = a[I,i]-Δt*(a[I,i]-a[I-δ(j,I),i])*A[i] over I ∈ slice(N,N[j],j)
-            end
-            # other direction, standard exit
-            i!=1 && @loop a[I,i] = A[i] over I ∈ slice(N,N[j],j)
-            # check flux
-            if i==1
-                ∮u = (∮(a,j)-A[j]*area[j])/∮(ones((N...,j)),j) # mass flux imbalance in domain
-                @loop a[I,i] -= ∮u over I ∈ slice(N,N[j],j)    # correct flux
-            end
-        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
-
-"""
-    Integrate the flux of `a` in the normal direction `j` over the exit
-    in that direction.
-"""
-function ∮(a::AbstractArray{T},j) where T
-    N,n = size_u(a)
-    return ∮(a,N,N[j],j)
-end
-"""
-    Integrate the flux of `a` in the normal direction `j` over a slice located
-    at `s``. This excluded ghosts cells.
-"""
-function ∮(a,N,s,j) 
-    sm = 0.0
-    for I ∈ slice(N.-1,s,j,2) # remove ghosts
-        sm += a[I,j] 
-    end
-    return sm
 end