differences between windows and Mac results

[ajpeliz]

Hi,
Same code, same versions getting different results in different platforms. Actually it runs in the MacOS and blows up in windows.
What should I look at? has anyone experienced something similar?
thanks

[simone-silvestri]

Hello @ajpeliz. There should be no difference between machines. Can you post a snippet of code that reproduces your issue (as minimal as possible) so we can test it?

[glwagner]

Is this with Julia 1.11?

[ajpeliz]

This ia a MSc Geophysis project I am developing with students. The aim is a minimal model of the North Atlantic and we are starting with very simple conditions. I developed my code based on your "hands on workshop". It runs Ok on my Mac laptop, but blows top on students computers (on Windows).
SmallAtlantic2_jl.txt

[ajpeliz]

1.11.4 No dia 15/04/2025, às 17:11, Gregory L. Wagner ***@***.***> escreveu: Is this with Julia 1.11? — Reply to this email directly, view it on GitHub<#4402 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/A5L6CIQSBSW24GKJTPUAJD32ZUVTJAVCNFSM6AAAAAB3F6Y3PKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTEOBUGQZDGNQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[simone-silvestri]

Ok, we have had some consistency problems with julia 1.11.4, can you try with julia 1.10.9 would it blow up?

[ajpeliz]

We finally could test with julia 1.10.9, and yes it runs!
Thank you.

[glwagner]

Hi @ajpeliz if you can please copy/paste your code into here, we can use it to diagnose the problem and hopefully solve it. We have to be able to run the code ourselves in order to do this work though.

[glwagner]

can you please format your code in markdown? ChatGPT can explain how to use triple backticks, as an example: https://chatgpt.com/share/68007b8f-8da8-8002-b039-d16eae9a41ca

[glwagner]

its still formatted as plain text rather than code

[ajpeliz]

using Oceananigans
using GLMakie
using Oceananigans.Units
#
# TO DO: 
# test a single gyre in f-plane
# test a single gyre in beta-plane
# test a double gyre in beta-plane
#

#Grid size
Nx=101  # lines
Ny=60   # columns
Nz=10
Lz=1000
lons=-Nx:0;
lats=0:Ny;

# model time
Δt=12minutes
stop_time=60days


# Set z-grid streatching 
# grids may be equally spaced on the vertical 
# or we can have more resolution at the surface 
# ana uneven dsitribution of the vertical grid is named 
# grid streatching
#
#
# Streach thhe grid!
refinement = 4 # controls spacing near surface (higher means finer spaced)
stretching = 12  # controls rate of stretching at bottom
# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz
# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement
# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
# Generating function of the cell interfaces  (depth of the surfaces deviding the cells)
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

# SET MODEL GRID
grid=LatitudeLongitudeGrid(size=(Nx,Ny,Nz), 
                           longitude = (-100,0),
                           latitude = (0,60),
                           z = z_faces, 
                           topology = (Bounded,Bounded,Bounded))


# deafult linear eq (based on a thermal expansion and salinity contractions coefs
buoyancy = SeawaterBuoyancy()
#
#  MIXING !!!
#
# Set horizontal and vertical viscosity (ν) and diffusivity (κ)
Az=HorizontalScalarDiffusivity(ν=1e4,κ=1e3)
Kz=VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(),ν=0.01,κ=1e-3)
# simple vertical mixing to force convection (hydrostatic will not do convection on its own)
convective_adjustment=ConvectiveAdjustmentVerticalDiffusivity(convective_κz= 1.0)
closure=(Az,Kz,convective_adjustment)
#closure=(Az,Kz)

# In case of boundary conditions check Simone Silvester hands on session at 1h:25min
# we can load fields from a .jld2 file with the same size of the grid
# imagine we need wind stress tau:
# fileOBC=jldopen("nameofthefile.jld2");
# tau=fileOBC["tau"]

# Alternatively we need to make tau field here!
u₁₀ = 10                             # m s⁻¹, average wind velocity 10 meters above the ocean
cᴰ = 2.5e-3                          # dimensionless drag coefficient
ρₒ = 1026.0                          # kg m⁻³, average density at the surface of the world ocean
ρₐ = 1.225                           # kg m⁻³, average density of air at sea-level
τx = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²
#τy =   ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) # m² s⁻²

# for the moment we go:
tau=τx*[cos(j*π/Ny) for i=1:Nx,j=1:Ny]
#tau=τx*[cos(j*2*π/Ny) for i=1:Nx,j=1:Ny]


# Imposing the Wind Stress
# we do it via a upper boundary condition
u_top_bc = FluxBoundaryCondition(-tau)
v_top_bc = FluxBoundaryCondition(tau*0.0)

# Condition for surface heat flux
# we can impose a static value or a restoring
# See silvester 1h:33 hands on workshop
# here we start with a simple value
sst=ones(Nx,Ny)*20.0 
T_top_bc = ValueBoundaryCondition(sst)

u_bcs=FieldBoundaryConditions(top=u_top_bc)
v_bcs=FieldBoundaryConditions(top=v_top_bc)
T_bcs=FieldBoundaryConditions(top=T_top_bc)
boundary_conditions=(u = u_bcs, v = v_bcs, T = T_bcs)

# Set model
model=HydrostaticFreeSurfaceModel(; grid,
		        		closure,
                                       #coriolis = FPlane(f=1e-4),
                                        coriolis=HydrostaticSphericalCoriolis(),
       		        		boundary_conditions,
		        		buoyancy,
		        		tracers = (:T, :S))

# initial fields;
Tini=ones(Nx,Ny,Nz);
for k=1:Nz 
 Tini[:,:,k].=k;
end

#set!(model, T = 15.0, S=35.0)
set!(model, T = Tini, S=35.0)

