External BCs based on biot savart

Initial commit with some of the functions needed to apply a vorticity-based BC on the velocity using a "multipole" (really multigrid) method.

src/Flow.jl

@@ -130,3 +130,42 @@ end
     end
     return s
 end
+
+function ml_ω!(ml::MultiLevelPoisson,a::Flow)
+    # vorticity on finest grid
+    @inside ml.levels[1].z[I] = curl(3,I,a.u)
+    # pool values at each level
+    for l ∈ 2:length(ml.levels)
+        restrict!(ml.levels[l].z,ml.levels[l-1].z)
+    end
+end
+
+@inline @fastmath biotsavart(x,j,ω,I,dx) = (p = x-dx*(SA[Tuple(I)...] .-0.5); ω[I]*p[j]/(p'*p))
+
+function u_ω(i,x,ml)
+    # initialize at coarsest level
+    ui = zero(eltype(x)); j = i%2+1
+    l = 1 #length(ml.levels)
+    R = inside(ml.levels[l].z)
+    Imax,dx,ω = 0,0,0
+
+    # loop levels
+    while l>=1
+        # set grid scale and index nearest to x
+        ω = ml.levels[l].z
+        dx = 2^(l-1)
+        Imax = CI(round.(Int,x/dx .+0.5)...)
+
+        # get contributions other than Imax
+        for I in R
+            I != Imax && (ui += biotsavart(x,j,ω,I,dx))
+        end
+
+        # move "up" one level near Imax
+        l=l-1
+        R = up(Imax)
+    end
+
+    # add Imax contribution
+    return ui + biotsavart(x,j,ω,Imax,dx)
+end

src/Metrics.jl

@@ -4,10 +4,6 @@ using StaticArrays
 @inline fSV(f,n) = SA[ntuple(f,n)...]
 @inline @fastmath fsum(f,n) = sum(ntuple(f,n))
 norm2(x) = √(x'*x)
-@fastmath function permute(f,i)
-    j,k = i%3+1,(i+1)%3+1
-    f(j,k)-f(k,j)
-end
 ×(a,b) = fSV(i->permute((j,k)->a[j]*b[k],i),3)
 
 """

src/MultiLevelPoisson.jl

@@ -44,15 +44,16 @@ struct MultiLevelPoisson{T,S<:AbstractArray{T},V<:AbstractArray{T}} <: AbstractP
     L::V
     z::S
     levels :: Vector{Poisson{T,S,V}}
+    maxdepth :: Int
     n :: Vector{Int16}
-    function MultiLevelPoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T},maxlevels=4) where T
+    function MultiLevelPoisson(x::AbstractArray{T},L::AbstractArray{T},z::AbstractArray{T},maxdepth=4) where T
         levels = Poisson[Poisson(x,L,z)]
-        while all(size(levels[end].x) .|> divisible) && length(levels) <= maxlevels
+        while all(size(levels[end].x) .|> divisible)
             push!(levels,restrictML(levels[end]))
         end
         text = "MultiLevelPoisson requires size=a2ⁿ, where n>2"
         @assert (length(levels)>2) text
-        new{T,typeof(x),typeof(L)}(x,L,z,levels,[])
+        new{T,typeof(x),typeof(L)}(x,L,z,levels,min(maxdepth,length(levels)),[])
     end
 end
 function update!(ml::MultiLevelPoisson)
@@ -70,7 +71,7 @@ function Vcycle!(ml::MultiLevelPoisson;l=1)
     restrict!(coarse.r,fine.r)
     fill!(coarse.x,0.)
     # solve coarse (with recursion if possible)
-    l+1<length(ml.levels) && Vcycle!(ml,l=l+1)
+    l+1<ml.maxdepth && Vcycle!(ml,l=l+1)
     smooth!(coarse)
     # correct fine
     prolongate!(fine.ϵ,coarse.x)

src/util.jl

@@ -174,4 +174,10 @@ 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
+end
+
+@fastmath function permute(f,i)
+    j,k = i%3+1,(i+1)%3+1
+    f(j,k)-f(k,j)
+end
+curl(i,I,u) = permute((j,k)->∂(j,CI(I,k),u), i)