From 5fcbcfa20de354a499309860f406eb69e6139263 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Mon, 20 Oct 2025 18:08:26 -0700 Subject: [PATCH 01/16] first commit for array averagedspecifiedtimes --- src/OutputWriters/windowed_time_average.jl | 58 +++++++++++++++++++--- src/Utils/prettytime.jl | 1 + src/Utils/schedules.jl | 2 + src/Utils/times_and_datetimes.jl | 8 +++ 4 files changed, 62 insertions(+), 7 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index c0c98d2692..73d8617a24 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -91,14 +91,12 @@ function (sch::AveragedTimeInterval)(model) return scheduled end initialize_schedule!(sch::AveragedTimeInterval, clock) = nothing -outside_window(sch::AveragedTimeInterval, clock) = clock.time <= next_actuation_time(sch) - sch.window +outside_window(sch::AveragedTimeInterval, clock) = clock.time <= next_actuation_time(sch) - sch.window end_of_window(sch::AveragedTimeInterval, clock) = clock.time >= next_actuation_time(sch) TimeInterval(sch::AveragedTimeInterval) = TimeInterval(sch.interval) Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window=sch.window, stride=sch.stride) - - """ mutable struct AveragedSpecifiedTimes <: AbstractSchedule @@ -144,6 +142,11 @@ function end_of_window(schedule::AveragedSpecifiedTimes, clock) return clock.time >= next_time end +TimeInterval(sch::AveragedSpecifiedTimes) = TimeInterval(sch.specified_times.times) +Base.copy(sch::AveragedSpecifiedTimes) = AveragedSpecifiedTimes(copy(sch.specified_times); window=sch.window, stride=sch.stride) + +next_actuation_time(sch::AveragedSpecifiedTimes) = Oceananigans.Utils.next_actuation_time(sch.specified_times) + ##### ##### WindowedTimeAverage ##### @@ -261,6 +264,40 @@ function advance_time_average!(wta::WindowedTimeAverage, model) return nothing end +function advance_time_average!(wta::SpecifiedWindowedTimeAverage, model) + + unscheduled = model.clock.iteration == 0 || outside_window(wta.schedule, model.clock) + if !(unscheduled) + if !(wta.schedule.collecting) + # Zero out result to begin new accumulation window + wta.result .= 0 + + # Begin collecting window-averaged increments + wta.schedule.collecting = true + + wta.window_start_time = next_actuation_time(wta.schedule) - wta.schedule.window + wta.previous_collection_time = wta.window_start_time + wta.window_start_iteration = model.clock.iteration - 1 + # @info "t $(prettytime(model.clock.time)), next actuation time: $(prettytime(next_actuation_time(wta.schedule))), window $(prettytime(wta.schedule.window))" + end + + if end_of_window(wta.schedule, model.clock) + accumulate_result!(wta, model) + # Save averaging start time and the initial data collection time + wta.schedule.collecting = false + wta.schedule.specified_times.previous_actuation += 1 + + elseif mod(model.clock.iteration - wta.window_start_iteration, stride(wta)) == 0 + accumulate_result!(wta, model) + else + # Off stride, so do nothing. + end + + end + return nothing +end + + # So it can be used as a Diagnostic run_diagnostic!(wta::WindowedTimeAverage, model) = advance_time_average!(wta, model) @@ -271,8 +308,14 @@ Base.summary(schedule::AveragedTimeInterval) = string("AveragedTimeInterval(", "stride=", schedule.stride, ", ", "interval=", prettytime(schedule.interval), ")") +Base.summary(schedule::AveragedSpecifiedTimes) = string("AveragedSpecifiedTimes(", + "window=", prettytime(schedule.window), ", ", + "stride=", schedule.stride, ", ", + "times=", schedule.specified_times, ")") + show_averaging_schedule(schedule) = "" show_averaging_schedule(schedule::AveragedTimeInterval) = string(" averaged on ", summary(schedule)) +show_averaging_schedule(schedule::AveragedSpecifiedTimes) = string(" averaged on ", summary(schedule)) output_averaging_schedule(output::WindowedTimeAverage) = output.schedule @@ -282,6 +325,8 @@ output_averaging_schedule(output::WindowedTimeAverage) = output.schedule time_average_outputs(schedule, outputs, model) = schedule, outputs # fallback +const AveragedTimeSchedule = Union{AveragedTimeInterval, AveragedSpecifiedTimes} + """ time_average_outputs(schedule::AveragedTimeInterval, outputs, model, field_slicer) @@ -290,17 +335,16 @@ Wrap each `output` in a `WindowedTimeAverage` on the time-averaged `schedule` an Returns the `TimeInterval` associated with `schedule` and a `NamedTuple` or `Dict` of the wrapped outputs. """ -function time_average_outputs(schedule::AveragedTimeInterval, outputs::Dict, model) +function time_average_outputs(schedule::AveragedTimeSchedule, outputs::Dict, model) averaged_outputs = Dict(name => WindowedTimeAverage(output, model; schedule=copy(schedule)) for (name, output) in outputs) return TimeInterval(schedule), averaged_outputs end -function time_average_outputs(schedule::AveragedTimeInterval, outputs::NamedTuple, model) +function time_average_outputs(schedule::AveragedTimeSchedule, outputs::NamedTuple, model) averaged_outputs = NamedTuple(name => WindowedTimeAverage(outputs[name], model; schedule=copy(schedule)) for name in keys(outputs)) return TimeInterval(schedule), averaged_outputs -end - +end \ No newline at end of file diff --git a/src/Utils/prettytime.jl b/src/Utils/prettytime.jl index ac71dbd676..f2e3e20f85 100644 --- a/src/Utils/prettytime.jl +++ b/src/Utils/prettytime.jl @@ -64,3 +64,4 @@ function prettytimeunits(t, longform=true) end prettytime(dt::AbstractTime) = "$dt" +prettytime(t::Array) = prettytime.(t) diff --git a/src/Utils/schedules.jl b/src/Utils/schedules.jl index d7c72763c6..88bab3269c 100644 --- a/src/Utils/schedules.jl +++ b/src/Utils/schedules.jl @@ -234,6 +234,8 @@ function specified_times_str(st) return string(str, "]") end +Base.copy(st::SpecifiedTimes) = SpecifiedTimes(copy(st.times), st.previous_actuation) + ##### ##### ConsecutiveIterations ##### diff --git a/src/Utils/times_and_datetimes.jl b/src/Utils/times_and_datetimes.jl index 7521c64115..a2c0df07f8 100644 --- a/src/Utils/times_and_datetimes.jl +++ b/src/Utils/times_and_datetimes.jl @@ -29,11 +29,19 @@ end @inline add_time_interval(base::AbstractTime, interval::Number, count=1) = base + seconds_to_nanosecond(interval * count) @inline add_time_interval(base::AbstractTime, interval::Period, count=1) = base + count * interval +@inline add_time_interval(base::Number, interval::Array{<:Number}, count=1) = interval[count] + function period_type(interval::Number) FT = Oceananigans.defaults.FloatType return FT end +function period_type(interval::Array{<:Number}) + FT = Oceananigans.defaults.FloatType + return Array{FT, 1} +end + period_type(interval::Dates.Period) = typeof(interval) time_type(interval::Number) = typeof(interval) time_type(interval::Dates.Period) = Dates.DateTime +time_type(interval::Array{<:Number}) = eltype(interval) From 5457bc51d2557cbb921086ff01372a766cb4f593 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 16:39:37 -0700 Subject: [PATCH 02/16] allow window to be arrays --- src/OutputWriters/windowed_time_average.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 73d8617a24..6f76b78612 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -102,18 +102,23 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window A schedule for averaging over windows that precede SpecifiedTimes. """ -mutable struct AveragedSpecifiedTimes <: AbstractSchedule +mutable struct AveragedSpecifiedTimes{W} <: AbstractSchedule specified_times :: SpecifiedTimes - window :: Float64 + window :: W stride :: Int collecting :: Bool end +const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{Vector{Float64}} + AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) = AveragedSpecifiedTimes(specified_times, window, stride, false) AveragedSpecifiedTimes(times; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); kw...) +get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1] +get_next_window(schedule::AveragedSpecifiedTimes) = schedule.window + function (schedule::AveragedSpecifiedTimes)(model) time = model.clock.time @@ -121,7 +126,7 @@ function (schedule::AveragedSpecifiedTimes)(model) next > length(schedule.specified_times.times) && return false next_time = schedule.specified_times.times[next] - window = schedule.window + window = get_next_window(schedule) schedule.collecting || time >= next_time - window end @@ -132,7 +137,8 @@ function outside_window(schedule::AveragedSpecifiedTimes, clock) next = schedule.specified_times.previous_actuation + 1 next > length(schedule.specified_times.times) && return true next_time = schedule.specified_times.times[next] - return clock.time < next_time - schedule.window + window = get_next_window(schedule) + return clock.time < next_time - window end function end_of_window(schedule::AveragedSpecifiedTimes, clock) @@ -171,7 +177,7 @@ stride(wta::SpecifiedWindowedTimeAverage) = wta.schedule.stride WindowedTimeAverage(operand, model=nothing; schedule) Returns an object for computing running averages of `operand` over `schedule.window` and -recurring on `schedule.interval`, where `schedule` is an `AveragedTimeInterval`. +recurring on `schedule.interval`, where `schedule` is an `AveragedTimeInterval` or `AveragedSpecifiedTimes`. During the collection period, averages are computed every `schedule.stride` iteration. `operand` may be a `Oceananigans.Field` or a function that returns an array or scalar. @@ -275,7 +281,7 @@ function advance_time_average!(wta::SpecifiedWindowedTimeAverage, model) # Begin collecting window-averaged increments wta.schedule.collecting = true - wta.window_start_time = next_actuation_time(wta.schedule) - wta.schedule.window + wta.window_start_time = next_actuation_time(wta.schedule) - get_next_window(wta.schedule) wta.previous_collection_time = wta.window_start_time wta.window_start_iteration = model.clock.iteration - 1 # @info "t $(prettytime(model.clock.time)), next actuation time: $(prettytime(next_actuation_time(wta.schedule))), window $(prettytime(wta.schedule.window))" From 96872aa73e57d83e42f488c934ff5e304b3718ad Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 17:52:25 -0700 Subject: [PATCH 03/16] Enhance AveragedSpecifiedTimes to support vector windows and add validation for overlapping windows --- src/OutputWriters/windowed_time_average.jl | 52 +++++++++++++++++++++- 1 file changed, 50 insertions(+), 2 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 6f76b78612..c174e820d8 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -102,7 +102,7 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window A schedule for averaging over windows that precede SpecifiedTimes. """ -mutable struct AveragedSpecifiedTimes{W} <: AbstractSchedule +mutable struct AveragedSpecifiedTimes{W <: Union{Float64, Vector{Float64}}} <: AbstractSchedule specified_times :: SpecifiedTimes window :: W stride :: Int @@ -114,7 +114,55 @@ const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{Vector{Float6 AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) = AveragedSpecifiedTimes(specified_times, window, stride, false) -AveragedSpecifiedTimes(times; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); kw...) +AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, window; kw...) +# AveragedSpecifiedTimes(times; window::Float64; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) + +function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) + length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) + perm = sortperm(times) + sorted_times = times[perm] + sorted_window = window[perm] + time_diff = diff(vcat(0, sorted_times)) + + @info "timediff", time_diff + @info "sortedwindow", sorted_window + + any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) + + return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) +end + +function AveragedSpecifiedTimes(times, window::Float64; kw...) + perm = sortperm(times) + sorted_times = times[perm] + time_diff = diff(vcat(0, sorted_times)) + + any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) + + return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) +end + +# function AveragedSpecifiedTimes(times; window, kw...) +# perm = sortperm(times) +# sorted_times = times[perm] +# time_diff = diff(vcat(0, sorted_times)) + +# if window isa Vector{Float64} +# length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) + +# sorted_window = window[perm] +# @info "timediff", time_diff +# @info "sortedwindow", sorted_window + +# any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) +# elseif window isa Number +# any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) +# else +# throw(ArgumentError("window must be a Float64 or a Vector{Float64}, got $(typeof(window))")) +# end + +# return AveragedSpecifiedTimes(SpecifiedTimes(times); window=window, kw...) +# end get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1] get_next_window(schedule::AveragedSpecifiedTimes) = schedule.window From 61d29ac359fe2b37231e899a01da140c92655550 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 18:02:59 -0700 Subject: [PATCH 04/16] Refactor AveragedSpecifiedTimes to ensure non-overlapping windows and improve time sorting --- src/OutputWriters/windowed_time_average.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index c174e820d8..a557f60420 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -129,12 +129,11 @@ function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) - return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) + return AveragedSpecifiedTimes(SpecifiedTimes(times); window=sorted_window, kw...) end function AveragedSpecifiedTimes(times, window::Float64; kw...) - perm = sortperm(times) - sorted_times = times[perm] + sorted_times = sort(times) time_diff = diff(vcat(0, sorted_times)) any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) From c20c9b09e550a21970d7fc7960b8d9eba9d1c7f4 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 18:03:25 -0700 Subject: [PATCH 05/16] Fix AveragedSpecifiedTimes to use sorted times for averaging windows --- src/OutputWriters/windowed_time_average.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index a557f60420..c7de95978f 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -129,7 +129,7 @@ function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) - return AveragedSpecifiedTimes(SpecifiedTimes(times); window=sorted_window, kw...) + return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) end function AveragedSpecifiedTimes(times, window::Float64; kw...) From 35c41a517021ea1f0a8751d7f9c2611717425dea Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 18:09:27 -0700 Subject: [PATCH 06/16] Remove debug logging from AveragedSpecifiedTimes function --- src/OutputWriters/windowed_time_average.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index c7de95978f..d1f67f1855 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -115,7 +115,6 @@ AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) = AveragedSpecifiedTimes(specified_times, window, stride, false) AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, window; kw...) -# AveragedSpecifiedTimes(times; window::Float64; kw...) = AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) @@ -124,9 +123,6 @@ function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) sorted_window = window[perm] time_diff = diff(vcat(0, sorted_times)) - @info "timediff", time_diff - @info "sortedwindow", sorted_window - any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) From 6d168250410165f3cc61b319cf87415ec7d6a769 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 21 Oct 2025 18:14:40 -0700 Subject: [PATCH 07/16] update commented out method --- src/OutputWriters/windowed_time_average.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index d1f67f1855..9cd1e25ecf 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -150,13 +150,13 @@ end # @info "sortedwindow", sorted_window # any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) +# return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) # elseif window isa Number # any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) +# return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) # else # throw(ArgumentError("window must be a Float64 or a Vector{Float64}, got $(typeof(window))")) # end - -# return AveragedSpecifiedTimes(SpecifiedTimes(times); window=window, kw...) # end get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1] From 6ea30f5eb0580d3d7c931e6f743f4791fafcfdea Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Wed, 22 Oct 2025 16:02:52 -0700 Subject: [PATCH 08/16] Refactor AveragedSpecifiedTimes to support varying window types and improve overlap validation --- src/OutputWriters/windowed_time_average.jl | 48 +++++++++------------- 1 file changed, 20 insertions(+), 28 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 9cd1e25ecf..146972a4b4 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -2,6 +2,7 @@ using Oceananigans.Diagnostics: AbstractDiagnostic using Oceananigans.OutputWriters: fetch_output using Oceananigans.Utils: AbstractSchedule, prettytime using Oceananigans.TimeSteppers: Clock +using Dates: Period import Oceananigans: run_diagnostic! import Oceananigans.Utils: TimeInterval, SpecifiedTimes @@ -102,63 +103,54 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window A schedule for averaging over windows that precede SpecifiedTimes. """ -mutable struct AveragedSpecifiedTimes{W <: Union{Float64, Vector{Float64}}} <: AbstractSchedule +mutable struct AveragedSpecifiedTimes{W} <: AbstractSchedule specified_times :: SpecifiedTimes window :: W stride :: Int collecting :: Bool end -const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{Vector{Float64}} +const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{<:Vector} AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) = AveragedSpecifiedTimes(specified_times, window, stride, false) AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, window; kw...) -function AveragedSpecifiedTimes(times, window::Vector{Float64}; kw...) +function determine_epsilon(eltype) + if eltype <: AbstractFloat + return eps(eltype) + elseif eltype <: Period + return Second(0) + else + return 0 + end +end + +function AveragedSpecifiedTimes(times, window::Vector; kw...) length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) perm = sortperm(times) sorted_times = times[perm] sorted_window = window[perm] time_diff = diff(vcat(0, sorted_times)) - any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) + epsilon = determine_epsilon(eltype(window)) + any(time_diff .- sorted_window .< -epsilon) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) end -function AveragedSpecifiedTimes(times, window::Float64; kw...) +function AveragedSpecifiedTimes(times, window::Union{<:Number, <:Period}; kw...) sorted_times = sort(times) time_diff = diff(vcat(0, sorted_times)) - any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) + epsilon = determine_epsilon(typeof(window)) + + any(time_diff .- window .< -epsilon) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) end -# function AveragedSpecifiedTimes(times; window, kw...) -# perm = sortperm(times) -# sorted_times = times[perm] -# time_diff = diff(vcat(0, sorted_times)) - -# if window isa Vector{Float64} -# length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) - -# sorted_window = window[perm] -# @info "timediff", time_diff -# @info "sortedwindow", sorted_window - -# any(time_diff .- sorted_window .< -eps(eltype(window))) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) -# return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) -# elseif window isa Number -# any(time_diff .- window .< -eps(typeof(window))) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) -# return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) -# else -# throw(ArgumentError("window must be a Float64 or a Vector{Float64}, got $(typeof(window))")) -# end -# end - get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1] get_next_window(schedule::AveragedSpecifiedTimes) = schedule.window From 56e6e55bb1d8a17f68244639a68d4c237310d17a Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Wed, 22 Oct 2025 16:05:22 -0700 Subject: [PATCH 09/16] Refactor AveragedTimeInterval struct to support generic types for interval and first actuation time --- src/OutputWriters/windowed_time_average.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 146972a4b4..34034f496a 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -13,11 +13,11 @@ import Oceananigans.Fields: location, indices, set! Container for parameters that configure and handle time-averaged output. """ -mutable struct AveragedTimeInterval <: AbstractSchedule - interval :: Float64 - window :: Float64 +mutable struct AveragedTimeInterval{I, T} <: AbstractSchedule + interval :: I + window :: I stride :: Int - first_actuation_time :: Float64 + first_actuation_time :: T actuations :: Int collecting :: Bool end From b3f1541e7ce9b62029a2dfcc17277624dfaed24e Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Wed, 22 Oct 2025 16:21:37 -0700 Subject: [PATCH 10/16] make averagedspecifiedtimes concrete --- src/OutputWriters/windowed_time_average.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 34034f496a..814b965fb9 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -103,8 +103,8 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window A schedule for averaging over windows that precede SpecifiedTimes. """ -mutable struct AveragedSpecifiedTimes{W} <: AbstractSchedule - specified_times :: SpecifiedTimes +mutable struct AveragedSpecifiedTimes{S<:SpecifiedTimes, W} <: AbstractSchedule + specified_times :: S window :: W stride :: Int collecting :: Bool From 0efdbc82ccfe6b32d40a52245ac357ab0354c2b8 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Thu, 23 Oct 2025 19:01:21 -0700 Subject: [PATCH 11/16] change dispatch method --- src/OutputWriters/windowed_time_average.jl | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 814b965fb9..a15ab91b1a 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -110,22 +110,16 @@ mutable struct AveragedSpecifiedTimes{S<:SpecifiedTimes, W} <: AbstractSchedule collecting :: Bool end -const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{<:Vector} +const VaryingWindowAveragedSpecifiedTimes = AveragedSpecifiedTimes{<:Any, <:Vector} AveragedSpecifiedTimes(specified_times::SpecifiedTimes; window, stride=1) = AveragedSpecifiedTimes(specified_times, window, stride, false) AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, window; kw...) -function determine_epsilon(eltype) - if eltype <: AbstractFloat - return eps(eltype) - elseif eltype <: Period - return Second(0) - else - return 0 - end -end +determine_epsilon(eltype) = 0 +determine_epsilon(::Type{T}) where T <: AbstractFloat = eps(T) +determine_epsilon(::Period) = Second(0) function AveragedSpecifiedTimes(times, window::Vector; kw...) length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) @@ -152,7 +146,7 @@ function AveragedSpecifiedTimes(times, window::Union{<:Number, <:Period}; kw...) end get_next_window(schedule::VaryingWindowAveragedSpecifiedTimes) = schedule.window[schedule.specified_times.previous_actuation + 1] -get_next_window(schedule::AveragedSpecifiedTimes) = schedule.window +get_next_window(schedule) = schedule.window function (schedule::AveragedSpecifiedTimes)(model) time = model.clock.time From 3ab8e0e54b971be3cea337a3fb6b06b5b94fb5c0 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Thu, 23 Oct 2025 19:04:18 -0700 Subject: [PATCH 12/16] change to fallback --- src/OutputWriters/windowed_time_average.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index a15ab91b1a..90cb6387d4 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -134,7 +134,7 @@ function AveragedSpecifiedTimes(times, window::Vector; kw...) return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) end -function AveragedSpecifiedTimes(times, window::Union{<:Number, <:Period}; kw...) +function AveragedSpecifiedTimes(times, window; kw...) sorted_times = sort(times) time_diff = diff(vcat(0, sorted_times)) From 5e769ba88a4c8a87ed5b230a12bd314abd9ec958 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Mon, 3 Nov 2025 16:55:39 -0500 Subject: [PATCH 13/16] refactor: update AveragedSpecifiedTimes struct to allow more flexible type parameters --- src/OutputWriters/windowed_time_average.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 90cb6387d4..7718c24ddb 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -103,7 +103,7 @@ Base.copy(sch::AveragedTimeInterval) = AveragedTimeInterval(sch.interval, window A schedule for averaging over windows that precede SpecifiedTimes. """ -mutable struct AveragedSpecifiedTimes{S<:SpecifiedTimes, W} <: AbstractSchedule +mutable struct AveragedSpecifiedTimes{S, W} <: AbstractSchedule specified_times :: S window :: W stride :: Int From c2132c37daf93f487585d73ef585b41d3bd067c3 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 4 Nov 2025 09:37:09 -0500 Subject: [PATCH 14/16] md files removed --- docs/src/models/background_fields.md | 107 --- docs/src/models/boundary_conditions.md | 620 ------------------ .../models/buoyancy_and_equation_of_state.md | 251 ------- docs/src/models/coriolis.md | 112 ---- docs/src/models/forcing_functions.md | 366 ----------- docs/src/models/lagrangian_particles.md | 136 ---- docs/src/models/models_overview.md | 208 ------ docs/src/models/turbulence_closures.md | 103 --- docs/src/simulations/callbacks.md | 93 --- docs/src/simulations/checkpointing.md | 61 -- docs/src/simulations/output_writers.md | 246 ------- 11 files changed, 2303 deletions(-) delete mode 100644 docs/src/models/background_fields.md delete mode 100644 docs/src/models/boundary_conditions.md delete mode 100644 docs/src/models/buoyancy_and_equation_of_state.md delete mode 100644 docs/src/models/coriolis.md delete mode 100644 docs/src/models/forcing_functions.md delete mode 100644 docs/src/models/lagrangian_particles.md delete mode 100644 docs/src/models/models_overview.md delete mode 100644 docs/src/models/turbulence_closures.md delete mode 100644 docs/src/simulations/callbacks.md delete mode 100644 docs/src/simulations/checkpointing.md delete mode 100644 docs/src/simulations/output_writers.md diff --git a/docs/src/models/background_fields.md b/docs/src/models/background_fields.md deleted file mode 100644 index ae55119b0b..0000000000 --- a/docs/src/models/background_fields.md +++ /dev/null @@ -1,107 +0,0 @@ -# Background fields - -`BackgroundField`s are velocity and tracer fields around which the resolved -velocity and tracer fields evolve. Only the _advective_ terms associated with -the interaction between background and resolved fields are included. -For example, tracer advection is described by - -```math -\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) \, , -``` - -where ``\boldsymbol{v}`` is the resolved velocity field and ``c`` is the resolved -tracer field corresponding to `model.tracers.c`. - -When a background field ``C`` is provided, the tracer advection term becomes - -```math -\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) - + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} C \right ) \, . -``` - -When both a background field velocity field ``\boldsymbol{U}`` and a background tracer field ``C`` -are provided, then the tracer advection term becomes - -```math -\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) - + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} C \right ) - + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{U} c \right ) \, . -``` - -Notice that the term ``\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{U} C \right )`` -is neglected: only the terms describing the advection of resolved tracer by the background -velocity field and the advection of background tracer by the resolved velocity field are included. -An analogous statement holds for the advection of background momentum by the resolved -velocity field. -Other possible terms associated with the Coriolis force, buoyancy, turbulence closures, -and surface waves acting on background fields are neglected. - -!!! compat "Model compatibility" - `BackgroundFields` are only supported by [`NonhydrostaticModel`](@ref). - -## Specifying background fields - -`BackgroundField`s are defined by functions of ``(x, y, z, t)`` and optional parameters. A -simple example is - -```jldoctest -using Oceananigans - -U(x, y, z, t) = 0.2 * z - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) - -model = NonhydrostaticModel(grid = grid, background_fields = (u=U,)) - -model.background_fields.velocities.u - -# output -FunctionField located at (Face, Center, Center) -├── func: U (generic function with 1 method) -├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo -├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days) -└── parameters: nothing -``` - -`BackgroundField`s are specified by passing them to the kwarg `background_fields` -in the `NonhydrostaticModel` constructor. The kwarg `background_fields` expects -a `NamedTuple` of fields, which are internally sorted into `velocities` and `tracers`, -wrapped in `FunctionField`s, and assigned their appropriate locations. - -`BackgroundField`s with parameters require using the `BackgroundField` wrapper: - -```jldoctest moar_background -using Oceananigans - -parameters = (α=3.14, N=1.0, f=0.1) - -# Background fields are defined via function of x, y, z, t, and optional parameters -U(x, y, z, t, α) = α * z -B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z - -U_field = BackgroundField(U, parameters=parameters.α) -B_field = BackgroundField(B, parameters=parameters) - -# output -BackgroundField{typeof(B), @NamedTuple{α::Float64, N::Float64, f::Float64}} -├── func: B (generic function with 1 method) -└── parameters: (α = 3.14, N = 1.0, f = 0.1) -``` - -When inserted into `NonhydrostaticModel`, we get - -```jldoctest moar_background -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) - -model = NonhydrostaticModel(grid = grid, background_fields = (u=U_field, b=B_field), - tracers=:b, buoyancy=BuoyancyTracer()) - -model.background_fields.tracers.b - -# output -FunctionField located at (Center, Center, Center) -├── func: B (generic function with 1 method) -├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo -├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days) -└── parameters: (α = 3.14, N = 1.0, f = 0.1) -``` diff --git a/docs/src/models/boundary_conditions.md b/docs/src/models/boundary_conditions.md deleted file mode 100644 index 2f03f15e69..0000000000 --- a/docs/src/models/boundary_conditions.md +++ /dev/null @@ -1,620 +0,0 @@ -# [Boundary conditions](@id model_step_bcs) - -Boundary conditions are intimately related to the grid topology, and only -need to be considered in directions with `Bounded` topology or across immersed boundaries. -In `Bounded` directions, tracer and momentum fluxes are conservative or "zero flux" -by default. Non-default boundary conditions are therefore required to specify non-zero fluxes -of tracers and momentum across `Bounded` directions, and across immersed boundaries -when using `ImmersedBoundaryGrid`. - -See [Numerical implementation of boundary conditions](@ref numerical_bcs) for more details. - -### Example: no-slip conditions on every boundary - -```@meta -DocTestSetup = quote - using Oceananigans - - using Random - Random.seed!(1234) -end -``` - -```jldoctest bcsintro -julia> using Oceananigans - -julia> grid = RectilinearGrid(size=(16, 16, 16), x=(0, 2π), y=(0, 1), z=(0, 1), topology=(Periodic, Bounded, Bounded)) -16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo -├── Periodic x ∈ [0.0, 6.28319) regularly spaced with Δx=0.392699 -├── Bounded y ∈ [0.0, 1.0] regularly spaced with Δy=0.0625 -└── Bounded z ∈ [0.0, 1.0] regularly spaced with Δz=0.0625 - -julia> no_slip_bc = ValueBoundaryCondition(0.0) -ValueBoundaryCondition: 0.0 -``` - -A "no-slip" [`BoundaryCondition`](@ref) specifies that velocity components tangential to `Bounded` -directions decay to `0` at the boundary, leading to a viscous loss of momentum. - -```jldoctest bcsintro -julia> no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc); - -julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=no_slip_field_bcs, v=no_slip_field_bcs, w=no_slip_field_bcs)) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -└── coriolis: Nothing - -julia> model.velocities.u.boundary_conditions -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: PeriodicBoundaryCondition -├── east: PeriodicBoundaryCondition -├── south: ValueBoundaryCondition: 0.0 -├── north: ValueBoundaryCondition: 0.0 -├── bottom: ValueBoundaryCondition: 0.0 -├── top: ValueBoundaryCondition: 0.0 -└── immersed: Nothing -``` - -Boundary conditions are passed to `FieldBoundaryCondition` to build boundary conditions for each -field individually, and then onto the model constructor (here `NonhydrotaticModel`) via the -keyword argument `boundary_conditions`. -The model constructor then "interprets" the input and builds appropriate boundary conditions -for the grid `topology`, given the user-specified `no_slip` default boundary condition for `Bounded` -directions. In the above example, note that the `west` and `east` boundary conditions are `PeriodicBoundaryCondition` -because the `x`-topology of the grid is `Periodic`. - -### Example: specifying boundary conditions on individual boundaries - -To specify no-slip boundary conditions on every `Bounded` direction _except_ -the surface, we write - -```jldoctest bcsintro -julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)); - -julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs)); - -julia> model.velocities.u.boundary_conditions -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: PeriodicBoundaryCondition -├── east: PeriodicBoundaryCondition -├── south: ValueBoundaryCondition: 0.0 -├── north: ValueBoundaryCondition: 0.0 -├── bottom: ValueBoundaryCondition: 0.0 -├── top: FluxBoundaryCondition: Nothing -└── immersed: Nothing - -julia> model.velocities.v.boundary_conditions -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: PeriodicBoundaryCondition -├── east: PeriodicBoundaryCondition -├── south: OpenBoundaryCondition{Nothing}: Nothing -├── north: OpenBoundaryCondition{Nothing}: Nothing -├── bottom: ValueBoundaryCondition: 0.0 -├── top: FluxBoundaryCondition: Nothing -└── immersed: Nothing -``` - -Now both `u` and `v` have `FluxBoundaryCondition(nothing)` at the `top` boundary, which is `Oceananigans` lingo -for "no-flux boundary condition". - -## Boundary condition classifications - -There are three primary boundary condition classifications: - -1. `FluxBoundaryCondition` specifies fluxes directly. - - Some applications of `FluxBoundaryCondition` are: - * surface momentum fluxes due to wind, or "wind stress"; - * linear or quadratic bottom drag; - * surface temperature fluxes due to heating or cooling; - * surface salinity fluxes due to precipitation and evaporation; - * relaxation boundary conditions that restores a field to some boundary distribution - over a given time-scale. - -2. `ValueBoundaryCondition` (Dirichlet) specifies the value of a field on - the given boundary, which when used in combination with a turbulence closure - results in a flux across the boundary. - - _Note_: Do not use `ValueBoundaryCondition` on a wall-normal velocity component - (see the note below about `ImpenetrableBoundaryCondition`). - - Some applications of `ValueBoundaryCondition` are: - * no-slip boundary condition for wall-tangential velocity components via `ValueBoundaryCondition(0)`; - * surface temperature distribution, where heat fluxes in and out of the domain - at a rate controlled by the near-surface temperature gradient and the temperature diffusivity; - * constant velocity tangential to a boundary as in a driven-cavity flow (for example), - where the top boundary is moving. Momentum will flux into the domain do the difference - between the top boundary velocity and the interior velocity, and the prescribed viscosity. - -3. `GradientBoundaryCondition` (Neumann) specifies the gradient of a field on a boundary. - For example, if there is a known `diffusivity`, we can express `FluxBoundaryCondition(flux)` - using `GradientBoundaryCondition(-flux / diffusivity)` (aka "Neumann" boundary condition). - -In addition to these primary boundary conditions, `ImpenetrableBoundaryCondition` applies to velocity -components in wall-normal directions. - -!!! warn "`ImpenetrableBoundaryCondition`" - `ImpenetrableBoundaryCondition` is internally enforced for fields created inside the model constructor. - As a result, `ImpenetrableBoundaryCondition` is only used for _additional_ velocity components - that are not evolved by a model, such as a velocity component used for (`AdvectiveForcing`)[@ref]. - -Finally, note that `Periodic` boundary conditions are internally enforced for `Periodic` directions, -and `DefaultBoundaryCondition`s may exist before boundary conditions are "materialized" by a model. - -## Default boundary conditions - -The default boundary condition in `Bounded` directions is no-flux, or `FluxBoundaryCondition(nothing)`. -The default boundary condition can be changed by passing a positional argument to `FieldBoundaryConditions`, -as in - -```jldoctest -julia> no_slip_bc = ValueBoundaryCondition(0.0) -ValueBoundaryCondition: 0.0 - -julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -├── east: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -├── south: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -├── north: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -├── bottom: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -├── top: FluxBoundaryCondition: Nothing -└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) -``` - -## Boundary condition structures - -Oceananigans uses a hierarchical structure to express boundary conditions: - -1. Each boundary of each field has one [`BoundaryCondition`](@ref) -2. Each field has seven [`BoundaryCondition`](@ref) (`west`, `east`, `south`, `north`, `bottom`, `top` and - `immersed`) -3. A set of `FieldBoundaryConditions`, up to one for each field, are grouped into a `NamedTuple` and passed - to the model constructor. - -## Specifying boundary conditions for a model - -Boundary conditions are defined at model construction time by passing a `NamedTuple` of `FieldBoundaryConditions` -specifying non-default boundary conditions for fields such as velocities and tracers. - -Fields for which boundary conditions are not specified are assigned a default boundary conditions. - -A few illustrations are provided below. See the examples for -further illustrations of boundary condition specification. - -## Creating individual boundary conditions with `BoundaryCondition` - -Boundary conditions may be specified with constants, functions, or arrays. -In this section we illustrate usage of the different [`BoundaryCondition`](@ref) constructors. - -### 1. Constant `Value` (Dirchlet) boundary condition - -```jldoctest bcs -julia> constant_T_bc = ValueBoundaryCondition(20.0) -ValueBoundaryCondition: 20.0 -``` - -A constant [`Value`](@ref) boundary condition can be used to specify constant tracer (such as temperature), -or a constant _tangential_ velocity component at a boundary. Note that boundary conditions on the -_normal_ velocity component must use the [`Open`](@ref) boundary condition type. - -Finally, note that `ValueBoundaryCondition(condition)` is an alias for `BoundaryCondition(Value, condition)`. - -### 2. Constant `Flux` boundary condition - -```jldoctest -julia> ρ₀ = 1027; # Reference density [kg/m³] - -julia> τₓ = 0.08; # Wind stress [N/m²] - -julia> wind_stress_bc = FluxBoundaryCondition(-τₓ/ρ₀) -FluxBoundaryCondition: -7.78968e-5 -``` - -A constant [`Flux`](@ref) boundary condition can be imposed on tracers and tangential velocity components -that can be used, for example, to specify cooling, heating, evaporation, or wind stress at the ocean surface. - -!!! info "The flux convention in Oceananigans" - `Oceananigans` uses the convention that positive fluxes produce transport in the - _positive_ direction (east, north, and up for ``x``, ``y``, ``z``). - This means, for example, that a _negative_ flux of momentum or velocity at a _top_ - boundary, such as in the above example, produces currents in the _positive_ direction, - because it prescribes a downwards flux of momentum into the domain from the top. - Likewise, a _positive_ temperature flux at the top boundary - causes _cooling_, because it transports heat _upwards_, out of the domain. - Conversely, a positive flux at a _bottom_ boundary acts to increase the interior - values of a quantity. - -### 3. Spatially- and temporally-varying flux - -Boundary conditions may be specified by functions, - -```jldoctest -julia> @inline surface_flux(x, y, t) = cos(2π * x) * cos(t); - -julia> top_tracer_bc = FluxBoundaryCondition(surface_flux) -FluxBoundaryCondition: ContinuousBoundaryFunction surface_flux at (Nothing, Nothing, Nothing) -``` - -!!! info "Boundary condition functions" - By default, a function boundary condition is called with the signature - ```julia - f(ξ, η, t) - ``` - where `t` is time and `ξ, η` are spatial coordinates that vary along the boundary: - * `f(y, z, t)` on `x`-boundaries; - * `f(x, z, t)` on `y`-boundaries; - * `f(x, y, t)` on `z`-boundaries. - Alternative function signatures are specified by keyword arguments to - `BoundaryCondition`, as illustrated in subsequent examples. - -### 4. Spatially- and temporally-varying flux with parameters - -Boundary condition functions may be 'parameterized', - -```jldoctest -julia> @inline wind_stress(x, y, t, p) = - p.τ * cos(p.k * x) * cos(p.ω * t); # function with parameters - -julia> top_u_bc = FluxBoundaryCondition(wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4)) -FluxBoundaryCondition: ContinuousBoundaryFunction wind_stress at (Nothing, Nothing, Nothing) -``` - -!!! info "Boundary condition functions with parameters" - The keyword argument `parameters` above specifies that `wind_stress` is called - with the signature `wind_stress(x, y, t, parameters)`. In principle, `parameters` is arbitrary. - However, relatively simple objects such as floating point numbers or `NamedTuple`s must be used - when running on the GPU. - -### 5. 'Field-dependent' boundary conditions - -Boundary conditions may also depend on model fields. For example, a linear drag boundary condition -is implemented with - -```jldoctest -julia> @inline linear_drag(x, y, t, u) = - 0.2 * u -linear_drag (generic function with 1 method) - -julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, field_dependencies=:u) -FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) -``` - -`field_dependencies` specifies the name of the dependent fields either with a `Symbol` or `Tuple` of `Symbol`s. - -### 6. 'Field-dependent' boundary conditions with parameters - -When boundary conditions depends on fields _and_ parameters, their functions take the form - -```jldoctest -julia> @inline quadratic_drag(x, y, t, u, v, drag_coeff) = - drag_coeff * u * sqrt(u^2 + v^2) -quadratic_drag (generic function with 1 method) - -julia> u_bottom_bc = FluxBoundaryCondition(quadratic_drag, field_dependencies=(:u, :v), parameters=1e-3) -FluxBoundaryCondition: ContinuousBoundaryFunction quadratic_drag at (Nothing, Nothing, Nothing) -``` - -Put differently, `ξ, η, t` come first in the function signature, followed by field dependencies, -followed by `parameters` is `!isnothing(parameters)`. - -### 7. Discrete-form boundary condition with parameters - -Discrete field data may also be accessed directly from boundary condition functions -using the `discrete_form`. For example: - -```jldoctest -@inline filtered_drag(i, j, grid, clock, model_fields) = - @inbounds - 0.05 * (model_fields.u[i-1, j, 1] + 2 * model_fields.u[i, j, 1] + model_fields.u[i-1, j, 1]) - -u_bottom_bc = FluxBoundaryCondition(filtered_drag, discrete_form=true) - -# output -FluxBoundaryCondition: DiscreteBoundaryFunction with filtered_drag -``` - -!!! info "The 'discrete form' for boundary condition functions" - The argument `discrete_form=true` indicates to [`BoundaryCondition`](@ref) that `filtered_drag` - uses the 'discrete form'. Boundary condition functions that use the 'discrete form' - are called with the signature - ```julia - f(i, j, grid, clock, model_fields) - ``` - where `i, j` are grid indices that vary along the boundary, `grid` is `model.grid`, - `clock` is the `model.clock`, and `model_fields` is a `NamedTuple` - containing `u, v, w` and the fields in `model.tracers`. - The signature is similar for ``x`` and ``y`` boundary conditions expect that `i, j` is replaced - with `j, k` and `i, k` respectively. - -### 8. Discrete-form boundary condition with parameters - -```jldoctest -julia> Cd = 0.2; # drag coefficient - -julia> @inline linear_drag(i, j, grid, clock, model_fields, Cd) = @inbounds - Cd * model_fields.u[i, j, 1]; - -julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, discrete_form=true, parameters=Cd) -FluxBoundaryCondition: DiscreteBoundaryFunction linear_drag with parameters 0.2 -``` - -!!! info "Inlining and avoiding bounds-checking in boundary condition functions" - Boundary condition functions should be decorated with `@inline` when running on CPUs for performance reasons. - On the GPU, all functions are force-inlined by default. - In addition, the annotation `@inbounds` should be used when accessing the elements of an array - in a boundary condition function (such as `model_fields.u[i, j, 1]` in the above example). - Using `@inbounds` will avoid a relatively expensive check that the index `i, j, 1` is 'in bounds'. - -### 9. A random, spatially-varying, constant-in-time temperature flux specified by an array - -```jldoctest -julia> Nx = Ny = 16; # Number of grid points. - -julia> Q = randn(Nx, Ny); # temperature flux - -julia> white_noise_T_bc = FluxBoundaryCondition(Q) -FluxBoundaryCondition: 16×16 Matrix{Float64} -``` - -When running on the GPU, `Q` must be converted to a `CuArray`. - -## Building boundary conditions on a field - -To create a set of [`FieldBoundaryConditions`](@ref) for a temperature field, -we write - -```jldoctest -julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0), - bottom = GradientBoundaryCondition(0.01)) -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── bottom: GradientBoundaryCondition: 0.01 -├── top: ValueBoundaryCondition: 20.0 -└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -``` - -If the grid is, e.g., horizontally-periodic, then each horizontal `DefaultBoundaryCondition` -is converted to `PeriodicBoundaryCondition` inside the model's constructor, before assigning the -boundary conditions to temperature `T`. - -In general, boundary condition defaults are inferred from the field location and `topology(grid)`. - -## Specifying model boundary conditions - -To specify non-default boundary conditions, a named tuple of [`FieldBoundaryConditions`](@ref) objects is -passed to the keyword argument `boundary_conditions` in the [`NonhydrostaticModel`](@ref) constructor. -The keys of `boundary_conditions` indicate the field to which the boundary condition is applied. -Below, non-default boundary conditions are imposed on the ``u``-velocity and tracer ``c``. - -```jldoctest -julia> topology = (Periodic, Periodic, Bounded); - -julia> grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology); - -julia> u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(+0.1), - bottom = ValueBoundaryCondition(-0.1)); - -julia> c_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0), - bottom = GradientBoundaryCondition(0.01)); - -julia> model = NonhydrostaticModel(grid=grid, boundary_conditions=(u=u_bcs, c=c_bcs), tracers=:c) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: c -├── closure: Nothing -├── buoyancy: Nothing -└── coriolis: Nothing - -julia> model.velocities.u -16×16×16 Field{Face, Center, Center} on RectilinearGrid on CPU -├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── boundary conditions: FieldBoundaryConditions -│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Value, top: Value, immersed: Nothing -└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19 - └── max=0.0, min=0.0, mean=0.0 - -julia> model.tracers.c -16×16×16 Field{Center, Center, Center} on RectilinearGrid on CPU -├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── boundary conditions: FieldBoundaryConditions -│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Gradient, top: Value, immersed: Nothing -└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19 - └── max=0.0, min=0.0, mean=0.0 -``` - -Notice that the specified non-default boundary conditions have been applied at -top and bottom of both `model.velocities.u` and `model.tracers.c`. - -## Immersed boundary conditions - -Immersed boundary conditions are supported experimentally. A no-slip boundary condition is specified -with - -```@setup immersed_bc -using Oceananigans - -# Save original stderr -original_stderr = stderr - -# Redirect stderr to /dev/null (Unix) or "nul" (Windows) -redirect_stderr(devnull) -``` - -```jldoctest immersed_bc -# Generate a simple ImmersedBoundaryGrid -hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2) -underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded)) -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill)) - -# Create a no-slip boundary condition for velocity fields. -# Note that the no-slip boundary condition is _only_ applied on immersed boundaries. -velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0)) -model = NonhydrostaticModel(; grid, boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs)) - -# Inspect the boundary condition on the vertical velocity: -model.velocities.w.boundary_conditions.immersed - -# output -┌ Warning: The FFT-based pressure_solver for NonhydrostaticModels on ImmersedBoundaryGrid -│ is approximate and will probably produce velocity fields that are divergent -│ adjacent to the immersed boundary. An experimental but improved pressure_solver -│ is available which may be used by writing -│ -│ using Oceananigans.Solvers: ConjugateGradientPoissonSolver -│ pressure_solver = ConjugateGradientPoissonSolver(grid) -│ -│ Please report issues to https://github.com/CliMA/Oceananigans.jl/issues. -└ @ Oceananigans.Models.NonhydrostaticModels ~/Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:55 -ImmersedBoundaryCondition: -├── west: ValueBoundaryCondition: 0.0 -├── east: ValueBoundaryCondition: 0.0 -├── south: ValueBoundaryCondition: 0.0 -├── north: ValueBoundaryCondition: 0.0 -├── bottom: Nothing -└── top: Nothing -``` - -!!! warning "`NonhydrostaticModel` on `ImmersedBoundaryGrid`" - The pressure solver for `NonhydrostaticModel` is approximate, and is unable to produce - a velocity field that is simultaneously divergence-free while also satisfying impenetrability - on the immersed boundary. As a result, simulated dynamics with `NonhydrostaticModel` can - exhibit egregiously unphysical errors and should be interpreted with caution. - -```jldoctest immersed_bc -hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2) -underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded)) -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill)) - -# Create a no-slip boundary condition for velocity fields. -# Note that the no-slip boundary condition is _only_ applied on immersed boundaries. -velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0)) -model = NonhydrostaticModel(; grid, boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs)) - -# output -┌ Warning: The FFT-based pressure_solver for NonhydrostaticModels on ImmersedBoundaryGrid -│ is approximate and will probably produce velocity fields that are divergent -│ adjacent to the immersed boundary. An experimental but improved pressure_solver -│ is available which may be used by writing -│ -│ using Oceananigans.Solvers: ConjugateGradientPoissonSolver -│ pressure_solver = ConjugateGradientPoissonSolver(grid) -│ -│ Please report issues to https://github.com/CliMA/Oceananigans.jl/issues. -└ @ Oceananigans.Models.NonhydrostaticModels ~/Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:55 -NonhydrostaticModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0) -├── grid: 32×32×16 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -└── coriolis: Nothing -``` - -An `ImmersedBoundaryCondition` encapsulates boundary conditions on each potential boundary-facet -of a boundary-adjacent cell. Boundary conditions on specific faces of immersed-boundary-adjacent -cells may also be specified by manually building an `ImmersedBoundaryCondition`: - -```jldoctest immersed_bc -bottom_drag_bc = ImmersedBoundaryCondition(bottom=ValueBoundaryCondition(0)) - -# output -ImmersedBoundaryCondition: -├── west: Nothing -├── east: Nothing -├── south: Nothing -├── north: Nothing -├── bottom: ValueBoundaryCondition: 0 -└── top: Nothing -``` - -The `ImmersedBoundaryCondition` may then be incorporated into the boundary conditions for a -`Field` by prescribing it to the `immersed` boundary label, - -```jldoctest immersed_bc -velocity_bcs = FieldBoundaryConditions(immersed=bottom_drag_bc) - -# output -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Value, top=Nothing -``` - -!!! warning "`ImmersedBoundaryCondition`" - `ImmersedBoundaryCondition` is experimental. - Therefore, one should use it only when a finer level of control over the boundary conditions - at the immersed boundary is required, and the user is familiar with the implementation of boundary - conditions on staggered grids. For all other cases , using the `immersed` argument of - `FieldBoundaryConditions` is preferred. - -A boundary condition that depends on the fields may be prescribed using the `immersed` -keyword argument in [`FieldBoundaryConditions`](@ref). -We illustrate field-dependent boundary conditions with an example that imposes linear bottom drag -on `u` on both the bottom facets of cells adjacent to an immersed boundary, _and_ the bottom boundary -of the underlying grid. - -First we create the boundary condition for the grid's bottom: - -```jldoctest immersed_bc -@inline linear_drag(x, y, t, u) = - 0.2 * u -drag_u = FluxBoundaryCondition(linear_drag, field_dependencies=:u) - -# output -FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) -``` - -Next, we create the immersed boundary condition by adding the argument `z` to `linear_drag` -and imposing drag only on "bottom" facets of cells that neighbor immersed cells: - -```jldoctest immersed_bc -@inline immersed_linear_drag(x, y, z, t, u) = - 0.2 * u -immersed_drag_u = FluxBoundaryCondition(immersed_linear_drag, field_dependencies=:u) - -u_immersed_bc = ImmersedBoundaryCondition(bottom = immersed_drag_u) - -# output -ImmersedBoundaryCondition: -├── west: Nothing -├── east: Nothing -├── south: Nothing -├── north: Nothing -├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction immersed_linear_drag at (Nothing, Nothing, Nothing) -└── top: Nothing -``` - -Finally, we combine the two: - -```jldoctest immersed_bc -u_bcs = FieldBoundaryConditions(bottom = drag_u, immersed = u_immersed_bc) - -# output -Oceananigans.FieldBoundaryConditions, with boundary conditions -├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) -├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) -└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Flux, top=Nothing -``` - -!!! warning "Positional argument requirements" - Note the difference between the arguments required for the function within the `bottom` boundary - condition versus the arguments for the function within the `immersed` boundary condition. E.g., - `x, y, t` in `linear_drag()` versus `x, y, z, t` in `immersed_linear_drag()`. - -```@setup immersed_bc -# Restore original stderr -redirect_stderr(original_stderr) -``` diff --git a/docs/src/models/buoyancy_and_equation_of_state.md b/docs/src/models/buoyancy_and_equation_of_state.md deleted file mode 100644 index a0dd6340d9..0000000000 --- a/docs/src/models/buoyancy_and_equation_of_state.md +++ /dev/null @@ -1,251 +0,0 @@ -# Buoyancy models and equations of state - -The buoyancy option selects how buoyancy is treated in `NonhydrostaticModel`s and -`HydrostaticFreeSurfaceModel`s (`ShallowWaterModel`s do not have that option given the physics of -the model). There are currently three alternatives: - -1. No buoyancy (and no gravity). -2. Evolve buoyancy as a tracer. -3. _Seawater buoyancy_: evolve temperature ``T`` and salinity ``S`` as tracers with a value for the gravitational - acceleration ``g`` and an equation of state of your choosing. - -## No buoyancy - -To turn off buoyancy (and gravity) you can simply pass `buoyancy = nothing` to the model -constructor. For example to create a `NonhydrostaticModel`: - - -```@meta -DocTestSetup = quote - using Oceananigans -end -``` - - -```jldoctest buoyancy -julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)); - -julia> model = NonhydrostaticModel(; grid, buoyancy=nothing) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -└── coriolis: Nothing -``` - -The option `buoyancy = nothing` is the default for [`NonhydrostaticModel`](@ref), so omitting the -`buoyancy` keyword argument from the `NonhydrostaticModel` constructor yields the same: - -```jldoctest buoyancy -julia> model = NonhydrostaticModel(; grid) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -└── coriolis: Nothing -``` - -The same is true for `HydrostaticFreeSurfaceModel`, - -```jldoctest buoyancy -julia> model = HydrostaticFreeSurfaceModel(; grid) -HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: QuasiAdamsBashforth2TimeStepper -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² -│ └── solver: FFTImplicitFreeSurfaceSolver -├── advection scheme: -│ └── momentum: VectorInvariant -├── vertical_coordinate: ZCoordinate -└── coriolis: Nothing -``` - -## Buoyancy as a tracer - -Both `NonhydrostaticModel` and `HydrostaticFreeSurfaceModel` support evolving -a buoyancy tracer by including `:b` in `tracers` and specifying `buoyancy = BuoyancyTracer()`: - -```jldoctest buoyancy -julia> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: b -├── closure: Nothing -├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() -└── coriolis: Nothing -``` - -Similarly for a `HydrostaticFreeSurfaceModel` with buoyancy as a tracer: - -```jldoctest buoyancy -julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) -HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: QuasiAdamsBashforth2TimeStepper -├── tracers: b -├── closure: Nothing -├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() -├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² -│ └── solver: FFTImplicitFreeSurfaceSolver -├── advection scheme: -│ ├── momentum: VectorInvariant -│ └── b: Centered(order=2) -├── vertical_coordinate: ZCoordinate -└── coriolis: Nothing -``` - -## Seawater buoyancy - -`NonhydrostaticModel` and `HydrostaticFreeSurfaceModel` support modeling the buoyancy of seawater -as a function of the gravitational acceleration, the conservative temperature ``T``, and the absolute -salinity ``S``. The relationship between ``T``, ``S``, the geopotential height, and the density -perturbation from a reference value is called the `equation_of_state`. - -Specifying `buoyancy = SeawaterBuoyancy()` returns a buoyancy model with a linear equation of state, -[Earth standard](https://en.wikipedia.org/wiki/Standard_gravity) `gravitational_acceleration = 9.80665` (in -S.I. units ``\text{m}\,\text{s}^{-2}``) and requires to add `:T` and `:S` as tracers: - -```jldoctest buoyancy -julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: (T, S) -├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() -└── coriolis: Nothing -``` - -and the same is true for `HydrostaticFreeSurfaceModel`, - -```jldoctest buoyancy -julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) -HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: QuasiAdamsBashforth2TimeStepper -├── tracers: (T, S) -├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() -├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² -│ └── solver: FFTImplicitFreeSurfaceSolver -├── advection scheme: -│ ├── momentum: VectorInvariant -│ ├── T: Centered(order=2) -│ └── S: Centered(order=2) -├── vertical_coordinate: ZCoordinate -└── coriolis: Nothing -``` - -To model flows near the surface of Europa where `gravitational_acceleration = 1.3` ``\text{m}\,\text{s}^{-2}``, -we might alternatively specify - -```jldoctest buoyancy -julia> buoyancy = SeawaterBuoyancy(gravitational_acceleration=1.3) -SeawaterBuoyancy{Float64}: -├── gravitational_acceleration: 1.3 -└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) - -julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=(:T, :S)) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: (T, S) -├── closure: Nothing -├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() -└── coriolis: Nothing -``` - -for example. - -### Linear equation of state - -To specify the thermal expansion and haline contraction coefficients -``\alpha = 2 \times 10^{-3} \; \text{K}^{-1}`` and ``\beta = 5 \times 10^{-4} \text{psu}^{-1}``, - -```jldoctest -julia> buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-3, haline_contraction=5e-4)) -SeawaterBuoyancy{Float64}: -├── gravitational_acceleration: 9.80665 -└── equation_of_state: LinearEquationOfState(thermal_expansion=0.002, haline_contraction=0.0005) -``` - -### Idealized nonlinear equations of state - -Instead of a linear equation of state, six idealized (second-order) nonlinear equations of state -as described by [Roquet15Idealized](@citet) may be used. These equations of state are provided -via the [SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package. - -```jldoctest buoyancy -julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials - -julia> eos = RoquetEquationOfState(:Freezing) -BoussinesqEquationOfState{Float64}: -├── seawater_polynomial: SecondOrderSeawaterPolynomial{Float64} -└── reference_density: 1024.6 - -julia> eos.seawater_polynomial # the density anomaly -ρ' = 0.7718 Sᴬ - 0.0491 Θ - 0.005027 Θ² - 2.5681e-5 Θ Z + 0.0 Sᴬ² + 0.0 Sᴬ Z + 0.0 Sᴬ Θ - -julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos) -SeawaterBuoyancy{Float64}: -├── gravitational_acceleration: 9.80665 -└── equation_of_state: BoussinesqEquationOfState{Float64} -``` - -### TEOS-10 equation of state - -A high-accuracy 55-term polynomial approximation to the TEOS-10 equation of state suitable for use in -Boussinesq models as described by [Roquet15TEOS](@citet) is implemented in the -[SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package and may be used. - -```jldoctest buoyancy -julia> using SeawaterPolynomials.TEOS10 - -julia> eos = TEOS10EquationOfState() -BoussinesqEquationOfState{Float64}: -├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64} -└── reference_density: 1020.0 -``` - -## The direction of gravitational acceleration - -To simulate gravitational accelerations that don't align with the vertical (`z`) coordinate, -we use `BuoyancyForce(formulation; gravity_unit_vector)`, wherein the buoyancy `formulation` -can be `BuoyancyTracer`, `SeawaterBuoyancy`, etc, in addition to the `gravity_unit_vector`. -For example, - -```jldoctest buoyancy -julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)); - -julia> θ = 45; # degrees - -julia> g̃ = (0, sind(θ), cosd(θ)); - -julia> buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃) -BuoyancyForce: -├── formulation: BuoyancyTracer -└── gravity_unit_vector: (0.0, 0.707107, 0.707107) - -julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=:b) -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: b -├── closure: Nothing -├── buoyancy: BuoyancyTracer with ĝ = (0.0, 0.707107, 0.707107) -└── coriolis: Nothing -``` diff --git a/docs/src/models/coriolis.md b/docs/src/models/coriolis.md deleted file mode 100644 index 013dddac75..0000000000 --- a/docs/src/models/coriolis.md +++ /dev/null @@ -1,112 +0,0 @@ -# Coriolis - -The Coriolis option determines whether the fluid experiences the effect of the Coriolis force, or rotation. Currently -three options are available: no rotation, ``f``-plane, and ``\beta``-plane. - -!!! info "Coriolis vs. rotation" - If you are wondering why this option is called "Coriolis" it is because rotational effects could include the - Coriolis and centripetal forces, both of which arise in non-inertial reference frames. But here the model only - considers the Coriolis force. - -## No rotation - -By default there is no rotation. This can be made explicit by passing `coriolis = nothing` to a model constructor. - -## Traditional ``f``-plane - -To set up an ``f``-plane with, for example, Coriolis parameter ``f = 10^{-4} \text{s}^{-1}`` - -```@meta -DocTestSetup = quote - using Oceananigans -end -``` - -```jldoctest -julia> coriolis = FPlane(f=1e-4) -FPlane{Float64}(f=0.0001) -``` - -An ``f``-plane can also be specified at some latitude on a spherical planet with a planetary rotation rate. For example, -to specify an ``f``-plane at a latitude of ``\varphi = 45°\text{N}`` on Earth which has a rotation rate of -``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` - -```jldoctest -julia> coriolis = FPlane(rotation_rate=7.292115e-5, latitude=45) -FPlane{Float64}(f=0.000103126) -``` - -in which case the value of ``f`` is given by ``2\Omega\sin\varphi``. - -## Coriolis term for constant rotation in a Cartesian coordinate system - -One can use `ConstantCartesianCoriolis` to set up a Coriolis acceleration term where the Coriolis parameter -is constant and the rotation axis is arbitrary. For example, with -``\boldsymbol{f} = (0, f_y, f_z) = (0, 2, 1) \times 10^{-4} \text{s}^{-1}``, - -```jldoctest -julia> coriolis = ConstantCartesianCoriolis(fx=0, fy=2e-4, fz=1e-4) -ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04 -``` - -Or alternatively, the same result can be achieved by specifying the magnitude of the Coriolis -frequency `f` and the `rotation_axis`. So another way to get a Coriolis acceleration with the same -values is: - -```jldoctest -julia> rotation_axis = (0, 2e-4, 1e-4)./√(2e-4^2 + 1e-4^2) # rotation_axis has to be a unit vector -(0.0, 0.8944271909999159, 0.4472135954999579) - -julia> coriolis = ConstantCartesianCoriolis(f=√(2e-4^2+1e-4^2), rotation_axis=rotation_axis) -ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04 -``` - -An ``f``-plane with non-traditional Coriolis terms can also be specified at some latitude on a spherical planet -with a planetary rotation rate. For example, to specify an ``f``-plane at a latitude of ``\varphi = 45°\text{N}`` -on Earth which has a rotation rate of ``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` - -```jldoctest -julia> coriolis = ConstantCartesianCoriolis(rotation_rate=7.292115e-5, latitude=45) -ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 1.03e-04, fz = 1.03e-04 -``` - -in which case ``f_z = 2\Omega\sin\varphi`` and ``f_y = 2\Omega\cos\varphi``. - -## Traditional ``\beta``-plane - -To set up a ``\beta``-plane the background rotation rate ``f_0`` and the ``\beta`` parameter must be specified. For example, -a ``\beta``-plane with ``f_0 = 10^{-4} \text{s}^{-1}`` and ``\beta = 1.5 \times 10^{-11} \text{s}^{-1}\text{m}^{-1}`` can be -set up with - -```jldoctest -julia> coriolis = BetaPlane(f₀=1e-4, β=1.5e-11) -BetaPlane{Float64}(f₀=0.0001, β=1.5e-11) -``` - -Alternatively, a ``\beta``-plane can also be set up at some latitude on a spherical planet with a planetary rotation rate -and planetary radius. For example, to specify a ``\beta``-plane at a latitude of ``\varphi = 10^\circ{S}`` on Earth -which has a rotation rate of ``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` and a radius of ``R = 6,371 \text{km}`` - -```jldoctest -julia> coriolis = BetaPlane(rotation_rate=7.292115e-5, latitude=-10, radius=6371e3) -BetaPlane{Float64}(f₀=-2.53252e-5, β=2.25438e-11) -``` - -in which case ``f_0 = 2\Omega\sin\varphi`` and ``\beta = 2\Omega\cos\varphi / R``. - -## Non-traditional ``\beta``-plane - -A non-traditional ``\beta``-plane requires either 5 parameters (by default Earth's radius and -rotation rate are used): - -```jldoctest -julia> NonTraditionalBetaPlane(fz=1e-4, fy=2e-4, β=4e-11, γ=-8e-11) -NonTraditionalBetaPlane{Float64}(fz = 1.00e-04, fy = 2.00e-04, β = 4.00e-11, γ = -8.00e-11, R = 6.37e+06) -``` - -or the rotation rate, radius, and latitude: - -```jldoctest -julia> NonTraditionalBetaPlane(rotation_rate=5.31e-5, radius=252.1e3, latitude=10) -NonTraditionalBetaPlane{Float64}(fz = 1.84e-05, fy = 1.05e-04, β = 4.15e-10, γ = -1.46e-10, R = 2.52e+05) -``` diff --git a/docs/src/models/forcing_functions.md b/docs/src/models/forcing_functions.md deleted file mode 100644 index 504f18539b..0000000000 --- a/docs/src/models/forcing_functions.md +++ /dev/null @@ -1,366 +0,0 @@ -# [Forcing functions](@id forcing_functions) - -"Forcings" are user-defined terms appended to right-hand side of -the momentum or tracer evolution equations. In `Oceananigans`, momentum -and tracer forcings are defined via julia functions. `Oceananigans` includes -an interface for implementing forcing functions that depend on spatial coordinates, -time, model velocity and tracer fields, and external parameters. - -```@meta -DocTestSetup = quote - using Oceananigans -end -``` - -Forcings are added to models by passing a `NamedTuple` of functions or forcing objects -to the `forcing` keyword argument in `NonhydrostaticModel`'s constructor. -By default, momentum and tracer forcing functions are assumed to be functions of -`x, y, z, t`. A basic example is - -```jldoctest -u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,)) - -model.forcing.u - -# output -ContinuousForcing{Nothing} at (Face, Center, Center) -├── func: u_forcing (generic function with 1 method) -├── parameters: nothing -└── field dependencies: () -``` - -More general forcing functions are built via the `Forcing` constructor -described below. `Oceananigans` also provides two convenience types: - - * `Relaxation` for damping terms that restore a field to a - target distribution outside of a masked region of space. `Relaxation` can be - used to implement sponge layers near the boundaries of a domain. - * `AdvectiveForcing` for advecting individual quantities by a separate or - "slip" velocity relative to both the prognostic model velocity field and any - `BackgroundField` velocity field. - -## The `Forcing` constructor - -The `Forcing` constructor provides an interface for specifying forcing functions that - -1. Depend on external parameters; and -2. Depend on model fields at the `x, y, z` location that forcing is applied; and/or -3. Require access to discrete model data. - -### Forcing functions with external parameters - -Most forcings involve external, changeable parameters. -Here are two examples of `forcing_func`tions that depend on -_(i)_ a single scalar parameter `s`, and _(ii)_ a `NamedTuple` of parameters, `p`: - -```jldoctest parameterized_forcing -# Forcing that depends on a scalar parameter `s` -u_forcing_func(x, y, z, t, s) = s * z - -u_forcing = Forcing(u_forcing_func, parameters=0.1) - -# Forcing that depends on a `NamedTuple` of parameters `p` -T_forcing_func(x, y, z, t, p) = - p.μ * exp(z / p.λ) * cos(p.k * x) * sin(p.ω * t) - -T_forcing = Forcing(T_forcing_func, parameters=(μ=1, λ=0.5, k=2π, ω=4π)) - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing, T=T_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) - -model.forcing.T - -# output -ContinuousForcing{@NamedTuple{μ::Int64, λ::Float64, k::Float64, ω::Float64}} at (Center, Center, Center) -├── func: T_forcing_func (generic function with 1 method) -├── parameters: (μ = 1, λ = 0.5, k = 6.283185307179586, ω = 12.566370614359172) -└── field dependencies: () -``` - -```jldoctest parameterized_forcing -model.forcing.u - -# output -ContinuousForcing{Float64} at (Face, Center, Center) -├── func: u_forcing_func (generic function with 1 method) -├── parameters: 0.1 -└── field dependencies: () -``` - -In this example, the objects passed to the `parameters` keyword in the construction of -`u_forcing` and `T_forcing` -- a floating point number for `u_forcing`, and a `NamedTuple` -of parameters for `T_forcing` -- are passed on to `u_forcing_func` and `T_forcing_func` when -they are called during time-stepping. The object passed to `parameters` is in principle arbitrary. -However, if using the GPU, then `typeof(parameters)` may be restricted by the requirements -of GPU-compiliability. - -### Forcing functions that depend on model fields - -Forcing functions may depend on model fields (velocity, tracers or auxiliary fields) evaluated at the `x, y, z` where forcing is applied. -Here's a somewhat non-sensical example: - -```jldoctest field_dependent_forcing -# Forcing that depends on the velocity fields `u`, `v`, and `w` -w_forcing_func(x, y, z, t, u, v, w) = - (u^2 + v^2 + w^2) / 2 - -w_forcing = Forcing(w_forcing_func, field_dependencies=(:u, :v, :w)) - -# Forcing that depends on salinity `S` and a scalar parameter -S_forcing_func(x, y, z, t, S, μ) = - μ * S - -S_forcing = Forcing(S_forcing_func, parameters=0.01, field_dependencies=:S) - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(grid=grid, forcing=(w=w_forcing, S=S_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) - -model.forcing.w - -# output -ContinuousForcing{Nothing} at (Center, Center, Face) -├── func: w_forcing_func (generic function with 1 method) -├── parameters: nothing -└── field dependencies: (:u, :v, :w) -``` - -```jldoctest field_dependent_forcing -model.forcing.S - -# output -ContinuousForcing{Float64} at (Center, Center, Center) -├── func: S_forcing_func (generic function with 1 method) -├── parameters: 0.01 -└── field dependencies: (:S,) -``` - -The `field_dependencies` arguments follow `x, y, z, t` in the forcing `func`tion in -the order they are specified in `Forcing`. -If both `field_dependencies` and `parameters` are specified, then the `field_dependencies` -arguments follow `x, y, z, t`, and `parameters` follow `field_dependencies`. - -Model fields that arise in the arguments of continuous `Forcing` `func`tions are -automatically interpolated to the staggered grid location at which the forcing is applied. - -### "Discrete form" forcing functions - -"Discrete form" forcing functions are either called with the signature - -```julia -func(i, j, k, grid, clock, model_fields) -``` - -or the parameterized form - -```julia -func(i, j, k, grid, clock, model_fields, parameters) -``` - -Discrete form forcing functions can access the entirety of model field -data through the argument `model_fields`. The object `model_fields` is a `NamedTuple` -whose properties include the velocity fields `model_fields.u`, `model_fields.v`, -`model_fields.w` and all fields in `model.tracers`. - -Using discrete forcing functions may require understanding the -staggered arrangement of velocity fields and tracers in `Oceananigans`. -Here's a slightly non-sensical example in which the vertical derivative of a buoyancy -tracer is used as a time-scale for damping the u-velocity field: - -```jldoctest discrete_forcing -# A damping term that depends on a "local average": -local_average(i, j, k, grid, c) = @inbounds (c[i, j, k] + c[i-1, j, k] + c[i+1, j, k] + - c[i, j-1, k] + c[i, j+1, k] + - c[i, j, k-1] + c[i, j, k+1]) / 7 - -b_forcing_func(i, j, k, grid, clock, model_fields) = - local_average(i, j, k, grid, model_fields.b) - -b_forcing = Forcing(b_forcing_func, discrete_form=true) - -# A term that damps the local velocity field in the presence of stratification -using Oceananigans.Operators: ∂zᶠᶜᶠ, ℑxzᶠᵃᶜ - -function u_forcing_func(i, j, k, grid, clock, model_fields, ε) - # The vertical derivative of buoyancy, interpolated to the u-velocity location: - N² = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶠᶜᶠ, model_fields.b) - - # Set to zero in unstable stratification where N² < 0: - N² = max(N², zero(typeof(N²))) - - return @inbounds - ε * sqrt(N²) * model_fields.u[i, j, k] -end - -u_forcing = Forcing(u_forcing_func, discrete_form=true, parameters=1e-3) - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(grid=grid, tracers=:b, buoyancy=BuoyancyTracer(), forcing=(u=u_forcing, b=b_forcing)) - -model.forcing.b - -# output -DiscreteForcing{Nothing} -├── func: b_forcing_func (generic function with 1 method) -└── parameters: nothing -``` - -```jldoctest discrete_forcing -model.forcing.u - -# output -DiscreteForcing{Float64} -├── func: u_forcing_func (generic function with 1 method) -└── parameters: 0.001 -``` - -The annotation `@inbounds` is crucial for performance when accessing array indices -of the fields in `model_fields`. - -## `Relaxation` - -`Relaxation` defines a special forcing function that restores a field at a specified `rate` to -a `target` distribution, within a region uncovered by a `mask`ing function. -`Relaxation` is useful for implementing sponge layers, as shown in the second example. - -The following code constructs a model in which all components -of the velocity field are damped to zero everywhere on a time-scale of 1000 seconds, or ~17 minutes: - -```jldoctest -damping = Relaxation(rate = 1/1000) - -grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) -model = NonhydrostaticModel(grid=grid, forcing=(u=damping, v=damping, w=damping)) - -model.forcing.w - -# output -ContinuousForcing{Nothing} at (Center, Center, Face) -├── func: Relaxation(rate=0.001, mask=1, target=0) -├── parameters: nothing -└── field dependencies: (:w,) -``` - -The constructor for `Relaxation` accepts the keyword arguments `mask`, and `target`, -which specify a `mask(x, y, z)` function that multiplies the forcing, and a `target(x, y, z)` -distribution for the quantity in question. By default, `mask` uncovered the whole domain -and `target` restores the field in question to 0 - -We illustrate usage of `mask` and `target` by implementing a sponge layer that relaxes -velocity fields to zero and restores temperature to a linear gradient in the bottom -1/10th of the domain: - -```jldoctest sponge_layer -grid = RectilinearGrid(size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(-1, 0)) - - damping_rate = 1/100 # relax fields on a 100 second time-scale -temperature_gradient = 0.001 # ⁰C m⁻¹ - surface_temperature = 20 # ⁰C (at z=0) - -target_temperature = LinearTarget{:z}(intercept=surface_temperature, gradient=temperature_gradient) - bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10) - -uvw_sponge = Relaxation(rate=damping_rate, mask=bottom_mask) - T_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=target_temperature) - -model = NonhydrostaticModel(grid=grid, forcing=(u=uvw_sponge, v=uvw_sponge, w=uvw_sponge, T=T_sponge), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) - -model.forcing.u - -# output -ContinuousForcing{Nothing} at (Face, Center, Center) -├── func: Relaxation(rate=0.01, mask=exp(-(z + 1.0)^2 / (2 * 0.1^2)), target=0) -├── parameters: nothing -└── field dependencies: (:u,) -``` - -```jldoctest sponge_layer -model.forcing.T - -# output -ContinuousForcing{Nothing} at (Center, Center, Center) -├── func: Relaxation(rate=0.01, mask=exp(-(z + 1.0)^2 / (2 * 0.1^2)), target=20.0 + 0.001 * z) -├── parameters: nothing -└── field dependencies: (:T,) -``` - -## `AdvectiveForcing` - -`AdvectiveForcing` defines a forcing function that represents advection by -a separate or "slip" velocity relative to the prognostic model velocity field. -`AdvectiveForcing` is implemented with native Oceananigans advection operators, -which means that tracers advected by the "flux form" advection term -``𝛁⋅𝐮_{\rm slip} c``. Caution is advised when ``𝐮_{\rm slip}`` is not divergence free. - -As an example, consider a model for sediment settling at a constant rate: - -```jldoctest -using Oceananigans - -r_sediment = 1e-4 # [m] "Fine sand" -ρ_sediment = 1200 # kg m⁻³ -ρ_ocean = 1026 # kg m⁻³ -Δb = 9.81 * (ρ_ocean - ρ_sediment) / ρ_ocean # m s⁻² -ν_molecular = 1.05e-6 # m² s⁻¹ -w_sediment = 2/9 * Δb / ν_molecular * r_sediment^2 # m s⁻¹ - -sinking = AdvectiveForcing(w=w_sediment) - -# output -AdvectiveForcing: -├── u: ZeroField{Int64} -├── v: ZeroField{Int64} -└── w: ConstantField(-0.00352102) -``` - -The three keyword arguments specify the `u`, `v`, and `w` components of the separate -slip velocity field. The default for each `u, v, w` is `ZeroField`. - -Next we consider a dynamically-evolving slip velocity. For this we use `ZFaceField` -with appropriate boundary conditions as our slip velocity: - -```jldoctest sinking -using Oceananigans -using Oceananigans.BoundaryConditions: ImpenetrableBoundaryCondition - -grid = RectilinearGrid(size=(32, 32, 32), x=(-10, 10), y=(-10, 10), z=(-4, 4), - topology=(Periodic, Periodic, Bounded)) - -no_penetration = ImpenetrableBoundaryCondition() -slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), - top=no_penetration, bottom=no_penetration) - -w_slip = ZFaceField(grid, boundary_conditions=slip_bcs) -sinking = AdvectiveForcing(w=w_slip) - -# output -AdvectiveForcing: -├── u: ZeroField{Int64} -├── v: ZeroField{Int64} -└── w: 32×32×33 Field{Center, Center, Face} on RectilinearGrid on CPU -``` - -To compute the slip velocity, we must add a `Callback`to `simulations.callback` that -computes `w_slip` ever iteration: - -```jldoctest sinking -using Oceananigans.BoundaryConditions: fill_halo_regions! - -model = NonhydrostaticModel(; grid, tracers=(:b, :P), forcing=(; P=sinking)) -simulation = Simulation(model; Δt=1, stop_iteration=100) - -# Build abstract operation for slip velocity -b_particle = - 1e-4 # relative buoyancy depends on reference density and initial buoyancy condition -b = model.tracers.b -R = 1e-3 # [m] mean particle radius -ν = 1.05e-6 # [m² s⁻¹] molecular kinematic viscosity of water -w_slip_op = 2/9 * (b - b_particle) / ν * R^2 # Stokes terminal velocity - -function compute_slip_velocity!(sim) - w_slip .= w_slip_op - fill_halo_regions!(w_slip) - return nothing -end - -simulation.callbacks[:slip] = Callback(compute_slip_velocity!) - -# output -Callback of compute_slip_velocity! on IterationInterval(1) -``` diff --git a/docs/src/models/lagrangian_particles.md b/docs/src/models/lagrangian_particles.md deleted file mode 100644 index e2b99571c6..0000000000 --- a/docs/src/models/lagrangian_particles.md +++ /dev/null @@ -1,136 +0,0 @@ -# Lagrangian particle tracking - -Models can keep track of the location and properties of neutrally buoyant particles. Particles are -advected with the flow field using forward Euler time-stepping at every model iteration. - -## Simple particles - -If you just need to keep of particle locations ``(x, y, z)`` then you can construct some Lagrangian particles -using the regular `LagrangianParticles` constructor - -```@meta -DocTestSetup = quote - using Oceananigans -end -``` - -```jldoctest particles -grid = RectilinearGrid(size=(10, 10, 10), extent=(1, 1, 1)) - -Nparticles = 10 - -x₀ = zeros(Nparticles) - -y₀ = rand(Nparticles) - -z₀ = -0.5 * ones(Nparticles) - -lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀) - -# output -10 LagrangianParticles with eltype Particle: -├── 3 properties: (:x, :y, :z) -├── particle-wall restitution coefficient: 1.0 -├── 0 tracked fields: () -└── dynamics: no_dynamics -``` - -then pass it to a model constructor - -```jldoctest particles -model = NonhydrostaticModel(grid=grid, particles=lagrangian_particles) - -# output -NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 10×10×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo -├── timestepper: RungeKutta3TimeStepper -├── advection scheme: Centered(order=2) -├── tracers: () -├── closure: Nothing -├── buoyancy: Nothing -├── coriolis: Nothing -└── particles: 10 LagrangianParticles with eltype Particle and properties (:x, :y, :z) -``` - -!!! warn "Lagrangian particles on GPUs" - Remember to use `CuArray` instead of regular `Array` when storing particle locations and properties on the GPU. - -## Custom particles - -If you want to keep track of custom properties, such as the species or DNA of a Lagrangian particle -representing a microbe in an agent-based model, then you can create your own custom particle type -and pass a `StructArray` to the `LagrangianParticles` constructor. - -```jldoctest particles -using Oceananigans -using StructArrays - -struct LagrangianMicrobe{T, S, D} - x :: T - y :: T - z :: T - species :: S - dna :: D -end - -Nparticles = 3 - -x₀ = zeros(Nparticles) - -y₀ = rand(Nparticles) - -z₀ = -0.5 * ones(Nparticles) - -species = [:rock, :paper, :scissors] - -dna = ["TATACCCC", "CCTAGGAC", "CGATTTAA"] - -particles = StructArray{LagrangianMicrobe}((x₀, y₀, z₀, species, dna)); - -lagrangian_particles = LagrangianParticles(particles) - -# output -3 LagrangianParticles with eltype LagrangianMicrobe: -├── 5 properties: (:x, :y, :z, :species, :dna) -├── particle-wall restitution coefficient: 1.0 -├── 0 tracked fields: () -└── dynamics: no_dynamics -``` - -!!! warn "Custom properties on GPUs" - Not all data types can be passed to GPU kernels. If you intend to advect particles on the GPU make sure - particle properties consist of only simple data types. The symbols and strings in this example won't - work on the GPU. - -## Writing particle properties to disk - -Particle properties can be written to disk using JLD2 or NetCDF. - -When writing to JLD2 you can pass `model.particles` as part of the named tuple of outputs. - -```@setup particles -using Oceananigans -using NCDatasets -grid = RectilinearGrid(size=(10, 10, 10), extent=(1, 1, 1)) -Nparticles = 3 -x₀ = zeros(Nparticles) -y₀ = rand(Nparticles) -z₀ = -0.5 * ones(Nparticles) -lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀) -model = NonhydrostaticModel(; grid, particles=lagrangian_particles) -``` - -```@example particles -JLD2Writer(model, (; model.particles), filename="particles", schedule=TimeInterval(15)) -``` - -When writing to NetCDF you should write particles to a separate file as the NetCDF dimensions differ for -particle trajectories. You can just pass `model.particles` straight to `NetCDFWriter`: - -```@example particles -NetCDFWriter(model, (; model.particles), filename="particles.nc", schedule=TimeInterval(15)) -``` - -!!! warn "Outputting custom particle properties to NetCDF" - NetCDF does not support arbitrary data types. If you need to write custom particle properties to disk - that are not supported by NetCDF then you should use JLD2 (which should support almost any Julia data type). diff --git a/docs/src/models/models_overview.md b/docs/src/models/models_overview.md deleted file mode 100644 index 75ce544750..0000000000 --- a/docs/src/models/models_overview.md +++ /dev/null @@ -1,208 +0,0 @@ -# Models - -In general in Oceananigans, the "model" object serves two main purposes: - * models store the configuration of a set of discrete equations. The discrete equations imply rules for evolving - prognostic variables, and computing diagnostic variables from the prognostic state. - * models provide a container for the prognostic and diagnostic state of those discrete equations at a particular time. - -## Two Oceananigans models for ocean simulations - -In addition to defining the abstract concept of a "model" that can be used with [Simulation](@ref), -Oceananigans provides two mature model implementations for simulating ocean-flavored fluid dynamics. -Both of these integrate the Navier-Stokes equations within the Boussinesq approximation -(we call these the "Boussinesq equations" for short): the [NonhydrostaticModel](@ref) and -the [HydrostaticFreeSurfaceModel](@ref). - -The [NonhydrostaticModel](@ref) integrates the Boussinesq equations _without_ making the hydrostatic approximation, -and therefore possessing a prognostic vertical momentum equation. The NonhydrostaticModel is useful for simulations -that resolve three-dimensional turbulence, such as large eddy simulations on [RectilinearGrid](@ref) with grid -spacings of O(1 m), as well as direct numerical simulation. The NonhydrostaticModel may also be used for -idealized classroom problems, as in the [two-dimensional turbulence example](@ref "Two dimensional turbulence example"). - -The [HydrostaticFreeSurfaceModel](@ref) integrates the hydrostatic or "primitive" Boussinesq equations -with a free surface on its top boundary. The hydrostatic approximation allosw the HydrostaticFreeSurfaceModel -to achieve much higher efficiency in simulations on curvilinear grids used for large-scale regional or global -simulations such as [LatitudeLongitudeGrid](@ref), [TripolarGrid](@ref), [ConformalCubedSphereGrid](@ref), -and other [OrthogonalSphericalShellGrid](@ref)s such as [RotatedLatitudeLongitudeGrid](@ref Oceananigans.OrthogonalSphericalShellGrids.RotatedLatitudeLongitudeGrid). -Because they span larger domains, simulations with the HydrostaticFreeSurfaceModel also usually involve coarser -grid spacings of O(30 m) up to O(100 km). Such coarse-grained simulations are usually paired with more elaborate -turbulence closures or "parameterizations" than small-scale simulations with NonhydrostaticModel, such as the -vertical mixing schemes [CATKEVerticalDiffusivity](@ref), [RiBasedVerticalDiffusivity](@ref), and -[TKEDissipationVerticalDiffusivity](@ref), and the mesoscale turbulence closure -[IsopycnalSkewSymmetricDiffusivity](@ref) (a.k.a. "Gent-McWilliams plus Redi"). - -A third, experimental [ShallowWaterModel](@ref) solves the shallow water equations. - -### Configuring NonhydrostaticModel with keyword arguments - -To illustrate the specification of discrete equations, consider first the docstring for [NonhydrostaticModel](@ref), - -```@docs -NonhydrostaticModel -``` - -The configuration operations for NonhydrostaticModel include "discretization choices", such as the `advection` scheme, -as well as aspects of the continuous underlying equations, such as the formulation of the buoyancy force. - -For our first example, we build the default `NonhydrostaticModel` (which is quite boring): - -```@example -using Oceananigans -grid = RectilinearGrid(size=(8, 8, 8), extent=(8, 8, 8)) -nh = NonhydrostaticModel(; grid) -``` - -The default `NonhydrostaticModel` has no tracers, no buoyancy force, no Coriolis force, and a second-order advection scheme. -We next consider a slightly more exciting NonhydrostaticModel configured with a WENO advection scheme, -the temperature/salinity-based [SeawaterBuoyancy](@ref), a boundary condition on the zonal momentum, -and a passive tracer forced by a cooked-up surface flux called "c": - -```@example first_model -using Oceananigans - -grid = RectilinearGrid(size=(128, 128), halo=(5, 5), x=(0, 256), z=(-128, 0), - topology = (Periodic, Flat, Bounded)) - -# Numerical method and physics choices -advection = WENO(order=9) # ninth‑order upwind for momentum and tracers -buoyancy = BuoyancyTracer() -coriolis = FPlane(f=1e-4) - -τx = - 8e-5 # m² s⁻² (τₓ < 0 ⟹ eastward wind stress) -u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(τx)) - -@inline Jc(x, t, Lx) = cos(2π / Lx * x) -c_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Jc, parameters=grid.Lx)) - -model = NonhydrostaticModel(; grid, advection, buoyancy, coriolis, - tracers = (:b, :c), - boundary_conditions = (; u=u_bcs, c=c_bcs)) -``` - -### Mutation of the model state - -In addition to providing an interface for configuring equations, models also -store the prognostic and diagnostic state associated with the solution to those equations. -Models thus also provide an interface for "setting" or fixing the prognostic state, which is typically -invoked to determine the initial conditions of a simulation. -To illustrate this we consider setting the above model to a stably-stratified and noisy condition: - -```@example first_model -N² = 1e-5 -bᵢ(x, z) = N² * z + 1e-6 * randn() -uᵢ(x, z) = 1e-3 * randn() -set!(model, b=bᵢ, u=uᵢ, w=uᵢ) - -model.tracers.b -``` - -Invoking `set!` above determine the model tracer `b` and the velocity components `u` and `w`. -`set!` also computes the diagnostic state of a model, which in the case of `NonhydrostaticModel` includes -the nonhydrostatic component of pressure, - -```@example first_model -model.pressures.pNHS -``` - -### Evolving models in time - -Model may be integrated or "stepped forward" in time by calling `time_step!(model, Δt)`, where `Δt` is the time step -and thus advancing the `model.clock`: - -```@example first_model -time_step!(model, 1) -model.clock -``` - -However, users are strongly encouraged to use the [`Simulation`](@ref) interface to manage time-stepping -along with other activities, like monitoring progress, writing output to disk, and more. - -```@example first_model -simulation = Simulation(model, Δt=1, stop_iteration=10) -run!(simulation) - -simulation -``` - -## Using the HydrostaticFreeSurfaceModel - -The HydrostaticFreeSurfaceModel has a similar interface as the NonhydrostaticModel, - -```@example -using Oceananigans -grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)) -model = HydrostaticFreeSurfaceModel(; grid) # default free surface, no tracers -``` - -The full array of keyword arguments used to configure a HydrostaticFreeSurfaceModel are detailed -in the docstring for [HydrostaticFreeSurfaceModel](@ref), - -```@docs -HydrostaticFreeSurfaceModel -``` - -A bit more involved HydrostaticFreeSurfaceModel example: - -```@example second_model -using Oceananigans -using SeawaterPolynomials: TEOS10EquationOfState - -grid = LatitudeLongitudeGrid(size = (180, 80, 10), - longitude = (0, 360), - latitude = (-80, 80), - z = (-1000, 0), - halo = (6, 6, 3)) - -momentum_advection = WENOVectorInvariant() -coriolis = HydrostaticSphericalCoriolis() -equation_of_state = TEOS10EquationOfState() -buoyancy = SeawaterBuoyancy(; equation_of_state) -closure = CATKEVerticalDiffusivity() - -# Generate a zonal wind stress that mimics Earth's mean winds -# with westerlies in mid-latitudes and easterlies near equator and poles -function zonal_wind_stress(λ, φ, t) - # Parameters - τ₀ = 1e-4 # Maximum wind stress magnitude (N/m²) - φ₀ = 30 # Latitude of maximum westerlies (degrees) - dφ = 10 - - # Approximate wind stress pattern - return - τ₀ * (+ exp(-(φ - φ₀)^2 / 2dφ^2) - - exp(-(φ + φ₀)^2 / 2dφ^2) - - 0.3 * exp(-φ^2 / dφ^2)) -end - -u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(zonal_wind_stress)) - -model = HydrostaticFreeSurfaceModel(; grid, momentum_advection, coriolis, closure, buoyancy, - boundary_conditions = (; u=u_bcs), tracers=(:T, :S, :e)) -``` - -Mutating the state of the HydrostaticFreeSurfaceModel works similarly as for the NonhydrostaticModel --- -except that the vertical velocity cannot be `set!`, because vertical velocity is not -prognostic in the hydrostatic equations. - -```@example second_model -using SeawaterPolynomials - -N² = 1e-5 -T₀ = 20 -S₀ = 35 -eos = model.buoyancy.formulation.equation_of_state -α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos) -g = model.buoyancy.formulation.gravitational_acceleration -dTdz = N² / (α * g) -Tᵢ(λ, φ, z) = T₀ + dTdz * z + 1e-3 * T₀ * randn() -uᵢ(λ, φ, z) = 1e-3 * randn() -set!(model, T=Tᵢ, S=S₀, u=uᵢ, v=uᵢ) - -model.tracers.T -``` - -## Where to go next - -- See [Quick start](@ref quick_start) for a compact, end-to-end workflow -- See the Examples gallery for longer tutorials covering specific cases, including large eddy simulation, Kelvin–Helmholtz instability and baroclinic instability. -- Other pages in the Models section: in‑depth pages on buoyancy, forcing, boundary conditions, closures, diagnostics, and output. -- Physics: governing equations and numerical forms for [`NonhydrostaticModel`](@ref) and [`HydrostaticFreeSurfaceModel`](@ref hydrostatic_free_surface_model). diff --git a/docs/src/models/turbulence_closures.md b/docs/src/models/turbulence_closures.md deleted file mode 100644 index 29c2605cfd..0000000000 --- a/docs/src/models/turbulence_closures.md +++ /dev/null @@ -1,103 +0,0 @@ -# [Turbulent diffusivity closures and large eddy simulation models](@id turbulence_closures) - -A turbulent diffusivity closure representing the effects of viscous dissipation and diffusion can -be passed via the `closure` keyword. - -See [turbulence closures](@ref numerical_closures) and [large eddy simulation](@ref numerical_les) for -more details on turbulent diffusivity closures. - -## Constant isotropic diffusivity - -To use constant isotropic values for the viscosity ``\nu`` and diffusivity ``\kappa`` you can use `ScalarDiffusivity`: - -```jldoctest -julia> using Oceananigans.TurbulenceClosures - -julia> closure = ScalarDiffusivity(ν=1e-2, κ=1e-2) -ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.01, κ=0.01) -``` - -## Constant anisotropic diffusivity - -To specify constant values for the horizontal and vertical viscosities, ``\nu_h`` and ``\nu_z``, -and horizontal and vertical diffusivities, ``\kappa_h`` and ``\kappa_z``, you can use -`HorizontalScalarDiffusivity()` and `VerticalScalarDiffusivity()`, e.g., - -```jldoctest -julia> using Oceananigans.TurbulenceClosures - -julia> horizontal_closure = HorizontalScalarDiffusivity(ν=1e-3, κ=2e-3) -HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.001, κ=0.002) - -julia> vertical_closure = VerticalScalarDiffusivity(ν=1e-3, κ=2e-3) -VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.001, κ=0.002) -``` - -After that you can set, e.g., `closure = (horizontal_closure, vertical_closure)` when constructing -the model so that all components will be taken into account when calculating the diffusivity term. -Note that `VerticalScalarDiffusivity` and `HorizontalScalarDiffusivity` are implemented using [different -schemes](https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#horizontal-dissipation) -with different conservation properties. - -## Tracer-specific diffusivities - -You can also set different diffusivities for each tracer in your simulation by passing -a [NamedTuple](https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple) as the argument ``\kappa``: - -```jldoctest -julia> using Oceananigans.TurbulenceClosures - -julia> ScalarDiffusivity(ν=1e-6, κ=(S=1e-7, T=1e-10)) -ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-6, κ=(S=1.0e-7, T=1.0e-10)) -``` - -The example above sets a viscosity of `1e-6`, a diffusivity for a tracer called `T` of `1e-7`, -and a diffusivity for a tracer called `S` of `1e-10`. Specifying diffusivities this way is also valid -for `HorizontalScalarDiffusivity` and `VerticalScalarDiffusivity`. If this method is used, diffusivities -for all tracers need to be specified. - - -## Smagorinsky-Lilly - -To use the default Smagorinsky-Lilly LES closure, we write - -```jldoctest -julia> using Oceananigans.TurbulenceClosures - -julia> closure = SmagorinskyLilly() -Smagorinsky closure with -├── coefficient = LillyCoefficient(smagorinsky = 0.16, reduction_factor = 1.0) -└── Pr = 1.0 -``` - -The parameters `C`, `Cb`, and `Pr` may alternatively be specified explicitly. -For more details see [`SmagorinskyLilly`](@ref). - -## Anisotropic minimum dissipation - -To use the constant anisotropic minimum dissipation (AMD) LES closure, - -```jldoctest -julia> using Oceananigans.TurbulenceClosures - -julia> closure = AnisotropicMinimumDissipation() -AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with: - Poincaré constant for momentum eddy viscosity Cν: 0.3333333333333333 - Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.3333333333333333 - Buoyancy modification multiplier Cb: nothing -``` - -no parameters are required although they may be specified. By default, the background viscosity and diffusivity -are assumed to be the molecular values for seawater. For more details see [`AnisotropicMinimumDissipation`](@ref). - -## Convective Adjustment Vertical Diffusivity--Viscosity - -To use the a convective adjustment scheme that applies enhanced values for vertical diffusivity ``\kappa_z`` and/or -viscosity ``\nu_z``, anytime and anywhere the background stratification becomes unstable. - -```jldoctest -julia> using Oceananigans - -julia> closure = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0, background_κz = 1e-3) -ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.001 convective_κz=1.0 background_νz=0.0 convective_νz=0.0) -``` diff --git a/docs/src/simulations/callbacks.md b/docs/src/simulations/callbacks.md deleted file mode 100644 index 7b677e185c..0000000000 --- a/docs/src/simulations/callbacks.md +++ /dev/null @@ -1,93 +0,0 @@ -# [Callbacks](@id callbacks) - -A [`Callback`](@ref) can be used to execute an arbitrary user-defined function on the -simulation at user-defined times. - -For example, we can specify a callback which displays the run time every 2 iterations: -```@meta -DocTestSetup = quote - using Oceananigans -end -``` - -```@example checkpointing -using Oceananigans - -model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) - -simulation = Simulation(model, Δt=1, stop_iteration=10) - -show_time(sim) = @info "Time is $(prettytime(sim.model.clock.time))" - -simulation.callbacks[:total_A] = Callback(show_time, IterationInterval(2)) - -simulation -``` - -Now when simulation runs the simulation the callback is called. - -```@example checkpointing -run!(simulation) -``` - -We can also use the convenience [`add_callback!`](@ref): - -```@example checkpointing -add_callback!(simulation, show_time, name=:total_A_via_convenience, IterationInterval(2)) - -simulation -``` - -The keyword argument `callsite` determines the moment at which the callback is executed. -By default, `callsite = TimeStepCallsite()`, indicating execution _after_ the completion of -a timestep. The other options are `callsite = TendencyCallsite()` that executes the callback -after the tendencies are computed but _before_ taking a timestep and `callsite = UpdateStateCallsite()` -that executes the callback within `update_state!`, after auxiliary variables have been computed -(for multi-stage time-steppers, `update_state!` may be called multiple times per timestep). - -As an example of a callback with `callsite = TendencyCallsite()` , we show below how we can -manually add to the tendency field of one of the velocity components. Here we've chosen -the `:u` field using parameters: - -```@example checkpointing -using Oceananigans -using Oceananigans: TendencyCallsite - -model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) - -simulation = Simulation(model, Δt=1, stop_iteration=10) - -function modify_tendency!(model, params) - model.timestepper.Gⁿ[params.c] .+= params.δ - return nothing -end - -simulation.callbacks[:modify_u] = Callback(modify_tendency!, IterationInterval(1), - callsite = TendencyCallsite(), - parameters = (c = :u, δ = 1)) - -run!(simulation) -``` - -Above there is no forcing at all, but due to the callback the ``u``-velocity is increased. - -```@example checkpointing -@info model.velocities.u -``` - -!!! note "Example only for illustration purposes" - The above is a redundant example since it could be implemented better with a simple forcing function. - We include it here though for illustration purposes of how one can use callbacks. - -## Functions - -Callback functions can only take one or two parameters `sim` - a simulation, or `model` for state callbacks, and optionally may also accept a NamedTuple of parameters. - -## Scheduling - -Callbacks rely on the schedules described in [Scheduling callbacks and output writers](@ref callback_schedules). -The most common choices are: - - [`IterationInterval`](@ref): runs every `n` iterations - - [`TimeInterval`](@ref): runs every `n` seconds of model time - - [`SpecifiedTimes`](@ref): runs at the specified times - - [`WallTimeInterval`](@ref): runs every `n` seconds of wall time diff --git a/docs/src/simulations/checkpointing.md b/docs/src/simulations/checkpointing.md deleted file mode 100644 index af4ce81ddd..0000000000 --- a/docs/src/simulations/checkpointing.md +++ /dev/null @@ -1,61 +0,0 @@ -# [Checkpointing](@id checkpointing) - -A [`Checkpointer`](@ref) can be used to serialize the entire model state to a file from which the model -can be restored at any time. This is useful if you'd like to periodically checkpoint when running -long simulations in case of crashes or hitting cluster time limits, but also if you'd like to restore -from a checkpoint and try out multiple scenarios. - -For example, to periodically checkpoint the model state to disk every 1,000,000 seconds of simulation -time to files of the form `model_checkpoint_iteration12500.jld2` where `12500` is the iteration -number (automatically filled in). - -Here's an example where we checkpoint every 5 iterations. This is far more often than appropriate for -typical applications: we only do it here for illustration purposes. - -```@repl checkpointing -using Oceananigans - -model = NonhydrostaticModel(grid=RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1))) - -simulation = Simulation(model, Δt=1, stop_iteration=8) - -simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(5), prefix="model_checkpoint") -``` - -Again, for illustration purposes of this example, we also add another callback so we can see the iteration -of the simulation - -```@repl checkpointing -show_iteration(sim) = @info "iteration: $(iteration(sim)); time: $(prettytime(sim.model.clock.time))" -add_callback!(simulation, show_iteration, name=:info, IterationInterval(1)) -``` - -Now let's run - -```@repl checkpointing -run!(simulation) -``` - -The default options should provide checkpoint files that are easy to restore from (in most cases). -For more advanced options and features, see [`Checkpointer`](@ref). - -## Picking up a simulation from a checkpoint file - -Picking up a simulation from a checkpoint requires the original script that was used to generate -the checkpoint data. Change the first instance of [`run!`](@ref) in the script to take `pickup=true`. - -When `pickup=true` is provided to `run!` then it finds the latest checkpoint file in the current working -directory, loads prognostic fields and their tendencies from file, resets the model clock and iteration, -to the clock time and iteration that the checkpoint corresponds to, and updates the model auxiliary state. -After that, the time-stepping loop. In this simple example, although the simulation run up to iteration 8, -the latest checkpoint is associated with iteration 5. - -```@repl checkpointing -simulation.stop_iteration = 12 - -run!(simulation, pickup=true) -``` - -Use `pickup=iteration`, where `iteration` is an `Integer`, to pick up from a specific iteration. -Or, use `pickup=filepath`, where `filepath` is a string, to pickup from a specific file located -at `filepath`. diff --git a/docs/src/simulations/output_writers.md b/docs/src/simulations/output_writers.md deleted file mode 100644 index a9c5b8ac80..0000000000 --- a/docs/src/simulations/output_writers.md +++ /dev/null @@ -1,246 +0,0 @@ -# [Output writers](@id output_writers) - -`AbstractOutputWriter`s save data to disk. -`Oceananigans` provides three ways to write output: - -1. [`NetCDFWriter`](@ref) for output of arrays and scalars that uses [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) -2. [`JLD2Writer`](@ref) for arbitrary julia data structures that uses [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) -3. [`Checkpointer`](@ref) that automatically saves as much model data as possible, using [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) - -The `Checkpointer` is discussed in detail on a separate [section](@ref checkpointing) of the documentation. - -## Basic usage - -[`NetCDFWriter`](@ref) and [`JLD2Writer`](@ref) require four inputs: - -1. The `model` from which output data is sourced (required to initialize the `OutputWriter`). -2. A key-value pairing of output "names" and "output" objects. `JLD2Writer` accepts `NamedTuple`s and `Dict`s; - `NetCDFWriter` accepts `Dict`s with string-valued keys. Output objects are either `AbstractField`s or - functions that return data when called via `func(model)`. -3. A `schedule` on which output is written. `TimeInterval`, `IterationInterval`, `WallTimeInterval` schedule - periodic output according to the simulation time, simulation interval, or "wall time" (the physical time - according to a clock on your wall). A fourth `schedule` called `AveragedTimeInterval` specifies - periodic output that is time-averaged over a `window` prior to being written. See [callback schedules](@ref callback_schedules) - for a deeper overview of schedule types and how to combine them. -4. The `filename` and `dir`ectory. - -Other important keyword arguments are - -* `indices` for outputting subregions, two- and one-dimensional slices of fields. Specifies the indices to write to disk with a `Tuple` of `Colon`, `UnitRange`,or `Int` elements. For example, `indices = (:, :, 1)` implies outputing ``x-y``-slices of the bottom-most index (`k=1`). Defaults to `(:, :, :)`, i.e., "all indices". -* `with_halos :: Boolean`: whether to output the halos (`true`) or only the interior points (`false`; default). - -* `array_type` for specifying the type of the array that holds outputted field data. The default is - `Array{Float64}`, or arrays of single-precision floating point numbers. - -Once an `OutputWriter` is created, it can be used to write output by adding it the -ordered dictionary `simulation.output_writers`. prior to calling `run!(simulation)`. - -More specific detail about the `NetCDFWriter` and `JLD2Writer` is given below. - -!!! tip "Time step alignment and output writing" - Oceananigans simulations will shorten the time step as needed to align model output with each - output writer's schedule. - -## NetCDF output writer - -Model data can be saved to NetCDF files along with associated metadata. The NetCDF output writer is generally used by -passing it a dictionary of (label, field) pairs and any indices for slicing if you don't want to save the full 3D field. - -### Examples - -Saving the u velocity field and temperature fields as full 3D fields, surface 2D slices, and -1D columns to separate NetCDF files: - -```@example netcdf1 -using Oceananigans -using NCDatasets - -grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1)) - -model = NonhydrostaticModel(grid=grid, tracers=:c) - -simulation = Simulation(model, Δt=12, stop_time=3600) - -fields = Dict("u" => model.velocities.u, "c" => model.tracers.c) - -simulation.output_writers[:field_writer] = - NetCDFWriter(model, fields, filename="more_fields.nc", schedule=TimeInterval(60)) -``` - -```@example netcdf1 -simulation.output_writers[:surface_slice_writer] = - NetCDFWriter(model, fields, filename="another_surface_xy_slice.nc", - schedule=TimeInterval(60), indices=(:, :, grid.Nz)) -``` - -```@example netcdf1 -simulation.output_writers[:averaged_profile_writer] = - NetCDFWriter(model, fields, - filename = "another_averaged_z_profile.nc", - schedule = AveragedTimeInterval(60, window=20), - indices = (1, 1, :)) -``` - -`NetCDFWriter` also accepts output functions that write scalars and arrays to disk, -provided that their `dimensions` are provided: - -```@example -using Oceananigans -using NCDatasets - -Nx, Ny, Nz = 16, 16, 16 - -grid = RectilinearGrid(size=(Nx, Ny, Nz), extent=(1, 2, 3)) - -model = NonhydrostaticModel(grid=grid) - -simulation = Simulation(model, Δt=1.25, stop_iteration=3) - -f(model) = model.clock.time^2; # scalar output - -g(model) = model.clock.time .* exp.(znodes(grid, Center())) # single-column profile output (vector) - -xC, yF = xnodes(grid, Center()), ynodes(grid, Face()) - -XC = [xC[i] for i in 1:Nx, j in 1:Ny] -YF = [yF[j] for i in 1:Nx, j in 1:Ny] - -h(model) = @. model.clock.time * sin(XC) * cos(YF) # x-y slice output (2D array) - -outputs = Dict("scalar" => f, "profile" => g, "slice" => h) - -dims = Dict("scalar" => (), "profile" => ("z_aac",), "slice" => ("x_caa", "y_aca")) - -output_attributes = Dict( - "scalar" => Dict("longname" => "Some scalar", "units" => "bananas"), - "profile" => Dict("longname" => "Some vertical profile", "units" => "watermelons"), - "slice" => Dict("longname" => "Some slice", "units" => "mushrooms")) - -global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7) - -simulation.output_writers[:things] = - NetCDFWriter(model, outputs, - schedule=IterationInterval(1), filename="things.nc", dimensions=dims, verbose=true, - global_attributes=global_attributes, output_attributes=output_attributes) -``` - -`NetCDFWriter` can also be configured for `outputs` that are interpolated or regridded -to a different grid than `model.grid`. To use this functionality, include the keyword argument -`grid = output_grid`. - -```@example -using Oceananigans -using Oceananigans.Fields: interpolate! -using NCDatasets - -grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1)); -model = NonhydrostaticModel(; grid) - -coarse_grid = RectilinearGrid(size=(grid.Nx, grid.Ny, grid.Nz÷2), extent=(grid.Lx, grid.Ly, grid.Lz)) -coarse_u = Field{Face, Center, Center}(coarse_grid) - -interpolate_u(model) = interpolate!(coarse_u, model.velocities.u) -outputs = (; u = interpolate_u) - -output_writer = NetCDFWriter(model, outputs; - grid = coarse_grid, - filename = "coarse_u.nc", - schedule = IterationInterval(1)) -``` - -See [`NetCDFWriter`](@ref) for more information. - -## JLD2 output writer - -JLD2 is a fast HDF5 compatible file format written in pure Julia. -JLD2 files can be opened in Julia with the [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) package -and in Python with the [h5py](https://www.h5py.org/) package. - -The `JLD2Writer` receives either a `Dict`ionary or `NamedTuple` containing -`name, output` pairs. The `name` can be a symbol or string. The `output` must either be -an `AbstractField` or a function called with `func(model)` that returns arbitrary output. -Whenever output needs to be written, the functions will be called and the output -of the function will be saved to the JLD2 file. - -### Examples - -Write out 3D fields for u, v, w, and a tracer c, along with a horizontal average: - -```@example jld2_output_writer -using Oceananigans -using Oceananigans.Utils: hour, minute - -model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)), tracers=(:c,)) -simulation = Simulation(model, Δt=12, stop_time=1hour) - -function init_save_some_metadata!(file, model) - file["author"] = "Chim Riggles" - file["parameters/coriolis_parameter"] = 1e-4 - file["parameters/density"] = 1027 - return nothing -end - -c_avg = Field(Average(model.tracers.c, dims=(1, 2))) - -# Note that model.velocities is NamedTuple -simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities, - filename = "some_more_data.jld2", - schedule = TimeInterval(20minute), - init = init_save_some_metadata!) -``` - -and a time- and horizontal-average of tracer `c` every 20 minutes of simulation time -to a file called `some_more_averaged_data.jld2` - -```@example jld2_output_writer -simulation.output_writers[:avg_c] = JLD2Writer(model, (; c=c_avg), - filename = "some_more_averaged_data.jld2", - schedule = AveragedTimeInterval(20minute, window=5minute)) -``` - -See [`JLD2Writer`](@ref) for more information. - -## Time-averaged output - -Time-averaged output is specified by setting the `schedule` keyword argument for either `NetCDFWriter` or -`JLD2Writer` to [`AveragedTimeInterval`](@ref). - -With `AveragedTimeInterval`, the time-average of ``a`` is taken as a left Riemann sum corresponding to - -```math -\langle a \rangle = \frac{1}{T} \int_{t_i-T}^{t_i} a \, \mathrm{d} t \, , -``` - -where ``\langle a \rangle`` is the time-average of ``a``, ``T`` is the time-`window` for averaging specified by -the `window` keyword argument to `AveragedTimeInterval`, and the ``t_i`` are discrete times separated by the -time `interval`. The ``t_i`` specify both the end of the averaging window and the time at which output is written. - -### Example - -Building an `AveragedTimeInterval` that averages over a 1 day window, every 4 days, - -```jldoctest averaged_time_interval -using Oceananigans -using Oceananigans.Units - -schedule = AveragedTimeInterval(4days, window=1day) - -# output -AveragedTimeInterval(window=1 day, stride=1, interval=4 days) -``` - -An `AveragedTimeInterval` schedule directs an output writer -to time-average its outputs before writing them to disk: - -```@example averaged_time_interval -using Oceananigans -using Oceananigans.Units - -model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) - -simulation = Simulation(model, Δt=10minutes, stop_time=30days) - -simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities, - filename = "even_more_averaged_velocity_data.jld2", - schedule = AveragedTimeInterval(4days, window=1day, stride=2)) -``` From 307f8b5c6e08c4e9d67493a66569693e966a9084 Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 4 Nov 2025 09:48:49 -0500 Subject: [PATCH 15/16] Revert "md files removed" This reverts commit c2132c37daf93f487585d73ef585b41d3bd067c3. --- docs/src/models/background_fields.md | 107 +++ docs/src/models/boundary_conditions.md | 620 ++++++++++++++++++ .../models/buoyancy_and_equation_of_state.md | 251 +++++++ docs/src/models/coriolis.md | 112 ++++ docs/src/models/forcing_functions.md | 366 +++++++++++ docs/src/models/lagrangian_particles.md | 136 ++++ docs/src/models/models_overview.md | 208 ++++++ docs/src/models/turbulence_closures.md | 103 +++ docs/src/simulations/callbacks.md | 93 +++ docs/src/simulations/checkpointing.md | 61 ++ docs/src/simulations/output_writers.md | 246 +++++++ 11 files changed, 2303 insertions(+) create mode 100644 docs/src/models/background_fields.md create mode 100644 docs/src/models/boundary_conditions.md create mode 100644 docs/src/models/buoyancy_and_equation_of_state.md create mode 100644 docs/src/models/coriolis.md create mode 100644 docs/src/models/forcing_functions.md create mode 100644 docs/src/models/lagrangian_particles.md create mode 100644 docs/src/models/models_overview.md create mode 100644 docs/src/models/turbulence_closures.md create mode 100644 docs/src/simulations/callbacks.md create mode 100644 docs/src/simulations/checkpointing.md create mode 100644 docs/src/simulations/output_writers.md diff --git a/docs/src/models/background_fields.md b/docs/src/models/background_fields.md new file mode 100644 index 0000000000..ae55119b0b --- /dev/null +++ b/docs/src/models/background_fields.md @@ -0,0 +1,107 @@ +# Background fields + +`BackgroundField`s are velocity and tracer fields around which the resolved +velocity and tracer fields evolve. Only the _advective_ terms associated with +the interaction between background and resolved fields are included. +For example, tracer advection is described by + +```math +\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) \, , +``` + +where ``\boldsymbol{v}`` is the resolved velocity field and ``c`` is the resolved +tracer field corresponding to `model.tracers.c`. + +When a background field ``C`` is provided, the tracer advection term becomes + +```math +\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) + + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} C \right ) \, . +``` + +When both a background field velocity field ``\boldsymbol{U}`` and a background tracer field ``C`` +are provided, then the tracer advection term becomes + +```math +\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) + + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} C \right ) + + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{U} c \right ) \, . +``` + +Notice that the term ``\boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{U} C \right )`` +is neglected: only the terms describing the advection of resolved tracer by the background +velocity field and the advection of background tracer by the resolved velocity field are included. +An analogous statement holds for the advection of background momentum by the resolved +velocity field. +Other possible terms associated with the Coriolis force, buoyancy, turbulence closures, +and surface waves acting on background fields are neglected. + +!!! compat "Model compatibility" + `BackgroundFields` are only supported by [`NonhydrostaticModel`](@ref). + +## Specifying background fields + +`BackgroundField`s are defined by functions of ``(x, y, z, t)`` and optional parameters. A +simple example is + +```jldoctest +using Oceananigans + +U(x, y, z, t) = 0.2 * z + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) + +model = NonhydrostaticModel(grid = grid, background_fields = (u=U,)) + +model.background_fields.velocities.u + +# output +FunctionField located at (Face, Center, Center) +├── func: U (generic function with 1 method) +├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo +├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days) +└── parameters: nothing +``` + +`BackgroundField`s are specified by passing them to the kwarg `background_fields` +in the `NonhydrostaticModel` constructor. The kwarg `background_fields` expects +a `NamedTuple` of fields, which are internally sorted into `velocities` and `tracers`, +wrapped in `FunctionField`s, and assigned their appropriate locations. + +`BackgroundField`s with parameters require using the `BackgroundField` wrapper: + +```jldoctest moar_background +using Oceananigans + +parameters = (α=3.14, N=1.0, f=0.1) + +# Background fields are defined via function of x, y, z, t, and optional parameters +U(x, y, z, t, α) = α * z +B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z + +U_field = BackgroundField(U, parameters=parameters.α) +B_field = BackgroundField(B, parameters=parameters) + +# output +BackgroundField{typeof(B), @NamedTuple{α::Float64, N::Float64, f::Float64}} +├── func: B (generic function with 1 method) +└── parameters: (α = 3.14, N = 1.0, f = 0.1) +``` + +When inserted into `NonhydrostaticModel`, we get + +```jldoctest moar_background +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) + +model = NonhydrostaticModel(grid = grid, background_fields = (u=U_field, b=B_field), + tracers=:b, buoyancy=BuoyancyTracer()) + +model.background_fields.tracers.b + +# output +FunctionField located at (Center, Center, Center) +├── func: B (generic function with 1 method) +├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×1 halo +├── clock: Clock{Float64, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days) +└── parameters: (α = 3.14, N = 1.0, f = 0.1) +``` diff --git a/docs/src/models/boundary_conditions.md b/docs/src/models/boundary_conditions.md new file mode 100644 index 0000000000..2f03f15e69 --- /dev/null +++ b/docs/src/models/boundary_conditions.md @@ -0,0 +1,620 @@ +# [Boundary conditions](@id model_step_bcs) + +Boundary conditions are intimately related to the grid topology, and only +need to be considered in directions with `Bounded` topology or across immersed boundaries. +In `Bounded` directions, tracer and momentum fluxes are conservative or "zero flux" +by default. Non-default boundary conditions are therefore required to specify non-zero fluxes +of tracers and momentum across `Bounded` directions, and across immersed boundaries +when using `ImmersedBoundaryGrid`. + +See [Numerical implementation of boundary conditions](@ref numerical_bcs) for more details. + +### Example: no-slip conditions on every boundary + +```@meta +DocTestSetup = quote + using Oceananigans + + using Random + Random.seed!(1234) +end +``` + +```jldoctest bcsintro +julia> using Oceananigans + +julia> grid = RectilinearGrid(size=(16, 16, 16), x=(0, 2π), y=(0, 1), z=(0, 1), topology=(Periodic, Bounded, Bounded)) +16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo +├── Periodic x ∈ [0.0, 6.28319) regularly spaced with Δx=0.392699 +├── Bounded y ∈ [0.0, 1.0] regularly spaced with Δy=0.0625 +└── Bounded z ∈ [0.0, 1.0] regularly spaced with Δz=0.0625 + +julia> no_slip_bc = ValueBoundaryCondition(0.0) +ValueBoundaryCondition: 0.0 +``` + +A "no-slip" [`BoundaryCondition`](@ref) specifies that velocity components tangential to `Bounded` +directions decay to `0` at the boundary, leading to a viscous loss of momentum. + +```jldoctest bcsintro +julia> no_slip_field_bcs = FieldBoundaryConditions(no_slip_bc); + +julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=no_slip_field_bcs, v=no_slip_field_bcs, w=no_slip_field_bcs)) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +└── coriolis: Nothing + +julia> model.velocities.u.boundary_conditions +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: PeriodicBoundaryCondition +├── east: PeriodicBoundaryCondition +├── south: ValueBoundaryCondition: 0.0 +├── north: ValueBoundaryCondition: 0.0 +├── bottom: ValueBoundaryCondition: 0.0 +├── top: ValueBoundaryCondition: 0.0 +└── immersed: Nothing +``` + +Boundary conditions are passed to `FieldBoundaryCondition` to build boundary conditions for each +field individually, and then onto the model constructor (here `NonhydrotaticModel`) via the +keyword argument `boundary_conditions`. +The model constructor then "interprets" the input and builds appropriate boundary conditions +for the grid `topology`, given the user-specified `no_slip` default boundary condition for `Bounded` +directions. In the above example, note that the `west` and `east` boundary conditions are `PeriodicBoundaryCondition` +because the `x`-topology of the grid is `Periodic`. + +### Example: specifying boundary conditions on individual boundaries + +To specify no-slip boundary conditions on every `Bounded` direction _except_ +the surface, we write + +```jldoctest bcsintro +julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)); + +julia> model = NonhydrostaticModel(; grid, boundary_conditions=(u=free_slip_surface_bcs, v=free_slip_surface_bcs, w=no_slip_field_bcs)); + +julia> model.velocities.u.boundary_conditions +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: PeriodicBoundaryCondition +├── east: PeriodicBoundaryCondition +├── south: ValueBoundaryCondition: 0.0 +├── north: ValueBoundaryCondition: 0.0 +├── bottom: ValueBoundaryCondition: 0.0 +├── top: FluxBoundaryCondition: Nothing +└── immersed: Nothing + +julia> model.velocities.v.boundary_conditions +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: PeriodicBoundaryCondition +├── east: PeriodicBoundaryCondition +├── south: OpenBoundaryCondition{Nothing}: Nothing +├── north: OpenBoundaryCondition{Nothing}: Nothing +├── bottom: ValueBoundaryCondition: 0.0 +├── top: FluxBoundaryCondition: Nothing +└── immersed: Nothing +``` + +Now both `u` and `v` have `FluxBoundaryCondition(nothing)` at the `top` boundary, which is `Oceananigans` lingo +for "no-flux boundary condition". + +## Boundary condition classifications + +There are three primary boundary condition classifications: + +1. `FluxBoundaryCondition` specifies fluxes directly. + + Some applications of `FluxBoundaryCondition` are: + * surface momentum fluxes due to wind, or "wind stress"; + * linear or quadratic bottom drag; + * surface temperature fluxes due to heating or cooling; + * surface salinity fluxes due to precipitation and evaporation; + * relaxation boundary conditions that restores a field to some boundary distribution + over a given time-scale. + +2. `ValueBoundaryCondition` (Dirichlet) specifies the value of a field on + the given boundary, which when used in combination with a turbulence closure + results in a flux across the boundary. + + _Note_: Do not use `ValueBoundaryCondition` on a wall-normal velocity component + (see the note below about `ImpenetrableBoundaryCondition`). + + Some applications of `ValueBoundaryCondition` are: + * no-slip boundary condition for wall-tangential velocity components via `ValueBoundaryCondition(0)`; + * surface temperature distribution, where heat fluxes in and out of the domain + at a rate controlled by the near-surface temperature gradient and the temperature diffusivity; + * constant velocity tangential to a boundary as in a driven-cavity flow (for example), + where the top boundary is moving. Momentum will flux into the domain do the difference + between the top boundary velocity and the interior velocity, and the prescribed viscosity. + +3. `GradientBoundaryCondition` (Neumann) specifies the gradient of a field on a boundary. + For example, if there is a known `diffusivity`, we can express `FluxBoundaryCondition(flux)` + using `GradientBoundaryCondition(-flux / diffusivity)` (aka "Neumann" boundary condition). + +In addition to these primary boundary conditions, `ImpenetrableBoundaryCondition` applies to velocity +components in wall-normal directions. + +!!! warn "`ImpenetrableBoundaryCondition`" + `ImpenetrableBoundaryCondition` is internally enforced for fields created inside the model constructor. + As a result, `ImpenetrableBoundaryCondition` is only used for _additional_ velocity components + that are not evolved by a model, such as a velocity component used for (`AdvectiveForcing`)[@ref]. + +Finally, note that `Periodic` boundary conditions are internally enforced for `Periodic` directions, +and `DefaultBoundaryCondition`s may exist before boundary conditions are "materialized" by a model. + +## Default boundary conditions + +The default boundary condition in `Bounded` directions is no-flux, or `FluxBoundaryCondition(nothing)`. +The default boundary condition can be changed by passing a positional argument to `FieldBoundaryConditions`, +as in + +```jldoctest +julia> no_slip_bc = ValueBoundaryCondition(0.0) +ValueBoundaryCondition: 0.0 + +julia> free_slip_surface_bcs = FieldBoundaryConditions(no_slip_bc, top=FluxBoundaryCondition(nothing)) +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +├── east: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +├── south: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +├── north: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +├── bottom: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +├── top: FluxBoundaryCondition: Nothing +└── immersed: DefaultBoundaryCondition (ValueBoundaryCondition: 0.0) +``` + +## Boundary condition structures + +Oceananigans uses a hierarchical structure to express boundary conditions: + +1. Each boundary of each field has one [`BoundaryCondition`](@ref) +2. Each field has seven [`BoundaryCondition`](@ref) (`west`, `east`, `south`, `north`, `bottom`, `top` and + `immersed`) +3. A set of `FieldBoundaryConditions`, up to one for each field, are grouped into a `NamedTuple` and passed + to the model constructor. + +## Specifying boundary conditions for a model + +Boundary conditions are defined at model construction time by passing a `NamedTuple` of `FieldBoundaryConditions` +specifying non-default boundary conditions for fields such as velocities and tracers. + +Fields for which boundary conditions are not specified are assigned a default boundary conditions. + +A few illustrations are provided below. See the examples for +further illustrations of boundary condition specification. + +## Creating individual boundary conditions with `BoundaryCondition` + +Boundary conditions may be specified with constants, functions, or arrays. +In this section we illustrate usage of the different [`BoundaryCondition`](@ref) constructors. + +### 1. Constant `Value` (Dirchlet) boundary condition + +```jldoctest bcs +julia> constant_T_bc = ValueBoundaryCondition(20.0) +ValueBoundaryCondition: 20.0 +``` + +A constant [`Value`](@ref) boundary condition can be used to specify constant tracer (such as temperature), +or a constant _tangential_ velocity component at a boundary. Note that boundary conditions on the +_normal_ velocity component must use the [`Open`](@ref) boundary condition type. + +Finally, note that `ValueBoundaryCondition(condition)` is an alias for `BoundaryCondition(Value, condition)`. + +### 2. Constant `Flux` boundary condition + +```jldoctest +julia> ρ₀ = 1027; # Reference density [kg/m³] + +julia> τₓ = 0.08; # Wind stress [N/m²] + +julia> wind_stress_bc = FluxBoundaryCondition(-τₓ/ρ₀) +FluxBoundaryCondition: -7.78968e-5 +``` + +A constant [`Flux`](@ref) boundary condition can be imposed on tracers and tangential velocity components +that can be used, for example, to specify cooling, heating, evaporation, or wind stress at the ocean surface. + +!!! info "The flux convention in Oceananigans" + `Oceananigans` uses the convention that positive fluxes produce transport in the + _positive_ direction (east, north, and up for ``x``, ``y``, ``z``). + This means, for example, that a _negative_ flux of momentum or velocity at a _top_ + boundary, such as in the above example, produces currents in the _positive_ direction, + because it prescribes a downwards flux of momentum into the domain from the top. + Likewise, a _positive_ temperature flux at the top boundary + causes _cooling_, because it transports heat _upwards_, out of the domain. + Conversely, a positive flux at a _bottom_ boundary acts to increase the interior + values of a quantity. + +### 3. Spatially- and temporally-varying flux + +Boundary conditions may be specified by functions, + +```jldoctest +julia> @inline surface_flux(x, y, t) = cos(2π * x) * cos(t); + +julia> top_tracer_bc = FluxBoundaryCondition(surface_flux) +FluxBoundaryCondition: ContinuousBoundaryFunction surface_flux at (Nothing, Nothing, Nothing) +``` + +!!! info "Boundary condition functions" + By default, a function boundary condition is called with the signature + ```julia + f(ξ, η, t) + ``` + where `t` is time and `ξ, η` are spatial coordinates that vary along the boundary: + * `f(y, z, t)` on `x`-boundaries; + * `f(x, z, t)` on `y`-boundaries; + * `f(x, y, t)` on `z`-boundaries. + Alternative function signatures are specified by keyword arguments to + `BoundaryCondition`, as illustrated in subsequent examples. + +### 4. Spatially- and temporally-varying flux with parameters + +Boundary condition functions may be 'parameterized', + +```jldoctest +julia> @inline wind_stress(x, y, t, p) = - p.τ * cos(p.k * x) * cos(p.ω * t); # function with parameters + +julia> top_u_bc = FluxBoundaryCondition(wind_stress, parameters=(k=4π, ω=3.0, τ=1e-4)) +FluxBoundaryCondition: ContinuousBoundaryFunction wind_stress at (Nothing, Nothing, Nothing) +``` + +!!! info "Boundary condition functions with parameters" + The keyword argument `parameters` above specifies that `wind_stress` is called + with the signature `wind_stress(x, y, t, parameters)`. In principle, `parameters` is arbitrary. + However, relatively simple objects such as floating point numbers or `NamedTuple`s must be used + when running on the GPU. + +### 5. 'Field-dependent' boundary conditions + +Boundary conditions may also depend on model fields. For example, a linear drag boundary condition +is implemented with + +```jldoctest +julia> @inline linear_drag(x, y, t, u) = - 0.2 * u +linear_drag (generic function with 1 method) + +julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, field_dependencies=:u) +FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) +``` + +`field_dependencies` specifies the name of the dependent fields either with a `Symbol` or `Tuple` of `Symbol`s. + +### 6. 'Field-dependent' boundary conditions with parameters + +When boundary conditions depends on fields _and_ parameters, their functions take the form + +```jldoctest +julia> @inline quadratic_drag(x, y, t, u, v, drag_coeff) = - drag_coeff * u * sqrt(u^2 + v^2) +quadratic_drag (generic function with 1 method) + +julia> u_bottom_bc = FluxBoundaryCondition(quadratic_drag, field_dependencies=(:u, :v), parameters=1e-3) +FluxBoundaryCondition: ContinuousBoundaryFunction quadratic_drag at (Nothing, Nothing, Nothing) +``` + +Put differently, `ξ, η, t` come first in the function signature, followed by field dependencies, +followed by `parameters` is `!isnothing(parameters)`. + +### 7. Discrete-form boundary condition with parameters + +Discrete field data may also be accessed directly from boundary condition functions +using the `discrete_form`. For example: + +```jldoctest +@inline filtered_drag(i, j, grid, clock, model_fields) = + @inbounds - 0.05 * (model_fields.u[i-1, j, 1] + 2 * model_fields.u[i, j, 1] + model_fields.u[i-1, j, 1]) + +u_bottom_bc = FluxBoundaryCondition(filtered_drag, discrete_form=true) + +# output +FluxBoundaryCondition: DiscreteBoundaryFunction with filtered_drag +``` + +!!! info "The 'discrete form' for boundary condition functions" + The argument `discrete_form=true` indicates to [`BoundaryCondition`](@ref) that `filtered_drag` + uses the 'discrete form'. Boundary condition functions that use the 'discrete form' + are called with the signature + ```julia + f(i, j, grid, clock, model_fields) + ``` + where `i, j` are grid indices that vary along the boundary, `grid` is `model.grid`, + `clock` is the `model.clock`, and `model_fields` is a `NamedTuple` + containing `u, v, w` and the fields in `model.tracers`. + The signature is similar for ``x`` and ``y`` boundary conditions expect that `i, j` is replaced + with `j, k` and `i, k` respectively. + +### 8. Discrete-form boundary condition with parameters + +```jldoctest +julia> Cd = 0.2; # drag coefficient + +julia> @inline linear_drag(i, j, grid, clock, model_fields, Cd) = @inbounds - Cd * model_fields.u[i, j, 1]; + +julia> u_bottom_bc = FluxBoundaryCondition(linear_drag, discrete_form=true, parameters=Cd) +FluxBoundaryCondition: DiscreteBoundaryFunction linear_drag with parameters 0.2 +``` + +!!! info "Inlining and avoiding bounds-checking in boundary condition functions" + Boundary condition functions should be decorated with `@inline` when running on CPUs for performance reasons. + On the GPU, all functions are force-inlined by default. + In addition, the annotation `@inbounds` should be used when accessing the elements of an array + in a boundary condition function (such as `model_fields.u[i, j, 1]` in the above example). + Using `@inbounds` will avoid a relatively expensive check that the index `i, j, 1` is 'in bounds'. + +### 9. A random, spatially-varying, constant-in-time temperature flux specified by an array + +```jldoctest +julia> Nx = Ny = 16; # Number of grid points. + +julia> Q = randn(Nx, Ny); # temperature flux + +julia> white_noise_T_bc = FluxBoundaryCondition(Q) +FluxBoundaryCondition: 16×16 Matrix{Float64} +``` + +When running on the GPU, `Q` must be converted to a `CuArray`. + +## Building boundary conditions on a field + +To create a set of [`FieldBoundaryConditions`](@ref) for a temperature field, +we write + +```jldoctest +julia> T_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0), + bottom = GradientBoundaryCondition(0.01)) +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── bottom: GradientBoundaryCondition: 0.01 +├── top: ValueBoundaryCondition: 20.0 +└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +``` + +If the grid is, e.g., horizontally-periodic, then each horizontal `DefaultBoundaryCondition` +is converted to `PeriodicBoundaryCondition` inside the model's constructor, before assigning the +boundary conditions to temperature `T`. + +In general, boundary condition defaults are inferred from the field location and `topology(grid)`. + +## Specifying model boundary conditions + +To specify non-default boundary conditions, a named tuple of [`FieldBoundaryConditions`](@ref) objects is +passed to the keyword argument `boundary_conditions` in the [`NonhydrostaticModel`](@ref) constructor. +The keys of `boundary_conditions` indicate the field to which the boundary condition is applied. +Below, non-default boundary conditions are imposed on the ``u``-velocity and tracer ``c``. + +```jldoctest +julia> topology = (Periodic, Periodic, Bounded); + +julia> grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=topology); + +julia> u_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(+0.1), + bottom = ValueBoundaryCondition(-0.1)); + +julia> c_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(20.0), + bottom = GradientBoundaryCondition(0.01)); + +julia> model = NonhydrostaticModel(grid=grid, boundary_conditions=(u=u_bcs, c=c_bcs), tracers=:c) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: c +├── closure: Nothing +├── buoyancy: Nothing +└── coriolis: Nothing + +julia> model.velocities.u +16×16×16 Field{Face, Center, Center} on RectilinearGrid on CPU +├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── boundary conditions: FieldBoundaryConditions +│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Value, top: Value, immersed: Nothing +└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19 + └── max=0.0, min=0.0, mean=0.0 + +julia> model.tracers.c +16×16×16 Field{Center, Center, Center} on RectilinearGrid on CPU +├── grid: 16×16×16 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── boundary conditions: FieldBoundaryConditions +│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Gradient, top: Value, immersed: Nothing +└── data: 22×22×22 OffsetArray(::Array{Float64, 3}, -2:19, -2:19, -2:19) with eltype Float64 with indices -2:19×-2:19×-2:19 + └── max=0.0, min=0.0, mean=0.0 +``` + +Notice that the specified non-default boundary conditions have been applied at +top and bottom of both `model.velocities.u` and `model.tracers.c`. + +## Immersed boundary conditions + +Immersed boundary conditions are supported experimentally. A no-slip boundary condition is specified +with + +```@setup immersed_bc +using Oceananigans + +# Save original stderr +original_stderr = stderr + +# Redirect stderr to /dev/null (Unix) or "nul" (Windows) +redirect_stderr(devnull) +``` + +```jldoctest immersed_bc +# Generate a simple ImmersedBoundaryGrid +hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2) +underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded)) +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill)) + +# Create a no-slip boundary condition for velocity fields. +# Note that the no-slip boundary condition is _only_ applied on immersed boundaries. +velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0)) +model = NonhydrostaticModel(; grid, boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs)) + +# Inspect the boundary condition on the vertical velocity: +model.velocities.w.boundary_conditions.immersed + +# output +┌ Warning: The FFT-based pressure_solver for NonhydrostaticModels on ImmersedBoundaryGrid +│ is approximate and will probably produce velocity fields that are divergent +│ adjacent to the immersed boundary. An experimental but improved pressure_solver +│ is available which may be used by writing +│ +│ using Oceananigans.Solvers: ConjugateGradientPoissonSolver +│ pressure_solver = ConjugateGradientPoissonSolver(grid) +│ +│ Please report issues to https://github.com/CliMA/Oceananigans.jl/issues. +└ @ Oceananigans.Models.NonhydrostaticModels ~/Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:55 +ImmersedBoundaryCondition: +├── west: ValueBoundaryCondition: 0.0 +├── east: ValueBoundaryCondition: 0.0 +├── south: ValueBoundaryCondition: 0.0 +├── north: ValueBoundaryCondition: 0.0 +├── bottom: Nothing +└── top: Nothing +``` + +!!! warning "`NonhydrostaticModel` on `ImmersedBoundaryGrid`" + The pressure solver for `NonhydrostaticModel` is approximate, and is unable to produce + a velocity field that is simultaneously divergence-free while also satisfying impenetrability + on the immersed boundary. As a result, simulated dynamics with `NonhydrostaticModel` can + exhibit egregiously unphysical errors and should be interpreted with caution. + +```jldoctest immersed_bc +hill(x, y) = 0.1 + 0.1 * exp(-x^2 - y^2) +underlying_grid = RectilinearGrid(size=(32, 32, 16), x=(-3, 3), y=(-3, 3), z=(0, 1), topology=(Periodic, Periodic, Bounded)) +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hill)) + +# Create a no-slip boundary condition for velocity fields. +# Note that the no-slip boundary condition is _only_ applied on immersed boundaries. +velocity_bcs = FieldBoundaryConditions(immersed=ValueBoundaryCondition(0)) +model = NonhydrostaticModel(; grid, boundary_conditions=(u=velocity_bcs, v=velocity_bcs, w=velocity_bcs)) + +# output +┌ Warning: The FFT-based pressure_solver for NonhydrostaticModels on ImmersedBoundaryGrid +│ is approximate and will probably produce velocity fields that are divergent +│ adjacent to the immersed boundary. An experimental but improved pressure_solver +│ is available which may be used by writing +│ +│ using Oceananigans.Solvers: ConjugateGradientPoissonSolver +│ pressure_solver = ConjugateGradientPoissonSolver(grid) +│ +│ Please report issues to https://github.com/CliMA/Oceananigans.jl/issues. +└ @ Oceananigans.Models.NonhydrostaticModels ~/Oceananigans.jl/src/Models/NonhydrostaticModels/NonhydrostaticModels.jl:55 +NonhydrostaticModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0) +├── grid: 32×32×16 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +└── coriolis: Nothing +``` + +An `ImmersedBoundaryCondition` encapsulates boundary conditions on each potential boundary-facet +of a boundary-adjacent cell. Boundary conditions on specific faces of immersed-boundary-adjacent +cells may also be specified by manually building an `ImmersedBoundaryCondition`: + +```jldoctest immersed_bc +bottom_drag_bc = ImmersedBoundaryCondition(bottom=ValueBoundaryCondition(0)) + +# output +ImmersedBoundaryCondition: +├── west: Nothing +├── east: Nothing +├── south: Nothing +├── north: Nothing +├── bottom: ValueBoundaryCondition: 0 +└── top: Nothing +``` + +The `ImmersedBoundaryCondition` may then be incorporated into the boundary conditions for a +`Field` by prescribing it to the `immersed` boundary label, + +```jldoctest immersed_bc +velocity_bcs = FieldBoundaryConditions(immersed=bottom_drag_bc) + +# output +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Value, top=Nothing +``` + +!!! warning "`ImmersedBoundaryCondition`" + `ImmersedBoundaryCondition` is experimental. + Therefore, one should use it only when a finer level of control over the boundary conditions + at the immersed boundary is required, and the user is familiar with the implementation of boundary + conditions on staggered grids. For all other cases , using the `immersed` argument of + `FieldBoundaryConditions` is preferred. + +A boundary condition that depends on the fields may be prescribed using the `immersed` +keyword argument in [`FieldBoundaryConditions`](@ref). +We illustrate field-dependent boundary conditions with an example that imposes linear bottom drag +on `u` on both the bottom facets of cells adjacent to an immersed boundary, _and_ the bottom boundary +of the underlying grid. + +First we create the boundary condition for the grid's bottom: + +```jldoctest immersed_bc +@inline linear_drag(x, y, t, u) = - 0.2 * u +drag_u = FluxBoundaryCondition(linear_drag, field_dependencies=:u) + +# output +FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) +``` + +Next, we create the immersed boundary condition by adding the argument `z` to `linear_drag` +and imposing drag only on "bottom" facets of cells that neighbor immersed cells: + +```jldoctest immersed_bc +@inline immersed_linear_drag(x, y, z, t, u) = - 0.2 * u +immersed_drag_u = FluxBoundaryCondition(immersed_linear_drag, field_dependencies=:u) + +u_immersed_bc = ImmersedBoundaryCondition(bottom = immersed_drag_u) + +# output +ImmersedBoundaryCondition: +├── west: Nothing +├── east: Nothing +├── south: Nothing +├── north: Nothing +├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction immersed_linear_drag at (Nothing, Nothing, Nothing) +└── top: Nothing +``` + +Finally, we combine the two: + +```jldoctest immersed_bc +u_bcs = FieldBoundaryConditions(bottom = drag_u, immersed = u_immersed_bc) + +# output +Oceananigans.FieldBoundaryConditions, with boundary conditions +├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction linear_drag at (Nothing, Nothing, Nothing) +├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing) +└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Flux, top=Nothing +``` + +!!! warning "Positional argument requirements" + Note the difference between the arguments required for the function within the `bottom` boundary + condition versus the arguments for the function within the `immersed` boundary condition. E.g., + `x, y, t` in `linear_drag()` versus `x, y, z, t` in `immersed_linear_drag()`. + +```@setup immersed_bc +# Restore original stderr +redirect_stderr(original_stderr) +``` diff --git a/docs/src/models/buoyancy_and_equation_of_state.md b/docs/src/models/buoyancy_and_equation_of_state.md new file mode 100644 index 0000000000..a0dd6340d9 --- /dev/null +++ b/docs/src/models/buoyancy_and_equation_of_state.md @@ -0,0 +1,251 @@ +# Buoyancy models and equations of state + +The buoyancy option selects how buoyancy is treated in `NonhydrostaticModel`s and +`HydrostaticFreeSurfaceModel`s (`ShallowWaterModel`s do not have that option given the physics of +the model). There are currently three alternatives: + +1. No buoyancy (and no gravity). +2. Evolve buoyancy as a tracer. +3. _Seawater buoyancy_: evolve temperature ``T`` and salinity ``S`` as tracers with a value for the gravitational + acceleration ``g`` and an equation of state of your choosing. + +## No buoyancy + +To turn off buoyancy (and gravity) you can simply pass `buoyancy = nothing` to the model +constructor. For example to create a `NonhydrostaticModel`: + + +```@meta +DocTestSetup = quote + using Oceananigans +end +``` + + +```jldoctest buoyancy +julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)); + +julia> model = NonhydrostaticModel(; grid, buoyancy=nothing) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +└── coriolis: Nothing +``` + +The option `buoyancy = nothing` is the default for [`NonhydrostaticModel`](@ref), so omitting the +`buoyancy` keyword argument from the `NonhydrostaticModel` constructor yields the same: + +```jldoctest buoyancy +julia> model = NonhydrostaticModel(; grid) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +└── coriolis: Nothing +``` + +The same is true for `HydrostaticFreeSurfaceModel`, + +```jldoctest buoyancy +julia> model = HydrostaticFreeSurfaceModel(; grid) +HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: QuasiAdamsBashforth2TimeStepper +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² +│ └── solver: FFTImplicitFreeSurfaceSolver +├── advection scheme: +│ └── momentum: VectorInvariant +├── vertical_coordinate: ZCoordinate +└── coriolis: Nothing +``` + +## Buoyancy as a tracer + +Both `NonhydrostaticModel` and `HydrostaticFreeSurfaceModel` support evolving +a buoyancy tracer by including `:b` in `tracers` and specifying `buoyancy = BuoyancyTracer()`: + +```jldoctest buoyancy +julia> model = NonhydrostaticModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: b +├── closure: Nothing +├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() +└── coriolis: Nothing +``` + +Similarly for a `HydrostaticFreeSurfaceModel` with buoyancy as a tracer: + +```jldoctest buoyancy +julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=BuoyancyTracer(), tracers=:b) +HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: QuasiAdamsBashforth2TimeStepper +├── tracers: b +├── closure: Nothing +├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection() +├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² +│ └── solver: FFTImplicitFreeSurfaceSolver +├── advection scheme: +│ ├── momentum: VectorInvariant +│ └── b: Centered(order=2) +├── vertical_coordinate: ZCoordinate +└── coriolis: Nothing +``` + +## Seawater buoyancy + +`NonhydrostaticModel` and `HydrostaticFreeSurfaceModel` support modeling the buoyancy of seawater +as a function of the gravitational acceleration, the conservative temperature ``T``, and the absolute +salinity ``S``. The relationship between ``T``, ``S``, the geopotential height, and the density +perturbation from a reference value is called the `equation_of_state`. + +Specifying `buoyancy = SeawaterBuoyancy()` returns a buoyancy model with a linear equation of state, +[Earth standard](https://en.wikipedia.org/wiki/Standard_gravity) `gravitational_acceleration = 9.80665` (in +S.I. units ``\text{m}\,\text{s}^{-2}``) and requires to add `:T` and `:S` as tracers: + +```jldoctest buoyancy +julia> model = NonhydrostaticModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: (T, S) +├── closure: Nothing +├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() +└── coriolis: Nothing +``` + +and the same is true for `HydrostaticFreeSurfaceModel`, + +```jldoctest buoyancy +julia> model = HydrostaticFreeSurfaceModel(; grid, buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) +HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: QuasiAdamsBashforth2TimeStepper +├── tracers: (T, S) +├── closure: Nothing +├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() +├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻² +│ └── solver: FFTImplicitFreeSurfaceSolver +├── advection scheme: +│ ├── momentum: VectorInvariant +│ ├── T: Centered(order=2) +│ └── S: Centered(order=2) +├── vertical_coordinate: ZCoordinate +└── coriolis: Nothing +``` + +To model flows near the surface of Europa where `gravitational_acceleration = 1.3` ``\text{m}\,\text{s}^{-2}``, +we might alternatively specify + +```jldoctest buoyancy +julia> buoyancy = SeawaterBuoyancy(gravitational_acceleration=1.3) +SeawaterBuoyancy{Float64}: +├── gravitational_acceleration: 1.3 +└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) + +julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=(:T, :S)) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: (T, S) +├── closure: Nothing +├── buoyancy: SeawaterBuoyancy with g=1.3 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection() +└── coriolis: Nothing +``` + +for example. + +### Linear equation of state + +To specify the thermal expansion and haline contraction coefficients +``\alpha = 2 \times 10^{-3} \; \text{K}^{-1}`` and ``\beta = 5 \times 10^{-4} \text{psu}^{-1}``, + +```jldoctest +julia> buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-3, haline_contraction=5e-4)) +SeawaterBuoyancy{Float64}: +├── gravitational_acceleration: 9.80665 +└── equation_of_state: LinearEquationOfState(thermal_expansion=0.002, haline_contraction=0.0005) +``` + +### Idealized nonlinear equations of state + +Instead of a linear equation of state, six idealized (second-order) nonlinear equations of state +as described by [Roquet15Idealized](@citet) may be used. These equations of state are provided +via the [SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package. + +```jldoctest buoyancy +julia> using SeawaterPolynomials.SecondOrderSeawaterPolynomials + +julia> eos = RoquetEquationOfState(:Freezing) +BoussinesqEquationOfState{Float64}: +├── seawater_polynomial: SecondOrderSeawaterPolynomial{Float64} +└── reference_density: 1024.6 + +julia> eos.seawater_polynomial # the density anomaly +ρ' = 0.7718 Sᴬ - 0.0491 Θ - 0.005027 Θ² - 2.5681e-5 Θ Z + 0.0 Sᴬ² + 0.0 Sᴬ Z + 0.0 Sᴬ Θ + +julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos) +SeawaterBuoyancy{Float64}: +├── gravitational_acceleration: 9.80665 +└── equation_of_state: BoussinesqEquationOfState{Float64} +``` + +### TEOS-10 equation of state + +A high-accuracy 55-term polynomial approximation to the TEOS-10 equation of state suitable for use in +Boussinesq models as described by [Roquet15TEOS](@citet) is implemented in the +[SeawaterPolynomials.jl](https://github.com/CliMA/SeawaterPolynomials.jl) package and may be used. + +```jldoctest buoyancy +julia> using SeawaterPolynomials.TEOS10 + +julia> eos = TEOS10EquationOfState() +BoussinesqEquationOfState{Float64}: +├── seawater_polynomial: TEOS10SeawaterPolynomial{Float64} +└── reference_density: 1020.0 +``` + +## The direction of gravitational acceleration + +To simulate gravitational accelerations that don't align with the vertical (`z`) coordinate, +we use `BuoyancyForce(formulation; gravity_unit_vector)`, wherein the buoyancy `formulation` +can be `BuoyancyTracer`, `SeawaterBuoyancy`, etc, in addition to the `gravity_unit_vector`. +For example, + +```jldoctest buoyancy +julia> grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)); + +julia> θ = 45; # degrees + +julia> g̃ = (0, sind(θ), cosd(θ)); + +julia> buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃) +BuoyancyForce: +├── formulation: BuoyancyTracer +└── gravity_unit_vector: (0.0, 0.707107, 0.707107) + +julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=:b) +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: b +├── closure: Nothing +├── buoyancy: BuoyancyTracer with ĝ = (0.0, 0.707107, 0.707107) +└── coriolis: Nothing +``` diff --git a/docs/src/models/coriolis.md b/docs/src/models/coriolis.md new file mode 100644 index 0000000000..013dddac75 --- /dev/null +++ b/docs/src/models/coriolis.md @@ -0,0 +1,112 @@ +# Coriolis + +The Coriolis option determines whether the fluid experiences the effect of the Coriolis force, or rotation. Currently +three options are available: no rotation, ``f``-plane, and ``\beta``-plane. + +!!! info "Coriolis vs. rotation" + If you are wondering why this option is called "Coriolis" it is because rotational effects could include the + Coriolis and centripetal forces, both of which arise in non-inertial reference frames. But here the model only + considers the Coriolis force. + +## No rotation + +By default there is no rotation. This can be made explicit by passing `coriolis = nothing` to a model constructor. + +## Traditional ``f``-plane + +To set up an ``f``-plane with, for example, Coriolis parameter ``f = 10^{-4} \text{s}^{-1}`` + +```@meta +DocTestSetup = quote + using Oceananigans +end +``` + +```jldoctest +julia> coriolis = FPlane(f=1e-4) +FPlane{Float64}(f=0.0001) +``` + +An ``f``-plane can also be specified at some latitude on a spherical planet with a planetary rotation rate. For example, +to specify an ``f``-plane at a latitude of ``\varphi = 45°\text{N}`` on Earth which has a rotation rate of +``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` + +```jldoctest +julia> coriolis = FPlane(rotation_rate=7.292115e-5, latitude=45) +FPlane{Float64}(f=0.000103126) +``` + +in which case the value of ``f`` is given by ``2\Omega\sin\varphi``. + +## Coriolis term for constant rotation in a Cartesian coordinate system + +One can use `ConstantCartesianCoriolis` to set up a Coriolis acceleration term where the Coriolis parameter +is constant and the rotation axis is arbitrary. For example, with +``\boldsymbol{f} = (0, f_y, f_z) = (0, 2, 1) \times 10^{-4} \text{s}^{-1}``, + +```jldoctest +julia> coriolis = ConstantCartesianCoriolis(fx=0, fy=2e-4, fz=1e-4) +ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04 +``` + +Or alternatively, the same result can be achieved by specifying the magnitude of the Coriolis +frequency `f` and the `rotation_axis`. So another way to get a Coriolis acceleration with the same +values is: + +```jldoctest +julia> rotation_axis = (0, 2e-4, 1e-4)./√(2e-4^2 + 1e-4^2) # rotation_axis has to be a unit vector +(0.0, 0.8944271909999159, 0.4472135954999579) + +julia> coriolis = ConstantCartesianCoriolis(f=√(2e-4^2+1e-4^2), rotation_axis=rotation_axis) +ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04 +``` + +An ``f``-plane with non-traditional Coriolis terms can also be specified at some latitude on a spherical planet +with a planetary rotation rate. For example, to specify an ``f``-plane at a latitude of ``\varphi = 45°\text{N}`` +on Earth which has a rotation rate of ``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` + +```jldoctest +julia> coriolis = ConstantCartesianCoriolis(rotation_rate=7.292115e-5, latitude=45) +ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 1.03e-04, fz = 1.03e-04 +``` + +in which case ``f_z = 2\Omega\sin\varphi`` and ``f_y = 2\Omega\cos\varphi``. + +## Traditional ``\beta``-plane + +To set up a ``\beta``-plane the background rotation rate ``f_0`` and the ``\beta`` parameter must be specified. For example, +a ``\beta``-plane with ``f_0 = 10^{-4} \text{s}^{-1}`` and ``\beta = 1.5 \times 10^{-11} \text{s}^{-1}\text{m}^{-1}`` can be +set up with + +```jldoctest +julia> coriolis = BetaPlane(f₀=1e-4, β=1.5e-11) +BetaPlane{Float64}(f₀=0.0001, β=1.5e-11) +``` + +Alternatively, a ``\beta``-plane can also be set up at some latitude on a spherical planet with a planetary rotation rate +and planetary radius. For example, to specify a ``\beta``-plane at a latitude of ``\varphi = 10^\circ{S}`` on Earth +which has a rotation rate of ``\Omega = 7.292115 \times 10^{-5} \text{s}^{-1}`` and a radius of ``R = 6,371 \text{km}`` + +```jldoctest +julia> coriolis = BetaPlane(rotation_rate=7.292115e-5, latitude=-10, radius=6371e3) +BetaPlane{Float64}(f₀=-2.53252e-5, β=2.25438e-11) +``` + +in which case ``f_0 = 2\Omega\sin\varphi`` and ``\beta = 2\Omega\cos\varphi / R``. + +## Non-traditional ``\beta``-plane + +A non-traditional ``\beta``-plane requires either 5 parameters (by default Earth's radius and +rotation rate are used): + +```jldoctest +julia> NonTraditionalBetaPlane(fz=1e-4, fy=2e-4, β=4e-11, γ=-8e-11) +NonTraditionalBetaPlane{Float64}(fz = 1.00e-04, fy = 2.00e-04, β = 4.00e-11, γ = -8.00e-11, R = 6.37e+06) +``` + +or the rotation rate, radius, and latitude: + +```jldoctest +julia> NonTraditionalBetaPlane(rotation_rate=5.31e-5, radius=252.1e3, latitude=10) +NonTraditionalBetaPlane{Float64}(fz = 1.84e-05, fy = 1.05e-04, β = 4.15e-10, γ = -1.46e-10, R = 2.52e+05) +``` diff --git a/docs/src/models/forcing_functions.md b/docs/src/models/forcing_functions.md new file mode 100644 index 0000000000..504f18539b --- /dev/null +++ b/docs/src/models/forcing_functions.md @@ -0,0 +1,366 @@ +# [Forcing functions](@id forcing_functions) + +"Forcings" are user-defined terms appended to right-hand side of +the momentum or tracer evolution equations. In `Oceananigans`, momentum +and tracer forcings are defined via julia functions. `Oceananigans` includes +an interface for implementing forcing functions that depend on spatial coordinates, +time, model velocity and tracer fields, and external parameters. + +```@meta +DocTestSetup = quote + using Oceananigans +end +``` + +Forcings are added to models by passing a `NamedTuple` of functions or forcing objects +to the `forcing` keyword argument in `NonhydrostaticModel`'s constructor. +By default, momentum and tracer forcing functions are assumed to be functions of +`x, y, z, t`. A basic example is + +```jldoctest +u_forcing(x, y, z, t) = exp(z) * cos(x) * sin(t) + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) +model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing,)) + +model.forcing.u + +# output +ContinuousForcing{Nothing} at (Face, Center, Center) +├── func: u_forcing (generic function with 1 method) +├── parameters: nothing +└── field dependencies: () +``` + +More general forcing functions are built via the `Forcing` constructor +described below. `Oceananigans` also provides two convenience types: + + * `Relaxation` for damping terms that restore a field to a + target distribution outside of a masked region of space. `Relaxation` can be + used to implement sponge layers near the boundaries of a domain. + * `AdvectiveForcing` for advecting individual quantities by a separate or + "slip" velocity relative to both the prognostic model velocity field and any + `BackgroundField` velocity field. + +## The `Forcing` constructor + +The `Forcing` constructor provides an interface for specifying forcing functions that + +1. Depend on external parameters; and +2. Depend on model fields at the `x, y, z` location that forcing is applied; and/or +3. Require access to discrete model data. + +### Forcing functions with external parameters + +Most forcings involve external, changeable parameters. +Here are two examples of `forcing_func`tions that depend on +_(i)_ a single scalar parameter `s`, and _(ii)_ a `NamedTuple` of parameters, `p`: + +```jldoctest parameterized_forcing +# Forcing that depends on a scalar parameter `s` +u_forcing_func(x, y, z, t, s) = s * z + +u_forcing = Forcing(u_forcing_func, parameters=0.1) + +# Forcing that depends on a `NamedTuple` of parameters `p` +T_forcing_func(x, y, z, t, p) = - p.μ * exp(z / p.λ) * cos(p.k * x) * sin(p.ω * t) + +T_forcing = Forcing(T_forcing_func, parameters=(μ=1, λ=0.5, k=2π, ω=4π)) + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) +model = NonhydrostaticModel(grid=grid, forcing=(u=u_forcing, T=T_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) + +model.forcing.T + +# output +ContinuousForcing{@NamedTuple{μ::Int64, λ::Float64, k::Float64, ω::Float64}} at (Center, Center, Center) +├── func: T_forcing_func (generic function with 1 method) +├── parameters: (μ = 1, λ = 0.5, k = 6.283185307179586, ω = 12.566370614359172) +└── field dependencies: () +``` + +```jldoctest parameterized_forcing +model.forcing.u + +# output +ContinuousForcing{Float64} at (Face, Center, Center) +├── func: u_forcing_func (generic function with 1 method) +├── parameters: 0.1 +└── field dependencies: () +``` + +In this example, the objects passed to the `parameters` keyword in the construction of +`u_forcing` and `T_forcing` -- a floating point number for `u_forcing`, and a `NamedTuple` +of parameters for `T_forcing` -- are passed on to `u_forcing_func` and `T_forcing_func` when +they are called during time-stepping. The object passed to `parameters` is in principle arbitrary. +However, if using the GPU, then `typeof(parameters)` may be restricted by the requirements +of GPU-compiliability. + +### Forcing functions that depend on model fields + +Forcing functions may depend on model fields (velocity, tracers or auxiliary fields) evaluated at the `x, y, z` where forcing is applied. +Here's a somewhat non-sensical example: + +```jldoctest field_dependent_forcing +# Forcing that depends on the velocity fields `u`, `v`, and `w` +w_forcing_func(x, y, z, t, u, v, w) = - (u^2 + v^2 + w^2) / 2 + +w_forcing = Forcing(w_forcing_func, field_dependencies=(:u, :v, :w)) + +# Forcing that depends on salinity `S` and a scalar parameter +S_forcing_func(x, y, z, t, S, μ) = - μ * S + +S_forcing = Forcing(S_forcing_func, parameters=0.01, field_dependencies=:S) + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) +model = NonhydrostaticModel(grid=grid, forcing=(w=w_forcing, S=S_forcing), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) + +model.forcing.w + +# output +ContinuousForcing{Nothing} at (Center, Center, Face) +├── func: w_forcing_func (generic function with 1 method) +├── parameters: nothing +└── field dependencies: (:u, :v, :w) +``` + +```jldoctest field_dependent_forcing +model.forcing.S + +# output +ContinuousForcing{Float64} at (Center, Center, Center) +├── func: S_forcing_func (generic function with 1 method) +├── parameters: 0.01 +└── field dependencies: (:S,) +``` + +The `field_dependencies` arguments follow `x, y, z, t` in the forcing `func`tion in +the order they are specified in `Forcing`. +If both `field_dependencies` and `parameters` are specified, then the `field_dependencies` +arguments follow `x, y, z, t`, and `parameters` follow `field_dependencies`. + +Model fields that arise in the arguments of continuous `Forcing` `func`tions are +automatically interpolated to the staggered grid location at which the forcing is applied. + +### "Discrete form" forcing functions + +"Discrete form" forcing functions are either called with the signature + +```julia +func(i, j, k, grid, clock, model_fields) +``` + +or the parameterized form + +```julia +func(i, j, k, grid, clock, model_fields, parameters) +``` + +Discrete form forcing functions can access the entirety of model field +data through the argument `model_fields`. The object `model_fields` is a `NamedTuple` +whose properties include the velocity fields `model_fields.u`, `model_fields.v`, +`model_fields.w` and all fields in `model.tracers`. + +Using discrete forcing functions may require understanding the +staggered arrangement of velocity fields and tracers in `Oceananigans`. +Here's a slightly non-sensical example in which the vertical derivative of a buoyancy +tracer is used as a time-scale for damping the u-velocity field: + +```jldoctest discrete_forcing +# A damping term that depends on a "local average": +local_average(i, j, k, grid, c) = @inbounds (c[i, j, k] + c[i-1, j, k] + c[i+1, j, k] + + c[i, j-1, k] + c[i, j+1, k] + + c[i, j, k-1] + c[i, j, k+1]) / 7 + +b_forcing_func(i, j, k, grid, clock, model_fields) = - local_average(i, j, k, grid, model_fields.b) + +b_forcing = Forcing(b_forcing_func, discrete_form=true) + +# A term that damps the local velocity field in the presence of stratification +using Oceananigans.Operators: ∂zᶠᶜᶠ, ℑxzᶠᵃᶜ + +function u_forcing_func(i, j, k, grid, clock, model_fields, ε) + # The vertical derivative of buoyancy, interpolated to the u-velocity location: + N² = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶠᶜᶠ, model_fields.b) + + # Set to zero in unstable stratification where N² < 0: + N² = max(N², zero(typeof(N²))) + + return @inbounds - ε * sqrt(N²) * model_fields.u[i, j, k] +end + +u_forcing = Forcing(u_forcing_func, discrete_form=true, parameters=1e-3) + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) +model = NonhydrostaticModel(grid=grid, tracers=:b, buoyancy=BuoyancyTracer(), forcing=(u=u_forcing, b=b_forcing)) + +model.forcing.b + +# output +DiscreteForcing{Nothing} +├── func: b_forcing_func (generic function with 1 method) +└── parameters: nothing +``` + +```jldoctest discrete_forcing +model.forcing.u + +# output +DiscreteForcing{Float64} +├── func: u_forcing_func (generic function with 1 method) +└── parameters: 0.001 +``` + +The annotation `@inbounds` is crucial for performance when accessing array indices +of the fields in `model_fields`. + +## `Relaxation` + +`Relaxation` defines a special forcing function that restores a field at a specified `rate` to +a `target` distribution, within a region uncovered by a `mask`ing function. +`Relaxation` is useful for implementing sponge layers, as shown in the second example. + +The following code constructs a model in which all components +of the velocity field are damped to zero everywhere on a time-scale of 1000 seconds, or ~17 minutes: + +```jldoctest +damping = Relaxation(rate = 1/1000) + +grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)) +model = NonhydrostaticModel(grid=grid, forcing=(u=damping, v=damping, w=damping)) + +model.forcing.w + +# output +ContinuousForcing{Nothing} at (Center, Center, Face) +├── func: Relaxation(rate=0.001, mask=1, target=0) +├── parameters: nothing +└── field dependencies: (:w,) +``` + +The constructor for `Relaxation` accepts the keyword arguments `mask`, and `target`, +which specify a `mask(x, y, z)` function that multiplies the forcing, and a `target(x, y, z)` +distribution for the quantity in question. By default, `mask` uncovered the whole domain +and `target` restores the field in question to 0 + +We illustrate usage of `mask` and `target` by implementing a sponge layer that relaxes +velocity fields to zero and restores temperature to a linear gradient in the bottom +1/10th of the domain: + +```jldoctest sponge_layer +grid = RectilinearGrid(size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(-1, 0)) + + damping_rate = 1/100 # relax fields on a 100 second time-scale +temperature_gradient = 0.001 # ⁰C m⁻¹ + surface_temperature = 20 # ⁰C (at z=0) + +target_temperature = LinearTarget{:z}(intercept=surface_temperature, gradient=temperature_gradient) + bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10) + +uvw_sponge = Relaxation(rate=damping_rate, mask=bottom_mask) + T_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=target_temperature) + +model = NonhydrostaticModel(grid=grid, forcing=(u=uvw_sponge, v=uvw_sponge, w=uvw_sponge, T=T_sponge), buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)) + +model.forcing.u + +# output +ContinuousForcing{Nothing} at (Face, Center, Center) +├── func: Relaxation(rate=0.01, mask=exp(-(z + 1.0)^2 / (2 * 0.1^2)), target=0) +├── parameters: nothing +└── field dependencies: (:u,) +``` + +```jldoctest sponge_layer +model.forcing.T + +# output +ContinuousForcing{Nothing} at (Center, Center, Center) +├── func: Relaxation(rate=0.01, mask=exp(-(z + 1.0)^2 / (2 * 0.1^2)), target=20.0 + 0.001 * z) +├── parameters: nothing +└── field dependencies: (:T,) +``` + +## `AdvectiveForcing` + +`AdvectiveForcing` defines a forcing function that represents advection by +a separate or "slip" velocity relative to the prognostic model velocity field. +`AdvectiveForcing` is implemented with native Oceananigans advection operators, +which means that tracers advected by the "flux form" advection term +``𝛁⋅𝐮_{\rm slip} c``. Caution is advised when ``𝐮_{\rm slip}`` is not divergence free. + +As an example, consider a model for sediment settling at a constant rate: + +```jldoctest +using Oceananigans + +r_sediment = 1e-4 # [m] "Fine sand" +ρ_sediment = 1200 # kg m⁻³ +ρ_ocean = 1026 # kg m⁻³ +Δb = 9.81 * (ρ_ocean - ρ_sediment) / ρ_ocean # m s⁻² +ν_molecular = 1.05e-6 # m² s⁻¹ +w_sediment = 2/9 * Δb / ν_molecular * r_sediment^2 # m s⁻¹ + +sinking = AdvectiveForcing(w=w_sediment) + +# output +AdvectiveForcing: +├── u: ZeroField{Int64} +├── v: ZeroField{Int64} +└── w: ConstantField(-0.00352102) +``` + +The three keyword arguments specify the `u`, `v`, and `w` components of the separate +slip velocity field. The default for each `u, v, w` is `ZeroField`. + +Next we consider a dynamically-evolving slip velocity. For this we use `ZFaceField` +with appropriate boundary conditions as our slip velocity: + +```jldoctest sinking +using Oceananigans +using Oceananigans.BoundaryConditions: ImpenetrableBoundaryCondition + +grid = RectilinearGrid(size=(32, 32, 32), x=(-10, 10), y=(-10, 10), z=(-4, 4), + topology=(Periodic, Periodic, Bounded)) + +no_penetration = ImpenetrableBoundaryCondition() +slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face), + top=no_penetration, bottom=no_penetration) + +w_slip = ZFaceField(grid, boundary_conditions=slip_bcs) +sinking = AdvectiveForcing(w=w_slip) + +# output +AdvectiveForcing: +├── u: ZeroField{Int64} +├── v: ZeroField{Int64} +└── w: 32×32×33 Field{Center, Center, Face} on RectilinearGrid on CPU +``` + +To compute the slip velocity, we must add a `Callback`to `simulations.callback` that +computes `w_slip` ever iteration: + +```jldoctest sinking +using Oceananigans.BoundaryConditions: fill_halo_regions! + +model = NonhydrostaticModel(; grid, tracers=(:b, :P), forcing=(; P=sinking)) +simulation = Simulation(model; Δt=1, stop_iteration=100) + +# Build abstract operation for slip velocity +b_particle = - 1e-4 # relative buoyancy depends on reference density and initial buoyancy condition +b = model.tracers.b +R = 1e-3 # [m] mean particle radius +ν = 1.05e-6 # [m² s⁻¹] molecular kinematic viscosity of water +w_slip_op = 2/9 * (b - b_particle) / ν * R^2 # Stokes terminal velocity + +function compute_slip_velocity!(sim) + w_slip .= w_slip_op + fill_halo_regions!(w_slip) + return nothing +end + +simulation.callbacks[:slip] = Callback(compute_slip_velocity!) + +# output +Callback of compute_slip_velocity! on IterationInterval(1) +``` diff --git a/docs/src/models/lagrangian_particles.md b/docs/src/models/lagrangian_particles.md new file mode 100644 index 0000000000..e2b99571c6 --- /dev/null +++ b/docs/src/models/lagrangian_particles.md @@ -0,0 +1,136 @@ +# Lagrangian particle tracking + +Models can keep track of the location and properties of neutrally buoyant particles. Particles are +advected with the flow field using forward Euler time-stepping at every model iteration. + +## Simple particles + +If you just need to keep of particle locations ``(x, y, z)`` then you can construct some Lagrangian particles +using the regular `LagrangianParticles` constructor + +```@meta +DocTestSetup = quote + using Oceananigans +end +``` + +```jldoctest particles +grid = RectilinearGrid(size=(10, 10, 10), extent=(1, 1, 1)) + +Nparticles = 10 + +x₀ = zeros(Nparticles) + +y₀ = rand(Nparticles) + +z₀ = -0.5 * ones(Nparticles) + +lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀) + +# output +10 LagrangianParticles with eltype Particle: +├── 3 properties: (:x, :y, :z) +├── particle-wall restitution coefficient: 1.0 +├── 0 tracked fields: () +└── dynamics: no_dynamics +``` + +then pass it to a model constructor + +```jldoctest particles +model = NonhydrostaticModel(grid=grid, particles=lagrangian_particles) + +# output +NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) +├── grid: 10×10×10 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── timestepper: RungeKutta3TimeStepper +├── advection scheme: Centered(order=2) +├── tracers: () +├── closure: Nothing +├── buoyancy: Nothing +├── coriolis: Nothing +└── particles: 10 LagrangianParticles with eltype Particle and properties (:x, :y, :z) +``` + +!!! warn "Lagrangian particles on GPUs" + Remember to use `CuArray` instead of regular `Array` when storing particle locations and properties on the GPU. + +## Custom particles + +If you want to keep track of custom properties, such as the species or DNA of a Lagrangian particle +representing a microbe in an agent-based model, then you can create your own custom particle type +and pass a `StructArray` to the `LagrangianParticles` constructor. + +```jldoctest particles +using Oceananigans +using StructArrays + +struct LagrangianMicrobe{T, S, D} + x :: T + y :: T + z :: T + species :: S + dna :: D +end + +Nparticles = 3 + +x₀ = zeros(Nparticles) + +y₀ = rand(Nparticles) + +z₀ = -0.5 * ones(Nparticles) + +species = [:rock, :paper, :scissors] + +dna = ["TATACCCC", "CCTAGGAC", "CGATTTAA"] + +particles = StructArray{LagrangianMicrobe}((x₀, y₀, z₀, species, dna)); + +lagrangian_particles = LagrangianParticles(particles) + +# output +3 LagrangianParticles with eltype LagrangianMicrobe: +├── 5 properties: (:x, :y, :z, :species, :dna) +├── particle-wall restitution coefficient: 1.0 +├── 0 tracked fields: () +└── dynamics: no_dynamics +``` + +!!! warn "Custom properties on GPUs" + Not all data types can be passed to GPU kernels. If you intend to advect particles on the GPU make sure + particle properties consist of only simple data types. The symbols and strings in this example won't + work on the GPU. + +## Writing particle properties to disk + +Particle properties can be written to disk using JLD2 or NetCDF. + +When writing to JLD2 you can pass `model.particles` as part of the named tuple of outputs. + +```@setup particles +using Oceananigans +using NCDatasets +grid = RectilinearGrid(size=(10, 10, 10), extent=(1, 1, 1)) +Nparticles = 3 +x₀ = zeros(Nparticles) +y₀ = rand(Nparticles) +z₀ = -0.5 * ones(Nparticles) +lagrangian_particles = LagrangianParticles(x=x₀, y=y₀, z=z₀) +model = NonhydrostaticModel(; grid, particles=lagrangian_particles) +``` + +```@example particles +JLD2Writer(model, (; model.particles), filename="particles", schedule=TimeInterval(15)) +``` + +When writing to NetCDF you should write particles to a separate file as the NetCDF dimensions differ for +particle trajectories. You can just pass `model.particles` straight to `NetCDFWriter`: + +```@example particles +NetCDFWriter(model, (; model.particles), filename="particles.nc", schedule=TimeInterval(15)) +``` + +!!! warn "Outputting custom particle properties to NetCDF" + NetCDF does not support arbitrary data types. If you need to write custom particle properties to disk + that are not supported by NetCDF then you should use JLD2 (which should support almost any Julia data type). diff --git a/docs/src/models/models_overview.md b/docs/src/models/models_overview.md new file mode 100644 index 0000000000..75ce544750 --- /dev/null +++ b/docs/src/models/models_overview.md @@ -0,0 +1,208 @@ +# Models + +In general in Oceananigans, the "model" object serves two main purposes: + * models store the configuration of a set of discrete equations. The discrete equations imply rules for evolving + prognostic variables, and computing diagnostic variables from the prognostic state. + * models provide a container for the prognostic and diagnostic state of those discrete equations at a particular time. + +## Two Oceananigans models for ocean simulations + +In addition to defining the abstract concept of a "model" that can be used with [Simulation](@ref), +Oceananigans provides two mature model implementations for simulating ocean-flavored fluid dynamics. +Both of these integrate the Navier-Stokes equations within the Boussinesq approximation +(we call these the "Boussinesq equations" for short): the [NonhydrostaticModel](@ref) and +the [HydrostaticFreeSurfaceModel](@ref). + +The [NonhydrostaticModel](@ref) integrates the Boussinesq equations _without_ making the hydrostatic approximation, +and therefore possessing a prognostic vertical momentum equation. The NonhydrostaticModel is useful for simulations +that resolve three-dimensional turbulence, such as large eddy simulations on [RectilinearGrid](@ref) with grid +spacings of O(1 m), as well as direct numerical simulation. The NonhydrostaticModel may also be used for +idealized classroom problems, as in the [two-dimensional turbulence example](@ref "Two dimensional turbulence example"). + +The [HydrostaticFreeSurfaceModel](@ref) integrates the hydrostatic or "primitive" Boussinesq equations +with a free surface on its top boundary. The hydrostatic approximation allosw the HydrostaticFreeSurfaceModel +to achieve much higher efficiency in simulations on curvilinear grids used for large-scale regional or global +simulations such as [LatitudeLongitudeGrid](@ref), [TripolarGrid](@ref), [ConformalCubedSphereGrid](@ref), +and other [OrthogonalSphericalShellGrid](@ref)s such as [RotatedLatitudeLongitudeGrid](@ref Oceananigans.OrthogonalSphericalShellGrids.RotatedLatitudeLongitudeGrid). +Because they span larger domains, simulations with the HydrostaticFreeSurfaceModel also usually involve coarser +grid spacings of O(30 m) up to O(100 km). Such coarse-grained simulations are usually paired with more elaborate +turbulence closures or "parameterizations" than small-scale simulations with NonhydrostaticModel, such as the +vertical mixing schemes [CATKEVerticalDiffusivity](@ref), [RiBasedVerticalDiffusivity](@ref), and +[TKEDissipationVerticalDiffusivity](@ref), and the mesoscale turbulence closure +[IsopycnalSkewSymmetricDiffusivity](@ref) (a.k.a. "Gent-McWilliams plus Redi"). + +A third, experimental [ShallowWaterModel](@ref) solves the shallow water equations. + +### Configuring NonhydrostaticModel with keyword arguments + +To illustrate the specification of discrete equations, consider first the docstring for [NonhydrostaticModel](@ref), + +```@docs +NonhydrostaticModel +``` + +The configuration operations for NonhydrostaticModel include "discretization choices", such as the `advection` scheme, +as well as aspects of the continuous underlying equations, such as the formulation of the buoyancy force. + +For our first example, we build the default `NonhydrostaticModel` (which is quite boring): + +```@example +using Oceananigans +grid = RectilinearGrid(size=(8, 8, 8), extent=(8, 8, 8)) +nh = NonhydrostaticModel(; grid) +``` + +The default `NonhydrostaticModel` has no tracers, no buoyancy force, no Coriolis force, and a second-order advection scheme. +We next consider a slightly more exciting NonhydrostaticModel configured with a WENO advection scheme, +the temperature/salinity-based [SeawaterBuoyancy](@ref), a boundary condition on the zonal momentum, +and a passive tracer forced by a cooked-up surface flux called "c": + +```@example first_model +using Oceananigans + +grid = RectilinearGrid(size=(128, 128), halo=(5, 5), x=(0, 256), z=(-128, 0), + topology = (Periodic, Flat, Bounded)) + +# Numerical method and physics choices +advection = WENO(order=9) # ninth‑order upwind for momentum and tracers +buoyancy = BuoyancyTracer() +coriolis = FPlane(f=1e-4) + +τx = - 8e-5 # m² s⁻² (τₓ < 0 ⟹ eastward wind stress) +u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(τx)) + +@inline Jc(x, t, Lx) = cos(2π / Lx * x) +c_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Jc, parameters=grid.Lx)) + +model = NonhydrostaticModel(; grid, advection, buoyancy, coriolis, + tracers = (:b, :c), + boundary_conditions = (; u=u_bcs, c=c_bcs)) +``` + +### Mutation of the model state + +In addition to providing an interface for configuring equations, models also +store the prognostic and diagnostic state associated with the solution to those equations. +Models thus also provide an interface for "setting" or fixing the prognostic state, which is typically +invoked to determine the initial conditions of a simulation. +To illustrate this we consider setting the above model to a stably-stratified and noisy condition: + +```@example first_model +N² = 1e-5 +bᵢ(x, z) = N² * z + 1e-6 * randn() +uᵢ(x, z) = 1e-3 * randn() +set!(model, b=bᵢ, u=uᵢ, w=uᵢ) + +model.tracers.b +``` + +Invoking `set!` above determine the model tracer `b` and the velocity components `u` and `w`. +`set!` also computes the diagnostic state of a model, which in the case of `NonhydrostaticModel` includes +the nonhydrostatic component of pressure, + +```@example first_model +model.pressures.pNHS +``` + +### Evolving models in time + +Model may be integrated or "stepped forward" in time by calling `time_step!(model, Δt)`, where `Δt` is the time step +and thus advancing the `model.clock`: + +```@example first_model +time_step!(model, 1) +model.clock +``` + +However, users are strongly encouraged to use the [`Simulation`](@ref) interface to manage time-stepping +along with other activities, like monitoring progress, writing output to disk, and more. + +```@example first_model +simulation = Simulation(model, Δt=1, stop_iteration=10) +run!(simulation) + +simulation +``` + +## Using the HydrostaticFreeSurfaceModel + +The HydrostaticFreeSurfaceModel has a similar interface as the NonhydrostaticModel, + +```@example +using Oceananigans +grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1)) +model = HydrostaticFreeSurfaceModel(; grid) # default free surface, no tracers +``` + +The full array of keyword arguments used to configure a HydrostaticFreeSurfaceModel are detailed +in the docstring for [HydrostaticFreeSurfaceModel](@ref), + +```@docs +HydrostaticFreeSurfaceModel +``` + +A bit more involved HydrostaticFreeSurfaceModel example: + +```@example second_model +using Oceananigans +using SeawaterPolynomials: TEOS10EquationOfState + +grid = LatitudeLongitudeGrid(size = (180, 80, 10), + longitude = (0, 360), + latitude = (-80, 80), + z = (-1000, 0), + halo = (6, 6, 3)) + +momentum_advection = WENOVectorInvariant() +coriolis = HydrostaticSphericalCoriolis() +equation_of_state = TEOS10EquationOfState() +buoyancy = SeawaterBuoyancy(; equation_of_state) +closure = CATKEVerticalDiffusivity() + +# Generate a zonal wind stress that mimics Earth's mean winds +# with westerlies in mid-latitudes and easterlies near equator and poles +function zonal_wind_stress(λ, φ, t) + # Parameters + τ₀ = 1e-4 # Maximum wind stress magnitude (N/m²) + φ₀ = 30 # Latitude of maximum westerlies (degrees) + dφ = 10 + + # Approximate wind stress pattern + return - τ₀ * (+ exp(-(φ - φ₀)^2 / 2dφ^2) + - exp(-(φ + φ₀)^2 / 2dφ^2) + - 0.3 * exp(-φ^2 / dφ^2)) +end + +u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(zonal_wind_stress)) + +model = HydrostaticFreeSurfaceModel(; grid, momentum_advection, coriolis, closure, buoyancy, + boundary_conditions = (; u=u_bcs), tracers=(:T, :S, :e)) +``` + +Mutating the state of the HydrostaticFreeSurfaceModel works similarly as for the NonhydrostaticModel --- +except that the vertical velocity cannot be `set!`, because vertical velocity is not +prognostic in the hydrostatic equations. + +```@example second_model +using SeawaterPolynomials + +N² = 1e-5 +T₀ = 20 +S₀ = 35 +eos = model.buoyancy.formulation.equation_of_state +α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos) +g = model.buoyancy.formulation.gravitational_acceleration +dTdz = N² / (α * g) +Tᵢ(λ, φ, z) = T₀ + dTdz * z + 1e-3 * T₀ * randn() +uᵢ(λ, φ, z) = 1e-3 * randn() +set!(model, T=Tᵢ, S=S₀, u=uᵢ, v=uᵢ) + +model.tracers.T +``` + +## Where to go next + +- See [Quick start](@ref quick_start) for a compact, end-to-end workflow +- See the Examples gallery for longer tutorials covering specific cases, including large eddy simulation, Kelvin–Helmholtz instability and baroclinic instability. +- Other pages in the Models section: in‑depth pages on buoyancy, forcing, boundary conditions, closures, diagnostics, and output. +- Physics: governing equations and numerical forms for [`NonhydrostaticModel`](@ref) and [`HydrostaticFreeSurfaceModel`](@ref hydrostatic_free_surface_model). diff --git a/docs/src/models/turbulence_closures.md b/docs/src/models/turbulence_closures.md new file mode 100644 index 0000000000..29c2605cfd --- /dev/null +++ b/docs/src/models/turbulence_closures.md @@ -0,0 +1,103 @@ +# [Turbulent diffusivity closures and large eddy simulation models](@id turbulence_closures) + +A turbulent diffusivity closure representing the effects of viscous dissipation and diffusion can +be passed via the `closure` keyword. + +See [turbulence closures](@ref numerical_closures) and [large eddy simulation](@ref numerical_les) for +more details on turbulent diffusivity closures. + +## Constant isotropic diffusivity + +To use constant isotropic values for the viscosity ``\nu`` and diffusivity ``\kappa`` you can use `ScalarDiffusivity`: + +```jldoctest +julia> using Oceananigans.TurbulenceClosures + +julia> closure = ScalarDiffusivity(ν=1e-2, κ=1e-2) +ScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.01, κ=0.01) +``` + +## Constant anisotropic diffusivity + +To specify constant values for the horizontal and vertical viscosities, ``\nu_h`` and ``\nu_z``, +and horizontal and vertical diffusivities, ``\kappa_h`` and ``\kappa_z``, you can use +`HorizontalScalarDiffusivity()` and `VerticalScalarDiffusivity()`, e.g., + +```jldoctest +julia> using Oceananigans.TurbulenceClosures + +julia> horizontal_closure = HorizontalScalarDiffusivity(ν=1e-3, κ=2e-3) +HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.001, κ=0.002) + +julia> vertical_closure = VerticalScalarDiffusivity(ν=1e-3, κ=2e-3) +VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.001, κ=0.002) +``` + +After that you can set, e.g., `closure = (horizontal_closure, vertical_closure)` when constructing +the model so that all components will be taken into account when calculating the diffusivity term. +Note that `VerticalScalarDiffusivity` and `HorizontalScalarDiffusivity` are implemented using [different +schemes](https://mitgcm.readthedocs.io/en/latest/algorithm/algorithm.html#horizontal-dissipation) +with different conservation properties. + +## Tracer-specific diffusivities + +You can also set different diffusivities for each tracer in your simulation by passing +a [NamedTuple](https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple) as the argument ``\kappa``: + +```jldoctest +julia> using Oceananigans.TurbulenceClosures + +julia> ScalarDiffusivity(ν=1e-6, κ=(S=1e-7, T=1e-10)) +ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-6, κ=(S=1.0e-7, T=1.0e-10)) +``` + +The example above sets a viscosity of `1e-6`, a diffusivity for a tracer called `T` of `1e-7`, +and a diffusivity for a tracer called `S` of `1e-10`. Specifying diffusivities this way is also valid +for `HorizontalScalarDiffusivity` and `VerticalScalarDiffusivity`. If this method is used, diffusivities +for all tracers need to be specified. + + +## Smagorinsky-Lilly + +To use the default Smagorinsky-Lilly LES closure, we write + +```jldoctest +julia> using Oceananigans.TurbulenceClosures + +julia> closure = SmagorinskyLilly() +Smagorinsky closure with +├── coefficient = LillyCoefficient(smagorinsky = 0.16, reduction_factor = 1.0) +└── Pr = 1.0 +``` + +The parameters `C`, `Cb`, and `Pr` may alternatively be specified explicitly. +For more details see [`SmagorinskyLilly`](@ref). + +## Anisotropic minimum dissipation + +To use the constant anisotropic minimum dissipation (AMD) LES closure, + +```jldoctest +julia> using Oceananigans.TurbulenceClosures + +julia> closure = AnisotropicMinimumDissipation() +AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with: + Poincaré constant for momentum eddy viscosity Cν: 0.3333333333333333 + Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.3333333333333333 + Buoyancy modification multiplier Cb: nothing +``` + +no parameters are required although they may be specified. By default, the background viscosity and diffusivity +are assumed to be the molecular values for seawater. For more details see [`AnisotropicMinimumDissipation`](@ref). + +## Convective Adjustment Vertical Diffusivity--Viscosity + +To use the a convective adjustment scheme that applies enhanced values for vertical diffusivity ``\kappa_z`` and/or +viscosity ``\nu_z``, anytime and anywhere the background stratification becomes unstable. + +```jldoctest +julia> using Oceananigans + +julia> closure = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0, background_κz = 1e-3) +ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.001 convective_κz=1.0 background_νz=0.0 convective_νz=0.0) +``` diff --git a/docs/src/simulations/callbacks.md b/docs/src/simulations/callbacks.md new file mode 100644 index 0000000000..7b677e185c --- /dev/null +++ b/docs/src/simulations/callbacks.md @@ -0,0 +1,93 @@ +# [Callbacks](@id callbacks) + +A [`Callback`](@ref) can be used to execute an arbitrary user-defined function on the +simulation at user-defined times. + +For example, we can specify a callback which displays the run time every 2 iterations: +```@meta +DocTestSetup = quote + using Oceananigans +end +``` + +```@example checkpointing +using Oceananigans + +model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=1, stop_iteration=10) + +show_time(sim) = @info "Time is $(prettytime(sim.model.clock.time))" + +simulation.callbacks[:total_A] = Callback(show_time, IterationInterval(2)) + +simulation +``` + +Now when simulation runs the simulation the callback is called. + +```@example checkpointing +run!(simulation) +``` + +We can also use the convenience [`add_callback!`](@ref): + +```@example checkpointing +add_callback!(simulation, show_time, name=:total_A_via_convenience, IterationInterval(2)) + +simulation +``` + +The keyword argument `callsite` determines the moment at which the callback is executed. +By default, `callsite = TimeStepCallsite()`, indicating execution _after_ the completion of +a timestep. The other options are `callsite = TendencyCallsite()` that executes the callback +after the tendencies are computed but _before_ taking a timestep and `callsite = UpdateStateCallsite()` +that executes the callback within `update_state!`, after auxiliary variables have been computed +(for multi-stage time-steppers, `update_state!` may be called multiple times per timestep). + +As an example of a callback with `callsite = TendencyCallsite()` , we show below how we can +manually add to the tendency field of one of the velocity components. Here we've chosen +the `:u` field using parameters: + +```@example checkpointing +using Oceananigans +using Oceananigans: TendencyCallsite + +model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=1, stop_iteration=10) + +function modify_tendency!(model, params) + model.timestepper.Gⁿ[params.c] .+= params.δ + return nothing +end + +simulation.callbacks[:modify_u] = Callback(modify_tendency!, IterationInterval(1), + callsite = TendencyCallsite(), + parameters = (c = :u, δ = 1)) + +run!(simulation) +``` + +Above there is no forcing at all, but due to the callback the ``u``-velocity is increased. + +```@example checkpointing +@info model.velocities.u +``` + +!!! note "Example only for illustration purposes" + The above is a redundant example since it could be implemented better with a simple forcing function. + We include it here though for illustration purposes of how one can use callbacks. + +## Functions + +Callback functions can only take one or two parameters `sim` - a simulation, or `model` for state callbacks, and optionally may also accept a NamedTuple of parameters. + +## Scheduling + +Callbacks rely on the schedules described in [Scheduling callbacks and output writers](@ref callback_schedules). +The most common choices are: + - [`IterationInterval`](@ref): runs every `n` iterations + - [`TimeInterval`](@ref): runs every `n` seconds of model time + - [`SpecifiedTimes`](@ref): runs at the specified times + - [`WallTimeInterval`](@ref): runs every `n` seconds of wall time diff --git a/docs/src/simulations/checkpointing.md b/docs/src/simulations/checkpointing.md new file mode 100644 index 0000000000..af4ce81ddd --- /dev/null +++ b/docs/src/simulations/checkpointing.md @@ -0,0 +1,61 @@ +# [Checkpointing](@id checkpointing) + +A [`Checkpointer`](@ref) can be used to serialize the entire model state to a file from which the model +can be restored at any time. This is useful if you'd like to periodically checkpoint when running +long simulations in case of crashes or hitting cluster time limits, but also if you'd like to restore +from a checkpoint and try out multiple scenarios. + +For example, to periodically checkpoint the model state to disk every 1,000,000 seconds of simulation +time to files of the form `model_checkpoint_iteration12500.jld2` where `12500` is the iteration +number (automatically filled in). + +Here's an example where we checkpoint every 5 iterations. This is far more often than appropriate for +typical applications: we only do it here for illustration purposes. + +```@repl checkpointing +using Oceananigans + +model = NonhydrostaticModel(grid=RectilinearGrid(size=(8, 8, 8), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=1, stop_iteration=8) + +simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(5), prefix="model_checkpoint") +``` + +Again, for illustration purposes of this example, we also add another callback so we can see the iteration +of the simulation + +```@repl checkpointing +show_iteration(sim) = @info "iteration: $(iteration(sim)); time: $(prettytime(sim.model.clock.time))" +add_callback!(simulation, show_iteration, name=:info, IterationInterval(1)) +``` + +Now let's run + +```@repl checkpointing +run!(simulation) +``` + +The default options should provide checkpoint files that are easy to restore from (in most cases). +For more advanced options and features, see [`Checkpointer`](@ref). + +## Picking up a simulation from a checkpoint file + +Picking up a simulation from a checkpoint requires the original script that was used to generate +the checkpoint data. Change the first instance of [`run!`](@ref) in the script to take `pickup=true`. + +When `pickup=true` is provided to `run!` then it finds the latest checkpoint file in the current working +directory, loads prognostic fields and their tendencies from file, resets the model clock and iteration, +to the clock time and iteration that the checkpoint corresponds to, and updates the model auxiliary state. +After that, the time-stepping loop. In this simple example, although the simulation run up to iteration 8, +the latest checkpoint is associated with iteration 5. + +```@repl checkpointing +simulation.stop_iteration = 12 + +run!(simulation, pickup=true) +``` + +Use `pickup=iteration`, where `iteration` is an `Integer`, to pick up from a specific iteration. +Or, use `pickup=filepath`, where `filepath` is a string, to pickup from a specific file located +at `filepath`. diff --git a/docs/src/simulations/output_writers.md b/docs/src/simulations/output_writers.md new file mode 100644 index 0000000000..a9c5b8ac80 --- /dev/null +++ b/docs/src/simulations/output_writers.md @@ -0,0 +1,246 @@ +# [Output writers](@id output_writers) + +`AbstractOutputWriter`s save data to disk. +`Oceananigans` provides three ways to write output: + +1. [`NetCDFWriter`](@ref) for output of arrays and scalars that uses [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) +2. [`JLD2Writer`](@ref) for arbitrary julia data structures that uses [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) +3. [`Checkpointer`](@ref) that automatically saves as much model data as possible, using [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) + +The `Checkpointer` is discussed in detail on a separate [section](@ref checkpointing) of the documentation. + +## Basic usage + +[`NetCDFWriter`](@ref) and [`JLD2Writer`](@ref) require four inputs: + +1. The `model` from which output data is sourced (required to initialize the `OutputWriter`). +2. A key-value pairing of output "names" and "output" objects. `JLD2Writer` accepts `NamedTuple`s and `Dict`s; + `NetCDFWriter` accepts `Dict`s with string-valued keys. Output objects are either `AbstractField`s or + functions that return data when called via `func(model)`. +3. A `schedule` on which output is written. `TimeInterval`, `IterationInterval`, `WallTimeInterval` schedule + periodic output according to the simulation time, simulation interval, or "wall time" (the physical time + according to a clock on your wall). A fourth `schedule` called `AveragedTimeInterval` specifies + periodic output that is time-averaged over a `window` prior to being written. See [callback schedules](@ref callback_schedules) + for a deeper overview of schedule types and how to combine them. +4. The `filename` and `dir`ectory. + +Other important keyword arguments are + +* `indices` for outputting subregions, two- and one-dimensional slices of fields. Specifies the indices to write to disk with a `Tuple` of `Colon`, `UnitRange`,or `Int` elements. For example, `indices = (:, :, 1)` implies outputing ``x-y``-slices of the bottom-most index (`k=1`). Defaults to `(:, :, :)`, i.e., "all indices". +* `with_halos :: Boolean`: whether to output the halos (`true`) or only the interior points (`false`; default). + +* `array_type` for specifying the type of the array that holds outputted field data. The default is + `Array{Float64}`, or arrays of single-precision floating point numbers. + +Once an `OutputWriter` is created, it can be used to write output by adding it the +ordered dictionary `simulation.output_writers`. prior to calling `run!(simulation)`. + +More specific detail about the `NetCDFWriter` and `JLD2Writer` is given below. + +!!! tip "Time step alignment and output writing" + Oceananigans simulations will shorten the time step as needed to align model output with each + output writer's schedule. + +## NetCDF output writer + +Model data can be saved to NetCDF files along with associated metadata. The NetCDF output writer is generally used by +passing it a dictionary of (label, field) pairs and any indices for slicing if you don't want to save the full 3D field. + +### Examples + +Saving the u velocity field and temperature fields as full 3D fields, surface 2D slices, and +1D columns to separate NetCDF files: + +```@example netcdf1 +using Oceananigans +using NCDatasets + +grid = RectilinearGrid(size=(16, 16, 16), extent=(1, 1, 1)) + +model = NonhydrostaticModel(grid=grid, tracers=:c) + +simulation = Simulation(model, Δt=12, stop_time=3600) + +fields = Dict("u" => model.velocities.u, "c" => model.tracers.c) + +simulation.output_writers[:field_writer] = + NetCDFWriter(model, fields, filename="more_fields.nc", schedule=TimeInterval(60)) +``` + +```@example netcdf1 +simulation.output_writers[:surface_slice_writer] = + NetCDFWriter(model, fields, filename="another_surface_xy_slice.nc", + schedule=TimeInterval(60), indices=(:, :, grid.Nz)) +``` + +```@example netcdf1 +simulation.output_writers[:averaged_profile_writer] = + NetCDFWriter(model, fields, + filename = "another_averaged_z_profile.nc", + schedule = AveragedTimeInterval(60, window=20), + indices = (1, 1, :)) +``` + +`NetCDFWriter` also accepts output functions that write scalars and arrays to disk, +provided that their `dimensions` are provided: + +```@example +using Oceananigans +using NCDatasets + +Nx, Ny, Nz = 16, 16, 16 + +grid = RectilinearGrid(size=(Nx, Ny, Nz), extent=(1, 2, 3)) + +model = NonhydrostaticModel(grid=grid) + +simulation = Simulation(model, Δt=1.25, stop_iteration=3) + +f(model) = model.clock.time^2; # scalar output + +g(model) = model.clock.time .* exp.(znodes(grid, Center())) # single-column profile output (vector) + +xC, yF = xnodes(grid, Center()), ynodes(grid, Face()) + +XC = [xC[i] for i in 1:Nx, j in 1:Ny] +YF = [yF[j] for i in 1:Nx, j in 1:Ny] + +h(model) = @. model.clock.time * sin(XC) * cos(YF) # x-y slice output (2D array) + +outputs = Dict("scalar" => f, "profile" => g, "slice" => h) + +dims = Dict("scalar" => (), "profile" => ("z_aac",), "slice" => ("x_caa", "y_aca")) + +output_attributes = Dict( + "scalar" => Dict("longname" => "Some scalar", "units" => "bananas"), + "profile" => Dict("longname" => "Some vertical profile", "units" => "watermelons"), + "slice" => Dict("longname" => "Some slice", "units" => "mushrooms")) + +global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7) + +simulation.output_writers[:things] = + NetCDFWriter(model, outputs, + schedule=IterationInterval(1), filename="things.nc", dimensions=dims, verbose=true, + global_attributes=global_attributes, output_attributes=output_attributes) +``` + +`NetCDFWriter` can also be configured for `outputs` that are interpolated or regridded +to a different grid than `model.grid`. To use this functionality, include the keyword argument +`grid = output_grid`. + +```@example +using Oceananigans +using Oceananigans.Fields: interpolate! +using NCDatasets + +grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 1)); +model = NonhydrostaticModel(; grid) + +coarse_grid = RectilinearGrid(size=(grid.Nx, grid.Ny, grid.Nz÷2), extent=(grid.Lx, grid.Ly, grid.Lz)) +coarse_u = Field{Face, Center, Center}(coarse_grid) + +interpolate_u(model) = interpolate!(coarse_u, model.velocities.u) +outputs = (; u = interpolate_u) + +output_writer = NetCDFWriter(model, outputs; + grid = coarse_grid, + filename = "coarse_u.nc", + schedule = IterationInterval(1)) +``` + +See [`NetCDFWriter`](@ref) for more information. + +## JLD2 output writer + +JLD2 is a fast HDF5 compatible file format written in pure Julia. +JLD2 files can be opened in Julia with the [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) package +and in Python with the [h5py](https://www.h5py.org/) package. + +The `JLD2Writer` receives either a `Dict`ionary or `NamedTuple` containing +`name, output` pairs. The `name` can be a symbol or string. The `output` must either be +an `AbstractField` or a function called with `func(model)` that returns arbitrary output. +Whenever output needs to be written, the functions will be called and the output +of the function will be saved to the JLD2 file. + +### Examples + +Write out 3D fields for u, v, w, and a tracer c, along with a horizontal average: + +```@example jld2_output_writer +using Oceananigans +using Oceananigans.Utils: hour, minute + +model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1)), tracers=(:c,)) +simulation = Simulation(model, Δt=12, stop_time=1hour) + +function init_save_some_metadata!(file, model) + file["author"] = "Chim Riggles" + file["parameters/coriolis_parameter"] = 1e-4 + file["parameters/density"] = 1027 + return nothing +end + +c_avg = Field(Average(model.tracers.c, dims=(1, 2))) + +# Note that model.velocities is NamedTuple +simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities, + filename = "some_more_data.jld2", + schedule = TimeInterval(20minute), + init = init_save_some_metadata!) +``` + +and a time- and horizontal-average of tracer `c` every 20 minutes of simulation time +to a file called `some_more_averaged_data.jld2` + +```@example jld2_output_writer +simulation.output_writers[:avg_c] = JLD2Writer(model, (; c=c_avg), + filename = "some_more_averaged_data.jld2", + schedule = AveragedTimeInterval(20minute, window=5minute)) +``` + +See [`JLD2Writer`](@ref) for more information. + +## Time-averaged output + +Time-averaged output is specified by setting the `schedule` keyword argument for either `NetCDFWriter` or +`JLD2Writer` to [`AveragedTimeInterval`](@ref). + +With `AveragedTimeInterval`, the time-average of ``a`` is taken as a left Riemann sum corresponding to + +```math +\langle a \rangle = \frac{1}{T} \int_{t_i-T}^{t_i} a \, \mathrm{d} t \, , +``` + +where ``\langle a \rangle`` is the time-average of ``a``, ``T`` is the time-`window` for averaging specified by +the `window` keyword argument to `AveragedTimeInterval`, and the ``t_i`` are discrete times separated by the +time `interval`. The ``t_i`` specify both the end of the averaging window and the time at which output is written. + +### Example + +Building an `AveragedTimeInterval` that averages over a 1 day window, every 4 days, + +```jldoctest averaged_time_interval +using Oceananigans +using Oceananigans.Units + +schedule = AveragedTimeInterval(4days, window=1day) + +# output +AveragedTimeInterval(window=1 day, stride=1, interval=4 days) +``` + +An `AveragedTimeInterval` schedule directs an output writer +to time-average its outputs before writing them to disk: + +```@example averaged_time_interval +using Oceananigans +using Oceananigans.Units + +model = NonhydrostaticModel(grid=RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=10minutes, stop_time=30days) + +simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities, + filename = "even_more_averaged_velocity_data.jld2", + schedule = AveragedTimeInterval(4days, window=1day, stride=2)) +``` From a71018d54ae01034e3dc7c1b3924e4b2f5dd9dca Mon Sep 17 00:00:00 2001 From: Xin Kai Lee Date: Tue, 4 Nov 2025 19:58:01 -0500 Subject: [PATCH 16/16] add compatiblity for AveragedSpecifiedTimes with datetime, remove checks for now --- src/OutputWriters/windowed_time_average.jl | 41 ++++++++++++---------- src/Utils/times_and_datetimes.jl | 8 +++-- 2 files changed, 29 insertions(+), 20 deletions(-) diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 7718c24ddb..bc5142ce46 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -2,7 +2,7 @@ using Oceananigans.Diagnostics: AbstractDiagnostic using Oceananigans.OutputWriters: fetch_output using Oceananigans.Utils: AbstractSchedule, prettytime using Oceananigans.TimeSteppers: Clock -using Dates: Period +using Dates: Period, Second, value import Oceananigans: run_diagnostic! import Oceananigans.Utils: TimeInterval, SpecifiedTimes @@ -119,28 +119,30 @@ AveragedSpecifiedTimes(times; window, kw...) = AveragedSpecifiedTimes(times, win determine_epsilon(eltype) = 0 determine_epsilon(::Type{T}) where T <: AbstractFloat = eps(T) -determine_epsilon(::Period) = Second(0) +determine_epsilon(::Type{<:Period}) = Second(0) function AveragedSpecifiedTimes(times, window::Vector; kw...) length(window) == length(times) || throw(ArgumentError("When providing a vector of windows, its length $(length(window)) must match the number of specified times $(length(times)).")) perm = sortperm(times) sorted_times = times[perm] sorted_window = window[perm] - time_diff = diff(vcat(0, sorted_times)) + # time_diff = diff(vcat(0, sorted_times)) + # time_diff = diff(sorted_times) - epsilon = determine_epsilon(eltype(window)) - any(time_diff .- sorted_window .< -epsilon) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) + # epsilon = determine_epsilon(eltype(window)) + # any(time_diff .- sorted_window .< -epsilon) && throw(ArgumentError("Averaging windows overlap. Ensure that for each specified time tᵢ, tᵢ - windowᵢ ≥ tᵢ₋₁.")) return AveragedSpecifiedTimes(SpecifiedTimes(sorted_times); window=sorted_window, kw...) end function AveragedSpecifiedTimes(times, window; kw...) sorted_times = sort(times) - time_diff = diff(vcat(0, sorted_times)) + # time_diff = diff(vcat(0, sorted_times)) + # time_diff = diff(sorted_times) - epsilon = determine_epsilon(typeof(window)) + # epsilon = determine_epsilon(typeof(window)) - any(time_diff .- window .< -epsilon) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) + # any(time_diff .- window .< -epsilon) && throw(ArgumentError("Averaging window $window is too large and causes overlapping windows. Ensure that for each specified time tᵢ, tᵢ - window ≥ tᵢ₋₁.")) return AveragedSpecifiedTimes(SpecifiedTimes(times); window, kw...) end @@ -167,7 +169,7 @@ function outside_window(schedule::AveragedSpecifiedTimes, clock) next > length(schedule.specified_times.times) && return true next_time = schedule.specified_times.times[next] window = get_next_window(schedule) - return clock.time < next_time - window + return clock.time <= next_time - window end function end_of_window(schedule::AveragedSpecifiedTimes, clock) @@ -186,18 +188,18 @@ next_actuation_time(sch::AveragedSpecifiedTimes) = Oceananigans.Utils.next_actua ##### WindowedTimeAverage ##### -mutable struct WindowedTimeAverage{OP, R, S} <: AbstractDiagnostic +mutable struct WindowedTimeAverage{OP, R, T, S} <: AbstractDiagnostic result :: R operand :: OP - window_start_time :: Float64 + window_start_time :: T window_start_iteration :: Int - previous_collection_time :: Float64 + previous_collection_time :: T schedule :: S fetch_operand :: Bool end -const IntervalWindowedTimeAverage = WindowedTimeAverage{<:Any, <:Any, <:AveragedTimeInterval} -const SpecifiedWindowedTimeAverage = WindowedTimeAverage{<:Any, <:Any, <:AveragedSpecifiedTimes} +const IntervalWindowedTimeAverage = WindowedTimeAverage{<:Any, <:Any, <:Any, <:AveragedTimeInterval} +const SpecifiedWindowedTimeAverage = WindowedTimeAverage{<:Any, <:Any, <:Any, <:AveragedSpecifiedTimes} stride(wta::IntervalWindowedTimeAverage) = wta.schedule.stride stride(wta::SpecifiedWindowedTimeAverage) = wta.schedule.stride @@ -224,7 +226,7 @@ function WindowedTimeAverage(operand, model=nothing; schedule, fetch_operand=tru result .= operand end - return WindowedTimeAverage(result, operand, 0.0, 0, 0.0, schedule, fetch_operand) + return WindowedTimeAverage(result, operand, model.clock.time, 0, model.clock.time, schedule, fetch_operand) end # Time-averaging doesn't change spatial location @@ -251,12 +253,15 @@ function accumulate_result!(wta, model) return accumulate_result!(wta, model.clock, integrand) end +period_to_number(p::Period) = value(p) +period_to_number(n::Number) = n + function accumulate_result!(wta, clock::Clock, integrand=wta.operand) # Time increment: - Δt = clock.time - wta.previous_collection_time + Δt = period_to_number(clock.time - wta.previous_collection_time) # Time intervals: - T_current = clock.time - wta.window_start_time - T_previous = wta.previous_collection_time - wta.window_start_time + T_current = period_to_number(clock.time - wta.window_start_time) + T_previous = period_to_number(wta.previous_collection_time - wta.window_start_time) # Accumulate left Riemann sum @. wta.result = (wta.result * T_previous + integrand * Δt) / T_current diff --git a/src/Utils/times_and_datetimes.jl b/src/Utils/times_and_datetimes.jl index a2c0df07f8..89bcea6d9a 100644 --- a/src/Utils/times_and_datetimes.jl +++ b/src/Utils/times_and_datetimes.jl @@ -29,7 +29,9 @@ end @inline add_time_interval(base::AbstractTime, interval::Number, count=1) = base + seconds_to_nanosecond(interval * count) @inline add_time_interval(base::AbstractTime, interval::Period, count=1) = base + count * interval -@inline add_time_interval(base::Number, interval::Array{<:Number}, count=1) = interval[count] +@inline add_time_interval(base, interval::Array{<:Number}, count=1) = interval[count] +@inline add_time_interval(base, interval::Array{<:Dates.Period}, count=1) = seconds_to_nanosecond(interval[count]) +@inline add_time_interval(base, interval::Array{Dates.DateTime}, count=1) = interval[count] function period_type(interval::Number) FT = Oceananigans.defaults.FloatType @@ -42,6 +44,8 @@ function period_type(interval::Array{<:Number}) end period_type(interval::Dates.Period) = typeof(interval) +period_type(interval) = typeof(interval) + time_type(interval::Number) = typeof(interval) time_type(interval::Dates.Period) = Dates.DateTime -time_type(interval::Array{<:Number}) = eltype(interval) +time_type(interval::Array) = eltype(interval) \ No newline at end of file