Regional polar simulation with prototype IceOceanModel

.gitignore

@@ -29,3 +29,7 @@ docs/site/
 *.swp
 *.svg
 *.gif
+*.tif
+*.zip
+*.nc
+*.mat

experiments/README.md

@@ -0,0 +1,21 @@
+# Experiments to ClimaOcean
+
+This repository holds ClimaOcean "experiments".
+An experiment has a few components:
+
+1. Script that builds and runs an `Oceananigans.Simulation`. More specifically, a "script"
+   is a file (such as `script.jl`) that can be run at the command line via
+
+```bash
+$ julia --project script.jl`,
+```
+
+or from the Julia REPL via
+
+```julia-repl
+julia> include("script.jl")
+```
+
+2. An "environment" - aa script that builds and runs a `Simulation`
+
+

experiments/antarctic_circumpolar_simulation/LocalPreferences.toml

@@ -0,0 +1 @@
+[CUDA_Runtime_jll]

experiments/antarctic_circumpolar_simulation/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
+
+[extras]
+CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"

experiments/antarctic_circumpolar_simulation/analyze_antarctic_circumpolar_simulation.jl

@@ -0,0 +1,59 @@
+using Oceananigans
+using GLMakie
+
+filename = "antarctic_circumpolar_surface.jld2"
+
+ut = FieldTimeSeries(filename, "u")
+Tt = FieldTimeSeries(filename, "T")
+St = FieldTimeSeries(filename, "S")
+
+times = Tt.times
+Nt = length(times)
+
+function set_zero_to_NaN!(u)
+    up = parent(u)
+    up[up .== 0] .= NaN
+    return nothing
+end
+
+for n = 1:Nt
+    un = ut[n]
+    Tn = Tt[n]
+    Sn = St[n]
+    set_zero_to_NaN!(Tn)
+    set_zero_to_NaN!(Sn)
+    if n > 1
+        set_zero_to_NaN!(un)
+    end
+end
+
+set_theme!(Theme(fontsize=32))
+fig = Figure(resolution=(2600, 2000))
+
+axT = Axis(fig[1, 1], xlabel="Longitude (degrees)", ylabel="Latitude (degrees)")
+axS = Axis(fig[2, 1], xlabel="Longitude (degrees)", ylabel="Latitude (degrees)")
+axu = Axis(fig[3, 1], xlabel="Longitude (degrees)", ylabel="Latitude (degrees)")
+
+slider = Slider(fig[4, 1], range=1:Nt, startvalue=1)
+n = slider.value
+
+Tn = @lift interior(Tt[$n], :, :, 1)
+Sn = @lift interior(St[$n], :, :, 1)
+un = @lift interior(ut[$n], :, :, 1)
+
+λc, φc, zc = nodes(Tt)
+λu, φc, zc = nodes(ut)
+
+hm = heatmap!(axT, λc, φc, Tn, colorrange=(-1, 30), nan_color=:gray, colormap=:thermal)
+Colorbar(fig[1, 2], hm, label="Temperature (ᵒC)")
+
+hm = heatmap!(axS, λc, φc, Sn, colorrange=(31, 35), nan_color=:gray, colormap=:haline)
+Colorbar(fig[2, 2], hm, label="Salinity (psu)")
+
+hm = heatmap!(axu, λu, φc, un, colorrange=(-0.2, 0.2),  nan_color=:gray, colormap=:redblue)
+Colorbar(fig[3, 2], hm, label="Zonal velocity (m s⁻¹)")
+
+
+
+display(fig)
+

experiments/antarctic_circumpolar_simulation/antarctic_circumpolar_ice_ocean_simulation.jl

@@ -0,0 +1,131 @@
+using Oceananigans
+using Oceananigans.Units
+using Oceananigans.Coriolis: ActiveCellEnstrophyConservingScheme
+using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
+using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
+using JLD2
+using Printf
+
+include("ice_ocean_model.jl")
+
+
+
+filename = "antarctic_circumpolar_grid_initial_conditions.jld2"
+file = jldopen(filename)
+Tᵢ = Array{Float64, 3}(file["T"])
+Sᵢ = Array{Float64, 3}(file["S"])
+Tᵢ = reverse(Tᵢ, dims=3)
+Sᵢ = reverse(Sᵢ, dims=3)
+longitude = file["longitude"]
+latitude = file["latitude"]
+z = file["z"]
+zb = file["bottom_height"]
+close(file)
+
+arch = GPU()
+Nx = 6 * 360
+Ny = length(latitude) - 1
+Nz = length(z) - 1
+
+grid = LatitudeLongitudeGrid(arch; longitude, latitude, z,
+                             size = (Nx, Ny, Nz),
+                             halo = (7, 7, 7),
+                             topology = (Periodic, Bounded, Bounded))
+
+@show grid
+
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(zb))
+
+equation_of_state = TEOS10EquationOfState()
+
+model = HydrostaticFreeSurfaceModel(; grid,
+                                    tracers = (:T, :S, :e),
+                                    buoyancy = SeawaterBuoyancy(; equation_of_state),
+                                    coriolis = HydrostaticSphericalCoriolis(scheme = ActiveCellEnstrophyConservingScheme()),
+                                    #coriolis = HydrostaticSphericalCoriolis(),
+                                    free_surface = SplitExplicitFreeSurface(; grid, cfl=0.5),
+                                    #momentum_advection = VectorInvariant(),
+                                    momentum_advection = VectorInvariant(vorticity_scheme = WENO(order=9),
+                                                                         divergence_scheme = WENO(order=9),
+                                                                         vertical_scheme = WENO(order=9)),
+                                    tracer_advection = WENO(order=7),
+                                    #tracer_advection = WENO(),
+                                    closure = CATKEVerticalDiffusivity())
+
+@show model
+
+set!(model, T=Tᵢ, S=Sᵢ)
+
+ocean_simulation = Simulation(model, Δt=1minutes, stop_time=30days)
+
+start_time = Ref(time_ns())
+
+function progress(sim)
+    elapsed = 1e-9 * (time_ns() - start_time[])
+
+    msg1 = string("Iter: ", iteration(sim),
+                  ", time: ", prettytime(sim),
+                  ", wall time: ", prettytime(elapsed))
+
+    u, v, w = sim.model.velocities
+    msg2 = @sprintf(", max|u|: (%.2e, %.2e, %.2e) m s⁻¹",
+                    maximum(abs, u),
+                    maximum(abs, v),
+                    maximum(abs, w))
+    @info msg1 * msg2
+
+    start_time[] = time_ns()
+
+    return nothing
+end
+
+ocean_simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
+
+output_dir = "." #/nobackup1/glwagner/"
+outputs = merge(model.velocities, model.tracers)
+Nz = size(grid, 3)
+ocean_simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs,
+                                                          schedule = TimeInterval(1day),
+                                                          indices = (:, :, Nz),
+                                                          dir = output_dir,
+                                                          filename = "antarctic_circumpolar_surface.jld2",
+                                                          overwrite_existing = true)
+
+
+#####
+##### Ice simulation
+#####
+
+ice_grid = LatitudeLongitudeGrid(arch; longitude, latitude,
+                                 size = (Nx, Ny),
+                                 halo = halo_size(grid),
+                                 topology = (Periodic, Bounded, Flat))
+
+land = zb .>= 0
+ice_grid = ImmersedBoundaryGrid(ice_grid, GridFittedBoundary(land))
+
+Nz = size(grid, 3)
+So = model.tracers.S
+ocean_surface_salinity = view(So, :, :, Nz)
+bottom_bc = IceWaterThermalEquilibrium(ConstantField(30)) #ocean_surface_salinity)
+
+ice_ocean_heat_flux = Field{Center, Center, Nothing}(ice_grid)
+
+ice_model = SlabSeaIceModel(ice_grid;
+                            velocities = nothing,
+                            advection = nothing, #WENO(),
+                            ice_consolidation_thickness = 0.05,
+                            ice_salinity = 4,
+                            internal_thermal_flux = ConductiveFlux(conductivity=2),
+                            top_thermal_flux = ConstantField(0), # W m⁻²
+                            top_thermal_boundary_condition = PrescribedTemperature(0),
+                            bottom_thermal_boundary_condition = bottom_bc,
+                            bottom_thermal_flux = ice_ocean_heat_flux)
+
+ice_simulation = Simulation(ice_model, Δt=20minutes, verbose=false)
+
+coupled_model = IceOceanModel(ice_simulation, ocean_simulation)
+coupled_simulation = Simulation(coupled_model, Δt=1minutes, stop_time=30days)
+
+run!(coupled_simulation)
+

experiments/antarctic_circumpolar_simulation/antarctic_circumpolar_simulation.jl

@@ -0,0 +1,89 @@
+using Oceananigans
+using Oceananigans.Units
+using Oceananigans.Coriolis: ActiveCellEnstrophyConservingScheme
+using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
+using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
+using JLD2
+using Printf
+
+filename = "antarctic_circumpolar_grid_initial_conditions.jld2"
+file = jldopen(filename)
+Tᵢ = Array{Float64, 3}(file["T"])
+Sᵢ = Array{Float64, 3}(file["S"])
+Tᵢ = reverse(Tᵢ, dims=3)
+Sᵢ = reverse(Sᵢ, dims=3)
+longitude = file["longitude"]
+latitude = file["latitude"]
+z = file["z"]
+zb = file["bottom_height"]
+close(file)
+
+arch = GPU()
+Nx = 6 * 360
+Ny = length(latitude) - 1
+Nz = length(z) - 1
+halo = (5, 5, 5)
+
+grid = LatitudeLongitudeGrid(arch; longitude, latitude, z, halo,
+                             size = (Nx, Ny, Nz),
+                             topology = (Periodic, Bounded, Bounded))
+
+@show grid
+
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(zb))
+
+equation_of_state = TEOS10EquationOfState()
+
+model = HydrostaticFreeSurfaceModel(; grid,
+                                    tracers = (:T, :S, :e),
+                                    buoyancy = SeawaterBuoyancy(; equation_of_state),
+                                    coriolis = HydrostaticSphericalCoriolis(scheme = ActiveCellEnstrophyConservingScheme()),
+                                    free_surface = SplitExplicitFreeSurface(; grid, cfl=0.5),
+                                    momentum_advection = VectorInvariant(vorticity_scheme = WENO(),
+                                                                         divergence_scheme = WENO(),
+                                                                         vertical_scheme = WENO()),
+                                    tracer_advection = WENO(),
+                                    closure = CATKEVerticalDiffusivity())
+
+@show model
+
+set!(model, T=Tᵢ, S=Sᵢ)
+
+simulation = Simulation(model, Δt=1minutes, stop_time=30days)
+
+start_time = Ref(time_ns())
+
+function progress(sim)
+    elapsed = 1e-9 * (time_ns() - start_time[])
+    start_time[] = time_ns()
+
+    msg1 = string("Iter: ", iteration(sim),
+                  ", time: ", prettytime(sim),
+                  ", wall time: ", prettytime(elapsed))
+
+    u, v, w = sim.model.velocities
+
+    msg2 = @sprintf(", max|u|: (%.2e, %.2e, %.2e) m s⁻¹",
+                    maximum(abs, parent(u)),
+                    maximum(abs, parent(v)),
+                    maximum(abs, parent(w)))
+
+    @info msg1 * msg2
+
+    return nothing
+end
+
+simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
+
+output_dir = "." #/nobackup1/glwagner/"
+outputs = merge(model.velocities, model.tracers)
+Nz = size(grid, 3)
+simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs,
+                                                    schedule = TimeInterval(1day),
+                                                    indices = (:, :, Nz),
+                                                    dir = output_dir,
+                                                    filename = "antarctic_circumpolar_surface.jld2",
+                                                    overwrite_existing = true)
+
+run!(simulation)
+