diff --git a/validation/near_global_lat_lon/analyze_bathymetry_forcing.jl b/validation/near_global_lat_lon/analyze_bathymetry_forcing.jl deleted file mode 100644 index 7998cfd683..0000000000 --- a/validation/near_global_lat_lon/analyze_bathymetry_forcing.jl +++ /dev/null @@ -1,52 +0,0 @@ -using JLD2 -using Printf -using GLMakie - -Nx = 128 -Ny = 60 - -bathymetry_path = "bathy_128x60var4.bin" -east_west_stress_path = "off_TAUXvar1.bin" -north_south_stress_path = "off_TAUY.bin" -sea_surface_temperature_path="sst25_128x60x12.bin" - -Nbytes = sizeof(Float32) * Nx * Ny -bathymetry = reshape(bswap.(reinterpret(Float32, read(bathymetry_path, Nbytes))), (Nx, Ny)) - -τˣ = reshape(bswap.(reinterpret(Float32, read(east_west_stress_path, 12Nbytes))), (Nx, Ny, 12)) -τʸ = reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, 12Nbytes))), (Nx, Ny, 12)) -target_sea_surface_temperature = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, 12Nbytes))), (Nx, Ny, 12)) - -τˣ = τˣ[:, :, 1] -τʸ = τʸ[:, :, 1] -target_sea_surface_temperature = target_sea_surface_temperature[:, :, 1] - -# bathymetry = Array{Float64, 2}(bathymetry) -# τˣ = Array{Float64, 2}(τˣ) -# τʸ = Array{Float64, 2}(τʸ) -# target_sea_surface_temperature = Array{Float64, 2}(target_sea_surface_temperature) - -bathymetry_file = jldopen("earth_bathymetry_128_60.jld2", "a+") -bathymetry_file["bathymetry"] = bathymetry -close(bathymetry_file) - -println("Bathymetry: min= ", minimum(bathymetry)," , max= ", maximum(bathymetry)) -println("τˣ: min= ", minimum(τˣ)," , max= ", maximum(τˣ)) -println("τʸ: min= ", minimum(τʸ)," , max= ", maximum(τʸ)) -println("SST★: min= ", minimum(target_sea_surface_temperature)," , max= ",maximum(target_sea_surface_temperature)) - -fig = Figure(resolution = (1920, 1080)) -ax = fig[1, 1] = LScene(fig, title = "Bathymetry") -heatmap!(ax, bathymetry) - -ax = fig[2, 1] = LScene(fig, title = "Bathymetry") -heatmap!(ax, target_sea_surface_temperature) - -ax = fig[1, 2] = LScene(fig, title = "East-west wind stress") -heatmap!(ax, τˣ) - -ax = fig[2, 2] = LScene(fig, title = "North-south wind stress") -heatmap!(ax, τʸ) - -display(fig) - diff --git a/validation/near_global_lat_lon/annual_cycle_earth_bathymetry.jl b/validation/near_global_lat_lon/annual_cycle_earth_bathymetry.jl deleted file mode 100644 index d853402906..0000000000 --- a/validation/near_global_lat_lon/annual_cycle_earth_bathymetry.jl +++ /dev/null @@ -1,334 +0,0 @@ -using Statistics -using JLD2 -using Printf -using Oceananigans -using Oceananigans.Units -using Oceananigans.Utils - -using Oceananigans.Coriolis: HydrostaticSphericalCoriolis -using Oceananigans.Architectures: arch_array -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom -using Oceananigans.TurbulenceClosures -using Oceananigans.Advection: VelocityStencil -using CUDA: @allowscalar -using Oceananigans.Operators: Δzᵃᵃᶜ - -include("cyclic_interpolate_utils.jl") - -##### -##### Grid -##### - -latitude = (-84.375, 84.375) -Δφ = latitude[2] - latitude[1] - -# 2.8125 degree resolution -Nx = 128 -Ny = 60 -Nz = 18 - -output_prefix = "annual_cycle_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_temp" - -arch = CPU() -reference_density = 1035 - -##### -##### Load forcing files roughly from CORE2 paper -##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html) -##### - -using DataDeps - -path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/" - -dh = DataDep("near_global_lat_lon_3_degrees", - "Forcing data for global latitude longitude simulation", - [path * "bathymetry_lat_lon_128x60_FP32.bin", - path * "sea_surface_temperature_25_128x60x12.jld2", - path * "tau_x_128x60x12.jld2", - path * "tau_y_128x60x12.jld2"] -) - -DataDeps.register(dh) - -datadep"near_global_lat_lon_3_degrees" - -##### -##### Load forcing files roughly from CORE2 paper -##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html) -##### - -filename = [:sea_surface_temperature_25_128x60x12, :tau_x_128x60x12, :tau_y_128x60x12] - -for name in filename - datadep_path = @datadep_str "near_global_lat_lon_3_degrees/" * string(name) * ".jld2" - file = Symbol(:file_, name) - @eval $file = jldopen($datadep_path) -end - -bathymetry_data = Array{Float32}(undef, Nx*Ny) -bathymetry_path = @datadep_str "near_global_lat_lon_3_degrees/bathymetry_lat_lon_128x60_FP32.bin" -read!(bathymetry_path, bathymetry_data) - -bathymetry_data = bswap.(bathymetry_data) |> Array{Float64} -bathymetry_data = reshape(bathymetry_data, Nx, Ny) - -τˣ = zeros(Nx, Ny, Nmonths) -τʸ = zeros(Nx, Ny, Nmonths) -T★ = zeros(Nx, Ny, Nmonths) - -for month in 1:Nmonths - τˣ[:, :, month] = file_tau_x_128x60x12["tau_x/$month"] ./ reference_density - τʸ[:, :, month] = file_tau_y_128x60x12["tau_y/$month"] ./ reference_density - T★[:, :, month] = file_sea_surface_temperature_25_128x60x12["sst25/$month"] -end - -bathymetry = arch_array(arch, bathymetry_data) - -H = 3600.0 -# H = - minimum(bathymetry) - -# Uncomment for a flat bottom: -# bathymetry = - H .* (bathymetry .< -10) - -# A spherical domain -@show underlying_grid = LatitudeLongitudeGrid(arch, - size = (Nx, Ny, Nz), - longitude = (-180, 180), - latitude = latitude, - halo = (5, 5, 5), - z = (-H, 0), - precompute_metrics = true) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -τˣ = arch_array(arch, - τˣ) -τʸ = arch_array(arch, - τʸ) - -target_sea_surface_temperature = T★ = arch_array(arch, T★) - -##### -##### Physics and model setup -##### - -νh = 1e+5 -νz = 1e+1 -κh = 1e+3 -κz = 1e-4 - -vertical_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), - ν = νz, κ = κz) - -horizontal_closure = HorizontalScalarDiffusivity(ν = νh, κ = κh) - -convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0) - -gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = 1000.0, - κ_symmetric = 900.0, - slope_limiter = FluxTapering(1e-2)) - -##### -##### Boundary conditions / constant-in-time surface forcing -##### - -Δz_top = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.underlying_grid) -Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.underlying_grid) - -@inline function surface_temperature_relaxation(i, j, grid, clock, fields, p) - time = clock.time - - n₁ = current_time_index(time) - n₂ = next_time_index(time) - - @inbounds begin - T★₁ = p.T★[i, j, n₁] - T★₂ = p.T★[i, j, n₂] - T_surface = fields.T[i, j, grid.Nz] - end - - T★ = cyclic_interpolate(T★₁, T★₂, time) - - return p.λ * (T_surface - T★) -end - -T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation, - discrete_form = true, - parameters = (λ = Δz_top/3days, T★ = target_sea_surface_temperature)) - -@inline function wind_stress(i, j, grid, clock, fields, τ) - time = clock.time - n₁ = current_time_index(time) - n₂ = next_time_index(time) - - @inbounds begin - τ₁ = τ[i, j, n₁] - τ₂ = τ[i, j, n₂] - end - - return cyclic_interpolate(τ₁, τ₂, time) -end - -u_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τˣ) -v_wind_stress_bc = FluxBoundaryCondition(wind_stress, discrete_form = true, parameters = τʸ) - -@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1] -@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1] - -# Linear bottom drag: -μ = Δz_bottom / 10days - -u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ) -v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ) - -u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc) -v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc) -T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc) - -free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver) - -equation_of_state=LinearEquationOfState(thermal_expansion=2e-4) - -model = HydrostaticFreeSurfaceModel(grid = grid, - free_surface = free_surface, - momentum_advection = WENO(vector_invariant=VelocityStencil()), - tracer_advection = WENO(), - coriolis = HydrostaticSphericalCoriolis(), - boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs), - buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=30), - tracers = :T, - closure = (vertical_closure, convective_adjustment, gent_mcwilliams_diffusivity)) - -##### -##### Initial condition: -##### - -u, v, w = model.velocities -η = model.free_surface.η -T = model.tracers.T -T .= -1 - -##### -##### Simulation setup -##### - -Δt = 20minutes - -simulation = Simulation(model, Δt = Δt, stop_time = 5years) - -start_time = [time_ns()] - -function progress(sim) - wall_time = (time_ns() - start_time[1]) * 1e-9 - - η = model.free_surface.η - u = model.velocities.u - @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, max(|w|): %.2e ms⁻¹, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, - maximum(abs, u), maximum(abs, w), - prettytime(wall_time)) - - start_time[1] = time_ns() - - return nothing -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) - -u, v, w = model.velocities - -T = model.tracers.T -η = model.free_surface.η - -save_interval = 5days - -simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, T, η), - schedule = TimeInterval(save_interval), - filename = output_prefix * "_surface", - indices = (:, :, grid.Nz), - overwrite_existing = true) - -simulation.output_writers[:atlantic] = JLD2OutputWriter(model, (; u, v, T, η), - schedule = TimeInterval(save_interval), - filename = output_prefix * "_atlantic", - indices = (60, :, :), - overwrite_existing = true) - -simulation.output_writers[:checkpointer] = Checkpointer(model, - schedule = TimeInterval(1year), - prefix = output_prefix * "_checkpoint", - overwrite_existing = true) - -# Let's goo! -@info "Running with Δt = $(prettytime(simulation.Δt))" - -run!(simulation) - -@info """ - - Simulation took $(prettytime(simulation.run_wall_time)) - Free surface: $(typeof(model.free_surface).name.wrapper) - Time step: $(prettytime(Δt)) -""" - -#### -#### Visualize solution -#### - -using GLMakie, JLD2 - -output_prefix = "annual_cycle_global_lat_lon_128_60_18_temp" - -surface_file = jldopen(output_prefix * "_surface.jld2") - -iterations = parse.(Int, keys(surface_file["timeseries/t"])) - -iter = Observable(0) - -ηi(iter) = surface_file["timeseries/η/" * string(iter)][:, :, 1] -ui(iter) = surface_file["timeseries/u/" * string(iter)][:, :, 1] -vi(iter) = surface_file["timeseries/v/" * string(iter)][:, :, 1] -Ti(iter) = surface_file["timeseries/T/" * string(iter)][:, :, 1] -ti(iter) = string(surface_file["timeseries/t/" * string(iter)] / day) - -η = @lift ηi($iter) -u = @lift ui($iter) -v = @lift vi($iter) -T = @lift Ti($iter) - -max_η = 2 -min_η = - max_η -max_u = 0.2 -min_u = - max_u -max_T = 32 -min_T = 0 - -fig = Figure(resolution = (1200, 900)) - -ax = Axis(fig[1, 1], title="Free surface displacement (m)") -hm = GLMakie.heatmap!(ax, η, colorrange=(min_η, max_η), colormap=:balance) -cb = Colorbar(fig[1, 2], hm) - -ax = Axis(fig[2, 1], title="Sea surface temperature (ᵒC)") -hm = GLMakie.heatmap!(ax, T, colorrange=(min_T, max_T), colormap=:thermal) -cb = Colorbar(fig[2, 2], hm) - -ax = Axis(fig[1, 3], title="East-west surface velocity (m s⁻¹)") -hm = GLMakie.heatmap!(ax, u, colorrange=(min_u, max_u), colormap=:balance) -cb = Colorbar(fig[1, 4], hm) - -ax = Axis(fig[2, 3], title="North-south surface velocity (m s⁻¹)") -hm = GLMakie.heatmap!(ax, v, colorrange=(min_u, max_u), colormap=:balance) -cb = Colorbar(fig[2, 4], hm) - -title_str = @lift "Earth day = " * ti($iter) -ax_t = fig[0, :] = Label(fig, title_str) - -GLMakie.record(fig, output_prefix * ".mp4", iterations, framerate=8) do i - @info "Plotting iteration $i of $(iterations[end])..." - iter[] = i -end - -display(fig) - -close(surface_file) diff --git a/validation/near_global_lat_lon/annual_cycle_forcing.jl b/validation/near_global_lat_lon/annual_cycle_forcing.jl deleted file mode 100644 index c9f9ab375b..0000000000 --- a/validation/near_global_lat_lon/annual_cycle_forcing.jl +++ /dev/null @@ -1,125 +0,0 @@ -using Printf -using GLMakie -using Oceananigans -using Oceananigans.Units - -include("cyclic_interpolate_utils.jl") - -using DataDeps - -path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/" - -dh = DataDep("near_global_lat_lon", - "Forcing data for global latitude longitude simulation", - [path * "bathymetry_lat_lon_128x60_FP32.bin", - path * "sea_surface_temperature_25_128x60x12.jld2", - path * "tau_x_128x60x12.jld2", - path * "tau_y_128x60x12.jld2"] -) - -DataDeps.register(dh) - -datadep"near_global_lat_lon" - -# 2.8125 degree resolution -Nx = 128 -Ny = 60 -Nz = 18 -reference_density = 1035 - -##### -##### Load forcing files roughly from CORE2 paper -##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html) -##### - -filename = [:sea_surface_temperature_25_128x60x12, :tau_x_128x60x12, :tau_y_128x60x12] - -for name in filename - datadep_path = @datadep_str "near_global_lat_lon/" * string(name) * ".jld2" - file = Symbol(:file_, name) - @eval $file = jldopen($datadep_path) -end - -bathymetry_data = Array{Float32}(undef, Nx*Ny) -bathymetry_path = @datadep_str "near_global_lat_lon/bathymetry_lat_lon_128x60_FP32.bin" -read!(bathymetry_path, bathymetry_data) - -bathymetry_data = bswap.(bathymetry_data) |> Array{Float64} -bathymetry_data = reshape(bathymetry_data, Nx, Ny) - -τˣ = zeros(Nx, Ny, Nmonths) -τʸ = zeros(Nx, Ny, Nmonths) -T★ = zeros(Nx, Ny, Nmonths) - -for month in 1:Nmonths - τˣ[:, :, month] = file_tau_x_128x60x12["tau_x/$month"] ./ reference_density - τʸ[:, :, month] = file_tau_y_128x60x12["tau_y/$month"] ./ reference_density - T★[:, :, month] = file_sea_surface_temperature_25_128x60x12["sst25/$month"] -end - -times = 0:1days:24*30days -discrete_times = 0:30days:11*30days - -#= -i = 80 -j = 12 # Southern ocean? -Nt = length(times) -τ_ij = zeros(Nt) - -for (m, t) in enumerate(times) - n₁ = current_time_index(t) - n₂ = next_time_index(t) - τ_ij[m] = cyclic_interpolate(τˣ[i, j, n₁], τˣ[i, j, n₂], t) -end - -fig = Figure(resolution = (1200, 600)) - -ax = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Target SST at (i, j) = ($i, $j)") -scatter!(ax, discrete_times ./ day, τˣ[i, j, :], markersize=20, marker=:utriangle, color=:red) -lines!(ax, times ./ day, τ_ij) - -display(fig) -=# - -##### -##### Analyze forcing -##### - -fig = Figure(resolution = (1200, 800)) - -t = Node(0.0) - -T★_continuous(t) = cyclic_interpolate(T★, t) -τˣ_continuous(t) = cyclic_interpolate(τˣ, t) -τʸ_continuous(t) = cyclic_interpolate(τʸ, t) - -T★_t = @lift T★_continuous($t) -τˣ_t = @lift τˣ_continuous($t) -τʸ_t = @lift τʸ_continuous($t) - -max_τ = maximum(abs, τˣ) -max_T = maximum(T★) -min_T = minimum(T★) - -fig = Figure(resolution = (1200, 600)) - -ax_τˣ = Axis(fig[1, 1], title="East-west wind stress (m² s⁻²)") -hm_τˣ = heatmap!(ax_τˣ, τˣ_t, colorrange=(-max_τ, max_τ), colormap=:balance) -cb_τˣ = Colorbar(fig[1, 2], hm_τˣ) - -ax_τʸ = Axis(fig[2, 1], title="North-south wind stress (m² s⁻²)") -hm_τʸ = heatmap!(ax_τʸ, τʸ_t, colorrange=(-max_τ, max_τ), colormap=:balance) -cb_τʸ = Colorbar(fig[2, 2], hm_τʸ) - -ax_T = Axis(fig[3, 1], title="Sea surface temperature (ᵒC)") -hm_T = heatmap!(ax_T, T★_t, colorrange=(min_T, max_T), colormap=:thermal) -cb_T = Colorbar(fig[3, 2], hm_T) - -title_str = @lift "Earth day = " * prettytime($t) -ax_t = fig[0, :] = Label(fig, title_str) - -display(fig) - -record(fig, "annual_cycle_forcing.mp4", times, framerate=8) do tn - t[] = tn -end diff --git a/validation/near_global_lat_lon/constant_forcing_earth_bathymetry.jl b/validation/near_global_lat_lon/constant_forcing_earth_bathymetry.jl deleted file mode 100644 index 093f50f571..0000000000 --- a/validation/near_global_lat_lon/constant_forcing_earth_bathymetry.jl +++ /dev/null @@ -1,287 +0,0 @@ -using Statistics -using JLD2 -using Printf -using GLMakie -using CUDA -using Oceananigans -using Oceananigans.Coriolis: HydrostaticSphericalCoriolis -using Oceananigans.Architectures: arch_array -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom -using Oceananigans.TurbulenceClosures: HorizontallyCurvilinearAnisotropicDiffusivity, - Explicit, - VerticallyImplicit -using Oceananigans.Units -using Oceananigans.Operators: Δzᵃᵃᶜ - -##### -##### Grid -##### - -latitude = (-84.375, 84.375) -Δφ = latitude[2] - latitude[1] - -# 2.8125 degree resolution -Nx = 128 -Ny = 60 -Nz = 18 - -arch = CPU() -reference_density = 1035 - -##### -##### Load forcing files roughly from CORE2 paper -##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html) -##### - -bathymetry_path = "bathy_128x60var4.bin" -east_west_stress_path = "off_TAUXvar1.bin" -north_south_stress_path = "off_TAUY.bin" -sea_surface_temperature_path="sst25_128x60x12.bin" - -bytes = sizeof(Float32) * Nx * Ny -bathymetry = reshape(bswap.(reinterpret(Float32, read(bathymetry_path, bytes))), (Nx, Ny)) -τˣ = - reshape(bswap.(reinterpret(Float32, read(east_west_stress_path, bytes))), (Nx, Ny)) ./ reference_density -τʸ = - reshape(bswap.(reinterpret(Float32, read(north_south_stress_path, bytes))), (Nx, Ny)) ./ reference_density -target_sea_surface_temperature = reshape(bswap.(reinterpret(Float32, read(sea_surface_temperature_path, bytes))), (Nx, Ny)) - -bathymetry = arch_array(arch, bathymetry) -τˣ = arch_array(arch, τˣ) -τʸ = arch_array(arch, τʸ) -target_sea_surface_temperature = T★ = arch_array(arch, target_sea_surface_temperature) - -H = 3600.0 -# bathymetry = - H .* (bathymetry .< -10) -# H = - minimum(bathymetry) - -# A spherical domain -@show underlying_grid = LatitudeLongitudeGrid(arch, - size = (Nx, Ny, Nz), - longitude = (-180, 180), - latitude = latitude, - halo = (3, 3, 3), - z = (-H, 0)) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -##### -##### Physics and model setup -##### - -νh = 1e+5 -νz = 1e-2 -κh = 1e+3 -κz = 1e-4 - -# Fully explicit -background_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz, time_discretization=Explicit()) -convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0, time_discretization=Explicit()) - -# Vertically-implicit, horizontally-explicit -#background_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz, time_discretization=VerticallyImplicit()) -#convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0, time_discretization=VerticallyImplicit()) - -##### -##### Boundary conditions / constant-in-time surface forcing -##### - -Δz_top = CUDA.@allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.grid) -Δz_bottom = CUDA.@allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.grid) - -@inline surface_temperature_relaxation(i, j, grid, clock, fields, p) = @inbounds p.λ * (fields.T[i, j, grid.Nz] - p.T★[i, j]) - -T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation, - discrete_form = true, - parameters = (λ = Δz_top/3days, T★ = target_sea_surface_temperature)) - -u_wind_stress_bc = FluxBoundaryCondition(τˣ) -v_wind_stress_bc = FluxBoundaryCondition(τʸ) - -@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1] -@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1] - -# Linear bottom drag: -μ = 1 / 10days * Δz_bottom - -u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ) -v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ) - -u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc) -v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc) -T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc) - -equation_of_state = LinearEquationOfState(thermal_expansion=2e-4) - -model = HydrostaticFreeSurfaceModel(grid = grid, - momentum_advection = VectorInvariant(), - tracer_advection = WENO(), - coriolis = HydrostaticSphericalCoriolis(), - boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs), - buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=true), - tracers = (:T, :S), - closure = (background_diffusivity, convective_adjustment)) - -#= -##### -##### Visualize the setup -##### - -setup_fig = Figure(resolution = (1200, 900)) - -free_surface = model.free_surface -∫Ax = free_surface.implicit_step_solver.vertically_integrated_lateral_areas.xᶠᶜᶜ -∫Ay = free_surface.implicit_step_solver.vertically_integrated_lateral_areas.yᶜᶠᶜ - -h_cpu = Array(bathymetry) -τˣ_cpu = Array(τˣ) -τʸ_cpu = Array(τʸ) -T★_cpu = Array(T★) -∫Ax_cpu = Array(parent(∫Ax))[:, :, 1] -∫Ay_cpu = Array(parent(∫Ay))[:, :, 1] - -ax_b = Axis(setup_fig[1, 1], title="Bathymetry (m)") -hm_b = heatmap!(ax_b, h_cpu) -cb_b = Colorbar(setup_fig[1, 2], hm_b) - -ax_T★ = Axis(setup_fig[1, 3], title="Target sea surface temperature (ᵒC)") -hm_T★ = heatmap!(ax_T★, T★_cpu) -cb_T★ = Colorbar(setup_fig[1, 4], hm_T★) - -ax_x = Axis(setup_fig[2, 1], title="Vertically-integrated x-face area (m³)") -hm_x = heatmap!(ax_x, ∫Ax_cpu) -cb_x = Colorbar(setup_fig[2, 2], hm_x) - -ax_y = Axis(setup_fig[2, 3], title="Vertically-integrated y-face area (m³)") -hm_y = heatmap!(ax_y, ∫Ay_cpu) -cb_y = Colorbar(setup_fig[2, 4], hm_y) - -ax_τˣ = Axis(setup_fig[3, 1], title="East-west wind stress (m² s⁻²)") -hm_τˣ = heatmap!(ax_τˣ, τˣ_cpu) -cb_τˣ = Colorbar(setup_fig[3, 2], hm_τˣ) - -ax_τʸ = Axis(setup_fig[3, 3], title="North-south wind stress (m² s⁻²)") -hm_τʸ = heatmap!(ax_τʸ, τʸ_cpu) -cb_τʸ = Colorbar(setup_fig[3, 4], hm_τʸ) - -display(setup_fig) -=# - -##### -##### Initial condition: -##### - -u, v, w = model.velocities -η = model.free_surface.η - -T = model.tracers.T -T .= 5 - -S = model.tracers.S -S .= 30 - -##### -##### Simulation setup -##### - -# Time-scale for gravity wave propagation across the smallest grid cell -g = model.free_surface.gravitational_acceleration -gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed - -#minimum_Δx = abs(grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ[1:grid.Ny])) * deg2rad(grid.Δλ)) -#minimum_Δy = abs(grid.radius * deg2rad(grid.Δφ)) -#wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed - -if model.free_surface isa ExplicitFreeSurface - Δt = 60seconds #0.2 * minimum_Δx / gravity_wave_speed -else - Δt = 20minutes -end - -simulation = Simulation(model, Δt = Δt, stop_time = 1day) - -start_time = [time_ns()] - -function progress(sim) - wall_time = (time_ns() - start_time[1]) * 1e-9 - - η = model.free_surface.η - - if model.free_surface isa ExplicitFreeSurface - @info @sprintf("Time: % 12s, iteration: %d, max(|η|): %.2e m, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, - maximum(abs, η), - prettytime(wall_time)) - else - @info @sprintf("Time: % 12s, iteration: %d, free surface iterations: %d, max(|η|): %.2e m, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, - sim.model.free_surface.implicit_step_solver.preconditioned_conjugate_gradient_solver.iteration, - maximum(abs, η), - prettytime(wall_time)) - end - - start_time[1] = time_ns() - - return nothing -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) - -output_fields = merge(model.velocities, model.tracers, (; η=model.free_surface.η)) -output_prefix = "global_lat_lon_$(grid.Nx)_$(grid.Ny)_$(grid.Nz)" - -simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields, - schedule = TimeInterval(10days), - filename = output_prefix, - overwrite_existing = true) - -# Let's goo! -@info "Running with Δt = $(prettytime(simulation.Δt))" - -run!(simulation) - -@info """ - - Background diffusivity: $background_diffusivity - Free surface: $(typeof(model.free_surface).name.wrapper) - Time step: $(prettytime(Δt)) - -""" -# Simulation took $(prettytime(simulation.run_wall_time)) -# Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale)) - -##### -##### Visualize solution -##### - -η_cpu = Array(interior(model.free_surface.η))[:, :, 1] -T_cpu = Array(interior(model.tracers.T))[:, :, end] -u_cpu = Array(interior(model.velocities.u))[:, :, end] -v_cpu = Array(interior(model.velocities.v))[:, :, end] - -max_η = maximum(abs, η_cpu) -max_u = maximum(abs, u_cpu) -max_v = maximum(abs, v_cpu) - -max_T = maximum(T_cpu) -min_T = minimum(T_cpu) - -solution_fig = Figure(resolution = (1200, 600)) - -ax_η = Axis(solution_fig[1, 1], title="Free surface displacement (m)") -hm_η = heatmap!(ax_η, η_cpu, colorrange=(-max_η, max_η), colormap=:balance) -cb_η = Colorbar(solution_fig[1, 2], hm_η) - -ax_T = Axis(solution_fig[2, 1], title="Sea surface temperature (ᵒC)") -hm_T = heatmap!(ax_T, T_cpu, colorrange=(min_T, max_T), colormap=:thermal) -cb_T = Colorbar(solution_fig[2, 2], hm_T) - -ax_u = Axis(solution_fig[1, 3], title="East-West velocity (m s⁻¹)") -hm_u = heatmap!(ax_u, u_cpu, colorrange=(-max_u, max_u), colormap=:balance) -cb_u = Colorbar(solution_fig[1, 4], hm_u) - -ax_v = Axis(solution_fig[2, 3], title="North-South velocity (m s⁻¹)") -hm_v = heatmap!(ax_v, v_cpu, colorrange=(-max_v, max_v), colormap=:balance) -cb_v = Colorbar(solution_fig[2, 4], hm_v) - -display(solution_fig) diff --git a/validation/near_global_lat_lon/cyclic_interpolate_utils.jl b/validation/near_global_lat_lon/cyclic_interpolate_utils.jl deleted file mode 100644 index 781ea0f978..0000000000 --- a/validation/near_global_lat_lon/cyclic_interpolate_utils.jl +++ /dev/null @@ -1,30 +0,0 @@ -const thirty_days = 30days -const Nmonths = 12 - -#@inline current_time_index(time, interval=2592000, length=12) = mod(unsafe_trunc(Int32, time / interval), length) + 1 -#@inline next_time_index(time, interval=2592000, length=12) = mod(unsafe_trunc(Int32, time / interval) + 1, length) + 1 -#@inline cyclic_interpolate(u₁::Number, u₂, time, interval=2592000) = u₁ + mod(time / interval, 1) * (u₂ - u₁) - -@inline current_time_index(time) = mod(unsafe_trunc(Int32, time / 2592000), 12) + 1 -@inline next_time_index(time) = mod(unsafe_trunc(Int32, time / 2592000) + 1, 12) + 1 -@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / 2592000, 1) * (u₂ - u₁) - -@inline function cyclic_interpolate(τ::AbstractArray, time, interval=2592000, length=12) - n₁ = current_time_index(time, interval, length) - n₂ = next_time_index(time, interval, length) - return cyclic_interpolate.(view(τ, :, :, n₁), view(τ, :, :, n₂), time, interval) -end - -@inline function interpolate_fluxes(old_array, Nx_old, Ny_old, Nx_new, Ny_new) - old_grid = LatitudeLongitudeGrid(size = (Nx_old, Ny_old, 1), latitude = (-80, 80), longitude = (-180, 180), z = (0, 1)) - new_grid = LatitudeLongitudeGrid(size = (Nx_new, Ny_new, 1), latitude = (-80, 80), longitude = (-180, 180), z = (0, 1)) - - old_field = Field{Center, Center, Center}(old_grid) - set!(old_field, old_array) - new_array = zeros(Nx_new, Ny_new) - - for i in 1:Nx_new, j in 1:Ny_new - new_array[i, j] = interpolate(old_field, new_grid.λᶜᵃᵃ[i], new_grid.φᵃᶜᵃ[j], old_grid.zᵃᵃᶜ[1]) - end - return new_array -end \ No newline at end of file diff --git a/validation/near_global_lat_lon/freely_decaying_barotropic_turbulence.jl b/validation/near_global_lat_lon/freely_decaying_barotropic_turbulence.jl deleted file mode 100644 index 75ca9c7c04..0000000000 --- a/validation/near_global_lat_lon/freely_decaying_barotropic_turbulence.jl +++ /dev/null @@ -1,225 +0,0 @@ -# # Freely decaying barotropic turbulence on a latitude-longitude strip - -using Oceananigans -using Oceananigans.Grids - -using Oceananigans.BoundaryConditions: fill_halo_regions! - -using Oceananigans.Coriolis: HydrostaticSphericalCoriolis - -using Oceananigans.Models.HydrostaticFreeSurfaceModels: - HydrostaticFreeSurfaceModel, - VerticalVorticityField, - VectorInvariant, - ExplicitFreeSurface, - ImplicitFreeSurface - -using Oceananigans.TurbulenceClosures - -using Oceananigans.Utils: prettytime, hours, day, days, years, year -using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval, IterationInterval -using Oceananigans.AbstractOperations: KernelFunctionOperation - -using Statistics -using JLD2 -using Printf -using CUDA - -##### -##### Grid -##### - -precompute = true - -latitude = (-80, 80) -Δφ = latitude[2] - latitude[1] - -resolution = 1/3 # degree -Nx = round(Int, 360 / resolution) -Ny = round(Int, Δφ / resolution) - -# A spherical domain -@show grid = LatitudeLongitudeGrid(GPU(), size = (Nx, Ny, 1), - longitude = (-180, 180), - latitude = latitude, - halo = (2, 2, 2), - z = (-100, 0), - precompute_metrics = precompute) - -##### -##### Physics and model setup -##### - -free_surface = ExplicitFreeSurface(gravitational_acceleration=1.0) - -CUDA.allowscalar(true) - -equatorial_Δx = grid.radius * deg2rad(mean(grid.Δλᶜᵃᵃ)) -diffusive_time_scale = 120days - -@show const νh₂ = equatorial_Δx^2 / diffusive_time_scale -@show const νh₄ = 1e-5 * equatorial_Δx^4 / diffusive_time_scale - -closure = ScalarBiharmonicDiffusivity(ν=νh₄) - -coriolis = HydrostaticSphericalCoriolis() -Ω = coriolis.rotation_rate / 20 -coriolis = HydrostaticSphericalCoriolis(rotation_rate=Ω) - -model = HydrostaticFreeSurfaceModel(grid = grid, - momentum_advection = VectorInvariant(), - free_surface = free_surface, - coriolis = coriolis, - tracers = nothing, - buoyancy = nothing, - closure = closure) - -##### -##### Initial condition: two streamfunction -##### - -# Random noise -ψ★ = Field(Face, Face, Center, model.architecture, model.grid) -set!(ψ★, (x, y, z) -> rand()) -fill_halo_regions!(ψ★) - -# Zonal wind -step(x, d, c) = 1/2 * (1 + tanh((x - c) / d)) -polar_mask(y) = step(y, -5, 60) * step(y, 5, -60) -zonal_ψ(y) = (cosd(4y)^3 + 0.5 * exp(-y^2 / 200)) * polar_mask(y) - -ψ̄ = Field((Face, Face, Center), model.grid) -set!(ψ̄, (x, y, z) -> zonal_ψ(y)) -fill_halo_regions!(ψ̄) - -ψ_total = 40 * ψ★ + Ny * ψ̄ - -u, v, w = model.velocities -η = model.free_surface.η - -if !isnothing(model.coriolis) - using Oceananigans.Coriolis: fᶠᶠᵃ - f = KernelFunctionOperation{Face, Face, Center}(fᶠᶠᵃ, model.grid, parameters=model.coriolis) - g = model.free_surface.gravitational_acceleration - η .= f * ψ_total / g -end - -u .= - ∂y(ψ_total) -v .= + ∂x(ψ_total) - -##### -##### Shenanigans for rescaling the velocity field to -##### 1. Have a magnitude (ish) that's a fixed fraction of -##### the surface gravity wave speed; -##### 2. Zero volume mean on the curvilinear LatitudeLongitudeGrid. -##### - -# Time-scale for gravity wave propagation across the smallest grid cell -g = model.free_surface.gravitational_acceleration -gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed - -minimum_Δx = grid.radius * cosd(maximum(abs, grid.φᵃᶜᵃ)) * deg2rad(maximum(abs, grid.Δλᶜᵃᵃ)) -minimum_Δy = grid.radius * deg2rad(minimum(abs, grid.Δφᵃᶜᵃ)) -wave_propagation_time_scale = min(minimum_Δx, minimum_Δy) / gravity_wave_speed - -@info "Max speeds prior to rescaling:" -@show max_u = maximum(u) -@show max_v = maximum(v) -max_speed_ish = sqrt(max_u^2 + max_v^2) - -target_speed = 0.5 * gravity_wave_speed -u .*= target_speed / max_speed_ish -v .*= target_speed / max_speed_ish - -# Zero out mean motion -using Oceananigans.AbstractOperations: volume - -u_cpu = XFaceField(CPU(), grid) -v_cpu = YFaceField(CPU(), grid) -set!(u_cpu, u) -set!(v_cpu, v) - -@show max_u = maximum(u) -@show max_v = maximum(v) - -u_dV = u_cpu * volume -u_reduced = Field(Average(u_dV, dims=(1, 2, 3))) -compute!(u_reduced) -mean!(u_reduced, u_dV) -integrated_u = u_reduced[1, 1, 1] - -v_dV = v_cpu * volume -v_reduced = Field(Average(v_dV, dims=(1, 2, 3))) -mean!(v_reduced, v_dV) -integrated_v = v_reduced[1, 1, 1] - -# Calculate total volume -u_cpu .= 1 -v_cpu .= 1 -compute!(u_reduced) -compute!(v_reduced) - -u_volume = u_reduced[1, 1, 1] -v_volume = v_reduced[1, 1, 1] - -@info "Max speeds prior zeroing out volume mean:" -@show maximum(u) -@show maximum(v) - -u .-= integrated_u / u_volume -v .-= integrated_v / v_volume - -@info "Initial max speeds:" -@show maximum(u) -@show maximum(v) - -##### -##### Simulation setup -##### - -ζ = VerticalVorticityField(model) -compute!(ζ) -Δt = 0.2 * minimum_Δx / gravity_wave_speed - -@info """ - Maximum vertical vorticity: $(maximum(ζ)) - Inverse maximum vertical vorticity: $(prettytime(1/maximum(ζ))) - Minimum wave propagation time scale: $(prettytime(wave_propagation_time_scale)) - Time step: $(prettytime(Δt)) -""" - -mutable struct Progress; interval_start_time::Float64; end - -function (p::Progress)(sim) - wall_time = (time_ns() - p.interval_start_time) * 1e-9 - - compute!(ζ) - - @info @sprintf("Time: %s, iteration: %d, max(|ζ|): %.2e s⁻¹, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, - maximum(abs, ζ), - prettytime(wall_time)) - - p.interval_start_time = time_ns() - - return nothing -end - -simulation = Simulation(model, - Δt = Δt, - stop_time = 100days, - iteration_interval = 1000, - progress = Progress(time_ns())) - -output_fields = merge(model.velocities, (η=model.free_surface.η, ζ=ζ)) - -output_prefix = "rotating_freely_decaying_barotropic_turbulence_Nx$(grid.Nx)_Ny$(grid.Ny)" - -simulation.output_writers[:fields] = JLD2OutputWriter(model, (ζ = ζ,), - schedule = TimeInterval(10day), - filename = output_prefix, - overwrite_existing = true) - -# Let's goo! -run!(simulation) \ No newline at end of file diff --git a/validation/near_global_lat_lon/near_global_quarter_degree.jl b/validation/near_global_lat_lon/near_global_quarter_degree.jl deleted file mode 100644 index b0a5107f85..0000000000 --- a/validation/near_global_lat_lon/near_global_quarter_degree.jl +++ /dev/null @@ -1,311 +0,0 @@ -using Statistics -using JLD2 -using Printf -using Oceananigans -using Oceananigans.Units - -using Oceananigans.Fields: interpolate, Field -using Oceananigans.Advection: VelocityStencil -using Oceananigans.Architectures: arch_array -using Oceananigans.Coriolis: HydrostaticSphericalCoriolis -using Oceananigans.BoundaryConditions -using Oceananigans.ImmersedBoundaries: inactive_node, peripheral_node -using CUDA: @allowscalar -using Oceananigans.Operators -using Oceananigans.Operators: Δzᵃᵃᶜ -using Oceananigans: prognostic_fields - -@inline function visualize(field, lev, dims) - (dims == 1) && (idx = (lev, :, :)) - (dims == 2) && (idx = (:, lev, :)) - (dims == 3) && (idx = (:, :, lev)) - - r = deepcopy(Array(interior(field)))[idx...] - r[r.==0] .= NaN - return r -end - -##### -##### Grid -##### - -arch = GPU() -reference_density = 1029 - -latitude = (-75, 75) - -# 0.25 degree resolution -Nx = 1440 -Ny = 600 -Nz = 48 - -const Nyears = 1 -const Nmonths = 12 -const thirty_days = 30days - -output_prefix = "near_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_fine" -pickup_file = false - -##### -##### Load forcing files and inital conditions from ECCO version 4 -##### https://ecco.jpl.nasa.gov/drive/files -##### Bathymetry is interpolated from ETOPO1 https://www.ngdc.noaa.gov/mgg/global/ -##### - -using DataDeps - -path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/ss/new_hydrostatic_data_after_cleared_bugs/quarter_degree_near_global_input_data/" - -datanames = ["z_faces-50-levels", - "bathymetry-1440x600", - "temp-1440x600-latitude-75", - "salt-1440x600-latitude-75", - "tau_x-1440x600-latitude-75", - "tau_y-1440x600-latitude-75", - "initial_conditions"] - -dh = DataDep("quarter_degree_near_global_lat_lon", - "Forcing data for global latitude longitude simulation", - [path * data * ".jld2" for data in datanames] -) - -DataDeps.register(dh) - -datadep"quarter_degree_near_global_lat_lon" - -files = [:file_z_faces, :file_bathymetry, :file_temp, :file_salt, :file_tau_x, :file_tau_y, :file_init] -for (data, file) in zip(datanames, files) - datadep_path = @datadep_str "quarter_degree_near_global_lat_lon/" * data * ".jld2" - @eval $file = jldopen($datadep_path) -end - -bathymetry = file_bathymetry["bathymetry"] - -τˣ = zeros(Nx, Ny, Nmonths) -τʸ = zeros(Nx, Ny, Nmonths) -T★ = zeros(Nx, Ny, Nmonths) -S★ = zeros(Nx, Ny, Nmonths) - -# Files contain 1 year (1992) of 12 monthly averages -τˣ = file_tau_x["field"] ./ reference_density -τʸ = file_tau_y["field"] ./ reference_density -T★ = file_temp["field"] -S★ = file_salt["field"] - -# Remember the convention!! On the surface a negative flux increases a positive decreases -bathymetry = arch_array(arch, bathymetry) - -# Stretched faces taken from ECCO Version 4 (50 levels in the vertical) -z_faces = file_z_faces["z_faces"][3:end] - -# A spherical domain -@show underlying_grid = LatitudeLongitudeGrid(arch, - size = (Nx, Ny, Nz), - longitude = (-180, 180), - latitude = latitude, - halo = (5, 5, 5), - z = z_faces, - precompute_metrics = true) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -τˣ = arch_array(arch, - τˣ) -τʸ = arch_array(arch, - τʸ) - -target_sea_surface_temperature = T★ = arch_array(arch, T★) -target_sea_surface_salinity = S★ = arch_array(arch, S★) - -##### -##### Physics and model setup -##### - -νz = 5e-3 -κz = 1e-4 - -convective_adjustment = RiBasedDiffusivity() -vertical_diffusivity = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), ν=νz, κ=κz) - -tracer_advection = WENO(underlying_grid) -momentum_advection = VectorInvariant(vorticity_scheme = WENO(), - divergence_scheme = WENO(), - vertical_scheme = WENO(underlying_grid)) - -##### -##### Boundary conditions / time-dependent fluxes -##### - -@inline current_time_index(time, tot_months) = mod(unsafe_trunc(Int32, time / thirty_days), tot_months) + 1 -@inline next_time_index(time, tot_months) = mod(unsafe_trunc(Int32, time / thirty_days) + 1, tot_months) + 1 -@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / thirty_days, 1) * (u₂ - u₁) - -Δz_top = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.underlying_grid) -Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.underlying_grid) - -@inline function surface_wind_stress(i, j, grid, clock, fields, τ) - time = clock.time - n₁ = current_time_index(time, Nmonths) - n₂ = next_time_index(time, Nmonths) - - @inbounds begin - τ₁ = τ[i, j, n₁] - τ₂ = τ[i, j, n₂] - end - - return cyclic_interpolate(τ₁, τ₂, time) -end - -u_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τˣ) -v_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τʸ) - -# Linear bottom drag: -μ = 0.001 # m s⁻¹ - -@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1] -@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1] -@inline u_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, k] -@inline v_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, k] - -u_immersed_bc = ImmersedBoundaryCondition(bottom = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form = true, parameters = μ)) -v_immersed_bc = ImmersedBoundaryCondition(bottom = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form = true, parameters = μ)) - -u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ) -v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ) - -@inline function surface_temperature_relaxation(i, j, grid, clock, fields, p) - time = clock.time - - n₁ = current_time_index(time, Nmonths) - n₂ = next_time_index(time, Nmonths) - - @inbounds begin - T★₁ = p.T★[i, j, n₁] - T★₂ = p.T★[i, j, n₂] - T_surface = fields.T[i, j, grid.Nz] - end - - T★ = cyclic_interpolate(T★₁, T★₂, time) - - return p.λ * (T_surface - T★) -end - -@inline function surface_salinity_relaxation(i, j, grid, clock, fields, p) - time = clock.time - - n₁ = current_time_index(time, Nmonths) - n₂ = next_time_index(time, Nmonths) - - @inbounds begin - S★₁ = p.S★[i, j, n₁] - S★₂ = p.S★[i, j, n₂] - S_surface = fields.S[i, j, grid.Nz] - end - - S★ = cyclic_interpolate(S★₁, S★₂, time) - - return p.λ * (S_surface - S★) -end - -T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation, - discrete_form = true, - parameters = (λ = Δz_top / 7days, T★ = target_sea_surface_temperature)) - -S_surface_relaxation_bc = FluxBoundaryCondition(surface_salinity_relaxation, - discrete_form = true, - parameters = (λ = Δz_top / 7days, S★ = target_sea_surface_salinity)) - -u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc, immersed = u_immersed_bc) -v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc, immersed = v_immersed_bc) -T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc) -S_bcs = FieldBoundaryConditions(top = S_surface_relaxation_bc) - -free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver) - -buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState()) - -model = HydrostaticFreeSurfaceModel(; grid, - free_surface, - momentum_advection, tracer_advection, - coriolis = HydrostaticSphericalCoriolis(), - buoyancy, - tracers = (:T, :S), - closure = (vertical_diffusivity, convective_adjustment), - boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs)) - -##### -##### Initial condition: -##### - -u, v, w = model.velocities -η = model.free_surface.η -T = model.tracers.T -S = model.tracers.S - -@info "Reading initial conditions" -T_init = file_init["T"] -S_init = file_init["S"] - -set!(model, T=T_init, S=S_init) -fill_halo_regions!(T) -fill_halo_regions!(S) - -@info "model initialized" - -##### -##### Simulation setup -##### - -Δt = 6minutes # probably we can go to 10min or 15min? - -simulation = Simulation(model, Δt = Δt, stop_time = Nyears*years) - -start_time = [time_ns()] - -using Oceananigans.Utils - -function progress(sim) - wall_time = (time_ns() - start_time[1]) * 1e-9 - - u = sim.model.velocities.u - η = sim.model.free_surface.η - - @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, maximum(abs, u), - prettytime(wall_time)) - - start_time[1] = time_ns() - - return nothing -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) - -u, v, w = model.velocities -T = model.tracers.T -S = model.tracers.S -η = model.free_surface.η - -output_fields = (; u, v, T, S, η) -save_interval = 5days - -# simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, T, S, η), -# schedule = TimeInterval(save_interval), -# filename = output_prefix * "_surface", -# indices = (:, :, grid.Nz), -# overwrite_existing = true) - -# simulation.output_writers[:checkpointer] = Checkpointer(model, -# schedule = TimeInterval(1year), -# prefix = output_prefix * "_checkpoint", -# overwrite_existing = true) - -# Let's goo! -@info "Running with Δt = $(prettytime(simulation.Δt))" - -run!(simulation, pickup = pickup_file) - -@info """ - Simulation took $(prettytime(simulation.run_wall_time)) - Free surface: $(typeof(model.free_surface).name.wrapper) - Time step: $(prettytime(Δt)) -""" diff --git a/validation/near_global_lat_lon/near_global_quarter_degree_one_level.jl b/validation/near_global_lat_lon/near_global_quarter_degree_one_level.jl deleted file mode 100644 index 545ee35214..0000000000 --- a/validation/near_global_lat_lon/near_global_quarter_degree_one_level.jl +++ /dev/null @@ -1,228 +0,0 @@ -using Statistics -using JLD2 -using Printf -using Plots -using Oceananigans -using Oceananigans.Units - -using Oceananigans.Fields: interpolate, Field -using Oceananigans.Architectures: arch_array -using Oceananigans.Coriolis: HydrostaticSphericalCoriolis -using Oceananigans.BoundaryConditions -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, inactive_node, peripheral_node -using CUDA: @allowscalar, device! -using Oceananigans.Operators -using Oceananigans.Operators: Δzᵃᵃᶜ -using Oceananigans: prognostic_fields - -using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField - -@inline function visualize(field, lev, dims) - (dims == 1) && (idx = (lev, :, :)) - (dims == 2) && (idx = (:, lev, :)) - (dims == 3) && (idx = (:, :, lev)) - - r = deepcopy(Array(interior(field)))[idx...] - r[ r.==0 ] .= NaN - return r -end - -##### -##### Grid -##### - -arch = GPU() -reference_density = 1029 - -latitude = (-75, 75) - -# 0.25 degree resolution -Nx = 1440 -Ny = 600 -Nz = 1 - -const Nyears = 1 -const Nmonths = 12 -const thirty_days = 30days - -output_prefix = "near_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_fine" -pickup_file = false - -##### -##### Load forcing files and inital conditions from ECCO version 4 -##### https://ecco.jpl.nasa.gov/drive/files -##### Bathymetry is interpolated from ETOPO1 https://www.ngdc.noaa.gov/mgg/global/ -##### - -using DataDeps - -path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/ss/new_hydrostatic_data_after_cleared_bugs/quarter_degree_near_global_input_data/" - -datanames = ["z_faces-50-levels", - "bathymetry-1440x600", - "temp-1440x600-latitude-75", - "salt-1440x600-latitude-75", - "tau_x-1440x600-latitude-75", - "tau_y-1440x600-latitude-75", - "initial_conditions"] - -dh = DataDep("quarter_degree_near_global_lat_lon", - "Forcing data for global latitude longitude simulation", - [path * data * ".jld2" for data in datanames] -) - -DataDeps.register(dh) - -datadep"quarter_degree_near_global_lat_lon" - -files = [:file_z_faces, :file_bathymetry, :file_temp, :file_salt, :file_tau_x, :file_tau_y, :file_init] -for (data, file) in zip(datanames, files) - datadep_path = @datadep_str "quarter_degree_near_global_lat_lon/" * data * ".jld2" - @eval $file = jldopen($datadep_path) -end - -bathymetry = file_bathymetry["bathymetry"] -bathymetry[bathymetry .< 0] .= -10e3 - -τˣ = zeros(Nx, Ny, Nmonths) -τʸ = zeros(Nx, Ny, Nmonths) - -# Files contain 1 year (1992) of 12 monthly averages -τˣ = file_tau_x["field"] ./ reference_density -τʸ = file_tau_y["field"] ./ reference_density - -# Remember the convention!! On the surface a negative flux increases a positive decreases -bathymetry = arch_array(arch, bathymetry) - -# A spherical domain -@show underlying_grid = LatitudeLongitudeGrid(arch, - size = (Nx, Ny, Nz), - longitude = (-180, 180), - latitude = latitude, - halo = (3, 3, 3), - z = (-10e3, 0), - #z = z_faces, - precompute_metrics = true) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -τˣ = arch_array(arch, - τˣ) -τʸ = arch_array(arch, - τʸ) - -##### -##### Boundary conditions / time-dependent fluxes -##### - -@inline current_time_index(time, tot_months) = mod(unsafe_trunc(Int32, time / thirty_days), tot_months) + 1 -@inline next_time_index(time, tot_months) = mod(unsafe_trunc(Int32, time / thirty_days) + 1, tot_months) + 1 -@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / thirty_days, 1) * (u₂ - u₁) - -Δz_top = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.underlying_grid) -Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.underlying_grid) - -@inline function surface_wind_stress(i, j, grid, clock, fields, τ) - time = clock.time - n₁ = current_time_index(time, Nmonths) - n₂ = next_time_index(time, Nmonths) - - @inbounds begin - τ₁ = τ[i, j, n₁] - τ₂ = τ[i, j, n₂] - end - - return cyclic_interpolate(τ₁, τ₂, time) -end - -u_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τˣ) -v_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τʸ) - -@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1] -@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1] - -# Linear bottom drag: -μ = 0.001 # ms⁻¹ - -u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ) -v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ) - -u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc) -v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc) - -free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver) - -buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState()) - -model = HydrostaticFreeSurfaceModel(grid = grid, - free_surface = free_surface, - momentum_advection = VectorInvariant(), - coriolis = HydrostaticSphericalCoriolis(), - buoyancy = nothing, - boundary_conditions = (u=u_bcs, v=v_bcs)) - -##### -##### Initial condition: -##### - -u, v, w = model.velocities -η = model.free_surface.η - -@info "model initialized" - -##### -##### Simulation setup -##### - -ζ = VerticalVorticityField(model) -compute!(ζ) - -Δt = 6minutes # for initialization, then we can go up to 6 minutes? -simulation = Simulation(model, Δt = Δt, stop_time = Nyears*years) -start_time = [time_ns()] - -using Oceananigans.Utils - -function progress(sim) - wall_time = (time_ns() - start_time[1]) * 1e-9 - - u = sim.model.velocities.u - η = sim.model.free_surface.η - - @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, wall time: %s", - prettytime(sim.model.clock.time), - sim.model.clock.iteration, maximum(abs, u), - prettytime(wall_time)) - - start_time[1] = time_ns() - - return nothing -end - -simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) - -u, v, w = model.velocities -η = model.free_surface.η - -output_fields = merge(model.velocities, (η=model.free_surface.η, ζ=ζ)) -save_interval = 1days - -simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, output_fields, #(; u, v, T, S, η), - schedule = TimeInterval(save_interval), - filename = output_prefix * "_surface", - indices = (:, :, grid.Nz), - overwrite_existing = true) - -simulation.output_writers[:checkpointer] = Checkpointer(model, - schedule = TimeInterval(1day), - prefix = output_prefix * "_checkpoint", - overwrite_existing = true) - -@info "Running with Δt = $(prettytime(simulation.Δt))" - -run!(simulation, pickup = pickup_file) - -@info """ - - Simulation took $(prettytime(simulation.run_wall_time)) - Free surface: $(typeof(model.free_surface).name.wrapper) - Time step: $(prettytime(Δt)) -""" diff --git a/validation/near_global_lat_lon/one_degree_setups/bathymetry-360x150-latitude-75.0.jld2 b/validation/near_global_lat_lon/one_degree_setups/bathymetry-360x150-latitude-75.0.jld2 deleted file mode 100644 index 6417fa0cbd..0000000000 Binary files a/validation/near_global_lat_lon/one_degree_setups/bathymetry-360x150-latitude-75.0.jld2 and /dev/null differ diff --git a/validation/near_global_lat_lon/one_degree_setups/idealized_near_global_one_degree.jl b/validation/near_global_lat_lon/one_degree_setups/idealized_near_global_one_degree.jl deleted file mode 100755 index e684bc1200..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/idealized_near_global_one_degree.jl +++ /dev/null @@ -1,197 +0,0 @@ -using Oceananigans -using Oceananigans.Units -using Oceananigans.Simulations: reset! -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom -using Printf -using JLD2 -using GLMakie -using SeawaterPolynomials.TEOS10 -using Oceananigans.Operators: Δx, Δy, ζ₃ᶠᶠᶜ - -using Oceananigans: fields -using Oceananigans.Grids: ynode -using Oceananigans.Utils: with_tracers -using Oceananigans.TurbulenceClosures: validate_closure - -arch = CPU() -filename = "idealized_near_global_one_degree" - -include("one_degree_artifacts.jl") -# bathymetry_path = download_bathymetry() # not needed because we uploaded to repo -bathymetry = jldopen(bathymetry_path)["bathymetry"] - -include("one_degree_interface_heights.jl") -z = one_degree_interface_heights() - -underlying_grid = LatitudeLongitudeGrid(arch; - size = (360, 150, 48), - halo = (4, 4, 4), - latitude = (-75, 75), - z, - longitude = (-180, 180), - precompute_metrics = true) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -@show grid - -##### -##### Closures -##### - -include("variable_biharmonic_diffusion_coefficient.jl") # defines VariableBiharmonicDiffusionCoefficient -vitd = VerticallyImplicitTimeDiscretization() -background_vertical_diffusivity = VerticalScalarDiffusivity(vitd, ν=1e-2, κ=1e-4) -dynamic_vertical_diffusivity = RiBasedVerticalDiffusivity() - -function idealized_one_degree_closure(; νh = (100kilometers)^2 / 1day, - κ_skew = 0, - κ_symmetric = κ_skew, - biharmonic_time_scale = 1day) - - ν₄ = VariableBiharmonicDiffusionCoefficient(biharmonic_time_scale) - biharmonic_viscosity = HorizontalScalarBiharmonicDiffusivity(ν=ν₄, discrete_form=true) - - gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew = κ_skew, - κ_symmetric = κ_symmetric, - slope_limiter = FluxTapering(1e-2)) - - horizontal_diffusivity = HorizontalScalarDiffusivity(ν=νh, κ=νh) - - closure_tuple = (gent_mcwilliams_diffusivity, - biharmonic_viscosity, - dynamic_vertical_diffusivity, - horizontal_diffusivity, - background_vertical_diffusivity) - - closure_tuple = with_tracers(tuple(:T), closure_tuple) - closure_tuple = validate_closure(closure_tuple) - - return closure_tuple -end - -free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver) -equation_of_state = LinearEquationOfState() -buoyancy = SeawaterBuoyancy(; equation_of_state, constant_salinity=35.0) - -@inline T_reference(φ) = max(0.0, 30.0 * cos(1.2 * π * φ / 180)) - -@inline function T_relaxation(i, j, grid, clock, fields, tᵣ) - φ = ynode(Center(), j, grid) - return @inbounds 1 / tᵣ * (fields.T[i, j, grid.Nz] - T_reference(φ)) -end - -T_top_bc = FluxBoundaryCondition(T_relaxation, discrete_form=true, parameters=30days) -T_bcs = FieldBoundaryConditions(top=T_top_bc) - -@inline surface_stress_x(φ, τ₀, τˢ, τᴺ) = τ₀ * (1 - exp(-φ^2 / 200)) - (τ₀ + τˢ) * exp(-(φ + 50)^2 / 200) - - (τ₀ + τᴺ) * exp(-(φ - 50)^2 / 200) - -@inline function surface_stress_x(i, j, grid, clock, fields, p) - φ = ynode(Center(), j, grid) - return surface_stress_x(φ, p.τ₀, p.τˢ, p.τᴺ) -end - -u_top_bc = FluxBoundaryCondition(surface_stress_x, discrete_form=true, parameters=(τ₀=6e-5, τˢ=2e-4, τᴺ=5e-5)) -u_bcs = FieldBoundaryConditions(top=u_top_bc) - -model = HydrostaticFreeSurfaceModel(; grid, free_surface, buoyancy, - momentum_advection = VectorInvariant(), - coriolis = HydrostaticSphericalCoriolis(), - tracers = :T, - closure = idealized_one_degree_closure(), - boundary_conditions = (u=u_bcs, T=T_bcs), - tracer_advection = WENO5(grid=underlying_grid)) - -@show model - -function T_initial(λ, φ, z) - H_stratification = 2000 - T_bottom = 0.0 - T_surface = T_reference(φ) - dTdz = (T_surface - T_bottom) / H_stratification - return max(T_bottom, T_surface + z * dTdz) -end - -set!(model, T=T_initial) - -simulation = Simulation(model; Δt=10minutes, stop_time=10years) - -wall_clock = Ref(time_ns()) -function progress(sim) - elapsed = 1e-9 * (time_ns() - wall_clock[]) - - u, v, w = sim.model.velocities - T = sim.model.tracers.T - η = sim.model.free_surface.η - umax = maximum(abs, u) - vmax = maximum(abs, v) - wmax = maximum(abs, w) - ηmax = maximum(abs, η) - Tmax = maximum(T) - Tmin = minimum(T) - - msg1 = @sprintf("Iteration: %d, time: %s, wall time: %s", iteration(sim), prettytime(sim), prettytime(elapsed)) - msg2 = @sprintf("├── max(u): (%.2e, %.2e, %.2e) m s⁻¹", umax, vmax, wmax) - msg3 = @sprintf("├── extrema(T): (%.2f, %.2f) ᵒC", Tmin, Tmax) - msg4 = @sprintf("└── max|η|: %.2e m", ηmax) - - @info string(msg1, "\n", msg2, "\n", msg3, "\n", msg4) - - wall_clock[] = time_ns() - - return nothing -end - -simulation.callbacks[:p] = Callback(progress, IterationInterval(10)) - -## Spin up simulation then reset. -#run!(simulation) -#reset!(simulation) -#simulation.stop_time = 30days -#model.closure = idealized_one_degree_closure(κ_skew=1e3) - -u, v, w = model.velocities -KE = @at (Center, Center, Center) u^2 + v^2 -ζ = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid, u, v) - -simulation.output_writers[:surface] = JLD2OutputWriter(model, merge(fields(model), (; KE, ζ)), - schedule = TimeInterval(1day), - filename = filename * "_surface.jld2", - indices = (:, :, grid.Nz), - overwrite_existing = true) - -run!(simulation) - -ut = FieldTimeSeries(filename * "_surface.jld2", "u") -vt = FieldTimeSeries(filename * "_surface.jld2", "v") -Tt = FieldTimeSeries(filename * "_surface.jld2", "T") -Kt = FieldTimeSeries(filename * "_surface.jld2", "KE") -Zt = FieldTimeSeries(filename * "_surface.jld2", "ζ") -t = ut.times -Nt = length(t) - -fig = Figure(resolution=(1800, 900)) -axk = Axis(fig[2, 2], xlabel="Longitude", ylabel="Latitude") -axz = Axis(fig[2, 3], xlabel="Longitude", ylabel="Latitude") -slider = Slider(fig[3, 2:3], range=1:Nt, startvalue=1) -n = slider.value - -title = @lift string("Surface kinetic energy at ", prettytime(t[$n])) -Label(fig[1, 1], title, tellwidth=false) - -KEⁿ = @lift interior(Kt[$n], :, :, 1) ./ 2 -ζⁿ = @lift interior(Zt[$n], :, :, 1) - -hmk = heatmap!(axk, KEⁿ, colorrange=(0, 0.1)) -hmz = heatmap!(axz, ζⁿ, colorrange=(-1e-4, 1e-4)) - -Colorbar(fig[2, 1], hmk, label="Surface kinetic energy (m² s⁻²)") -Colorbar(fig[2, 4], hmz, label="Surface vertical vorticity (s⁻¹)") - -display(fig) - -record(fig, filename * ".mp4", 1:Nt, framerate=12) do nn - n[] = nn -end - diff --git a/validation/near_global_lat_lon/one_degree_setups/one_degree_artifacts.jl b/validation/near_global_lat_lon/one_degree_setups/one_degree_artifacts.jl deleted file mode 100644 index e0b46a9d33..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/one_degree_artifacts.jl +++ /dev/null @@ -1,17 +0,0 @@ -using Downloads - -bathymetry_path = joinpath(@__DIR__, "bathymetry-360x150-latitude-75.0.jld2") -boundary_conditions_path = joinpath(@__DIR__, "boundary_conditions-1degree.jld2") -initial_condition_path = joinpath(@__DIR__, "initial_conditions-1degree.jld2") - -# TODO: convert to DataDeps - -download_bathymetry(path=bathymetry_path) = - Downloads.download("https://www.dropbox.com/s/axyzt88g0nr9dbc/bathymetry-360x150-latitude-75.0.jld2", path) - -download_boundary_conditions(path=boundary_conditions_path) = - Downloads.download("https://www.dropbox.com/s/7sbrq5fcnuhtvqp/boundary_conditions-1degree.jld2", path) - -download_initial_condition(path=initial_condition_path) = - Downloads.download("https://www.dropbox.com/s/hbbzeg301hwm2xd/initial_conditions-1degree.jld2", path) - diff --git a/validation/near_global_lat_lon/one_degree_setups/one_degree_interface_heights.jl b/validation/near_global_lat_lon/one_degree_setups/one_degree_interface_heights.jl deleted file mode 100644 index b71e508f5d..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/one_degree_interface_heights.jl +++ /dev/null @@ -1,51 +0,0 @@ -one_degree_interface_heights() = - [-5244.5 - -4834.0 - -4446.5 - -4082.0 - -3740.5 - -3422.0 - -3126.5 - -2854.0 - -2604.5 - -2378.0 - -2174.4 - -1993.6 - -1834.7 - -1695.6 - -1572.8 - -1461.4 - -1356.9 - -1255.5 - -1155.5 - -1056.3 - -958.0 - -861.4 - -767.5 - -677.3 - -592.2 - -513.3 - -441.7 - -378.2 - -323.2 - -276.7 - -238.3 - -207.2 - -182.3 - -162.5 - -146.5 - -133.0 - -121.3 - -110.5 - -100.1 - -90.0 - -80.0 - -70.0 - -60.0 - -50.0 - -40.0 - -30.0 - -20.0 - -10.0 - 0.0] - diff --git a/validation/near_global_lat_lon/one_degree_setups/plot_bathymetry.jl b/validation/near_global_lat_lon/one_degree_setups/plot_bathymetry.jl deleted file mode 100644 index 15195ecbc3..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/plot_bathymetry.jl +++ /dev/null @@ -1,15 +0,0 @@ -using GLMakie -using JLD2 - -filename = "bathymetry-360x150-latitude-75.0.jld2" -file = jldopen(filename) -h = file["bathymetry"] -close(file) -h[h .> 0] .= NaN - -fig = Figure(resolution=(1800, 800)) -ax = Axis(fig[1, 1], xlabel="Longitude", ylabel="Latitude") -hm = heatmap!(ax, h) -Colorbar(fig[1, 2], hm, label="Elevation (m)") -display(fig) - diff --git a/validation/near_global_lat_lon/one_degree_setups/plot_boundary_conditions.jl b/validation/near_global_lat_lon/one_degree_setups/plot_boundary_conditions.jl deleted file mode 100644 index 6bcae49c86..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/plot_boundary_conditions.jl +++ /dev/null @@ -1,103 +0,0 @@ -using JLD2 -using GLMakie - -include("one_degree_artifacts.jl") -# bathymetry_path = download_bathymetry() # not needed because we uploaded to repo -boundary_conditions_path = download_boundary_conditions() - -file = jldopen(bathymetry_path) -h = file["bathymetry"] -land = h .> 0 -wet_Nx = sum(h -> h < 0, h, dims=1) -close(file) - -ρ₀ = 1025 -φ = -74.5:74.5 -file = jldopen(boundary_conditions_path) - -τˣ = file["τˣ"] ./ ρ₀ -τʸ = file["τʸ"] ./ ρ₀ -Tˢ = file["Tˢ"] -Sˢ = file["Sˢ"] - -close(file) - -τˣ[isnan.(τˣ)] .= 0 -τʸ[isnan.(τʸ)] .= 0 - -τlim = 0.8 * max(maximum(abs, τˣ), maximum(abs, τʸ)) - -for n = 1:12 - view(τˣ, :, :, n)[land] .= NaN - view(τʸ, :, :, n)[land] .= NaN - view(Tˢ, :, :, n)[land] .= NaN - view(Sˢ, :, :, n)[land] .= NaN -end - -fig = Figure(resolution=(2000, 900)) -axx = Axis(fig[2, 2], xlabel="Longitude", ylabel="Latitude", title="Zonal surface stress") -axy = Axis(fig[2, 3], xlabel="Longitude", ylabel="Latitude", title="Meridonal surface stress") -axT = Axis(fig[3, 2], xlabel="Longitude", ylabel="Latitude") -axS = Axis(fig[3, 3], xlabel="Longitude", ylabel="Latitude") -slider = Slider(fig[4, 2:3], range=1:12, startvalue=1) -n = slider.value - -title = @lift string("Boundary conditions in month ", $n) -Label(fig[1, 2:3], title) - -τˣn = @lift view(τˣ, :, :, $n) -τʸn = @lift view(τʸ, :, :, $n) -Tn = @lift view(Tˢ, :, :, $n) -Sn = @lift view(Sˢ, :, :, $n) - -hmx = heatmap!(axx, τˣn, colorrange=(-τlim, τlim), colormap=:redblue) -hmy = heatmap!(axy, τʸn, colorrange=(-τlim, τlim), colormap=:redblue) -hmT = heatmap!(axT, Tn, colorrange=(-1, 31), colormap=:thermal) -hmS = heatmap!(axS, Sn, colorrange=(25, 38), colormap=:haline) - -Colorbar(fig[2, 4], hmx, label="Stress (m² s⁻²)") -Colorbar(fig[3, 1], hmT, label="Surface temperature (ᵒC)", flipaxis=false) -Colorbar(fig[3, 4], hmS, label="Surface salinity (psu)") - -display(fig) - -record(fig, "one_degree_boundary_conditions.mp4", 1:12, framerate=1) do nn - n[] = nn -end - -##### -##### Zonal means -##### - -filternan(t) = ifelse(isnan(t), 0.0, t) -zonal_mean_τˣ = sum(filternan, τˣ, dims=1) ./ wet_Nx -zonal_mean_τʸ = sum(filternan, τʸ, dims=1) ./ wet_Nx -zonal_mean_Tˢ = sum(filternan, Tˢ, dims=1) ./ wet_Nx -zonal_mean_Sˢ = sum(filternan, Sˢ, dims=1) ./ wet_Nx - -zonal_mean_τˣ = dropdims(zonal_mean_τˣ, dims=1)' -zonal_mean_τʸ = dropdims(zonal_mean_τʸ, dims=1)' -zonal_mean_Tˢ = dropdims(zonal_mean_Tˢ, dims=1)' -zonal_mean_Sˢ = dropdims(zonal_mean_Sˢ, dims=1)' - -fig = Figure(resolution=(2000, 900)) - -axx = Axis(fig[1, 2], xlabel="Month", ylabel="Latitude", title="Zonally-averaged zonal surface stress") -axy = Axis(fig[1, 3], xlabel="Month", ylabel="Latitude", title="Zonally-averaged meridonal surface stress") -axT = Axis(fig[2, 2], xlabel="Month", ylabel="Latitude") -axS = Axis(fig[2, 3], xlabel="Month", ylabel="Latitude") - -τlim = 1.5e-4 -hmx = heatmap!(axx, zonal_mean_τˣ, colorrange=(-τlim, τlim), colormap=:redblue) -hmy = heatmap!(axy, zonal_mean_τʸ, colorrange=(-τlim, τlim), colormap=:redblue) -hmT = heatmap!(axT, zonal_mean_Tˢ, colorrange=(-1, 31), colormap=:thermal) -hmS = heatmap!(axS, zonal_mean_Sˢ, colorrange=(25, 38), colormap=:haline) - -Colorbar(fig[1, 4], hmx, label="Stress (m² s⁻²)") -Colorbar(fig[2, 1], hmT, label="Surface temperature (ᵒC)", flipaxis=false) -Colorbar(fig[2, 4], hmS, label="Surface salinity (psu)") - -display(fig) - -save("zonally_averaged_one_degree_boundary_conditions.png", fig) - diff --git a/validation/near_global_lat_lon/one_degree_setups/plot_initial_condition.jl b/validation/near_global_lat_lon/one_degree_setups/plot_initial_condition.jl deleted file mode 100644 index fdc456003e..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/plot_initial_condition.jl +++ /dev/null @@ -1,44 +0,0 @@ -using JLD2 -using GLMakie - -include("one_degree_artifacts.jl") -include("one_degree_interface_heights.jl") -z_interfaces = one_degree_interface_heights() -z_centers = (z_interfaces[1:end-1] .+ z_interfaces[2:end]) ./ 2 - -filename = "bathymetry-360x150-latitude-75.0.jld2" -file = jldopen(filename) -h = file["bathymetry"] -close(file) -land = h .> 0 - -initial_condition_path = download_initial_condition() -file = jldopen(initial_condition_path) -Tᵢ = file["T"] -Nz = size(Tᵢ, 3) - -for k in 1:Nz - z = z_centers[k] - earth = h .> z - view(Tᵢ, :, :, k)[earth] .= NaN -end - -fig = Figure(resolution=(1800, 800)) -ax = Axis(fig[2, 1]) -slider = Slider(fig[3, 1:2], range=1:Nz, startvalue=Nz) -k = slider.value - -depth = @lift z_centers[$k] -title = @lift string("Temperature (ᵒC) at z = ", $depth, " m") -Label(fig[1, 1:2], title) - -Tᵏ = @lift view(Tᵢ, :, :, $k) - -hm = heatmap!(ax, Tᵏ, colorrange=(-2, 32), colormap=:thermal) -Colorbar(fig[2, 2], hm, label="Temperature (ᵒC)") - -display(fig) - -record(fig, "one_degree_initial_condition.mp4", 1:Nz, framerate=12) do kk - k[] = kk -end diff --git a/validation/near_global_lat_lon/one_degree_setups/variable_biharmonic_diffusion_coefficient.jl b/validation/near_global_lat_lon/one_degree_setups/variable_biharmonic_diffusion_coefficient.jl deleted file mode 100644 index 43b4c7ce7a..0000000000 --- a/validation/near_global_lat_lon/one_degree_setups/variable_biharmonic_diffusion_coefficient.jl +++ /dev/null @@ -1,15 +0,0 @@ -using Oceananigans.Utils: prettysummary - -struct VariableBiharmonicDiffusionCoefficient - time_scale :: Float64 -end - -@inline function (vbd::VariableBiharmonicDiffusionCoefficient)(i, j, k, grid, lx, ly, lz) - Δh⁻² = 1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2 - Δh⁴ = 1 / Δh⁻²^2 - return Δh⁴ / vbd.time_scale -end - -Base.summary(vbdc::VariableBiharmonicDiffusionCoefficient) = string("VariableBiharmonicDiffusionCoefficient(time_scale=", - prettytime(vbdc.time_scale), ")") - diff --git a/validation/near_global_lat_lon/visualize.jl b/validation/near_global_lat_lon/visualize.jl deleted file mode 100644 index 647c1f9120..0000000000 --- a/validation/near_global_lat_lon/visualize.jl +++ /dev/null @@ -1,43 +0,0 @@ -using Oceananigans.Grids -using Oceananigans.Utils: prettytime, hours, day, days, years - -using Statistics -using JLD2 -using Printf -using CairoMakie - -output_prefix = "near_global_lat_lon_1440_600_1_fine_surface" -filepath = output_prefix * ".jld2" -file = jldopen(filepath) - -Nx = file["grid/underlying_grid/Nx"] -Ny = file["grid/underlying_grid/Ny"] -Lλ = file["grid/underlying_grid/Lx"] -Lφ = file["grid/underlying_grid/Ly"] -Lz = file["grid/underlying_grid/Lz"] - -grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1), - longitude = (-180, 180), - latitude = (-Lφ/2, Lφ/2), - z = (-Lz, 0)) - -x, y, z = nodes((Center, Center, Center), grid) - -bottom = Float32.(file["grid/immersed_boundary/bottom_height"][4:end-3,3:end-3,1]) -bottom[ bottom .>0 ] .= NaN -bottom[ bottom .<0 ] .= 0.0 - -iter = Observable(0) -iters = parse.(Int, keys(file["timeseries/t"])) -ζ′ = @lift file["timeseries/ζ/" * string($iter)][:, :, 1] -title = @lift(@sprintf("Surface Vorticity in Hydrostatic Model at time = %s", prettytime(file["timeseries/t/" * string($iter)]))) -fig = Figure(resolution = (2000, 1000)) -ax = Axis(fig[1,1], xlabel = "longitude", ylabel = "latitude", title=title) -heatmap_plot = heatmap!(ax, ζ′, colormap=:balance, colorrange=(-1e-6, 1e-6), nan_color=:black) -Colorbar(fig[1,2], heatmap_plot , width=25) -display(fig) - -record(fig, output_prefix * ".mp4", iters[2:end], framerate=6) do i - @info "Plotting iteration $i of $(iters[end])..." - iter[] = i -end diff --git a/validation/near_global_lat_lon/visualize_annual_cycle_solution.jl b/validation/near_global_lat_lon/visualize_annual_cycle_solution.jl deleted file mode 100644 index 065adfc782..0000000000 --- a/validation/near_global_lat_lon/visualize_annual_cycle_solution.jl +++ /dev/null @@ -1,301 +0,0 @@ -using Statistics -using JLD2 -using Printf -using GLMakie -using Oceananigans -using Oceananigans.Units -using Oceananigans.Fields: fill_halo_regions!, ReducedField - -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom -using CUDA: @allowscalar - -##### -##### Grid -##### - -latitude = (-84.375, 84.375) -Δφ = latitude[2] - latitude[1] - -# 2.8125 degree resolution -Nx = 128 -Ny = 60 -Nz = 18 - -using DataDeps - -path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/" - -dh = DataDep("near_global_lat_lon", - "Forcing data for global latitude longitude simulation", path * "bathymetry_lat_lon_128x60_FP32.bin" -) - -DataDeps.register(dh) - -datadep"near_global_lat_lon" - -bathymetry_data = Array{Float32}(undef, Nx*Ny) -bathymetry_path = @datadep_str "near_global_lat_lon/bathymetry_lat_lon_128x60_FP32.bin" -read!(bathymetry_path, bathymetry_data) - -bathymetry_data = bswap.(bathymetry_data) |> Array{Float64} -bathymetry = reshape(bathymetry_data, Nx, Ny) - -Nmonths = 12 -bytes = sizeof(Float32) * Nx * Ny - -H = 3600.0 -#bathymetry = - H .* (bathymetry .< -10) - -# A spherical domain -underlying_grid = LatitudeLongitudeGrid(size = (Nx, Ny, Nz), - longitude = (-180, 180), - latitude = latitude, - halo = (3, 3, 3), - z = (-H, 0)) - -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry)) - -function geographic2cartesian(λ, φ; r=1) - λ = cat(λ, λ[1:1] .+ 360, dims=1) - - Nλ = length(λ) - Nφ = length(φ) - - λ = repeat(reshape(λ, Nλ, 1), 1, Nφ) - φ = repeat(reshape(φ, 1, Nφ), Nλ, 1) - - λ_azimuthal = λ .+ 180 # Convert to λ ∈ [0°, 360°] - φ_azimuthal = 90 .- φ # Convert to φ ∈ [0°, 180°] (0° at north pole) - - x = @. r * cosd(λ_azimuthal) * sind(φ_azimuthal) - y = @. r * sind(λ_azimuthal) * sind(φ_azimuthal) - z = @. r * cosd(φ_azimuthal) - - return x, y, z -end - -λu, ϕu, ru = nodes((Face, Center, Center), grid) -λv, ϕv, rv = nodes((Center, Face, Center), grid) -λc, ϕc, rc = nodes((Center, Center, Center), grid) - -xu, yu, zu = geographic2cartesian(λu, ϕu) -xv, yv, zv = geographic2cartesian(λv, ϕv) -xc, yc, zc = geographic2cartesian(λc, ϕc) - -output_prefix = "annual_cycle_global_lat_lon_128_60_18_temp" - -surface_file = jldopen(output_prefix * "_surface.jld2") -bottom_file = jldopen(output_prefix * "_bottom.jld2") - -iterations = parse.(Int, keys(surface_file["timeseries/t"])) - -iter = Node(0) - -# Continents -land_ccc = bathymetry .> - underlying_grid.Δzᵃᵃᶜ / 2 - -function mask_land_cfc(v) - land_cfc = cat(land_ccc, land_ccc[:, end:end], dims=2) - v[land_cfc] .= NaN - return v -end - -function mask_land(data) - data[land_ccc] .= NaN - return data -end - -ηi(iter) = mask_land(surface_file["timeseries/η/" * string(iter)][:, :, 1]) -ui(iter) = mask_land(surface_file["timeseries/u/" * string(iter)][:, :, 1]) -vi(iter) = mask_land_cfc(surface_file["timeseries/v/" * string(iter)][:, :, 1]) -Ti(iter) = mask_land(surface_file["timeseries/T/" * string(iter)][:, :, 1]) -ti(iter) = prettytime(surface_file["timeseries/t/" * string(iter)]) - -ubi(iter) = mask_land(bottom_file["timeseries/u/" * string(iter)][:, :, 1]) -vbi(iter) = mask_land_cfc(bottom_file["timeseries/v/" * string(iter)][:, :, 1]) - -uri = ReducedField(Face, Center, Nothing, CPU(), grid, dims=3) -vri = ReducedField(Center, Face, Nothing, CPU(), grid, dims=3) -spi = ReducedField(Center, Center, Nothing, CPU(), grid, dims=3) - -bathymetry = cat(bathymetry, bathymetry[1:1, :], dims=1) -land = bathymetry .> -10 -water = @. !(land) - -continents = 1.0 .* land -continents[water] .= NaN - -sp = @lift begin - uri .= ui($iter) - vri .= vi($iter) - fill_halo_regions!(uri) - fill_halo_regions!(vri) - spi .= sqrt(uri^2 + vri^2) - fill_halo_regions!(spi) - sp = spi[1:Nx+1, 1:Ny, 1] - # sp[land] .= NaN - return sp -end - -η = @lift ηi($iter) -u = @lift ui($iter) -v = @lift vi($iter) -T = @lift Ti($iter) - -T_sphere = @lift begin - T = cat(Ti($iter), Ti($iter)[1:1, :], dims=1) - T[land] .= NaN - return T -end - -v_sphere = @lift begin - v = cat(vi($iter)[:, 1:Ny], vi($iter)[1:1, 1:Ny], dims=1) - v[land] .= NaN - return v -end - -u_sphere = @lift begin - u = cat(ui($iter), ui($iter)[1:1, :], dims=1) - u[land] .= NaN - return u -end - -ub = @lift ubi($iter) -vb = @lift vbi($iter) - -max_η = 2 -min_η = - max_η -max_u = 0.1 -min_u = - max_u -max_T = 32 -min_T = 0 - -dλbg = 0.01 -dφbg = 0.01 -λbg = range(-180, stop=180 - dλbg, step=dλbg) -φbg = range(-89.5, stop=89.5 - dφbg, step=dφbg) -xbg, ybg, zbg = geographic2cartesian(λbg, φbg) - -sphere_continents!(ax_T) = surface!(ax_T, xbg, ybg, zbg, color=fill("#080", 1, 1)) - -##### -##### Sphere -##### - -sphere_fig = Figure(resolution = (150, 90)) - -graylim = 1.0 -α = 0.99 - -#surface!(ax, xc, yc, zc, color=sp, colormap=:blues, colorrange=(0, max_u), -# show_axis=false, shading=false, ssao=true) - -# Temperature -ax_T = sphere_fig[1, 1:3] = LScene(sphere_fig) - -surface!(ax_T, xc, yc, zc, color=T_sphere, colormap=:thermal, colorrange=(0, max_T), - show_axis=false, shading=false, ssao=true) - -sphere_continents!(ax_T) - -rotate_cam!(ax_T.scene, (π/8, π/6, 0)) - -# Meridional velocity -ax_v = sphere_fig[1, 4:6] = LScene(sphere_fig) - -surface!(ax_v, xc, yc, zc, color=u_sphere, colormap=:balance, colorrange=(min_u, max_u), - show_axis=false, shading=false, ssao=true) - -sphere_continents!(ax_v) - -rotate_cam!(ax_v.scene, (π/8, π/6, 0)) - -plot_title = @lift "t = " * ti($iter) -supertitle = sphere_fig[0, 3:4] = Label(sphere_fig, plot_title, fontsize=40) - -#v_title = sphere_fig[0, 5] = Label(sphere_fig, "Ocean meridional velocity", fontsize=20) -#T_title = sphere_fig[0, 2] = Label(sphere_fig, "Ocean temperature", fontsize=20) - -display(sphere_fig) - -save_interval = 5days -full_rotation_savepoints = round(Int, 10year / save_interval) -dθ = 2π / full_rotation_savepoints -dϕ = π/4 / full_rotation_savepoints - -function dϕi(i) - savepoint = i/iterations[end] * length(iterations) - mod_savepoint = mod1(savepoint, full_rotation_savepoints) - halfway_savepoint = 1/2 * length(iterations) - return mod_i < halfway_savepoint ? dϕ : - dϕ -end - -GLMakie.record(sphere_fig, output_prefix * "_sphere.mp4", iterations, framerate=8) do i - @info "Plotting iteration $i of $(iterations[end])..." - iter[] = i - rotate_cam!(ax_T.scene, cameracontrols(ax_T.scene), (dϕi(i), dθ, 0)) - rotate_cam!(ax_v.scene, cameracontrols(ax_v.scene), (dϕi(i), dθ, 0)) -end - -#= -##### -##### Meridional velocity -##### - -plane_fig = Figure(resolution = (1200, 900)) - -# kwargs = NamedTuple() -kwargs = (; interpolate=true) - -function continents!(ax) - Nλ = length(λc) - Nφ = length(ϕc) - λ = repeat(reshape(λc, Nλ, 1), 1, Nφ) - φ = repeat(reshape(ϕc, 1, Nφ), Nλ, 1) - surface!(ax, λ, φ, zeros(Nx, Ny) .- 0.02, color=fill("#6e7f80", 1, 1)) - return nothing -end - -ax = Axis(plane_fig[1, 1], title="Free surface displacement (m)") -continents!(ax) -hm = heatmap!(ax, λc, ϕc, η; colorrange=(min_η, max_η), colormap=:balance, kwargs...) -cb = Colorbar(plane_fig[1, 2], hm) - -ax = Axis(plane_fig[2, 1], title="Sea surface temperature (ᵒC)") -continents!(ax) -hm = heatmap!(ax, λc, ϕc, T; colorrange=(min_T, max_T), colormap=:thermal, kwargs...) -cb = Colorbar(plane_fig[2, 2], hm) - -ax = Axis(plane_fig[1, 3], title="East-west surface velocity (m s⁻¹)") -continents!(ax) -hm = heatmap!(ax, λc, ϕc, u; colorrange=(min_u, max_u), colormap=:balance, kwargs...) -cb = Colorbar(plane_fig[1, 4], hm) - -ax = Axis(plane_fig[2, 3], title="North-south surface velocity (m s⁻¹)") -continents!(ax) -hm = heatmap!(ax, λc, ϕv, v; colorrange=(min_u, max_u), colormap=:balance, kwargs...) -cb = Colorbar(plane_fig[2, 4], hm) - -ax = Axis(plane_fig[3, 1], title="East-west bottom velocity (m s⁻¹)") -continents!(ax) -hm = heatmap!(ax, λc, ϕc, ub; colorrange=(min_u, max_u), colormap=:balance, kwargs...) -cb = Colorbar(plane_fig[3, 2], hm) - -ax = Axis(plane_fig[3, 3], title="North-south bottom velocity (m s⁻¹)") -continents!(ax) -hm = heatmap!(ax, λc, ϕv, vb; colorrange=(min_u, max_u), colormap=:balance, kwargs...) -cb = Colorbar(plane_fig[3, 4], hm) - -title_str = @lift "Earth day = " * ti($iter) -ax_t = plane_fig[0, :] = Label(plane_fig, title_str) - -GLMakie.record(plane_fig, output_prefix * "_plane.mp4", iterations, framerate=8) do i - @info "Plotting iteration $i of $(iterations[end])..." - iter[] = i -end - -display(plane_fig) -=# - -close(surface_file) -close(bottom_file)