Non-determinism in examples, traced back to initial implicit residual in `ImplicitEquationJacobian`

[charleskawczynski]

Found (by luck) in #2225. I've since initialized the data to zeros, to help improve CI determinism, but if after

    ∂Yₜ∂Y = FieldMatrix(...)

we do

    fill!(parent(∂Yₜ∂Y[ᶜρ_name, ᶠ𝕄_name]), NaN)
    fill!(parent(∂Yₜ∂Y[ᶜ𝔼_name, ᶠ𝕄_name]), NaN)
    fill!(parent(∂Yₜ∂Y[ᶠ𝕄_name, ᶜρ_name]), NaN)
    fill!(parent(∂Yₜ∂Y[ᶠ𝕄_name, ᶜ𝔼_name]), NaN)
    fill!(parent(∂Yₜ∂Y[ᶠ𝕄_name, ᶠ𝕄_name]), NaN)

then the solution is (sometimes) driven to NaNs. This appears to be a bug because the solution should not depend on the initial state of the residual--it should be (properly) computed before it is used.

Here are two builds where I've observed this (I also observed this locally):

• https://buildkite.com/clima/climacore-ci/builds/5419#019582d3-de2a-490d-a2f1-a05109e36c87
• https://buildkite.com/clima/climacore-ci/builds/5419#01958518-d497-41c1-b2f2-5bacc4da043e

However, I just tried ∂Yₜ∂Y with NaNs, and the solution was NaN-free, so this seems to be non-deterministic.

Running on the CPU, with Julia 1.10, here is a "reproducer": ENV["TEST_NAME"] = "sphere/held_suarez_rhoe"; using Revise; include("examples/hybrid/driver.jl")

cc @dennisYatunin

[charleskawczynski]

Actually, I'm not 100% confident that the issue with this initialization, I used the DebugOnly module and ran with a NaN-callback, and this is where I did find NaNs, but I don't know if that's how it got into the solution. I'm going to leave this issue open for now, since it's still a bit odd. Let's see how the latest push in #2225 goes.