Skip to content

Commit

Permalink
Merge pull request #78 from jbisits/jib-initialconditionsfortripleper…
Browse files Browse the repository at this point in the history
…iodic

Initial conditions for triple periodic domain that is not just noise
  • Loading branch information
jbisits authored Nov 24, 2024
2 parents 9567634 + e7558a8 commit b8ab0f5
Show file tree
Hide file tree
Showing 13 changed files with 403 additions and 337 deletions.
2 changes: 1 addition & 1 deletion examples/four_stair_staircase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ set_staircase_initial_conditions!(sdns)
Δt = 1e-1
stop_time = 10 * 60 # seconds
save_schedule = 10 # seconds
output_path = joinpath(@__DIR__, "output")
output_path = joinpath(@__DIR__, "output_staircase")
simulation = SDNS_simulation_setup(sdns, stop_time, save_computed_output!;
output_path, Δt, save_schedule,)

Expand Down
9 changes: 3 additions & 6 deletions examples/single_interface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 = SingleInterfaceICs(eos, depth_of_interface, salinity, temperature, maintain_interface = true)
interface_ics = SingleInterfaceICs(eos, depth_of_interface, salinity, temperature, interface_smoothing = Tanh)
velocity_noise = VelocityNoise(0.0, 0.0, 1e-7)

