From 111c8f717d26a6ca550c529bd8e70fae7cf54239 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 6 Sep 2024 16:07:45 +0200 Subject: [PATCH 01/33] initial commit --- src/Flow.jl | 24 ++++++++++++------------ src/WaterLily.jl | 6 +++--- src/util.jl | 16 ++++++++++------ 3 files changed, 25 insertions(+), 21 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 9dbd39e6..a8d9b48e 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -65,11 +65,11 @@ using EllipsisNotation Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. """ -accelerate!(r,dt,g::Function,::Tuple,t=sum(dt)) = for i ∈ 1:last(size(r)) - r[..,i] .+= g(i,t) +accelerate!(r,dt,::Tuple,U::Function,t=sum(dt)) = for i ∈ 1:last(size(r)) + @loop r[I,i] += U(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(r) 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,g::Nothing,U::Function) = accelerate!(r,dt,(),(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t)) +accelerate!(r,dt,g::Function,U::Function) = accelerate!(r,dt,(),(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t)) accelerate!(r,dt,::Nothing,::Tuple) = nothing """ BCTuple(U,dt,N) @@ -112,7 +112,7 @@ struct Flow{D, T, Sf<:AbstractArray{T}, Vf<:AbstractArray{T}, Tf<:AbstractArray{ 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.) + BC!(u,U,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 @@ -153,18 +153,18 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp @fastmath function mom_step!(a::Flow{N},b::AbstractPoisson) where N a.u⁰ .= a.u; scale_u!(a,0) # predictor u → u' - U = BCTuple(a.U,@view(a.Δt[1:end-1]),N) + t = sum(@view(a.Δt[1:end-1])) 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) + BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t) + a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit + project!(a,b); BC!(a.u,a.U,a.exitBC,a.perdir,t) # corrector u → u¹ - U = BCTuple(a.U,a.Δt,N) + t = sum(a.Δt) 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) + BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) + project!(a,b,0.5); BC!(a.u,a.U,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)) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 811cea4b..0c501a31 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -39,7 +39,7 @@ 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)` + tuple `u_BC[i]=uᵢ, i=eachindex(dims)`, or a time and space-varying function `f(i,x,t)` - `L`: Simulation length scale. - `U`: Simulation velocity scale. - `Δt`: Initial time step. @@ -68,8 +68,8 @@ mutable struct Simulation 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λ + isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" + uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->u_BC(i,x,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) measure!(flow,body;ϵ) diff --git a/src/util.jl b/src/util.jl index d85f921f..4c72c1b3 100644 --- a/src/util.jl +++ b/src/util.jl @@ -166,7 +166,7 @@ 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=()) +function BC!(a,A,saveexit=false,perdir=(),t=0) N,n = size_u(a) for i ∈ 1:n, j ∈ 1:n if j in perdir @@ -175,9 +175,9 @@ 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,A,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,A,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) @@ -185,16 +185,20 @@ function BC!(a,A,saveexit=false,perdir=()) end end end +uBC(i,A,I,t) = A[i] +uBC(i,A::Function,I,t) = A(i,loc(i,I),t) + """ 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(@views(u[slice(N,2,1,2),1])) # inflow mass + @loop u[I,1] = u⁰[I,1]-U*Δt*(u⁰[I,1]-u⁰[I-δ(1,I),1]) over I ∈ exitR + ∮u = sum(@views(u[exitR,1]))/length(exitR)-U # mass flux imbalance @loop u[I,1] -= ∮u over I ∈ exitR # correct flux end """ From 81724f06f15305ad8822048fd544a7addf88af0c Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 6 Sep 2024 16:35:23 +0200 Subject: [PATCH 02/33] correct inflow mass flux --- src/util.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util.jl b/src/util.jl index 4c72c1b3..bb7cae38 100644 --- a/src/util.jl +++ b/src/util.jl @@ -196,9 +196,9 @@ Apply a 1D convection scheme to fill the ghost cell on the exit of the domain. function exitBC!(u,u⁰,Δt) N,_ = size_u(u) exitR = slice(N.-1,N[1],1,2) # exit slice excluding ghosts - U = sum(@views(u[slice(N,2,1,2),1])) # inflow mass + U = sum(@view(u[slice(N,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(@views(u[exitR,1]))/length(exitR)-U # mass flux imbalance + ∮u = sum(@view(u[exitR,1]))/length(exitR)-U # mass flux imbalance @loop u[I,1] -= ∮u over I ∈ exitR # correct flux end """ From 31086c227b1a956f213e3ad8c26f19cc8d54bc74 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Mon, 9 Sep 2024 17:57:09 +0200 Subject: [PATCH 03/33] fix accelerate! --- src/Flow.jl | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index a8d9b48e..c8c58350 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -70,14 +70,8 @@ accelerate!(r,dt,::Tuple,U::Function,t=sum(dt)) = for i ∈ 1:last(size(r)) end accelerate!(r,dt,g::Nothing,U::Function) = accelerate!(r,dt,(),(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t)) accelerate!(r,dt,g::Function,U::Function) = accelerate!(r,dt,(),(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t)) +accelerate!(r,dt,g::Function,::Tuple) = accelerate!(r,dt,(),(i,x,t)->g(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 """ Flow{D::Int, T::Float, Sf<:AbstractArray{T,D}, Vf<:AbstractArray{T,D+1}, Tf<:AbstractArray{T,D+2}} From e97d012a586e9ac51c776864d041b8921d4a86be Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 13 Sep 2024 10:54:38 +0200 Subject: [PATCH 04/33] fix test and accelerate! --- src/Flow.jl | 16 ++++++++-------- src/util.jl | 2 +- test/maintests.jl | 13 +++++++------ 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index c8c58350..3e41f875 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -65,13 +65,13 @@ using EllipsisNotation Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. """ -accelerate!(r,dt,::Tuple,U::Function,t=sum(dt)) = for i ∈ 1:last(size(r)) - @loop r[I,i] += U(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(r) +accelerate!(r,U::Function,t) = for i ∈ 1:last(size(r)) + @loop r[I,i] += U(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) end -accelerate!(r,dt,g::Nothing,U::Function) = accelerate!(r,dt,(),(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t)) -accelerate!(r,dt,g::Function,U::Function) = accelerate!(r,dt,(),(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t)) -accelerate!(r,dt,g::Function,::Tuple) = accelerate!(r,dt,(),(i,x,t)->g(i,t)) -accelerate!(r,dt,::Nothing,::Tuple) = nothing +accelerate!(r,t,g::Nothing,U::Function) = accelerate!(r,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t),t) +accelerate!(r,t,g::Function,U::Function) = accelerate!(r,(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t),t) +accelerate!(r,t,g::Function,::Tuple) = accelerate!(r,(i,x,t)->g(i,t),t) +accelerate!(r,t,::Nothing,::Tuple) = nothing """ Flow{D::Int, T::Float, Sf<:AbstractArray{T,D}, Vf<:AbstractArray{T,D+1}, Tf<:AbstractArray{T,D+2}} @@ -149,14 +149,14 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp # predictor u → u' t = sum(@view(a.Δt[1:end-1])) conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir) - accelerate!(a.f,@view(a.Δt[1:end-1]),a.g,a.U) + accelerate!(a.f,t,a.g,a.U) BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t) a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit project!(a,b); BC!(a.u,a.U,a.exitBC,a.perdir,t) # corrector u → u¹ t = sum(a.Δt) conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) - accelerate!(a.f,a.Δt,a.g,a.U) + accelerate!(a.f,t,a.g,a.U) BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) push!(a.Δt,CFL(a)) diff --git a/src/util.jl b/src/util.jl index bb7cae38..91b86f93 100644 --- a/src/util.jl +++ b/src/util.jl @@ -196,7 +196,7 @@ Apply a 1D convection scheme to fill the ghost cell on the exit of the domain. function exitBC!(u,u⁰,Δt) N,_ = size_u(u) exitR = slice(N.-1,N[1],1,2) # exit slice excluding ghosts - U = sum(@view(u[slice(N,2,1,2),1]))/length(exitR) # inflow mass flux + 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 diff --git a/test/maintests.jl b/test/maintests.jl index a32a0b54..cc961ab6 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -42,9 +42,13 @@ 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]) + A = [1.,2.,3]; Af(i,x,t) = i + @test all([WaterLily.uBC(i,A,CartesianIndex(1,1,1),0) for i in 1:3] .== A) + @test all([WaterLily.uBC(i,Af,CartesianIndex(1,1,1),0) for i in 1:3] .== A) + 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 @@ -154,9 +158,6 @@ 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 @@ -164,9 +165,9 @@ end @test all(a .== 0) WaterLily.accelerate!(a,[1],(i,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,t) -> i==1 ? t : 2*t,(i,x,t) -> i==1 ? -t : -2*t) @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) end # Impulsive flow in a box From 2766de4be0908ccbbeff3fd969a0ac759be66977 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 13 Sep 2024 11:32:31 +0200 Subject: [PATCH 05/33] fix maintest error --- test/maintests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/maintests.jl b/test/maintests.jl index cc961ab6..5d7eace5 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -161,13 +161,13 @@ end # 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,t) -> i==1 ? t : 2*t,()) @test all(a[:,:,1] .== 1) && all(a[:,:,2] .== 2) - WaterLily.accelerate!(a,[1],nothing,(i,x,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,x,t) -> i==1 ? -t : -2*t) + WaterLily.accelerate!(a,1,(i,t) -> i==1 ? t : 2*t,(i,x,t) -> i==1 ? -t : -2*t) @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) end # Impulsive flow in a box From d451f1cdfddf5f7c1c30b4e37d64921f2d58ede5 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 27 Feb 2025 15:04:00 +0100 Subject: [PATCH 06/33] fix allocation issues, merge body_force! and update test --- src/Flow.jl | 20 ++++++++++++++++---- src/WaterLily.jl | 8 ++++---- src/util.jl | 9 ++++----- test/maintests.jl | 27 ++++++++++++++++++++++++--- 4 files changed, 48 insertions(+), 16 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 3e41f875..b30b8e11 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -61,18 +61,28 @@ upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I using EllipsisNotation """ - accelerate!(r,dt,g) + accelerate!(r,t,g) Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. """ -accelerate!(r,U::Function,t) = for i ∈ 1:last(size(r)) - @loop r[I,i] += U(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) +accelerate!(r,g::Function,t) = for i ∈ 1:last(size(r)) + @loop r[I,i] += g(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) end accelerate!(r,t,g::Nothing,U::Function) = accelerate!(r,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t),t) accelerate!(r,t,g::Function,U::Function) = accelerate!(r,(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t),t) accelerate!(r,t,g::Function,::Tuple) = accelerate!(r,(i,x,t)->g(i,t),t) accelerate!(r,t,::Nothing,::Tuple) = nothing +""" + body_force!(r,force,t) + +Adds a body force to the RHS +- `force` is either a vector field or a function `f(i,x,t)` that returns the component + `i` of the vector field at time `t` and position `x`. +""" +body_force!(r,::Nothing,t=0) = nothing +body_force!(r,force::AbstractArray,t=0) = r .+= force +body_force!(r,force::Function,t) = accelerate!(r,force,t) """ Flow{D::Int, T::Float, Sf<:AbstractArray{T,D}, Vf<:AbstractArray{T,D+1}, Tf<:AbstractArray{T,D+2}} @@ -144,11 +154,12 @@ 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 +@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson;body_force=nothing) where N a.u⁰ .= a.u; scale_u!(a,0) # predictor u → u' t = sum(@view(a.Δt[1:end-1])) conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir) + body_force!(a.f,body_force,t) accelerate!(a.f,t,a.g,a.U) BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t) a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit @@ -156,6 +167,7 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp # corrector u → u¹ t = sum(a.Δt) conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) + body_force!(a.f,body_force,t) accelerate!(a.f,t,a.g,a.U) BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 0c501a31..b8a26084 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -94,17 +94,17 @@ 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. """ -function sim_step!(sim::Simulation,t_end;remeasure=true,max_steps=typemax(Int),verbose=false) +function sim_step!(sim::Simulation,t_end;remeasure=true,max_steps=typemax(Int),body_force=nothing,verbose=false) 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, body_force) verbose && println("tU/L=",round(sim_time(sim),digits=4), ", Δt=",round(sim.flow.Δt[end],digits=3)) end end -function sim_step!(sim::Simulation;remeasure=true) +function sim_step!(sim::Simulation;remeasure=true,body_force=nothing) remeasure && measure!(sim) - mom_step!(sim.flow,sim.pois) + mom_step!(sim.flow, sim.pois; body_force) end """ diff --git a/src/util.jl b/src/util.jl index 91b86f93..b0f01d48 100644 --- a/src/util.jl +++ b/src/util.jl @@ -166,7 +166,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=(),t=0) +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 @@ -175,9 +176,9 @@ function BC!(a,A,saveexit=false,perdir=(),t=0) else if i==j # Normal direction, Dirichlet for s ∈ (1,2) - @loop a[I,i] = uBC(i,A,I,t) 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] = uBC(i,A,I,t) 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) @@ -185,8 +186,6 @@ function BC!(a,A,saveexit=false,perdir=(),t=0) end end end -uBC(i,A,I,t) = A[i] -uBC(i,A::Function,I,t) = A(i,loc(i,I),t) """ exitBC!(u,u⁰,U,Δt) diff --git a/test/maintests.jl b/test/maintests.jl index 5d7eace5..76bfa98b 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -45,9 +45,16 @@ using ReadVTK, WriteVTK WaterLily.exitBC!(u,u,0) # conservative exit check @test GPUArrays.@allowscalar all(u[end,2:end-1, 1] .== U[1]) - A = [1.,2.,3]; Af(i,x,t) = i - @test all([WaterLily.uBC(i,A,CartesianIndex(1,1,1),0) for i in 1:3] .== A) - @test all([WaterLily.uBC(i,Af,CartesianIndex(1,1,1),0) for i in 1:3] .== A) + # 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 all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) + @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) + # test 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]) @@ -169,6 +176,20 @@ end @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) WaterLily.accelerate!(a,1,(i,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 + b = zeros(N,N,2) |> f; bf = ones(N,N,2) |> f + WaterLily.body_force!(b,nothing) + @test all(b .== 0) + WaterLily.body_force!(b, bf) + @test all(b .== 1) + apply!((i,x)->x[i], bf) + WaterLily.body_force!(b, bf) + a .= 0 # reset + WaterLily.accelerate!(a,(i,x,t)->x[i]+1.0,1.) + @test all(b .== a) + WaterLily.body_force!(b,(i,x,t)->x[i]+1.0,1.) + WaterLily.accelerate!(a,(i,x,t)->x[i]+1.0,1.) + @test all(b .== a) end # Impulsive flow in a box U = (2/3, -1/3) From 633b6d99ed537aafd0973d2b57a40df7349eb02d Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 27 Feb 2025 15:48:06 +0100 Subject: [PATCH 07/33] tidy interface and update test --- src/Flow.jl | 15 +++++++-------- test/maintests.jl | 8 ++++---- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index b30b8e11..bd348104 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -61,16 +61,13 @@ upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I using EllipsisNotation """ - accelerate!(r,t,g) + accelerate!(r,t,g,U) Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. """ -accelerate!(r,g::Function,t) = for i ∈ 1:last(size(r)) - @loop r[I,i] += g(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) -end -accelerate!(r,t,g::Nothing,U::Function) = accelerate!(r,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t),t) -accelerate!(r,t,g::Function,U::Function) = accelerate!(r,(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t),t) -accelerate!(r,t,g::Function,::Tuple) = accelerate!(r,(i,x,t)->g(i,t),t) +accelerate!(r,t,g::Nothing,U::Function) = body_force!(r,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t),t) +accelerate!(r,t,g::Function,U::Function) = body_force!(r,(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t),t) +accelerate!(r,t,g::Function,::Tuple) = body_force!(r,(i,x,t)->g(i,t),t) accelerate!(r,t,::Nothing,::Tuple) = nothing """ body_force!(r,force,t) @@ -82,7 +79,9 @@ Adds a body force to the RHS """ body_force!(r,::Nothing,t=0) = nothing body_force!(r,force::AbstractArray,t=0) = r .+= force -body_force!(r,force::Function,t) = accelerate!(r,force,t) +body_force!(r,force::Function,t) = for i ∈ 1:last(size(r)) + @loop r[I,i] += force(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) +end """ Flow{D::Int, T::Float, Sf<:AbstractArray{T,D}, Vf<:AbstractArray{T,D+1}, Tf<:AbstractArray{T,D+2}} diff --git a/test/maintests.jl b/test/maintests.jl index 76bfa98b..363b566d 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -178,17 +178,17 @@ end @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) # check applying body force b = zeros(N,N,2) |> f; bf = ones(N,N,2) |> f - WaterLily.body_force!(b,nothing) + WaterLily.body_force!(b, nothing) @test all(b .== 0) WaterLily.body_force!(b, bf) @test all(b .== 1) apply!((i,x)->x[i], bf) WaterLily.body_force!(b, bf) - a .= 0 # reset - WaterLily.accelerate!(a,(i,x,t)->x[i]+1.0,1.) + a .= 0 # reset and accelerate using a non-uniform velocity field + WaterLily.accelerate!(a,1.,nothing,(i,x,t)->t*(x[i]+1.0)) @test all(b .== a) WaterLily.body_force!(b,(i,x,t)->x[i]+1.0,1.) - WaterLily.accelerate!(a,(i,x,t)->x[i]+1.0,1.) + WaterLily.accelerate!(a,1.,nothing,(i,x,t)->t*(x[i]+1.0)) @test all(b .== a) end # Impulsive flow in a box From ab62bbc5e645fd5429befab19ae20c774056b69e Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 28 Feb 2025 16:07:51 +0100 Subject: [PATCH 08/33] try solve compatibility issue --- src/WaterLily.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index d3d40a15..9aa7d273 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,16 +69,15 @@ mutable struct Simulation <: AbstractSimulation 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" - uλ = (isnothing(uλ) && !isa(u_BC,Function)) ? (i,x)->u_BC[i] : uλ - if hasmethod(u_BC, Tuple{Int,Number}) # uniform case - @assert all(typeof.(ntuple(i->u_BC(i,zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = (i,x)->u_BC(i,zero(T)) - elseif hasmethod(u_BC, Tuple{Int,SVector,Number}) # non-uniform case - @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = (i,x)->u_BC(i,x,zero(T)) - end + hasmethod(u_BC,Tuple{Int,Number}) && (uBC(i,x,t) = u_BC(i,t)) + hasmethod(u_BC,Tuple{Int,SVector,Number}) && (uBC(i,x,t) = u_BC(i,x,t)) + # u_BC = hasmethod(u_BC,Tuple{Int,Number}) ? (tmp=u_BC; (i,x,t)->tmp(i,t)) : u_BC # map to non-uniform case + # @show hasmethod(u_BC,Tuple{Int,Number}) + # @show hasmethod(u_BC,Tuple{Int,SVector,Number}) + isa(u_BC,Function) && @assert all(typeof.(ntuple(i->uBC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" + uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->uBC(i,x,zero(T)),(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) + 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 From 05336b550c04831d06f333b2ad58a9d5d150c412 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Mon, 3 Mar 2025 12:52:22 +0100 Subject: [PATCH 09/33] fix Simulation constructor --- src/WaterLily.jl | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 9aa7d273..ece031f3 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,14 +69,10 @@ mutable struct Simulation <: AbstractSimulation 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" - hasmethod(u_BC,Tuple{Int,Number}) && (uBC(i,x,t) = u_BC(i,t)) - hasmethod(u_BC,Tuple{Int,SVector,Number}) && (uBC(i,x,t) = u_BC(i,x,t)) - # u_BC = hasmethod(u_BC,Tuple{Int,Number}) ? (tmp=u_BC; (i,x,t)->tmp(i,t)) : u_BC # map to non-uniform case - # @show hasmethod(u_BC,Tuple{Int,Number}) - # @show hasmethod(u_BC,Tuple{Int,SVector,Number}) - isa(u_BC,Function) && @assert all(typeof.(ntuple(i->uBC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->uBC(i,x,zero(T)),(i,x)->u_BC[i]) : uλ - U = isnothing(U) ? √sum(abs2,u_BC) : U # default if not specified + hasmethod(u_BC,Tuple{Int,Number}) ? (uBC(i,x,t) = u_BC(i,t)) : uBC = u_BC # copy the tuple or the function U(i,x,t) + isa(uBC,Function) && @assert all(typeof.(ntuple(i->uBC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" + uλ = isnothing(uλ) ? ifelse(isa(uBC,Function),(i,x)->uBC(i,x,zero(T)),(i,x)->uBC[i]) : uλ + U = isnothing(U) ? √sum(abs2,uBC) : U # default if not specified 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)) From 3653bbf1f50a2a38360f0fda9e941769b119b6ce Mon Sep 17 00:00:00 2001 From: marinlauber Date: Tue, 4 Mar 2025 12:07:09 +0100 Subject: [PATCH 10/33] fix BC time and constructor --- src/Flow.jl | 19 +++++++++---------- src/WaterLily.jl | 10 +++++----- test/maintests.jl | 2 +- 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index fc330bf8..19403eb6 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -154,23 +154,22 @@ 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{N},b::AbstractPoisson;body_force=nothing) where N - a.u⁰ .= a.u; scale_u!(a,0); t = sum(@view(a.Δt[1:end-1])) + 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) - body_force!(a.f,body_force,t) - accelerate!(a.f,t,a.g,a.U) - BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t) + body_force!(a.f,body_force,t₀) + accelerate!(a.f,t₀,a.g,a.U) + BDIM!(a); BC!(a.u,a.U,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.U,a.exitBC,a.perdir,t) + project!(a,b); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) # corrector u → u¹ @log "c" - t = sum(a.Δt) conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) - body_force!(a.f,body_force,t) - accelerate!(a.f,t,a.g,a.U) - BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) - project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t) + body_force!(a.f,body_force,t₁) + accelerate!(a.f,t₁,a.g,a.U) + BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) + project!(a,b,0.5); BC!(a.u,a.U,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)) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index ece031f3..da938957 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,11 +69,11 @@ mutable struct Simulation <: AbstractSimulation 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" - hasmethod(u_BC,Tuple{Int,Number}) ? (uBC(i,x,t) = u_BC(i,t)) : uBC = u_BC # copy the tuple or the function U(i,x,t) - isa(uBC,Function) && @assert all(typeof.(ntuple(i->uBC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = isnothing(uλ) ? ifelse(isa(uBC,Function),(i,x)->uBC(i,x,zero(T)),(i,x)->uBC[i]) : uλ - U = isnothing(U) ? √sum(abs2,uBC) : U # default if not specified - flow = Flow(dims,uBC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) + isa(u_BC,Function) && @assert hasmethod(u_BC,Tuple{Int,SVector,Number}) "`u_BC` must be a function of (i,x,t)" + isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" + uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->u_BC(i,x,zero(T)),(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) measure!(flow,body;ϵ) new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) end diff --git a/test/maintests.jl b/test/maintests.jl index 3c0d0188..0aa57304 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -326,7 +326,7 @@ 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) From 5226d1aa46be7ef01c38ed26d0a5061138bfebb7 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 4 Mar 2025 13:31:00 +0100 Subject: [PATCH 11/33] =?UTF-8?q?Small=20adjustements=20on=20how=20to=20cr?= =?UTF-8?q?eate=20u=5FBC=20and=20u=CE=BB=20in=20Simulation=20constructor.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/WaterLily.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index da938957..3022d21e 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,10 +69,10 @@ mutable struct Simulation <: AbstractSimulation 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 hasmethod(u_BC,Tuple{Int,SVector,Number}) "`u_BC` must be a function of (i,x,t)" + isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC as a function need to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" - uλ = isnothing(uλ) ? ifelse(isa(u_BC,Function),(i,x)->u_BC(i,x,zero(T)),(i,x)->u_BC[i]) : uλ - U = isnothing(U) ? √sum(abs2,u_BC) : U # default if not specified + isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,zero(T)) : uλ = (i,_)->u_BC[i]) + isnothing(U) && (U = √sum(abs2,u_BC)) flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) From 247c1a826785bc252eb368b30a503ccdc8f53944 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 4 Mar 2025 13:48:42 +0100 Subject: [PATCH 12/33] =?UTF-8?q?Reverted=20from=20zero(T)=20to=200=20in?= =?UTF-8?q?=20defining=20u=CE=BB=20cause=20of=20GPU=20nonbits=20error?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/WaterLily.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 3022d21e..ef9f524d 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -71,7 +71,7 @@ mutable struct Simulation <: AbstractSimulation @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC as a function need to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" - isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,zero(T)) : uλ = (i,_)->u_BC[i]) + isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) isnothing(U) && (U = √sum(abs2,u_BC)) flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) From f4f4697e8d6bc4d3ad074484cadd754d8ba8033e Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 4 Mar 2025 14:00:37 +0100 Subject: [PATCH 13/33] u_BC everywhere (instead of mixed with uBC) and same default Float32 type for Simulation and Flow (Flow was Float64 previously). --- src/Flow.jl | 4 ++-- src/util.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 19403eb6..4c3d8819 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -111,7 +111,7 @@ struct Flow{D, T, Sf<:AbstractArray{T}, Vf<:AbstractArray{T}, Tf<:AbstractArray{ 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 + uλ::Function=(i, x) -> 0., perdir=(), exitBC=false, T=Float32) where D Ng = N .+ 2 Nd = (Ng..., D) u = Array{T}(undef, Nd...) |> f; apply!(uλ, u); @@ -154,7 +154,7 @@ 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{N},b::AbstractPoisson;body_force=nothing) where N - a.u⁰ .= a.u; scale_u!(a,0); t₁ = sum(a.Δt); t₀ = t₁-a.Δt[end] + 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) diff --git a/src/util.jl b/src/util.jl index e669459a..6a55698b 100644 --- a/src/util.jl +++ b/src/util.jl @@ -190,7 +190,7 @@ boundary. For example `aₓ(x)=Aₓ ∀ x ∈ minmax(X)`. A zero Neumann conditi is applied to the tangential components. """ 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) +function BC!(a,u_BC::Function,saveexit=false,perdir=(),t=0) N,n = size_u(a) for i ∈ 1:n, j ∈ 1:n if j in perdir @@ -199,9 +199,9 @@ function BC!(a,uBC::Function,saveexit=false,perdir=(),t=0) else if i==j # Normal direction, Dirichlet for s ∈ (1,2) - @loop a[I,i] = uBC(i,loc(i,I),t) over I ∈ slice(N,s,j) + @loop a[I,i] = u_BC(i,loc(i,I),t) over I ∈ slice(N,s,j) end - (!saveexit || i>1) && (@loop a[I,i] = uBC(i,loc(i,I),t) over I ∈ slice(N,N[j],j)) # overwrite exit + (!saveexit || i>1) && (@loop a[I,i] = u_BC(i,loc(i,I),t) over I ∈ slice(N,N[j],j)) # overwrite exit else # Tangential directions, Neumann @loop a[I,i] = a[I+δ(j,I),i] over I ∈ slice(N,1,j) @loop a[I,i] = a[I-δ(j,I),i] over I ∈ slice(N,N[j],j) From be9e85f8a6a420904c7c9e8dfddf34e822588ee7 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 4 Mar 2025 16:58:29 +0100 Subject: [PATCH 14/33] Added a user defined function (udf) feature to pass functions into sim_step with arbitrary keyword arguments. The function can modify the ::Flow object and is called both in the predictor and corrector steps. Added tests for its implementation of an increasing body force, but using udf instead of body_force. All tests passing on Array and CuArray. --- src/Flow.jl | 19 ++++++++++++++----- src/WaterLily.jl | 12 +++++++----- test/maintests.jl | 17 +++++++++++++++-- 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 4c3d8819..3448420c 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -77,8 +77,8 @@ Adds a body force to the RHS - `force` is either a vector field or a function `f(i,x,t)` that returns the component `i` of the vector field at time `t` and position `x`. """ -body_force!(r,::Nothing,t=0) = nothing -body_force!(r,force::AbstractArray,t=0) = r .+= force +body_force!(_,::Nothing,_) = nothing +body_force!(r,force::AbstractArray,_) = r .+= force body_force!(r,force::Function,t) = for i ∈ 1:last(size(r)) @loop r[I,i] += force(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) end @@ -153,12 +153,12 @@ 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;body_force=nothing) where N +@fastmath function mom_step!(a::Flow{N},b::AbstractPoisson;body_force=nothing,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) - body_force!(a.f,body_force,t₀) + body_force!(a.f,body_force,t₀); udf!(a,udf,t₀; kwargs...) accelerate!(a.f,t₀,a.g,a.U) BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) # BC MUST be at t₁ a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit @@ -166,7 +166,7 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp # corrector u → u¹ @log "c" conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) - body_force!(a.f,body_force,t₁) + body_force!(a.f,body_force,t₁); udf!(a,udf,t₁; kwargs...) accelerate!(a.f,t₁,a.g,a.U) BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) @@ -185,3 +185,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!(_,::Nothing,_) = nothing +udf!(flow::Flow,force!::Function,t; kwargs...) = force!(flow,t; kwargs...) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index ef9f524d..307af650 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,7 +69,7 @@ mutable struct Simulation <: AbstractSimulation 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 first(methods(u_BC)).nargs==4 "u_BC as a function need to be defined as u_BC(i,x,t)" + isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) isnothing(U) && (U = √sum(abs2,u_BC)) @@ -95,18 +95,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),body_force=nothing,verbose=false) +function sim_step!(sim::AbstractSimulation,t_end;remeasure=true,max_steps=typemax(Int),body_force=nothing,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, body_force) + sim_step!(sim; remeasure, body_force, 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,body_force=nothing) +function sim_step!(sim::AbstractSimulation;remeasure=true,body_force=nothing,udf=nothing,kwargs...) remeasure && measure!(sim) - mom_step!(sim.flow, sim.pois; body_force) + mom_step!(sim.flow, sim.pois; body_force, udf, kwargs...) end """ diff --git a/test/maintests.jl b/test/maintests.jl index 0aa57304..e22405d9 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -299,20 +299,24 @@ 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 + !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 @@ -320,6 +324,15 @@ end 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.U[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 From 24aa6b71207589df9f685e883e78a68e09da59a8 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 4 Mar 2025 17:13:16 +0100 Subject: [PATCH 15/33] Fixed body_force tests. Remember to always runs all the tests... --- test/maintests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/maintests.jl b/test/maintests.jl index e22405d9..0d26ddf7 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -178,12 +178,12 @@ end @test all(a[:,:,1] .== 0) && all(a[:,:,2] .== 0) # check applying body force b = zeros(N,N,2) |> f; bf = ones(N,N,2) |> f - WaterLily.body_force!(b, nothing) + WaterLily.body_force!(b, nothing, 0) @test all(b .== 0) - WaterLily.body_force!(b, bf) + WaterLily.body_force!(b, bf, 0) @test all(b .== 1) apply!((i,x)->x[i], bf) - WaterLily.body_force!(b, bf) + WaterLily.body_force!(b, bf, 0) a .= 0 # reset and accelerate using a non-uniform velocity field WaterLily.accelerate!(a,1.,nothing,(i,x,t)->t*(x[i]+1.0)) @test all(b .== a) From a9326397b036368ad02af54cfdd59b02ae63c658 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 5 Mar 2025 13:44:55 +0100 Subject: [PATCH 16/33] removes body_force and defines acceleration g as a function of g(i,x,t), hence implementing the space-varying featured previously in charge of body_force. Tests have been mdofied accordingly, an all are locally passing for CPU and GPU. Need to update examples repo. --- src/Flow.jl | 30 ++++++++++-------------------- src/WaterLily.jl | 8 ++++---- test/maintests.jl | 21 +++++++++------------ 3 files changed, 23 insertions(+), 36 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 3448420c..483aca36 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -63,25 +63,15 @@ using EllipsisNotation """ accelerate!(r,t,g,U) -Add a uniform acceleration `gᵢ+dUᵢ/dt` at time `t=sum(dt)` to field `r`. +Adds a space-time acceleration field g(i,x,t) and velocity field U(i,x,t) such that `rᵢ += gᵢ(i,x,t)+dUᵢ(i,x,t)/dt` at time `t=sum(dt)` to vector field `r`. """ -accelerate!(r,t,g::Nothing,U::Function) = body_force!(r,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t),t) -accelerate!(r,t,g::Function,U::Function) = body_force!(r,(i,x,t)->g(i,t)+ForwardDiff.derivative(τ->U(i,x,τ),t),t) -accelerate!(r,t,g::Function,::Tuple) = body_force!(r,(i,x,t)->g(i,t),t) -accelerate!(r,t,::Nothing,::Tuple) = nothing -""" - body_force!(r,force,t) - -Adds a body force to the RHS - -- `force` is either a vector field or a function `f(i,x,t)` that returns the component - `i` of the vector field at time `t` and position `x`. -""" -body_force!(_,::Nothing,_) = nothing -body_force!(r,force::AbstractArray,_) = r .+= force -body_force!(r,force::Function,t) = for i ∈ 1:last(size(r)) - @loop r[I,i] += force(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) +accelerate!(_,_,::Nothing,::Union{Nothing,Tuple}) = nothing +accelerate!(r,t,f) = for i ∈ 1:last(size(r)) + @loop r[I,i] += f(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) end +accelerate!(r,t,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}} @@ -153,12 +143,12 @@ 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;body_force=nothing,udf=nothing,kwargs...) where 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) - body_force!(a.f,body_force,t₀); udf!(a,udf,t₀; kwargs...) + udf!(a,udf,t₀; kwargs...) accelerate!(a.f,t₀,a.g,a.U) BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) # BC MUST be at t₁ a.exitBC && exitBC!(a.u,a.u⁰,a.Δt[end]) # convective exit @@ -166,7 +156,7 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp # corrector u → u¹ @log "c" conv_diff!(a.f,a.u,a.σ,ν=a.ν,perdir=a.perdir) - body_force!(a.f,body_force,t₁); udf!(a,udf,t₁; kwargs...) + udf!(a,udf,t₁; kwargs...) accelerate!(a.f,t₁,a.g,a.U) BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 307af650..bf3abd4f 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -98,17 +98,17 @@ 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),body_force=nothing,udf=nothing,verbose=false,kwargs...) +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, body_force, udf, kwargs...) + 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,body_force=nothing,udf=nothing,kwargs...) +function sim_step!(sim::AbstractSimulation;remeasure=true,udf=nothing,kwargs...) remeasure && measure!(sim) - mom_step!(sim.flow, sim.pois; body_force, udf, kwargs...) + mom_step!(sim.flow, sim.pois; udf, kwargs...) end """ diff --git a/test/maintests.jl b/test/maintests.jl index 0d26ddf7..dec628a8 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -170,24 +170,21 @@ end N = 4; a = zeros(N,N,2) |> f 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,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,x,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 - b = zeros(N,N,2) |> f; bf = ones(N,N,2) |> f - WaterLily.body_force!(b, nothing, 0) - @test all(b .== 0) - WaterLily.body_force!(b, bf, 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) - apply!((i,x)->x[i], bf) - WaterLily.body_force!(b, bf, 0) a .= 0 # reset and accelerate using a non-uniform velocity field - WaterLily.accelerate!(a,1.,nothing,(i,x,t)->t*(x[i]+1.0)) + 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.body_force!(b,(i,x,t)->x[i]+1.0,1.) + 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 @@ -304,7 +301,7 @@ function acceleratingFlow(N;use_g=false,T=Float64,perdir=(1,),jerk=4,mem=Array) # 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 From 404ddabdb2ebcec3fe2d3c550101bf143acda03f Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 5 Mar 2025 14:07:54 +0100 Subject: [PATCH 17/33] Added test for when both g and U are present. --- test/maintests.jl | 260 +++++++++++++++++++++++----------------------- 1 file changed, 131 insertions(+), 129 deletions(-) diff --git a/test/maintests.jl b/test/maintests.jl index dec628a8..a65d969a 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -2,130 +2,130 @@ using GPUArrays using ReadVTK, WriteVTK @info "Test backends: $(join(arrays,", "))" -@testset "util.jl" begin - I = CartesianIndex(1,2,3,4) - @test I+δ(3,I) == CartesianIndex(1,2,4,4) - @test WaterLily.CI(I,5)==CartesianIndex(1,2,3,4,5) - @test WaterLily.CIj(3,I,5)==CartesianIndex(1,2,5,4) - @test WaterLily.CIj(2,CartesianIndex(16,16,16,3),14)==CartesianIndex(16,14,16,3) - - @test loc(3,CartesianIndex(3,4,5)) == SVector(3,4,4.5) .- 1.5 - I = CartesianIndex(rand(2:10,3)...) - @test loc(0,I) == SVector(I.I...) .- 1.5 - - ex,sym = :(a[I,i] = Math.add(p.b[I],func(I,q))),[] - WaterLily.grab!(sym,ex) - @test ex == :(a[I, i] = Math.add(b[I], func(I, q))) - @test sym == [:a, :I, :i, :(p.b), :q] - - for f ∈ arrays - p = zeros(4,5) |> f - apply!(x->x[1]+x[2]+3,p) # add 2×1.5 to move edge to origin - @test inside(p) == CartesianIndices((2:3,2:4)) - @test inside(p,buff=0) == CartesianIndices(p) - @test L₂(p) == 187 - - u = zeros(5,5,2) |> f - apply!((i,x)->x[i],u) - @test GPUArrays.@allowscalar [u[i,j,1].-(i-2) for i in 1:3, j in 1:3]==zeros(3,3) - - Ng, D, U = (6, 6), 2, (1.0, 0.5) - u = rand(Ng..., D) |> f # vector - σ = rand(Ng...) |> f # scalar - BC!(u, U) - @test GPUArrays.@allowscalar 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 GPUArrays.@allowscalar 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]) - - GPUArrays.@allowscalar u[end,:,1] .= 3 - BC!(u,U,true) # save exit values - @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3) - - 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 all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) - @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) - # test 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 - @test GPUArrays.@allowscalar all(σ[1, 2:end-1] .== σ[end-1, 2:end-1]) && all(σ[2:end-1, 1] .== σ[2:end-1, end-1]) - - u = rand(Ng..., D) |> f # vector - 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 interpolation - a = zeros(5,5,2) |> f; b = zeros(5,5) |> f - apply!((i,x)->x[i]+1.5,a); apply!(x->x[1]+1.5,b) # offset for start of grid - @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(2.5,1),a) .≈ [2.5,1.]) - @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(3.5,3),a) .≈ [3.5,3.]) - @test GPUArrays.@allowscalar WaterLily.interp(SVector(2.5,1),b) ≈ 2.5 - @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5 - end -end - -function Poisson_setup(poisson,N::NTuple{D};f=Array,T=Float32) where D - c = ones(T,N...,D) |> f; BC!(c, ntuple(zero,D)) - x = zeros(T,N) |> f; z = copy(x) - pois = poisson(x,c,z) - soln = map(I->T(I.I[1]),CartesianIndices(N)) |> f - I = first(inside(x)) - GPUArrays.@allowscalar @. soln -= soln[I] - z = mult!(pois,soln) - solver!(pois) - GPUArrays.@allowscalar @. x -= x[I] - return L₂(x-soln)/L₂(soln),pois -end - -@testset "Poisson.jl" begin - for f ∈ arrays - err,pois = Poisson_setup(Poisson,(5,5);f) - @test GPUArrays.@allowscalar parent(pois.D)==f(Float32[0 0 0 0 0; 0 -2 -3 -2 0; 0 -3 -4 -3 0; 0 -2 -3 -2 0; 0 0 0 0 0]) - @test GPUArrays.@allowscalar parent(pois.iD)≈f(Float32[0 0 0 0 0; 0 -1/2 -1/3 -1/2 0; 0 -1/3 -1/4 -1/3 0; 0 -1/2 -1/3 -1/2 0; 0 0 0 0 0]) - @test err < 1e-5 - err,pois = Poisson_setup(Poisson,(2^6+2,2^6+2);f) - @test err < 1e-6 - @test pois.n[] < 310 - err,pois = Poisson_setup(Poisson,(2^4+2,2^4+2,2^4+2);f) - @test err < 1e-6 - @test pois.n[] < 35 - end -end - -@testset "MultiLevelPoisson.jl" begin - I = CartesianIndex(4,3,2) - @test all(WaterLily.down(J)==I for J ∈ WaterLily.up(I)) - @test_throws AssertionError("MultiLevelPoisson requires size=a2ⁿ, where n>2") Poisson_setup(MultiLevelPoisson,(15+2,3^4+2)) - - err,pois = Poisson_setup(MultiLevelPoisson,(10,10)) - @test pois.levels[3].D == Float32[0 0 0 0; 0 -2 -2 0; 0 -2 -2 0; 0 0 0 0] - @test err < 1e-5 - - pois.levels[1].L[5:6,:,1].=0 - WaterLily.update!(pois) - @test pois.levels[3].D == Float32[0 0 0 0; 0 -1 -1 0; 0 -1 -1 0; 0 0 0 0] - - for f ∈ arrays - err,pois = Poisson_setup(MultiLevelPoisson,(2^6+2,2^6+2);f) - @test err < 1e-6 - @test pois.n[] ≤ 3 - err,pois = Poisson_setup(MultiLevelPoisson,(2^4+2,2^4+2,2^4+2);f) - @test err < 1e-6 - @test pois.n[] ≤ 3 - end -end +# @testset "util.jl" begin +# I = CartesianIndex(1,2,3,4) +# @test I+δ(3,I) == CartesianIndex(1,2,4,4) +# @test WaterLily.CI(I,5)==CartesianIndex(1,2,3,4,5) +# @test WaterLily.CIj(3,I,5)==CartesianIndex(1,2,5,4) +# @test WaterLily.CIj(2,CartesianIndex(16,16,16,3),14)==CartesianIndex(16,14,16,3) + +# @test loc(3,CartesianIndex(3,4,5)) == SVector(3,4,4.5) .- 1.5 +# I = CartesianIndex(rand(2:10,3)...) +# @test loc(0,I) == SVector(I.I...) .- 1.5 + +# ex,sym = :(a[I,i] = Math.add(p.b[I],func(I,q))),[] +# WaterLily.grab!(sym,ex) +# @test ex == :(a[I, i] = Math.add(b[I], func(I, q))) +# @test sym == [:a, :I, :i, :(p.b), :q] + +# for f ∈ arrays +# p = zeros(4,5) |> f +# apply!(x->x[1]+x[2]+3,p) # add 2×1.5 to move edge to origin +# @test inside(p) == CartesianIndices((2:3,2:4)) +# @test inside(p,buff=0) == CartesianIndices(p) +# @test L₂(p) == 187 + +# u = zeros(5,5,2) |> f +# apply!((i,x)->x[i],u) +# @test GPUArrays.@allowscalar [u[i,j,1].-(i-2) for i in 1:3, j in 1:3]==zeros(3,3) + +# Ng, D, U = (6, 6), 2, (1.0, 0.5) +# u = rand(Ng..., D) |> f # vector +# σ = rand(Ng...) |> f # scalar +# BC!(u, U) +# @test GPUArrays.@allowscalar 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 GPUArrays.@allowscalar 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]) + +# GPUArrays.@allowscalar u[end,:,1] .= 3 +# BC!(u,U,true) # save exit values +# @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3) + +# 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 all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) +# @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) +# # test 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 +# @test GPUArrays.@allowscalar all(σ[1, 2:end-1] .== σ[end-1, 2:end-1]) && all(σ[2:end-1, 1] .== σ[2:end-1, end-1]) + +# u = rand(Ng..., D) |> f # vector +# 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 interpolation +# a = zeros(5,5,2) |> f; b = zeros(5,5) |> f +# apply!((i,x)->x[i]+1.5,a); apply!(x->x[1]+1.5,b) # offset for start of grid +# @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(2.5,1),a) .≈ [2.5,1.]) +# @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(3.5,3),a) .≈ [3.5,3.]) +# @test GPUArrays.@allowscalar WaterLily.interp(SVector(2.5,1),b) ≈ 2.5 +# @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5 +# end +# end + +# function Poisson_setup(poisson,N::NTuple{D};f=Array,T=Float32) where D +# c = ones(T,N...,D) |> f; BC!(c, ntuple(zero,D)) +# x = zeros(T,N) |> f; z = copy(x) +# pois = poisson(x,c,z) +# soln = map(I->T(I.I[1]),CartesianIndices(N)) |> f +# I = first(inside(x)) +# GPUArrays.@allowscalar @. soln -= soln[I] +# z = mult!(pois,soln) +# solver!(pois) +# GPUArrays.@allowscalar @. x -= x[I] +# return L₂(x-soln)/L₂(soln),pois +# end + +# @testset "Poisson.jl" begin +# for f ∈ arrays +# err,pois = Poisson_setup(Poisson,(5,5);f) +# @test GPUArrays.@allowscalar parent(pois.D)==f(Float32[0 0 0 0 0; 0 -2 -3 -2 0; 0 -3 -4 -3 0; 0 -2 -3 -2 0; 0 0 0 0 0]) +# @test GPUArrays.@allowscalar parent(pois.iD)≈f(Float32[0 0 0 0 0; 0 -1/2 -1/3 -1/2 0; 0 -1/3 -1/4 -1/3 0; 0 -1/2 -1/3 -1/2 0; 0 0 0 0 0]) +# @test err < 1e-5 +# err,pois = Poisson_setup(Poisson,(2^6+2,2^6+2);f) +# @test err < 1e-6 +# @test pois.n[] < 310 +# err,pois = Poisson_setup(Poisson,(2^4+2,2^4+2,2^4+2);f) +# @test err < 1e-6 +# @test pois.n[] < 35 +# end +# end + +# @testset "MultiLevelPoisson.jl" begin +# I = CartesianIndex(4,3,2) +# @test all(WaterLily.down(J)==I for J ∈ WaterLily.up(I)) +# @test_throws AssertionError("MultiLevelPoisson requires size=a2ⁿ, where n>2") Poisson_setup(MultiLevelPoisson,(15+2,3^4+2)) + +# err,pois = Poisson_setup(MultiLevelPoisson,(10,10)) +# @test pois.levels[3].D == Float32[0 0 0 0; 0 -2 -2 0; 0 -2 -2 0; 0 0 0 0] +# @test err < 1e-5 + +# pois.levels[1].L[5:6,:,1].=0 +# WaterLily.update!(pois) +# @test pois.levels[3].D == Float32[0 0 0 0; 0 -1 -1 0; 0 -1 -1 0; 0 0 0 0] + +# for f ∈ arrays +# err,pois = Poisson_setup(MultiLevelPoisson,(2^6+2,2^6+2);f) +# @test err < 1e-6 +# @test pois.n[] ≤ 3 +# err,pois = Poisson_setup(MultiLevelPoisson,(2^4+2,2^4+2,2^4+2);f) +# @test err < 1e-6 +# @test pois.n[] ≤ 3 +# end +# end @testset "Flow.jl" begin # test than vanLeer behaves correctly @@ -170,17 +170,19 @@ end N = 4; a = zeros(N,N,2) |> f WaterLily.accelerate!(a,1,nothing,()) @test all(a .== 0) - WaterLily.accelerate!(a,1,(i,x,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,x,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,x,t) -> i==1 ? t : 2*t,(i,x,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) + WaterLily.accelerate!(b,0.,(i,x,t)->1,nothing) @test all(b .== 1) - a .= 0 # reset and accelerate using a non-uniform velocity field + 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) From 2e6a49bfcd2d4cbd25e38906dd97d6e640a26de9 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 5 Mar 2025 16:06:27 +0100 Subject: [PATCH 18/33] Added assert for g function definition in Simulation constructor. --- src/WaterLily.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index bf3abd4f..01cb96fd 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -69,6 +69,7 @@ mutable struct Simulation <: AbstractSimulation 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(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) From 05a989ddc77382b16f6d744d2a2961b0dc5a181c Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 5 Mar 2025 16:09:56 +0100 Subject: [PATCH 19/33] Reactivated tests (forgot to do it in last commit). --- test/maintests.jl | 248 +++++++++++++++++++++++----------------------- 1 file changed, 124 insertions(+), 124 deletions(-) diff --git a/test/maintests.jl b/test/maintests.jl index a65d969a..c280c84b 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -2,130 +2,130 @@ using GPUArrays using ReadVTK, WriteVTK @info "Test backends: $(join(arrays,", "))" -# @testset "util.jl" begin -# I = CartesianIndex(1,2,3,4) -# @test I+δ(3,I) == CartesianIndex(1,2,4,4) -# @test WaterLily.CI(I,5)==CartesianIndex(1,2,3,4,5) -# @test WaterLily.CIj(3,I,5)==CartesianIndex(1,2,5,4) -# @test WaterLily.CIj(2,CartesianIndex(16,16,16,3),14)==CartesianIndex(16,14,16,3) - -# @test loc(3,CartesianIndex(3,4,5)) == SVector(3,4,4.5) .- 1.5 -# I = CartesianIndex(rand(2:10,3)...) -# @test loc(0,I) == SVector(I.I...) .- 1.5 - -# ex,sym = :(a[I,i] = Math.add(p.b[I],func(I,q))),[] -# WaterLily.grab!(sym,ex) -# @test ex == :(a[I, i] = Math.add(b[I], func(I, q))) -# @test sym == [:a, :I, :i, :(p.b), :q] - -# for f ∈ arrays -# p = zeros(4,5) |> f -# apply!(x->x[1]+x[2]+3,p) # add 2×1.5 to move edge to origin -# @test inside(p) == CartesianIndices((2:3,2:4)) -# @test inside(p,buff=0) == CartesianIndices(p) -# @test L₂(p) == 187 - -# u = zeros(5,5,2) |> f -# apply!((i,x)->x[i],u) -# @test GPUArrays.@allowscalar [u[i,j,1].-(i-2) for i in 1:3, j in 1:3]==zeros(3,3) - -# Ng, D, U = (6, 6), 2, (1.0, 0.5) -# u = rand(Ng..., D) |> f # vector -# σ = rand(Ng...) |> f # scalar -# BC!(u, U) -# @test GPUArrays.@allowscalar 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 GPUArrays.@allowscalar 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]) - -# GPUArrays.@allowscalar u[end,:,1] .= 3 -# BC!(u,U,true) # save exit values -# @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3) - -# 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 all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) -# @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) -# # test 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 -# @test GPUArrays.@allowscalar all(σ[1, 2:end-1] .== σ[end-1, 2:end-1]) && all(σ[2:end-1, 1] .== σ[2:end-1, end-1]) - -# u = rand(Ng..., D) |> f # vector -# 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 interpolation -# a = zeros(5,5,2) |> f; b = zeros(5,5) |> f -# apply!((i,x)->x[i]+1.5,a); apply!(x->x[1]+1.5,b) # offset for start of grid -# @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(2.5,1),a) .≈ [2.5,1.]) -# @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(3.5,3),a) .≈ [3.5,3.]) -# @test GPUArrays.@allowscalar WaterLily.interp(SVector(2.5,1),b) ≈ 2.5 -# @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5 -# end -# end - -# function Poisson_setup(poisson,N::NTuple{D};f=Array,T=Float32) where D -# c = ones(T,N...,D) |> f; BC!(c, ntuple(zero,D)) -# x = zeros(T,N) |> f; z = copy(x) -# pois = poisson(x,c,z) -# soln = map(I->T(I.I[1]),CartesianIndices(N)) |> f -# I = first(inside(x)) -# GPUArrays.@allowscalar @. soln -= soln[I] -# z = mult!(pois,soln) -# solver!(pois) -# GPUArrays.@allowscalar @. x -= x[I] -# return L₂(x-soln)/L₂(soln),pois -# end - -# @testset "Poisson.jl" begin -# for f ∈ arrays -# err,pois = Poisson_setup(Poisson,(5,5);f) -# @test GPUArrays.@allowscalar parent(pois.D)==f(Float32[0 0 0 0 0; 0 -2 -3 -2 0; 0 -3 -4 -3 0; 0 -2 -3 -2 0; 0 0 0 0 0]) -# @test GPUArrays.@allowscalar parent(pois.iD)≈f(Float32[0 0 0 0 0; 0 -1/2 -1/3 -1/2 0; 0 -1/3 -1/4 -1/3 0; 0 -1/2 -1/3 -1/2 0; 0 0 0 0 0]) -# @test err < 1e-5 -# err,pois = Poisson_setup(Poisson,(2^6+2,2^6+2);f) -# @test err < 1e-6 -# @test pois.n[] < 310 -# err,pois = Poisson_setup(Poisson,(2^4+2,2^4+2,2^4+2);f) -# @test err < 1e-6 -# @test pois.n[] < 35 -# end -# end - -# @testset "MultiLevelPoisson.jl" begin -# I = CartesianIndex(4,3,2) -# @test all(WaterLily.down(J)==I for J ∈ WaterLily.up(I)) -# @test_throws AssertionError("MultiLevelPoisson requires size=a2ⁿ, where n>2") Poisson_setup(MultiLevelPoisson,(15+2,3^4+2)) - -# err,pois = Poisson_setup(MultiLevelPoisson,(10,10)) -# @test pois.levels[3].D == Float32[0 0 0 0; 0 -2 -2 0; 0 -2 -2 0; 0 0 0 0] -# @test err < 1e-5 - -# pois.levels[1].L[5:6,:,1].=0 -# WaterLily.update!(pois) -# @test pois.levels[3].D == Float32[0 0 0 0; 0 -1 -1 0; 0 -1 -1 0; 0 0 0 0] - -# for f ∈ arrays -# err,pois = Poisson_setup(MultiLevelPoisson,(2^6+2,2^6+2);f) -# @test err < 1e-6 -# @test pois.n[] ≤ 3 -# err,pois = Poisson_setup(MultiLevelPoisson,(2^4+2,2^4+2,2^4+2);f) -# @test err < 1e-6 -# @test pois.n[] ≤ 3 -# end -# end +@testset "util.jl" begin + I = CartesianIndex(1,2,3,4) + @test I+δ(3,I) == CartesianIndex(1,2,4,4) + @test WaterLily.CI(I,5)==CartesianIndex(1,2,3,4,5) + @test WaterLily.CIj(3,I,5)==CartesianIndex(1,2,5,4) + @test WaterLily.CIj(2,CartesianIndex(16,16,16,3),14)==CartesianIndex(16,14,16,3) + + @test loc(3,CartesianIndex(3,4,5)) == SVector(3,4,4.5) .- 1.5 + I = CartesianIndex(rand(2:10,3)...) + @test loc(0,I) == SVector(I.I...) .- 1.5 + + ex,sym = :(a[I,i] = Math.add(p.b[I],func(I,q))),[] + WaterLily.grab!(sym,ex) + @test ex == :(a[I, i] = Math.add(b[I], func(I, q))) + @test sym == [:a, :I, :i, :(p.b), :q] + + for f ∈ arrays + p = zeros(4,5) |> f + apply!(x->x[1]+x[2]+3,p) # add 2×1.5 to move edge to origin + @test inside(p) == CartesianIndices((2:3,2:4)) + @test inside(p,buff=0) == CartesianIndices(p) + @test L₂(p) == 187 + + u = zeros(5,5,2) |> f + apply!((i,x)->x[i],u) + @test GPUArrays.@allowscalar [u[i,j,1].-(i-2) for i in 1:3, j in 1:3]==zeros(3,3) + + Ng, D, U = (6, 6), 2, (1.0, 0.5) + u = rand(Ng..., D) |> f # vector + σ = rand(Ng...) |> f # scalar + BC!(u, U) + @test GPUArrays.@allowscalar 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 GPUArrays.@allowscalar 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]) + + GPUArrays.@allowscalar u[end,:,1] .= 3 + BC!(u,U,true) # save exit values + @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3) + + 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 all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) + @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) + # test 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 + @test GPUArrays.@allowscalar all(σ[1, 2:end-1] .== σ[end-1, 2:end-1]) && all(σ[2:end-1, 1] .== σ[2:end-1, end-1]) + + u = rand(Ng..., D) |> f # vector + 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 interpolation + a = zeros(5,5,2) |> f; b = zeros(5,5) |> f + apply!((i,x)->x[i]+1.5,a); apply!(x->x[1]+1.5,b) # offset for start of grid + @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(2.5,1),a) .≈ [2.5,1.]) + @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(3.5,3),a) .≈ [3.5,3.]) + @test GPUArrays.@allowscalar WaterLily.interp(SVector(2.5,1),b) ≈ 2.5 + @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5 + end +end + +function Poisson_setup(poisson,N::NTuple{D};f=Array,T=Float32) where D + c = ones(T,N...,D) |> f; BC!(c, ntuple(zero,D)) + x = zeros(T,N) |> f; z = copy(x) + pois = poisson(x,c,z) + soln = map(I->T(I.I[1]),CartesianIndices(N)) |> f + I = first(inside(x)) + GPUArrays.@allowscalar @. soln -= soln[I] + z = mult!(pois,soln) + solver!(pois) + GPUArrays.@allowscalar @. x -= x[I] + return L₂(x-soln)/L₂(soln),pois +end + +@testset "Poisson.jl" begin + for f ∈ arrays + err,pois = Poisson_setup(Poisson,(5,5);f) + @test GPUArrays.@allowscalar parent(pois.D)==f(Float32[0 0 0 0 0; 0 -2 -3 -2 0; 0 -3 -4 -3 0; 0 -2 -3 -2 0; 0 0 0 0 0]) + @test GPUArrays.@allowscalar parent(pois.iD)≈f(Float32[0 0 0 0 0; 0 -1/2 -1/3 -1/2 0; 0 -1/3 -1/4 -1/3 0; 0 -1/2 -1/3 -1/2 0; 0 0 0 0 0]) + @test err < 1e-5 + err,pois = Poisson_setup(Poisson,(2^6+2,2^6+2);f) + @test err < 1e-6 + @test pois.n[] < 310 + err,pois = Poisson_setup(Poisson,(2^4+2,2^4+2,2^4+2);f) + @test err < 1e-6 + @test pois.n[] < 35 + end +end + +@testset "MultiLevelPoisson.jl" begin + I = CartesianIndex(4,3,2) + @test all(WaterLily.down(J)==I for J ∈ WaterLily.up(I)) + @test_throws AssertionError("MultiLevelPoisson requires size=a2ⁿ, where n>2") Poisson_setup(MultiLevelPoisson,(15+2,3^4+2)) + + err,pois = Poisson_setup(MultiLevelPoisson,(10,10)) + @test pois.levels[3].D == Float32[0 0 0 0; 0 -2 -2 0; 0 -2 -2 0; 0 0 0 0] + @test err < 1e-5 + + pois.levels[1].L[5:6,:,1].=0 + WaterLily.update!(pois) + @test pois.levels[3].D == Float32[0 0 0 0; 0 -1 -1 0; 0 -1 -1 0; 0 0 0 0] + + for f ∈ arrays + err,pois = Poisson_setup(MultiLevelPoisson,(2^6+2,2^6+2);f) + @test err < 1e-6 + @test pois.n[] ≤ 3 + err,pois = Poisson_setup(MultiLevelPoisson,(2^4+2,2^4+2,2^4+2);f) + @test err < 1e-6 + @test pois.n[] ≤ 3 + end +end @testset "Flow.jl" begin # test than vanLeer behaves correctly From 35c1b23230f70d003a4c195d542a689294d42303 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 5 Mar 2025 16:57:40 +0100 Subject: [PATCH 20/33] Removed EllipsisNotation from Flow.jl (not needed anymore) and added it to utils.jl (required by interp). --- src/Flow.jl | 1 - src/util.jl | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Flow.jl b/src/Flow.jl index 483aca36..f26a39d7 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -59,7 +59,6 @@ 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,t,g,U) diff --git a/src/util.jl b/src/util.jl index 6a55698b..47559e9b 100644 --- a/src/util.jl +++ b/src/util.jl @@ -253,6 +253,7 @@ 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) From c9a4db2fbf00be2386fe12eb68a72fa02e8d501d Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Fri, 7 Mar 2025 13:36:38 +0100 Subject: [PATCH 21/33] Uses arg names instead of underscore because it breaks (unsure why not caught by tests). Also renamed the rate-of-strain tensor function to S. --- src/Flow.jl | 6 +++--- src/Metrics.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index f26a39d7..70b14fed 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -64,7 +64,7 @@ upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I Adds a space-time acceleration field g(i,x,t) and velocity field U(i,x,t) such that `rᵢ += gᵢ(i,x,t)+dUᵢ(i,x,t)/dt` at time `t=sum(dt)` to vector field `r`. """ -accelerate!(_,_,::Nothing,::Union{Nothing,Tuple}) = nothing +accelerate!(r,t,::Nothing,::Union{Nothing,Tuple}) = nothing accelerate!(r,t,f) = for i ∈ 1:last(size(r)) @loop r[I,i] += f(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) end @@ -181,5 +181,5 @@ end 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!(_,::Nothing,_) = nothing -udf!(flow::Flow,force!::Function,t; kwargs...) = force!(flow,t; kwargs...) +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..1646ab84 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) From 5cb8752a8fbe4943e8e3dd704341c1f5bc10f19a Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Fri, 7 Mar 2025 14:26:31 +0100 Subject: [PATCH 22/33] Fixed tests. --- test/maintests.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/maintests.jl b/test/maintests.jl index c280c84b..90db42d7 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -381,16 +381,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(WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ 0) + @test GPUArrays.@allowscalar all(WaterLily.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(WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[2 0; 0 2]) + @test GPUArrays.@allowscalar all(WaterLily.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(WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[0 2; 2 0]) + @test GPUArrays.@allowscalar all(WaterLily.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) From 3026334ac210d1cf102ebaf41bae7116f8b79901 Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Mon, 10 Mar 2025 11:55:49 +0100 Subject: [PATCH 23/33] Fixed viscous force tests. --- src/Metrics.jl | 2 +- test/maintests.jl | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Metrics.jl b/src/Metrics.jl index 1646ab84..b8f6afe7 100644 --- a/src/Metrics.jl +++ b/src/Metrics.jl @@ -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/test/maintests.jl b/test/maintests.jl index 90db42d7..ac23a3fd 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -381,16 +381,16 @@ import WaterLily: × # stress tensor u₂ = zeros(N,N,2) |> f u₃ = zeros(N,N,N,3) |> f - @test GPUArrays.@allowscalar all(WaterLily.S(CartesianIndex(N÷2,N÷2),u₂) .≈ 0) - @test GPUArrays.@allowscalar all(WaterLily.S(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.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[2 0; 0 2]) - @test GPUArrays.@allowscalar all(WaterLily.S(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.S(CartesianIndex(N÷2,N÷2),u₂) .≈ SA[0 2; 2 0]) - @test GPUArrays.@allowscalar all(WaterLily.S(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) From 341f046c26e8f5f68c494a89ea18db37c6eb248a Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Tue, 11 Mar 2025 14:35:18 +0100 Subject: [PATCH 24/33] Added type stable check for g function in Simulation constructor. Fixed accelerating flow test accordingly. --- src/WaterLily.jl | 1 + test/maintests.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 01cb96fd..4597eb85 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -70,6 +70,7 @@ mutable struct Simulation <: AbstractSimulation @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(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" + isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N}),zero(T)),N)).==T) "`g` is not type stable" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) diff --git a/test/maintests.jl b/test/maintests.jl index ac23a3fd..b3d02b17 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -303,7 +303,7 @@ function acceleratingFlow(N;use_g=false,T=Float64,perdir=(1,),jerk=4,mem=Array) # 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,x,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 From c80d19852b798a7958a212284991e7d78bc9488d Mon Sep 17 00:00:00 2001 From: marinlauber Date: Wed, 12 Mar 2025 12:31:06 +0100 Subject: [PATCH 25/33] improve tests for non-uniform BCs and improve error message --- src/WaterLily.jl | 2 +- test/maintests.jl | 29 ++++++++++++++++++++++++++--- 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 4597eb85..aea545c4 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -67,7 +67,7 @@ mutable struct Simulation <: AbstractSimulation Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) where N - @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` and `uλ` cannot be both specified as Function" + @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both" @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N}),zero(T)),N)).==T) "`g` is not type stable" diff --git a/test/maintests.jl b/test/maintests.jl index b3d02b17..8cb6be25 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -49,8 +49,8 @@ using ReadVTK, WriteVTK Ubc(i,x,t) = i==1 ? 1.0 : 0.5 v = rand(Ng..., D) |> f # vector BC!(v,Ubc,false); BC!(u,U,false) # make sure we apply the same - @test all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1]) - @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2]) + @test GPUArrays.@allowscalar all(v[1, :, 1] .== u[1, :, 1]) && all(v[2, :, 1] .== u[2, :, 1]) && all(v[end, :, 1] .== u[end, :, 1]) + @test GPUArrays.@allowscalar all(v[:, 1, 2] .== u[:, 1, 2]) && all(v[:, 2, 2] .== u[:, 2, 2]) && all(v[:, end, 2] .== u[:, end, 2]) # test exit bc GPUArrays.@allowscalar v[end,:,1] .= 3 BC!(v,Ubc,true) # save exit values @@ -65,6 +65,20 @@ using ReadVTK, WriteVTK BC!(u,U,true,(1,)) #saveexit has no effect here as x-periodic @test GPUArrays.@allowscalar all(u[1:2, :, 1] .== u[end-1:end, :, 1]) && all(u[1:2, :, 2] .== u[end-1:end, :, 2]) && all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) + # test non-uniform BCs + Ubc_1(i,x,t) = i==1 ? x[2] : x[1] + v .= 0; BC!(v,Ubc_1) + # the tangential BC change the value of the ghost cells on the other axis, so we cannot check it + @test GPUArrays.@allowscalar all(v[1,2:end-1,1] .≈ v[end,2:end-1,1]) + @test GPUArrays.@allowscalar all(v[2:end-1,1,2] .≈ v[2:end-1,end,2]) + # more complex + Ng, D = (8, 8, 8), 3 + u = zeros(Ng..., D) |> f # vector + Ubc_2(i,x,t) = i==1 ? cos(2π*x[1]/8) : i==2 ? sin(2π*x[2]/8) : tan(π*x[3]/16) + BC!(u,Ubc_2) + @test GPUArrays.@allowscalar all(u[1,:,:,1] .≈ cos(-1π/4)) && all(u[2,:,:,1] .≈ cos(0)) && all(u[end,:,:,1] .≈ cos(6π/4)) + @test GPUArrays.@allowscalar all(u[:,1,:,2] .≈ sin(-1π/4)) && all(u[:,2,:,2] .≈ sin(0)) && all(u[:,end,:,2] .≈ sin(6π/4)) + @test GPUArrays.@allowscalar all(u[:,:,1,3] .≈ tan(-1π/16)) && all(u[:,:,2,3] .≈ tan(0)) && all(u[:,:,end,3].-tan(6π/16).<1e-6) # test interpolation a = zeros(5,5,2) |> f; b = zeros(5,5) |> f @@ -334,7 +348,16 @@ end ) end end - +@testset "Boundary Layer Flow" begin + for f ∈ arrays + make_bl_flow(L=32;T=Float32) = Simulation((L,L), + (i,x,t)-> i==1 ? convert(Float32,4.0*(((x[2]+0.5)/2L)-((x[2]+0.5)/2L)^2)) : 0.f0, + L;ν=0.001,U=1,mem=f,T,exitBC=false) # fails with exitBC=true, but the profile is maintained + sim = make_bl_flow(32) + sim_step!(sim,10) + @test GPUArrays.@allowscalar all(sim.flow.u[1,:,1] .≈ sim.flow.u[end,:,1]) + end +end @testset "Circle in accelerating flow" begin for f ∈ arrays make_accel_circle(radius=32,H=16) = Simulation(radius.*(2H,2H), From 311b5b1941c61ba087ca81677cd08e3a5f411055 Mon Sep 17 00:00:00 2001 From: Gabriel Weymouth Date: Fri, 14 Mar 2025 17:04:05 +0100 Subject: [PATCH 26/33] Change tangential BCs and add rotating Ref test --- src/Flow.jl | 6 ++---- src/util.jl | 4 ++-- test/maintests.jl | 20 ++++++++++++++++++++ 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 70b14fed..ea4229c7 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -62,12 +62,10 @@ upperBoundary!(r,u,Φ,ν,i,j,N,::Val{true}) = @loop r[I-δ(j,I),i] -= Φ[CIj(j,I """ accelerate!(r,t,g,U) -Adds a space-time acceleration field g(i,x,t) and velocity field U(i,x,t) such that `rᵢ += gᵢ(i,x,t)+dUᵢ(i,x,t)/dt` at time `t=sum(dt)` to vector field `r`. +Accounts for applied and reference-frame acceleration using `rᵢ += g(i,x,t)+dU(i,x,t)/dt` """ accelerate!(r,t,::Nothing,::Union{Nothing,Tuple}) = nothing -accelerate!(r,t,f) = for i ∈ 1:last(size(r)) - @loop r[I,i] += f(i,loc(i,I,eltype(r)),t) over I ∈ CartesianIndices(Base.front(size(r))) -end +accelerate!(r,t,f::Function) = @loop r[Ii] += f(last(Ii),loc(Ii,eltype(r)),t) over Ii ∈ CartesianIndices(r) accelerate!(r,t,g::Function,::Union{Nothing,Tuple}) = accelerate!(r,t,g) accelerate!(r,t,::Nothing,U::Function) = accelerate!(r,t,(i,x,t)->ForwardDiff.derivative(τ->U(i,x,τ),t)) accelerate!(r,t,g::Function,U::Function) = accelerate!(r,t,(i,x,t)->g(i,x,t)+ForwardDiff.derivative(τ->U(i,x,τ),t)) diff --git a/src/util.jl b/src/util.jl index 47559e9b..9e3b86f0 100644 --- a/src/util.jl +++ b/src/util.jl @@ -203,8 +203,8 @@ function BC!(a,u_BC::Function,saveexit=false,perdir=(),t=0) end (!saveexit || i>1) && (@loop a[I,i] = u_BC(i,loc(i,I),t) over I ∈ slice(N,N[j],j)) # overwrite exit else # Tangential directions, Neumann - @loop a[I,i] = a[I+δ(j,I),i] over I ∈ slice(N,1,j) - @loop a[I,i] = a[I-δ(j,I),i] over I ∈ slice(N,N[j],j) + @loop a[I,i] = u_BC(i,loc(i,I),t)+a[I+δ(j,I),i]-u_BC(i,loc(i,I+δ(j,I)),t) over I ∈ slice(N,1,j) + @loop a[I,i] = u_BC(i,loc(i,I),t)+a[I-δ(j,I),i]-u_BC(i,loc(i,I-δ(j,I)),t) over I ∈ slice(N,N[j],j) end end end diff --git a/test/maintests.jl b/test/maintests.jl index 8cb6be25..3dc5f25a 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -358,6 +358,26 @@ end @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₂(simg.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), From 9a7753d5b007ba2c382a9eef2097d54a1b0a3a00 Mon Sep 17 00:00:00 2001 From: Gabriel Weymouth Date: Tue, 18 Mar 2025 09:56:35 +0100 Subject: [PATCH 27/33] fix typo in rotating reference test --- test/maintests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/maintests.jl b/test/maintests.jl index 3dc5f25a..f3f12039 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -375,7 +375,7 @@ 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₂(simg.flow.p)<3e-3 # should be zero + @test L₂(simg.flow.p)==L₂(sim.flow.p)<3e-3 # should be zero end @testset "Circle in accelerating flow" begin From 8c55457470a5e21c7d05b8e81c504534f2427ada Mon Sep 17 00:00:00 2001 From: marinlauber Date: Thu, 27 Mar 2025 12:13:20 +0100 Subject: [PATCH 28/33] check type with x,t not only t --- src/WaterLily.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index aea545c4..779c5932 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -70,9 +70,9 @@ mutable struct Simulation <: AbstractSimulation @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both" @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" - isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N}),zero(T)),N)).==T) "`g` is not type stable" + isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`g` is not type stable" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" - isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N}),zero(T)),N)).==T) "`u_BC` is not type stable" + isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`u_BC` is not type stable" isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) isnothing(U) && (U = √sum(abs2,u_BC)) flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) From 561db0dc86ca24efdfd371c5596a63a49d871ed5 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 28 Mar 2025 09:27:53 +0100 Subject: [PATCH 29/33] =?UTF-8?q?remove=20u=CE=BB=20from=20Simulation=20an?= =?UTF-8?q?d=20Flow?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Flow.jl | 5 +++-- src/WaterLily.jl | 6 +++--- test/maintests.jl | 4 ++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index ea4229c7..987b1d77 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -98,10 +98,11 @@ struct Flow{D, T, Sf<:AbstractArray{T}, Vf<:AbstractArray{T}, Tf<:AbstractArray{ 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=Float32) where D + perdir=(), exitBC=false, T=Float32) where D Ng = N .+ 2 Nd = (Ng..., D) - u = Array{T}(undef, Nd...) |> f; apply!(uλ, u); + u = Array{T}(undef, Nd...) |> f + isa(U, Function) ? apply!((i,x)->U(i,x,0.), u) : apply!((i,x)->U[i], u) BC!(u,U,exitBC,perdir); exitBC!(u,u,0.) u⁰ = copy(u); fv, p, σ = zeros(T, Nd) |> f, zeros(T, Ng) |> f, zeros(T, Ng) |> f diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 779c5932..3d16007a 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -67,15 +67,15 @@ mutable struct Simulation <: AbstractSimulation Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) where N - @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both" + # @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both" @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`g` is not type stable" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`u_BC` is not type stable" - isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) + # @assert !isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) isnothing(U) && (U = √sum(abs2,u_BC)) - flow = Flow(dims,u_BC;uλ,Δt,ν,g,T,f=mem,perdir,exitBC) + flow = Flow(dims,u_BC;Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) end diff --git a/test/maintests.jl b/test/maintests.jl index f3f12039..3e1339c8 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -266,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 @@ -274,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 From 957aed0b36fe4029a2ef70aaa632de0215bcb617 Mon Sep 17 00:00:00 2001 From: marinlauber Date: Fri, 28 Mar 2025 09:29:04 +0100 Subject: [PATCH 30/33] remove old checks in WaterLily.jl --- src/WaterLily.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 3d16007a..e601e23a 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -67,13 +67,11 @@ mutable struct Simulation <: AbstractSimulation Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) where N - # @assert !(isa(u_BC,Function) && isa(uλ,Function)) "`u_BC` will be used to generate `uλ=u_BC(t=0)` do not provide both" @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`g` is not type stable" isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`u_BC` is not type stable" - # @assert !isnothing(uλ) && (uλ = isa(u_BC, Function) ? uλ = (i,x)->u_BC(i,x,0) : uλ = (i,_)->u_BC[i]) isnothing(U) && (U = √sum(abs2,u_BC)) flow = Flow(dims,u_BC;Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) From bbbc15ad8416f3169cc1aaf5eaed2c10cccc28f4 Mon Sep 17 00:00:00 2001 From: weymouth Date: Sun, 30 Mar 2025 22:58:15 +0200 Subject: [PATCH 31/33] Name and doc changes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Changed name of `u_BC` to `Uλ`. Added function to check `g,Uλ`. Updated doc strings for `Simulation` and `Flow`. --- src/Flow.jl | 4 ++-- src/WaterLily.jl | 35 ++++++++++++++++++----------------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index 987b1d77..ff7130a2 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -91,10 +91,10 @@ 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 + U :: 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, diff --git a/src/WaterLily.jl b/src/WaterLily.jl index e601e23a..cc8bda89 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -30,26 +30,25 @@ 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, Uλ::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 and space-varying function `u_BC(i,x,t)` + - `Uλ`: Domain velocity applied to the initial, boundary and acceleration conditions. + Defined by `Uλ[i]::NTuple`, or function of space and time `Uλ(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. - `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,21 +62,23 @@ mutable struct Simulation <: AbstractSimulation flow :: Flow body :: AbstractBody pois :: AbstractPoisson - function Simulation(dims::NTuple{N}, u_BC, L::Number; + function Simulation(dims::NTuple{N}, Uλ, L::Number; Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), - uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), + exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) where N - @assert !(isnothing(U) && isa(u_BC,Function)) "`U` must be specified if `u_BC` is a Function" - isa(g,Function) && @assert first(methods(g)).nargs==4 "g::Function needs to be defined as g(i,x,t)" - isa(g,Function) && @assert all(typeof.(ntuple(i->g(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`g` is not type stable" - isa(u_BC,Function) && @assert first(methods(u_BC)).nargs==4 "u_BC::Function needs to be defined as u_BC(i,x,t)" - isa(u_BC,Function) && @assert all(typeof.(ntuple(i->u_BC(i,zeros(SVector{N,T}),zero(T)),N)).==T) "`u_BC` is not type stable" - isnothing(U) && (U = √sum(abs2,u_BC)) - flow = Flow(dims,u_BC;Δt,ν,g,T,f=mem,perdir,exitBC) + @assert !(isnothing(U) && isa(Uλ,Function)) "`U` must be specified if `u_BC` is a Function" + isnothing(U) && (U = √sum(abs2,Uλ)) + check_fn(g,N,T); check_fn(Uλ,N,T) + flow = Flow(dims,Uλ;Δt,ν,g,T,f=mem,perdir,exitBC) measure!(flow,body;ϵ) new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.p,flow.μ₀,flow.σ;perdir)) end end +function check_fn(f::Function,N,T) + @assert first(methods(f)).nargs==4 "$f needs to be defined as f(i,x,t)" + @assert all(typeof.(ntuple(i->f(i,zeros(SVector{N,T}),zero(T)),N)).==T) "$f is not type stable" +end +check_fn(f,N,T) = nothing time(sim::AbstractSimulation) = time(sim.flow) """ From 01d9a249aafd043a8d2ab89d5f5cf5406854de29 Mon Sep 17 00:00:00 2001 From: weymouth Date: Sun, 30 Mar 2025 23:02:24 +0200 Subject: [PATCH 32/33] =?UTF-8?q?Slightly=20better=20U=CE=BB=20doc=20strin?= =?UTF-8?q?g?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/WaterLily.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index cc8bda89..75aeb65a 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -39,8 +39,8 @@ abstract type AbstractSimulation end Constructor for a WaterLily.jl simulation: - `dims`: Simulation domain dimensions. - - `Uλ`: Domain velocity applied to the initial, boundary and acceleration conditions. - Defined by `Uλ[i]::NTuple`, or function of space and time `Uλ(i,x,t)`. + - `Uλ`: Domain velocity field applied to the initial, boundary and acceleration conditions. + Define a uniform field by `Uλ[i]::NTuple`, or supply a function of space and time `Uλ(i,x,t)`. - `L`: Simulation length scale. - `U`: Simulation velocity scale. Required if using `Uλ::Function`. - `Δt`: Initial time step. From 399061a77737f7c03a627ff114606b523d0d475b Mon Sep 17 00:00:00 2001 From: Bernat Font Date: Wed, 16 Apr 2025 19:56:38 +0200 Subject: [PATCH 33/33] =?UTF-8?q?Reverted=20behaviour=20to=20uBC=20and=20u?= =?UTF-8?q?=CE=BB=20for=20BCs=20and=20IC,=20respectively.=20Also=20cleaned?= =?UTF-8?q?=20up=20across=20Simulation=20and=20Flow.=20All=20tests=20passi?= =?UTF-8?q?ng=20locally=20on=20CPU=20and=20GPU.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Flow.jl | 27 ++++++++++++++------------- src/WaterLily.jl | 27 ++++++++++++--------------- src/util.jl | 23 +++++++++++++++++------ test/maintests.jl | 4 ++-- 4 files changed, 45 insertions(+), 36 deletions(-) diff --git a/src/Flow.jl b/src/Flow.jl index ff7130a2..4650a007 100644 --- a/src/Flow.jl +++ b/src/Flow.jl @@ -91,24 +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/function + 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} # 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, - perdir=(), exitBC=false, T=Float32) 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) + isnothing(uλ) && (uλ = ic_function(uBC)) u = Array{T}(undef, Nd...) |> f - isa(U, Function) ? apply!((i,x)->U(i,x,0.), u) : apply!((i,x)->U[i], u) - BC!(u,U,exitBC,perdir); exitBC!(u,u,0.) - u⁰ = copy(u); + 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 @@ -147,17 +148,17 @@ and the `AbstractPoisson` pressure solver to project the velocity onto an incomp @log "p" conv_diff!(a.f,a.u⁰,a.σ,ν=a.ν,perdir=a.perdir) udf!(a,udf,t₀; kwargs...) - accelerate!(a.f,t₀,a.g,a.U) - BDIM!(a); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) # BC MUST be at t₁ + 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.U,a.exitBC,a.perdir,t₁) + 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) udf!(a,udf,t₁; kwargs...) - accelerate!(a.f,t₁,a.g,a.U) - BDIM!(a); scale_u!(a,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) - project!(a,b,0.5); BC!(a.u,a.U,a.exitBC,a.perdir,t₁) + 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)) diff --git a/src/WaterLily.jl b/src/WaterLily.jl index 75aeb65a..24e721f0 100644 --- a/src/WaterLily.jl +++ b/src/WaterLily.jl @@ -30,8 +30,8 @@ include("Metrics.jl") abstract type AbstractSimulation end """ - Simulation(dims::NTuple, Uλ::Union{NTuple,Function}, L::Number; - U=norm2(Uλ), Δt=0.25, ν=0., ϵ=1, g=nothing, + 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) @@ -39,8 +39,8 @@ abstract type AbstractSimulation end Constructor for a WaterLily.jl simulation: - `dims`: Simulation domain dimensions. - - `Uλ`: Domain velocity field applied to the initial, boundary and acceleration conditions. - Define a uniform field by `Uλ[i]::NTuple`, or supply a function of space and time `Uλ(i,x,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. Required if using `Uλ::Function`. - `Δt`: Initial time step. @@ -48,6 +48,8 @@ Constructor for a WaterLily.jl simulation: - `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. - `body`: Immersed geometry. - `T`: Array element type. @@ -62,23 +64,18 @@ mutable struct Simulation <: AbstractSimulation flow :: Flow body :: AbstractBody pois :: AbstractPoisson - function Simulation(dims::NTuple{N}, Uλ, L::Number; + function Simulation(dims::NTuple{N}, uBC, L::Number; Δt=0.25, ν=0., g=nothing, U=nothing, ϵ=1, perdir=(), - exitBC=false, body::AbstractBody=NoBody(), + uλ=nothing, exitBC=false, body::AbstractBody=NoBody(), T=Float32, mem=Array) where N - @assert !(isnothing(U) && isa(Uλ,Function)) "`U` must be specified if `u_BC` is a Function" - isnothing(U) && (U = √sum(abs2,Uλ)) - check_fn(g,N,T); check_fn(Uλ,N,T) - flow = Flow(dims,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 end -function check_fn(f::Function,N,T) - @assert first(methods(f)).nargs==4 "$f needs to be defined as f(i,x,t)" - @assert all(typeof.(ntuple(i->f(i,zeros(SVector{N,T}),zero(T)),N)).==T) "$f is not type stable" -end -check_fn(f,N,T) = nothing time(sim::AbstractSimulation) = time(sim.flow) """ diff --git a/src/util.jl b/src/util.jl index 9e3b86f0..25c05189 100644 --- a/src/util.jl +++ b/src/util.jl @@ -190,7 +190,7 @@ boundary. For example `aₓ(x)=Aₓ ∀ x ∈ minmax(X)`. A zero Neumann conditi is applied to the tangential components. """ BC!(a,U,saveexit=false,perdir=(),t=0) = BC!(a,(i,x,t)->U[i],saveexit,perdir,t) -function BC!(a,u_BC::Function,saveexit=false,perdir=(),t=0) +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 @@ -199,12 +199,12 @@ function BC!(a,u_BC::Function,saveexit=false,perdir=(),t=0) else if i==j # Normal direction, Dirichlet for s ∈ (1,2) - @loop a[I,i] = u_BC(i,loc(i,I),t) 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] = u_BC(i,loc(i,I),t) 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] = u_BC(i,loc(i,I),t)+a[I+δ(j,I),i]-u_BC(i,loc(i,I+δ(j,I)),t) over I ∈ slice(N,1,j) - @loop a[I,i] = u_BC(i,loc(i,I),t)+a[I-δ(j,I),i]-u_BC(i,loc(i,I-δ(j,I)),t) over I ∈ slice(N,N[j],j) + @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 @@ -258,4 +258,15 @@ 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 3e1339c8..282b407a 100644 --- a/test/maintests.jl +++ b/test/maintests.jl @@ -332,7 +332,7 @@ end 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 @@ -341,7 +341,7 @@ end # 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.U[1] + 0.5*jerk*WaterLily.time(sim_udf)^2 + 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