diff --git a/Manifest.toml b/Manifest.toml index 041a784..5ec1f7d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -735,9 +735,9 @@ version = "1.2.0" [[deps.Oceananigans]] deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "4672af7242405313743af45168bfce3d87b84b2c" +git-tree-sha1 = "3a5d2b1d9a237731f886a657e72e4a52a7f6fb2b" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.90.11" +version = "0.90.14" [deps.Oceananigans.extensions] OceananigansEnzymeExt = "Enzyme" diff --git a/idealized/run_lesbrary_simulation.jl b/idealized/run_lesbrary_simulation.jl index 9652645..fbb3e24 100644 --- a/idealized/run_lesbrary_simulation.jl +++ b/idealized/run_lesbrary_simulation.jl @@ -32,22 +32,24 @@ using LESbrary.IdealizedExperiments: seventy_two_hour_suite_parameters # which slices of the simulation is saved) from the default 2 minutes to 1 minute # (to make pretty movies). -architecture = GPU() -#size = (64, 64, 64) -size = (128, 128, 128) +architecture = CPU() +#size = (32, 32, 32) +size = (64, 64, 64) +#size = (128, 128, 128) #size = (256, 256, 256) #size = (256, 256, 384) # case = :strong_wind snapshot_time_interval = 10minute -data_directory = "/home/greg/Projects/LESbrary.jl/data" +data_directory = "." #/home/greg/Projects/LESbrary.jl/data" cases = [ - :strong_wind, - :free_convection, - :weak_wind_strong_cooling, - :med_wind_med_cooling, - :strong_wind_weak_cooling, - :strong_wind_no_rotation + :strong_wind_and_sunny, + #:strong_wind, + #:free_convection, + #:weak_wind_strong_cooling, + #:med_wind_med_cooling, + #:strong_wind_weak_cooling, + #:strong_wind_no_rotation ] #= @@ -56,10 +58,10 @@ suites = [twelve_hour_suite_parameters, forty_eight_hour_suite_parameters] =# -#suite = six_hour_suite_parameters +suite = six_hour_suite_parameters #suite = twelve_hour_suite_parameters #suite = twenty_four_hour_suite_parameters -suite = forty_eight_hour_suite_parameters +#suite = forty_eight_hour_suite_parameters #suite = seventy_two_hour_suite_parameters @inline κˢ(x, y, z, t) = ifelse(z > -5, 1e-2, zero(z)) diff --git a/src/IdealizedExperiments/IdealizedExperiments.jl b/src/IdealizedExperiments/IdealizedExperiments.jl index a07855c..3a8fca2 100644 --- a/src/IdealizedExperiments/IdealizedExperiments.jl +++ b/src/IdealizedExperiments/IdealizedExperiments.jl @@ -21,6 +21,7 @@ six_hour_suite_parameters = Dict{Symbol, Any}( :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -1.2e-3, :buoyancy_flux => 4.0e-7, :f => 1e-4), :strong_wind => Dict{Symbol, Any}(:momentum_flux => -1.4e-3, :buoyancy_flux => 0.0, :f => 1e-4), :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -1.1e-3, :buoyancy_flux => 0.0, :f => 0.0), + :strong_wind_and_sunny => Dict{Symbol, Any}(:momentum_flux => -1.5e-3, :penetrating_buoyancy_flux => -5e-7, :f => 0), ) for (name, set) in six_hour_suite_parameters @@ -36,6 +37,7 @@ twelve_hour_suite_parameters = Dict{Symbol, Any}( :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -8.0e-4, :buoyancy_flux => 2.0e-7, :f => 1e-4), :strong_wind => Dict{Symbol, Any}(:momentum_flux => -9.0e-4, :buoyancy_flux => 0.0, :f => 1e-4), :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -6.0e-4, :buoyancy_flux => 0.0, :f => 0.0), + :strong_wind_and_sunny => Dict{Symbol, Any}(:momentum_flux => -7.0e-4, :penetrating_buoyancy_flux => -5e-7, :f => 0), ) for (name, set) in twelve_hour_suite_parameters @@ -44,21 +46,6 @@ for (name, set) in twelve_hour_suite_parameters set[:stokes_drift] = true end -eighteen_hour_suite_parameters = Dict{Symbol, Any}( - :free_convection => Dict{Symbol, Any}(:momentum_flux => 0.0, :buoyancy_flux => 3.6e-7, :f => 1e-4), - :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -3.5e-4, :buoyancy_flux => 3.0e-7, :f => 1e-4), - :med_wind_med_cooling => Dict{Symbol, Any}(:momentum_flux => -5.0e-4, :buoyancy_flux => 2.4e-7, :f => 1e-4), - :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -7.0e-4, :buoyancy_flux => 1.5e-7, :f => 1e-4), - :strong_wind => Dict{Symbol, Any}(:momentum_flux => -8.0e-4, :buoyancy_flux => 0.0, :f => 1e-4), - :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -4.2e-4, :buoyancy_flux => 0.0, :f => 0.0), -) - -for (name, set) in eighteen_hour_suite_parameters - set[:name] = string(name) - set[:stop_time] = 18hours - set[:stokes_drift] = true -end - twenty_four_hour_suite_parameters = Dict{Symbol, Any}( :free_convection => Dict{Symbol, Any}(:momentum_flux => 0.0, :buoyancy_flux => 2.4e-7, :f => 1e-4), :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -3.0e-4, :buoyancy_flux => 2.0e-7, :f => 1e-4), @@ -66,6 +53,7 @@ twenty_four_hour_suite_parameters = Dict{Symbol, Any}( :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -5.9e-4, :buoyancy_flux => 1.0e-7, :f => 1e-4), :strong_wind => Dict{Symbol, Any}(:momentum_flux => -6.8e-4, :buoyancy_flux => 0.0, :f => 1e-4), :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -3.0e-4, :buoyancy_flux => 0.0, :f => 0.0), + :strong_wind_and_sunny => Dict{Symbol, Any}(:momentum_flux => -4.0e-4, :penetrating_buoyancy_flux => -5e-7, :f => 0), ) for (name, set) in twenty_four_hour_suite_parameters @@ -74,21 +62,6 @@ for (name, set) in twenty_four_hour_suite_parameters set[:stokes_drift] = true end -thirty_six_hour_suite_parameters = Dict{Symbol, Any}( - :free_convection => Dict{Symbol, Any}(:momentum_flux => 0.0, :buoyancy_flux => 1.8e-7, :f => 1e-4), - :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -2.5e-4, :buoyancy_flux => 1.5e-7, :f => 1e-4), - :med_wind_med_cooling => Dict{Symbol, Any}(:momentum_flux => -3.9e-4, :buoyancy_flux => 1.2e-7, :f => 1e-4), - :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -4.8e-4, :buoyancy_flux => 7.5e-8, :f => 1e-4), - :strong_wind => Dict{Symbol, Any}(:momentum_flux => -5.6e-4, :buoyancy_flux => 0.0, :f => 1e-4), - :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -2.3e-4, :buoyancy_flux => 0.0, :f => 0.0), -) - -for (name, set) in thirty_six_hour_suite_parameters - set[:name] = string(name) - set[:stop_time] = 36hours - set[:stokes_drift] = true -end - forty_eight_hour_suite_parameters = Dict{Symbol, Any}( :free_convection => Dict{Symbol, Any}(:momentum_flux => 0.0, :buoyancy_flux => 1.2e-7, :f => 1e-4), :weak_wind_strong_cooling => Dict{Symbol, Any}(:momentum_flux => -2.0e-4, :buoyancy_flux => 1.0e-7, :f => 1e-4), @@ -96,6 +69,7 @@ forty_eight_hour_suite_parameters = Dict{Symbol, Any}( :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -3.8e-4, :buoyancy_flux => 5.0e-8, :f => 1e-4), :strong_wind => Dict{Symbol, Any}(:momentum_flux => -4.5e-4, :buoyancy_flux => 0.0, :f => 1e-4), :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -1.6e-4, :buoyancy_flux => 0.0, :f => 0.0), + :strong_wind_and_sunny => Dict{Symbol, Any}(:momentum_flux => -2.0e-4, :penetrating_buoyancy_flux => -5e-7, :f => 0), ) for (name, set) in forty_eight_hour_suite_parameters @@ -111,6 +85,7 @@ seventy_two_hour_suite_parameters = Dict{Symbol, Any}( :strong_wind_weak_cooling => Dict{Symbol, Any}(:momentum_flux => -3.4e-4, :buoyancy_flux => 3.8e-8, :f => 1e-4), :strong_wind => Dict{Symbol, Any}(:momentum_flux => -4.1e-4, :buoyancy_flux => 0.0, :f => 1e-4), :strong_wind_no_rotation => Dict{Symbol, Any}(:momentum_flux => -1.1e-4, :buoyancy_flux => 0.0, :f => 0.0), + :strong_wind_and_sunny => Dict{Symbol, Any}(:momentum_flux => -1.5e-4, :penetrating_buoyancy_flux => -5e-7, :f => 0), ) for (name, set) in seventy_two_hour_suite_parameters @@ -120,3 +95,4 @@ for (name, set) in seventy_two_hour_suite_parameters end end # module + diff --git a/src/IdealizedExperiments/three_layer_constant_fluxes.jl b/src/IdealizedExperiments/three_layer_constant_fluxes.jl index 11e4731..3468164 100644 --- a/src/IdealizedExperiments/three_layer_constant_fluxes.jl +++ b/src/IdealizedExperiments/three_layer_constant_fluxes.jl @@ -5,6 +5,7 @@ using OrderedCollections using Oceanostics using Oceananigans using Oceananigans.Units +using Oceananigans.Grids: zspacing using Oceananigans.Utils: WallTimeInterval using Oceananigans.Operators: Δzᶜᶜᶜ @@ -23,6 +24,36 @@ Logging.global_logger(OceananigansLogger()) @inline passive_tracer_forcing(x, y, z, t, p) = p.μ⁺ * exp(-(z - p.z₀)^2 / (2 * p.λ^2)) - p.μ⁻ +const c = Center() +const f = Face() + +struct TwoExponentialPenetratingFlux{FT} <: Function + surface_flux :: FT + first_fraction :: FT + first_decay_scale :: FT + second_decay_scale :: FT +end + +function TwoExponentialPenetratingFlux(surface_flux, + first_fraction = 0.6, + first_decay_scale = 1.0, + second_decay_scale = 20.0) + + return TwoExponentialPenetratingFlux(surface_flux, + first_fraction, + first_decay_scale, + second_decay_scale) +end + +function minus_penetrating_flux_divergence(grid; I₀, ϵ₁=0.6, λ₁=1.0, λ₂=20.0) + I_field = Field{Nothing, Nothing, Face}(grid) + I(z) = I₀ * (ϵ₁ * exp(z / λ₁) + (1 - ϵ₁) * exp(z / λ₂)) + set!(I_field, I) + dIdz = Field(-1 * ∂z(I_field)) + compute!(dIdz) + return dIdz +end + """ three_layer_constant_fluxes_simulation(; kw...) @@ -34,14 +65,15 @@ function three_layer_constant_fluxes_simulation(; name = "", size = (32, 32, 32), passive_tracers = false, - background_diffusivity = nothing, + subgrid_scale_closure = true, extent = (512meters, 512meters, 256meters), architecture = CPU(), stop_time = 0.1hours, initial_Δt = 1.0, - f = 1e-4, - buoyancy_flux = 1e-8, - momentum_flux = -1e-4, + f = 0.0, + buoyancy_flux = 0.0, + penetrating_buoyancy_flux = nothing, + momentum_flux = 0.0, thermocline_type = "linear", surface_layer_depth = 48.0, thermocline_width = 24.0, @@ -66,13 +98,14 @@ function three_layer_constant_fluxes_simulation(; Nx, Ny, Nz = size Lx, Ly, Lz = extent slice_depth = 8.0 - Qᵇ = buoyancy_flux - Qᵘ = momentum_flux + Jᵇ = buoyancy_flux + Iᵇ = penetrating_buoyancy_flux + τˣ = momentum_flux stop_hours = stop_time / hour ## Determine filepath prefix prefix = @sprintf("three_layer_constant_fluxes_%s_hr%d_Qu%.1e_Qb%.1e_f%.1e_Nh%d_Nz%d_", - thermocline_type, stop_hours, abs(Qᵘ), Qᵇ, f, Nx, Nz) + thermocline_type, stop_hours, abs(τˣ), Jᵇ, f, Nx, Nz) data_directory = joinpath(data_directory, prefix * name) # save data in /data/prefix @@ -116,15 +149,15 @@ function three_layer_constant_fluxes_simulation(; α = buoyancy.equation_of_state.thermal_expansion g = buoyancy.gravitational_acceleration - Qᶿ = Qᵇ / (α * g) + Jᶿ = Jᵇ / (α * g) dθdz_surface_layer = N²_surface_layer / (α * g) dθdz_thermocline = N²_thermocline / (α * g) dθdz_deep = N²_deep / (α * g) - θ_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᶿ), + θ_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᶿ), bottom = GradientBoundaryCondition(dθdz_deep)) - u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) + u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τˣ)) # Tracer forcing @@ -149,10 +182,7 @@ function three_layer_constant_fluxes_simulation(; μ⁺ = 1 / 6hour μ₀ = √(2π) * λ / grid.Lz * μ⁺ / 2 μ∞ = √(2π) * λ / grid.Lz * μ⁺ - - c₀_forcing = Forcing(passive_tracer_forcing, parameters=(z₀= 0.0, λ=λ, μ⁺=μ⁺, μ⁻=μ₀)) - c₁_forcing = Forcing(passive_tracer_forcing, parameters=(z₀=-48.0, λ=λ, μ⁺=μ⁺, μ⁻=μ∞)) - c₂_forcing = Forcing(passive_tracer_forcing, parameters=(z₀=-96.0, λ=λ, μ⁺=μ⁺, μ⁻=μ∞)) + c_forcing = Forcing(passive_tracer_forcing, parameters=(z₀=-48.0, λ=λ, μ⁺=μ⁺, μ⁻=μ∞)) # Sponge layer for u, v, w, and T gaussian_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10) @@ -161,6 +191,15 @@ function three_layer_constant_fluxes_simulation(; T_sponge = Relaxation(rate = 4/hour, target = LinearTarget{:z}(intercept = θ_deep - z_deep*dθdz_deep, gradient = dθdz_deep), mask = gaussian_mask) + + if penetrating_buoyancy_flux isa Number + Iᶿ = penetrating_buoyancy_flux / (α * g) + dIdz = minus_penetrating_flux_divergence(grid; I₀=Iᶿ) + T_penetrating_flux = Forcing(dIdz) + T_forcing = (T_sponge, T_penetrating_flux) + elseif isnothing(penetrating_buoyancy_flux) + T_forcing = T_sponge + end if stokes_drift && momentum_flux != 0.0 @info "Whipping up the Stokes drift..." @@ -181,26 +220,30 @@ function three_layer_constant_fluxes_simulation(; @info "Framing the model..." - #tracers = passive_tracers ? (:T, :c₀, :c₁, :c₂) : :T tracers = passive_tracers ? (:T, :c) : :T - if isnothing(background_diffusivity) - closure = nothing + if subgrid_scale_closure + Nz = Base.size(grid, 3) + Δz = zspacing(1, 1, Nz, grid, c, c, c) + C = SurfaceEnhancedModelConstant(Δz) + closure = AnisotropicMinimumDissipation(; C) + advection = CenteredSecondOrder() else - closure = ScalarDiffusivity(κ=background_diffusivity) + closure = nothing + advection = WENO(order=9) end - model = NonhydrostaticModel(; grid, buoyancy, tracers, stokes_drift, closure, + model = NonhydrostaticModel(; grid, buoyancy, tracers, stokes_drift, closure, advection, timestepper = :RungeKutta3, - advection = WENO(order=9), coriolis = FPlane(; f), boundary_conditions = (T=θ_bcs, u=u_bcs), - forcing = (u=u_sponge, v=v_sponge, w=w_sponge, T=T_sponge, - c=c₁_forcing)) - #c₀=c₀_forcing, c₁=c₁_forcing, c₂=c₂_forcing)) + forcing = (u=u_sponge, v=v_sponge, w=w_sponge, T=T_forcing, c=c_forcing)) # # Set Initial condition + @info "Build model:" + @info "$model" + @info "Setting initial conditions..." ## Noise with 8 m decay scale @@ -281,12 +324,11 @@ function three_layer_constant_fluxes_simulation(; cff_scratch = Field{Center, Face, Face}(model.grid) if statistics == "first_order" - primitive_statistics = first_through_second_order(model, b=b, p=p, w_scratch=ccf_scratch, c_scratch=ccc_scratch) - elseif statistics == "second_order" primitive_statistics = first_order_statistics(model, b=b, p=p, w_scratch=ccf_scratch, c_scratch=ccc_scratch) + elseif statistics == "second_order" + primitive_statistics = first_through_second_order(model, b=b, p=p, w_scratch=ccf_scratch, c_scratch=ccc_scratch) end - U = Field(primitive_statistics[:u]) V = Field(primitive_statistics[:v]) B = Field(primitive_statistics[:b]) @@ -297,28 +339,31 @@ function three_layer_constant_fluxes_simulation(; additional_statistics = Dict(:e => Average(TurbulentKineticEnergy(model, U=U, V=V), dims=(1, 2)), :Ri => ∂z(B) / (∂z(U)^2 + ∂z(V)^2)) - # We're changing to WENO(order=9) so there are no subfilter fluxes... + # If we change to WENO(order=9) so there are no subfilter fluxes... #subfilter_flux_statistics = merge(subfilter_momentum_fluxes(model), subfilter_tracer_fluxes(model)) #statistics_to_output = merge(primitive_statistics, subfilter_flux_statistics, additional_statistics) - statistics_to_output = merge(primitive_statistics, additional_statistics) + @info "Garnishing output writers..." + @info " - with fields: $(keys(fields_to_output))" + @info " - with statistics: $(keys(statistics_to_output))" global_attributes = OrderedDict() global_attributes[:LESbrary_jl_commit_SHA1] = execute(`git rev-parse HEAD`).stdout |> strip global_attributes[:name] = name global_attributes[:thermocline_type] = thermocline_type - global_attributes[:buoyancy_flux] = Qᵇ - global_attributes[:momentum_flux] = Qᵘ - global_attributes[:temperature_flux] = Qᶿ + global_attributes[:buoyancy_flux] = Jᵇ + global_attributes[:penetrating_buoyancy_flux] = penetrating_buoyancy_flux + global_attributes[:momentum_flux] = τˣ + global_attributes[:temperature_flux] = Jᶿ global_attributes[:coriolis_parameter] = f global_attributes[:thermal_expansion_coefficient] = α global_attributes[:gravitational_acceleration] = g - global_attributes[:boundary_condition_θ_top] = Qᶿ + global_attributes[:boundary_condition_θ_top] = Jᶿ global_attributes[:boundary_condition_θ_bottom] = dθdz_deep - global_attributes[:boundary_condition_u_top] = Qᵘ + global_attributes[:boundary_condition_u_top] = τˣ global_attributes[:boundary_condition_u_bottom] = 0.0 global_attributes[:surface_layer_depth] = surface_layer_depth global_attributes[:thermocline_width] = thermocline_width diff --git a/src/NearSurfaceTurbulenceModels.jl b/src/NearSurfaceTurbulenceModels.jl index d8e47fe..994d721 100644 --- a/src/NearSurfaceTurbulenceModels.jl +++ b/src/NearSurfaceTurbulenceModels.jl @@ -21,7 +21,7 @@ end """ SurfaceEnhancedModelConstant(Δz; FT=Float64, C₀=1/12, z₀=-Δz/2, - enhancement=2, decay_scale=8Δz) + enhancement=2, decay_scale=4Δz) Returns a callable object representing a spatially-variable model constant for an LES eddy diffusivity model with the surface-enhanced form @@ -30,7 +30,7 @@ for an LES eddy diffusivity model with the surface-enhanced form """ function SurfaceEnhancedModelConstant(Δz; FT=Float64, C₀=1/12, z₀=-Δz/2, - enhancement=2, decay_scale=8Δz) + enhancement=4, decay_scale=4Δz) return SurfaceEnhancedModelConstant{FT}(C₀, Δz, z₀, enhancement, decay_scale) end @@ -57,7 +57,7 @@ function SurfaceEnhancedModelConstant(filename::String) close(file) - Cᴬᴹᴰ = SurfaceEnhancedModelConstant(Δz; C₀=1/12, z₀=-Δz/2, enhancement=2, decay_scale=8Δz) + Cᴬᴹᴰ = SurfaceEnhancedModelConstant(Δz; C₀=1/12, z₀=-Δz/2, enhancement=4, decay_scale=4Δz) return Cᴬᴹᴰ end