Issue with TimeStepWizard

[Sumanshekhar17]

I am facing some issue with TimeStepWizard. It is saying that there is no method like TimeStepWizard , I have also faced this issue in previous version of Oceananigans but that time after updating Oceananigans resolved it. This time these are the information of version- Oceananigans = v0.61.3 julia = 1.6.1

I am pasting me whole code so that you can reproduce the error-

using Printf
using Plots
using Oceananigans
using Oceananigans.Utils
using Oceananigans.Units: minutes, hour, hours, day
using Oceananigans.Grids: nodes
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters: JLD2OutputWriter, FieldSlicer, TimeInterval
using Oceananigans.Diagnostics: accurate_cell_advection_timescale

#Defining Grid
#number of grid spacing in south,north and vertical direction
const Nx=256
const Ny=256
const Nz=64

#Length of grid in south, north and vertical direction
const Lx=4
const Ly=4
const Lz=0.1

const S = 1.6 # Stretching factor
hyperbolically_spaced_nodes(k) = -Lz-Lz*(tanh(S * ( (-(k-Nz-1) ) / Nz - 1)) / tanh(S))
computational_grid = VerticallyStretchedRectilinearGrid(size = (Nx, Ny, Nz), 
                                                               architecture = CPU(),
                                                               x = (0,Lx),
                                                               y = (0,Ly),
                                                               halo = (3, 3, 3),
                                                               z_faces = hyperbolically_spaced_nodes)

#Governing parameters 
#Coefficient of Thermal expansion
const alpha= 2e-4

#Coefficient of Salinity
const beta=8e-4

#diffusive viscocity
const v=1e-5 
#diffusivity
const k=2e-6 

const f=0   #coriolis parameter

g = 300

#Temperature(Tracer) boundary condition 
Amplitude = 950 # Amplitude (W/m²) of Heat flux

Q_cool = -143.1092 # (W/m²) {due to Longwave + Latent heat + Sensible Heat}

gausian(t) =  exp(-((t)^2)/(0.025*(86400)^2)) 

# These are the times during which peak of the heat flux can be observed
peak1 = 0.56day
peak2 = 1.56day
peak3 = 2.56day
peak4 = 3.56day
peak5 = 4.56day
peak6 = 5.56day

# Heat absorbed due to shortwave
Q_sh(t) = Amplitude*(gausian(t-peak1) + gausian(t-peak2) + gausian(t-peak3) + gausian(t-peak4)+ gausian(t-peak5)+ gausian(t-peak6) )

Q_net(t) = Q_sh(t) + Q_cool

function heatflux(t)
    if Q_net(t) <= 0
        bias = 1
        multiplier = 0
    else
        bias = 0
        multiplier = 1
    end

    return Q = bias*Q_cool
end

dTdz = 0
T_bcs = FieldBoundaryConditions( top = FluxBoundaryCondition(heatflux(t)),
                                 bottom = GradientBoundaryCondition(dTdz))

#Defining Heat Source Term as a forcing 
function ShortWavePenetration(x,y,z,t,params)

    if Q_net(t) <= 0
        bias = 1
        multiplier = 0
    else
        bias = 0
        multiplier = 1
    end

    Q_band1(z) = (params.Ι/params.λ1)*exp(-z/params.λ1) 
    Q_band2(z) = (params.Ι/params.λ2)*exp(-z/params.λ2)

    if  z >= computational_grid.zᵃᵃᶜ[Nz]
        cooling_source = multiplier*Q_cool/(params.ρ*params.Cp*computational_grid.Δzᵃᵃᶜ[Nz])
        return Q_sh(t) * (Q_band1(z) + Q_band2(z))/(params.ρ*params.Cp) + (cooling_source*computational_grid.Δzᵃᵃᶜ[Nz])

    else 
        return Q_sh(t) * (Q_band1(z) + Q_band2(z))/(params.ρ*params.Cp)
    end

end
ShortWavePenetration_parameters = (ρ = 1000,
                                   Cp = 4182,
                                   λ1 = 0.35,
                                   λ2 = 23.0,
                                   Ι = 0.58
                                    )
heat_source_term = Forcing(ShortWavePenetration,
                            parameters = ShortWavePenetration_parameters)

#Velocity Boundary conditions 
const Qᵘ=0
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))
v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))

#Buoyancy that depends on temperature and salinity 
buoyancy = SeawaterBuoyancy(gravitational_acceleration = g,equation_of_state=LinearEquationOfState(α=alpha, β=beta))

#Model instantiation 
using Oceananigans.Advection
using Oceananigans.TurbulenceClosures

model = NonhydrostaticModel(architecture = CPU(),
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            grid = computational_grid,
                            tracers = (:T, :S),
                            coriolis = FPlane(f=f),
                            buoyancy = buoyancy,
                            closure = SmagorinskyLilly(),
                            forcing = (T=heat_source_term,),
                            boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs))

#Setting Initial Conditions 
using CSV
using DataFrames
initial_temperature = Matrix(CSV.read("initial_temperature.csv", DataFrame))
itemp = zeros(Nx, Ny, Nz)
for x ∈ 1:Nx
    for y ∈ 1:Ny
        itemp[x,y,:] = initial_temperature
    end