## setup model
Expand All @@ -21,14 +21,11 @@ sdns = StaircaseDNS(model_setup, interface_ics, velocity_noise)
Δt = 1e-1
stop_time = 210 * 60 # seconds
save_schedule = 10 # seconds
output_path = joinpath(@__DIR__, "output")
simulation = SDNS_simulation_setup(sdns, stop_time, save_computed_output!,
save_vertical_velocities!;
output_path = joinpath(@__DIR__, "output_single_interface")
simulation = SDNS_simulation_setup(sdns, stop_time, save_computed_output!, save_vertical_velocities!;
Δt, save_schedule,
output_path, max_Δt = 5)
## Run
run!(simulation)

## Compute density ratio
compute_R_ρ!(simulation.output_writers[:computed_output].filepath,
simulation.output_writers[:tracers].filepath, eos)
10 changes: 4 additions & 6 deletions examples/single_interface_periodic_tanh_background.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,7 @@ model_setup = (;architecture, diffusivities, domain_extent, resolution, eos)
depth_of_interface = -0.5
salinity = [34.56, 34.70]
temperature = [-1.5, 0.5]
interface_ics = PeriodoicSingleInterfaceICs(eos, depth_of_interface, salinity, temperature,
BackgroundTanh())
interface_ics = PeriodoicSingleInterfaceICs(eos, depth_of_interface, salinity, background_state = BackgroundTanh())
tracer_noise = TracerNoise(1e-6, 1e-6)

## setup model
Expand All @@ -22,14 +21,13 @@ sdns = StaircaseDNS(model_setup, interface_ics, tracer_noise)
Δt = 1e-1
stop_time = 4 * 60 * 60 # seconds
save_schedule = 30 # seconds
output_path = joinpath(@__DIR__, "output_tanh_background")
simulation = SDNS_simulation_setup(sdns, stop_time, save_computed_output!,
save_vertical_velocities!;
output_path = joinpath(@__DIR__, "output_periodic")
simulation = SDNS_simulation_setup(sdns, stop_time, save_computed_output!, save_vertical_velocities!;
Δt, save_schedule,
output_path, max_Δt = 5)
## Run
run!(simulation)

## Compute density ratio
# Compute density ratio
compute_R_ρ!(simulation.output_writers[:computed_output].filepath,
simulation.output_writers[:tracers].filepath, eos)
12 changes: 9 additions & 3 deletions src/StaircaseShenanigans.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,16 +22,19 @@ import SeawaterPolynomials.SecondOrderSeawaterPolynomials: RoquetSeawaterPolynom
abstract type AbstractStaircaseModel end
abstract type AbstractInitialConditions end
abstract type AbstractNoise end
abstract type BackgroundFunction end
abstract type AbstractBackgroundFunction end
abstract type AbstractInterfaceSmoothing end

export StaircaseDNS, PeriodicStaircaseDNS, DNSModel, SDNS_simulation_setup
export StaircaseDNS, DNSModel, SDNS_simulation_setup

export STStaircaseInitialConditions, StaircaseICs, SmoothSTStaircaseInitialConditions,
STSingleInterfaceInitialConditions, SingleInterfaceICs,
PeriodicSTSingleInterfaceInitialConditions, PeriodoicSingleInterfaceICs,
set_staircase_initial_conditions!

export BackgroundTanh, BackgroundLinear, tanh_background, linear_background
export TanhInterfaceSmoothing, Tanh, NoSmoothing

export BackgroundTanh, BackgroundLinear, tanh_background, linear_background, NoBackground

export AbstractNoise, VelocityNoise, TracerNoise

Expand All @@ -46,11 +49,14 @@ export compute_R_ρ!
export animate_tracers, animate_density, visualise_initial_conditions, visualise_initial_density,
animate_tracers_anomaly, animate_density_anomaly

include("interface_smoothing.jl")
include("staircase_initial_conditions.jl")
include("single_interfaces_initial_conditions.jl")
include("staircase_noise.jl")
include("staircase_diagnostics.jl")
include("staircase_model.jl")
include("set_staircase_initial_conditions.jl")
include("staircase_background.jl")
include("staircase_restoring.jl")
include("alt_lineareos.jl")
include("makie_plotting_functions.jl")
Expand Down
24 changes: 24 additions & 0 deletions src/interface_smoothing.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
"""
struct TanhInterfaceSmoothing
Container to set a hyperbolic tangent over a single interface as an initial condition.
"""
struct TanhInterfaceSmoothing{T} <: AbstractInterfaceSmoothing
"Tracer in the deeper of the two layers seperated by an interface"
Cₗ :: T
"Change in tracer content over an interface"
ΔC :: T
"Steepness of the tanh change"
D :: T
"Location of the interface"
z_interface :: T
"Total z range of the two layers seperated by a given interface"
z_range :: T
end
const Tanh = TanhInterfaceSmoothing{T} where {T}
@inline (p::TanhInterfaceSmoothing)(x, y, z) = p.Cₗ - 0.5 * p.ΔC * (1 + tanh(p.D * (z - p.z_interface) / p.z_range))

Base.summary(p::Type{<:TanhInterfaceSmoothing}) = "tanh smoothing"

struct NoInterfaceSmoothing <: AbstractInterfaceSmoothing end
const NoSmoothing = NoInterfaceSmoothing() = nothing
Base.summary(ns::NoInterfaceSmoothing) = "no interface smoothing"
25 changes: 20 additions & 5 deletions src/set_staircase_initial_conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,9 @@ function set_staircase_initial_conditions!(model, ics::SmoothSTStaircaseInitialC

return nothing
end
function set_staircase_initial_conditions!(model, ics::SingleInterfaceICs)
set_staircase_initial_conditions!(model, ics::Union{SingleInterfaceICs, PeriodoicSingleInterfaceICs}) =
set_staircase_initial_conditions!(model, ics, ics.interface_smoothing)
function set_staircase_initial_conditions!(model, ics, interface_smoothing::Nothing)

depth_of_interface = ics.depth_of_interface
z = znodes(model.grid, Center())
Expand All @@ -51,9 +53,20 @@ function set_staircase_initial_conditions!(model, ics::SingleInterfaceICs)

return nothing
end
"Here the `BackgroundField` behaves as the `initial_condition` and noise is added to the
tracer fields to create an instability."
set_staircase_initial_conditions!(model, ics::PeriodoicSingleInterfaceICs) = nothing
function set_staircase_initial_conditions!(model, ics, interface_smoothing::Type{<:Tanh})

depth_of_interface = ics.depth_of_interface
Lz = model.grid.Lz
S = ics.salinity_values
T = ics.temperature_values

S₀(x, y, z) = Tanh(S[2], diff(S)[1], 100.0, depth_of_interface, abs(Lz))(x, y, z)
T₀(x, y, z) = Tanh(T[2], diff(T)[1], 100.0, depth_of_interface, abs(Lz))(x, y, z)

set!(model, S = S₀, T = T₀)

return nothing
end

"Fallback --- don't set any noise."
set_noise!(model, noise::Nothing) = nothing
Expand All @@ -80,8 +93,10 @@ function set_noise!(model, noise::TracerNoise)
S, T = model.tracers
S_noise = noise.S_magnitude * randn(size(S))
T_noise = noise.T_magnitude * randn(size(T))
S₀ = S + S_noise
T₀ = T + T_noise

set!(model, S = S_noise, T = T_noise)
set!(model, S = S₀, T = T₀)

end
"""
Expand Down
112 changes: 112 additions & 0 deletions src/single_interfaces_initial_conditions.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
abstract type AbstractSingleInterfaceInitialConditions <: AbstractInitialConditions end

Base.iterate(sics::AbstractSingleInterfaceInitialConditions, state = 1) =
state > length(fieldnames(sics)) ? nothing : (getfield(sics, state), state + 1)
"""
struct STSingleInterfaceInitialConditions
Initial conditions for a single interface (i.e. two layers with uniform `S` and `T` seperated
by a step change). The property `maintain_interface` is a `Boolean` which if set to `true` will
set [reentrant_boundary_conditions](@ref) so that the interface will be maintained (by not
letting the system run down).
"""
struct STSingleInterfaceInitialConditions{T, A, IS} <: AbstractSingleInterfaceInitialConditions
"The depth of the interface"
depth_of_interface :: T
"Salinity values in each layer"
salinity_values :: A
"Temperature values in each layer"
temperature_values :: A
"Initial smoothing function over interface"
interface_smoothing :: IS
"Initial R_ρ at the interface"
R_ρ :: T
end
function STSingleInterfaceInitialConditions(model, depth_of_interface, salinity, temperature;
interface_smoothing = nothing)

eos = model.buoyancy.model.equation_of_state

R_ρ = compute_R_ρ(salinity, temperature, depth_of_interface, eos)

return STSingleInterfaceInitialConditions(depth_of_interface, salinity, temperature,
interface_smoothing, R_ρ)

end
function STSingleInterfaceInitialConditions(eos::BoussinesqEquationOfState, depth_of_interface,
salinity, temperature; interface_smoothing = nothing)

R_ρ = compute_R_ρ(salinity, temperature, depth_of_interface, eos)

return STSingleInterfaceInitialConditions(depth_of_interface, salinity, temperature,
interface_smoothing, R_ρ)

end
const SingleInterfaceICs = STSingleInterfaceInitialConditions # alias
Base.summary(ics::SingleInterfaceICs) = "Single S-T interface at z = $(ics.depth_of_interface) on bounded z domain"

"""
struct PeriodicSTSingleInterfaceInitialConditions
Sets a `BackgroundField` according to `background_State` and uses a triply periodic domain
to evolve salinity and temperature anomalies about the background state.
"""
struct PeriodicSTSingleInterfaceInitialConditions{T, A, IS, BF} <: AbstractSingleInterfaceInitialConditions
"The depth of the interface"
depth_of_interface :: T
"Salinity values in each layer"
salinity_values :: A
"Temperature values in each layer"
temperature_values :: A
"Initial smoothing function over interface"
interface_smoothing :: IS
"Initial R_ρ at the interface"
R_ρ :: T
"`BackgroundField` about which anomaly is advected. Should be an `AbstractBackgroundFunction`."
background_state :: BF
end
function PeriodicSTSingleInterfaceInitialConditions(eos::BoussinesqEquationOfState, depth_of_interface,
salinity, temperature;
interface_smoothing = NoSmoothing(),
background_state = NoBackground())

R_ρ = compute_R_ρ(salinity, temperature, depth_of_interface, eos)

return PeriodicSTSingleInterfaceInitialConditions(depth_of_interface, salinity, temperature,
interface_smoothing, R_ρ, background_state)

end
const PeriodoicSingleInterfaceICs = PeriodicSTSingleInterfaceInitialConditions # alias
Base.summary(ics::PeriodoicSingleInterfaceICs) = "Single S-T interface at z = $(ics.depth_of_interface) on triply periodic domain"

function compute_R_ρ(salinity, temperature, depth_of_interfaces::Number, eos)

S_u = S_g = salinity[1]
S_l = S_f = salinity[2]
T_u = T_f = temperature[1]
T_l = T_g = temperature[2]

ρ_u = ρ(T_u, S_u, depth_of_interfaces, eos)
ρ_l = ρ(T_l, S_l, depth_of_interfaces, eos)
ρ_f = ρ(T_f, S_f, depth_of_interfaces, eos)
ρ_g = ρ(T_g, S_g, depth_of_interfaces, eos)

return (0.5 * (ρ_f - ρ_u) + 0.5 * (ρ_l - ρ_g)) / (0.5 * (ρ_f - ρ_l) + 0.5 * (ρ_u - ρ_g))
end

function Base.show(io::IO, sics::AbstractSingleInterfaceInitialConditions)
if sics isa PeriodicSTSingleInterfaceInitialConditions
println(io, "PeriodicSTSingleInterfaceInitialConditions")
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, "┣━ interface_smoothing: $(summary(sics.interface_smoothing))")
println(io, "┣━━━━ background_state: $(summary(sics.background_state))")
print(io, "┗━━━━━━━━━━━━━━━━━ R_ρ: $(round.(sics.R_ρ; digits = 2))")
elseif sics isa STSingleInterfaceInitialConditions
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, "┣━ interface_smoothing: $(summary(sics.interface_smoothing))")
print(io, "┗━━━━━━━━━━━━━━━━━ R_ρ: $(round.(sics.R_ρ; digits = 2))")
end
end
Loading

0 comments on commit b8ab0f5

Please sign in to comment.