Force Differences in WaterLily Version

[mingels28]

I have noticed a difference in the force calculation outputs in WaterLily 1.2.0 and older versions. I rolled back the version of WaterLily on my install until I observed the difference and is appears to be between v1.0.4 and v1.1.0. To clarify, the forces between v1.2.0 (pressure_force) and v1.1.0 are exactly the same, but going back further to v1.0.4 is where the differences occur. I have included a simple script for the forces on a circle and some plots demonstrating the differences in lift and drag (unscaled).

using WaterLily, CUDA
function circle(n,m;Re=250,U=1,mem=Array)
    radius, center = m/8, m/2
    body = AutoBody((x,t)->√sum(abs2, x .- center) - radius)
    Simulation((n,m), (U,0), radius; ν=U*radius/Re, body, mem=mem)
end

function get_force(sim,t)
    sim_step!(sim, t, remeasure=false, verbose=true)
    #v1.2.0 force routine
    # f = WaterLily.pressure_force(sim)     

    # force routine for versions ≤v1.0.1
    sz = size(sim.flow.p)
    df = ones(Float32, tuple(sz..., length(sz))) |> CuArray
    f  = WaterLily.∮nds(sim.flow.p,df,sim.body,t*sim.L/sim.U)

	return f
end

sim = circle(2^8, 2^8, Re=5000, mem=CuArray);
t = 0:0.01:200;
f = [get_force(sim,tᵢ) for tᵢ ∈ t];

using Plots
plot(t,[first.(f), last.(f)], 
    labels=permutedims(["drag","lift"]),
    xlabel="tU/L",
    ylabel="Force",
    title="v1.2.0",
    xlims=(0,t[end]),
    ylims=(-110,110))

[b-fg]

My guess is on #129, since there were a lot of changes there, and some related to the pressure solver (which I think is the only part that can be affecting this issue).

[weymouth]

Thanks for the issue @mingels28.

I suggest we first create a really quick test script to quantify the issue. The best would be something we can add to the testset in the future. A test which only does the first time step and outputs the drag force would be good since we can compare against the analytic potential flow added mass to determine which answer is right! If it turns out that this difference only shows up after a long time (like in your example) then that's interesting in itself.

Once we have that we can identify the exact PR that caused the change and what to do about it.

[weymouth]

I made this example which does not replicate your issue.

This is a 2D flat plate added mass calculation which takes one time step and should give the result of Cₐ=F/π*R^2 = 1. The example gives Cₐ=1.05918... (with R=96) for all tagged versions back to 1.0.4.

Note that I'm using the same force routine in every version to keep it consistent. So either the difference is the force function, or this is a divergence that builds up over a long-time simulation.

using WaterLily,StaticArrays
# consistent pressure force calculation
@inline function nds(body,x,t)
    d,n,_ = measure(body,x,t) # no fast in old versions
    n*WaterLily.kern(clamp(d,-1,1))
end
pressure_force(sim) = pressure_force(sim.flow,sim.body)
pressure_force(flow,body) = pressure_force(flow.p,flow.f,body,sum(flow.Δt)-flow.Δt[end]) # time missing in some versions
function pressure_force(p,df,body,t=0,T=promote_type(Float64,eltype(p)))
    df .= zero(eltype(p))
    WaterLily.@loop df[I,:] .= p[I]*nds(body,loc(0,I),t) over I ∈ inside(p)
    sum(T,df,dims=ntuple(i->i,ndims(p)))[:] |> Array
end

# simulated + measure force
function impulsive_force(;R=32,n=8R,T=Float32,thk=T(1.5))
    sim = Simulation((n,n),(0,0),R;U=1,T,body=AutoBody(
        (x,t)->√sum(abs2,max.(0,abs.(x .- n÷2)-SA[0,R-thk]))-thk, # 2D plate
        (x,t)->x-SA[0.5f0t^2,0] # unit acceleration
    ))
    sim_step!(sim)  # get t≈0, "inviscid" pressure 
    return pressure_force(sim)/(π*R^2) # Cₐ
end
impulsive_force(R=96) # gives [1.059,0] for v1.2.0 & v1.1.0 & v1.0.4