Constant Smagorinsky halo size

[ali-ramadhan]

@jm-c was right about Smagorinsky requiring a larger halo size than the rest of the code.

Running test_smag_divflux_finiteness(T) causes a BoundsError when trying to access index -1 which suggests we need halos of at least size 2.

We do support halos of arbitrary size in all three directions so it would be easy to increase the size of the halos if Smagorinsky is used as a closure. Just wanted to open this issue to see if this is what we wanted to do.

For now I'm skipping the test so it's showing up as broken.

Constant Smagorinsky: Error During Test at /home/alir/Oceananigans.jl/test/runtests.jl:466
  Test threw exception
  Expression: test_smag_divflux_finiteness(T)
  BoundsError: attempt to access OffsetArray(::Array{Float32,3}, 0:4, 0:4, 1:3) with eltype Float32 with indices 0:4×0:4×1:3 at index [3, -1, 2]
  Stacktrace:
   [1] throw_boundserror(::OffsetArray{Float32,3,Array{Float32,3}}, ::Tuple{Int64,Int64,Int64}) at ./abstractarray.jl:484
   [2] checkbounds at ./abstractarray.jl:449 [inlined]
   [3] getindex at /home/alir/.julia/packages/OffsetArrays/ruvC7/src/OffsetArrays.jl:130 [inlined]
   [4] ∂y_afa at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:8 [inlined]
   [5] ∂y_u(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::OffsetArray{Float32,3,Array{Float32,3}}) at /home/alir/Oceananigans.jl/src/closures/velocity_gradients.jl:14
   [6] Σ₁₂ at /home/alir/Oceananigans.jl/src/closures/velocity_gradients.jl:33 [inlined]
   [7] Σ₁₂² at /home/alir/Oceananigans.jl/src/closures/velocity_gradients.jl:44 [inlined]
   [8] Σ₁₂² at /home/alir/Oceananigans.jl/src/closures/velocity_gradients.jl:80 [inlined]
   [9] ▶x_caa at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:191 [inlined]
   [10] ▶y_aca at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:217 [inlined]
   [11] ▶xy_cca at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:275 [inlined]
   [12] ΣᵢⱼΣᵢⱼ_ccc at /home/alir/Oceananigans.jl/src/closures/constant_smagorinsky.jl:32 [inlined]
   [13] κ_ccc(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::ConstantSmagorinsky{Float32}, ::LinearEquationOfState{Float64}, ::Float64, ::OffsetArray{Float32,3,Array{Float32,3}}, ::OffsetArray{Float32,3,Array{Float32,3}}, ::OffsetArray{Float32,3,Array{Float32,3}}, ::OffsetArray{Float32,3,Array{Float32,3}}, ::OffsetArray{Float32,3,Array{Float32,3}}) at /home/alir/Oceananigans.jl/src/closures/constant_smagorinsky.jl:88
   [14] ▶y_afa(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::Function, ::ConstantSmagorinsky{Float32}, ::Vararg{Any,N} where N) at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:204
   [15] κ_∂y_ϕ(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::OffsetArray{Float32,3,Array{Float32,3}}, ::Function, ::ConstantSmagorinsky{Float32}, ::LinearEquationOfState{Float64}, ::Vararg{Any,N} where N) at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:482
   [16] ∂y_aca(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::Function, ::OffsetArray{Float32,3,Array{Float32,3}}, ::Vararg{Any,N} where N) at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:114
   [17] ∇_κ_∇ϕ(::Int64, ::Int64, ::Int64, ::RegularCartesianGrid{Float32,StepRangeLen{Float32,Float64,Float64}}, ::OffsetArray{Float32,3,Array{Float32,3}}, ::ConstantSmagorinsky{Float32}, ::LinearEquationOfState{Float64}, ::Vararg{Any,N} where N) at /home/alir/Oceananigans.jl/src/closures/closure_operators.jl:506
   [18] test_smag_divflux_finiteness(::Type) at /home/alir/Oceananigans.jl/test/test_turbulence_closures.jl:173
   [19] top-level scope at /home/alir/Oceananigans.jl/test/runtests.jl:466
   [20] top-level scope at /build/julia/src/julia-1.1.0/usr/share/julia/stdlib/v1.1/Test/src/Test.jl:1083
   [21] top-level scope at /home/alir/Oceananigans.jl/test/runtests.jl:464
   [22] top-level scope at /build/julia/src/julia-1.1.0/usr/share/julia/stdlib/v1.1/Test/src/Test.jl:1083
   [23] top-level scope at /home/alir/Oceananigans.jl/test/runtests.jl:427
   [24] top-level scope at /build/julia/src/julia-1.1.0/usr/share/julia/stdlib/v1.1/Test/src/Test.jl:1083
   [25] top-level scope at /home/alir/Oceananigans.jl/test/runtests.jl:20

[glwagner]

That's a fine strategy for now since Smagorinsky cannot be used for any practical purpose at the moment on master.

This will be good to keep in mind when discussing #259: the halo or stencil size is a property of the equation being solved. This motivates implementing the equation abstraction as a type.

[glwagner]

This issue is stale since we use a different strategy to calculate nonlinear viscosities now.

[ali-ramadhan]

Can the test_smag_divflux_finiteness test be enabled or should we remove it?

[glwagner]

It can, it needs either to be modified, or for the closures themselves to be modified to avoid spitting NaNs (constant Smagorinsky has the same issue as AMD in that fields cannot be constant, which is why this test was failing). I think it makes sense to modify the closure implementation to solve this problem, since I think the resulting behavior will actually be more correct (the NaN is a numerical issue rather than a mathematical feature of the model) and the computational cost of the if statements are hopefully in the noise.