end

initial_salinity = Matrix(CSV.read("initial_salinity.csv", DataFrame))
isal = zeros(Nx, Ny, Nz)
for x ∈ 1:Nx
    for y ∈ 1:Ny
        isal[x,y,:] = initial_salinity
    end
end

# `set!` the `model` fields using functions or constants:
set!(model, T = itemp, S = isal)

#Setting up a simulation 
using Oceananigans.Diagnostics: accurate_cell_advection_timescale
wizard = TimeStepWizard(cfl=0.5,Δt=0.1, max_change=1.1, max_Δt=1minutes,cell_advection_timescale = accurate_cell_advection_timescale)
start_time = time_ns() # so we can print the total elapsed wall time

# Print a progress message
progress_message(sim) =
    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
            sim.model.clock.iteration, prettytime(model.clock.time),
            prettytime(wizard), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))

simulation = Simulation(model,
                    Δt = wizard,
             stop_time = 6day,
    iteration_interval = 1,
              progress = progress_message
)
#Output
fields = Dict("u" => model.velocities.u,"v" => model.velocities.v,"w" => model.velocities.w, "T" => model.tracers.T)

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filepath="DevangSetup.nc",
                       schedule=TimeInterval(6) )
run!(simulation)

And the error message is following (I am pasting partially) -

MethodError: no method matching zero(::TimeStepWizard{Float64, typeof(accurate_cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\Dates\src\periods.jl:53
  zero(::Colorant) at C:\Users\My Account\.julia\packages\ColorTypes\6m8P7\src\traits.jl:477
  zero(::TaylorSeries.Taylor1) at C:\Users\My Account\.julia\packages\TaylorSeries\tveWm\src\arithmetic.jl:37
  ...

Stacktrace:
 [1] iszero(x::TimeStepWizard{Float64, typeof(accurate_cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
   @ Base .\number.jl:40
 [2] prettytime(t::TimeStepWizard{Float64, typeof(accurate_cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
   @ Oceananigans.Utils C:\Users\My Account\.julia\packages\Oceananigans\To7WB\src\Utils\pretty_time.jl:18

I am attaching the file of Initial data which is used in the code- initial_salinity.csv initial_temperature.csv

[francispoulin]

Are you sure this is exactly the same code that you ran?

I tried executing it as is and found an error on this line

T_bcs = FieldBoundaryConditions( top = FluxBoundaryCondition(heatflux(t)),
                                 bottom = GradientBoundaryCondition(dTdz))

The error is the following:

ERROR: UndefVarError: t not defined
Stacktrace:
 [1] top-level scope
   @ REPL[38]:1
 [2] top-level scope
   @ ~/.julia/packages/CUDA/fRSUT/src/initialization.jl:52

[glwagner]

@Sumanshekhar17 can you write up a shorter example that does not involve data? Hopefully, as short as possible to reproduce the error. The code you've posted is too long to efficiently debug.

[Sumanshekhar17]

Are you sure this is exactly the same code that you ran?

I tried executing it as is and found an error on this line

T_bcs = FieldBoundaryConditions( top = FluxBoundaryCondition(heatflux(t)),
                                 bottom = GradientBoundaryCondition(dTdz))

The error is the following:

ERROR: UndefVarError: t not defined
Stacktrace:
 [1] top-level scope
   @ REPL[38]:1
 [2] top-level scope
   @ ~/.julia/packages/CUDA/fRSUT/src/initialization.jl:52

@francispoulin sorry, I did one mistake there. Though I am pasting a short code here which can produce the same error -

using Printf
using Oceananigans
using Oceananigans.Units: minutes, hour, hours, day

const Nx=256
const Ny=256
const Nz=64

#Length of grid in south, north and vertical direction
const Lx=4
const Ly=4
const Lz=0.1

const S = 1.6 # Stretching factor
hyperbolically_spaced_nodes(k) = -Lz-Lz*(tanh(S * ( (-(k-Nz-1) ) / Nz - 1)) / tanh(S))
computational_grid = VerticallyStretchedRectilinearGrid(size = (Nx, Ny, Nz), 
                                                               architecture = CPU(),
                                                               x = (0,Lx),
                                                               y = (0,Ly),
                                                               halo = (3, 3, 3),
                                                               z_faces = hyperbolically_spaced_nodes)

#Coefficient of Thermal expansion
const alpha= 2e-4

#Coefficient of Salinity
const beta=8e-4
#heat flux
const Bo=3.6e-4
#coriolis parameter
const f=-0.5
#gravitational acceleration
const g=300

Q = Bo
const dTdz = 0 # K m⁻¹
T_bcs = FieldBoundaryConditions(
                                 top = FluxBoundaryCondition(Q),
                                 bottom = GradientBoundaryCondition(dTdz))

const Qᵘ=0
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))
v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))

buoyancy = SeawaterBuoyancy(gravitational_acceleration = g,equation_of_state=LinearEquationOfState(α=alpha, β=beta))

using Oceananigans.Advection
using Oceananigans.TurbulenceClosures

model = NonhydrostaticModel(architecture = CPU(),
                            advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            grid = computational_grid,
                            tracers = (:T, :S),
                            coriolis = FPlane(f=f),
                            buoyancy = buoyancy,
                            closure = SmagorinskyLilly(),
                            boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs))
using Oceananigans.Diagnostics: accurate_cell_advection_timescale
wizard = TimeStepWizard(cfl=0.5,Δt=0.1, max_change=1.1, max_Δt=1minutes,cell_advection_timescale = accurate_cell_advection_timescale)
start_time = time_ns() # so we can print the total elapsed wall time

# Print a progress message
progress_message(sim) =
    @printf("i: %04d, t: %s, Δt: %s, wmax = %.1e ms⁻¹, wall time: %s\n",
            sim.model.clock.iteration, prettytime(model.clock.time),
            prettytime(wizard), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))

simulation = Simulation(model,
                    Δt = wizard,
             stop_time = 12minutes,
    iteration_interval = 1,
              progress = progress_message
)
fields = Dict("u" => model.velocities.u,"v" => model.velocities.v,"w" => model.velocities.w, "T" => model.tracers.T)

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filepath="TestSetup.nc",
                       schedule=TimeInterval(6) )
run!(simulation)

[glwagner]

I think the problem is that prettytime(wizard) does not work.

MWE:

julia> using Oceananigans

julia> wizard = TimeStepWizard(cfl=0.1, Δt=1.0)
TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)}(0.1, Inf, 2.0, 0.5, Inf, 0.0, 1.0, Oceananigans.Utils.cell_advection_timescale, Oceananigans.Simulations.infinite_diffusion_timescale)

