3D simulations for several shape: flow u is NaN or 0

[f3fora]

I discovered this simulator and I am really appreciating it.

I have been playing with it and I found out that in several case, when dealing with 3D shape, like a cylinder, the sim.flow.u array is composed by only 0 and NaN.

Test it trying the donut example changing its sdf.

body = AutoBody() do xyz, t
    x, y, z = xyz - center
    # donut
    # norm2(SA[x, norm2(SA[y, z])-R]) - r 
    
    # cylinder
    d = SA[norm2(SA[y, z])-R, abs(x)-r]
    min(maximum(d), 0.0) + norm2(max.(d, 0.0))
end

Now, maybe I misunderstood something or it is just an unlucky choice of parameters, or it is a problem.

Edit:
The same problem arises with a 2D box.

[weymouth]

Can you put up a complete minimal example reproducing the problem?

[f3fora]

# disk_and_box_2D.jl
using WaterLily
using StaticArrays
using GLMakie

import CUDA
@assert CUDA.functional()

include("plot.jl")
function disk(p=6; Re=2e2, mem=Array, U=1)
    n = 2^p
    center = SA[n/2, n/2]
    R = n / 4
    ν = U * R / Re

    norm2(x) = √sum(abs2, x)
    body = AutoBody() do xyz, t
        norm2(xyz - center) - R
    end

    Simulation((2n, n), (U, 0), R; ν, body, mem), center
end

sim, center = disk(mem=CUDA.CuArray);

dat = sim.flow.σ[inside(sim.flow.σ)] |> Array;

makie_video!(sim, dat, ω_θ!, name="disk.mp4", duration=10, step=0.25) do obs, ext
    f = Figure()
    x, y, m = size(sim.flow.V)
    ax = Axis(f[1, 1], limits=(0, x - m, 0, y - m), aspect=DataAspect())

    contourf!(ax, obs, levels=ext)

    body_plot!(ax, sim)
    tightlimits!(ax)
    f
end

function box(p=6; Re=2e2, mem=Array, U=1)
    n = 2^p
    center = SA[n/2, n/2]
    R = n / 4
    r = n / 16
    ν = U * R / Re

    norm2(x) = √sum(abs2, x)
    body = AutoBody() do xyz, t
        d = abs.(xyz - center) - SA[R, r]
        norm2(max.(d, 0.0)) + min(maximum(d), 0.0)
    end

    Simulation((2n, n), (U, 0), R; ν, body, mem), center
end

sim, center = box(mem=CUDA.CuArray);

dat = sim.flow.σ[inside(sim.flow.σ)] |> Array;

makie_video!(sim, dat, ω_θ!, name="box.mp4", duration=10, step=0.25) do obs, ext
    f = Figure()
    x, y, m = size(sim.flow.V)
    ax = Axis(f[1, 1], limits=(0, x - m, 0, y - m), aspect=DataAspect())

    contourf!(ax, obs, levels=ext)

    body_plot!(ax, sim)
    tightlimits!(ax)
    f
end

# plot.jl
function ω_θ!(dat, sim, center=center)
    dt, a = sim.L / sim.U, sim.flow.σ
    @inside a[I] = WaterLily.curl(3, I, sim.flow.u) * dt
    copyto!(dat, a[inside(a)])
end

function makie_video!(makie_plot, sim, dat, obs_update!; remeasure=false, name="file.mp4", duration=1, step=0.1, framerate=30, compression=20)
    t₀ = round(sim_time(sim))
    t = range(t₀, t₀ + step; step)
    X = obs_update!(dat, sim)
    for tᵢ in t
        sim_step!(sim, tᵢ; remeasure)
        X = obs_update!(dat, sim)
    end
    obs = X |> Observable
    ext = lift(obs) do x
        i, a = extrema(x)
        LinRange(i, a, 20)
    end
    f = makie_plot(obs, ext)
    
    t = range(t₀, t₀ + duration; step)
    record(f, name, t; framerate, compression) do tᵢ
        sim_step!(sim, tᵢ; remeasure)
        X = obs_update!(dat, sim)
        obs[] = X
        println("simulation ", round(Int, (tᵢ - t₀) / duration * 100), "% complete")
        println("length: ", length(sim.flow.u), ", nan: ", count(isnan, sim.flow.u), ", zero: ", count(iszero, sim.flow.u))
    end
    return f
end

