positive tracers

[masonrogers14]

I am trying to use oceananigans to simulate concentration fields that should be nonnegative. However, even for the simple case of a constant velocity field advecting a Gaussian patch of passive tracer, the model develops steadily growing negative tracer values that eventually become comparable to the largest positive values. Here is a MWE:

using Oceananigans, Printf

#parameters---------------------------------------------------------------------
#gridding
Nx = 64
Ny = 64
Nz = 64
Lx = Nx * 1
Ly = Ny * 1
Lz = Nz * 1

#dimensional parameters
κ = 1e-6
ν = 1e-6

#simulation
sim_time = 86400.

#grid, boundary conditions, and fields------------------------------------------
#architecture
arch = Oceananigans.GPU()

#make grid
grid = RectilinearGrid(
    arch,
    size = (Nx, Ny, Nz),
    x = (-Lx/2, Lx/2), y = (-Ly/2, Ly/2), z = (-Lz, 0),
    topology = (Periodic, Periodic, Bounded),
    halo = (5, 5, 5)
)

#initialize model and simulation------------------------------------------------
#model
model = NonhydrostaticModel(; 
    grid,
    timestepper = :RungeKutta3,
    advection = WENO(grid),
    closure = ScalarDiffusivity(ν=ν, κ=κ),
    tracers = :c0,
    hydrostatic_pressure_anomaly = CenterField(grid)
)

#initial conditions
u_init = XFaceField(grid)
c_init = Field{Center, Center, Center}(grid)
set!(u_init, 0.1)
set!(c_init, (x, y, z) -> (
    exp(-(x^2 + y^2) / (Lx * Ly / 100)) * 
    exp(-z^2 / (Lz^2 / 100))
))
set!(model; u=u_init, c0=c_init)

#simulation
simulation = Simulation(model; Δt=0.5, stop_time=sim_time)
progress(sim) = @printf(
    "time: %s, max(U): %.4f, max(c): %.4e, min(c): %.4e \n",
    prettytime(sim),
    max(
        maximum(abs(sim.model.velocities.u)),
        maximum(abs(sim.model.velocities.v)),
        maximum(abs(sim.model.velocities.w))
    ),
    maximum(sim.model.tracers.c0),
    minimum(sim.model.tracers.c0)
)

#callbacks
add_callback!(simulation, progress, IterationInterval(120))

#run----------------------------------------------------------------------------
run!(simulation)

The behavior persists with any advection scheme other than first-order upwind biased, which is insufficient for more turbulent flows.

Any tips or suggestions? I heard a rumor that @simone-silvestri might know something useful...

[simone-silvestri]

Hi Mason, I would try bounding the advection scheme. We have a simple flux limiter that should preserve bounds in one-dimensions

advection = WENO(bounds = (0, 1))

It is not really bounds preserving in three-dimensions, but it should work.

[masonrogers14]

Thanks for the suggestion. Unfortunately I get a bit of a mysterious stacktrace when I try it:

