NaN error

[Sumanshekhar17]

Respected sir,
I want to simulate a small prototype of Open Ocean convection and did all necessary scaling according to that. Please have a look at my code, I don't know why my code isn't working. It is giving NaN error.

using Printf
using Plots
using JLD2

using Oceananigans
using Oceananigans.Utils

using Oceananigans.Grids: nodes
using Oceananigans.Diagnostics: FieldMaximum
using Oceananigans.OutputWriters: JLD2OutputWriter, FieldSlicer, TimeInterval

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

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

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))

#Coefficient of Thermal expansion
alpha= 2e-4

#Coefficient of Salinity
saline=0


Ra=1.8e8 #rayleigh number

Ro=0.24  #Horizontal Rossby number

v=1e-6 #diffusive viscocity

k=1e-6 #diffusivity 

l=Lx/2  #center of gaussian field

m=Ly/2  #center of gausian field


Bo=9e-5 #maximum surface flux

f=0.5386   #coriolis parameter

beta=4.039e-3

#Gaussian Surface forcing
heat_flux(x,y,t) = Bo*exp(-0.5*((x-l)^2+(y-m)^2))

dTdz = 0 # K m⁻¹

T_bcs = TracerBoundaryConditions(grid,
                                 top = BoundaryCondition(Flux, heat_flux),
                                 bottom = BoundaryCondition(Gradient, dTdz))


Qᵘ=0    #Zero flux boundary condition at the top surface

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ), bottom=ValueBoundaryCondition(0.0))
v_bcs = VVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ), bottom=ValueBoundaryCondition(0.0))
w_bcs = WVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(0.0), bottom=ValueBoundaryCondition(0.0))

buoyancy = SeawaterBuoyancy(gravitational_acceleration = 300,equation_of_state=LinearEquationOfState(α=alpha, β=saline))


#Incompressible model initiation 

using Oceananigans.Advection
using Oceananigans.TurbulenceClosures

model = IncompressibleModel(architecture = CPU(),

                            timestepper = :RungeKutta3,
                            grid = grid,
                            coriolis = BetaPlane(f₀=f, β=beta),
                            buoyancy = buoyancy,
                            closure = SmagorinskyLilly(),
                            boundary_conditions = (u=u_bcs, v=v_bcs, w=w_bcs, T=T_bcs))


# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

# Temperature initial condition: a stable density gradient with random noise superposed.
dTz = 81.2762

Tᵢ(x, y, z) = 30 + dTz * z + dTz * model.grid.Lz * 1e-6 * Ξ(z)

# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-8 * Ξ(z)

# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, v=uᵢ, w=uᵢ, T=Tᵢ)

wizard = TimeStepWizard(cfl=1.5, Δt=0.001, max_change=1.1, max_Δt=0.2minute)

wmax = FieldMaximum(abs, model.velocities.w)

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.Δt), wmax(sim.model),
            prettytime((time_ns() - start_time) * 1e-9))


simulation = Simulation(model, Δt=wizard, stop_time=15minutes, iteration_interval=10,
                        progress=progress_message)



simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                           prefix = "modified_open_ocean_convection2",
                         schedule = TimeInterval(0.2minute),
                            force = true)

run!(simulation)````

[glwagner]

NaNs often occur either because the time step is too small, or because the problem being posed is too extreme or not physical (such that a reasonable and stable time step cannot be found).

One issue is that the use of ValueBoundaryCondition is incorrect here:

w_bcs = WVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(0.0), bottom=ValueBoundaryCondition(0.0))

You are trying to specify an impenetrable boundary condition on the vertical velocity, which requires NormalFlowBoundaryCondition:

w_bcs = WVelocityBoundaryConditions(grid, top = NormalFlowBoundaryCondition(0.0), bottom=NormalFlowBoundaryCondition(0.0))

this is also the default for a grid with topology(grid) = (Periodic, Periodic, Bounded), so there's no need to specify it explicitly. In other words, try changing

boundary_conditions = (u=u_bcs, v=v_bcs, w=w_bcs, T=T_bcs))

to

boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs))

in your constructor for IncompressibleModel.

A good approach for developing a complex simulation is to start with something simple (or something that worked previously) and modify the script incrementally (slowly!), re-running and re-testing the script after each significant change. For this is useful to run either a two-dimensional version of the problem or a very low resolution simulation. I also suggest adding some lines at the bottom of your script to plot the solution so you can visually inspect it after a run.