#
# define the simulation
#
#



simulation = Simulation(model; Δt, stop_time)

progress(sim) = @info "Time : $(prettytime(sim.model.clock.time)), 
                       max(|u|): $(maximum(abs, sim.model.velocities.u)), 
                       max(|v|): $(maximum(abs, sim.model.velocities.v)), 
                       max(T) : $(maximum(model.tracers.T))"

simulation.callbacks[:progress] =Callback(progress, IterationInterval(100))

run!(simulation)


#
# Plots and figures
#


println(" ")
println(" ")
println(" Sim end ")
println(" ")
println(" ")
println(" Plotting !!! ")
println(" ")
println(" ")

# Extract surface fields
uu=Array(interior(model.velocities.u, :,:,grid.Nz));
vv=Array(interior(model.velocities.v, :,:,grid.Nz));
tt=Array(interior(model.tracers.T, :,:,grid.Nz));

# make Makie fig
fig = Figure(size = (500, 900));
#def axis
axis_kwargs = (xlabel="longitude",
               ylabel="latitude",
               aspect = AxisAspect(grid.Lx/grid.Ly),
               limits = ((-grid.Lx,0), (0, grid.Ly)));

# colormaping ranges (same as caxis)
taulims=(-0.001,0.001);
ulims  =(-0.5,0.5);
tlims  =(17,20);

# logica da divisão dos eixos ax=Axis(fig[2,1]) duas linhs uma coluna
ax_tau  = Axis(fig[1, 1]; title = "τₓ", axis_kwargs...);
ax_u    = Axis(fig[2, 1]; title = "Surface u", axis_kwargs...);
ax_v    = Axis(fig[3, 1]; title = "Surface v", axis_kwargs...);
ax_t    = Axis(fig[4, 1]; title = "Surface T", axis_kwargs...);

hm_tau = heatmap!(ax_tau, lons,lats,tau; colormap = :balance, colorrange = taulims);
Colorbar(fig[1, 2], hm_tau; label = "N m⁻2")

hm_uu = heatmap!(ax_u, lons,lats,uu; colormap = :balance, colorrange = ulims);
Colorbar(fig[2, 2], hm_uu; label = "m s⁻1")

hm_vv = heatmap!(ax_v, lons,lats,vv; colormap = :balance, colorrange = ulims);
Colorbar(fig[3, 2], hm_vv; label = "m s⁻1")

hm_tt = heatmap!(ax_t, lons,lats,tt; colormap = :balance, colorrange = tlims);
Colorbar(fig[4, 2], hm_tt; label = "C")

# compute depth integrated flow (a vertical heighted sum)
# from model velocities
u3d=Array(interior(model.velocities.u, :,:,:));
v3d=Array(interior(model.velocities.v, :,:,:));
Hz=ones(Nx,Ny,Nz); # will be the vertical size of cells)
U=ones(Nx,Ny,Nz);
V=ones(Nx,Ny,Nz);
Ui=ones(Nx,Ny);
Vi=ones(Nx,Ny);

for i in 1:Nx
 for j in 1:Ny
  for k in 1:Nz
   Hz[i,j,k]=abs(z_faces(k)-z_faces(k+1))
   U[i,j,k]=u3d[i,j,k]*Hz[i,j,k]
   V[i,j,k]=v3d[i,j,k]*Hz[i,j,k]
  end
 Ui[i,j]=sum(U[i,j,:])
 Vi[i,j]=sum(V[i,j,:])
 end
end


# make Makie fig
fig2= Figure(size = (500, 500));
#def axis
axis_kwargs = (xlabel="longitude",
               ylabel="latitude",
               aspect = AxisAspect(grid.Lx/grid.Ly),
               limits = ((-grid.Lx,0), (0, grid.Ly)));
# colormaping ranges (same as caxis)
ulims  =(-200,200);
# logica da divisão dos eixos ax=Axis(fig[2,1]) duas linhs uma coluna
ax_U  = Axis(fig2[1, 1]; title = "Ui", axis_kwargs...);
ax_V  = Axis(fig2[2, 1]; title = "Vi", axis_kwargs...);
hm_U = heatmap!(ax_U, lons,lats,Ui; colormap = :balance, colorrange = ulims);
Colorbar(fig2[1, 2], hm_U; label = "m2 /s")
hm_V = heatmap!(ax_V, lons,lats,Vi; colormap = :balance, colorrange = ulims);
Colorbar(fig2[2, 2], hm_V; label = "m2 s⁻1")

tmp=reverse(Ui, dims=2);   # We integrate Ui in the y (latitude direction)
tmp=cumsum(tmp, dims=2);
Us=reverse(tmp, dims=2);


# plotting Stream function of depth integrated flow map with projection
using GeoMakie

figgeo = Figure()
ax = GeoAxis(figgeo[1,1],
      source = "+proj=longlat +datum=WGS84",
      dest = "+proj=lcc +lon_0=-50 +lat_1=0 +lat_2=45")#,
    # limits=(-110,10),(-10,70))
sp = contour!(ax, lons[1:end-1], lats[1:end-1], Us, levels=15)
#ax.xticks = Makie.WilkinsonTicks(5; k_min = 3)
#ax.yticks = Makie.WilkinsonTicks(5; k_min = 3)

# display figures 
window = display(GLMakie.Screen(), fig)
window = display(GLMakie.Screen(), fig2)
window = display(GLMakie.Screen(), figgeo)

[glwagner]

Sorry to harp on this. To clarify, you need to use triple backticks: (```). I will do this for you. IF you want to learn more, I made another chatgpt explanation: https://chatgpt.com/share/6801af43-92a4-8002-9a47-667b3a98a954