Noisy dependence of the body position

[weymouth]

Simple test case demonstrating issue:

[updated example]

using WaterLily,StaticArrays

function circle(center;ϵ=1,Re=100,n=3*2^6,m=2^7,U=1,T=typeof(center))
    radius = m÷8
    body = AutoBody((x,t)->√sum(abs2, x .- center) - radius)
    Simulation((n,m), (U,0), radius; ν=U*radius/Re, body,T,ϵ)
end
function circle_ke(x;kwargs...)
    sim = circle(x;kwargs...)
    sim_step!(sim,0.1)
    # sum(I->WaterLily.ke(I,sim.flow.u),inside(sim.flow.p))
    # sum(I->WaterLily.ke(I,sim.flow.u,SA[1,0]),inside(sim.flow.p))
    sim.flow.u[83,66,2] # within the BDIM-ϵ region
end

using Plots,TypedTables 
data = map(range(64,66,200)) do x
    (x=x,double=circle_ke(Float64(x)),single=circle_ke(Float32(x)))
end |> Table
data2 = map(range(64,66,200)) do x
    (x=x,double=circle_ke(Float64(x),ϵ=2),single=circle_ke(Float32(x),ϵ=2))
end |> Table

plot(data.x,data.double,label="Float64, ϵ=1",c=:blue)
plot!(data.x,data.single,label="Float32, ϵ=1",ls=:dash,c=:lightblue)
plot!(data2.x,data2.double,label="Float64, ϵ=2",c=:red)
plot!(data2.x,data2.single,label="Float32, ϵ=2",ls=:dash,c=:orange)

[b-fg]

The wavy pattern probably has a wavenumber related to BDIM's epsilon right?

[weymouth]

The amplitude of that wave is related to epsilon. The wavenumber is set by the grid spacing.

[weymouth]

There might be a problem with the Waterlily.ke function because when I subtract the background flow using

sum(I->WaterLily.ke(I,sim.flow.u,SA[1,0]),inside(sim.flow.p))

I get a very different amount of noise

[weymouth]

To separate the two issues, here is the result just sampling the velocity sim.flow.u[91,65,2], which looks basically perfect. So perhaps this is purely and issue with metrics that accumulate over fields! 😱

[b-fg]

Well, that's a good catch. Then I think there is probably noise in the bodies or the boundaries, and that should be masked out when computing integrals values.

[b-fg]

I tested to exclude boundaries and internal body cells with no improvement.

[weymouth]

Looks like the noise is strongly dependent on the finite pressure tolerance. Here is the measurement of sum(sim.flow.u[2:end-1,1:end-1,2]) for different tolerances:

The current default tol=2e-4:

Ten times smaller tol=2e-5

The floor tol=10eps(T)

Same, but without the Float32 runs so you can see the Float64 levels...

[weymouth]

So it looks like there are three sources of "noise".

1. BDIM noise with wavelength=h and amplitude dependent on epsilon. Only really noticeable in some metrics (probably mostly BL related).
2. Pressure solver noise with amplitude dependent on the solver tolerance. Not clear if there is a "wavelength" of if it's basically broadband.
3. Finite precision noise. Maybe fixable with some tricks in the accumulation when using Float32? #134

[b-fg]

That PR basically has implemented sum(Float64,...) so maybe that improves the results regarding the finite precision noise when running with Float32.

Also, to make the derivative smooth, we could first do the integrals of the AD variable with Float64 (as in the PR), but then cast it to Float32, regardless of the simulation precision. That could be a way to automatically smooth the derivative if, as shown, the noise is around O(10^-8).

[weymouth]

As far as I can tell, nothing in that PR actually improves the results for Float32.