diff --git a/examples/four_stair_staircase.jl b/examples/four_stair_staircase.jl index 8ebd664..05db535 100644 --- a/examples/four_stair_staircase.jl +++ b/examples/four_stair_staircase.jl @@ -11,7 +11,7 @@ model = DNSModel(architecture, domain_extent, resolution, diffusivities) ## Initial conditions number_of_steps = 4 depth_of_steps = [-0.2, -0.4, -0.6, -0.8] -salinity = [34.57, 34.6, 34.63, 34.66, 34.69] +salinity = [34.57, 34.60, 34.63, 34.66, 34.69] temperature = [-1.5, -1.0, -0.5, 0.0, 0.5] step_ics = StepInitialConditions(model, number_of_steps, depth_of_steps, salinity, temperature) @@ -21,7 +21,7 @@ sdns = StaircaseDNS(model, step_ics) set_staircase_initial_conditions!(sdns) Δt = 1e-2 -stop_time = 2* 60 # seconds +stop_time = 2 * 60 # seconds save_schedule = 60 # seconds output_path = joinpath(@__DIR__, "output") simulation = SDNS_simulation_setup(sdns, Δt, stop_time, save_schedule, save_computed_output!; output_path) diff --git a/src/StaircaseShenanigans.jl b/src/StaircaseShenanigans.jl index 1b6a5dc..ce1c511 100644 --- a/src/StaircaseShenanigans.jl +++ b/src/StaircaseShenanigans.jl @@ -5,7 +5,7 @@ using Oceananigans: seawater_density using SeawaterPolynomials.TEOS10 using SeawaterPolynomials.SecondOrderSeawaterPolynomials using SeawaterPolynomials: BoussinesqEquationOfState -using SeawaterPolynomials: thermal_expansion, haline_contraction +using SeawaterPolynomials: thermal_expansion, haline_contraction, ρ using GibbsSeaWater: gsw_alpha, gsw_beta using NCDatasets, JLD2 diff --git a/src/alt_lineareos.jl b/src/alt_lineareos.jl index 3d2fab8..48c85c2 100644 --- a/src/alt_lineareos.jl +++ b/src/alt_lineareos.jl @@ -3,20 +3,20 @@ Return a linear `RoquetSeawaterPolynomial` with custom values for α and β which are computed using GibbsSeaWater.jl from `Θ` and `S` (`S` is practical salinity). """ -function CustomLinearRoquetSeawaterPolynomial(Θ, S, FT) +function CustomLinearRoquetSeawaterPolynomial(Θ, S, reference_density, FT) - α = gsw_alpha(S, Θ, 0) - β = gsw_beta(S, Θ, 0) + α = reference_density * gsw_alpha(S, Θ, 0) + β = reference_density * gsw_beta(S, Θ, 0) return SecondOrderSeawaterPolynomial{FT}(R₀₁₀ = α, R₁₀₀ = β) end "Extend `RoquetSeawaterPolynomial` for `:CustomLinear` coefficent set" -SecondOrderSeawaterPolynomials.RoquetSeawaterPolynomial(Θ, S, FT=Float64, coefficient_set=:CustomLinear) = - eval(Symbol(coefficient_set, :RoquetSeawaterPolynomial))(Θ, S, FT) +SecondOrderSeawaterPolynomials.RoquetSeawaterPolynomial(Θ, S, reference_density, FT=Float64, coefficient_set=:CustomLinear) = + eval(Symbol(coefficient_set, :RoquetSeawaterPolynomial))(Θ, S, reference_density, FT) "Constructor for `CustomLienarEquationOfState`. The `Θ` and `S` inputs are the values that α and β (the coefficients in the equation of state) will be calculated at." CustomLinearEquationOfState(Θ, S; reference_density = 1024.6) = - BoussinesqEquationOfState(RoquetSeawaterPolynomial(Θ, S), reference_density) + BoussinesqEquationOfState(RoquetSeawaterPolynomial(Θ, S, reference_density), reference_density) diff --git a/src/staircase_initial_conditions.jl b/src/staircase_initial_conditions.jl index 3379b1f..1ec9ec0 100644 --- a/src/staircase_initial_conditions.jl +++ b/src/staircase_initial_conditions.jl @@ -3,6 +3,7 @@ abstract type AbstractStaircaseInitialConditions <: AbstractInitialConditions en Base.iterate(sics::AbstractStaircaseInitialConditions, state = 1) = state > length(fieldnames(sics)) ? nothing : (getfield(sdns, state), state + 1) +"Container for initial conditions that have well mixed layers seperated by sharp step interfaces." struct StepInitialConditions{T} <: AbstractStaircaseInitialConditions "Number of step changes in the initial state" number_of_steps :: Int @@ -18,59 +19,65 @@ end function StepInitialConditions(model, number_of_steps, depth_of_steps, salinity, temperature) eos = model.buoyancy.model.equation_of_state - ΔT = diff(temperature) - ΔS = diff(salinity) - α, β = compute_α_and_β(salinity, temperature, depth_of_steps, eos) - - R_ρ = similar(depth_of_steps) - @. R_ρ = (β * ΔS) / (α * ΔT) + R_ρ = compute_R_ρ(salinity, temperature, depth_of_steps, eos) return StepInitialConditions(number_of_steps, depth_of_steps, salinity, temperature, R_ρ) end """ - function compute_α_and_β(salinity, temperature, eos) -Compute thermal expansion and haline contraction coefficients at the interface of the steps. -The coefficients are computed as α̂ = 0.5 * (α(Sᵢ, 0.5 * (Θᵢ+Θⱼ), 0) + α(Sⱼ, 0.5 * (Θᵢ+Θⱼ), 0)) -where j = i + 1. Still need to double check if there is a more accurate way to do this as -I think this is a slight simplication of the method I am meant to be using. + function compute_R_ρ(salinity, temperature, depth_of_steps, eos) +Compute the density ratio, ``R_{\rho}``, at a diffusive interface with a non-linear equation of state +as defined in [McDougall (1981)](https://www.sciencedirect.com/science/article/abs/pii/0079661181900021) +equation (1) on page 92. """ -function compute_α_and_β(salinity, temperature, depth_of_steps, eos) +function compute_R_ρ(salinity, temperature, depth_of_steps, eos) + + R_ρ = similar(depth_of_steps) - S1 = @view salinity[1:end-1] - S2 = @view salinity[2:end] - Sstepmean = (S1 .+ S2) / 2 + if isequal(is_linear_eos(eos.seawater_polynomial), "nonlineareos") - T1 = @view temperature[1:end-1] - T2 = @view temperature[2:end] - Tstepmean = (T1 .+ T2) / 2 + S_u = S_g = @view salinity[1:end-1] + S_l = S_f = @view salinity[2:end] + Θ_u = Θ_f = @view temperature[1:end-1] + Θ_l = Θ_g = @view temperature[2:end] - eos_vec = fill(eos, length(S1)) + eos_vec = fill(eos, length(S_u)) - α = 0.5 * (thermal_expansion.(Tstepmean, S1, depth_of_steps, eos_vec) + - thermal_expansion.(Tstepmean, S2, depth_of_steps, eos_vec)) + ρ_u = ρ.(Θ_u, S_u, depth_of_steps, eos_vec) + ρ_l = ρ.(Θ_l, S_l, depth_of_steps, eos_vec) + ρ_f = ρ.(Θ_f, S_f, depth_of_steps, eos_vec) + ρ_g = ρ.(Θ_g, S_g, depth_of_steps, eos_vec) - β = 0.5 * (haline_contraction.(T1, Sstepmean, depth_of_steps, eos_vec) + - haline_contraction.(T2, Sstepmean, depth_of_steps, eos_vec)) + R_ρ = @. (0.5 * (ρ_f - ρ_u) + 0.5 * (ρ_l - ρ_g)) / + (0.5 * (ρ_f - ρ_l) + 0.5 * (ρ_u - ρ_g)) - return α, β -end -# Would prefer to implement these when I get a chance. -# "Get the `geopotential_height` from the `grid` of a `model`." -# geopotential_height(grid) = KernelFunctionOperation{Center, Center, Face}(Zᶜᶜᶠ, grid) + else + + S_u = @view salinity[1:end-1] + S_l = @view salinity[2:end] + S_m = (S_u .+ S_l) / 2 -# function compute_gp_height(grid, depth_of_steps) + Θ_u = @view temperature[1:end-1] + Θ_l = @view temperature[2:end] + Θ_m = (Θ_u .+ Θ_l) / 2 -# gh = Field(geopotential_height(grid)) -# compute!(gh) + eos_vec = fill(eos, length(depth_of_steps)) -# step_gh_idx = [findfirst(interior(gh, 1, 1, :) .≥ d) for d ∈ depth_of_steps] -# step_gh_height = interior(gh, 1, 1, step_gh_idx) + α = thermal_expansion.(Θ_m, S_m, depth_of_steps, eos_vec) + β = haline_contraction.(Θ_m, S_m, depth_of_steps, eos_vec) -# return Array(step_gh_height) -# end + ΔT = diff(temperature) + ΔS = diff(salinity) + + R_ρ = @. (β * ΔS) / (α * ΔT) + + end + + return R_ρ +end +"Container for initial conditions that have well mixed layers seperated by smoothed step interfaces." struct SmoothStepInitialConditions{T} <: AbstractStaircaseInitialConditions "Number of step changes in the initial state" number_of_steps :: Int diff --git a/test/runtests.jl b/test/runtests.jl index 7d845c7..a68e5a2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,8 +31,9 @@ using Test for _ in 1:5 Θ = rand(range(-1.5, 20.5, length = 20)) S = rand(range(34.5, 34.7, length = 20)) - true_coefficients = gsw_alpha(S, Θ, 0), gsw_beta(S, Θ, 0) - custom_linear = CustomLinearEquationOfState(Θ, S) + ρ₀ = 1024.6 + true_coefficients = ρ₀ * gsw_alpha(S, Θ, 0), ρ₀ * gsw_beta(S, Θ, 0) + custom_linear = CustomLinearEquationOfState(Θ, S, reference_density = ρ₀) sp = custom_linear.seawater_polynomial linear_coefficients = sp.R₀₁₀, sp.R₁₀₀ other_coefficients = sp.R₁₀₁, sp.R₂₀₀, sp.R₀₂₀, sp.R₁₁₀