julia> prettytime(wizard)
ERROR: MethodError: no method matching zero(::TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:53
  zero(::ArrayInterface.NDIndex{N, I} where I<:Tuple{Vararg{Any, N}}) where N at /Users/gregorywagner/.julia/packages/ArrayInterface/VFy81/src/ndindex.jl:95
  zero(::LinearAlgebra.UniformScaling{T}) where T at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/uniformscaling.jl:136
  ...
Stacktrace:
 [1] iszero(x::TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
   @ Base ./number.jl:40
 [2] prettytime(t::TimeStepWizard{Float64, typeof(Oceananigans.Utils.cell_advection_timescale), typeof(Oceananigans.Simulations.infinite_diffusion_timescale)})
   @ Oceananigans.Utils ~/.julia/packages/Oceananigans/MRn50/src/Utils/pretty_time.jl:18
 [3] top-level scope
   @ REPL[3]:1

I'm not sure what prettytime(wizard) should return --- I think this is expected.

[francispoulin]

Thanks @Sumanshekhar17 for the simpler example and @glwagner has answered your question I believe.

Why don't you try it without prettytime(wizard) to see if that works?

[Sumanshekhar17]

It is strange It was working in the earlier versions.

[Sumanshekhar17]

@francispoulin and @glwagner Thank you for this. before closing this issue I wanted to ask weather I have written my boundary condtion as a function heatflux is correct or not -

function heatflux(x, y, z, t)
    if Q_net(t) <= 0
        bias = 1
        multiplier = 0
    else
        bias = 0
        multiplier = 1
    end

    return bias*Q_cool
end
Temp_bc = FluxBoundaryCondition(heatflux)
dTdz = 0
T_bcs = FieldBoundaryConditions( top = Temp_bc,
                                 bottom = GradientBoundaryCondition(dTdz))

I am constantly getting error stating that -

TaskFailedException

    nested task error: MethodError: no method matching heatflux(::Float64, ::Float64, ::Float64)
    Closest candidates are:
      heatflux(::Any, ::Any, ::Any, ::Any) at In[9]:1
    Stacktrace:
      [1] call
        @ C:\Users\My Account\.julia\packages\Cassette\r4kKQ\src\context.jl:456 [inlined]
      [2] fallback
        @ C:\Users\My Account\.julia\packages\Cassette\r4kKQ\src\context.jl:454 [inlined]
      [3] _overdub_fallback(::Any, ::Vararg{Any, N} where N)
        @ C:\Users\My Account\.julia\packages\Cassette\r4kKQ\src\overdub.jl:585 [inlined]
      [4] overdub
        @ C:\Users\My Account\.julia\packages\Cassette\r4kKQ\src\overdub.jl:585 [inlined]
      [5] overdub
        @ C:\Users\My Account\.julia\packages\Oceananigans\To7WB\src\BoundaryConditions\continuous_boundary_function.jl:122 [inlined]

[glwagner]

heatflux(x, y, t) should be a function of x, y, t only.

From the docs:

https://clima.github.io/OceananigansDocumentation/stable/model_setup/boundary_conditions/#.-Spatially-and-temporally-varying-flux

[tomchor]

@Sumanshekhar17 Has this helped? Can we close this issue?

[Sumanshekhar17]

Ohh Yess It helped me Thank you everyone .