diff --git a/src/Flow.jl b/src/Flow.jl index 1113a94a..4650a007 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -59,26 +59,16 @@ lowerBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop ( Φ[I] = ϕuP(j,CIj(j,CI(I,i),N[j]-2),CI(I,i),u,ϕ(i,CI(I,j),u)) -ν*∂(j,CI(I,i),u); r[I,i] += Φ[I]) over I ∈ slice(N,2,j,2) upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I,2)] over I ∈ slice(N,N[j],j,2) -using EllipsisNotation """ - accelerate!(r,dt,g) + accelerate!(r,t,g,U) -Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. +Accounts for applied and reference-frame acceleration using `rᵢ += g(i,x,t)+dU(i,x,t)/dt` """ -accelerate!(r,dt,g::Function,::Tuple,t=sum(dt)) = for i ∈ 1:last(size(r)) - r[..,i] .+= g(i,t) -end -accelerate!(r,dt,g::Nothing,U::Function) = accelerate!(r,dt,(i,t)->ForwardDiff.derivative(τ->U(i,τ),t),()) -accelerate!(r,dt,g::Function,U::Function) = accelerate!(r,dt,(i,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,τ),t),()) -accelerate!(r,dt,::Nothing,::Tuple) = nothing -""" - BCTuple(U,dt,N) - -Return BC tuple `U(i∈1:N, t=sum(dt))`. -""" -BCTuple(f::Function,dt,N,t=sum(dt))=ntuple(i->f(i,t),N) -BCTuple(f::Tuple,dt,N)=f - +accelerate!(r,t,::Nothing,::Union{Nothing,Tuple}) = nothing +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)) """ Flow{D::Int, T::Float, Sf<:AbstractArray{T,D}, Vf<:AbstractArray{T,D+1}, Tf<:AbstractArray{T,D+2}} @@ -101,23 +91,25 @@ struct Flow{D, T, Sf<:AbstractArray{T}, Vf<:AbstractArray{T}, Tf<:AbstractArray{ μ₀:: Vf # zeroth-moment vector μ₁:: Tf # first-moment tensor field # Non-fields - U :: Union{NTuple{D,Number},Function} # domain boundary values + uBC :: Union{NTuple{D,Number},Function} # domain boundary values/function Δt:: Vector{T} # time step (stored in CPU memory) ν :: T # kinematic viscosity - g :: Union{Function,Nothing} # (possibly time-varying) uniform acceleration field + g :: Union{Function,Nothing} # acceleration field funciton exitBC :: Bool # Convection exit perdir :: NTuple # tuple of periodic direction - function Flow(N::NTuple{D}, U; f=Array, Δt=0.25, ν=0., g=nothing, - uλ::Function=(i, x) -> 0., perdir=(), exitBC=false, T=Float64) where D + function Flow(N::NTuple{D}, uBC; f=Array, Δt=0.25, ν=0., g=nothing, + uλ=nothing, perdir=(), exitBC=false, T=Float32) where D Ng = N .+ 2 Nd = (Ng..., D) - u = Array{T}(undef, Nd...) |> f; apply!(uλ, u); - BC!(u,BCTuple(U,0.,D),exitBC,perdir); exitBC!(u,u,BCTuple(U,0.,D),0.) - u⁰ = copy(u); + isnothing(uλ) && (uλ = ic_function(uBC)) + u = Array{T}(undef, Nd...) |> f + isa(uλ, Function) ? apply!(uλ, u) : apply!((i,x)->uλ[i], u) + BC!(u,uBC,exitBC,perdir); exitBC!(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, ones(T, Nd) |> f, zeros(T, Ng..., D, D) |> f BC!(μ₀,ntuple(zero, D),false,perdir) - new{D,T,typeof(p),typeof(u),typeof(μ₁)}(u,u⁰,fv,p,σ,V,μ₀,μ₁,U,T[Δt],ν,g,exitBC,perdir) + new{D,T,typeof(p),typeof(u),typeof(μ₁)}(u,u⁰,fv,p,σ,V,μ₀,μ₁,uBC,T[Δt],ν,g,exitBC,perdir) end end @@ -150,21 +142,23 @@ end Integrate the `Flow` one time step using the [Boundary Data Immersion Method](https://eprints.soton.ac.uk/369635/) and the `AbstractPoisson` pressure solver to project the velocity onto an incompressible flow. """ -@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson) where N - a.u⁰ .= a.u; scale_u!(a,0); U = BCTuple(a.U,a.Δt,N) +@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson;udf=nothing,kwargs...) where N + a.u⁰ .= a.u; scale_u!(a,0); t₁ = sum(a.Δt); t₀ = t₁-a.Δt[end] # predictor u → u' @log "p" conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir) - accelerate!(a.f,@view(a.Δt[1:end-1]),a.g,a.U) - BDIM!(a); BC!(a.u,U,a.exitBC,a.perdir) - a.exitBC && exitBC!(a.u,a.u⁰,U,a.Δt[end]) # convective exit - project!(a,b); BC!(a.u,U,a.exitBC,a.perdir) + udf!(a,udf,t₀; kwargs...) + accelerate!(a.f,t₀,a.g,a.uBC) + BDIM!(a); BC!(a.u,a.uBC,a.exitBC,a.perdir,t₁) # BC MUST be at t₁ + a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit + project!(a,b); BC!(a.u,a.uBC,a.exitBC,a.perdir,t₁) # corrector u → u¹ @log "c" conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) - accelerate!(a.f,a.Δt,a.g,a.U) - BDIM!(a); scale_u!(a,0.5); BC!(a.u,U,a.exitBC,a.perdir) - project!(a,b,0.5); BC!(a.u,U,a.exitBC,a.perdir) + udf!(a,udf,t₁; kwargs...) + accelerate!(a.f,t₁,a.g,a.uBC) + BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.uBC,a.exitBC,a.perdir,t₁) + project!(a,b,0.5); BC!(a.u,a.uBC,a.exitBC,a.perdir,t₁) push!(a.Δt,CFL(a)) end scale_u!(a,scale) = @loop a.u[Ii] *= scale over Ii ∈ inside_u(size(a.p)) @@ -180,3 +174,12 @@ end end return s end + +""" + udf!(flow::Flow,udf::Function,t) + +User defined function using `udf::Function` to operate on `flow::Flow` during the predictor and corrector step, in sync with time `t`. +Keyword arguments must be passed to `sim_step!` for them to be carried over the actual function call. +""" +udf!(flow,::Nothing,t; kwargs...) = nothing +udf!(flow,force!::Function,t; kwargs...) = force!(flow,t; kwargs...) diff --git a/src/Metrics.jl b/src/Metrics.jl index c5011299..b8f6afe7 100644 --- a/src/Metrics.jl +++ b/src/Metrics.jl @@ -100,12 +100,12 @@ function pressure_force(p,df,body,t=0,T=promote_type(Float64,eltype(p))) end """ - ∇²u(I::CartesianIndex,u) + S(I::CartesianIndex,u) Rate-of-strain tensor. """ -∇²u(I::CartesianIndex{2},u) = @SMatrix [∂(i,j,I,u)+∂(j,i,I,u) for i ∈ 1:2, j ∈ 1:2] -∇²u(I::CartesianIndex{3},u) = @SMatrix [∂(i,j,I,u)+∂(j,i,I,u) for i ∈ 1:3, j ∈ 1:3] +S(I::CartesianIndex{2},u) = @SMatrix [0.5*(∂(i,j,I,u)+∂(j,i,I,u)) for i ∈ 1:2, j ∈ 1:2] +S(I::CartesianIndex{3},u) = @SMatrix [0.5*(∂(i,j,I,u)+∂(j,i,I,u)) for i ∈ 1:3, j ∈ 1:3] """ viscous_force(sim::Simulation) @@ -115,7 +115,7 @@ viscous_force(sim) = viscous_force(sim.flow,sim.body) viscous_force(flow,body) = viscous_force(flow.u,flow.ν,flow.f,body,time(flow)) function viscous_force(u,ν,df,body,t=0,T=promote_type(Float64,eltype(u))) df .= zero(eltype(u)) - @loop df[I,:] .= -ν*∇²u(I,u)*nds(body,loc(0,I,T),t) over I ∈ inside_u(u) + @loop df[I,:] .= -2ν*S(I,u)*nds(body,loc(0,I,T),t) over I ∈ inside_u(u) sum(T,df,dims=ntuple(i->i,ndims(u)-1))[:] |> Array end diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 202dbab5..24e721f0 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -30,26 +30,27 @@ include("Metrics.jl") abstract type AbstractSimulation end """ - Simulation(dims::NTuple, u_BC::Union{NTuple,Function}, L::Number; - U=norm2(u_BC), Δt=0.25, ν=0., ϵ=1, perdir=() - uλ::nothing, g=nothing, exitBC=false, + Simulation(dims::NTuple, uBC::Union{NTuple,Function}, L::Number; + U=norm2(Uλ), Δt=0.25, ν=0., ϵ=1, g=nothing, + perdir=(), exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) Constructor for a WaterLily.jl simulation: - `dims`: Simulation domain dimensions. - - `u_BC`: Simulation domain velocity boundary conditions, either a - tuple `u_BC[i]=uᵢ, i=eachindex(dims)`, or a time-varying function `f(i,t)` + - `uBC`: Velocity field applied to boundary and acceleration conditions. + Define a `Tuple` for constant BCs, or a `Function` for space and time varying BCs `uBC(i,x,t)`. - `L`: Simulation length scale. - - `U`: Simulation velocity scale. + - `U`: Simulation velocity scale. Required if using `Uλ::Function`. - `Δt`: Initial time step. - `ν`: Scaled viscosity (`Re=UL/ν`). - - `g`: Domain acceleration, `g(i,t)=duᵢ/dt` + - `g`: Domain acceleration, `g(i,x,t)=duᵢ/dt` - `ϵ`: BDIM kernel width. - `perdir`: Domain periodic boundary condition in the `(i,)` direction. + - `uλ`: Velocity field applied to the initial condition. + Define a Tuple for homogeneous (per direction) IC, or a `Function` for space varying IC `uλ(i,x)`. - `exitBC`: Convective exit boundary condition in the `i=1` direction. - - `uλ`: Function to generate the initial velocity field. - `body`: Immersed geometry. - `T`: Array element type. - `mem`: memory location. `Array`, `CuArray`, `ROCm` to run on CPU, NVIDIA, or AMD devices, respectively. @@ -63,16 +64,14 @@ mutable struct Simulation <: AbstractSimulation flow :: Flow body :: AbstractBody pois :: AbstractPoisson - function Simulation(dims::NTuple{N}, u_BC, L::Number; + function Simulation(dims::NTuple{N}, uBC, L::Number; Δ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 !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" - isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->u_BC(i,0.),(i,x)->u_BC[i]) : uλ - U = isnothing(U) ? √sum(abs2,u_BC) : U # default if not specified - flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) + @assert !(isnothing(U) && isa(uBC,Function)) "`U` (velocity scale) must be specified if boundary conditions `uBC` is a `Function`" + isnothing(U) && (U = √sum(abs2,uBC)) + check_fn(uBC,N,T,3); check_fn(g,N,T,3); check_fn(uλ,N,T,2) + flow = Flow(dims,uBC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) end @@ -94,18 +93,20 @@ sim_time(sim::AbstractSimulation) = time(sim)*sim.U/sim.L Integrate the simulation `sim` up to dimensionless time `t_end`. If `remeasure=true`, the body is remeasured at every time step. Can be set to `false` for static geometries to speed up simulation. +A user-defined function `udf` can be passed to arbitrarily modify the `::Flow` during the predictor and corrector steps. +If the `udf` user keyword arguments, these needs to be included in the `sim_step!` call as well. """ -function sim_step!(sim::AbstractSimulation,t_end;remeasure=true,max_steps=typemax(Int),verbose=false) +function sim_step!(sim::AbstractSimulation,t_end;remeasure=true,max_steps=typemax(Int),udf=nothing,verbose=false,kwargs...) steps₀ = length(sim.flow.Δt) while sim_time(sim) < t_end && length(sim.flow.Δt) - steps₀ < max_steps - sim_step!(sim; remeasure) + sim_step!(sim; remeasure, udf, kwargs...) verbose && println("tU/L=",round(sim_time(sim),digits=4), ", Δt=",round(sim.flow.Δt[end],digits=3)) end end -function sim_step!(sim::AbstractSimulation;remeasure=true) +function sim_step!(sim::AbstractSimulation;remeasure=true,udf=nothing,kwargs...) remeasure && measure!(sim) - mom_step!(sim.flow,sim.pois) + mom_step!(sim.flow, sim.pois; udf, kwargs...) end """ diff --git a/src/util.jl b/src/util.jl index 0c12c598..25c05189 100644 --- a/src/util.jl +++ b/src/util.jl @@ -189,7 +189,8 @@ 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,saveexit=false,perdir=()) +BC!(a,U,saveexit=false,perdir=(),t=0) = BC!(a,(i,x,t)->U[i],saveexit,perdir,t) +function BC!(a,uBC::Function,saveexit=false,perdir=(),t=0) N,n = size_u(a) for i ∈ 1:n, j ∈ 1:n if j in perdir @@ -198,26 +199,28 @@ function BC!(a,A,saveexit=false,perdir=()) else if i==j # Normal direction, Dirichlet for s ∈ (1,2) - @loop a[I,i] = A[i] over I ∈ slice(N,s,j) + @loop a[I,i] = uBC(i,loc(i,I),t) 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 + (!saveexit || i>1) && (@loop a[I,i] = uBC(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] = uBC(i,loc(i,I),t)+a[I+δ(j,I),i]-uBC(i,loc(i,I+δ(j,I)),t) over I ∈ slice(N,1,j) + @loop a[I,i] = uBC(i,loc(i,I),t)+a[I-δ(j,I),i]-uBC(i,loc(i,I-δ(j,I)),t) over I ∈ slice(N,N[j],j) end end end end + """ exitBC!(u,u⁰,U,Δt) Apply a 1D convection scheme to fill the ghost cell on the exit of the domain. """ -function exitBC!(u,u⁰,U,Δt) +function exitBC!(u,u⁰,Δt) N,_ = size_u(u) 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 + U = sum(@view(u[slice(N.-1,2,1,2),1]))/length(exitR) # inflow mass flux + @loop u[I,1] = u⁰[I,1]-U*Δt*(u⁰[I,1]-u⁰[I-δ(1,I),1]) over I ∈ exitR + ∮u = sum(@view(u[exitR,1]))/length(exitR)-U # mass flux imbalance @loop u[I,1] -= ∮u over I ∈ exitR # correct flux end """ @@ -250,8 +253,20 @@ function interp(x::SVector{D}, arr::AbstractArray{T,D}) where {D,T} end return s end +using EllipsisNotation function interp(x::SVector{D}, varr::AbstractArray) where {D} # Shift to align with each staggered grid component and interpolate @inline shift(i) = SVector{D}(ifelse(i==j,0.5,0.0) for j in 1:D) return SVector{D}(interp(x+shift(i),@view(varr[..,i])) for i in 1:D) -end \ No newline at end of file +end + +check_fn(f,N,T,nargs) = nothing +function check_fn(f::Function,N,T,nargs) + @assert first(methods(f)).nargs==nargs+1 "$f signature needs $nargs arguments" + @assert all(typeof.(ntuple(i->f(i,xtargs(Val{}(nargs),N,T)...),N)).==T) "$f is not type stable" +end +xtargs(::Val{2},N,T) = (zeros(SVector{N,T}),) +xtargs(::Val{3},N,T) = (zeros(SVector{N,T}),zero(T)) + +ic_function(uBC::Function) = (i,x)->uBC(i,x,0) +ic_function(uBC::Tuple) = (i,x)->uBC[i] \ No newline at end of file diff --git a/test/maintests.jl b/test/maintests.jl index d46c99eb..282b407a 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -42,9 +42,20 @@ using ReadVTK, WriteVTK BC!(u,U,true) # save exit values @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3) - WaterLily.exitBC!(u,u,U,0) # conservative exit check + WaterLily.exitBC!(u,u,0) # conservative exit check @test GPUArrays.@allowscalar all(u[end,2:end-1, 1] .== U[1]) + # test BC with function + 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 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 + @test GPUArrays.@allowscalar all(v[end, :, 1] .== 3) + BC!(u,U,true,(2,)) # periodic in y and save exit values @test GPUArrays.@allowscalar all(u[:, 1:2, 1] .== u[:, end-1:end, 1]) && all(u[:, 1:2, 1] .== u[:,end-1:end,1]) WaterLily.perBC!(σ,(1,2)) # periodic in two directions @@ -54,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 @@ -154,20 +179,30 @@ end Ip = WaterLily.CIj(1,I,length(f)-2); # make periodic @test ϕuP(1,Ip,I,f,1)==λ(f[Ip],f[I-δ(1,I)],f[I]) - @test all(WaterLily.BCTuple((1,2,3),[0],3).==WaterLily.BCTuple((i,t)->i,0,3)) - @test all(WaterLily.BCTuple((i,t)->t,[1.234],3).==ntuple(i->1.234,3)) - # check applying acceleration for f ∈ arrays N = 4; a = zeros(N,N,2) |> f - WaterLily.accelerate!(a,[1],nothing,()) + WaterLily.accelerate!(a,1,nothing,()) @test all(a .== 0) - WaterLily.accelerate!(a,[1],(i,t) -> i==1 ? t : 2*t,()) + WaterLily.accelerate!(a,1.,(i,x,t)->i==1 ? t : 2*t,()) @test all(a[:,:,1] .== 1) && all(a[:,:,2] .== 2) - WaterLily.accelerate!(a,[1],nothing,(i,t) -> i==1 ? -t : -2*t) + WaterLily.accelerate!(a,1.,nothing,(i,x,t) -> i==1 ? -t : -2*t) @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) - WaterLily.accelerate!(a,[1],(i,t) -> i==1 ? t : 2*t,(i,t) -> i==1 ? -t : -2*t) + WaterLily.accelerate!(a,1.,(i,x,t) -> i==1 ? t : 2*t,(i,x,t) -> i==1 ? -t : -2*t) @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) + # check applying body force (changes in x but not t) + b = zeros(N,N,2) |> f + WaterLily.accelerate!(b,0.,(i,x,t)->1,nothing) + @test all(b .== 1) + WaterLily.accelerate!(b,1.,(i,x,t)->0,(i,x,t)->t) + @test all(b .== 2) + a .= 0; b .= 1 # reset and accelerate using a non-uniform velocity field + WaterLily.accelerate!(a,0.,nothing,(i,x,t)->t*(x[i]+1.0)) + WaterLily.accelerate!(b,0,(i,x,t)->x[i],nothing) + @test all(b .== a) + WaterLily.accelerate!(b,1.,(i,x,t)->x[i]+1.0,nothing) + WaterLily.accelerate!(a,1.,nothing,(i,x,t)->t*(x[i]+1.0)) + @test all(b .== a) end # Impulsive flow in a box U = (2/3, -1/3) @@ -231,7 +266,7 @@ end function TGVsim(mem;perdir=(1,2),Re=1e8,T=typeof(Re)) # Define vortex size, velocity, viscosity - L = 64; κ=2π/L; ν = 1/(κ*Re); + L = 64; κ = T(2π/L); ν = T(1/(κ*Re)); # TGV vortex in 2D function TGV(i,xy,t,κ,ν) x,y = @. (xy)*κ # scaled coordinates @@ -239,7 +274,7 @@ function TGVsim(mem;perdir=(1,2),Re=1e8,T=typeof(Re)) return cos(x)*sin(y)*exp(-2*κ^2*ν*t) # u_y end # Initialize simulation - return Simulation((L,L),(0,0),L;U=1,uλ=(i,x)->TGV(i,x,0.0,κ,ν),ν,T,mem,perdir),TGV + return Simulation((L,L),(i,x,t)->TGV(i,x,t,κ,ν),L;U=1,ν,T,mem,perdir),TGV end @testset "Flow.jl periodic TGV" begin for f ∈ arrays @@ -277,34 +312,76 @@ end @test derivative(lift,2.0) ≈ (lift(2+h)-lift(2-h))/2h rtol=√h end -function acceleratingFlow(N;T=Float64,perdir=(1,),jerk=4,mem=Array) +function acceleratingFlow(N;use_g=false,T=Float64,perdir=(1,),jerk=4,mem=Array) # periodic in x, Neumann in y # assuming gravitational scale is 1 and Fr is 1, U scale is Fr*√gL UScale = √N # this is also initial U # constant jerk in x, zero acceleration in y - g(i,t) = i==1 ? t*jerk : 0 + g(i,x,t) = i==1 ? t*jerk : 0. + !use_g && (g = nothing) return WaterLily.Simulation( (N,N), (UScale,0.), N; ν=0.001,g,Δt=0.001,perdir,T,mem ),jerk end +gravity!(flow::Flow,t; jerk=4) = for i ∈ 1:last(size(flow.f)) + WaterLily.@loop flow.f[I,i] += i==1 ? t*jerk : 0 over I ∈ CartesianIndices(Base.front(size(flow.f))) +end @testset "Flow.jl with increasing body force" begin for f ∈ arrays N = 8 - sim,jerk = acceleratingFlow(N;mem=f) + sim,jerk = acceleratingFlow(N;use_g=true,mem=f) sim_step!(sim,1.0); u = sim.flow.u |> Array # Exact uₓ = uₓ₀ + ∫ a dt = uₓ₀ + ∫ jerk*t dt = uₓ₀ + 0.5*jerk*t^2 - uFinal = sim.flow.U[1] + 0.5*jerk*WaterLily.time(sim)^2 + uFinal = sim.flow.uBC[1] + 0.5*jerk*WaterLily.time(sim)^2 @test ( WaterLily.L₂(u[:,:,1].-uFinal) < 1e-4 && WaterLily.L₂(u[:,:,2].-0) < 1e-4 ) + + # Test with user defined function instead of acceleration + sim_udf,_ = acceleratingFlow(N;mem=f) + sim_step!(sim_udf,1.0; udf=gravity!, jerk=jerk); u_udf = sim_udf.flow.u |> Array + uFinal = sim_udf.flow.uBC[1] + 0.5*jerk*WaterLily.time(sim_udf)^2 + @test ( + WaterLily.L₂(u_udf[:,:,1].-uFinal) < 1e-4 && + WaterLily.L₂(u_udf[:,:,2].-0) < 1e-4 + ) + 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 make_accel_circle(radius=32,H=16) = Simulation(radius.*(2H,2H), - (i,t)-> i==1 ? t : zero(t), radius; U=1, mem=f, + (i,x,t)-> i==1 ? t : zero(t), radius; U=1, mem=f, body=AutoBody((x,t)->√sum(abs2,x .-H*radius)-radius)) sim = make_accel_circle(); sim_step!(sim) @test isapprox(WaterLily.pressure_force(sim)/(π*sim.L^2),[-1,0],atol=0.04) @@ -347,16 +424,16 @@ import WaterLily: × # stress tensor u₂ = zeros(N,N,2) |> f u₃ = zeros(N,N,N,3) |> f - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2),u₂) .≈ 0) - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ 0) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ 0) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ 0) apply!((i,x)->x[i],u₂) # uniform gradient apply!((i,x)->x[i],u₃) - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[2 0; 0 2]) - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ SA[2 0 0; 0 2 0; 0 0 2]) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[2 0; 0 2]) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ SA[2 0 0; 0 2 0; 0 0 2]) apply!((i,x)->x[i%2+1],u₂) # shear apply!((i,x)->x[i%3+1],u₃) - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[0 2; 2 0]) - @test GPUArrays.@allowscalar all(WaterLily.∇²u(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ SA[0 1 1; 1 0 1; 1 1 0]) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[0 2; 2 0]) + @test GPUArrays.@allowscalar all(2WaterLily.S(CartesianIndex(N÷2,N÷2,N÷2),u₃) .≈ SA[0 1 1; 1 0 1; 1 1 0]) # viscous force u₂ .= 0; u₃ .= 0 @test all(WaterLily.viscous_force(u₂,1.0,df₂,body) .≈ 0)