Merge pull request #57 from jbisits/jib-betterbfiledsparameters

More flexibility in setting background fields

examples/single_interface_periodic.jl

@@ -11,7 +11,7 @@ model_setup = (;architecture, diffusivities, domain_extent, resolution, eos)
 depth_of_interface = -0.5
 salinity = [34.58, 34.70]
 temperature = [-1.5, 0.5]
-interface_ics = PeriodoicSingleInterfaceICs(eos, depth_of_interface, salinity, temperature, tanh_background)
+interface_ics = PeriodoicSingleInterfaceICs(eos, depth_of_interface, salinity, temperature, BackgroundTanh())
 tracer_noise = TracerNoise(1e-6, 1e-6)
 
 ## setup model

src/StaircaseShenanigans.jl

@@ -5,6 +5,7 @@ using Oceananigans: seawater_density
 using Oceananigans.BoundaryConditions: update_boundary_condition!
 using Oceananigans.Operators: ℑzᵃᵃᶠ
 using Oceananigans.Fields: condition_operand
+import Oceananigans.BackgroundField
 using SeawaterPolynomials.TEOS10
 using SeawaterPolynomials.SecondOrderSeawaterPolynomials
 using SeawaterPolynomials: BoussinesqEquationOfState
@@ -29,7 +30,7 @@ export STStaircaseInitialConditions, StaircaseICs, SmoothSTStaircaseInitialCondi
        PeriodicSTSingleInterfaceInitialConditions, PeriodoicSingleInterfaceICs,
        set_staircase_initial_conditions!
 
-export tanh_background, linear_background
+export BackgroundTanh, BackgroundLinear, tanh_background, linear_background
 
 export AbstractNoise, VelocityNoise, TracerNoise
 

src/staircase_initial_conditions.jl

@@ -148,7 +148,14 @@ struct SmoothSTStaircaseInitialConditions{T} <: AbstractStaircaseInitialConditio
       smoothing_funciton :: Function
 end
 function Base.show(io::IO, sics::AbstractStaircaseInitialConditions)
-    if sics isa STSingleInterfaceInitialConditions
+    if sics isa PeriodicSTSingleInterfaceInitialConditions
+        println(io, "STSingleInterfaceInitialConditions")
+        println(io, "┣━ depth_of_interface: $(sics.depth_of_interface)")
+        println(io, "┣━━━━ salinity_values: $(sics.salinity_values)")
+        println(io, "┣━ temperature_values: $(sics.temperature_values)")
+        println(io, "┣━━━━━━━━━━━━━━━━ R_ρ: $(round.(sics.R_ρ; digits = 2))")
+        println(io, "┗━━━ background_state: $(typeof(sics.background_state))")
+    elseif sics isa STSingleInterfaceInitialConditions
         println(io, "STSingleInterfaceInitialConditions")
         println(io, "┣━ depth_of_interface: $(sics.depth_of_interface)")
         println(io, "┣━━━━ salinity_values: $(sics.salinity_values)")

src/staircase_model.jl

