From 633557e8cd97f1cba4da17c54519311f9c72625a Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 7 Sep 2023 10:16:55 +0200 Subject: [PATCH 1/6] initial Biot-Savart BC --- BiotSavart.jl | 73 ++++++++++++++++++ Biot_Savart_BC.jl | 165 +++++++++++++++++++++++++++++++++++++++++ Biot_Savart_Profile.jl | 96 ++++++++++++++++++++++++ Biot_Savart_Vortex.jl | 65 ++++++++++++++++ 4 files changed, 399 insertions(+) create mode 100644 BiotSavart.jl create mode 100644 Biot_Savart_BC.jl create mode 100644 Biot_Savart_Profile.jl create mode 100644 Biot_Savart_Vortex.jl diff --git a/BiotSavart.jl b/BiotSavart.jl new file mode 100644 index 00000000..55122654 --- /dev/null +++ b/BiotSavart.jl @@ -0,0 +1,73 @@ +using LinearAlgebra +""" +Biot-Savart integral +""" +function BiotSavart!(u,ω) + N,n = WaterLily.size_u(u) + u .= 0.0 + for i ∈ 1:n + j=i%2+1 # the only component not zero in the vorticity + for I ∈ WaterLily.inside_u(N,i) + @WaterLily.loop u[I,i] += K(i,I,j,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) + end + end +end + +""" +Biot-Savart integral +""" +function BiotSavart!(Is,u,ω) + N,n = WaterLily.size_u(u) + for i ∈ 1:n + j=i%2+1 # the only component not zero in the vorticity + for I ∈ Is + u[I,i] = 0.0 + @WaterLily.loop u[I,i] += K(i,I,j,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) + end + end +end + +""" +Biot-Savart kernel +""" +function K(i,I,j,J,ϵ=1e-6) + # face centered velocity at I due to vorticity at cell edge J + r = loc(i,I) .- loc(0,J) .+ 0.5 # cell edge vorticity + rⁿ = norm(r)^2 + return sign(i-j)*r[j]/(2π*rⁿ+ϵ^2) +end + +function error_at_slice(ue,u,dn,up) + N,n = WaterLily.size_u(u) + error = 0.0 + for i ∈ 1:n, s ∈ [up,dn] + for I ∈ WaterLily.slice(N,s,2) + error += abs.(u[I,i] - ue[I,i]) + end + end + return error / (2N[1]) +end + +function velocity_at_slice(ue,dn,up) + res = zeros(eltype(ue),size(u)) + N,n = WaterLily.size_u(u) + for i ∈ 1:n, s ∈ [up,dn] + for I ∈ WaterLily.slice(N,s,2) + res[I,i] = ue[I,i] + end + end + return res +end + +function ∮(a::AbstractArray{T},j) where T + N,n = WaterLily.size_u(a) + return ∮(a,N,N[j],j) +end + +function ∮(a,N,s,j) + sm = 0.0 + for I ∈ WaterLily.slice(N.-1,s,j,2) # remove ghosts + sm += a[I,j] + end + return sm +end diff --git a/Biot_Savart_BC.jl b/Biot_Savart_BC.jl new file mode 100644 index 00000000..85df0322 --- /dev/null +++ b/Biot_Savart_BC.jl @@ -0,0 +1,165 @@ +using WaterLily +using LinearAlgebra +using Plots; gr() +# include("examples/TwoD_plots.jl") +include("BiotSavart.jl") + +# function ∮(a::AbstractArray{T},j) where T +# N,n = WaterLily.size_u(a) +# return ∮(a,N,N[j],j) +# end + +# function ∮(a,N,s,j) +# sm = 0.0 +# for I ∈ WaterLily.slice(N.-1,s,j,2) # remove ghosts +# sm += a[I,j] +# end +# return sm +# end + +function BC_BiotSavart!(a,A,ω) + N,n = WaterLily.size_u(a) + @inside ω[I] = WaterLily.curl(3,I,a) + area = ntuple(i->prod(N.-2)/(N[i]-2),n) # are of domain (without ghosts) + for j ∈ 1:n, i ∈ 1:n + k = i%2+1 # the only component not zero in the vorticity + if i==j # Normal direction, Dirichlet + if i==1 + for s ∈ (1,2,N[j]) + for I ∈ WaterLily.slice(N,s,j) + a[I,i] = A[i] + @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) + end + end + # @WaterLily.loop a[I,i] = A[i] over I ∈ WaterLily.slice(N,N[j],j) + else + for s ∈ (1,2,N[j]) + for I ∈ WaterLily.slice(N,s,j) + a[I,i] = A[i] + @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) + end + end + end + # make global face flux zero + # ∮u = (∮(a,N,N[j],j)-∮(a,N,2,j))/∮(ones((N...,j)),j)# mass flux imbalance in domain + # @WaterLily.loop a[I,i] += ∮u/2 over I ∈ WaterLily.slice(N,2,j) + # @WaterLily.loop a[I,i] -= ∮u/2 over I ∈ WaterLily.slice(N,N[j],j) + else # Tangential directions, Neumanns + # @WaterLily.loop a[I,i] = a[I+δ(j,I),i] over I ∈ WaterLily.slice(N,1,j) + # @WaterLily.loop a[I,i] = a[I-δ(j,I),i] over I ∈ WaterLily.slice(N,N[j],j) + for s ∈ (1,N[j]) + for I ∈ WaterLily.slice(N,s,j) + a[I,i] = A[i] + @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) + end + end + end + end +end + +BC!(a,A,σ) = BC_BiotSavart!(a,A,σ) +# BC!(a,A,σ) = WaterLily.BC!(a,A) + +@fastmath function mom_step_inner!(a::Flow,b::AbstractPoisson) + a.u⁰ .= a.u; a.u .= 0 + # predictor u → u' + WaterLily.conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν) + WaterLily.BDIM!(a); + @inside a.σ[I] = WaterLily.curl(3,I,a.u) + BC!(a.u,a.U,a.σ) + WaterLily.project!(a,b) + @inside a.σ[I] = WaterLily.curl(3,I,a.u) + BC!(a.u,a.U,a.σ) + + # corrector u → u¹ + WaterLily.conv_diff!(a.f,a.u,a.σ,ν=a.ν) + WaterLily.BDIM!(a); a.u ./= 2; + @inside a.σ[I] = WaterLily.curl(3,I,a.u) + BC!(a.u,a.U,a.σ) + b.x ./=2 #scaled pressure guess, no necessary + WaterLily.project!(a,b); + @inside a.σ[I] = WaterLily.curl(3,I,a.u) + BC!(a.u,a.U,a.σ) + push!(a.Δt,WaterLily.CFL(a)) +end + + +# some definitons +U = 1 +Re = 250 +m, n = 2^6, 2^7 +println("$m x $n: ", prod((m,n))) +radius = 16 +center = 32 + 1.5 + +# make a circle body +body = AutoBody((x,t)->√sum(abs2, x .- center) - 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=Float32) + +# duration of the simulation +duration = 0.0 +step = 0.1 +t₀ = 0.0 + +# init plot +N,dim = WaterLily.size_u(sim.flow.u) +Px = 1; Py1 = 2; Py2=N[2] +l = @layout [ a ; b c ] +@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) +p1 = contourf(axes(sim.flow.σ,1),axes(sim.flow.σ,2),clamp.(sim.flow.σ'*sim.L/sim.U,-5,5), + linewidth=0,color=:RdBu_11,clims=(-5,5),levels=10, + aspect_ratio=:equal) +p2 = plot(sim.flow.u[Px,3:end-1,1],1:m-1,color=:blue,label="u",xlims=(-1,2)) +p3 = plot(1:n-1,sim.flow.u[3:end-1,Py1,1],color=:blue,label="u b",ylims=(-1,2)) +p = plot(p1,p2,p3,layout = l) +plot!(p[2],sim.flow.u[Px,3:end-1,2],1:m-1,color=:red,label="v") +plot!(p[3],1:n-1,sim.flow.u[3:end-1,Py1,2],color=:red,label="v b") +plot!(p[3],1:n-1,sim.flow.u[3:end-1,Py2,1],color=:blue,linestyle=:dot,label="u t") +plot!(p[3],1:n-1,-sim.flow.u[3:end-1,Py2,2],color=:red,linestyle=:dot,label="v t") + +@time anim = @animate 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_inner!(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) + p[1][1][:z] = clamp.(sim.flow.σ*sim.L/sim.U,-5,5) + + p[2][1][:x] = sim.flow.u[Px,3:end-1,1] + p[2][2][:x] = sim.flow.u[Px,3:end-1,2] + + p[3][1][:y] = sim.flow.u[3:end-1,Py1,1] + p[3][2][:y] = sim.flow.u[3:end-1,Py1,2] + p[3][3][:y] = sim.flow.u[3:end-1,Py2,1] + p[3][4][:y] = -sim.flow.u[3:end-1,Py2,2] + + # print time step + println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) +end +gif(anim, "Biot_Savart_Circle.gif") + +# N, = WaterLily.size_u(sim.flow.u) +# println("Inflow :", WaterLily.∮(sim.flow.u,N,2,2)) +# println("Outflow :", WaterLily.∮(sim.flow.u,N,N[2],2)) +# plot(sim.flow.u[:,2,2],label="v botton") +# plot!(-sim.flow.u[:,end,2],label="v top") +# plot!(1.0.-sim.flow.u[:,2,1],label="u bottom") +# plot!(1.0.-sim.flow.u[:,end,1],label="u top") + +# plot(sim.flow.u[1,:,1],label="u inlet") +# plot!(sim.flow.u[2,:,1],label="u inlet") + +# flood(sim.flow.u[:,:,1]) diff --git a/Biot_Savart_Profile.jl b/Biot_Savart_Profile.jl new file mode 100644 index 00000000..e3e7e1ca --- /dev/null +++ b/Biot_Savart_Profile.jl @@ -0,0 +1,96 @@ +using WaterLily +using Plots; gr() +include("BiotSavart.jl") + +# some definitons +U = 1 +Re = 250 +m, n = 2^8, 3*2^7 +radius, center = 32, 128+1.5 +Px, Py = 64, 28 + +# make a circle body +body = AutoBody((x,t)->√sum(abs2, x .- center) - 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) + +# duration of the simulation +duration = 10 +step = 0.1 +t₀ = 0.0 + +# N,_ = WaterLily.size_u(sim.flow.u) +u = copy(sim.flow.u) +N,_ = WaterLily.size_u(sim.flow.u) + +# init plot +l = @layout [ a ; b c ] +@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) +BiotSavart!(WaterLily.slice(N,Px,1),u,sim.flow.σ) +p1 = contourf(axes(sim.flow.σ,1),axes(sim.flow.σ,2),clamp.(sim.flow.σ*sim.L/sim.U,-5,5)', + linewidth=0,color=:RdBu_11,clims=(-5,5),levels=10, + aspect_ratio=:equal) +p2 = plot(sim.flow.u[Px,2:end-1,1],1:m,color=:blue,label="u",xlims=(-1,2)) +p3 = plot(1:n,sim.flow.u[2:end-1,Py,1],color=:blue,label="u b",ylims=(-1,2)) +p = plot(p1,p2,p3,layout = l) + +plot!(p[1],[Px,Px],[1,256],color=:black,linestyle=:dash,legend=:none) +plot!(p[1],[1,n],[Py,Py],color=:black,linestyle=:dash,legend=:none) +# plot!(p[1],[1,n],[228,228],color=:black,linestyle=:dash,legend=:none) + +plot!(p[2],sim.flow.u[Px,2:end-1,2],1:m,color=:red,label="v") +plot!(p[2],1.0.+u[Px,2:end-1,1],1:m,color=:blue,linestyle=:dashdot,label="u BS") +plot!(p[2], u[Px,2:end-1,2],1:m,linestyle=:dashdot,color=:red,label="v BS") + +BiotSavart!(WaterLily.slice(N,Py,2),u,sim.flow.σ) +BiotSavart!(WaterLily.slice(N,m-Py,2),u,sim.flow.σ) +plot!(p[3],1:n,sim.flow.u[2:end-1,Py,2],color=:red,label="v b") +plot!(p[3],1:n,1.0.+u[2:end-1,Py,1],color=:blue,linestyle=:dashdot,label="u b BS") +plot!(p[3],1:n, u[2:end-1,Py,2],color=:red,linestyle=:dashdot,label="v b BS") +plot!(p[3],1:n, u[2:end-1,m-Py,1],color=:blue,linestyle=:dot,label="u t BS") +plot!(p[3],1:n, u[2:end-1,m-Py,2],color=:red,linestyle=:dot,label="v t BS") + +@time anim = @animate 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!(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) + p[1][1][:z] = clamp.(sim.flow.σ*sim.L/sim.U,-5,5) + + BiotSavart!(WaterLily.slice(N,Px,1),u,sim.flow.σ) + p[2][1][:x] = sim.flow.u[Px,2:end-1,1] + p[2][2][:x] = sim.flow.u[Px,2:end-1,2] + p[2][3][:x] = 1.0.+u[Px,2:end-1,1] + p[2][4][:x] = u[Px,2:end-1,2] + # mass imbalance + ∮u = ∮(u,N,Px,1) + + BiotSavart!(WaterLily.slice(N,Py,2),u,sim.flow.σ) + BiotSavart!(WaterLily.slice(N,m-Py,2),u,sim.flow.σ) + p[3][1][:y] = sim.flow.u[2:end-1,Py,1] + p[3][2][:y] = sim.flow.u[2:end-1,Py,2] + p[3][3][:y] = 1.0.+u[2:end-1,Py,1] + p[3][4][:y] = u[2:end-1,Py,2] + p[3][5][:y] = 1.0.+u[2:end-1,m-Py,1] + p[3][6][:y] = u[2:end-1,m-Py,2] + # mass imbalance + ∮v = ∮(u,N,Py,2) + + # print time step + println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) + println("mass imbalance ",∮u,", ",∮v) +endo +gif(anim,"Biot_Savart_Profile.gif") diff --git a/Biot_Savart_Vortex.jl b/Biot_Savart_Vortex.jl new file mode 100644 index 00000000..6030b50d --- /dev/null +++ b/Biot_Savart_Vortex.jl @@ -0,0 +1,65 @@ +using WaterLily +using LinearAlgebra +include("examples/TwoD_plots.jl") +include("BiotSavart.jl") + +""" +RankineVortex(i,xy,center,R,Γ) +""" +function RankineVortex(i, xy, center, R=4, Γ=1) + xy = (xy .- 0.5 .- center) + x,y = xy + θ = atan(y,x) + r = norm(xy) + vθ =Γ/2π*(r<=R ? r/R^2 : 1/r) + v = [-vθ*sin(θ),vθ*cos(θ)] + return v[i] +end + + +# some definitons +U = 1 +Re = 250 +m, n = 2^6, 2^6 + +# make a simulation +sim = Simulation((n,m), (U,0), m; ν=U*m/Re, T=Float64) +u = copy(sim.flow.u) + +# make a Rankine vortex +f(i,x) = RankineVortex(i,x,(m/2,m/2),10, 1) + +# apply it to the flow +apply!(f, sim.flow.u) + +flood(sim.flow.u[:,:,1]; shift=(-0.5,-0.5)) +flood(sim.flow.u[:,:,2]; shift=(-0.5,-0.5)) + +# compute vorticity +@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) +flood(sim.flow.σ; shift=(-0.5,-0.5)) + +# compute velocity +BiotSavart!(u,sim.flow.σ) + +# error +BC!(sim.flow.u,zeros(2)) +BC!(u,zeros(2)) +println("L₂-norm error u-velocity ", WaterLily.L₂(u[:,:,1].-sim.flow.u[:,:,1])) +println("L₂-norm error v-velocity ", WaterLily.L₂(u[:,:,2].-sim.flow.u[:,:,2])) + +# check +flood(u[:,:,1].-sim.flow.u[:,:,1]; shift=(-0.5,-0.5)) +flood(u[:,:,2].-sim.flow.u[:,:,2]; shift=(-0.5,-0.5)) + + +# # show Kernel +# xs = collect(1:0.01:100) +# kernel = zeros(length(xs)) +# ϵ = 1e-6 +# for i in range(1,length(xs)) +# r = [xs[i],0.0] +# rⁿ = norm(r)^2 +# kernel[i] = r[1]/(2π*rⁿ+ϵ^2) +# end +# plot(xs, abs.(kernel), scale=:log10) From 6465360c88acdca4ec5735b9eec8e24c04b985d4 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 7 Sep 2023 10:50:01 +0200 Subject: [PATCH 2/6] initial convective exit commit --- Convective-Exit.jl | 106 +++++++++++++++++++++++++++++++++++++++++++++ src/Flow.jl | 8 ++-- src/util.jl | 53 +++++++++++++++++++++++ 3 files changed, 163 insertions(+), 4 deletions(-) create mode 100644 Convective-Exit.jl diff --git a/Convective-Exit.jl b/Convective-Exit.jl new file mode 100644 index 00000000..5bbc79fe --- /dev/null +++ b/Convective-Exit.jl @@ -0,0 +1,106 @@ +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)) + end +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) + + diff --git a/src/Flow.jl b/src/Flow.jl index bb2df7a8..41a15e64 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -110,12 +110,12 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp a.u⁰ .= a.u; a.u .= 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Δt!(a.u,a.U;Δt=a.Δt[end]) # convective exit + project!(a,b); BCΔt!(a.u,a.U) # 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); a.u ./= 2; BCΔt!(a.u,a.U) + project!(a,b,2); BCΔt!(a.u,a.U) push!(a.Δt,CFL(a)) end diff --git a/src/util.jl b/src/util.jl index e3eac80c..83f3fcf3 100644 --- a/src/util.jl +++ b/src/util.jl @@ -174,4 +174,57 @@ 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 \ No newline at end of file From 2d4787c19f1675d756acb4763539223a5999ea19 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 7 Sep 2023 10:52:01 +0200 Subject: [PATCH 3/6] removed old files --- BiotSavart.jl | 73 ------------------ Biot_Savart_BC.jl | 165 ----------------------------------------- Biot_Savart_Profile.jl | 96 ------------------------ Biot_Savart_Vortex.jl | 65 ---------------- 4 files changed, 399 deletions(-) delete mode 100644 BiotSavart.jl delete mode 100644 Biot_Savart_BC.jl delete mode 100644 Biot_Savart_Profile.jl delete mode 100644 Biot_Savart_Vortex.jl diff --git a/BiotSavart.jl b/BiotSavart.jl deleted file mode 100644 index 55122654..00000000 --- a/BiotSavart.jl +++ /dev/null @@ -1,73 +0,0 @@ -using LinearAlgebra -""" -Biot-Savart integral -""" -function BiotSavart!(u,ω) - N,n = WaterLily.size_u(u) - u .= 0.0 - for i ∈ 1:n - j=i%2+1 # the only component not zero in the vorticity - for I ∈ WaterLily.inside_u(N,i) - @WaterLily.loop u[I,i] += K(i,I,j,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) - end - end -end - -""" -Biot-Savart integral -""" -function BiotSavart!(Is,u,ω) - N,n = WaterLily.size_u(u) - for i ∈ 1:n - j=i%2+1 # the only component not zero in the vorticity - for I ∈ Is - u[I,i] = 0.0 - @WaterLily.loop u[I,i] += K(i,I,j,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) - end - end -end - -""" -Biot-Savart kernel -""" -function K(i,I,j,J,ϵ=1e-6) - # face centered velocity at I due to vorticity at cell edge J - r = loc(i,I) .- loc(0,J) .+ 0.5 # cell edge vorticity - rⁿ = norm(r)^2 - return sign(i-j)*r[j]/(2π*rⁿ+ϵ^2) -end - -function error_at_slice(ue,u,dn,up) - N,n = WaterLily.size_u(u) - error = 0.0 - for i ∈ 1:n, s ∈ [up,dn] - for I ∈ WaterLily.slice(N,s,2) - error += abs.(u[I,i] - ue[I,i]) - end - end - return error / (2N[1]) -end - -function velocity_at_slice(ue,dn,up) - res = zeros(eltype(ue),size(u)) - N,n = WaterLily.size_u(u) - for i ∈ 1:n, s ∈ [up,dn] - for I ∈ WaterLily.slice(N,s,2) - res[I,i] = ue[I,i] - end - end - return res -end - -function ∮(a::AbstractArray{T},j) where T - N,n = WaterLily.size_u(a) - return ∮(a,N,N[j],j) -end - -function ∮(a,N,s,j) - sm = 0.0 - for I ∈ WaterLily.slice(N.-1,s,j,2) # remove ghosts - sm += a[I,j] - end - return sm -end diff --git a/Biot_Savart_BC.jl b/Biot_Savart_BC.jl deleted file mode 100644 index 85df0322..00000000 --- a/Biot_Savart_BC.jl +++ /dev/null @@ -1,165 +0,0 @@ -using WaterLily -using LinearAlgebra -using Plots; gr() -# include("examples/TwoD_plots.jl") -include("BiotSavart.jl") - -# function ∮(a::AbstractArray{T},j) where T -# N,n = WaterLily.size_u(a) -# return ∮(a,N,N[j],j) -# end - -# function ∮(a,N,s,j) -# sm = 0.0 -# for I ∈ WaterLily.slice(N.-1,s,j,2) # remove ghosts -# sm += a[I,j] -# end -# return sm -# end - -function BC_BiotSavart!(a,A,ω) - N,n = WaterLily.size_u(a) - @inside ω[I] = WaterLily.curl(3,I,a) - area = ntuple(i->prod(N.-2)/(N[i]-2),n) # are of domain (without ghosts) - for j ∈ 1:n, i ∈ 1:n - k = i%2+1 # the only component not zero in the vorticity - if i==j # Normal direction, Dirichlet - if i==1 - for s ∈ (1,2,N[j]) - for I ∈ WaterLily.slice(N,s,j) - a[I,i] = A[i] - @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) - end - end - # @WaterLily.loop a[I,i] = A[i] over I ∈ WaterLily.slice(N,N[j],j) - else - for s ∈ (1,2,N[j]) - for I ∈ WaterLily.slice(N,s,j) - a[I,i] = A[i] - @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) - end - end - end - # make global face flux zero - # ∮u = (∮(a,N,N[j],j)-∮(a,N,2,j))/∮(ones((N...,j)),j)# mass flux imbalance in domain - # @WaterLily.loop a[I,i] += ∮u/2 over I ∈ WaterLily.slice(N,2,j) - # @WaterLily.loop a[I,i] -= ∮u/2 over I ∈ WaterLily.slice(N,N[j],j) - else # Tangential directions, Neumanns - # @WaterLily.loop a[I,i] = a[I+δ(j,I),i] over I ∈ WaterLily.slice(N,1,j) - # @WaterLily.loop a[I,i] = a[I-δ(j,I),i] over I ∈ WaterLily.slice(N,N[j],j) - for s ∈ (1,N[j]) - for I ∈ WaterLily.slice(N,s,j) - a[I,i] = A[i] - @WaterLily.loop a[I,i] += K(i,I,k,J)*ω[J] over J ∈ WaterLily.inside_u(N,i) - end - end - end - end -end - -BC!(a,A,σ) = BC_BiotSavart!(a,A,σ) -# BC!(a,A,σ) = WaterLily.BC!(a,A) - -@fastmath function mom_step_inner!(a::Flow,b::AbstractPoisson) - a.u⁰ .= a.u; a.u .= 0 - # predictor u → u' - WaterLily.conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν) - WaterLily.BDIM!(a); - @inside a.σ[I] = WaterLily.curl(3,I,a.u) - BC!(a.u,a.U,a.σ) - WaterLily.project!(a,b) - @inside a.σ[I] = WaterLily.curl(3,I,a.u) - BC!(a.u,a.U,a.σ) - - # corrector u → u¹ - WaterLily.conv_diff!(a.f,a.u,a.σ,ν=a.ν) - WaterLily.BDIM!(a); a.u ./= 2; - @inside a.σ[I] = WaterLily.curl(3,I,a.u) - BC!(a.u,a.U,a.σ) - b.x ./=2 #scaled pressure guess, no necessary - WaterLily.project!(a,b); - @inside a.σ[I] = WaterLily.curl(3,I,a.u) - BC!(a.u,a.U,a.σ) - push!(a.Δt,WaterLily.CFL(a)) -end - - -# some definitons -U = 1 -Re = 250 -m, n = 2^6, 2^7 -println("$m x $n: ", prod((m,n))) -radius = 16 -center = 32 + 1.5 - -# make a circle body -body = AutoBody((x,t)->√sum(abs2, x .- center) - 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=Float32) - -# duration of the simulation -duration = 0.0 -step = 0.1 -t₀ = 0.0 - -# init plot -N,dim = WaterLily.size_u(sim.flow.u) -Px = 1; Py1 = 2; Py2=N[2] -l = @layout [ a ; b c ] -@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) -p1 = contourf(axes(sim.flow.σ,1),axes(sim.flow.σ,2),clamp.(sim.flow.σ'*sim.L/sim.U,-5,5), - linewidth=0,color=:RdBu_11,clims=(-5,5),levels=10, - aspect_ratio=:equal) -p2 = plot(sim.flow.u[Px,3:end-1,1],1:m-1,color=:blue,label="u",xlims=(-1,2)) -p3 = plot(1:n-1,sim.flow.u[3:end-1,Py1,1],color=:blue,label="u b",ylims=(-1,2)) -p = plot(p1,p2,p3,layout = l) -plot!(p[2],sim.flow.u[Px,3:end-1,2],1:m-1,color=:red,label="v") -plot!(p[3],1:n-1,sim.flow.u[3:end-1,Py1,2],color=:red,label="v b") -plot!(p[3],1:n-1,sim.flow.u[3:end-1,Py2,1],color=:blue,linestyle=:dot,label="u t") -plot!(p[3],1:n-1,-sim.flow.u[3:end-1,Py2,2],color=:red,linestyle=:dot,label="v t") - -@time anim = @animate 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_inner!(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) - p[1][1][:z] = clamp.(sim.flow.σ*sim.L/sim.U,-5,5) - - p[2][1][:x] = sim.flow.u[Px,3:end-1,1] - p[2][2][:x] = sim.flow.u[Px,3:end-1,2] - - p[3][1][:y] = sim.flow.u[3:end-1,Py1,1] - p[3][2][:y] = sim.flow.u[3:end-1,Py1,2] - p[3][3][:y] = sim.flow.u[3:end-1,Py2,1] - p[3][4][:y] = -sim.flow.u[3:end-1,Py2,2] - - # print time step - println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) -end -gif(anim, "Biot_Savart_Circle.gif") - -# N, = WaterLily.size_u(sim.flow.u) -# println("Inflow :", WaterLily.∮(sim.flow.u,N,2,2)) -# println("Outflow :", WaterLily.∮(sim.flow.u,N,N[2],2)) -# plot(sim.flow.u[:,2,2],label="v botton") -# plot!(-sim.flow.u[:,end,2],label="v top") -# plot!(1.0.-sim.flow.u[:,2,1],label="u bottom") -# plot!(1.0.-sim.flow.u[:,end,1],label="u top") - -# plot(sim.flow.u[1,:,1],label="u inlet") -# plot!(sim.flow.u[2,:,1],label="u inlet") - -# flood(sim.flow.u[:,:,1]) diff --git a/Biot_Savart_Profile.jl b/Biot_Savart_Profile.jl deleted file mode 100644 index e3e7e1ca..00000000 --- a/Biot_Savart_Profile.jl +++ /dev/null @@ -1,96 +0,0 @@ -using WaterLily -using Plots; gr() -include("BiotSavart.jl") - -# some definitons -U = 1 -Re = 250 -m, n = 2^8, 3*2^7 -radius, center = 32, 128+1.5 -Px, Py = 64, 28 - -# make a circle body -body = AutoBody((x,t)->√sum(abs2, x .- center) - 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) - -# duration of the simulation -duration = 10 -step = 0.1 -t₀ = 0.0 - -# N,_ = WaterLily.size_u(sim.flow.u) -u = copy(sim.flow.u) -N,_ = WaterLily.size_u(sim.flow.u) - -# init plot -l = @layout [ a ; b c ] -@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) -BiotSavart!(WaterLily.slice(N,Px,1),u,sim.flow.σ) -p1 = contourf(axes(sim.flow.σ,1),axes(sim.flow.σ,2),clamp.(sim.flow.σ*sim.L/sim.U,-5,5)', - linewidth=0,color=:RdBu_11,clims=(-5,5),levels=10, - aspect_ratio=:equal) -p2 = plot(sim.flow.u[Px,2:end-1,1],1:m,color=:blue,label="u",xlims=(-1,2)) -p3 = plot(1:n,sim.flow.u[2:end-1,Py,1],color=:blue,label="u b",ylims=(-1,2)) -p = plot(p1,p2,p3,layout = l) - -plot!(p[1],[Px,Px],[1,256],color=:black,linestyle=:dash,legend=:none) -plot!(p[1],[1,n],[Py,Py],color=:black,linestyle=:dash,legend=:none) -# plot!(p[1],[1,n],[228,228],color=:black,linestyle=:dash,legend=:none) - -plot!(p[2],sim.flow.u[Px,2:end-1,2],1:m,color=:red,label="v") -plot!(p[2],1.0.+u[Px,2:end-1,1],1:m,color=:blue,linestyle=:dashdot,label="u BS") -plot!(p[2], u[Px,2:end-1,2],1:m,linestyle=:dashdot,color=:red,label="v BS") - -BiotSavart!(WaterLily.slice(N,Py,2),u,sim.flow.σ) -BiotSavart!(WaterLily.slice(N,m-Py,2),u,sim.flow.σ) -plot!(p[3],1:n,sim.flow.u[2:end-1,Py,2],color=:red,label="v b") -plot!(p[3],1:n,1.0.+u[2:end-1,Py,1],color=:blue,linestyle=:dashdot,label="u b BS") -plot!(p[3],1:n, u[2:end-1,Py,2],color=:red,linestyle=:dashdot,label="v b BS") -plot!(p[3],1:n, u[2:end-1,m-Py,1],color=:blue,linestyle=:dot,label="u t BS") -plot!(p[3],1:n, u[2:end-1,m-Py,2],color=:red,linestyle=:dot,label="v t BS") - -@time anim = @animate 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!(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) - p[1][1][:z] = clamp.(sim.flow.σ*sim.L/sim.U,-5,5) - - BiotSavart!(WaterLily.slice(N,Px,1),u,sim.flow.σ) - p[2][1][:x] = sim.flow.u[Px,2:end-1,1] - p[2][2][:x] = sim.flow.u[Px,2:end-1,2] - p[2][3][:x] = 1.0.+u[Px,2:end-1,1] - p[2][4][:x] = u[Px,2:end-1,2] - # mass imbalance - ∮u = ∮(u,N,Px,1) - - BiotSavart!(WaterLily.slice(N,Py,2),u,sim.flow.σ) - BiotSavart!(WaterLily.slice(N,m-Py,2),u,sim.flow.σ) - p[3][1][:y] = sim.flow.u[2:end-1,Py,1] - p[3][2][:y] = sim.flow.u[2:end-1,Py,2] - p[3][3][:y] = 1.0.+u[2:end-1,Py,1] - p[3][4][:y] = u[2:end-1,Py,2] - p[3][5][:y] = 1.0.+u[2:end-1,m-Py,1] - p[3][6][:y] = u[2:end-1,m-Py,2] - # mass imbalance - ∮v = ∮(u,N,Py,2) - - # print time step - println("tU/L=",round(tᵢ,digits=4),", Δt=",round(sim.flow.Δt[end],digits=3)) - println("mass imbalance ",∮u,", ",∮v) -endo -gif(anim,"Biot_Savart_Profile.gif") diff --git a/Biot_Savart_Vortex.jl b/Biot_Savart_Vortex.jl deleted file mode 100644 index 6030b50d..00000000 --- a/Biot_Savart_Vortex.jl +++ /dev/null @@ -1,65 +0,0 @@ -using WaterLily -using LinearAlgebra -include("examples/TwoD_plots.jl") -include("BiotSavart.jl") - -""" -RankineVortex(i,xy,center,R,Γ) -""" -function RankineVortex(i, xy, center, R=4, Γ=1) - xy = (xy .- 0.5 .- center) - x,y = xy - θ = atan(y,x) - r = norm(xy) - vθ =Γ/2π*(r<=R ? r/R^2 : 1/r) - v = [-vθ*sin(θ),vθ*cos(θ)] - return v[i] -end - - -# some definitons -U = 1 -Re = 250 -m, n = 2^6, 2^6 - -# make a simulation -sim = Simulation((n,m), (U,0), m; ν=U*m/Re, T=Float64) -u = copy(sim.flow.u) - -# make a Rankine vortex -f(i,x) = RankineVortex(i,x,(m/2,m/2),10, 1) - -# apply it to the flow -apply!(f, sim.flow.u) - -flood(sim.flow.u[:,:,1]; shift=(-0.5,-0.5)) -flood(sim.flow.u[:,:,2]; shift=(-0.5,-0.5)) - -# compute vorticity -@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u) -flood(sim.flow.σ; shift=(-0.5,-0.5)) - -# compute velocity -BiotSavart!(u,sim.flow.σ) - -# error -BC!(sim.flow.u,zeros(2)) -BC!(u,zeros(2)) -println("L₂-norm error u-velocity ", WaterLily.L₂(u[:,:,1].-sim.flow.u[:,:,1])) -println("L₂-norm error v-velocity ", WaterLily.L₂(u[:,:,2].-sim.flow.u[:,:,2])) - -# check -flood(u[:,:,1].-sim.flow.u[:,:,1]; shift=(-0.5,-0.5)) -flood(u[:,:,2].-sim.flow.u[:,:,2]; shift=(-0.5,-0.5)) - - -# # show Kernel -# xs = collect(1:0.01:100) -# kernel = zeros(length(xs)) -# ϵ = 1e-6 -# for i in range(1,length(xs)) -# r = [xs[i],0.0] -# rⁿ = norm(r)^2 -# kernel[i] = r[1]/(2π*rⁿ+ϵ^2) -# end -# plot(xs, abs.(kernel), scale=:log10) From 63b88be4ef2a50fe35d1ff2d5f58d078a83540e0 Mon Sep 17 00:00:00 2001 From: weymouth Date: Tue, 10 Oct 2023 22:04:33 +0200 Subject: [PATCH 4/6] 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 | 133 +++++++++++---------------------------------- Project.toml | 1 + src/Flow.jl | 21 ++++--- src/WaterLily.jl | 4 +- src/util.jl | 82 ++++++---------------------- 5 files changed, 63 insertions(+), 178 deletions(-) diff --git a/Convective-Exit.jl b/Convective-Exit.jl index 5bbc79fe..9207ff6c 100644 --- a/Convective-Exit.jl +++ b/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)) \ No newline at end of file diff --git a/Project.toml b/Project.toml index 725575ee..61ffe13e 100644 --- a/Project.toml +++ b/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" diff --git a/src/Flow.jl b/src/Flow.jl index 41a15e64..ddf48e9d 100644 --- a/src/Flow.jl +++ b/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) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 399d9c54..657f16f1 100644 --- a/src/WaterLily.jl +++ b/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 diff --git a/src/util.jl b/src/util.jl index 83f3fcf3..43b13e9a 100644 --- a/src/util.jl +++ b/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,13 +146,14 @@ 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) @@ -164,6 +161,14 @@ function BC!(a,A) 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 \ No newline at end of file From cdb354c41eca0f5a8f225a0d6357701386b6b3ba Mon Sep 17 00:00:00 2001 From: weymouth Date: Thu, 12 Oct 2023 22:09:19 +0200 Subject: [PATCH 5/6] Move to test --- Convective-Exit.jl | 35 ----------------------------------- Project.toml | 1 - test/runtests.jl | 28 +++++++++++++++++----------- 3 files changed, 17 insertions(+), 47 deletions(-) delete mode 100644 Convective-Exit.jl diff --git a/Convective-Exit.jl b/Convective-Exit.jl deleted file mode 100644 index 9207ff6c..00000000 --- a/Convective-Exit.jl +++ /dev/null @@ -1,35 +0,0 @@ -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 - -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)) \ No newline at end of file diff --git a/Project.toml b/Project.toml index 61ffe13e..725575ee 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,6 @@ 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" diff --git a/test/runtests.jl b/test/runtests.jl index ac4c6c9c..be7aaae5 100644 --- a/test/runtests.jl +++ b/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, exit=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, exit) end @testset "WaterLily.jl" begin - for mem ∈ arrays - sim = sphere_sim(32;mem); + for mem ∈ arrays, exit ∈ (true,false) + sim = sphere_sim(;mem,exit); @test sim_time(sim) == 0 sim_step!(sim,0.1,remeasure=false) @test length(sim.flow.Δt)-1 == length(sim.pois.n)÷2 From 62793e0c86933ba2685e160ed28a6f996694508f Mon Sep 17 00:00:00 2001 From: Gabriel Weymouth Date: Fri, 13 Oct 2023 17:02:57 +0100 Subject: [PATCH 6/6] Change exit to exitBC to avoid keyword conflict --- src/Flow.jl | 18 +++++++++--------- src/WaterLily.jl | 4 ++-- src/util.jl | 8 ++++---- test/runtests.jl | 8 ++++---- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index ddf48e9d..c79a567d 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -70,18 +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 - 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 + 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, exit); exitBC!(u,u,U,0.) + 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],ν,exit) + new{D,T,typeof(p),typeof(u),typeof(μ₁)}(u,u⁰,fv,p,σ,V,σᵥ,μ₀,μ₁,U,T[Δt],ν,exitBC) end end @@ -111,13 +111,13 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp 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,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) + 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); scale_u!(a,0.5); BC!(a.u,a.U,a.exit) - project!(a,b,0.5); BC!(a.u,a.U,a.exit) + 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)) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 657f16f1..74e6be9d 100644 --- a/src/WaterLily.jl +++ b/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],exit=false, + 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,exit) + 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 diff --git a/src/util.jl b/src/util.jl index 43b13e9a..8a4cf645 100644 --- a/src/util.jl +++ b/src/util.jl @@ -163,10 +163,10 @@ 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 + 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 """ diff --git a/test/runtests.jl b/test/runtests.jl index be7aaae5..ed1e6619 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -211,13 +211,13 @@ end end end -function sphere_sim(radius = 8; mem=Array, exit=false) +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, exit) + return Simulation(radius.*(6,4),(1,0),radius; body, ν=radius/250, T=Float32, mem, exitBC) end @testset "WaterLily.jl" begin - for mem ∈ arrays, exit ∈ (true,false) - sim = sphere_sim(;mem,exit); + 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