From 63b88be4ef2a50fe35d1ff2d5f58d078a83540e0 Mon Sep 17 00:00:00 2001 From: weymouth Date: Tue, 10 Oct 2023 22:04:33 +0200 Subject: [PATCH] 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