Your setup looks like a laboratory-scale simulation of convection. I think starting with uniform cooling (rather than spatially varying), omitting the effect of salinity, and using an IsotropicDiffusivity in two-dimensions (rather than starting out in three dimensions) might be a good path forward for working and building up your script.

[Sumanshekhar17]

sir, this code worked when I omitted the effect of rotation, But after adding the rotation effect it started creating problems.

[ali-ramadhan]

@Sumanshekhar17 We are more than happy to help with debugging Oceananigans.jl scripts and simulation setups, but you need to show some effort and willingness too. We cannot set up your simulation for you.

@glwagner provided some excellent advice on how you can develop a complex simulation, I suggest you follow it.

This is how we develop our simulations, slowly and iteratively. Most importantly: we check what it is doing by plotting the solution and making movies. Plotting and debugging simulations in 2D is easier.

If you can post some figures and movies of your solution we might be able to help.

[ali-ramadhan]

@Sumanshekhar17 Sorry if I was a little rude. We want to help but we need more information.

How did you come up with your values for f and beta? Maybe they're too large?

Making a movie of the solution without rotation then with rotation could be very helpful.

I'll turn this issue into a GitHub discussion so we can continue discussing there.

[Sumanshekhar17]

@ali-ramadhan sir, I have scaled down this value from a simulation originally simulated on 200(km)x200(km)x4(km) domain; it was from a research paper, and I am using Oceananigan to test it's scaled problem so that I can become comfortable with this solver. I initially tested by turning off the Coriolis effect. It worked -

open_ocean_convection.mp4

And sir, I am extremely sorry for not posting the complete code earlier. Here is the code after implementing the Coriolis effect-

using Printf
using Plots
using JLD2

using Oceananigans
using Oceananigans.Utils

using Oceananigans.Grids: nodes
using Oceananigans.Diagnostics: FieldMaximum
using Oceananigans.OutputWriters: JLD2OutputWriter, FieldSlicer, TimeInterval

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

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

grid = RegularCartesianGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))

#Coefficient of Thermal expansion
alpha= 2e-4

#Coefficient of Salinity
saline=0


Ra=1.8e8 #rayleigh number

Ro=0.24  #Horizontal Rossby number

v=1e-6 #diffusive viscocity

k=1e-6 #diffusivity 

l=Lx/2  # centre of the gaussian field

m=Ly/2  # centre of the gaussian field


Bo=9e-5 #maximum surface flux

f=0.5386   #coriolis parameter

beta=4.039e-3

#Surface forcing
heat_flux(x,y,t) = Bo*exp(-0.5*((x-l)^2+(y-m)^2))

dTdz = 0 # K m⁻¹

T_bcs = TracerBoundaryConditions(grid,
                                 top = BoundaryCondition(Flux, heat_flux),
                                 bottom = BoundaryCondition(Gradient, dTdz))


Qᵘ=0    #Zero flux boundary condition at the top surface

u_bcs = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ), bottom=ValueBoundaryCondition(0.0))
v_bcs = VVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ), bottom=ValueBoundaryCondition(0.0))
# w_bcs = WVelocityBoundaryConditions(grid, top = ValueBoundaryCondition(0.0), bottom=ValueBoundaryCondition(0.0))

buoyancy = SeawaterBuoyancy(gravitational_acceleration = 300,equation_of_state=LinearEquationOfState(α=alpha, β=saline))


#Incompressible model initiation 

using Oceananigans.Advection
using Oceananigans.TurbulenceClosures

model = IncompressibleModel(architecture = CPU(),
#                              advection = UpwindBiasedFifthOrder(),
                            timestepper = :RungeKutta3,
                            grid = grid,
                            coriolis = BetaPlane(f=f, β=beta),
                            buoyancy = buoyancy,
                            closure = SmagorinskyLilly(),
                            boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs))


# Random noise damped at top and bottom
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

# Temperature initial condition: a stable density gradient with random noise superposed.
dTz = 81.2762

Tᵢ(x, y, z) = 30 + dTz * z + dTz * model.grid.Lz * 1e-6 * Ξ(z)

# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-8 * Ξ(z)

# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, v=uᵢ, w=uᵢ, T=Tᵢ)

wizard = TimeStepWizard(cfl=0.1, Δt=0.2, max_change=1.1, max_Δt=0.2minute)

wmax = FieldMaximum(abs, model.velocities.w)

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.Δt), wmax(sim.model),
            prettytime((time_ns() - start_time) * 1e-9))