function body_plot!(ax, sim; levels=[0], R=inside(sim.flow.p))
    WaterLily.measure_sdf!(sim.flow.σ, sim.body, WaterLily.time(sim))
    contour!(ax, (sim.flow.σ[R] |> Array)'; levels)
end

[f3fora]

In the previous comment, I included two files. The latter is a modified version of the one for plotting in 3D with makie, but for 2D. The former contains two example. One is a 2D disk which behaves correctly, the other one is a rectangular box. As I said before, and also in the case of some 3D body, the flow u is either NaN or 0.

[weymouth]

A good minimum reproducible example should be as short as possible, while still reproducing the error. You should try to get rid of all the packages and functions which aren't related to the issue you're seeing.

Can you try to make an example like that?

[f3fora]

# box
using WaterLily
using StaticArrays

function ω_θ!(dat, sim)
    dt, a = sim.L / sim.U, sim.flow.σ
    @inside a[I] = WaterLily.curl(3, I, sim.flow.u) * dt
    copyto!(dat, a[inside(a)])
end

U = 1
n = 2^6
Re = 2e2
center = SA[n/2, n/2]
R = n / 4
r = n / 16
ν = U * R / Re

norm2(x) = √sum(abs2, x)
body = AutoBody() do xyz, t
    d = abs.(xyz - center) - SA[R, r]
    norm2(max.(d, 0.0)) + min(maximum(d), 0.0)
end

sim = Simulation((2n, n), (U, 0), R; ν, body, mem=Array)

dat = sim.flow.σ[inside(sim.flow.σ)] |> Array;

step = 0.1
duration = 2
t₀ = round(sim_time(sim))

for tᵢ in range(t₀, t₀ + duration; step)
    sim_step!(sim, tᵢ)
    ω_θ!(dat, sim)
    println("length: ", length(sim.flow.u), ", nan: ", count(isnan, sim.flow.u), ", zero: ", count(iszero, sim.flow.u))
end

[weymouth]

Sorry to be so Socratic, but is the problem really with ω_θ!? Do you actually require sim_step! for the error to appear?

I suggest creating the Simulation and then counting isnan within that structure. There should be none but if there are, which arrays are they in? That should isolate the problem nicely.

[f3fora]

Do you actually require sim_step! for the error to appear?

The NaN are propagated in sim_step!.

I suggest creating the Simulation and then counting isnan within that structure. There should be none but if there are, which arrays are they in? That should isolate the problem nicely.

Some NaNs appear initially in sim.flow.μ₀ and sim.flow.μ₁

using WaterLily
using StaticArrays

U = 1
n = 2^6
Re = 2e2
center = SA[n/2, n/2]
R = n / 4
r = n / 16
ν = U * R / Re

norm2(x) = √sum(abs2, x)
body = AutoBody() do xyz, t
    d = abs.(xyz - center) - SA[R, r]
    norm2(max.(d, 0.0)) + min(maximum(d), 0.0)
end

sim = Simulation((2n, n), (U, 0), R; ν, body, mem=Array)

for f in fieldnames(typeof(sim.flow))
    println(f, ":", getfield(sim.flow, f) .|> isnan |> any)
end

#=
sim.flow
u:false
u⁰:false
f:false
p:false
σ:false
V:false
σᵥ:false
μ₀:true
μ₁:true
U:false
Δt:false
ν:false
=#

[weymouth]

Ah - now we're getting down to it! But now that we know it's all in the AutoBody, we can reduce the problem much further to see that WaterLily isn't needed at all!

using ForwardDiff, StaticArrays

function sdf(x)
    r = abs.(x) .- 1
    √sum(abs2,max.(r, 0.0)) + min(maximum(r), 0.0)
end

n(x) = ForwardDiff.gradient(sdf,x)
n(SA[5,5]) # correctly gives the gradient (0.707,0.707)
n(SA[0,0]) # correctly gives NaN
n(SA[0.5,0]) # incorrectly gives NaN, should be (1,0)

The box SDF is continuous, but it's derivative isn't defined along some of the lines in the interior.

So it is correct that ForwardDiff returns NaN at some points. However, there is something else wrong as well since the entire interior is returning NaN. I think the problem is that the last line of the SDF combines discontinuous functions together (due to the max and min). These happen to cancel, but NaN-NaN doesn't equal zero in most programming languages.

Unfortunately, I don't have a good solution to this problem. So, I think your options are:

1. Avoid shapes without defined analytic derivatives.
2. Raise an issue with the ForwardDiff package and see if they can help suggest ways to rewrite the function to work better with ForwardDiff.
3. Raise an issue on the Julia discourse or slack and see if people have other suggestions for workarounds.

I am curious if you come up with something. So please let me know.

[weymouth]

Note that #77 has also "fixed" this issue. I use scare-quotes because it hasn't fixed the underlying problem with the function definition and so there are a LOT of undefined normals. AutoBody now zeros those out, meaning the simulation runs, but there may be issues with the solution accuracy.