Setting to `ECCOMetadata` fails on a coarse 2 degree grid

[glwagner]

This code:

using ClimaOcean
using Oceananigans

# 2 degree grid
grid = LatitudeLongitudeGrid(size = (180, 75, 1),
                             halo = (7, 7, 7),
                             z = (-10, 0),
                             latitude  = (-75, 75),
                             longitude = (0, 360))

bottom_height = regrid_bathymetry(grid;
                                  minimum_depth = 10,
                                  interpolation_passes = 5,
                                  major_basins = 3)

grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))

atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory())
ocean = ocean_simulation(grid)

set!(ocean.model; T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity))

fails with

ERROR: LoadError: DomainError with -1.244133350327736e21:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:686 [inlined]
  [3] s
    @ ~/.julia/packages/SeawaterPolynomials/F5A5Q/src/TEOS10.jl:83 [inlined]
  [4] thermal_sensitivity
    @ ~/.julia/packages/SeawaterPolynomials/F5A5Q/src/TEOS10.jl:261 [inlined]

which I think must have to do with the fact that the ECCO data is not correctly masked / inpainted when it is set.

[simone-silvestri]

Is this still an issue? I cannot reproduce it on main:

julia> using ClimaOcean

julia> using Oceananigans

       # 2 degree grid

julia> grid = LatitudeLongitudeGrid(size = (180, 75, 1),
                                    halo = (7, 7, 7),
                                    z = (-10, 0),
                                    latitude  = (-75, 75),
                                    longitude = (0, 360))
180×75×1 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 7×7×7 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [0.0, 360.0)  regularly spaced with Δλ=2.0
├── latitude:  Bounded  φ ∈ [-75.0, 75.0] regularly spaced with Δφ=2.0
└── z:         Bounded  z ∈ [-10.0, 0.0]  regularly spaced with Δz=10.0

julia> bottom_height = regrid_bathymetry(grid;
                                         minimum_depth = 10,
                                         interpolation_passes = 5,
                                         major_basins = 3)
[ Info: Regridding bathymetry from existing file /Users/simonesilvestri/.julia/scratchspaces/0376089a-ecfe-4b0e-a64f-9c555d74d754/Bathymetry/ETOPO_2022_v1_60s_N90W180_surface.nc.
180×75×1 Field{Center, Center, Nothing} reduced over dims = (3,) on LatitudeLongitudeGrid on CPU
├── grid: 180×75×1 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 7×7×7 halo and with precomputed metrics
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: ZeroFlux, north: ZeroFlux, bottom: Nothing, top: Nothing, immersed: ZeroFlux
└── data: 194×89×1 OffsetArray(::Array{Float64, 3}, -6:187, -6:82, 1:1) with eltype Float64 with indices -6:187×-6:82×1:1
    └── max=0.0, min=-9225.52, mean=-2516.54

julia> grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
180×75×1 ImmersedBoundaryGrid{Float64, Periodic, Bounded, Bounded} on CPU with 7×7×7 halo:
├── immersed_boundary: GridFittedBottom(mean(z)=-2516.54, min(z)=-9225.52, max(z)=0.0)
├── underlying_grid: 180×75×1 LatitudeLongitudeGrid{Float64, Periodic, Bounded, Bounded} on CPU with 7×7×7 halo and with precomputed metrics
├── longitude: Periodic λ ∈ [0.0, 360.0)  regularly spaced with Δλ=2.0
├── latitude:  Bounded  φ ∈ [-75.0, 75.0] regularly spaced with Δφ=2.0
└── z:         Bounded  z ∈ [-10.0, 0.0]  regularly spaced with Δz=10.0

julia> atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory())
PrescribedAtmosphere

julia> ocean = ocean_simulation(grid)
┌ Warning: This simulation will run forever as stop iteration = stop time = wall time limit = Inf.
└ @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/Hkk5J/src/Simulations/simulation.jl:55
Simulation of HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: Inf days
├── Stop iteration : Inf
├── Wall time limit: Inf
├── Callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries

julia> set!(ocean.model; T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity))
[ Info: In-painting ecco temperature
[ Info: In-painting ecco salinity

If I heatmap the initial conditions I get

julia> fig = Figure(); axT = Axis(fig[1, 1], title = "Initial temperature"); axS = Axis(fig[1, 2], title = "Initial Salinity")
Axis with 0 plots:

julia> heatmap!(axT, ocean.model.tracers.T); heatmap!(axS, ocean.model.tracers.S)
Heatmap{Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float32}}}