simulation = Simulation(model, Δt=wizard, stop_time=15minutes, iteration_interval=10,
                        progress=progress_message)



simulation.output_writers[:fields] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                           prefix = "modified_open_ocean_convection2",
                         schedule = TimeInterval(0.2minute),
                            force = true)

run!(simulation)





# load("open_ocean_convection1.jld2")
file = jldopen("modified_open_ocean_convection.jld2")

# u_snapshot = file["timeseries/w/0"][:, :,60]
# v_snapshot = file["timeseries/v/0"][:, :, 60]
# speed_snapshot = sqrt.(u_snapshot.*u_snapshot + v_snapshot.*v_snapshot)

# Coordinate arrays
xC, yC, zC = file["grid/xC"][1:256],file["grid/yC"][1:256],file["grid/zC"][1:65]
Lx, Ly, Lz = file["grid/Lx"],file["grid/Ly"],file["grid/Lz"]

# Extract a vector of iterations
iterations = parse.(Int, keys(file["timeseries/t"]))


@info "Making a neat movie of verticle velocity and Temperature..."

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]
    u_snapshot = file["timeseries/u/$iteration"][:, :, 60]
    v_snapshot = file["timeseries/v/$iteration"][:, :, 60]
    w_snapshot = file["timeseries/w/$iteration"][:, 128, :]
    speed_snapshot = sqrt.(u_snapshot.*u_snapshot + v_snapshot.*v_snapshot)
  #  T_snapshot = file["timeseries/T/$iteration"][:, 128, :]
   

    ulims = 5.0e-4
    
    ulevels = range(-ulims, stop=ulims, length=50)
    
    slims = 0.004
    
    slevels = range(0, stop=slims, length=50)
    
    Tlims = 2.5e3
  
    Tlevels = range(0, stop=Tlims, length=50)

    kwargs1 = (xlabel="x", ylabel="z", aspectratio=50, linewidth=0, colorbar=true,
              xlims=(0, Lx), ylims=(-Lz,0))
    
    kwargs2 = (xlabel="x", ylabel="y", aspectratio=1, linewidth=0, colorbar=true,
              xlims=(0, Lx), ylims=(0, Ly))

    u_plot = contourf(xC, zC, clamp.(w_snapshot', -ulims, ulims);
                        color = :balance,
                        clims=(-ulims, ulims), levels=ulevels,
                        kwargs1...)

    s_plot = contourf(xC, yC, clamp.(speed_snapshot', -slims,slims );
                       color = :thermal,
                       clims=(0, slims), levels=slevels,
                       kwargs2...)
#title
    u_title = @sprintf("vertical velocity (m s⁻¹), t = %s", prettytime(t))
    s_title = "Horizontal speed"
    plot( u_plot, s_plot, title= [u_title s_title], layout=(1, 2), size=(1200, 600))
end

gif(anim, "./modified_open_ocean_convection.gif", fps = 1) # hide````

I will check once again the parameters though I have checked them many times.  Please forgive me If I made any mistake or disturbed the decorum of this community, I am learning.
Thank you!

[Sumanshekhar17]

Dear @glwagner and @ali-ramadhan sir should I upload my calculations also where I have scaled down the values?

[Sumanshekhar17]

dear sir,
my code ran. Thanks, @ali-ramadhan sir for highlighting the beta term. Actually, the beta plane approximation is not appropriate for periodic boundary condition and for such a small domain. by taking Fplane approximation it worked.

successful_animation.mp4

above is the animation of the same.

[ali-ramadhan]

Ah sorry I missed your reply yesterday. Indeed a beta-plane might not make sense with periodic boundary conditions!

Glad you were able to visualize the solution. Not sure if this is what you were expecting but the movies look pretty cool. I'm guessing the vertical velocity seems to form columns since the fluid is unstratified?

Also, what advection scheme are you using? I can see some numerical artifacts in the right movie of the horizontal speed. You can try passing a higher-order advection scheme to the IncompressibleModel constructor like

model = IncompressibleModel(..., advection = UpwindBiasedFifthOrder(), ...)

which should improve the fidelity of the solution. The best advection scheme is probably WENO5() but it's pretty expensive on the CPU so I usually only use it on the GPU. UpwindBiasedFifthOrder() is pretty good to use on the CPU.

EDIT: Ah I see you have UpwindBiasedFifthOrder() in your latest script but it seems to be commented out.

[glwagner]

In fact, we probably want to throw an error if someone tries to use BetaPlane when y is Periodic...