glwagner
commented
3 years ago A few other users have run into this bug too! I think @ali-ramadhan has a solution but the PR hasn't been merged yet. @mukund-gupta may know more.
Closed Sumanshekhar17 closed 2 years ago
glwagner
commented
3 years ago A few other users have run into this bug too! I think @ali-ramadhan has a solution but the PR hasn't been merged yet. @mukund-gupta may know more.
mukund-gupta
commented
3 years ago Hey! Yes, I got the same error. Working with ali/unclog-docs branch got rid of that problem for me, but I guess it would be better to wait for the PR to be merged.
Sumanshekhar17
commented
3 years ago @mukund-gupta, thank you for the information. In addition to that, please point me to the discussion where you raised that issue. It could be worth reading. And I am new to git-hub, so I face difficulty in understanding ali/unclog-docs branch
glwagner
commented
3 years ago @Sumanshekhar17 to use the code in ali/unclog-docs, you can use the package manager:
using Pkg
pkg"add Oceananigans#ali/unclog-docs"This will reinstall the version of Oceananigans.jl contained on the branch ali/unclog-docs.
The documentation for the julia package manager may also be useful:
Sumanshekhar17
commented
3 years ago @glwagner, I got rid of this UnderVarError but now I am facing domain error while applying Chebyshev spacing in z.
But at the same time, when I apply any linear profile, It works.
linear function like -(2k-1/2Nz)
I get the error message after applying the set() function.-
julia> # `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, v=uᵢ, w=uᵢ, T=Tᵢ)
ERROR: TaskFailedException:
DomainError with -1.1554673348527535e-7:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.If we see the error message, it says to replace x^y, but there is no arithmetic like that in my code.
glwagner
commented
3 years ago Can you paste the whole script?
Sumanshekhar17
commented
3 years ago #############################################################################
using Printf
using Pkg
pkg"add Oceananigans#ali/unclog-docs"
using Oceananigans
using Oceananigans.Utils
using Oceananigans.Units: minute, minutes, hour
using Oceananigans.Grids: nodes
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters: JLD2OutputWriter, FieldSlicer, TimeInterval
using Statistics
# using CUDA
#number of grid spacing in south, north and vertical direction and in oceananigans
#they call this as the size of one grid in that direction
const Nx=256
const Ny=256
const Nz=65
#Length of grid in south, north and vertical direction
const Lx=4
const Ly=4
const Lz=0.1
#Vetrical temperature gradient
const dTz = 70
#scaled gravitational acceleration
const g=300
const R0=1
const T0=30
const Factor_T =1e-6
const Factor_V=1e-8
const CFL=0.6
const maximum_Δt= 0.2minute
const max_allowable_change_in_Δt= 1.1
const initial_timestep = 0.05
Name_of_simulation = "ocean_convection_Fplane_Vertical_strecthed_grid"
# stretching_function= cos(π / 2 * (z - 1)) , won't work as z is only defined after grid function.
const stop_time_info= 12minute #Stop the simulation once this much model clock time has passed.
const iteration_interval_info= 10
#How often to update the time step, check stop criteria, and call progress function (in number of iterations).
#Coefficient of Thermal expansion
const alpha= 2e-4
#Coefficient of Salinity
const saline=0
#diffusive viscocity
const v=1e-5
#diffusivity
const k=2e-6
const l=Lx/2 #center of gaussian field
const m=Ly/2 #center of gausian field
const Bo=3.6e-4 #maximum surface flux
const f=-0.5 #coriolis parameter
#Constructing the Vertically stretched grid with Chebyshev spacing in z
computational_grid = VerticallyStretchedRectilinearGrid(size = (Nx, Ny, Nz),
x = (0, Lx),
y = (0, Ly),
halo = (3, 3, 3),
zF = k -> cos((2k - 1) / 2Nz))
# computational_grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
show(computational_grid)
##xC,yC,zC are not defined for vertically stretched grid
computational_grid.zᵃᵃᶜ
computational_grid.xᶜᵃᵃ
computational_grid.yᵃᶜᵃ
###############################################################
#checking the value of all the governing parameters.
Raf= (Bo*Lz^4)/(v*k^2) #Raleigh number
Ro= ((Bo/(-f)^3)^0.5)/Lz #Rossby number in vertical direction
Ro_l= (Bo/(-f)^3*R0^2)^(1/4) #Rossby number in horizontal direction
N= (g*alpha*dTz)^0.5 #vertical stratification
N_cap= N/(-f) #normalised vertical stratification
################################################################
#no Beta plane approximation for periodic boundary condition
#surface forcing function for top boundary condition
Q(x,y,t) = Bo*exp((-0.5)*(((x-l)/0.4)^2+((y-m)/0.4)^2))
##################################################################
const dTdz = 0 # K m⁻¹
T_bcs = TracerBoundaryConditions(computational_grid,
top = FluxBoundaryCondition(Q),
bottom = GradientBoundaryCondition(dTdz))
const Qᵘ=0 #Zero flux boundary condition at the top surface
u_bcs = UVelocityBoundaryConditions(computational_grid, top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))
v_bcs = VVelocityBoundaryConditions(computational_grid, top = FluxBoundaryCondition(Qᵘ), bottom=ValueBoundaryCondition(0.0))
buoyancy = SeawaterBuoyancy(gravitational_acceleration = g,equation_of_state=LinearEquationOfState(α=alpha, β=saline))
#Incompressible model initiation
using Oceananigans.Advection
using Oceananigans.TurbulenceClosures
model = IncompressibleModel(architecture = CPU(),
advection = UpwindBiasedFifthOrder(),
timestepper = :RungeKutta3,
grid = computational_grid,
coriolis = FPlane(f=f),
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.
Tᵢ(x, y, z) = T0 + dTz * z + dTz * model.grid.Lz * Factor_T * Ξ(z)
# Velocity initial condition: random noise scaled by the friction velocity.
uᵢ(x, y, z) = sqrt(abs(Qᵘ)) * Factor_V * Ξ(z)
# `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, v=uᵢ, w=uᵢ, T=Tᵢ)
wizard = TimeStepWizard(cfl= CFL, Δt = initial_timestep, max_change = max_allowable_change_in_Δt, max_Δt = maximum_Δt)
# A type for calculating adaptive time steps based on capping the CFL number at `cfl`.
# On calling `update_Δt!(wizard, model)`, the `TimeStepWizard` computes a time-step such that
# ``cfl = max(u/Δx, v/Δy, w/Δz) Δt``, where ``max(u/Δx, v/Δy, w/Δz)`` is the maximum ratio
# between model velocity and along-velocity grid spacing anywhere on the model grid. The new
# `Δt` is constrained to change by a multiplicative factor no more than `max_change` or no
# less than `min_change` from the previous `Δt`, and to be no greater in absolute magnitude
# than `max_Δt` and no less than `min_Δt`.
# wmax = FieldMaximum(abs, model.velocities.w) ##not working due to update, It has been renamed into other func.
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), maximum(abs, sim.model.velocities.w),
prettytime((time_ns() - start_time) * 1e-9))
simulation = Simulation(model, Δt=wizard, stop_time= stop_time_info, iteration_interval = iteration_interval_info,
progress=progress_message)
simulation.output_writers[:fields] =
JLD2OutputWriter(model, merge(model.velocities, model.tracers),
prefix = Name_of_simulation,
schedule = TimeInterval(0.2minute),
force = true)
run!(simulation)
using JLD2
using Plots
# load("modified_open_ocean_convection_Fplane_GPU.jld2")
file = jldopen("ocean_convection_Fplane_Vertical_strecthed_grid.jld2")
# 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 = 0.0025
ulevels = range(-ulims, stop=ulims, length=50)
slims = 0.025
slevels = range(0, stop=slims, length=50)
Tlims = 2.5e3
Tlevels = range(0, stop=Tlims, length=70)
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, "./ocean_convection_Fplane_Vertical_strecthed_grid.gif", fps = 1) #animation functionI think this is resolved. Please re-open if not.
Hi all, I have tried the Vertical stretching function in my code, and I tried Constructing the Vertically stretched grid with Chebyshev spacing in the z direction.-
But when I set up the model, I got an UndefVarError stating this-
And I think that this is because of the Vertical stretching function, but I couldn't get what is Id here? I haven't used Id anywhere. Am I heading in the right direction?