Update to R_rho + correctly set coefficients in linear eos

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)

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
 

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)

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

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₁₁₀