Diffusivities not being calculated for model before first time-step

[tomchor]

Maybe I'm missing something, but it seems like the viscosity isn't being calculated when the model is built:

julia> using Oceananigans

julia> grid = RegularRectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (4, 4, 4)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (0.25, 0.25, 0.25)

julia> closure = ConstantSmagorinsky(ν=1e-6, κ=1e-7)
SmagorinskyLilly: C=0.23, Cb=1.0, Pr=1.0, ν=1.0e-6, κ=1.0e-7

julia> model = NonhydrostaticModel(grid=grid, closure=closure)
NonhydrostaticModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
├── tracers: (:T, :S)
├── closure: Oceananigans.TurbulenceClosures.SmagorinskyLilly{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Float64, NamedTuple{(:T, :S), Tuple{Float64, Float64}}, NamedTuple{(:T, :S), Tuple{Float64, Float64}}}
├── buoyancy: SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}
└── coriolis: Nothing

julia> interior(model.diffusivities.νₑ)
4×4×4 view(::Array{Float64, 3}, 2:5, 2:5, 2:5) with eltype Float64:
[:, :, 1] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

[:, :, 2] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

[:, :, 3] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

[:, :, 4] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

And this seems to be impacting the calculation of the diffusivity as well:

julia> κe = ComputedField(model.diffusivities.κₑ.S)
ComputedField located at (Center, Center, Center) of BinaryOperation at (Center, Center, Center)
├── data: OffsetArrays.OffsetArray{Float64, 3, Array{Float64, 3}}, size: (4, 4, 4)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=4, Ny=4, Nz=4)
├── operand: BinaryOperation at (Center, Center, Center)
└── status: time=0.0

julia> compute!(κe)

julia> interior(κe)
4×4×4 view(::Array{Float64, 3}, 2:5, 2:5, 2:5) with eltype Float64:
[:, :, 1] =
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7

[:, :, 2] =
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7

[:, :, 3] =
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7

[:, :, 4] =
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7
 -9.0e-7  -9.0e-7  -9.0e-7  -9.0e-7

Is this expected behavior? I find it kinda odd and it actually threw me off in a calculation that I'm doing for a project.

[glwagner]

It's calculated during update_state! and therefore set!:

https://github.com/CliMA/Oceananigans.jl/blob/6df0b32a3277772557f112eb2462159fbd1f53a2/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl#L49

It's also calculated in run! before starting a time-stepping loop in case set! is not used:

https://github.com/CliMA/Oceananigans.jl/blob/6df0b32a3277772557f112eb2462159fbd1f53a2/src/Simulations/run.jl#L142-L143

During time-stepping, it's called after a time step is complete.

When the model is created the velocities and tracers are initialized to 0 so I think the diffusivity is zero then (typically...)? I believe this is the case in your example --- isn't it correct that the diffusivities are zero?

To cover the case that explicit tracer and velocity fields are supplied to a model with non-zero values, we could call update_state! in the constructor --- not a bad idea at all. That just requires changing this line:

https://github.com/CliMA/Oceananigans.jl/blob/6df0b32a3277772557f112eb2462159fbd1f53a2/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl#L168

to

    model = NonhydrostaticModel(architecture, grid, clock, advection, buoyancy, coriolis, stokes_drift,
                                forcing, closure, background_fields, particles, velocities, tracers,
                                pressures, diffusivity_fields, timestepper, pressure_solver, immersed_boundary)

    update_state!(model)

    return model

[tomchor]

When the model is created the velocities and tracers are initialized to 0 so I think the diffusivity is zero then (typically...)? I believe this is the case in your example --- isn't it correct that the diffusivities are zero?

Not in this case because it ends up not adding the molecular viscosity.

To cover the case that explicit tracer and velocity fields are supplied to a model with non-zero values, we could call update_state! in the constructor --- not a bad idea at all. That just requires changing this line:

I think this makes things more intuitive, no? Are there any downsides besides a slightly longer building time for model?

[glwagner]

I think this makes things more intuitive, no? Are there any downsides besides a slightly longer building time for model?

There's no downside. update_state! has to be efficient for the model to run so I don't think there's a performance issue. I think it's more "correct", since without that call the auxiliary state may be wrong initially.