From 0e3a926be4e2caab36023c7a4122b956ef13763d Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Mon, 13 Nov 2023 10:08:19 -0700 Subject: [PATCH] Clean up examples dir --- examples/analyze_neverworld_simulation.jl | 43 --- examples/animate_neverworld_simulation.jl | 152 ----------- examples/compute_mixed_layer_depth.jl | 27 -- examples/melting_baroclinicity.jl | 244 ------------------ .../quarter_degree_near_global_simulation.jl | 96 ------- examples/quarter_degree_omip_run.jl | 36 --- examples/run_neverworld_simulation.jl | 135 ---------- ...ed_perfect_one_degree_model_calibration.jl | 0 .../gm_one_degree_model_calibration.jl | 0 ...spect_one_degree_near_global_simulation.jl | 0 .../one_degree_near_global_simulation.jl | 0 .../perfect_one_degree_model_calibration.jl | 0 12 files changed, 733 deletions(-) delete mode 100644 examples/analyze_neverworld_simulation.jl delete mode 100644 examples/animate_neverworld_simulation.jl delete mode 100644 examples/compute_mixed_layer_depth.jl delete mode 100644 examples/melting_baroclinicity.jl delete mode 100644 examples/quarter_degree_near_global_simulation.jl delete mode 100644 examples/quarter_degree_omip_run.jl delete mode 100644 examples/run_neverworld_simulation.jl rename {examples => experiments/one_degree_simulations}/distributed_perfect_one_degree_model_calibration.jl (100%) rename {examples => experiments/one_degree_simulations}/gm_one_degree_model_calibration.jl (100%) rename {examples => experiments/one_degree_simulations}/inspect_one_degree_near_global_simulation.jl (100%) rename {examples => experiments/one_degree_simulations}/one_degree_near_global_simulation.jl (100%) rename {examples => experiments/one_degree_simulations}/perfect_one_degree_model_calibration.jl (100%) diff --git a/examples/analyze_neverworld_simulation.jl b/examples/analyze_neverworld_simulation.jl deleted file mode 100644 index f76a0035..00000000 --- a/examples/analyze_neverworld_simulation.jl +++ /dev/null @@ -1,43 +0,0 @@ -using Oceananigans -using Oceananigans.Units -using Printf -using JLD2 - -dir = "archive" -backend = OnDisk() -bt = FieldTimeSeries("$dir/neverworld_xyz.jld2", "b"; backend) -t = bt.times -Nt = length(t) -grid = bt.grid - -Nyears = 20 -t_end = 200 * 360 * day # 200 "years" -t_start = (200 - Nyears) * 360 * day # 200 "years" - -t = bt.times -n_stop = searchsortedfirst(t, t_stop) -n_start = searchsortedfirst(t, t_start) - 1 - -for season = 1:4 - n_average_start = n_start + season - 1 - n_average = n_average_start:4:n_stop - - for name in ("b", "u", "v", "w", "e", "κᶜ") - ψt = FieldTimeSeries("$dir/neverworld_xyz.jld2", name; backend) - LX, LY, LZ = location(ψt) - ψ_avg = Field{LX, LY, LZ}(grid) - - for n in n_average - @info name season n - ψn = ψt[n] - ψ_avg .+= ψn / length(n_average) - end - - filename = string("climatology_", name, "_season_", season, ".jld2") - file = jldopen(filename, "a+") - file[name] = ψ_avg - file["season"] = season - close(file) - end -end - diff --git a/examples/animate_neverworld_simulation.jl b/examples/animate_neverworld_simulation.jl deleted file mode 100644 index 822ff075..00000000 --- a/examples/animate_neverworld_simulation.jl +++ /dev/null @@ -1,152 +0,0 @@ -using Oceananigans -using Oceananigans.Units -using GLMakie -using Printf -using Statistics - -# Some plotting parameters -ζlim = 6e-5 - -# Plot a limited number of years (here, 10) -Nyears = 3 -t_end = 200 * 360 * day # 200 "years" -t_start = (200 - Nyears) * 360 * day # 200 "years" - -# Load OnDisk time-series without loading massive amount of data -backend = OnDisk() -bxyt = FieldTimeSeries("neverworld_xy.jld2", "b"; backend) - -# Determine which iterations to load into memory -t̃ = bxyt.times -n_end = searchsortedfirst(t̃, t_end) -n_start = searchsortedfirst(t̃, t_start) - 1 - -t = times = t̃[n_start:n_end] -Nt = length(t) - -bxyt = FieldTimeSeries("neverworld_xy.jld2", "b"; times) -uxyt = FieldTimeSeries("neverworld_xy.jld2", "u"; times) -vxyt = FieldTimeSeries("neverworld_xy.jld2", "v"; times) -exyt = FieldTimeSeries("neverworld_xy.jld2", "e"; times) -κxyt = FieldTimeSeries("neverworld_xy.jld2", "κᶜ"; times) - -byzt = FieldTimeSeries("neverworld_zonal_average.jld2", "b"; times) -eyzt = FieldTimeSeries("neverworld_zonal_average.jld2", "e"; times) -κyzt = FieldTimeSeries("neverworld_zonal_average.jld2", "κᶜ"; times) - -# Hack to set immersed regions to NaN -for ft in (bxyt, uxyt, exyt, κxyt, byzt, eyzt, κyzt) - fp = parent(ft) - fp[fp .== 0] .= NaN -end - -fig = Figure(resolution=(2000, 1800)) - -axbxy = Axis(fig[2, 1], xlabel="Longitude", ylabel="Latitude", title="b(x, y)") -axζxy = Axis(fig[2, 2], xlabel="Longitude", ylabel="Latitude", title="ζ(x, y)") -axexy = Axis(fig[2, 3], xlabel="Longitude", ylabel="Latitude", title="e(x, y)") - -axNyz = Axis(fig[3, 1], xlabel="Longitude", ylabel="z (m)", title="N²(y, z)") -axeyz = Axis(fig[3, 2], xlabel="Longitude", ylabel="z (m)", title="e(y, z)") -axκyz = Axis(fig[3, 3], xlabel="Longitude", ylabel="z (m)", title="κᶜ(y, z)") - -slider = Slider(fig[5, 1:3], range=1:Nt, startvalue=Nt) -n = slider.value - -title = @lift begin - t_day = t[$n] / day - t_year = floor(Int, t_day / 360) - yearday = t_day % 360 - return @sprintf("CATKE Neverworld at % 3d years, % 3d days", t_year, yearday) -end - -Label(fig[0, 1:3], title, fontsize=36) - -using Oceananigans.Operators: ζ₃ᶠᶠᶜ -grid = bxyt.grid -grid = grid.underlying_grid -Nx, Ny, Nz = size(grid) -uxy = Field{Face, Center, Center}(grid, indices=(:, :, 1)) -vxy = Field{Center, Face, Center}(grid, indices=(:, :, 1)) -ζop = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid, uxy, vxy) -ζxy = Field(ζop, indices=(:, :, 1)) - -ζxyn = @lift begin - uxyn = uxyt[$n] - vxyn = vxyt[$n] - interior(uxy) .= interior(uxyn) - interior(vxy) .= interior(vxyn) - compute!(ζxy) - return interior(ζxy, :, :, 1) -end - -byz = Field{Nothing, Center, Center}(grid) -N²yz = Field(∂z(byz)) - -N²yzn = @lift begin - byzn = byzt[$n] - interior(byz) .= interior(byzn) - compute!(N²yz) - return interior(N²yz, 1, :, :) -end - -bxyn = @lift interior(bxyt[$n], :, :, 1) -exyn = @lift interior(exyt[$n], :, :, 1) - -byzn = @lift interior(byzt[$n], 1, :, :) -eyzn = @lift interior(eyzt[$n], 1, :, :) -κyzn = @lift interior(κyzt[$n], 1, :, :) - -κlims = @lift begin - κ_mean = mean(filter(!isnan, interior(κyzt[$n]))) - (κ_mean/2, 3κ_mean) -end - -eyzlims = @lift begin - e_mean = mean(filter(!isnan, interior(eyzt[$n]))) - (e_mean/2, 3e_mean) -end - -exylims = @lift begin - e_mean = mean(filter(!isnan, interior(eyzt[$n]))) - (e_mean/10, e_mean) -end - -x, y, z = nodes(bxyt) -xf, yf, zf = nodes(grid, Face(), Face(), Face()) -cbkw = (vertical=false, flipaxis=true, tellwidth=false) - -hm = heatmap!(axbxy, x, y, bxyn, colorrange=(-0.05, 0.05)) -Colorbar(fig[1, 1], hm; label="Buoyancy", cbkw...) - -hm = heatmap!(axζxy, x, y, ζxyn, colormap=:balance, colorrange=(-ζlim, ζlim)) -Colorbar(fig[1, 2], hm; label="Relative vertical vorticity (s⁻¹)", cbkw...) - -hm = heatmap!(axexy, x, y, exyn, colormap=:solar, colorrange=exylims) -Colorbar(fig[1, 3], hm; label="Turbulent kinetic energy (m² s⁻²)", cbkw...) - -cbkw = (vertical=false, flipaxis=false, tellwidth=false) - -hm = heatmap!(axNyz, y, zf, N²yzn, nan_color=:gray, colorrange=(1e-6, 1e-4)) -Colorbar(fig[4, 1], hm; label="Buoyancy gradient (s⁻²)", cbkw...) -contour!(axNyz, y, z, byzn, levels=30, color=:white, linewidth=3) - -hm = heatmap!(axeyz, y, z, eyzn, nan_color=:gray, colormap=:solar, colorrange=eyzlims) -Colorbar(fig[4, 2], hm; label="Turbulent kinetic energy (m² s⁻²)", cbkw...) -contour!(axeyz, y, z, byzn, levels=30, color=:white, linewidth=3) - -hm = heatmap!(axκyz, y, zf, κyzn, nan_color=:gray, colormap=:solar, colorrange=κlims) -Colorbar(fig[4, 3], hm; label="Tracer eddy diffusivity (m² s⁻¹)", cbkw...) -contour!(axκyz, y, z, byzn, levels=30, color=:white, linewidth=3) - -ylims!(axNyz, -4000, 0) -ylims!(axeyz, -4000, 0) -ylims!(axκyz, -4000, 0) - -display(fig) - -record(fig, "catke_neverworld.mp4", 1:Nt, framerate=8) do nn - @info "Recording frame $nn of $Nt..." - n[] = nn -end - diff --git a/examples/compute_mixed_layer_depth.jl b/examples/compute_mixed_layer_depth.jl deleted file mode 100644 index 07defabf..00000000 --- a/examples/compute_mixed_layer_depth.jl +++ /dev/null @@ -1,27 +0,0 @@ -using ClimaOcean.NearGlobalSimulations: one_degree_near_global_simulation -using ClimaOcean.Diagnostics -using Oceananigans -using Oceananigans - -using GLMakie - -# Build the simulation -simulation = one_degree_near_global_simulation(CPU()) -model = simulation.model - -h = MixedLayerDepthField(model) -compute!(h) -x, y, z = nodes(model.tracers.T) - -fig = Figure(resolution=(1800, 1200)) -ax = Axis(fig[1, 1], xlabel="Longitude", ylabel="Latitude") -hm = heatmap!(ax, x, y, interior(h, :, :, 1), colorrange=(-1, 100)) - -Colorbar(fig[1, 2], hm; - vertical = true, - tellheight = false, - flipaxis = true, - label = "Mixed layer depth (m)") - -display(fig) - diff --git a/examples/melting_baroclinicity.jl b/examples/melting_baroclinicity.jl deleted file mode 100644 index fda05588..00000000 --- a/examples/melting_baroclinicity.jl +++ /dev/null @@ -1,244 +0,0 @@ -using Oceananigans -using Oceananigans.Architectures: arch_array -using Oceananigans.Fields: ZeroField, ConstantField -using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity -using Oceananigans.Units -using Oceananigans.Utils: prettysummary - -using SeawaterPolynomials: TEOS10EquationOfState, thermal_expansion, haline_contraction - -using ClimaSeaIce -using ClimaSeaIce: melting_temperature -using ClimaSeaIce.ThermalBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium - -using Printf -using GLMakie -using Statistics - -include("ice_ocean_model.jl") - -arch = GPU() -Nx = Ny = 256 -Nz = 10 -Lz = 400 -x = y = (-50kilometers, 50kilometers) -halo = (4, 4, 4) -topology = (Periodic, Bounded, Bounded) - -ice_grid = RectilinearGrid(arch; x, y, - size = (Nx, Ny), - topology = (topology[1], topology[2], Flat), - halo = halo[1:2]) - -ocean_grid = RectilinearGrid(arch; topology, halo, x, y, - size = (Nx, Ny, Nz), - z = (-Lz, 0)) - -# Top boundary conditions: -# - outgoing radiative fluxes emitted from surface -# - incoming shortwave radiation starting after 40 days - -ice_ocean_heat_flux = Field{Center, Center, Nothing}(ice_grid) -top_ocean_heat_flux = Qᵀ = Field{Center, Center, Nothing}(ice_grid) -top_salt_flux = Qˢ = Field{Center, Center, Nothing}(ice_grid) -# top_salt_flux = Qˢ = arch_array(arch, zeros(Nx, Ny)) - -# Generate a zero-dimensional grid for a single column slab model - -boundary_conditions = (T = FieldBoundaryConditions(top=FluxBoundaryCondition(Qᵀ)), - S = FieldBoundaryConditions(top=FluxBoundaryCondition(Qˢ))) - -equation_of_state = TEOS10EquationOfState() -buoyancy = SeawaterBuoyancy(; equation_of_state) -horizontal_biharmonic_diffusivity = HorizontalScalarBiharmonicDiffusivity(κ=5e6) - -ocean_model = HydrostaticFreeSurfaceModel(; buoyancy, boundary_conditions, - grid = ocean_grid, - momentum_advection = WENO(), - tracer_advection = WENO(), - #closure = (horizontal_biharmonic_diffusivity, CATKEVerticalDiffusivity()), - closure = CATKEVerticalDiffusivity(), - coriolis = FPlane(f=1.4e-4), - tracers = (:T, :S, :e)) - -Nz = size(ocean_grid, 3) -So = ocean_model.tracers.S -ocean_surface_salinity = view(So, :, :, Nz) -bottom_bc = IceWaterThermalEquilibrium(ConstantField(30)) #ocean_surface_salinity) - -u, v, w = ocean_model.velocities -ocean_surface_velocities = (u = view(u, :, :, Nz), #interior(u, :, :, Nz), - v = view(v, :, :, Nz), #interior(v, :, :, Nz), - w = ZeroField()) - -ice_model = SlabSeaIceModel(ice_grid; - velocities = ocean_surface_velocities, - advection = nothing, #WENO(), - ice_consolidation_thickness = 0.05, - ice_salinity = 4, - internal_thermal_flux = ConductiveFlux(conductivity=2), - #top_thermal_flux = ConstantField(-100), # W m⁻² - 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) - -ocean_simulation = Simulation(ocean_model; Δt=20minutes, verbose=false) -ice_simulation = Simulation(ice_model, Δt=20minutes, verbose=false) - -# Initial condition -S₀ = 30 -T₀ = melting_temperature(ice_model.phase_transitions.liquidus, S₀) + 2.0 - -N²S = 1e-6 -β = haline_contraction(T₀, S₀, 0, equation_of_state) -g = ocean_model.buoyancy.model.gravitational_acceleration -dSdz = - g * β * N²S - -uᵢ(x, y, z) = 0.0 -Tᵢ(x, y, z) = T₀ # + 0.1 * randn() -Sᵢ(x, y, z) = S₀ + dSdz * z #+ 0.1 * randn() - -function hᵢ(x, y) - if sqrt(x^2 + y^2) < 20kilometers - #return 1 + 0.1 * rand() - return 2 - else - return 0 - end -end - -set!(ocean_model, u=uᵢ, S=Sᵢ, T=T₀) -set!(ice_model, h=hᵢ) - -coupled_model = IceOceanModel(ice_simulation, ocean_simulation) -coupled_simulation = Simulation(coupled_model, Δt=1minutes, stop_time=20days) - -S = ocean_model.tracers.S -by = - g * β * ∂y(S) - -function progress(sim) - h = sim.model.ice.model.ice_thickness - S = sim.model.ocean.model.tracers.S - T = sim.model.ocean.model.tracers.T - u = sim.model.ocean.model.velocities.u - msg1 = @sprintf("Iter: % 6d, time: % 12s", iteration(sim), prettytime(sim)) - msg2 = @sprintf(", max(h): %.2f", maximum(h)) - msg3 = @sprintf(", min(S): %.2f", minimum(S)) - msg4 = @sprintf(", extrema(T): (%.2f, %.2f)", minimum(T), maximum(T)) - msg5 = @sprintf(", max|∂y b|: %.2e", maximum(abs, by)) - msg6 = @sprintf(", max|u|: %.2e", maximum(abs, u)) - @info msg1 * msg2 * msg3 * msg4 * msg5 * msg6 - return nothing -end - -coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) - -h = ice_model.ice_thickness -T = ocean_model.tracers.T -S = ocean_model.tracers.S -u, v, w = ocean_model.velocities -η = ocean_model.free_surface.η - -ht = [] -Tt = [] -Ft = [] -Qt = [] -St = [] -ut = [] -vt = [] -ηt = [] -ζt = [] -tt = [] - -ζ = Field(∂x(v) - ∂y(u)) - -function saveoutput(sim) - compute!(ζ) - hn = Array(interior(h, :, :, 1)) - Fn = Array(interior(Qˢ, :, :, 1)) - Qn = Array(interior(Qᵀ, :, :, 1)) - Tn = Array(interior(T, :, :, Nz)) - Sn = Array(interior(S, :, :, Nz)) - un = Array(interior(u, :, :, Nz)) - vn = Array(interior(v, :, :, Nz)) - ηn = Array(interior(η, :, :, 1)) - ζn = Array(interior(ζ, :, :, Nz)) - push!(ht, hn) - push!(Ft, Fn) - push!(Qt, Qn) - push!(Tt, Tn) - push!(St, Sn) - push!(ut, un) - push!(vt, vn) - push!(ηt, ηn) - push!(ζt, ζn) - push!(tt, time(sim)) -end - -coupled_simulation.callbacks[:output] = Callback(saveoutput, IterationInterval(10)) - -run!(coupled_simulation) - -##### -##### Viz -##### - -set_theme!(Theme(fontsize=24)) - -x = xnodes(ocean_grid, Center()) -y = ynodes(ocean_grid, Center()) - -fig = Figure(resolution=(2400, 700)) - -axh = Axis(fig[1, 1], xlabel="x (km)", ylabel="y (km)", title="Ice thickness") -axT = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)", title="Ocean surface temperature") -axS = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)", title="Ocean surface salinity") -axZ = Axis(fig[1, 4], xlabel="x (km)", ylabel="y (km)", title="Ocean vorticity") - -Nt = length(tt) -slider = Slider(fig[2, 1:4], range=1:Nt, startvalue=Nt) -n = slider.value - -title = @lift string("Melt-driven baroclinic instability after ", prettytime(tt[$n])) -Label(fig[0, 1:3], title) - -hn = @lift ht[$n] -Fn = @lift Ft[$n] -Tn = @lift Tt[$n] -Sn = @lift St[$n] -un = @lift ut[$n] -vn = @lift vt[$n] -ηn = @lift ηt[$n] -ζn = @lift ζt[$n] -Un = @lift mean(ut[$n], dims=1)[:] - -x = x ./ 1e3 -y = y ./ 1e3 - -Stop = view(S, :, :, Nz) -Smax = maximum(Stop) -Smin = minimum(Stop) - -compute!(ζ) -ζtop = view(ζ, :, :, Nz) -ζmax = maximum(abs, ζtop) -ζlim = 2e-4 #ζmax / 2 - -heatmap!(axh, x, y, hn, colorrange=(0, 1), colormap=:grays) -heatmap!(axT, x, y, Tn, colormap=:thermal) -heatmap!(axS, x, y, Sn, colorrange = (29, 30), colormap=:haline) -heatmap!(axZ, x, y, ζn, colorrange=(-ζlim, ζlim), colormap=:redblue) - -#heatmap!(axZ, x, y, Tn, colormap=:thermal) -#heatmap!(axF, x, y, Fn) - -display(fig) - -#= -record(fig, "salty_baroclinic_ice_cube.mp4", 1:Nt, framerate=48) do nn - @info string(nn) - n[] = nn -end -=# - diff --git a/examples/quarter_degree_near_global_simulation.jl b/examples/quarter_degree_near_global_simulation.jl deleted file mode 100644 index b687792a..00000000 --- a/examples/quarter_degree_near_global_simulation.jl +++ /dev/null @@ -1,96 +0,0 @@ -using ClimaOcean.NearGlobalSimulations: quarter_degree_near_global_simulation - -using Oceananigans -using Oceananigans.Units -using Oceananigans.Utils: WallTimeInterval -using Oceananigans.BuoyancyModels: buoyancy -using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: - MixingLength, TurbulentKineticEnergyEquation, CATKEVerticalDiffusivity - -using ParameterEstimocean.Parameters: closure_with_parameters -using JLD2 - -##### -##### Boundary layer turbulence closure options -##### - -# "Ri-based" --- uses calibrated defaults in Oceananigans -ri_based = RiBasedVerticalDiffusivity() - -# Choose closure -boundary_layer_turbulence_closure = ri_based - -@show boundary_layer_turbulence_closure - -##### -##### Build the simulation -##### - -start_time = 0days -stop_time = start_time + 10years - -simulation = quarter_degree_near_global_simulation(; start_time, stop_time, boundary_layer_turbulence_closure) - -# Define output -slices_save_interval = 1day -fields_save_interval = 30days -Nx, Ny, Nz = size(simulation.model.grid) - -dir = "/storage2/simone/quarter-degree-data/" -closure_name = typeof(boundary_layer_turbulence_closure).name.wrapper -output_prefix = "near_global_$(Nx)_$(Ny)_$(Nz)_$closure_name" - -simulation.output_writers[:checkpointer] = Checkpointer(simulation.model; dir, - prefix = output_prefix * "_checkpointer", - schedule = WallTimeInterval(10minutes), - cleanup = true, - overwrite_existing = true) - -model = simulation.model - -simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers); dir, - schedule = TimeInterval(slices_save_interval), - filename = output_prefix * "_fields", - with_halos = true, - overwrite_existing = true) - -slice_indices = [(:, :, Nz), (:, :, Nz-10)] -output_names = [:surface, :near_surface] -for n = 1:2 - indices = slice_indices[n] - - outputs = Dict() - - for name in keys(model.tracers) - c = model.tracers[name] - outputs[name] = Field(c; indices) - end - - outputs[:u] = Field(model.velocities.u; indices) - outputs[:v] = Field(model.velocities.v; indices) - outputs[:w] = Field(model.velocities.w; indices) - outputs[:η] = model.free_surface.η - outputs[:ζ] = VerticalVorticityField(model.grid, model.velocities; indices) - - name = output_names[n] - simulation.output_writers[name] = JLD2OutputWriter(model, outputs; dir, - schedule = TimeInterval(slices_save_interval), - filename = output_prefix * "_fields_$name", - with_halos = true, - overwrite_existing = true) -end - -@info "Running a simulation with Δt = $(prettytime(simulation.Δt))" - -spinup = 11days -simulation.Δt = 1minute -simulation.stop_time = time(simulation) + spinup -run!(simulation) - -simulation.stop_time = start_time + 2.2years - spinup -simulation.Δt = 10minutes -run!(simulation) - -@info "Simulation took $(prettytime(simulation.run_wall_time))." - diff --git a/examples/quarter_degree_omip_run.jl b/examples/quarter_degree_omip_run.jl deleted file mode 100644 index 8a5afb84..00000000 --- a/examples/quarter_degree_omip_run.jl +++ /dev/null @@ -1,36 +0,0 @@ -using GLMakie -using Oceananigans -using ClimaOcean.Bathymetry: regrid_bathymetry -using ClimaOcean.VerticalGrids: stretched_vertical_faces, PowerLawStretching - -##### -##### Quarter degree near-global grid -##### - -arch = CPU() - -z_faces = stretched_vertical_faces(surface_layer_Δz = 8, - surface_layer_height = 128, - stretching = PowerLawStretching(1.02), - minimum_depth = 6000) - -grid = LatitudeLongitudeGrid(arch; - size = (4 * 360, 4 * 160, 1), - latitude = (-80, 80), - longitude = (-180, 180), - z = z_faces, - halo = (4, 4, 4)) - -h = regrid_bathymetry(grid, height_above_water=1, minimum_depth=5) - -λ, φ, z = nodes(h) - -land = interior(h) .> 0 -interior(h)[land] .= NaN - -fig = Figure(resolution=(2400, 1200)) -ax = Axis(fig[1, 1], xlabel="Longitude", ylabel="Latitude") -heatmap!(ax, λ, φ, interior(h, :, :, 1), nan_color=:white, colorrange=(-5000, 0)) - -display(fig) - diff --git a/examples/run_neverworld_simulation.jl b/examples/run_neverworld_simulation.jl deleted file mode 100644 index 0acddbe3..00000000 --- a/examples/run_neverworld_simulation.jl +++ /dev/null @@ -1,135 +0,0 @@ -using Oceananigans -using Oceananigans.Units -using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity -using Oceananigans.ImmersedBoundaries: PartialCellBottom, GridFittedBottom -using ClimaOcean.IdealizedSimulations: neverworld_simulation -using ClimaOcean.VerticalGrids: stretched_vertical_faces, PowerLawStretching -using Printf -using CUDA - -closure = CATKEVerticalDiffusivity(minimum_turbulent_kinetic_energy = 1e-6, - minimum_convective_buoyancy_flux = 1e-11) - -# closure = RiBasedVerticalDiffusivity() - -z = stretched_vertical_faces(surface_layer_Δz = 8, - surface_layer_height = 64, - stretching = PowerLawStretching(1.02), - maximum_Δz = 400.0, - minimum_depth = 4000) - -simulation = neverworld_simulation(GPU(); z, - #ImmersedBoundaryType = GridFittedBottom, - ImmersedBoundaryType = PartialCellBottom, - horizontal_resolution = 1/8, - longitude = (0, 60), - latitude = (-70, 0), - time_step = 10minutes, - stop_time = 4 * 360days, - closure) - -model = simulation.model -grid = model.grid - -@show grid -@show model - -start_time = Ref(time_ns()) -previous_model_time = Ref(time(simulation)) - -function progress(sim) - b = sim.model.tracers.b - e = sim.model.tracers.e - u, v, w = sim.model.velocities - - msg = @sprintf("Iter: %d, time: %s, extrema(b): (%6.2e, %6.2e)", - iteration(sim), prettytime(sim), minimum(b), maximum(b)) - - msg *= @sprintf(", max(e): %6.2e", maximum(e)) - - msg *= @sprintf(", max|u|: %6.2e, max|w|: %6.2e", - maximum(maximum(abs, q) for q in (u, v, w)), maximum(abs, w)) - - try - κᶜ = sim.model.diffusivity_fields.κᶜ - msg *= @sprintf(", max(κᶜ): %6.2e", maximum(κᶜ)) - catch - end - - elapsed = 1e-9 * (time_ns() - start_time[]) - elapsed_model_time = time(sim) - previous_model_time[] - SYPD = (elapsed_model_time/360days) / (elapsed/day) - - msg *= @sprintf(", wall time: %s, SYPD: %.1f", prettytime(elapsed), SYPD) - start_time[] = time_ns() - previous_model_time[] = time(sim) - - @info msg - - return nothing -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) - -# Set up output -Nx, Ny, Nz = size(grid) -Δt = simulation.Δt -Δt_minutes = round(Int, Δt / minutes) -ib_str = grid.immersed_boundary isa PartialCellBottom ? "partial_cells" : "full_cells" -output_suffix = "$(Nx)_$(Ny)_$(Nz)_dt$(Δt_minutes)_$(ib_str).jld2" -output_dir = "." #/nobackup1/glwagner/" -fine_output_frequency = 1day -i = round(Int, Nx/10) # index for yz-sliced output - -z = znodes(grid, Face(), with_halos=true) - -K = CUDA.@allowscalar [Nz, - searchsortedfirst(z, -100), - searchsortedfirst(z, -400)] - -i = round(Int, Nx/10) # index for yz-sliced output - -diffusivity_fields = (; κᶜ = model.diffusivity_fields.κᶜ) -outputs = merge(model.velocities, model.tracers, diffusivity_fields) -zonally_averaged_outputs = NamedTuple(n => Average(outputs[n], dims=1) for n in keys(outputs)) - -simulation.output_writers[:yz] = JLD2OutputWriter(model, outputs; - schedule = TimeInterval(fine_output_frequency), - filename = joinpath(output_dir, "neverworld_yz_" * output_suffix), - indices = (i, :, :), - with_halos = true, - overwrite_existing = true) - -simulation.output_writers[:zonal] = JLD2OutputWriter(model, zonally_averaged_outputs; - schedule = TimeInterval(fine_output_frequency), - filename = joinpath(output_dir, "neverworld_zonal_average_" * output_suffix), - with_halos = true, - overwrite_existing = true) - -for (n, k) in enumerate(K) - name = Symbol(:xy, n) - simulation.output_writers[name] = JLD2OutputWriter(model, outputs; - schedule = TimeInterval(fine_output_frequency), - filename = joinpath(output_dir, "neverworld_xy$(n)_" * output_suffix), - indices = (:, :, Nz), - with_halos = true, - overwrite_existing = true) -end - -simulation.output_writers[:xyz] = JLD2OutputWriter(model, outputs; - schedule = TimeInterval(90days), - filename = joinpath(output_dir, "neverworld_xyz_" * output_suffix), - with_halos = true, - overwrite_existing = true) - -simulation.output_writers[:checkpointer] = Checkpointer(model, - schedule = TimeInterval(360days), - dir = output_dir, - prefix = "neverworld_$(Nx)_$(Ny)_$(Nz)_checkpoint", - cleanup = true) - -@info "Running..." -@show simulation - -run!(simulation) - diff --git a/examples/distributed_perfect_one_degree_model_calibration.jl b/experiments/one_degree_simulations/distributed_perfect_one_degree_model_calibration.jl similarity index 100% rename from examples/distributed_perfect_one_degree_model_calibration.jl rename to experiments/one_degree_simulations/distributed_perfect_one_degree_model_calibration.jl diff --git a/examples/gm_one_degree_model_calibration.jl b/experiments/one_degree_simulations/gm_one_degree_model_calibration.jl similarity index 100% rename from examples/gm_one_degree_model_calibration.jl rename to experiments/one_degree_simulations/gm_one_degree_model_calibration.jl diff --git a/examples/inspect_one_degree_near_global_simulation.jl b/experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl similarity index 100% rename from examples/inspect_one_degree_near_global_simulation.jl rename to experiments/one_degree_simulations/inspect_one_degree_near_global_simulation.jl diff --git a/examples/one_degree_near_global_simulation.jl b/experiments/one_degree_simulations/one_degree_near_global_simulation.jl similarity index 100% rename from examples/one_degree_near_global_simulation.jl rename to experiments/one_degree_simulations/one_degree_near_global_simulation.jl diff --git a/examples/perfect_one_degree_model_calibration.jl b/experiments/one_degree_simulations/perfect_one_degree_model_calibration.jl similarity index 100% rename from examples/perfect_one_degree_model_calibration.jl rename to experiments/one_degree_simulations/perfect_one_degree_model_calibration.jl