@@ -53,8 +53,7 @@ function StaircaseDNS(model_setup::NamedTuple, initial_conditions::PeriodoicSing
 
     architecture, diffusivities, domain_extent, resolution, eos = model_setup
     z_topology = Periodic
-    τ = diffusivities.κ.S / diffusivities.κ.T
-    background_fields = S_and_T_background_fields(initial_conditions, domain_extent.Lz, τ)
+    background_fields = S_and_T_background_fields(initial_conditions, domain_extent.Lz)
     model = DNSModel(architecture, domain_extent, resolution, diffusivities, eos;
                      z_topology, background_fields)
 

src/staircase_restoring.jl

@@ -1,3 +1,140 @@
+abstract type BackgroundFunction end
+
+Oceananigans.BackgroundField(bf::BackgroundFunction) =
+    BackgroundField(bf.func, parameters = bf.parameters)
+"""
+    mutable struct BackgroundTanh{F, T, P}
+Container for a tanh background field.
+"""
+mutable struct BackgroundTanh{F, T} <: BackgroundFunction
+    "tanh function"
+          func :: F
+    "Scale the steepness of the `tanh` change"
+         scale :: T
+    "Parameters for the tanh background field"
+    parameters :: Any
+end
+BackgroundTanh() = BackgroundTanh(tanh_background, 100, NamedTuple())
+BackgroundTanh(scale) = BackgroundTanh(tanh_background, scale, NamedTuple())
+function Base.show(io, bt::BackgroundTanh)
+    println(io, "BackgroundTanh")
+    println(io, "┣━━━ function: $(bt.func)")
+    println(io, "┣━━━━━━ scale: $(bt.scale)")
+      print(io, "┗━ parameters: $(bt.parameters)")
+end
+"""
+    mutable struct BackgroundLinear{F, P}
+Container for a linear background field.
+"""
+mutable struct BackgroundLinear{F} <: BackgroundFunction
+    "Linear function"
+          func :: F
+    "Parameters for the tanh background field"
+    parameters :: Any
+end
+BackgroundLinear() = BackgroundLinear(linear_background, NamedTuple())
+function Base.show(io, bl::BackgroundLinear)
+    println(io, "BackgroundLinear")
+    println(io, "┣━━━ function: $(bl.func)")
+      print(io, "┗━ parameters: $(bl.parameters)")
+end
+"""
+    function S_and_T_background_fields(initial_conditions)
+Set background fields for the `S` and `T` tracer fields where the domain is triply periodic.
+"""
+function S_and_T_background_fields(ics::PeriodicSTSingleInterfaceInitialConditions, Lz)
+
+    get_parameters!(ics, :salinity_values, Lz)
+    S_background = BackgroundField(ics.background_state)
+    get_parameters!(ics, :temperature_values, Lz)
+    T_background = BackgroundField(ics.background_state)
+
+    return (S = S_background, T = T_background)
+end
+"Sets a background state that is hyperbolic tangent. The scaling `τ` is set by
+the diffusivity ratio κₛ / κₜ ."
+  tanh_background(x, y, z, t, p) = p.Cₗ - 0.5 * p.ΔC * (1  + tanh(p.D * (z - p.z_interface) / p.Lz))
+linear_background(x, y, z, t, p) = p.Cᵤ - p.ΔC * z / p.Lz
+
+function get_parameters!(ics::PeriodicSTSingleInterfaceInitialConditions, tracer::Symbol, Lz)
+
+    z_interface = ics.depth_of_interface
+    C = getproperty(ics, tracer)
+    ΔC = diff(C)[1]
+    Cᵤ, Cₗ = C
+    update_parameters!(ics.background_state, ΔC, Cᵤ, Cₗ, abs(Lz), z_interface)
+
+    return nothing
+end
+
+function update_parameters!(backgound_state::BackgroundTanh, ΔC, Cᵤ, Cₗ, Lz, z_interface)
+
+    backgound_state.parameters = (; ΔC, Cᵤ, Cₗ, Lz, z_interface, D = backgound_state.scale)
+    return nothing
+end
+function update_parameters!(backgound_state::BackgroundLinear, ΔC, Cᵤ, Lz, z_interface)
+
+    backgound_state.parameters = (; ΔC, Cᵤ, Lz)
+    return nothing
+end
+
+# The following functions are to be used as `BoundaryConditions` so that tracers can
+# re-enter the domain with the initial gradient added effectively allowing the gradient to
+# be maintained. Another option is to add background tracer fields and only evolve the anomaly
+
+"""
+    T_reentry(i, j, grid, clock, model_fields, ΔT)
+Discrete boundary condition to set the temperature on a vertical boundary to tracer value at
+the **top** of the domain and adds ΔT. That is a reentrant T condition to be used on the
+**bottom** of the domain.
+"""
+@inline T_reentry(i, j, grid, clock, model_fields, ΔT) =
+    @inbounds ℑzᵃᵃᶠ(i, j, grid.Nz+1, model_fields.T) + ΔT
+"""
+    S_reentry(i, j, grid, clock, model_fields, ΔS)
+As for [T_reentry](@ref) but using salinity tracer instead.
+"""
+@inline S_reentry(i, j, grid, clock, model_fields, ΔS) =
+    @inbounds ℑzᵃᵃᶠ(i, j, grid.Nz+1, model_fields.S) + ΔS
+
+"Access velocity field at the _bottom_ of the domain for use as `BoundaryCondition`."
+@inline _w_bottom(i, j, grid, clock, model_fields) = @inbounds model_fields.w[i, j, 1]
+"Access velocity field at the _top_ of the domain for use as `BoundaryCondition`."
+@inline _w_top(i, j, grid, clock, model_fields) = @inbounds model_fields.w[i, j, grid.Nz+1]
+
+"Interpolated velocity between top and bottom (vertical) faces of domain."
+@inline w_top_bottom_interpolate(i, j, grid, clock, model_fields) =
+    @inbounds 0.5 * (model_fields.w[i, j, 1] + model_fields.w[i, j, grid.Nz+1])
+
+""""
+    function reentrant_boundary_conditions(ics::SingleInterfaceICs)
+Setup boundary conditions to maintain a gradient (ΔC) between the two layers of a
+`SingleInterface` model. The boundary conditions for the tracers are re-entrant from top to
+bottom with ΔC added to maintain a ΔC difference between the top and bottom layer. The
+vertical velocitiy field also has modified vertical boundary conditions where the velocity
+on the bottom face is set to the velocity on the top face and vice versa.
+"""
+function reentrant_boundary_conditions(ics::SingleInterfaceICs)
+
+    bcs = if ics.maintain_interface
+            ΔT = diff(ics.temperature_values)[1]
+            ΔS = diff(ics.salinity_values)[1]
+            T_bottom_reentry = ValueBoundaryCondition(T_reentry, discrete_form=true, parameters = ΔT)
+            S_bottom_reentry = ValueBoundaryCondition(S_reentry, discrete_form=true, parameters = ΔS)
+            T_bcs = FieldBoundaryConditions(bottom = T_bottom_reentry)
+            S_bcs = FieldBoundaryConditions(bottom = S_bottom_reentry)
+
+            w_top = OpenBoundaryCondition(_w_bottom, discrete_form=true)
+            w_bottom = OpenBoundaryCondition(_w_top, discrete_form=true)
+            w_bcs = FieldBoundaryConditions(top = w_top, bottom = w_bottom)
+
+            (T=T_bcs, S=S_bcs, w=w_bcs)
+        else
+            NamedTuple()
+        end
+
+    return bcs
+end
 """
     struct OuterStairMask{T}
 Container for location of first and last stair for `mask` in `Relaxation`. This need not be
@@ -88,85 +225,3 @@ struct ExponentialTarget{T}
 end
 
 @inline (p::ExponentialTarget)(x, y, z, t) = p.A * exp(-p.λ * z)
-
-# The following functions are to be used as `BoundaryConditions` so that tracers can
-# re-enter the domain with the initial gradient added effectively allowing the gradient to
-# be maintained. Another option is to add background tracer fields and only evolve the anomaly
-
-"""
-    T_reentry(i, j, grid, clock, model_fields, ΔT)
-Discrete boundary condition to set the temperature on a vertical boundary to tracer value at
-the **top** of the domain and adds ΔT. That is a reentrant T condition to be used on the
-**bottom** of the domain.
-"""
-@inline T_reentry(i, j, grid, clock, model_fields, ΔT) =
-    @inbounds ℑzᵃᵃᶠ(i, j, grid.Nz+1, model_fields.T) + ΔT
-"""
-    S_reentry(i, j, grid, clock, model_fields, ΔS)
-As for [T_reentry](@ref) but using salinity tracer instead.
-"""
-@inline S_reentry(i, j, grid, clock, model_fields, ΔS) =
-    @inbounds ℑzᵃᵃᶠ(i, j, grid.Nz+1, model_fields.S) + ΔS
-
-"Access velocity field at the _bottom_ of the domain for use as `BoundaryCondition`."
-@inline _w_bottom(i, j, grid, clock, model_fields) = @inbounds model_fields.w[i, j, 1]
-"Access velocity field at the _top_ of the domain for use as `BoundaryCondition`."
-@inline _w_top(i, j, grid, clock, model_fields) = @inbounds model_fields.w[i, j, grid.Nz+1]
-
-"Interpolated velocity between top and bottom (vertical) faces of domain."
-@inline w_top_bottom_interpolate(i, j, grid, clock, model_fields) =
-    @inbounds 0.5 * (model_fields.w[i, j, 1] + model_fields.w[i, j, grid.Nz+1])
-
-""""
-    function reentrant_boundary_conditions(ics::SingleInterfaceICs)
-Setup boundary conditions to maintain a gradient (ΔC) between the two layers of a
-`SingleInterface` model. The boundary conditions for the tracers are re-entrant from top to
-bottom with ΔC added to maintain a ΔC difference between the top and bottom layer. The
-vertical velocitiy field also has modified vertical boundary conditions where the velocity
-on the bottom face is set to the velocity on the top face and vice versa.
-"""
-function reentrant_boundary_conditions(ics::SingleInterfaceICs)
-
-    bcs = if ics.maintain_interface
-            ΔT = diff(ics.temperature_values)[1]
-            ΔS = diff(ics.salinity_values)[1]
-            T_bottom_reentry = ValueBoundaryCondition(T_reentry, discrete_form=true, parameters = ΔT)
-            S_bottom_reentry = ValueBoundaryCondition(S_reentry, discrete_form=true, parameters = ΔS)
-            T_bcs = FieldBoundaryConditions(bottom = T_bottom_reentry)
-            S_bcs = FieldBoundaryConditions(bottom = S_bottom_reentry)
-
-            w_top = OpenBoundaryCondition(_w_bottom, discrete_form=true)
-            w_bottom = OpenBoundaryCondition(_w_top, discrete_form=true)
-            w_bcs = FieldBoundaryConditions(top = w_top, bottom = w_bottom)
-
-            (T=T_bcs, S=S_bcs, w=w_bcs)
-        else
-            NamedTuple()
-        end
-
-    return bcs
-end
-
-"""
-    function S_and_T_background_fields(initial_conditions)
-Set background fields for the `S` and `T` tracer fields where the domain is triply periodic.
-"""
-function S_and_T_background_fields(ics::PeriodicSTSingleInterfaceInitialConditions, Lz, τ)
-
-    z_interface = ics.depth_of_interface
-    ΔT = diff(ics.temperature_values)[1]
-    Tᵤ, Tₗ = ics.temperature_values
-    T_parameters = (Cᵤ = Tᵤ, Cₗ = Tₗ, ΔC = ΔT, Lz = abs(Lz), z_interface, τ)
-
-    ΔS = diff(ics.salinity_values)[1]
-    Sᵤ, Sₗ = ics.salinity_values
-    S_parameters = (Cᵤ = Sᵤ, Cₗ = Sₗ, ΔC = ΔS, Lz = abs(Lz), z_interface, τ)
-    S_background = BackgroundField(ics.background_state, parameters=S_parameters)
-    T_background = BackgroundField(ics.background_state, parameters=T_parameters)
-
-    return (S = S_background, T = T_background)
-end
-"Sets a background state that is hyperbolic tangent. The scaling `τ` is set by
-the diffusivity ratio κₛ / κₜ ."
-tanh_background(x, y, z, t, p) =  p.Cₗ - 0.5 * p.ΔC * (1  + tanh(round(1 / p.τ) * (z - p.z_interface) / p.Lz))
-linear_background(x, y, z, t, p) = p.Cᵤ - p.ΔC * z / p.Lz