Skip to content

Commit

Permalink
Add solar heating cases
Browse files Browse the repository at this point in the history
  • Loading branch information
@glwagner
glwagner committed May 22, 2024
1 parent 2f4825b commit be00b8c
Show file tree
Hide file tree
Showing 5 changed files with 103 additions and 80 deletions.
4 changes: 2 additions & 2 deletions Manifest.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
26 changes: 14 additions & 12 deletions idealized/run_lesbrary_simulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
]

#=
Expand All @@ -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))
Expand Down
36 changes: 6 additions & 30 deletions src/IdealizedExperiments/IdealizedExperiments.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -44,28 +46,14 @@ 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),
:med_wind_med_cooling => Dict{Symbol, Any}(:momentum_flux => -4.5e-4, :buoyancy_flux => 1.6e-7, :f => 1e-4),
: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
Expand All @@ -74,28 +62,14 @@ 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),
:med_wind_med_cooling => Dict{Symbol, Any}(:momentum_flux => -3.4e-4, :buoyancy_flux => 8.0e-8, :f => 1e-4),
: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
Expand All @@ -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
Expand All @@ -120,3 +95,4 @@ for (name, set) in seventy_two_hour_suite_parameters
end

end # module

111 changes: 78 additions & 33 deletions src/IdealizedExperiments/three_layer_constant_fluxes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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ᶜᶜᶜ

Expand All @@ -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...)
Expand All @@ -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,
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand All @@ -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..."
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand All @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/NearSurfaceTurbulenceModels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit be00b8c

Please sign in to comment.