Skip to content

Running ocean_simulation with closure being a tuple of closures which includes CATKE throws initialization error due to tracer :e missing. #214

@xkykai

Description

@xkykai

Here's a MWE:

using Oceananigans
using ClimaOcean
using OrthogonalSphericalShellGrids
using CUDA

z_faces = exponential_z_faces(Nz=40, depth=6000)

Nx = 360
Ny = 180
Nz = length(z_faces) - 1

grid = TripolarGrid(CPU();
                    size = (Nx, Ny, Nz),
                    halo = (7, 7, 7),
                    z = z_faces, 
                    north_poles_latitude = 55,
                    first_pole_longitude = 70)

bottom_height = regrid_bathymetry(grid; 
                                  minimum_depth = 10,
                                  interpolation_passes = 5,
                                  major_basins = 3,
                                  dir = "./scratch")

CUDA.@allowscalar interior(bottom_height)[isnan.(bottom_height)] .= 0;
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))

closure = (Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew = 1000.0, κ_symmetric = 1000.0), 
           ClimaOcean.OceanSimulations.default_ocean_closure())
ocean = ocean_simulation(grid; closure)

which throws

ERROR: ArgumentError: Tracers must contain :e to represent turbulent kinetic energy for `CATKEVerticalDiffusivity`.
Stacktrace:
  [1] with_tracers
    @ C:\Users\xinle\.julia\packages\Oceananigans\IsF1o\src\TurbulenceClosures\turbulence_closure_implementations\TKEBasedVerticalDiffusivities\catke_vertical_diffusivity.jl:123 [inlined]
  [2] #22
    @ .\none:0 [inlined]
  [3] iterate
    @ .\generator.jl:47 [inlined]
  [4] collect_to!(dest::Vector{…}, itr::Base.Generator{…}, offs::Int64, st::Int64)
    @ Base .\array.jl:892
  [5] collect_to_with_first!
    @ .\array.jl:870 [inlined]
  [6] collect
    @ .\array.jl:844 [inlined]
  [7] _totuple
    @ .\tuple.jl:425 [inlined]
  [8] Tuple
    @ .\tuple.jl:391 [inlined]
  [9] with_tracers(tracers::Tuple{…}, closure_tuple::Tuple{…})
    @ Oceananigans.TurbulenceClosures C:\Users\xinle\.julia\packages\Oceananigans\IsF1o\src\TurbulenceClosures\closure_tuples.jl:74
 [10] HydrostaticFreeSurfaceModel(; grid::ImmersedBoundaryGrid{…}, clock::Clock{…}, momentum_advection::VectorInvariant{…}, tracer_advection::FluxFormAdvection{…}, buoyancy::SeawaterBuoyancy{…}, coriolis::HydrostaticSphericalCoriolis{…}, free_surface::SplitExplicitFreeSurface{…}, tracers::Tuple{…}, forcing::@NamedTuple{}, closure::Tuple{…}, boundary_conditions::@NamedTuple{}, particles::Nothing, biogeochemistry::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::@NamedTuple{})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels C:\Users\xinle\.julia\packages\Oceananigans\IsF1o\src\Models\HydrostaticFreeSurfaceModels\hydrostatic_free_surface_model.jl:175
 [11] ocean_simulation(grid::ImmersedBoundaryGrid{…}; Δt::Float64, closure::Tuple{…}, free_surface::SplitExplicitFreeSurface{…}, reference_density::Int64, rotation_rate::Float64, gravitational_acceleration::Float64, bottom_drag_coefficient::ClimaOcean.OceanSimulations.Default{…}, forcing::@NamedTuple{}, coriolis::HydrostaticSphericalCoriolis{…}, momentum_advection::VectorInvariant{…}, tracer_advection::FluxFormAdvection{…}, verbose::Bool)
    @ ClimaOcean.OceanSimulations C:\Users\xinle\.julia\packages\ClimaOcean\l5TH1\src\OceanSimulations\OceanSimulations.jl:122
 [12] top-level scope
    @ c:\Users\xinle\MIT\NN_global_simulation\global_simulation copy.jl:30
Some type information was truncated. Use `show(err)` to see complete types.

The issue is https://github.com/CliMA/ClimaOcean.jl/blob/cd2677b91b4f1a6469fda9ded13f937333760997/src/OceanSimulations/OceanSimulations.jl#L116C1-L120C8 does not allow for closure being a tuple.

Suggested change could be either checking for whether closure is a tuple eg.

    tracers = (:T, :S)
    if closure isa Tuple
        for clo in closure
            if clo isa CATKEVerticalDiffusivity
                tracers = tuple(tracers..., :e)
                tracer_advection = (; T = tracer_advection, S = tracer_advection, e = nothing)
            end
        end
    else
        if closure isa CATKEVerticalDiffusivity
            tracers = tuple(tracers..., :e)
            tracer_advection = (; T = tracer_advection, S = tracer_advection, e = nothing)
        end
    end

or allowing the user to freely define tracers as a function argument.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions