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