[ Info: Initializing simulation...
ERROR: LoadError: InvalidIRError: compiling MethodInstance for Oceananigans.Models.NonhydrostaticModels.gpu_compute_Gc!(::KernelAbstractions.CompilerMetadata{…}, ::OffsetArrays.OffsetArray{…}, ::RectilinearGrid{…}, ::Nothing, ::Tuple{…}) resulted in invalid LLVM IR
Reason: unsupported dynamic function invocation (call to div_Uc)
Stacktrace:
 [1] tracer_tendency
   @ ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl:265
 [2] macro expansion
   @ ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl:188
 [3] gpu_compute_Gc!
   @ ~/.julia/packages/KernelAbstractions/iW1Rw/src/macros.jl:97
 [4] gpu_compute_Gc!
   @ ./none:0
Hint: catch this exception as `err` and call `code_typed(err; interactive = true)` to introspect the erronous code with Cthulhu.jl
Stacktrace:
  [1] check_ir(job::GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams}, args::LLVM.Module)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/validation.jl:147
  [2] macro expansion
    @ ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:382 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/TimerOutputs/NRdsv/src/TimerOutput.jl:253 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:381 [inlined]
  [5] emit_llvm(job::GPUCompiler.CompilerJob; toplevel::Bool, libraries::Bool, optimize::Bool, cleanup::Bool, validate::Bool, only_entry::Bool)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/utils.jl:108
  [6] emit_llvm
    @ ~/.julia/packages/GPUCompiler/2CW9L/src/utils.jl:106 [inlined]
  [7] codegen(output::Symbol, job::GPUCompiler.CompilerJob; toplevel::Bool, libraries::Bool, optimize::Bool, cleanup::Bool, validate::Bool, strip::Bool, only_entry::Bool, parent_job::Nothing)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:100
  [8] codegen
    @ ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:82 [inlined]
  [9] compile(target::Symbol, job::GPUCompiler.CompilerJob; kwargs::@Kwargs{})
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:79
 [10] compile
    @ ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:74 [inlined]
 [11] #1145
    @ ~/.julia/packages/CUDA/2kjXI/src/compiler/compilation.jl:250 [inlined]
 [12] JuliaContext(f::CUDA.var"#1145#1148"{GPUCompiler.CompilerJob{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams}}; kwargs::@Kwargs{})
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:34
 [13] JuliaContext(f::Function)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/driver.jl:25
 [14] compile(job::GPUCompiler.CompilerJob)
    @ CUDA ~/.julia/packages/CUDA/2kjXI/src/compiler/compilation.jl:249
 [15] actual_compilation(cache::Dict{Any, CUDA.CuFunction}, src::Core.MethodInstance, world::UInt64, cfg::GPUCompiler.CompilerConfig{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams}, compiler::typeof(CUDA.compile), linker::typeof(CUDA.link))
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/execution.jl:237
 [16] cached_compilation(cache::Dict{Any, CUDA.CuFunction}, src::Core.MethodInstance, cfg::GPUCompiler.CompilerConfig{GPUCompiler.PTXCompilerTarget, CUDA.CUDACompilerParams}, compiler::Function, linker::Function)
    @ GPUCompiler ~/.julia/packages/GPUCompiler/2CW9L/src/execution.jl:151
 [17] macro expansion
    @ ~/.julia/packages/CUDA/2kjXI/src/compiler/execution.jl:380 [inlined]
 [18] macro expansion
    @ ./lock.jl:273 [inlined]
 [19] cufunction(f::typeof(Oceananigans.Models.NonhydrostaticModels.gpu_compute_Gc!), tt::Type{Tuple{KernelAbstractions.CompilerMetadata{…}, OffsetArrays.OffsetArray{…}, RectilinearGrid{…}, Nothing, Tuple{…}}}; kwargs::@Kwargs{always_inline::Bool, maxthreads::Int64})
    @ CUDA ~/.julia/packages/CUDA/2kjXI/src/compiler/execution.jl:375
 [20] macro expansion
    @ ~/.julia/packages/CUDA/2kjXI/src/compiler/execution.jl:112 [inlined]
 [21] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…}; ndrange::Nothing, workgroupsize::Nothing)
    @ CUDA.CUDAKernels ~/.julia/packages/CUDA/2kjXI/src/CUDAKernels.jl:103
 [22] (::KernelAbstractions.Kernel{…})(::Field{…}, ::Vararg{…})
    @ CUDA.CUDAKernels ~/.julia/packages/CUDA/2kjXI/src/CUDAKernels.jl:89
 [23] _launch!(::GPU, ::RectilinearGrid{…}, ::Symbol, ::typeof(Oceananigans.Models.NonhydrostaticModels.compute_Gc!), ::Field{…}, ::RectilinearGrid{…}, ::Vararg{…}; exclude_periphery::Bool, reduced_dimensions::Tuple{}, active_cells_map::Nothing, only_local_halos::Bool, async::Bool)
    @ Oceananigans.Utils ~/.julia/packages/Oceananigans/p4IUv/src/Utils/kernel_launching.jl:286
 [24] _launch!
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Utils/kernel_launching.jl:268 [inlined]
 [25] launch!
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Utils/kernel_launching.jl:251 [inlined]
 [26] compute_interior_tendency_contributions!(model::NonhydrostaticModel{…}, kernel_parameters::Symbol; active_cells_map::Nothing)
    @ Oceananigans.Models.NonhydrostaticModels ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl:133
 [27] compute_interior_tendency_contributions!
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl:57 [inlined]
 [28] compute_tendencies!(model::NonhydrostaticModel{…}, callbacks::Vector{…})
    @ Oceananigans.Models.NonhydrostaticModels ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl:35
 [29] #apply_regionally!#56
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Utils/multi_region_transformation.jl:121 [inlined]
 [30] apply_regionally!
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Utils/multi_region_transformation.jl:118 [inlined]
 [31] macro expansion
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Utils/multi_region_transformation.jl:206 [inlined]
 [32] update_state!(model::NonhydrostaticModel{…}, callbacks::Vector{…}; compute_tendencies::Bool)
    @ Oceananigans.Models.NonhydrostaticModels ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl:52
 [33] update_state! (repeats 2 times)
    @ ~/.julia/packages/Oceananigans/p4IUv/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl:20 [inlined]
 [34] initialize!(sim::Simulation{NonhydrostaticModel{…}, Float64, Float64, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/p4IUv/src/Simulations/run.jl:191
 [35] time_step!(sim::Simulation{NonhydrostaticModel{…}, Float64, Float64, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/p4IUv/src/Simulations/run.jl:112
 [36] run!(sim::Simulation{NonhydrostaticModel{…}, Float64, Float64, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}}; pickup::Bool)
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/p4IUv/src/Simulations/run.jl:97
 [37] run!(sim::Simulation{NonhydrostaticModel{…}, Float64, Float64, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}, OrderedCollections.OrderedDict{…}})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/p4IUv/src/Simulations/run.jl:85
 [38] top-level scope
    @ /nobackup/users/masonr/inertial/ht_xinkai.jl:71
 [39] include(fname::String)
    @ Main ./sysimg.jl:38
 [40] top-level scope
    @ REPL[1]:1
in expression starting at /nobackup/users/masonr/inertial/ht_xinkai.jl:71
Some type information was truncated. Use `show(err)` to see complete types.

[glwagner]

@masonrogers14 can you provide your code or a simple example that illustrates it?

[masonrogers14]

Here it is (same as MWE above but with the advection scheme definition suggested by Simone):

using Oceananigans, Printf

#parameters---------------------------------------------------------------------
#gridding
Nx = 64
Ny = 64
Nz = 64
Lx = Nx * 1
Ly = Ny * 1
Lz = Nz * 1

#dimensional parameters
κ = 1e-6
ν = 1e-6

#simulation
sim_time = 86400.

#grid, boundary conditions, and fields------------------------------------------
#architecture
arch = Oceananigans.GPU()

#make grid
grid = RectilinearGrid(
    arch,
    size = (Nx, Ny, Nz),
    x = (-Lx/2, Lx/2), y = (-Ly/2, Ly/2), z = (-Lz, 0),
    topology = (Periodic, Periodic, Bounded),
    halo = (5, 5, 5)
)

#initialize model and simulation------------------------------------------------
#model
model = NonhydrostaticModel(; 
    grid,
    timestepper = :RungeKutta3,
    advection = WENO(bounds=(0, 1)),
    closure = ScalarDiffusivity(ν=ν, κ=κ),
    tracers = :c0,
    hydrostatic_pressure_anomaly = CenterField(grid)
)

#initial conditions
u_init = XFaceField(grid)
c_init = Field{Center, Center, Center}(grid)
set!(u_init, 0.1)
set!(c_init, (x, y, z) -> (
    exp(-(x^2 + y^2) / (Lx * Ly / 100)) * 
    exp(-z^2 / (Lz^2 / 100))
))
set!(model; u=u_init, c0=c_init)

#simulation
simulation = Simulation(model; Δt=0.5, stop_time=sim_time)
progress(sim) = @printf(
    "time: %s, max(U): %.4f, max(c): %.4e, min(c): %.4e \n",
    prettytime(sim),
    max(
        maximum(abs(sim.model.velocities.u)),
        maximum(abs(sim.model.velocities.v)),
        maximum(abs(sim.model.velocities.w))
    ),
    maximum(sim.model.tracers.c0),
    minimum(sim.model.tracers.c0)
)

#callbacks
add_callback!(simulation, progress, IterationInterval(120))

#run----------------------------------------------------------------------------
run!(simulation)