Merge pull request #131 from CliMA/glw/infinite-impulsive-storm

Stokes drift for constant momentum fluxes

Project.toml

@@ -7,9 +7,6 @@ version = "0.1.0"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
-GeoData = "9b6fcbb8-86d6-11e9-1ce7-23a6bb139a78"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
@@ -20,15 +17,15 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 Oceanostics = "d0ccf422-c8fb-49b5-a76d-74acdde946ac"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Oceananigans = "^0.68.2"
-Oceanostics = "^0.6.3"
 julia = "^1.6"
 
 [extras]

examples/actually_run_three_layer_constant_fluxes.jl

@@ -0,0 +1,109 @@
+# Example script for generating and analyzing LESbrary data
+
+using Oceananigans
+using Oceananigans.Units
+
+using CUDA
+
+using LESbrary.IdealizedExperiments: three_layer_constant_fluxes_simulation
+using LESbrary.IdealizedExperiments: six_hour_suite_parameters
+using LESbrary.IdealizedExperiments: eighteen_hour_suite_parameters
+using LESbrary.IdealizedExperiments: twelve_hour_suite_parameters
+using LESbrary.IdealizedExperiments: twenty_four_hour_suite_parameters
+using LESbrary.IdealizedExperiments: thirty_six_hour_suite_parameters
+using LESbrary.IdealizedExperiments: forty_eight_hour_suite_parameters
+using LESbrary.IdealizedExperiments: seventy_two_hour_suite_parameters
+
+# LESbrary parameters
+# ===================
+#
+# Canonical LESbrary parameters are organized by the _intended duration_ of the simulation
+# --- two days for the "two day suite", four days for the "four day suite", et cetera.
+# The strength of the forcing is tuned roughly so that the boundary layer deepens to roughly
+# 128 meters, half the depth of a 256 meter domain.
+#
+# Each suite has five cases:
+#
+# * :free_convection
+# * :weak_wind_strong_cooling
+# * :strong_wind_weak_cooling
+# * :strong_wind
+# * :strong_wind_no_rotation
+#
+# In addition to selecting the architecture, size, and case to run, we can also tweak
+# certain parameters below.
+
+configuration = (;
+    architecture = GPU(),
+    size = (448, 448, 512),
+    #size = (384, 384, 384),
+    #size = (64, 64, 64),
+    snapshot_time_interval = 10minutes,
+    passive_tracers = false,
+    jld2_output = true,
+    checkpoint = true,
+    time_averaged_statistics = true,
+    data_directory = "/nobackup/users/glwagner"
+)
+
+cases = [
+         :free_convection,
+         :weak_wind_strong_cooling, 
+         :med_wind_med_cooling,     
+         :strong_wind_weak_cooling, 
+         :strong_wind,
+         :strong_wind_no_rotation,
+        ]
+
+#for case in cases
+
+#case = :free_convection
+#case = :strong_wind_no_rotation
+#case = :weak_wind_strong_cooling
+#case = :med_wind_med_cooling
+#case = :strong_wind_weak_cooling
+case = :strong_wind
+
+    parameters = six_hour_suite_parameters[case]
+    #parameters = twelve_hour_suite_parameters[case]
+    #parameters = eighteen_hour_suite_parameters[case]
+    #parameters = seventy_two_hour_suite_parameters[case]
+
+    println("Running $case")
+    for (k, v) in parameters
+        println("    - ", k, ": ", v)
+    end
+
+    simulation = three_layer_constant_fluxes_simulation(; configuration..., parameters...)
+
+    if !isnothing(simulation.model.stokes_drift)
+        stokes_drift = simulation.model.stokes_drift
+        grid = simulation.model.grid
+        @show uˢ₀ = CUDA.@allowscalar stokes_drift.uˢ[1, 1, grid.Nz]
+
+        a★ = stokes_drift.air_friction_velocity
+        ρʷ = stokes_drift.water_density
+        ρᵃ = stokes_drift.air_density
+        u★ = a★ * sqrt(ρᵃ / ρʷ)
+        @show La = sqrt(u★ / uˢ₀)
+    end
+
+    run!(simulation)
+#end
+
+#=
+#case = :free_convection
+#case = :weak_wind_strong_cooling
+#case = :med_wind_med_cooling
+#case = :strong_wind_weak_cooling
+#case = :strong_wind
+case = :strong_wind_no_rotation
+
+parameters = thirty_six_hour_suite_parameters[case]
+@info "Running case $case with $parameters..."
+
+start_time = time_ns()
+simulation = three_layer_constant_fluxes_simulation(; configuration..., parameters...)
+run!(simulation)
+elapsed = 1e-9 * (time_ns() - start_time)
+=#

