Merge branch 'add-nonuniform-BCs' into turbulence_modelling

Author: b-fg

src/Flow.jl

@@ -62,12 +62,10 @@ upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I
 """
     accelerate!(r,t,g,U)
 
-Adds a space-time acceleration field g(i,x,t) and velocity field U(i,x,t) such that `rᵢ += gᵢ(i,x,t)+dUᵢ(i,x,t)/dt` at time `t=sum(dt)` to vector field `r`.
+Accounts for applied and reference-frame acceleration using `rᵢ += g(i,x,t)+dU(i,x,t)/dt`
 """
 accelerate!(r,t,::Nothing,::Union{Nothing,Tuple}) = nothing
-accelerate!(r,t,f) = for i ∈ 1:last(size(r))
-    @loop r[I,i] += f(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r)))
-end
+accelerate!(r,t,f::Function) = @loop r[Ii] += f(last(Ii),loc(Ii,eltype(r)),t) over Ii ∈ CartesianIndices(r)
 accelerate!(r,t,g::Function,::Union{Nothing,Tuple}) = accelerate!(r,t,g)
 accelerate!(r,t,::Nothing,U::Function) = accelerate!(r,t,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t))
 accelerate!(r,t,g::Function,U::Function) = accelerate!(r,t,(i,x,t)->g(i,x,t)+ForwardDiff.derivative(τ->U(i,x,τ),t))

src/WaterLily.jl

@@ -67,7 +67,7 @@ mutable struct Simulation <: AbstractSimulation
                         Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(),
                         uλ=nothing, exitBC=false, body::AbstractBody=NoBody(),
                         T=Float32, mem=Array) where N
-        @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` and `uλ` cannot be both specified as Function"
+        @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both"
         @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function"
         isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)"
         isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N}),zero(T)),N)).==T) "`g` is not type stable"

src/util.jl

@@ -203,8 +203,8 @@ function BC!(a,u_BC::Function,saveexit=false,perdir=(),t=0)
                 end
                 (!saveexit || i>1) && (@loop a[I,i] = u_BC(i,loc(i,I),t) 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)
+                @loop a[I,i] = u_BC(i,loc(i,I),t)+a[I+δ(j,I),i]-u_BC(i,loc(i,I+δ(j,I)),t) over I ∈ slice(N,1,j)
+                @loop a[I,i] = u_BC(i,loc(i,I),t)+a[I-δ(j,I),i]-u_BC(i,loc(i,I-δ(j,I)),t) over I ∈ slice(N,N[j],j)
             end
         end
     end

test/maintests.jl

@@ -49,8 +49,8 @@ using ReadVTK, WriteVTK
         Ubc(i,x,t) = i==1 ? 1.0 : 0.5
         v = rand(Ng..., D) |> f # vector
         BC!(v,Ubc,false); BC!(u,U,false) # make sure we apply the same
-        @test all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1])
-        @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2])
+        @test GPUArrays.@allowscalar all(v[1, :, 1] .== u[1, :, 1]) && all(v[2, :, 1] .== u[2, :, 1]) && all(v[end, :, 1] .== u[end, :, 1])
+        @test GPUArrays.@allowscalar all(v[:, 1, 2] .== u[:, 1, 2]) && all(v[:, 2, 2] .== u[:, 2, 2]) && all(v[:, end, 2] .== u[:, end, 2])
         # test exit bc
         GPUArrays.@allowscalar v[end,:,1] .= 3
         BC!(v,Ubc,true) # save exit values
@@ -65,6 +65,20 @@ using ReadVTK, WriteVTK
         BC!(u,U,true,(1,)) #saveexit has no effect here as x-periodic
         @test GPUArrays.@allowscalar all(u[1:2, :, 1] .== u[end-1:end, :, 1]) && all(u[1:2, :, 2] .== u[end-1:end, :, 2]) &&
                            all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2])
+        # test non-uniform BCs
+        Ubc_1(i,x,t) = i==1 ? x[2] : x[1]
+        v .= 0; BC!(v,Ubc_1)
+        # the tangential BC change the value of the ghost cells on the other axis, so we cannot check it
+        @test GPUArrays.@allowscalar all(v[1,2:end-1,1] .≈ v[end,2:end-1,1])
+        @test GPUArrays.@allowscalar all(v[2:end-1,1,2] .≈ v[2:end-1,end,2])
+        # more complex
+        Ng, D = (8, 8, 8), 3
+        u = zeros(Ng..., D) |> f # vector
+        Ubc_2(i,x,t) = i==1 ? cos(2π*x[1]/8) : i==2 ? sin(2π*x[2]/8) : tan(π*x[3]/16)
+        BC!(u,Ubc_2)
+        @test GPUArrays.@allowscalar all(u[1,:,:,1] .≈ cos(-1π/4))  && all(u[2,:,:,1] .≈ cos(0)) && all(u[end,:,:,1] .≈ cos(6π/4))
+        @test GPUArrays.@allowscalar all(u[:,1,:,2] .≈ sin(-1π/4))  && all(u[:,2,:,2] .≈ sin(0)) && all(u[:,end,:,2] .≈ sin(6π/4))
+        @test GPUArrays.@allowscalar all(u[:,:,1,3] .≈ tan(-1π/16)) && all(u[:,:,2,3] .≈ tan(0)) && all(u[:,:,end,3].-tan(6π/16).<1e-6)
 
         # test interpolation
         a = zeros(5,5,2) |> f; b = zeros(5,5) |> f
@@ -334,6 +348,35 @@ end
         )
     end
 end
+@testset "Boundary Layer Flow" begin
+    for f ∈ arrays
+        make_bl_flow(L=32;T=Float32) = Simulation((L,L),
+            (i,x,t)-> i==1 ? convert(Float32,4.0*(((x[2]+0.5)/2L)-((x[2]+0.5)/2L)^2)) : 0.f0,
+            L;ν=0.001,U=1,mem=f,T,exitBC=false) # fails with exitBC=true, but the profile is maintained
+        sim = make_bl_flow(32)
+        sim_step!(sim,10)
+        @test GPUArrays.@allowscalar all(sim.flow.u[1,:,1] .≈ sim.flow.u[end,:,1])
+    end
+end
+
+@testset "Rotating reference frame" begin
+    function rotating_reference(N,x₀::SVector{2,T},ω::T,mem=Array) where T
+        function velocity(i,x,t)
+            s,c = sincos(ω*t); y = ω*(x-x₀)
+            i==1 ? s*y[1]+c*y[2] : -c*y[1]+s*y[2]
+        end
+        coriolis(i,x,t) = i==1 ? 2ω*velocity(2,x,t) : -2ω*velocity(1,x,t)
+        centrifugal(i,x,t) = ω^2*(x-x₀)[i]
+        g(i,x,t) = coriolis(i,x,t)+centrifugal(i,x,t)
+        udf(a::Flow,t) = WaterLily.@loop a.f[Ii] += g(last(Ii),loc(Ii,eltype(a.f)),t) over Ii in CartesianIndices(a.f)
+        simg = Simulation((N,N),velocity,N; g, U=1, T, mem) # use g
+        simg,Simulation((N,N),velocity,N; U=1, T, mem),udf
+    end
+    L = 4
+    simg,sim,udf = rotating_reference(2L,SA_F64[L,L],1/L,Array)
+    sim_step!(simg);sim_step!(sim;udf)
+    @test L₂(simg.flow.p)==L₂(sim.flow.p)<3e-3 # should be zero
+end
 
 @testset "Circle in accelerating flow" begin
     for f ∈ arrays