Climatemachine.jl: [Ocean] ocean_gyre.jl experiment produces NaNs

[glwagner]

Description

The experiments/ocean_gyre.jl experiment produces NaNs:

julia> include("ocean_gyre.jl")
[ Info: CUDAdrv.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)
[ Info: Initializing ocean_gyre
┌ Info: Starting ocean_gyre
│ dt = 5.92728e+02
│ timeend = 2592000.00
│ number of steps = 4373
└ norm(Q) = 6.9665869213940251e+08
┌ Info: Update
│ simtime = 7705.47 / 2592000.00
│ runtime = 00:01:01
└ norm(Q) = NaN
┌ Info: Update
│ simtime = 17781.84 / 2592000.00
│ runtime = 00:02:04
└ norm(Q) = NaN
┌ Info: Update
│ simtime = 27858.22 / 2592000.00
│ runtime = 00:03:06
└ norm(Q) = NaN
┌ Info: Update
│ simtime = 37934.60 / 2592000.00
│ runtime = 00:04:08
└ norm(Q) = NaN

Steps to reproduce

Steps to reproduce the behavior:

First I added CLIMA to my julia environment

pkg> https://github.com/climate-machine/CLIMA.git

at the same time, I cloned it onto my laptop. When I cloned it, the output of git log -1 was

~/Projects/CLIMA$ git log -1
commit cf990ffb758b663b33bfeebc448b0596127b4da7 (HEAD -> master, origin/master, origin/HEAD)
Merge: 454aed78 ccefa332
Author: bors[bot] <26634292+bors[bot]@users.noreply.github.com>
Date:   Wed Mar 18 22:30:19 2020 +0000

    Merge #857

    857: Rename some tests r=charleskawczynski a=charleskawczynski



    Co-authored-by: Charles Kawczynski <[email protected]>

I was able to precompile:

julia> using CLIMA
[ Info: Precompiling CLIMA [777c4786-024f-11e9-21a3-85d5d4106250]
┌ Warning: Replacing docs for `CLIMA.DGmethods.knl_nodal_update_aux! :: Union{}` in module `CLIMA.DGmethods`
└ @ Base.Docs docs/Docs.jl:223
[ Info: CUDAdrv.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)

I then navigated to /experiments, opened julia, and ran:

julia> include("ocean_gyre.jl")
[ Info: CUDAdrv.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)
[ Info: Initializing ocean_gyre
┌ Info: Starting ocean_gyre
│ dt              = 5.92728e+02
│ timeend         = 2592000.00
│ number of steps = 4373
└ norm(Q)         = 6.9665869213940251e+08
┌ Info: Update
│ simtime =  7705.47 / 2592000.00
│ runtime = 00:01:01
└ norm(Q) = NaN
┌ Info: Update
│ simtime = 17781.84 / 2592000.00
│ runtime = 00:02:04
└ norm(Q) = NaN

Expected behavior

I expected the model to run without producing NaNs.

However, I also could not deduce what norm(Q) is I'm not sure what to expect, exactly. The script does not produce any plots.

For CLIMA developers

• [x] The bug can be reproduced
• [x] Labels have been added and developers assigned
• [ ] The bug can be fixed (if not, please explain why not in a comment below)
• [ ] There is a timeline to fix the bug
• [ ] The bug has been fixed (Please link the PR)

[jm-c]

@glwagner there have been discussion about this issue (before today daily zoom session,
during and after) so it's a little bit misleading not to mention all these exchanges.
Since it was working fine before (not sure when it got broken but certainly was working on March 04), I would be tempted to ask the person who made code changes (and we could go back in time to check exactly which PR merge did it) to check and fix this problem.

[jm-c]

to be more specific: with git tag "545c51f17b3c0d67a49b174dab9af9405ee4f33c" (from March 06. 13:16)
I can run: test/Ocean/Hydrostatic_Boussinesq/simple_box.jl
but not: experiments/OceanBoxGCM/ocean_gyre.jl (which seems to have been broken at that time).
If someone know the latest git tag that is still working, this could help.

[jm-c]

May be it's the time-step ?
With this tag: 2959529adf93d652141d9d5b33f504d4a36e5f86
(from March 11, after merging PR #802) I can run: experiments/OceanBoxGCM/ocean_gyre.jl
and it prints: dt = 1.39467e+02
And with this tag: bff4887b2a17c1f0f82bfb8811f726a7cd6c11d7
(from March 12, after merging PR 833) when I try to run ocean_gyre.jl,
it prints dt = 5.92728e+02 and got NaN very soon.

[glwagner]

Thanks @jm-c. A link to the source code at that time is:

https://github.com/climate-machine/CLIMA/tree/545c51f17b3c0d67a49b174dab9af9405ee4f33c

@blallen is it possible to set the time step size in ocean_gyre.jl? There doesn't appear to be a reference to the time step size in the file currently so I'm not sure where to start.

[blallen]

in the solver_config constructor you can pass the optional argument ode_dt to set the time step manually

[glwagner]

setting ode_dt = 1e2 in CLIMA.setup_solver fixed this problem for me.

How is the time-step determined if it is not set by the user?

This is the script that worked for me:

using CLIMA, CLIMA.HydrostaticBoussinesq, CLIMA.Mesh.Filters

using CLIMA.PlanetParameters: grav
using CLIMA.Mesh.Grids: polynomialorder

CLIMA.init()

# Box dimensions in meters
Lx = 4e6    # m 
Ly = 4e6    # m 
Lz = 1000   # m 

# Problem resolution
polynomial_order = 4 
number_of_elements = (20, 20, 20) 

problem = OceanGyre{Float64}(Lx, Ly, Lz) 
model = HydrostaticBoussinesqModel{Float64}(problem, cʰ=sqrt(grav*Lz))
driver_configuration = CLIMA.OceanBoxGCMConfiguration("ocean_gyre", polynomial_order, 
                                                      number_of_elements, model)

grid = driver_configuration.grid

nodal_filters = (vert_filter = CutoffFilter(grid, polynomialorder(grid) - 1), 
                  exp_filter = ExponentialFilter(grid, 1, 8)) 

day = 86400.0
start_time = 0.0 
stop_time = 30day

solver_configuration = CLIMA.setup_solver(start_time, stop_time, driver_configuration,
        init_on_cpu = true,
             ode_dt = 1e2, 
    ode_solver_type = CLIMA.IMEXSolverType(linear_model=LinearHBModel),
          modeldata = nodal_filters
)

result = CLIMA.invoke!(solver_configuration)

Some output:

julia> include("ocean_gyre.jl")
[ Info: CUDAdrv.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)
[ Info: Initializing ocean_gyre
┌ Info: Starting ocean_gyre
│ dt              = 1.00000e+02
│ timeend         = 2592000.00
│ number of steps = 25920
└ norm(Q)         = 6.9665869213940251e+08
┌ Info: Update
│ simtime =  1600.00 / 2592000.00
│ runtime = 00:01:01
└ norm(Q) = 6.9677165893633997e+08
┌ Info: Update
│ simtime =  3400.00 / 2592000.00
│ runtime = 00:02:02
└ norm(Q) = 6.9692173530689240e+08

[blallen]

hmmm that means the calculate_dt functions in src/Ocean/Courant.jl aren't calculating the right time step...

the rest of the physics should be good tho

[christophernhill]

@glwagner etc... on CPU RHEL 8 system (tartarus.mit.edu) I get similar, but on GPU I get to NaN fairly quickly.... - hmmmm??? I think GPU run is longer, so maybe just stability.

[cnh@TARTARUS CLIMA]$ diff -b foo.jl foo-cpu.jl 
6c6
< CLIMA.init()
---
> CLIMA.init(disable_gpu=true)

[cnh@TARTARUS CLIMA]$ julia --project=@. foo.jl 
[TARTARUS.MIT.EDU:27185] mca_base_component_repository_open: unable to open mca_btl_uct: /usr/lib64/openmpi/lib/openmpi/mca_btl_uct.so: undefined symbol: uct_release_md_resource_list (ignored)
[1584672363.081868] [TARTARUS:27185:0]         parser.c:1600 UCX  WARN  unused env variable: UCX_MEM_MMAP_RELOC (set UCX_WARN_UNUSED_ENV_VARS=n to suppress this warning)
[ Info: Initializing ocean_gyre
┌ Info: Starting ocean_gyre
│ dt              = 1.00000e+02
│ timeend         = 2592000.00
│ number of steps = 25920
└ norm(Q)         = 6.9665869213936841e+08
┌ Info: Update
│ simtime = 90600.00 / 2592000.00
│ runtime = 00:00:59
└ norm(Q) = 7.0270326470991445e+08
┌ Info: Update
│ simtime = 215000.00 / 2592000.00
│ runtime = 00:01:59
└ norm(Q) = 7.0789351226975501e+08
┌ Info: Update
│ simtime = 339000.00 / 2592000.00
│ runtime = 00:03:00
└ norm(Q) = NaN

[cnh@TARTARUS CLIMA]$ julia --project=@. foo-cpu.jl 
[TARTARUS.MIT.EDU:27085] mca_base_component_repository_open: unable to open mca_btl_uct: /usr/lib64/openmpi/lib/openmpi/mca_btl_uct.so: undefined symbol: uct_release_md_resource_list (ignored)
[1584671843.709311] [TARTARUS:27085:0]         parser.c:1600 UCX  WARN  unused env variable: UCX_MEM_MMAP_RELOC (set UCX_WARN_UNUSED_ENV_VARS=n to suppress this warning)
[ Info: Initializing ocean_gyre
┌ Info: Starting ocean_gyre
│ dt              = 1.00000e+02
│ timeend         = 2592000.00
│ number of steps = 25920
└ norm(Q)         = 6.9665869213933420e+08
┌ Info: Update
│ simtime =  1400.00 / 2592000.00
│ runtime = 00:01:02
└ norm(Q) = 6.9675608628497553e+08
┌ Info: Update
│ simtime =  3100.00 / 2592000.00
│ runtime = 00:02:05
└ norm(Q) = 6.9689613501868594e+08
┌ Info: Update
│ simtime =  4800.00 / 2592000.00
│ runtime = 00:03:06
└ norm(Q) = 6.9704189442832410e+08
┌ Info: Update
│ simtime =  6500.00 / 2592000.00
│ runtime = 00:04:08
└ norm(Q) = 6.9718788767513585e+08
┌ Info: Update
│ simtime =  8200.00 / 2592000.00
│ runtime = 00:05:10
└ norm(Q) = 6.9733318861567092e+08
┌ Info: Update
│ simtime =  9900.00 / 2592000.00
│ runtime = 00:06:12
└ norm(Q) = 6.9747744194268930e+08

[cnh@TARTARUS CLIMA]$ cat foo.jl
using CLIMA, CLIMA.HydrostaticBoussinesq, CLIMA.Mesh.Filters

using CLIMA.PlanetParameters: grav
using CLIMA.Mesh.Grids: polynomialorder

CLIMA.init()

# Box dimensions in meters
Lx = 4e6    # m 
Ly = 4e6    # m 
Lz = 1000   # m 

# Problem resolution
polynomial_order = 4 
number_of_elements = (20, 20, 20) 

problem = OceanGyre{Float64}(Lx, Ly, Lz) 
model = HydrostaticBoussinesqModel{Float64}(problem, cʰ=sqrt(grav*Lz))
driver_configuration = CLIMA.OceanBoxGCMConfiguration("ocean_gyre", polynomial_order, 
                                                      number_of_elements, model)

grid = driver_configuration.grid

nodal_filters = (vert_filter = CutoffFilter(grid, polynomialorder(grid) - 1), 
                  exp_filter = ExponentialFilter(grid, 1, 8)) 

day = 86400.0
start_time = 0.0 
stop_time = 30day

solver_configuration = CLIMA.setup_solver(start_time, stop_time, driver_configuration,
        init_on_cpu = true,
             ode_dt = 1e2, 
    ode_solver_type = CLIMA.IMEXSolverType(linear_model=LinearHBModel),
          modeldata = nodal_filters
)

result = CLIMA.invoke!(solver_configuration)

[jm-c]

I forgot to mention earlier: all my short test were done on CPU.
And re-run (also on CPU) current master branch: ocean_gyre.jl with explicit setting of ode_dt to 120.
and able to finish 1.day without NaN.
┌ Info: Finished
│ norm(Q) = 7.0322019181584179e+08
│ norm(Q) / norm(Q₀) = 1.0094185283991741e+00
└ norm(Q) - norm(Q₀) = 6.5614996764464378e+06
A bit unfortunate that this standard experiment does not produce any output file.

[jm-c]

Sorry, I made terrible mistake in previous comments but fixed now.

[glwagner]

@christophernhill oh, good point.

Thanks @jm-c. Yes, we need to be able to analyze the solution to diagnose the issue.

As far as I can tell this example is trivial, with initial conditions determined by a function in src/Ocean/HydrostaticBoussinesq/SimpleBoxProblem.jl:

function ocean_init_state!(p::OceanGyre, Q, A, coords, t)
  @inbounds z = coords[3]
  @inbounds H = p.H

  Q.u = @SVector [0,0]
  Q.η = 0
  Q.θ = 9 + 8z/H

  return nothing
end

Edit: actually, there appear to be inhomogeneous boundary conditions:

@inline velocity_flux(p::OceanGyre, y, ρ) = - (p.τₒ / ρ) * cos(y * π / p.Lʸ)

@inline function temperature_flux(p::OceanGyre, y, θ)
  θʳ =  p.θᴱ * (1 - y / p.Lʸ)
  return p.λʳ * (θʳ - θ)
end

[blallen]

We fixed this!