examples/free_convection.jl

@@ -3,61 +3,41 @@
 # This script runs a simulation of convection driven by cooling at the
 # surface of an idealized, stratified, rotating ocean surface boundary layer.
 
-using LESbrary, Printf, Statistics, Oceananigans, Oceananigans.Units
+using LESbrary
+using Printf
+using Statistics
+using Oceananigans
+using Oceananigans.Units
+using LESbrary.Utils: SimulationProgressMessenger
 
 # Domain
 
-grid = RectilinearGrid(CPU(), size=(32, 32, 32), x=(0, 128), y=(0, 128), z=(-64, 0))
+Nx = Ny = Nz = 128
+grid = RectilinearGrid(GPU(), size=(Nx, Ny, Nz), x=(0, 512), y=(0, 512), z=(-256, 0))
 
 # Buoyancy and boundary conditions
 
 Qᵇ = 1e-7
 N² = 1e-5
 
-buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4), constant_salinity=35.0)
-
-α = buoyancy.equation_of_state.α
-g = buoyancy.gravitational_acceleration
-
-## Compute temperature flux and gradient from buoyancy flux and gradient
-Qᵀ = Qᵇ / (α * g)
-dTdz = N² / (α * g)
-
-T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ),
-                                bottom = GradientBoundaryCondition(dTdz))
-
-# LES Model
-
-# Wall-aware AMD model constant which is 'enhanced' near the upper boundary.
-# Necessary to obtain a smooth temperature distribution.
-using LESbrary.NearSurfaceTurbulenceModels: SurfaceEnhancedModelConstant
-
-Δz = grid.Δzᵃᵃᶜ
-Cᴬᴹᴰ = SurfaceEnhancedModelConstant(Δz, C₀ = 1/12, enhancement = 7, decay_scale = 4Δz)
-
-# Instantiate Oceananigans.NonhydrostaticModel
-
-using Oceananigans
-using Oceananigans.Advection: WENO5
+b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))
 
-model = NonhydrostaticModel(; grid, buoyancy,
+model = NonhydrostaticModel(; grid,
                             timestepper = :RungeKutta3,
-                            advection = WENO5(),
-                            tracers = :T,
-                            coriolis = FPlane(f=1e-4),
-                            closure = AnisotropicMinimumDissipation(C=Cᴬᴹᴰ),
-                            boundary_conditions = (; T=T_bcs))
+                            buoyancy = BuoyancyTracer(),
+                            advection = WENO(),
+                            tracers = :b,
+                            closure = AnisotropicMinimumDissipation(),
+                            boundary_conditions = (; b=b_bcs))
 
 # # Initial condition
 
 Ξ(z) = rand() * exp(z / 8)
-Tᵢ(x, y, z) = dTdz * z + 1e-6 * Ξ(z) * dTdz * grid.Lz
-set!(model, T=Tᵢ)
+bᵢ(x, y, z) = N² * z + 1e-6 * Ξ(z) * N² * grid.Lz
+set!(model, b=bᵢ)
 
 # # Prepare the simulation
 
-using LESbrary.Utils: SimulationProgressMessenger
-
 # Adaptive time-stepping
 
 simulation = Simulation(model, Δt=2.0, stop_time=8hour)
@@ -68,31 +48,18 @@ simulation.callbacks[:progress] = Callback(SimulationProgressMessenger(wizard),
 
 # Prepare Output
 
-using Oceananigans.Utils: GiB
-using Oceananigans.OutputWriters
-
-prefix = @sprintf("free_convection_Qb%.1e_Nsq%.1e_Nh%d_Nz%d", Qᵇ, N², grid.Nx, grid.Nz)
-
-data_directory = joinpath(@__DIR__, "..", "data", prefix) # save data in /data/prefix
-
-# Copy this file into the directory with data
-mkpath(data_directory)
-cp(@__FILE__, joinpath(data_directory, basename(@__FILE__)), force=true)
+filename = @sprintf("free_convection_Qb%.1e_Nsq%.1e_Nh%d_Nz%d", Qᵇ, N², grid.Nx, grid.Nz)
 
 simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                                                       schedule = TimeInterval(4hour),
-                                                      prefix = prefix * "_fields",
-                                                      dir = data_directory,
-                                                      max_filesize = 2GiB,
-                                                      force = true)
+                                                      filename = filename * "_fields",
+                                                      overwrite_existing = true)
 
 simulation.output_writers[:slices] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                       schedule = AveragedTimeInterval(1hour, window=15minute),
-                                                      prefix = prefix * "_slices",
-                                                      field_slicer = FieldSlicer(j=floor(Int, grid.Ny/2)),
-                                                      dir = data_directory,
-                                                      max_filesize = 2GiB,
-                                                      force = true)
+                                                      filename = filename * "_slices",
+                                                      indices = (:, floor(Int, grid.Ny/2), :),
+                                                      overwrite_existing = true)
 
 # Horizontally-averaged turbulence statistics
 turbulence_statistics = LESbrary.TurbulenceStatistics.first_through_second_order(model)
@@ -101,100 +68,10 @@ tke_budget_statistics = LESbrary.TurbulenceStatistics.turbulent_kinetic_energy_b
 simulation.output_writers[:statistics] =
     JLD2OutputWriter(model, merge(turbulence_statistics, tke_budget_statistics),
                      schedule = AveragedTimeInterval(1hour, window=15minute),
-                     prefix = prefix * "_statistics",
-                     dir = data_directory,
-                     force = true)
+                     filename = filename * "_statistics",
+                     overwrite_existing = true)
 
 # # Run
 
-LESbrary.Utils.print_banner(simulation)
-
 run!(simulation)
 
-# # Load and plot turbulence statistics
-
-using JLD2, Plots
-
-## Some plot parameters
-linewidth = 3
-ylim = (-64, 0)
-plot_size = (1000, 500)
-zC = znodes(Center, grid)
-zF = znodes(Face, grid)
-
-## Load data
-file = jldopen(simulation.output_writers[:statistics].filepath)
-
-iterations = parse.(Int, keys(file["timeseries/t"]))
-iter = iterations[end] # plot final iteration
-
-## Temperature
-T = file["timeseries/T/$iter"][1, 1, :]
-
-## Velocity variances
-w² = file["timeseries/ww/$iter"][1, 1, :]
-e  = file["timeseries/e/$iter"][1, 1, :]
-
-## Terms in the TKE budget
-buoyancy_flux =   file["timeseries/tke_buoyancy_flux/$iter"][1, 1, :]
-  dissipation = - file["timeseries/tke_dissipation/$iter"][1, 1, :]
-
- pressure_flux = - file["timeseries/tke_pressure_flux/$iter"][1, 1, :]
-advective_flux = - file["timeseries/tke_advective_flux/$iter"][1, 1, :]
-
-transport = zeros(grid.Nz)
-transport = (pressure_flux[2:end] .+ advective_flux[2:end]
-             .- pressure_flux[1:end-1] .- advective_flux[1:end-1]) / Δz
-
-## For mixing length calculation
-wT = file["timeseries/wT/$iter"][1, 1, 2:end-1]
-
-close(file)
-
-## Post-process the data to determine the mixing length
-
-## Mixing length, computed at cell interfaces and omitting boundaries
-Tz = @. (T[2:end] - T[1:end-1]) / Δz
-bz = @. α * g * Tz
-eᶠ = @. (e[1:end-1] + e[2:end]) / 2
-
-## Mixing length model: wT ∝ - ℓᵀ √e ∂z T ⟹  ℓᵀ = wT / (√e ∂z T)
-ℓ_measured = @. - wT / (√(eᶠ) * Tz)
-ℓ_estimated = @. min(-zF[2:end-1], sqrt(eᶠ / max(0, bz)))
-
-# Plot data
-
-temperature = plot(T, zC, size = plot_size,
-                     linewidth = linewidth,
-                        xlabel = "Temperature (ᵒC)",
-                        ylabel = "z (m)",
-                          ylim = ylim,
-                         label = nothing)
-
-variances = plot(e, zC, size = plot_size,
-                     linewidth = linewidth,
-                        xlabel = "Velocity variances (m² s⁻²)",
-                        ylabel = "z (m)",
-                          ylim = ylim,
-                         label = "(u² + v² + w²) / 2")
-
-plot!(variances, 1/2 .* w², zF, linewidth = linewidth,
-                                    label = "w² / 2")
-
-budget = plot([buoyancy_flux dissipation transport], zC, size = plot_size,
-              linewidth = linewidth,
-                 xlabel = "TKE budget terms",
-                 ylabel = "z (m)",
-                   ylim = ylim,
-                  label = ["buoyancy flux" "dissipation" "kinetic energy transport"])
-
-mixing_length = plot([ℓ_measured ℓ_estimated], zF[2:end-1], size = plot_size,
-                                                       linewidth = linewidth,
-                                                          xlabel = "Mixing length (m)",
-                                                          ylabel = "z (m)",
-                                                            xlim = (-5, 20),
-                                                            ylim = ylim,
-                                                           label = ["measured" "estimated"])
-
-plot(temperature, variances, budget, mixing_length, layout=(1, 4))
-