diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 4657501702..faee3cd8f7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -15,8 +15,9 @@ steps: - "wget -N -P $SVERDRUP_HOME https://julialang-s3.julialang.org/bin/linux/x64/$JULIA_MINOR_VERSION/julia-$JULIA_VERSION-linux-x86_64.tar.gz" - "tar xf $SVERDRUP_HOME/julia-$JULIA_VERSION-linux-x86_64.tar.gz -C $SVERDRUP_HOME" - # Instantiate and precompile + # Instantiate, build, and precompile - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia --color=yes --project -e 'using Pkg; Pkg.instantiate(; verbose=true)'" + - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia --color=yes --project -e 'using Pkg; Pkg.build()'" - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia --color=yes --project -e 'using Pkg; Pkg.precompile()'" - "$SVERDRUP_HOME/julia-$JULIA_VERSION/bin/julia --color=yes --project -e 'using Pkg; Pkg.status()'" diff --git a/docs/src/library.md b/docs/src/library.md index 84ddb5671f..79e7178453 100644 --- a/docs/src/library.md +++ b/docs/src/library.md @@ -168,6 +168,7 @@ Pages = [ "OutputWriters/output_writer_utils.jl", "OutputWriters/jld2_output_writer.jl", "OutputWriters/netcdf_output_writer.jl", + "OutputWriters/windowed_time_average.jl", "OutputWriters/checkpointer.jl" ] ``` diff --git a/docs/src/model_setup/checkpointing.md b/docs/src/model_setup/checkpointing.md index 5270ce918a..a7b1ae4438 100644 --- a/docs/src/model_setup/checkpointing.md +++ b/docs/src/model_setup/checkpointing.md @@ -14,12 +14,12 @@ end ``` ```jldoctest -julia> using Oceananigans.OutputWriters; +julia> using Oceananigans.OutputWriters, Oceananigans.Utils julia> model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))); -julia> checkpointer = Checkpointer(model, time_interval=1000000, prefix="model_checkpoint") -Checkpointer{Nothing,Int64,Array{Symbol,1}}(nothing, 1000000, 0.0, ".", "model_checkpoint", [:architecture, :grid, :clock, :coriolis, :buoyancy, :closure, :velocities, :tracers, :timestepper], false, false) +julia> checkpointer = Checkpointer(model, schedule=TimeInterval(5years), prefix="model_checkpoint") +Checkpointer{TimeInterval,Array{Symbol,1}}(TimeInterval(1.5768e8, 0.0), ".", "model_checkpoint", [:architecture, :grid, :clock, :coriolis, :buoyancy, :closure, :velocities, :tracers, :timestepper], false, false) ``` The default options should provide checkpoint files that are easy to restore from in most cases. For more advanced diff --git a/docs/src/model_setup/output_writers.md b/docs/src/model_setup/output_writers.md index 59c9559efd..90e6f1d68d 100644 --- a/docs/src/model_setup/output_writers.md +++ b/docs/src/model_setup/output_writers.md @@ -1,20 +1,52 @@ # Output writers -Saving model data to disk can be done in a flexible manner using output writers. The two main output writers currently -implemented are a NetCDF output writer (relying on [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl)) -and a JLD2 output writer (relying on [JLD2.jl](https://github.com/JuliaIO/JLD2.jl)). +`AbstractOutputWriter`s save data to disk. +`Oceananigans` provides three ways to write output: -Output writers are stored as a list of output writers in `simulation.output_writers`. Output writers can be specified -at model creation time or be specified at any later time and appended (or assigned with a key value pair) to -`simulation.output_writers`. +1. [`NetCDFOutputWriter`](@ref) for output of arrays and scalars that uses [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) +2. [`JLD2OutputWriter`](@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 on a separate documentation page. + +## Basic usage + +[`NetCDFOutputWriter`](@ref) and [`JLD2OutputWriter`](@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. `JLD2OutputWriter` accepts `NamedTuple`s and `Dict`s; + `NetCDFOutputWriter` 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. +4. The filename and directory. Currently `NetCDFOutputWriter` accepts one `filepath` argument, while + `JLD2OutputWriter` accepts a filename `prefix` and `dir`ectory. + +Other important keyword arguments are + +* `field_slicer::FieldSlicer` for outputting subregions, two- and one-dimensional slices of fields. + By default a `FieldSlicer` is used to remove halo regions from fields so that only the physical + portion of model data is saved to disk. + +* `array_type` for specifying the type of the array that holds outputted field data. The default is + `Array{Float32}`, 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 `NetCDFOutputWriter` and `JLD2OutputWriter` is given below. ## 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. -The following example shows how to construct NetCDF output writers for two different kinds of outputs (3D fields and -slices) along with output attributes +### Examples + +Saving the u velocity field and temperature fields, the full 3D fields and surface 2D slices +to separate NetCDF files: ```jldoctest netcdf1 using Oceananigans, Oceananigans.OutputWriters @@ -28,26 +60,48 @@ simulation = Simulation(model, Δt=12, stop_time=3600); fields = Dict("u" => model.velocities.u, "T" => model.tracers.T); simulation.output_writers[:field_writer] = - NetCDFOutputWriter(model, fields, filepath="output_fields.nc", time_interval=60) + NetCDFOutputWriter(model, fields, filepath="more_fields.nc", schedule=TimeInterval(60)) # output -NetCDFOutputWriter (time_interval=60): output_fields.nc +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: more_fields.nc ├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0) -└── 2 outputs: ["T", "u"] +├── 2 outputs: ["T", "u"] +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +└── array type: Array{Float32} ``` - ```jldoctest netcdf1 simulation.output_writers[:surface_slice_writer] = - NetCDFOutputWriter(model, fields, filepath="output_surface_xy_slice.nc", - time_interval=60, field_slicer=FieldSlicer(k=grid.Nz)) + NetCDFOutputWriter(model, fields, filepath="another_surface_xy_slice.nc", + schedule=TimeInterval(60), field_slicer=FieldSlicer(k=grid.Nz)) # output -NetCDFOutputWriter (time_interval=60): output_surface_xy_slice.nc +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: another_surface_xy_slice.nc ├── dimensions: zC(1), zF(1), xC(16), yF(16), xF(16), yC(16), time(0) -└── 2 outputs: ["T", "u"] +├── 2 outputs: ["T", "u"] +├── field slicer: FieldSlicer(:, :, 16, with_halos=false) +└── array type: Array{Float32} ``` -Writing a scalar, profile, and slice to NetCDF: +```jldoctest netcdf1 +simulation.output_writers[:averaged_profile_writer] = + NetCDFOutputWriter(model, fields, + filepath = "another_averaged_z_profile.nc", + schedule = AveragedTimeInterval(60, window=20), + field_slicer = FieldSlicer(i=1, j=1)) + +# output +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: another_averaged_z_profile.nc +├── dimensions: zC(16), zF(17), xC(1), yF(1), xF(1), yC(1), time(0) +├── 2 outputs: ["T", "u"] averaged on AveragedTimeInterval(window=20 seconds, stride=1, interval=1 minute) +├── field slicer: FieldSlicer(1, 1, :, with_halos=false) +└── array type: Array{Float32} +``` + +`NetCDFOutputWriter` also accepts output functions that write scalars and arrays to disk, +provided that their `dimensions` are provided: ```jldoctest using Oceananigans, Oceananigans.OutputWriters @@ -59,9 +113,7 @@ model = IncompressibleModel(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(Cell, grid)); # vector/profile output - h(model) = model.clock.time .* ( sin.(xnodes(Cell, grid, reshape=true)[:, :, 1]) .* cos.(ynodes(Face, grid, reshape=true)[:, :, 1])); # xy slice output @@ -77,54 +129,139 @@ output_attributes = Dict( global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7); -simulation.output_writers[:stuff] = +simulation.output_writers[:things] = NetCDFOutputWriter(model, outputs, - iteration_interval=1, filepath="stuff.nc", dimensions=dims, verbose=true, + schedule=IterationInterval(1), filepath="some_things.nc", dimensions=dims, verbose=true, global_attributes=global_attributes, output_attributes=output_attributes) # output -NetCDFOutputWriter (iteration_interval=1): stuff.nc +NetCDFOutputWriter scheduled on IterationInterval(1): +├── filepath: some_things.nc ├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0) -└── 3 outputs: ["profile", "slice", "scalar"] +├── 3 outputs: ["profile", "slice", "scalar"] +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +└── array type: Array{Float32} ``` -See [`NetCDFOutputWriter`](@ref) for more details and options. +See [`NetCDFOutputWriter`](@ref) for more information. ## JLD2 output writer -JLD2 is a an HDF5 compatible file format written in pure Julia and is generally pretty fast. JLD2 files can be opened in -Python with the [h5py](https://www.h5py.org/) package. +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 JLD2 output writer is generally used by passing it a dictionary or named tuple of (label, function) pairs where the -functions have a single input `model`. Whenever output needs to be written, the functions will be called and the output -of the function will be saved to the JLD2 file. For example, to write out 3D fields for w and T and a horizontal average -of T every 1 hour of simulation time to a file called `some_data.jld2` +The `JLD2OutputWriter` 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. -```julia -using Oceananigans -using Oceananigans.OutputWriters +### Examples + +Write out 3D fields for w and T and a horizontal average: +```jldoctest jld2_output_writer +using Oceananigans, Oceananigans.OutputWriters, Oceananigans.Fields using Oceananigans.Utils: hour, minute -model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))) +model = IncompressibleModel(grid=RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))) + simulation = Simulation(model, Δt=12, stop_time=1hour) -function init_save_some_metadata(file, model) +function init_save_some_metadata!(file, model) file["author"] = "Chim Riggles" file["parameters/coriolis_parameter"] = 1e-4 file["parameters/density"] = 1027 + return nothing end T_avg = AveragedField(model.tracers.T, dims=(1, 2)) -outputs = Dict( - :w => model -> model.velocities.u, - :T => model -> model.tracers.T, - :T_avg => model -> T_avg(model) -) +# Note that model.velocities is NamedTuple +simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, + prefix = "some_more_data", + schedule = TimeInterval(20minute), + init = init_save_some_metadata!) + +# output +JLD2OutputWriter scheduled on TimeInterval(20 minutes): +├── filepath: ./some_more_data.jld2 +├── 3 outputs: (:u, :v, :w) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` + +and a time- and horizontal-average of temperature `T` every 1 hour of simulation time +to a file called `some_averaged_data.jld2` + +```jldoctest jld2_output_writer +simulation.output_writers[:avg_T] = JLD2OutputWriter(model, (T=T_avg,), + prefix = "some_more_averaged_data", + schedule = AveragedTimeInterval(20minute, window=5minute)) +# output +JLD2OutputWriter scheduled on TimeInterval(20 minutes): +├── filepath: ./some_more_averaged_data.jld2 +├── 1 outputs: (:T,) averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=20 minutes) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` + +See [`JLD2OutputWriter`](@ref) for more information. -jld2_writer = JLD2OutputWriter(model, outputs, init=init_save_some_metadata, interval=20minute, prefix="some_data") +## Time-averaged output -push!(simulation.output_writers, jld2_writer) +Time-averaged output is specified by setting the `schedule` keyword argument for either `NetCDFOutputWriter` or +`JLD2OutputWriter` 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 year window, every 4 years, + +```jldoctest averaged_time_interval +using Oceananigans.OutputWriters: AveragedTimeInterval +using Oceananigans.Utils: year, years + +schedule = AveragedTimeInterval(4years, window=1year) -See [`JLD2OutputWriter`](@ref) for more details and options. +# output +AveragedTimeInterval(window=1 year, stride=1, interval=4 years) +``` + +An `AveragedTimeInterval` schedule directs an output writer +to time-average its outputs before writing them to disk: + +```jldoctest averaged_time_interval +using Oceananigans +using Oceananigans.OutputWriters: JLD2OutputWriter +using Oceananigans.Utils: minutes + +model = IncompressibleModel(grid=RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=10minutes, stop_time=30years) + +simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, + prefix = "even_more_averaged_velocity_data", + schedule = AveragedTimeInterval(4years, window=1year, stride=2)) + +# output +JLD2OutputWriter scheduled on TimeInterval(4 years): +├── filepath: ./even_more_averaged_velocity_data.jld2 +├── 3 outputs: (:u, :v, :w) averaged on AveragedTimeInterval(window=1 year, stride=2, interval=4 years) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` diff --git a/examples/eady_turbulence.jl b/examples/eady_turbulence.jl index b6c4e5fdc1..641f6e3453 100644 --- a/examples/eady_turbulence.jl +++ b/examples/eady_turbulence.jl @@ -287,12 +287,12 @@ nothing # hide # we create a `JLD2OutputWriter` that saves `ζ` and `δ` and add it to # `simulation`. -using Oceananigans.OutputWriters: JLD2OutputWriter +using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval simulation.output_writers[:fields] = JLD2OutputWriter(model, (ζ=ζ, δ=δ), - time_interval = 2hour, - prefix = "eady_turbulence", - force = true) + schedule = TimeInterval(2hour), + prefix = "eady_turbulence", + force = true) nothing # hide # All that's left is to press the big red button: diff --git a/examples/internal_wave.jl b/examples/internal_wave.jl index 5b01dea73d..a70e5e8b0c 100644 --- a/examples/internal_wave.jl +++ b/examples/internal_wave.jl @@ -118,12 +118,12 @@ simulation = Simulation(model, Δt = 0.02 * 2π/ω, stop_iteration = 100) # and add an output writer that saves the vertical velocity field every two iterations: -using Oceananigans.OutputWriters: JLD2OutputWriter +using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, - iteration_interval = 2, - prefix = "internal_wave", - force = true) + schedule = IterationInterval(2), + prefix = "internal_wave", + force = true) # With initial conditions set and an output writer at the ready, we run the simulation diff --git a/examples/langmuir_turbulence.jl b/examples/langmuir_turbulence.jl index a5103ffd7d..8c57a09e0e 100644 --- a/examples/langmuir_turbulence.jl +++ b/examples/langmuir_turbulence.jl @@ -215,7 +215,7 @@ using Oceananigans.Utils: minute fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,)) simulation.output_writers[:fields] = - JLD2OutputWriter(model, fields_to_output, time_interval = 2minute, + JLD2OutputWriter(model, fields_to_output, schedule=TimeInterval(2minute), prefix = "langmuir_turbulence", force = true) nothing # hide diff --git a/examples/ocean_wind_mixing_and_convection.jl b/examples/ocean_wind_mixing_and_convection.jl index 4275aad3c7..7113d149ee 100644 --- a/examples/ocean_wind_mixing_and_convection.jl +++ b/examples/ocean_wind_mixing_and_convection.jl @@ -129,7 +129,7 @@ nothing # hide ## Instantiate a JLD2OutputWriter to write fields. We will add it to the simulation before ## running it. -field_writer = JLD2OutputWriter(model, fields_to_output; time_interval=hour/4, +field_writer = JLD2OutputWriter(model, fields_to_output; schedule=TimeInterval(hour/4), prefix="ocean_wind_mixing_and_convection", force=true) # ## Running the simulation diff --git a/examples/one_dimensional_diffusion.jl b/examples/one_dimensional_diffusion.jl index 2faf9b3de8..78ca461c66 100644 --- a/examples/one_dimensional_diffusion.jl +++ b/examples/one_dimensional_diffusion.jl @@ -87,11 +87,11 @@ plot!(p, interior(T)[1, 1, :], z, linewidth=2, label=tracer_label(model.clock.ti # Interesting! Next, we add an output writer that saves the temperature field # in JLD2 files, and run the simulation for longer so that we can animate the results. -using Oceananigans.OutputWriters: JLD2OutputWriter +using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval simulation.output_writers[:temperature] = JLD2OutputWriter(model, model.tracers, prefix = "one_dimensional_diffusion", - iteration_interval = 100, force = true) + schedule=IterationInterval(100), force = true) ## Run simulation for 10,000 more iterations simulation.stop_iteration += 10000 diff --git a/examples/two_dimensional_turbulence.jl b/examples/two_dimensional_turbulence.jl index 7a02fbb842..f14864daa0 100644 --- a/examples/two_dimensional_turbulence.jl +++ b/examples/two_dimensional_turbulence.jl @@ -60,7 +60,7 @@ simulation = Simulation(model, Δt=0.1, stop_iteration=2000) using Oceananigans.OutputWriters simulation.output_writers[:fields] = - JLD2OutputWriter(model, (ω = ω,), iteration_interval = 20, + JLD2OutputWriter(model, (ω = ω,), schedule=IterationInterval(20), prefix = "2d_turbulence_vorticity", force = true) # ## Running the simulation diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl index 2a9477d35d..cde101d8bf 100644 --- a/src/Diagnostics/Diagnostics.jl +++ b/src/Diagnostics/Diagnostics.jl @@ -1,17 +1,19 @@ module Diagnostics export - TimeSeries, FieldMaximum, - CFL, AdvectiveCFL, DiffusiveCFL, NaNChecker, - run_diagnostic! + NaNChecker, + FieldMaximum, + CFL, AdvectiveCFL, DiffusiveCFL, + run_diagnostic!, + TimeInterval, IterationInterval, WallTimeInterval using Oceananigans using Oceananigans.Operators +using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval using Oceananigans: AbstractDiagnostic include("nan_checker.jl") -include("time_series.jl") include("field_maximum.jl") include("cfl.jl") diff --git a/src/Diagnostics/cfl.jl b/src/Diagnostics/cfl.jl index 75ae19fe40..123b8d63c7 100644 --- a/src/Diagnostics/cfl.jl +++ b/src/Diagnostics/cfl.jl @@ -7,7 +7,7 @@ using Oceananigans.TurbulenceClosures: cell_diffusion_timescale An object for computing the Courant-Freidrichs-Lewy (CFL) number. """ struct CFL{D, S} - Δt :: D + Δt :: D timescale :: S end diff --git a/src/Diagnostics/nan_checker.jl b/src/Diagnostics/nan_checker.jl index 2a13ec57df..5dc75e98cf 100644 --- a/src/Diagnostics/nan_checker.jl +++ b/src/Diagnostics/nan_checker.jl @@ -3,21 +3,18 @@ A diagnostic that checks for `NaN` values and aborts the simulation if any are found. """ -struct NaNChecker{F} <: AbstractDiagnostic - iteration_interval :: Int - fields :: F +struct NaNChecker{T, F} <: AbstractDiagnostic + schedule :: T + fields :: F end """ - NaNChecker(model; iteration_interval, fields) + NaNChecker(; schedule, fields) -Construct a `NaNChecker` for `model`. `fields` should be a `Dict{Symbol,Field}`. An -`iteration_interval` should be passed to indicate how often to check for NaNs (in -number of iterations). +Returns a `NaNChecker` that checks for `NaN` anywhere within `fields` +when `schedule` actuates. """ -function NaNChecker(model; iteration_interval, fields) - return NaNChecker(iteration_interval, fields) -end +NaNChecker(model=nothing; schedule, fields) = NaNChecker(schedule, fields) function run_diagnostic!(nc::NaNChecker, model) for (name, field) in nc.fields diff --git a/src/Diagnostics/time_series.jl b/src/Diagnostics/time_series.jl deleted file mode 100644 index 19438199bc..0000000000 --- a/src/Diagnostics/time_series.jl +++ /dev/null @@ -1,130 +0,0 @@ -""" - TimeSeries{D, II, TI, T, TT} <: AbstractDiagnostic - -A diagnostic for collecting and storing time series. -""" -struct TimeSeries{D, II, TI, T, TT} <: AbstractDiagnostic - diagnostic :: D - iteration_interval :: II - time_interval :: TI - data :: T - time :: Vector{TT} -end - -""" - TimeSeries(diagnostic, model; iteration_interval=nothing, time_interval=nothing) - -A `TimeSeries` `Diagnostic` that records a time series of `diagnostic(model)`. - -Example -======= - -```jldoctest timeseries1 -julia> using Oceananigans, Oceananigans.Diagnostics - -julia> model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))); - -julia> sim = Simulation(model, Δt=1, stop_iteration=3); - -julia> max_u = TimeSeries(FieldMaximum(abs, model.velocities.u), model; iteration_interval=1); - -julia> sim.diagnostics[:max_u] = max_u; - -julia> run!(sim); -[ Info: Simulation is stopping. Model iteration 3 has hit or exceeded simulation stop iteration 3. -``` - -```jldoctest timeseries1 -julia> max_u.data -4-element Array{Float64,1}: - 0.0 - 0.0 - 0.0 - 0.0 -``` -""" -function TimeSeries(diagnostic, model; iteration_interval=nothing, time_interval=nothing) - TD = typeof(diagnostic(model)) - TT = typeof(model.clock.time) - return TimeSeries(diagnostic, iteration_interval, time_interval, TD[], TT[]) -end - -@inline (time_series::TimeSeries)(model) = time_series.diagnostic(model) - -function run_diagnostic!(diag::TimeSeries, model) - push!(diag.data, diag(model)) - push!(diag.time, model.clock.time) - return nothing -end - -""" - TimeSeries(diagnostics::NamedTuple, model; iteration_interval=nothing, time_interval=nothing) - -A `TimeSeries` `Diagnostic` that records a `NamedTuple` of time series of -`diag(model)` for each `diag` in `diagnostics`. - -Example -======= - -```jldoctest timeseries2 -julia> using Oceananigans, Oceananigans.Diagnostics - -julia> model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))); - -julia> Δt = 1.0; - -julia> cfl = TimeSeries((adv=AdvectiveCFL(Δt), diff=DiffusiveCFL(Δt)), model; iteration_interval=1); - -julia> sim = Simulation(model, Δt=Δt, stop_iteration=3); - -julia> sim.diagnostics[:cfl] = cfl; - -julia> run!(sim); -[ Info: Simulation is stopping. Model iteration 3 has hit or exceeded simulation stop iteration 3. -``` - -```jldoctest timeseries2 -julia> cfl.data -(adv = [0.0, 0.0, 0.0, 0.0], diff = [0.0002688, 0.0002688, 0.0002688, 0.0002688]) -``` - -``` jldoctest timeseries2 -julia> cfl.adv -4-element Array{Float64,1}: - 0.0 - 0.0 - 0.0 - 0.0 -``` -""" -function TimeSeries(diagnostics::NamedTuple, model; iteration_interval=nothing, time_interval=nothing) - TT = typeof(model.clock.time) - TDs = Tuple(typeof(diag(model)) for diag in diagnostics) - data = NamedTuple{propertynames(diagnostics)}(Tuple(T[] for T in TDs)) - return TimeSeries(diagnostics, iteration_interval, time_interval, data, TT[]) -end - -function (time_series::TimeSeries{<:NamedTuple})(model) - names = propertynames(time_series.diagnostic) - return NamedTuple{names}(Tuple(diag(model) for diag in time_series.diagnostics)) -end - -function run_diagnostic!(diag::TimeSeries{<:NamedTuple}, model) - ntuple(Val(length(diag.diagnostic))) do i - Base.@_inline_meta - push!(diag.data[i], diag.diagnostic[i](model)) - end - push!(diag.time, model.clock.time) - return nothing -end - -Base.getproperty(ts::TimeSeries{<:NamedTuple}, name::Symbol) = get_time_series_field(ts, name) - -"Returns a field of time series, or a field of `time_series.data` when possible." -function get_time_series_field(time_series, name) - if name ∈ propertynames(time_series) - return getfield(time_series, name) - else - return getproperty(time_series.data, name) - end -end diff --git a/src/OutputWriters/OutputWriters.jl b/src/OutputWriters/OutputWriters.jl index 3fbb121a0c..32ba5d0e7d 100644 --- a/src/OutputWriters/OutputWriters.jl +++ b/src/OutputWriters/OutputWriters.jl @@ -6,7 +6,8 @@ export JLD2OutputWriter, NetCDFOutputWriter, Checkpointer, restore_from_checkpoint, - WindowedTimeAverage + WindowedTimeAverage, + TimeInterval, IterationInterval, WallTimeInterval, AveragedTimeInterval using CUDA @@ -14,6 +15,7 @@ using Oceananigans.Architectures using Oceananigans.Grids using Oceananigans.Fields using Oceananigans.Models +using Oceananigans.Utils: TimeInterval, IterationInterval, WallTimeInterval using Oceananigans: AbstractOutputWriter using Oceananigans.Fields: OffsetArray @@ -25,6 +27,7 @@ include("output_writer_utils.jl") include("field_slicer.jl") include("fetch_output.jl") include("windowed_time_average.jl") +include("time_average_outputs.jl") include("jld2_output_writer.jl") include("netcdf_output_writer.jl") include("checkpointer.jl") diff --git a/src/OutputWriters/checkpointer.jl b/src/OutputWriters/checkpointer.jl index bc6a795912..4da8fb78cb 100644 --- a/src/OutputWriters/checkpointer.jl +++ b/src/OutputWriters/checkpointer.jl @@ -5,55 +5,64 @@ using Oceananigans.Fields: offset_data An output writer for checkpointing models to a JLD2 file from which models can be restored. """ -mutable struct Checkpointer{I, T, P} <: AbstractOutputWriter - iteration_interval :: I - time_interval :: T - previous :: Float64 - dir :: String - prefix :: String - properties :: P - force :: Bool - verbose :: Bool +mutable struct Checkpointer{T, P} <: AbstractOutputWriter + schedule :: T + dir :: String + prefix :: String + properties :: P + force :: Bool + verbose :: Bool end """ - Checkpointer(model; iteration_interval=nothing, time_interval=nothing, dir=".", - prefix="checkpoint", force=false, verbose=false, + Checkpointer(model; schedule, + dir = ".", + prefix = "checkpoint", + force = false, + verbose = false, properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis, - :buoyancy, :closure, :velocities, :tracers, :timestepper]) + :buoyancy, :closure, :velocities, :tracers, :timestepper] + ) Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as -specified by `iteration_interval` or `time_interval`. The `model.clock.iteration` is included +specified by `schedule`. The `model.clock.iteration` is included in the filename to distinguish between multiple checkpoint files. Note that extra model `properties` can be safely specified, but removing crucial properties such as `:velocities` will make restoring from the checkpoint impossible. The checkpoint file is generated by serializing model properties to JLD2. However, -functions cannot be serialized to disk (at least not with JLD2). So if a model property +functions cannot be serialized to disk with JLD2. So if a model property contains a reference somewhere in its hierarchy it will not be included in the checkpoint file (and you will have to manually restore them). Keyword arguments ================= -- `iteration_interval`: Save output every `n` model iterations. -- `time_interval`: Save output every `t` units of model clock time. +- `schedule` (required): Schedule that determines when to checkpoint. + - `dir`: Directory to save output to. Default: "." (current working directory). + - `prefix`: Descriptive filename prefixed to all output files. Default: "checkpoint". + - `force`: Remove existing files if their filenames conflict. Default: `false`. + - `verbose`: Log what the output writer is doing with statistics on compute/write times - and file sizes. Default: `false`. + and file sizes. Default: `false`. + - `properties`: List of model properties to checkpoint. Some are required. """ -function Checkpointer(model; iteration_interval=nothing, time_interval=nothing, - dir=".", prefix="checkpoint", force=false, verbose=false, +function Checkpointer(model; schedule, + dir = ".", + prefix = "checkpoint", + force = false, + verbose = false, properties = [:architecture, :grid, :clock, :coriolis, - :buoyancy, :closure, :velocities, :tracers, :timestepper]) - - validate_intervals(iteration_interval, time_interval) + :buoyancy, :closure, :velocities, :tracers, :timestepper] + ) # Certain properties are required for `restore_from_checkpoint` to work. required_properties = (:grid, :architecture, :velocities, :tracers, :timestepper) + for rp in required_properties if rp ∉ properties @warn "$rp is required for checkpointing. It will be added to checkpointed properties" @@ -72,7 +81,8 @@ function Checkpointer(model; iteration_interval=nothing, time_interval=nothing, end mkpath(dir) - return Checkpointer(iteration_interval, time_interval, 0.0, dir, prefix, properties, force, verbose) + + return Checkpointer(schedule, dir, prefix, properties, force, verbose) end function write_output!(c::Checkpointer, model) diff --git a/src/OutputWriters/field_slicer.jl b/src/OutputWriters/field_slicer.jl index e078ee9ccd..618dc64a37 100644 --- a/src/OutputWriters/field_slicer.jl +++ b/src/OutputWriters/field_slicer.jl @@ -1,3 +1,5 @@ +import Oceananigans: short_show + """ struct FieldSlicer{I, J, K, W} @@ -88,3 +90,11 @@ function slice_parent(slicer, field) end slice_parent(::Nothing, field) = parent(field) + +show_index(i) = string(i) +show_index(i::Colon) = ":" + +short_show(fs::FieldSlicer) = string("FieldSlicer(", + "$(show_index(fs.i)), ", + "$(show_index(fs.j)), ", + "$(show_index(fs.k)), with_halos=$(fs.with_halos))") diff --git a/src/OutputWriters/jld2_output_writer.jl b/src/OutputWriters/jld2_output_writer.jl index 0c93268ca8..ca710d7f3b 100644 --- a/src/OutputWriters/jld2_output_writer.jl +++ b/src/OutputWriters/jld2_output_writer.jl @@ -1,22 +1,21 @@ using Printf using JLD2 using Oceananigans.Utils +using Oceananigans.Utils: TimeInterval, pretty_filesize """ JLD2OutputWriter{I, T, O, IF, IN, KW} <: AbstractOutputWriter An output writer for writing to JLD2 files. """ -mutable struct JLD2OutputWriter{I, T, O, FS, D, IF, IN, KW} <: AbstractOutputWriter +mutable struct JLD2OutputWriter{O, T, FS, D, IF, IN, KW} <: AbstractOutputWriter filepath :: String outputs :: O - iteration_interval :: I - time_interval :: T + schedule :: T field_slicer :: FS array_type :: D init :: IF including :: IN - previous :: Float64 part :: Int max_filesize :: Float64 force :: Bool @@ -27,21 +26,17 @@ end noinit(args...) = nothing """ - JLD2OutputWriter(model, outputs; prefix, - dir = ".", - iteration_interval = nothing, - time_interval = nothing, - time_averaging_window = nothing, - time_averaging_stride = 1, - field_slicer = FieldSlicer(), - array_type = Array{Float32}, - max_filesize = Inf, - force = false, - init = noinit, - including = [:grid, :coriolis, :buoyancy, :closure], - verbose = false, - part = 1, - jld2_kw = Dict{Symbol, Any}()) + JLD2OutputWriter(model, outputs; prefix, schedule, + dir = ".", + field_slicer = FieldSlicer(), + array_type = Array{Float32}, + max_filesize = Inf, + force = false, + init = noinit, + including = [:grid, :coriolis, :buoyancy, :closure], + verbose = false, + part = 1, + jld2_kw = Dict{Symbol, Any}()) Construct a `JLD2OutputWriter` for an Oceananigans `model` that writes `label, output` pairs in `outputs` to a JLD2 file. @@ -56,27 +51,15 @@ Keyword arguments ## Filenaming - - `prefix`: Descriptive filename prefixed to all output files. + - `prefix` (required): Descriptive filename prefixed to all output files. - `dir`: Directory to save output to. Default: "." (current working directory). ## Output frequency and time-averaging - - `iteration_interval`: Save output every `iteration_interval` model iterations. + - `schedule` (required): `AbstractSchedule` that determines when output is saved. - - `time_interval`: Save output every `time_interval` units of `model.clock.time`. - - - `time_averaging_window`: Specifies a time window over which each member of `output` is averaged before - being saved. For this each member of output is converted to - `Oceananigans.Diagnostics.WindowedTimeAverage`. - Default `nothing` indicates no averaging. - - - `time_averaging_stride`: Specifies a iteration 'stride' between the calculation of each `output` during - time-averaging. Longer strides means that output is calculated less frequently, - and that the resulting time-average is faster to compute, but less accurate. - Default: 1. - ## Slicing and type conversion prior to output - `field_slicer`: An object for slicing field output in ``(x, y, z)``, including omitting halos. @@ -113,44 +96,77 @@ Keyword arguments Default: 1. - `jld2_kw`: Dict of kwargs to be passed to `jldopen` when data is written. -""" -function JLD2OutputWriter(model, outputs; prefix, - dir = ".", - iteration_interval = nothing, - time_interval = nothing, - time_averaging_window = nothing, - time_averaging_stride = 1, - field_slicer = FieldSlicer(), - array_type = Array{Float32}, - max_filesize = Inf, - force = false, - init = noinit, - including = [:grid, :coriolis, :buoyancy, :closure], - verbose = false, - part = 1, - jld2_kw = Dict{Symbol, Any}()) - - validate_intervals(iteration_interval, time_interval) - - # Convert each output to WindowedTimeAverage if time_averaging_window is specified - if !isnothing(time_averaging_window) - - !isnothing(iteration_interval) && error("Cannot specify iteration_interval with time_averaging_window.") - - output_names = Tuple(keys(outputs)) - - averaged_output = - Tuple( - WindowedTimeAverage(outputs[name], model; time_interval = time_interval, - time_window = time_averaging_window, - stride = time_averaging_stride, - field_slicer = field_slicer) - for name in output_names - ) - - outputs = NamedTuple{output_names}(averaged_output) - end +Example +======= + +Write out 3D fields for w and T and a horizontal average: + +```jldoctest jld2_output_writer +using Oceananigans, Oceananigans.OutputWriters, Oceananigans.Fields +using Oceananigans.Utils: hour, minute + +model = IncompressibleModel(grid=RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))) +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 + +T_avg = AveragedField(model.tracers.T, dims=(1, 2)) + +# Note that model.velocities is NamedTuple +simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, + prefix = "some_data", + schedule = TimeInterval(20minute), + init = init_save_some_metadata!) + +# output +JLD2OutputWriter scheduled on TimeInterval(20 minutes): +├── filepath: ./some_data.jld2 +├── 3 outputs: (:u, :v, :w) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` + +and a time- and horizontal-average of temperature `T` every 1 hour of simulation time +to a file called `some_averaged_data.jld2` + +```jldoctest jld2_output_writer +simulation.output_writers[:avg_T] = JLD2OutputWriter(model, (T=T_avg,), + prefix = "some_averaged_data", + schedule = AveragedTimeInterval(20minute, window=5minute)) + +# output +JLD2OutputWriter scheduled on TimeInterval(20 minutes): +├── filepath: ./some_averaged_data.jld2 +├── 1 outputs: (:T,) averaged on AveragedTimeInterval(window=5 minutes, stride=1, interval=20 minutes) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` +""" +function JLD2OutputWriter(model, outputs; prefix, schedule, + dir = ".", + field_slicer = FieldSlicer(), + array_type = Array{Float32}, + max_filesize = Inf, + force = false, + init = noinit, + including = [:grid, :coriolis, :buoyancy, :closure], + verbose = false, + part = 1, + jld2_kw = Dict{Symbol, Any}()) + + # Convert each output to WindowedTimeAverage if schedule::AveragedTimeWindow is specified + schedule, outputs = time_average_outputs(schedule, outputs, model, field_slicer) + mkpath(dir) filepath = joinpath(dir, prefix * ".jld2") force && isfile(filepath) && rm(filepath, force=true) @@ -160,20 +176,9 @@ function JLD2OutputWriter(model, outputs; prefix, saveproperties!(file, model, including) end - return JLD2OutputWriter(filepath, - outputs, - iteration_interval, - time_interval, - field_slicer, - array_type, - init, - including, - 0.0, - part, - max_filesize, - force, - verbose, - jld2_kw) + return JLD2OutputWriter(filepath, outputs, schedule, field_slicer, + array_type, init, including, part, max_filesize, + force, verbose, jld2_kw) end function write_output!(writer::JLD2OutputWriter, model) @@ -253,3 +258,16 @@ function start_next_file(model, writer::JLD2OutputWriter) return nothing end + +function Base.show(io::IO, ow::JLD2OutputWriter) + + averaging_schedule = output_averaging_schedule(ow) + + print(io, "JLD2OutputWriter scheduled on $(show_schedule(ow.schedule)):", '\n', + "├── filepath: $(ow.filepath)", '\n', + "├── $(length(ow.outputs)) outputs: $(keys(ow.outputs))", show_averaging_schedule(averaging_schedule), '\n', + "├── field slicer: $(short_show(ow.field_slicer))", '\n', + "├── array type: ", show_array_type(ow.array_type), '\n', + "├── including: ", ow.including, '\n', + "└── max filesize: ", pretty_filesize(ow.max_filesize)) +end diff --git a/src/OutputWriters/netcdf_output_writer.jl b/src/OutputWriters/netcdf_output_writer.jl index a4f72376a7..8e1b259632 100644 --- a/src/OutputWriters/netcdf_output_writer.jl +++ b/src/OutputWriters/netcdf_output_writer.jl @@ -1,20 +1,26 @@ using NCDatasets using Oceananigans.Fields +using Oceananigans.Utils: show_schedule using Dates: now using Oceananigans.Grids: topology, halo_size -using Oceananigans.Utils: validate_intervals, versioninfo_with_gpu, oceananigans_versioninfo +using Oceananigans.Utils: versioninfo_with_gpu, oceananigans_versioninfo -xdim(::Type{Face}) = "xF" -ydim(::Type{Face}) = "yF" -zdim(::Type{Face}) = "zF" +xdim(::Type{Face}) = ("xF",) +ydim(::Type{Face}) = ("yF",) +zdim(::Type{Face}) = ("zF",) -xdim(::Type{Cell}) = "xC" -ydim(::Type{Cell}) = "yC" -zdim(::Type{Cell}) = "zC" +xdim(::Type{Cell}) = ("xC",) +ydim(::Type{Cell}) = ("yC",) +zdim(::Type{Cell}) = ("zC",) -netcdf_spatial_dimensions(::Field{LX, LY, LZ}) where {LX, LY, LZ} = xdim(LX), ydim(LY), zdim(LZ) +xdim(::Type{Nothing}) = () +ydim(::Type{Nothing}) = () +zdim(::Type{Nothing}) = () + +netcdf_spatial_dimensions(::AbstractField{LX, LY, LZ}) where {LX, LY, LZ} = + tuple(xdim(LX)..., ydim(LY)..., zdim(LZ)...) const default_dimension_attributes = Dict( "xC" => Dict("longname" => "Locations of the cell centers in the x-direction.", "units" => "m"), @@ -34,38 +40,76 @@ const default_output_attributes = Dict( "S" => Dict("longname" => "Absolute salinity", "units" => "g/kg") ) +add_schedule_metadata!(attributes, schedule) = nothing + +function add_schedule_metadata!(global_attributes, schedule::IterationInterval) + global_attributes["schedule"] = "IterationInterval" + global_attributes["interval"] = schedule.interval + global_attributes["output iteration interval"] = + "Output was saved every $(schedule.interval) iteration(s)." + + return nothing +end + +function add_schedule_metadata!(global_attributes, schedule::TimeInterval) + global_attributes["schedule"] = "TimeInterval" + global_attributes["interval"] = schedule.interval + global_attributes["output time interval"] = + "Output was saved every $(prettytime(schedule.interval))." + + return nothing +end + +function add_schedule_metadata!(global_attributes, schedule::WallTimeInterval) + global_attributes["schedule"] = "WallTimeInterval" + global_attributes["interval"] = schedule.interval + global_attributes["output time interval"] = + "Output was saved every $(prettytime(schedule.interval))." + + return nothing +end + +function add_schedule_metadata!(global_attributes, schedule::AveragedTimeInterval) + add_schedule_metadata!(global_attributes, TimeInterval(schedule)) + + global_attributes["time_averaging_window"] = schedule.window + global_attributes["time averaging window"] = + "Output was time averaged with a window size of $(prettytime(schedule.window))" + + global_attributes["time_averaging_stride"] = schedule.stride + global_attributes["time averaging stride"] = + "Output was time averaged with a stride of $(schedule.stride) iteration(s) within the time averaging window." + + return nothing +end + """ NetCDFOutputWriter{D, O, I, T, S} <: AbstractOutputWriter An output writer for writing to NetCDF files. """ -mutable struct NetCDFOutputWriter{D, O, I, T, S, A} <: AbstractOutputWriter - filepath :: String - dataset :: D - outputs :: O - iteration_interval :: I - time_interval :: T - mode :: String - field_slicer :: S - array_type :: A - previous :: Float64 - verbose :: Bool +mutable struct NetCDFOutputWriter{D, O, T, S, A} <: AbstractOutputWriter + filepath :: String + dataset :: D + outputs :: O + schedule :: T + mode :: String + field_slicer :: S + array_type :: A + previous :: Float64 + verbose :: Bool end """ -function NetCDFOutputWriter(model, outputs; filepath, - iteration_interval = nothing, - time_interval = nothing, - time_averaging_window = nothing, - time_averaging_stride = 1, - array_type = Array{Float32}, - field_slicer = FieldSlicer(), - global_attributes = Dict(), - output_attributes = Dict(), - dimensions = Dict(), - mode = "c", - compression = 0, - verbose = false) +function NetCDFOutputWriter(model, outputs; filepath, schedule + array_type = Array{Float32}, + field_slicer = FieldSlicer(), + global_attributes = Dict(), + output_attributes = Dict(), + dimensions = Dict(), + mode = "c", + compression = 0, + verbose = false) Construct a `NetCDFOutputWriter` that writes `(label, output)` pairs in `outputs` (which should be a `Dict`) to a NetCDF file, where `label` is a string that labels the output and `output` is @@ -77,48 +121,32 @@ Keyword arguments ================= - `filepath` (required): Filepath to save output to. -- `iteration_interval`: Save output every `n` model iterations. - -- `time_interval`: Save output every `t` units of model clock time. - -- `time_averaging_window`: Specifies a time window over which each member of `output` is averaged before - being saved. For this each member of output is converted to `Oceananigans.Diagnostics.WindowedTimeAverage`. - Default `nothing` indicates no averaging. - -- `time_averaging_stride`: Specifies a iteration 'stride' between the calculation of each `output` during - time-averaging. Longer strides means that output is calculated less frequently, and that the resulting - time-average is faster to compute, but less accurate. Default: 1. +- `schedule` (required): `AbstractSchedule` that determines when output is saved. - `array_type`: The array type to which output arrays are converted to prior to saving. - Default: Array{Float32}. + Default: Array{Float32}. - `field_slicer`: An object for slicing field output in ``(x, y, z)``, including omitting halos. - Has no effect on output that is not a field. `field_slicer = nothing` means - no slicing occurs, so that all field data, including halo regions, is saved. - Default: FieldSlicer(), which slices halo regions. + Has no effect on output that is not a field. `field_slicer = nothing` means + no slicing occurs, so that all field data, including halo regions, is saved. + Default: `FieldSlicer()`, which slices halo regions. - `global_attributes`: Dict of model properties to save with every file (deafult: `Dict()`) - `output_attributes`: Dict of attributes to be saved with each field variable (reasonable - defaults are provided for velocities, buoyancy, temperature, and salinity). + defaults are provided for velocities, buoyancy, temperature, and salinity; + otherwise `output_attributes` *must* be user-provided). -- `dimensions`: A `Dict` of dimension tuples to apply to outputs (useful for function outputs - as field dimensions can be inferred). +- `dimensions`: A `Dict` of dimension tuples to apply to outputs (required for function outputs) - `with_halos`: Include the halo regions in the grid coordinates and output fields - (default: `false`). + (default: `false`). - `mode`: "a" (for append) and "c" (for clobber or create). Default: "c". See NCDatasets.jl - documentation for more information on the `mode` option. + documentation for more information on the `mode` option. - `compression`: Determines the compression level of data (0-9, default 0) -- `slice_kwargs`: `dimname = Union{OrdinalRange, Integer}` will slice the dimension `dimname`. - All other keywords are ignored. E.g. `xC = 3:10` will only produce output along the dimension - `xC` between indices 3 and 10 for all fields with `xC` as one of their dimensions. `xC = 1` - is treated like `xC = 1:1`. Multiple dimensions can be sliced in one call. Not providing slices - writes output over the entire domain (including halo regions if `with_halos=true`). - Examples ======== Saving the u velocity field and temperature fields, the full 3D fields and surface 2D slices @@ -136,26 +164,49 @@ simulation = Simulation(model, Δt=12, stop_time=3600); fields = Dict("u" => model.velocities.u, "T" => model.tracers.T); simulation.output_writers[:field_writer] = - NetCDFOutputWriter(model, fields, filepath="fields.nc", time_interval=60) + NetCDFOutputWriter(model, fields, filepath="fields.nc", schedule=TimeInterval(60)) # output -NetCDFOutputWriter (time_interval=60): fields.nc +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: fields.nc ├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0) -└── 2 outputs: ["T", "u"] +├── 2 outputs: ["T", "u"] +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +└── array type: Array{Float32} ``` ```jldoctest netcdf1 simulation.output_writers[:surface_slice_writer] = NetCDFOutputWriter(model, fields, filepath="surface_xy_slice.nc", - time_interval=60, field_slicer=FieldSlicer(k=grid.Nz)) + schedule=TimeInterval(60), field_slicer=FieldSlicer(k=grid.Nz)) # output -NetCDFOutputWriter (time_interval=60): surface_xy_slice.nc +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: surface_xy_slice.nc ├── dimensions: zC(1), zF(1), xC(16), yF(16), xF(16), yC(16), time(0) -└── 2 outputs: ["T", "u"] +├── 2 outputs: ["T", "u"] +├── field slicer: FieldSlicer(:, :, 16, with_halos=false) +└── array type: Array{Float32} ``` -Writing a scalar, profile, and slice to NetCDF: +```jldoctest netcdf1 +simulation.output_writers[:averaged_profile_writer] = + NetCDFOutputWriter(model, fields, + filepath = "averaged_z_profile.nc", + schedule = AveragedTimeInterval(60, window=20), + field_slicer = FieldSlicer(i=1, j=1)) + +# output +NetCDFOutputWriter scheduled on TimeInterval(1 minute): +├── filepath: averaged_z_profile.nc +├── dimensions: zC(16), zF(17), xC(1), yF(1), xF(1), yC(1), time(0) +├── 2 outputs: ["T", "u"] averaged on AveragedTimeInterval(window=20 seconds, stride=1, interval=1 minute) +├── field slicer: FieldSlicer(1, 1, :, with_halos=false) +└── array type: Array{Float32} +``` + +`NetCDFOutputWriter` also accepts output functions that write scalars and arrays to disk, +provided that their `dimensions` are provided: ```jldoctest using Oceananigans, Oceananigans.OutputWriters @@ -187,41 +238,27 @@ global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7); simulation.output_writers[:things] = NetCDFOutputWriter(model, outputs, - iteration_interval=1, filepath="things.nc", dimensions=dims, verbose=true, + schedule=IterationInterval(1), filepath="things.nc", dimensions=dims, verbose=true, global_attributes=global_attributes, output_attributes=output_attributes) # output -NetCDFOutputWriter (iteration_interval=1): things.nc +NetCDFOutputWriter scheduled on IterationInterval(1): +├── filepath: things.nc ├── dimensions: zC(16), zF(17), xC(16), yF(16), xF(16), yC(16), time(0) -└── 3 outputs: ["profile", "slice", "scalar"] +├── 3 outputs: ["profile", "slice", "scalar"] +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +└── array type: Array{Float32} ``` """ -function NetCDFOutputWriter(model, outputs; filepath, - iteration_interval = nothing, - time_interval = nothing, - time_averaging_window = nothing, - time_averaging_stride = 1, - array_type = Array{Float32}, - field_slicer = FieldSlicer(), - global_attributes = Dict(), - output_attributes = Dict(), - dimensions = Dict(), - mode = "c", - compression = 0, - verbose = false) - - validate_intervals(iteration_interval, time_interval) - - # Convert each output to WindowedTimeAverage if time_averaging_window is specified - if !isnothing(time_averaging_window) - - !isnothing(iteration_interval) && error("Cannot specify iteration_interval with time_averaging_window.") - - outputs = Dict(name => WindowedTimeAverage(outputs[name], model, time_interval = time_interval, - time_window = time_averaging_window, stride = time_averaging_stride, - field_slicer = field_slicer) - for name in keys(outputs)) - end +function NetCDFOutputWriter(model, outputs; filepath, schedule, + array_type = Array{Float32}, + field_slicer = FieldSlicer(), + global_attributes = Dict(), + output_attributes = Dict(), + dimensions = Dict(), + mode = "c", + compression = 0, + verbose = false) # Ensure we can add any kind of metadata to the global attributes later by converting to pairs of type {Any, Any}. global_attributes = Dict{Any, Any}(k => v for (k, v) in global_attributes) @@ -231,35 +268,12 @@ function NetCDFOutputWriter(model, outputs; filepath, global_attributes["Julia"] = "This file was generated using " * versioninfo_with_gpu() global_attributes["Oceananigans"] = "This file was generated using " * oceananigans_versioninfo() - # Add useful metadata about intervals and time averaging - if isnothing(time_averaging_window) - if !isnothing(iteration_interval) - global_attributes["iteration_interval"] = iteration_interval - global_attributes["output iteration interval"] = - "Output was saved every $iteration_interval iteration(s)." - end - - if !isnothing(time_interval) - global_attributes["time_interval"] = time_interval - global_attributes["output time interval"] = - "Output was saved every $(prettytime(time_interval))." - end - else - global_attributes["time_interval"] = time_interval - global_attributes["output time interval"] = - "Output was time averaged and saved every $(prettytime(time_interval))." - - global_attributes["time_averaging_window"] = time_averaging_window - global_attributes["time averaging window"] = - "Output was time averaged with a window size of $(prettytime(time_averaging_window))" - - if time_averaging_stride != 1 - global_attributes["time_averaging_stride"] = time_averaging_stride - global_attributes["time averaging stride"] = - "Output was time averaged with a stride of $time_averaging_stride iteration(s) within the time averaging window." - end - end + add_schedule_metadata!(global_attributes, schedule) + # Convert schedule to TimeInterval and each output to WindowedTimeAverage if + # schedule::AveragedTimeInterval + schedule, outputs = time_average_outputs(schedule, outputs, model, field_slicer) + grid = model.grid Nx, Ny, Nz = size(grid) Hx, Hy, Hz = halo_size(grid) @@ -296,29 +310,45 @@ function NetCDFOutputWriter(model, outputs; filepath, end end - # Initiates empty variables for fields from the user-supplied `outputs`. for (name, output) in outputs - if output isa Field - defVar(dataset, name, eltype(array_type), (netcdf_spatial_dimensions(output)..., "time"), - compression=compression, attrib=output_attributes[name]) - elseif output isa WindowedTimeAverage && output.operand isa Field - defVar(dataset, name, eltype(array_type), (netcdf_spatial_dimensions(output.operand)..., "time"), - compression=compression, attrib=output_attributes[name]) - else - name ∉ keys(dimensions) && error("Custom output $name needs dimensions!") - - defVar(dataset, name, eltype(array_type), (dimensions[name]..., "time"), - compression=compression, attrib=output_attributes[name]) - end + attributes = try output_attributes[name]; catch; Dict(); end + define_output_variable!(dataset, output, name, array_type, compression, attributes, dimensions) end sync(dataset) end - return NetCDFOutputWriter(filepath, dataset, outputs, iteration_interval, time_interval, - mode, field_slicer, array_type, 0.0, verbose) + return NetCDFOutputWriter(filepath, dataset, outputs, schedule, mode, field_slicer, array_type, 0.0, verbose) +end + +##### +##### Variable definition +##### + +""" Defines empty variables for 'custom' user-supplied `output`. """ +function define_output_variable!(dataset, output, name, array_type, compression, output_attributes, dimensions) + name ∉ keys(dimensions) && error("Custom output $name needs dimensions!") + + defVar(dataset, name, eltype(array_type), (dimensions[name]..., "time"), + compression=compression, attrib=output_attributes) + + return nothing end +""" Defines empty field variable. """ +define_output_variable!(dataset, output::AbstractField, name, array_type, compression, output_attributes, dimensions) = + defVar(dataset, name, eltype(array_type), + (netcdf_spatial_dimensions(output)..., "time"), + compression=compression, attrib=output_attributes) + +""" Defines empty field variable for `WindowedTimeAverage`s over fields. """ +define_output_variable!(dataset, output::WindowedTimeAverage{<:AbstractField}, args...) = + define_output_variable!(dataset, output.operand, args...) + +##### +##### Write output +##### + Base.open(ow::NetCDFOutputWriter) = Dataset(ow.filepath, "a") Base.close(ow::NetCDFOutputWriter) = close(ow.dataset) @@ -347,12 +377,7 @@ function write_output!(ow::NetCDFOutputWriter, model) verbose && (t0′ = time_ns()) data = fetch_and_convert_output(output, model, ow) - - if output isa AveragedField - data = dropdims(data, dims=output.dims) - elseif output isa WindowedTimeAverage && output.operand isa AveragedField - data = dropdims(data, dims=output.operand.dims) - end + data = drop_averaged_dims(output, data) colons = Tuple(Colon() for _ in 1:ndims(data)) ds[name][colons..., time_index] = data @@ -360,7 +385,7 @@ function write_output!(ow::NetCDFOutputWriter, model) if verbose # Time after computing this output. t1′ = time_ns() - @info "Computing $name done: time=$(prettytime((t1′-t0′)/1e9))" + @info "Computing $name done: time=$(prettytime((t1′-t0′) / 1e9))" end end @@ -378,16 +403,24 @@ function write_output!(ow::NetCDFOutputWriter, model) return nothing end -function Base.show(io::IO, ow::NetCDFOutputWriter) - freq_int = isnothing(ow.iteration_interval) && isnothing(ow.time_interval) ? "" : - isnothing(ow.iteration_interval) ? "(time_interval=$(ow.time_interval))" : - isnothing(ow.time_interval) ? "(iteration_interval=$(ow.iteration_interval))" : - "(iteration_interval=$(ow.iteration_interval), time_interval=$(ow.time_interval))" +drop_averaged_dims(output, data) = data # fallback +drop_averaged_dims(output::AveragedField, data) = dropdims(data, dims=output.dims) +drop_averaged_dims(output::WindowedTimeAverage{<:AveragedField}, data) = dropdims(data, dims=output.operand.dims) +##### +##### Show +##### + +function Base.show(io::IO, ow::NetCDFOutputWriter) dims = join([dim * "(" * string(length(ow.dataset[dim])) * "), " for dim in keys(ow.dataset.dim)])[1:end-2] - print(io, "NetCDFOutputWriter $freq_int: $(ow.filepath)\n", - "├── dimensions: $dims\n", - "└── $(length(ow.outputs)) outputs: $(keys(ow.outputs))") + averaging_schedule = output_averaging_schedule(ow) + + print(io, "NetCDFOutputWriter scheduled on $(show_schedule(ow.schedule)):", '\n', + "├── filepath: $(ow.filepath)", '\n', + "├── dimensions: $dims", '\n', + "├── $(length(ow.outputs)) outputs: $(keys(ow.outputs))", show_averaging_schedule(averaging_schedule), '\n', + "├── field slicer: $(short_show(ow.field_slicer))", '\n', + "└── array type: ", show_array_type(ow.array_type)) end diff --git a/src/OutputWriters/output_writer_utils.jl b/src/OutputWriters/output_writer_utils.jl index 1f8cd70bab..fa3d1204eb 100644 --- a/src/OutputWriters/output_writer_utils.jl +++ b/src/OutputWriters/output_writer_utils.jl @@ -106,3 +106,13 @@ function has_reference(has_type, obj) return typeof(obj) <: has_type end end + +""" Returns the schedule for output averaging determined by the first output value. """ +function output_averaging_schedule(ow::AbstractOutputWriter) + first_output = first(values(ow.outputs)) + return output_averaging_schedule(first_output) +end + +output_averaging_schedule(output) = nothing # fallback + +show_array_type(a::Type{Array{T}}) where T = "Array{$T}" diff --git a/src/OutputWriters/time_average_outputs.jl b/src/OutputWriters/time_average_outputs.jl new file mode 100644 index 0000000000..8a6c057222 --- /dev/null +++ b/src/OutputWriters/time_average_outputs.jl @@ -0,0 +1,27 @@ +time_average_outputs(schedule, outputs, model, field_slicer) = schedule, outputs # fallback + +""" + time_average_outputs(schedule::AveragedTimeInterval, outputs, model, field_slicer) + +Wrap each `output` in a `WindowedTimeAverage` on the time-averaged `schedule` and with `field_slicer`. + +Returns the `TimeInterval` associated with `schedule` and a `NamedTuple` or `Dict` of the wrapped +outputs. +""" +function time_average_outputs(schedule::AveragedTimeInterval, outputs::Dict, model, field_slicer) + averaged_outputs = Dict(name => WindowedTimeAverage(output, model; schedule=schedule, field_slicer=field_slicer) + for (name, output) in outputs) + + return TimeInterval(schedule), averaged_outputs +end + +function time_average_outputs(schedule::AveragedTimeInterval, outputs::NamedTuple, model, field_slicer) + output_names = Tuple(keys(outputs)) + + output_values = Tuple(WindowedTimeAverage(outputs[name], model; schedule=schedule, field_slicer=field_slicer) + for name in output_names) + + averaged_outputs = NamedTuple{output_names}(output_values) + + return TimeInterval(schedule), averaged_outputs +end diff --git a/src/OutputWriters/windowed_time_average.jl b/src/OutputWriters/windowed_time_average.jl index 3acab11794..4bc5af16a9 100644 --- a/src/OutputWriters/windowed_time_average.jl +++ b/src/OutputWriters/windowed_time_average.jl @@ -1,68 +1,125 @@ using Oceananigans.Diagnostics: AbstractDiagnostic +using Oceananigans.OutputWriters: fetch_output +using Oceananigans.Utils: AbstractSchedule, prettytime -import Oceananigans.Utils: time_to_run -import Oceananigans.Fields: AbstractField, compute! -import Oceananigans.OutputWriters: fetch_output +import Oceananigans.Utils: TimeInterval, show_schedule import Oceananigans.Diagnostics: run_diagnostic! """ - WindowedTimeAverage{R, OP, FS} <: AbstractDiagnostic + mutable struct AveragedTimeInterval <: AbstractSchedule -An object for computing 'windowed' time averages, or moving time-averages -of a `operand` over a specified `time_window`, collected on `time_interval`. +Container for parameters that configure and handle time-averaged output. """ -mutable struct WindowedTimeAverage{R, OP, FS} <: AbstractDiagnostic - result :: R - operand :: OP - time_window :: Float64 - time_interval :: Float64 +mutable struct AveragedTimeInterval <: AbstractSchedule + interval :: Float64 + window :: Float64 stride :: Int - field_slicer :: FS - window_start_time :: Float64 - window_start_iteration :: Int - previous_collection_time :: Float64 previous_interval_stop_time :: Float64 collecting :: Bool end """ - WindowedTimeAverage(operand; time_window, - time_interval, - stride = 1, - field_slicer = FieldSlicer()) + AveragedTimeInterval(interval; window=interval, stride=1) + +Returns a `schedule` that specifies periodic time-averaging of output. +The time `window` specifies the extent of the time-average, which +reoccurs every `interval`. + +`output` is computed and accumulated into the average every `stride` iterations +during the averaging window. For example, `stride=1` computs output every iteration, +whereas `stride=2` computes output every other iteration. Time-averages with +longer `stride`s are faster to compute, but less accurate. + +The time-average of ``a`` is a left Riemann sum corresponding to + +`` ⟨a⟩ = 1/T \\int_{tᵢ-T}^T a \\mathrm{d} t ,`` + +where ``⟨a⟩`` is the time-average of ``a``, ``T`` is the time-window for averaging, +and the ``tᵢ`` are discrete times separated by the time `interval`. The ``tᵢ`` specify +both the end of the averaging window and the time at which output is written. + +Example +======= + +```jldoctest averaged_time_interval +using Oceananigans.OutputWriters: AveragedTimeInterval +using Oceananigans.Utils: year, years + +schedule = AveragedTimeInterval(4years, window=1year) + +# output +AveragedTimeInterval(window=1 year, stride=1, interval=4 years) +``` + +An `AveragedTimeInterval` schedule directs an output writer +to time-average its outputs before writing them to disk: + +```jldoctest averaged_time_interval +using Oceananigans +using Oceananigans.OutputWriters: JLD2OutputWriter +using Oceananigans.Utils: minutes + +model = IncompressibleModel(grid=RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))) + +simulation = Simulation(model, Δt=10minutes, stop_time=30years) + +simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, + prefix = "averaged_velocity_data", + schedule = AveragedTimeInterval(4years, window=1year, stride=2)) + +# output +JLD2OutputWriter scheduled on TimeInterval(4 years): +├── filepath: ./averaged_velocity_data.jld2 +├── 3 outputs: (:u, :v, :w) averaged on AveragedTimeInterval(window=1 year, stride=2, interval=4 years) +├── field slicer: FieldSlicer(:, :, :, with_halos=false) +├── array type: Array{Float32} +├── including: [:grid, :coriolis, :buoyancy, :closure] +└── max filesize: Inf YiB +``` +""" +AveragedTimeInterval(interval; window=interval, stride=1) = + AveragedTimeInterval(Float64(interval), Float64(window), stride, 0.0, false) + +# Determines whether or not to call run_diagnostic +(schedule::AveragedTimeInterval)(model) = + schedule.collecting || model.clock.time >= schedule.previous_interval_stop_time + schedule.interval - schedule.window + +TimeInterval(schedule::AveragedTimeInterval) = TimeInterval(schedule.interval) + +""" + WindowedTimeAverage{OP, R, FS} <: AbstractDiagnostic + +An object for computing 'windowed' time averages, or moving time-averages +of a `operand` over a specified `window`, collected on `interval`. +""" +mutable struct WindowedTimeAverage{OP, R, FS} <: AbstractDiagnostic + result :: R + operand :: OP + window_start_time :: Float64 + window_start_iteration :: Int + previous_collection_time :: Float64 + field_slicer :: FS + schedule :: AveragedTimeInterval +end + +""" + WindowedTimeAverage(operand, model=nothing; schedule, field_slicer=FieldSlicer()) -Returns an object for computing running averages of `operand` over `time_window`, -recurring on `time_interval`. During the collection period, averages are computed -every `stride` iteration. +Returns an object for computing running averages of `operand` over `schedule.window` and +recurring on `schedule.interval`, where `schedule` is an `AveragedTimeInterval`. +During the collection period, averages are computed every `schedule.stride` iteration. -`operand` may be a `Oceananigans.Field` or a subtype of it. +`operand` may be a `Oceananigans.Field` or a function that returns an array or scalar. Calling `wta(model)` for `wta::WindowedTimeAverage` object returns `wta.result`. """ -function WindowedTimeAverage(operand, model=nothing; time_window, time_interval, stride=1, - field_slicer = FieldSlicer()) - +function WindowedTimeAverage(operand, model=nothing; schedule, field_slicer=FieldSlicer()) + output = fetch_output(operand, model, field_slicer) result = similar(output) # convert views to arrays result .= output # initialize `result` with initial output - return WindowedTimeAverage(result, - operand, - Float64(time_window), - Float64(time_interval), - stride, - field_slicer, - 0.0, 0, 0.0, 0.0, false) -end - -function time_to_run(clock, wta::WindowedTimeAverage) - if (wta.collecting || - clock.time >= wta.previous_interval_stop_time + wta.time_interval - wta.time_window) - - return true - else - return false - end + return WindowedTimeAverage(result, operand, 0.0, 0, 0.0, field_slicer, schedule) end function accumulate_result!(wta, model) @@ -76,6 +133,7 @@ function accumulate_result!(wta, model) # Accumulate left Riemann sum integrand = fetch_output(wta.operand, model, wta.field_slicer) + @. wta.result = (wta.result * T_previous + integrand * Δt) / T_current # Save time of integrand collection @@ -87,24 +145,28 @@ end function run_diagnostic!(wta::WindowedTimeAverage, model) if model.clock.iteration == 0 # initialize previous interval stop time - wta.previous_interval_stop_time = model.clock.time + wta.schedule.previous_interval_stop_time = model.clock.time end # Don't start collecting if we are *only* "initializing" run_diagnostic! at the beginning # of a Simulation. # - # Note: this can be false at the zeroth iteration if time_interval == time_window (which + # Note: this can be false at the zeroth iteration if interval == window (which # implies we are always collecting) - initializing = model.clock.iteration == 0 && - model.clock.time < wta.previous_interval_stop_time + wta.time_interval - wta.time_window - - if !(wta.collecting) && !(initializing) - # run_diagnostic! has been called, but we are not currently collecting data. + unscheduled = model.clock.iteration == 0 && + model.clock.time < wta.schedule.previous_interval_stop_time + wta.schedule.interval - wta.schedule.window + + if unscheduled + # This is an "unscheduled" call to run_diagnostic! --- which occurs when run_diagnostic! + # is called at the beginning of a run (and schedule.interval != schedule.window). + # In this case we do nothing. + elseif !(wta.schedule.collecting) + # run_diagnostic! has been called on schedule, but we are not currently collecting data. # Initialize data collection: # Start averaging period - wta.collecting = true + wta.schedule.collecting = true # Zero out result wta.result .= 0 @@ -114,18 +176,17 @@ function run_diagnostic!(wta::WindowedTimeAverage, model) wta.window_start_iteration = model.clock.iteration wta.previous_collection_time = model.clock.time - elseif model.clock.time >= wta.previous_interval_stop_time + wta.time_interval + elseif model.clock.time >= wta.schedule.previous_interval_stop_time + wta.schedule.interval # Output is imminent. Finalize averages and cease data collection. accumulate_result!(wta, model) # Averaging period is complete. - wta.collecting = false + wta.schedule.collecting = false - # Reset the "previous" interval time, - # subtracting a sliver that presents window overshoot from accumulating. - wta.previous_interval_stop_time = model.clock.time - rem(model.clock.time, wta.time_interval) + # Reset the "previous" interval time, subtracting a sliver that presents overshoot from accumulating. + wta.schedule.previous_interval_stop_time = model.clock.time - rem(model.clock.time, wta.schedule.interval) - elseif mod(model.clock.iteration - wta.window_start_iteration, wta.stride) == 0 + elseif mod(model.clock.iteration - wta.window_start_iteration, wta.schedule.stride) == 0 # Collect data as usual accumulate_result!(wta, model) end @@ -136,9 +197,21 @@ end function (wta::WindowedTimeAverage)(model) # For the paranoid - wta.collecting && + wta.schedule.collecting && model.clock.iteration > 0 && @warn "Returning a WindowedTimeAverage before the collection period is complete." return wta.result end + +Base.show(io::IO, schedule::AveragedTimeInterval) = print(io, short_show(schedule)) + +short_show(schedule::AveragedTimeInterval) = string("AveragedTimeInterval(", + "window=", prettytime(schedule.window), ", ", + "stride=", schedule.stride, ", ", + "interval=", prettytime(schedule.interval), ")") + +show_averaging_schedule(schedule) = "" +show_averaging_schedule(schedule::AveragedTimeInterval) = string(" averaged on ", short_show(schedule)) + +output_averaging_schedule(output::WindowedTimeAverage) = output.schedule diff --git a/src/Simulations/run.jl b/src/Simulations/run.jl index 516e585ae1..bf90fc4d66 100644 --- a/src/Simulations/run.jl +++ b/src/Simulations/run.jl @@ -1,3 +1,4 @@ +using Oceananigans.Utils: initialize_schedule! using Oceananigans.OutputWriters: WindowedTimeAverage using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper, RungeKutta3TimeStepper @@ -70,13 +71,19 @@ function run!(sim) model = sim.model clock = model.clock - [open(writer) for writer in values(sim.output_writers)] + # Initialization + for writer in values(sim.output_writers) + open(writer) + initialize_schedule!(writer.schedule) + add_dependencies!(sim.diagnostics, writer) + end - [add_dependencies!(sim.diagnostics, writer) for writer in values(sim.output_writers)] + [initialize_schedule!(diag.schedule) for diag in values(sim.diagnostics)] while !stop(sim) time_before = time() + # Evaluate all diagnostics and write output at first iteration if clock.iteration == 0 [run_diagnostic!(diag, sim.model) for diag in values(sim.diagnostics)] [write_output!(out, sim.model) for out in values(sim.output_writers)] @@ -86,8 +93,8 @@ function run!(sim) euler = clock.iteration == 0 || (sim.Δt isa TimeStepWizard && n == 1) ab2_or_rk3_time_step!(model, get_Δt(sim.Δt), euler=euler) - [time_to_run(clock, diag) && run_diagnostic!(diag, sim.model) for diag in values(sim.diagnostics)] - [time_to_run(clock, out) && write_output!(out, sim.model) for out in values(sim.output_writers)] + [ diag.schedule(model) && run_diagnostic!(diag, sim.model) for diag in values(sim.diagnostics) ] + [ writer.schedule(model) && write_output!(writer, sim.model) for writer in values(sim.output_writers) ] end sim.progress(sim) diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 7ef76b58eb..633a8558fb 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -43,5 +43,6 @@ include("output_writer_diagnostic_utils.jl") include("ordered_dict_show.jl") include("with_tracers.jl") include("versioninfo.jl") +include("schedules.jl") end diff --git a/src/Utils/output_writer_diagnostic_utils.jl b/src/Utils/output_writer_diagnostic_utils.jl index 8528c5a5d8..4ee0941afe 100644 --- a/src/Utils/output_writer_diagnostic_utils.jl +++ b/src/Utils/output_writer_diagnostic_utils.jl @@ -27,39 +27,3 @@ function push!(container::DiagOrWriterDict, elems...) end return nothing end - -""" - validate_intervals(iteration_interval, time_interval) - -Warn the user if iteration and time intervals are both `nothing`. -""" -function validate_intervals(iteration_interval, time_interval) - isnothing(iteration_interval) && isnothing(time_interval) && @warn "You have not specified any intervals." - return nothing -end - -has_iteration_interval(obj) = :iteration_interval in propertynames(obj) && obj.iteration_interval != nothing -has_time_interval(obj) = :time_interval in propertynames(obj) && obj.time_interval != nothing - -iteration_interval_is_ripe(clock, obj) = has_iteration_interval(obj) && clock.iteration % obj.iteration_interval == 0 - -function time_interval_is_ripe(clock, obj) - if has_time_interval(obj) && clock.time >= obj.previous + obj.time_interval - obj.previous = clock.time - rem(clock.time, obj.time_interval) - return true - else - return false - end -end - -function time_to_run(clock, diag) - iteration_interval_is_ripe(clock, diag) && return true - time_interval_is_ripe(clock, diag) && return true - - # If the output writer does not specify any intervals, - # it is unable to write output and we throw an error as a convenience. - has_iteration_interval(diag) || has_time_interval(diag) || - error("$(typeof(diag)) must have an iteration or time interval specified!") - - return false -end diff --git a/src/Utils/schedules.jl b/src/Utils/schedules.jl new file mode 100644 index 0000000000..13a123d86a --- /dev/null +++ b/src/Utils/schedules.jl @@ -0,0 +1,100 @@ +""" + AbstractSchedule + +Supertype for objects that schedule `OutputWriter`s and `Diagnostics`. +Schedules must define a function `Schedule(model)` that returns true or +false. +""" +abstract type AbstractSchedule end + +initialize_schedule!(abstract_schedule) = nothing # fallback + +##### +##### TimeInterval +##### + +""" + struct TimeInterval <: AbstractSchedule + +Callable `TimeInterval` schedule for periodic output or diagnostic evaluation +according to `model.clock.time`. +""" +mutable struct TimeInterval <: AbstractSchedule + interval :: Float64 + previous_actuation_time :: Float64 +end + +""" + TimeInterval(interval) + +Returns a callable `TimeInterval` that schedules periodic output or diagnostic evaluation +on a `interval` of simulation time, as kept by `model.clock`. +""" +TimeInterval(interval) = TimeInterval(Float64(interval), 0.0) + +function (schedule::TimeInterval)(model) + time = model.clock.time + + if time >= schedule.previous_actuation_time + schedule.interval + # Shave overshoot off previous_actuation_time to prevent overshoot from accumulating + schedule.previous_actuation_time = time - rem(time, schedule.interval) + return true + else + return false + end + +end + +##### +##### IterationInterval +##### + +""" + IterationInterval(interval) + +Returns a callable IterationInterval that schedules periodic output or diagnostic evaluation +over `interval`, and therefore when `model.clock.iteration % schedule.interval == 0`. +""" +struct IterationInterval <: AbstractSchedule + interval :: Int +end + +(schedule::IterationInterval)(model) = model.clock.iteration % schedule.interval == 0 + +""" + mutable struct WallTimeInterval <: AbstractSchedule + +""" +mutable struct WallTimeInterval <: AbstractSchedule + interval :: Float64 + previous_actuation_time :: Float64 +end + +""" + WallTimeInterval(interval) + +Returns a callable WallTimeInterval that schedules periodic output or diagnostic evaluation +on a `interval` of "wall time" while a simulation runs. The "wall time" +is the actual real world time, as kept by an actual or hypothetical clock hanging +on your wall. +""" +WallTimeInterval(interval) = WallTimeInterval(Float64(interval), time_ns() * 1e-9) + +initialize_schedule!(schedule::WallTimeInterval) = schedule.previous_actuation_time = time_ns() * 1e-9 + +function (schedule::WallTimeInterval)(model) + wall_time = time_ns() * 1e-9 + + if wall_time >= schedule.previous_actuation_time + schedule.interval + # Shave overshoot off previous_actuation_time to prevent overshoot from accumulating + schedule.previous_actuation_time = wall_time - rem(wall_time, schedule.interval) + return true + else + return false + end + +end + +show_schedule(schedule) = string(schedule) +show_schedule(schedule::IterationInterval) = string("IterationInterval(", schedule.interval, ")") +show_schedule(schedule::TimeInterval) = string("TimeInterval(", prettytime(schedule.interval), ")") diff --git a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl index 0c8b74b082..6931182dca 100644 --- a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl +++ b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl @@ -53,7 +53,7 @@ function run_ocean_large_eddy_simulation_regression_test(arch, closure) simulation.stop_iteration = spinup_steps-test_steps run!(simulation) - checkpointer = Checkpointer(model, iteration_interval = test_steps, prefix = name, + checkpointer = Checkpointer(model, schedule = IterationInterval(test_steps), prefix = name, dir = joinpath(dirname(@__FILE__), "data")) simulation.output_writers[:checkpointer] = checkpointer diff --git a/test/regression_tests/rayleigh_benard_regression_test.jl b/test/regression_tests/rayleigh_benard_regression_test.jl index e02939ed85..15e2a0ba4c 100644 --- a/test/regression_tests/rayleigh_benard_regression_test.jl +++ b/test/regression_tests/rayleigh_benard_regression_test.jl @@ -56,7 +56,7 @@ function run_rayleigh_benard_regression_test(arch) prefix = "rayleigh_benard" - checkpointer = Checkpointer(model, iteration_interval=test_steps, prefix=prefix, + checkpointer = Checkpointer(model, schedule=IterationInterval(test_steps), prefix=prefix, dir=joinpath(dirname(@__FILE__), "data")) ##### diff --git a/test/regression_tests/thermal_bubble_regression_test.jl b/test/regression_tests/thermal_bubble_regression_test.jl index e9d58296fe..c8285a98eb 100644 --- a/test/regression_tests/thermal_bubble_regression_test.jl +++ b/test/regression_tests/thermal_bubble_regression_test.jl @@ -33,7 +33,7 @@ function run_thermal_bubble_regression_test(arch) "T" => model.tracers.T, "S" => model.tracers.S) - nc_writer = NetCDFOutputWriter(model, outputs, filename=regression_data_filepath, iteration_interval=10) + nc_writer = NetCDFOutputWriter(model, outputs, filename=regression_data_filepath, schedule=IterationInterval(10)) push!(simulation.output_writers, nc_writer) =# diff --git a/test/runtests.jl b/test/runtests.jl index ea2998b537..df42cb6c81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,7 +70,7 @@ closures = ( CUDA.allowscalar(true) -include("runtests_utils.jl") +include("utils_for_runtests.jl") group = get(ENV, "TEST_GROUP", :all) |> Symbol diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl index 14eace41cd..c71f888572 100644 --- a/test/test_diagnostics.jl +++ b/test/test_diagnostics.jl @@ -1,9 +1,11 @@ +using Oceananigans.Diagnostics + function nan_checker_aborts_simulation(arch, FT) grid = RegularCartesianGrid(size=(4, 2, 1), extent=(1, 1, 1)) model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT) # It checks for NaNs in w by default. - nc = NaNChecker(model; iteration_interval=1, fields=Dict(:w => model.velocities.w.data.parent)) + nc = NaNChecker(model; schedule=IterationInterval(1), fields=Dict(:w => model.velocities.w.data.parent)) push!(model.diagnostics, nc) model.velocities.w[3, 2, 1] = NaN @@ -63,48 +65,26 @@ end get_iteration(model) = model.clock.iteration get_time(model) = model.clock.time -function timeseries_diagnostic_works(arch, FT) - model = TestModel(arch, FT) - Δt = FT(1e-16) - simulation = Simulation(model, Δt=Δt, stop_iteration=1) - iter_diag = TimeSeries(get_iteration, model, iteration_interval=1) - push!(simulation.diagnostics, iter_diag) - run!(simulation) - return iter_diag.time[end] == Δt && iter_diag.data[end] == 1 -end - -function timeseries_diagnostic_tuples(arch, FT) - model = TestModel(arch, FT) - Δt = FT(1e-16) - simulation = Simulation(model, Δt=Δt, stop_iteration=2) - timeseries = TimeSeries((iters=get_iteration, itertimes=get_time), model, iteration_interval=2) - simulation.diagnostics[:timeseries] = timeseries - run!(simulation) - return timeseries.iters[end] == 2 && timeseries.itertimes[end] == 2Δt -end - function diagnostics_getindex(arch, FT) model = TestModel(arch, FT) simulation = Simulation(model, Δt=0, stop_iteration=0) - iter_timeseries = TimeSeries(get_iteration, model) - time_timeseries = TimeSeries(get_time, model) - simulation.diagnostics[:iters] = iter_timeseries - simulation.diagnostics[:times] = time_timeseries - return simulation.diagnostics[2] == time_timeseries + nc = NaNChecker(model; schedule=IterationInterval(1), fields=Dict(:w => model.velocities.w.data.parent)) + simulation.diagnostics[:nc] = nc + return simulation.diagnostics[1] == nc end function diagnostics_setindex(arch, FT) model = TestModel(arch, FT) simulation = Simulation(model, Δt=0, stop_iteration=0) - iter_timeseries = TimeSeries(get_iteration, model) - time_timeseries = TimeSeries(get_time, model) - max_abs_u_timeseries = TimeSeries(FieldMaximum(abs, model.velocities.u), model, iteration_interval=1) + nc1 = NaNChecker(model; schedule=IterationInterval(1), fields=Dict(:w => model.velocities.w.data.parent)) + nc2 = NaNChecker(model; schedule=IterationInterval(2), fields=Dict(:u => model.velocities.u.data.parent)) + nc3 = NaNChecker(model; schedule=IterationInterval(3), fields=Dict(:v => model.velocities.v.data.parent)) - push!(simulation.diagnostics, iter_timeseries, time_timeseries) - simulation.diagnostics[2] = max_abs_u_timeseries + push!(simulation.diagnostics, nc1, nc2) + simulation.diagnostics[2] = nc3 - return simulation.diagnostics[:diag2] == max_abs_u_timeseries + return simulation.diagnostics[:diag2] == nc3 end @testset "Diagnostics" begin @@ -126,8 +106,6 @@ end @test diffusive_cfl_diagnostic_is_correct(arch, FT) @test advective_cfl_diagnostic_is_correct(arch, FT) @test max_abs_field_diagnostic_is_correct(arch, FT) - @test timeseries_diagnostic_works(arch, FT) - @test timeseries_diagnostic_tuples(arch, FT) @test diagnostics_getindex(arch, FT) @test diagnostics_setindex(arch, FT) end diff --git a/test/test_output_writers.jl b/test/test_output_writers.jl index df6eaa346f..acb13ae432 100644 --- a/test/test_output_writers.jl +++ b/test/test_output_writers.jl @@ -3,6 +3,7 @@ using NCDatasets using Oceananigans.BoundaryConditions: BoundaryFunction, PBC, FBC, ZFBC using Oceananigans.Diagnostics using Oceananigans.Fields +using Oceananigans.OutputWriters function instantiate_windowed_time_average(model) @@ -13,7 +14,7 @@ function instantiate_windowed_time_average(model) u₀ = similar(interior(u)) u₀ .= interior(u) - wta = WindowedTimeAverage(model.velocities.u, time_window=1.0, time_interval=10.0) + wta = WindowedTimeAverage(model.velocities.u, schedule=AveragedTimeInterval(10, window=1)) return all(wta(model) .== u₀) end @@ -25,7 +26,7 @@ function time_step_with_windowed_time_average(model) set!(model, u=0, v=0, w=0, T=0, S=0) - wta = WindowedTimeAverage(model.velocities.u, time_window=2.0, time_interval=4.0) + wta = WindowedTimeAverage(model.velocities.u, schedule=AveragedTimeInterval(4, window=2)) simulation = Simulation(model, Δt=1.0, stop_time=4.0) simulation.diagnostics[:u_avg] = wta @@ -62,7 +63,7 @@ function run_thermal_bubble_netcdf_tests(arch) "S" => model.tracers.S) nc_filepath = "test_dump_$(typeof(arch)).nc" - nc_writer = NetCDFOutputWriter(model, outputs, filepath=nc_filepath, iteration_interval=10, verbose=true) + nc_writer = NetCDFOutputWriter(model, outputs, filepath=nc_filepath, schedule=IterationInterval(10), verbose=true) push!(simulation.output_writers, nc_writer) i_slice = 1:10 @@ -72,7 +73,7 @@ function run_thermal_bubble_netcdf_tests(arch) j_slice = j_slice:j_slice # So we can correctly index with it for later tests. nc_sliced_filepath = "test_dump_sliced_$(typeof(arch)).nc" - nc_sliced_writer = NetCDFOutputWriter(model, outputs, filepath=nc_sliced_filepath, iteration_interval=10, + nc_sliced_writer = NetCDFOutputWriter(model, outputs, filepath=nc_sliced_filepath, schedule=IterationInterval(10), field_slicer=field_slicer, verbose=true) push!(simulation.output_writers, nc_sliced_writer) @@ -84,7 +85,8 @@ function run_thermal_bubble_netcdf_tests(arch) @test haskey(ds3.attrib, "date") && !isnothing(ds3.attrib["date"]) @test haskey(ds3.attrib, "Julia") && !isnothing(ds3.attrib["Julia"]) @test haskey(ds3.attrib, "Oceananigans") && !isnothing(ds3.attrib["Oceananigans"]) - @test haskey(ds3.attrib, "iteration_interval") && ds3.attrib["iteration_interval"] == 10 + @test haskey(ds3.attrib, "schedule") && ds3.attrib["schedule"] == "IterationInterval" + @test haskey(ds3.attrib, "interval") && ds3.attrib["interval"] == 10 @test haskey(ds3.attrib, "output iteration interval") && !isnothing(ds3.attrib["output iteration interval"]) @test eltype(ds3["time"]) == eltype(model.clock.time) @@ -142,7 +144,8 @@ function run_thermal_bubble_netcdf_tests(arch) @test haskey(ds2.attrib, "date") && !isnothing(ds2.attrib["date"]) @test haskey(ds2.attrib, "Julia") && !isnothing(ds2.attrib["Julia"]) @test haskey(ds2.attrib, "Oceananigans") && !isnothing(ds2.attrib["Oceananigans"]) - @test haskey(ds2.attrib, "iteration_interval") && ds2.attrib["iteration_interval"] == 10 + @test haskey(ds2.attrib, "schedule") && ds2.attrib["schedule"] == "IterationInterval" + @test haskey(ds2.attrib, "interval") && ds2.attrib["interval"] == 10 @test haskey(ds2.attrib, "output iteration interval") && !isnothing(ds2.attrib["output iteration interval"]) @test eltype(ds2["time"]) == eltype(model.clock.time) @@ -222,7 +225,7 @@ function run_thermal_bubble_netcdf_tests_with_halos(arch) ) nc_filepath = "test_dump_with_halos_$(typeof(arch)).nc" - nc_writer = NetCDFOutputWriter(model, outputs, filepath=nc_filepath, iteration_interval=10, + nc_writer = NetCDFOutputWriter(model, outputs, filepath=nc_filepath, schedule=IterationInterval(10), field_slicer=FieldSlicer(with_halos=true)) push!(simulation.output_writers, nc_writer) @@ -233,7 +236,8 @@ function run_thermal_bubble_netcdf_tests_with_halos(arch) @test haskey(ds.attrib, "date") && !isnothing(ds.attrib["date"]) @test haskey(ds.attrib, "Julia") && !isnothing(ds.attrib["Julia"]) @test haskey(ds.attrib, "Oceananigans") && !isnothing(ds.attrib["Oceananigans"]) - @test haskey(ds.attrib, "iteration_interval") && ds.attrib["iteration_interval"] == 10 + @test haskey(ds.attrib, "schedule") && ds.attrib["schedule"] == "IterationInterval" + @test haskey(ds.attrib, "interval") && ds.attrib["interval"] == 10 @test haskey(ds.attrib, "output iteration interval") && !isnothing(ds.attrib["output iteration interval"]) @test eltype(ds["time"]) == eltype(model.clock.time) @@ -319,13 +323,11 @@ function run_netcdf_function_output_tests(arch) global_attributes = Dict("location" => "Bay of Fundy", "onions" => 7) nc_filepath = "test_function_outputs_$(typeof(arch)).nc" + simulation.output_writers[:food] = NetCDFOutputWriter(model, outputs; filepath=nc_filepath, - time_interval=Δt, dimensions=dims, array_type=Array{Float64}, verbose=true, - global_attributes=global_attributes, output_attributes=output_attributes) - - # We should have no problem with time_interval=Δt=1.25 as 1.25 can be represented exactly with - # a floating-point number. + schedule=TimeInterval(Δt), dimensions=dims, array_type=Array{Float64}, verbose=true, + global_attributes=global_attributes, output_attributes=output_attributes) run!(simulation) @@ -334,7 +336,8 @@ function run_netcdf_function_output_tests(arch) @test haskey(ds.attrib, "date") && !isnothing(ds.attrib["date"]) @test haskey(ds.attrib, "Julia") && !isnothing(ds.attrib["Julia"]) @test haskey(ds.attrib, "Oceananigans") && !isnothing(ds.attrib["Oceananigans"]) - @test haskey(ds.attrib, "time_interval") && !isnothing(ds.attrib["time_interval"]) + @test haskey(ds.attrib, "schedule") && !isnothing(ds.attrib["schedule"]) + @test haskey(ds.attrib, "interval") && !isnothing(ds.attrib["interval"]) @test haskey(ds.attrib, "output time interval") && !isnothing(ds.attrib["output time interval"]) @test eltype(ds["time"]) == eltype(model.clock.time) @@ -413,8 +416,8 @@ function run_netcdf_function_output_tests(arch) simulation.output_writers[:food] = NetCDFOutputWriter(model, outputs; filepath=nc_filepath, mode="a", - iteration_interval=1, array_type=Array{Float64}, dimensions=dims, verbose=true, - global_attributes=global_attributes, output_attributes=output_attributes) + schedule=IterationInterval(1), array_type=Array{Float64}, dimensions=dims, verbose=true, + global_attributes=global_attributes, output_attributes=output_attributes) run!(simulation) @@ -457,7 +460,7 @@ function jld2_field_output(model) simulation = Simulation(model, Δt=1.0, stop_iteration=1) simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, - time_interval = 1.0, + schedule = TimeInterval(1), dir = ".", prefix = "test", force = true) @@ -497,7 +500,7 @@ function jld2_sliced_field_output(model) simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities, - time_interval = 1.0, + schedule = TimeInterval(1), field_slicer = FieldSlicer(i=1:2, j=1:3, k=:), dir = ".", prefix = "test", @@ -526,7 +529,7 @@ function run_jld2_file_splitting_tests(arch) file["boundary_conditions/fake"] = π end - ow = JLD2OutputWriter(model, (u=model.velocities.u,); dir=".", prefix="test", iteration_interval=1, + ow = JLD2OutputWriter(model, (u=model.velocities.u,); dir=".", prefix="test", schedule=IterationInterval(1), init=fake_bc_init, including=[:grid], field_slicer=nothing, array_type=Array{Float64}, max_filesize=200KiB, force=true) @@ -587,7 +590,7 @@ function run_thermal_bubble_checkpointer_tests(arch) run!(true_simulation) checkpointed_simulation = Simulation(checkpointed_model, Δt=Δt, stop_iteration=5) - checkpointer = Checkpointer(checkpointed_model, iteration_interval=5, force=true) + checkpointer = Checkpointer(checkpointed_model, schedule=IterationInterval(5), force=true) push!(checkpointed_simulation.output_writers, checkpointer) # Checkpoint should be saved as "checkpoint5.jld" after the 5th iteration. @@ -635,7 +638,7 @@ function run_checkpoint_with_function_bcs_tests(arch) model = IncompressibleModel(architecture=arch, grid=grid, boundary_conditions=(u=u_bcs, T=T_bcs)) set!(model, u=π/2, v=ℯ, T=Base.MathConstants.γ, S=Base.MathConstants.φ) - checkpointer = Checkpointer(model) + checkpointer = Checkpointer(model, schedule=IterationInterval(1)) write_output!(checkpointer, model) model = nothing @@ -705,7 +708,7 @@ function run_cross_architecture_checkpointer_tests(arch1, arch2) model = IncompressibleModel(architecture=arch1, grid=grid) set!(model, u=π/2, v=ℯ, T=Base.MathConstants.γ, S=Base.MathConstants.φ) - checkpointer = Checkpointer(model) + checkpointer = Checkpointer(model, schedule=IterationInterval(1)) write_output!(checkpointer, model) model = nothing @@ -753,20 +756,20 @@ end function run_dependency_adding_tests(model) - windowed_time_average = WindowedTimeAverage(model.velocities.u, time_window=2.0, time_interval=4.0) + windowed_time_average = WindowedTimeAverage(model.velocities.u, schedule=AveragedTimeInterval(4, window=2)) output = Dict("time_average" => windowed_time_average) attributes = Dict("time_average" => Dict("longname" => "A time average", "units" => "arbitrary")) dimensions = Dict("time_average" => ("xF", "yC", "zC")) # JLD2 dependencies test - jld2_output_writer = JLD2OutputWriter(model, output, time_interval=4.0, dir=".", prefix="test", force=true) + jld2_output_writer = JLD2OutputWriter(model, output, schedule=TimeInterval(4), dir=".", prefix="test", force=true) @test dependencies_added_correctly!(model, windowed_time_average, jld2_output_writer) # NetCDF dependency test netcdf_output_writer = NetCDFOutputWriter(model, output, - time_interval = 4.0, + schedule = TimeInterval(4), filepath = "test.nc", output_attributes = attributes, dimensions = dimensions) @@ -785,16 +788,15 @@ function run_windowed_time_averaging_simulation_tests!(model) simulation = Simulation(model, Δt=1.0, stop_iteration=0) jld2_output_writer = JLD2OutputWriter(model, model.velocities, - time_interval = π, - time_averaging_window = 1.0, - prefix = "test", - force = true) + schedule = AveragedTimeInterval(π, window=1), + prefix = "test", + force = true) + # https://github.com/Alexander-Barth/NCDatasets.jl/issues/105 nc_filepath = "windowed_time_average_test1.nc" nc_outputs = Dict(string(name) => field for (name, field) in pairs(model.velocities)) nc_output_writer = NetCDFOutputWriter(model, nc_outputs, filepath=nc_filepath, - # https://github.com/Alexander-Barth/NCDatasets.jl/issues/105 - time_interval = Float64(π), time_averaging_window = 1.0) + schedule = AveragedTimeInterval(π, window=1)) jld2_outputs_are_time_averaged = Tuple(typeof(out) <: WindowedTimeAverage for out in jld2_output_writer.outputs) nc_outputs_are_time_averaged = Tuple(typeof(out) <: WindowedTimeAverage for out in values(nc_output_writer.outputs)) @@ -812,8 +814,8 @@ function run_windowed_time_averaging_simulation_tests!(model) jld2_u_windowed_time_average = simulation.output_writers[:jld2].outputs.u nc_w_windowed_time_average = simulation.output_writers[:nc].outputs["w"] - @test !(jld2_u_windowed_time_average.collecting) - @test !(nc_w_windowed_time_average.collecting) + @test !(jld2_u_windowed_time_average.schedule.collecting) + @test !(nc_w_windowed_time_average.schedule.collecting) # Test that time-averaging is finalized prior to output even when averaging over # time_window is not fully realized. For this, step forward to a time at which @@ -822,40 +824,38 @@ function run_windowed_time_averaging_simulation_tests!(model) simulation.stop_iteration = 2 run!(simulation) # model.clock.time = 3.0, just before output but after average-collection. - @test jld2_u_windowed_time_average.collecting - @test nc_w_windowed_time_average.collecting + @test jld2_u_windowed_time_average.schedule.collecting + @test nc_w_windowed_time_average.schedule.collecting # Step forward such that time_window is not reached, but output will occur. simulation.Δt = π - 3 + 0.01 # ≈ 0.15 < 1.0 simulation.stop_iteration = 3 run!(simulation) # model.clock.time ≈ 3.15, after output - @test jld2_u_windowed_time_average.previous_interval_stop_time == - model.clock.time - rem(model.clock.time, jld2_u_windowed_time_average.time_interval) + @test jld2_u_windowed_time_average.schedule.previous_interval_stop_time == + model.clock.time - rem(model.clock.time, jld2_u_windowed_time_average.schedule.interval) - @test nc_w_windowed_time_average.previous_interval_stop_time == - model.clock.time - rem(model.clock.time, nc_w_windowed_time_average.time_interval) + @test nc_w_windowed_time_average.schedule.previous_interval_stop_time == + model.clock.time - rem(model.clock.time, nc_w_windowed_time_average.schedule.interval) # Test that collection does start when a simulation is initialized and # time_interval == time_averaging_window model.clock.iteration = model.clock.time = 0 simulation.output_writers[:jld2] = JLD2OutputWriter(model, model.velocities, - time_interval = π, - time_averaging_window = π, - prefix = "test", - force = true) + schedule = AveragedTimeInterval(π, window=π), + prefix = "test", + force = true) nc_filepath = "windowed_time_average_test2.nc" nc_outputs = Dict(string(name) => field for (name, field) in pairs(model.velocities)) simulation.output_writers[:nc] = NetCDFOutputWriter(model, nc_outputs, filepath=nc_filepath, - # https://github.com/Alexander-Barth/NCDatasets.jl/issues/105 - time_interval = Float64(π), time_averaging_window = Float64(π)) + schedule=AveragedTimeInterval(π, window=π)) run!(simulation) - @test simulation.output_writers[:jld2].outputs.u.collecting - @test simulation.output_writers[:nc].outputs["w"].collecting + @test simulation.output_writers[:jld2].outputs.u.schedule.collecting + @test simulation.output_writers[:nc].outputs["w"].schedule.collecting return nothing end @@ -870,6 +870,9 @@ function jld2_time_averaging_of_horizontal_averages(model) w = (x, y, z) -> 0, T = (x, y, z) -> 4) + u, v, w = model.velocities + T, S = model.tracers + simulation = Simulation(model, Δt=1.0, stop_iteration=5) u, v, w = model.velocities @@ -879,18 +882,17 @@ function jld2_time_averaging_of_horizontal_averages(model) uv = AveragedField(u * v, dims=(1, 2)), wT = AveragedField(w * T, dims=(1, 2))) - simulation.output_writers[:velocities] = JLD2OutputWriter(model, average_fluxes, - time_interval = 4.0, - time_averaging_window = 2.0, - dir = ".", - prefix = "test", - force = true) + simulation.output_writers[:fluxes] = JLD2OutputWriter(model, average_fluxes, + schedule = AveragedTimeInterval(4, window=2), + dir = ".", + prefix = "test", + force = true) run!(simulation) file = jldopen("test.jld2") - # Data is saved with halos by default + # Data is saved without halos by default wu = file["timeseries/wu/4"][1, 1, 3] uv = file["timeseries/uv/4"][1, 1, 3] wT = file["timeseries/wT/4"][1, 1, 3] @@ -937,16 +939,16 @@ function run_netcdf_time_averaging_tests(arch) horizontal_average_nc_filepath = "decay_averaged_field_test.nc" simulation.output_writers[:horizontal_average] = - NetCDFOutputWriter(model, nc_outputs, filepath=horizontal_average_nc_filepath, time_interval=10Δt, + NetCDFOutputWriter(model, nc_outputs, filepath=horizontal_average_nc_filepath, schedule=TimeInterval(10Δt), dimensions=nc_dimensions, array_type=Array{Float64}, verbose=true) time_average_nc_filepath = "decay_windowed_time_average_test.nc" window = 6Δt stride = 2 simulation.output_writers[:time_average] = - NetCDFOutputWriter(model, nc_outputs, filepath=time_average_nc_filepath, array_type=Array{Float64}, - time_interval=10Δt, time_averaging_window=window, time_averaging_stride=stride, - dimensions=nc_dimensions, verbose=true) + NetCDFOutputWriter(model, nc_outputs, filepath=time_average_nc_filepath, array_type=Array{Float64}, + schedule=AveragedTimeInterval(10Δt, window=window, stride=stride), + dimensions=nc_dimensions, verbose=true) run!(simulation) @@ -979,10 +981,9 @@ function run_netcdf_time_averaging_tests(arch) ds = NCDataset(time_average_nc_filepath) - attribute_names = ( - "time_interval", "output time interval", - "time_averaging_window", "time averaging window", - "time_averaging_stride", "time averaging stride") + attribute_names = ("schedule", "interval", "output time interval", + "time_averaging_window", "time averaging window", + "time_averaging_stride", "time averaging stride") for name in attribute_names @test haskey(ds.attrib, name) && !isnothing(ds.attrib[name]) @@ -1015,7 +1016,9 @@ end for arch in archs # Some tests can reuse this same grid and model. - grid = RegularCartesianGrid(size=(4, 4, 4), extent=(1, 1, 1)) + grid = RegularCartesianGrid(size=(4, 4, 4), extent=(1, 1, 1), + topology=(Periodic, Periodic, Bounded)) + model = IncompressibleModel(architecture=arch, grid=grid) @testset "WindowedTimeAverage" begin diff --git a/test/runtests_utils.jl b/test/utils_for_runtests.jl similarity index 100% rename from test/runtests_utils.jl rename to test/utils_for_runtests.jl diff --git a/verification/convergence_tests/ConvergenceTests/DoublyPeriodicTaylorGreen.jl b/verification/convergence_tests/ConvergenceTests/DoublyPeriodicTaylorGreen.jl index 0b0a8d9373..c38a289849 100644 --- a/verification/convergence_tests/ConvergenceTests/DoublyPeriodicTaylorGreen.jl +++ b/verification/convergence_tests/ConvergenceTests/DoublyPeriodicTaylorGreen.jl @@ -43,7 +43,7 @@ function setup_simulation(; Nx, Δt, stop_iteration, U=1, architecture=CPU(), di simulation.output_writers[:fields] = JLD2OutputWriter(model, model.velocities; dir = dir, force = true, field_slicer = nothing, prefix = @sprintf("taylor_green_Nx%d_Δt%.1e", Nx, Δt), - time_interval = stop_iteration * Δt / 10) + schedule = TimeInterval(stop_iteration * Δt / 10)) return simulation end diff --git a/verification/convergence_tests/ConvergenceTests/ForcedFlowFixedSlip.jl b/verification/convergence_tests/ConvergenceTests/ForcedFlowFixedSlip.jl index 88cdffcadc..a03086af46 100644 --- a/verification/convergence_tests/ConvergenceTests/ForcedFlowFixedSlip.jl +++ b/verification/convergence_tests/ConvergenceTests/ForcedFlowFixedSlip.jl @@ -73,7 +73,7 @@ function setup_xy_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= dir = dir, force = true, prefix = @sprintf("forced_fixed_slip_xy_Nx%d_Δt%.1e", Nx, Δt), field_slicer = nothing, - time_interval = stop_iteration * Δt / 2) + schedule = TimeInterval(stop_iteration * Δt / 2)) return simulation end @@ -119,7 +119,7 @@ function setup_xz_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= dir = dir, force = true, prefix = @sprintf("forced_fixed_slip_xz_Nx%d_Δt%.1e", Nx, Δt), field_slicer = nothing, - time_interval = stop_iteration * Δt / 2) + schedule = TimeInterval(stop_iteration * Δt / 2)) return simulation end diff --git a/verification/convergence_tests/ConvergenceTests/ForcedFlowFreeSlip.jl b/verification/convergence_tests/ConvergenceTests/ForcedFlowFreeSlip.jl index f11d4967f1..9ea1dfe4ef 100644 --- a/verification/convergence_tests/ConvergenceTests/ForcedFlowFreeSlip.jl +++ b/verification/convergence_tests/ConvergenceTests/ForcedFlowFreeSlip.jl @@ -59,7 +59,7 @@ function setup_xz_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= dir = dir, force = true, field_slicer = nothing, prefix = @sprintf("forced_free_slip_xz_Nx%d_Δt%.1e", Nx, Δt), - time_interval = stop_iteration * Δt / 2) + schedule = TimeInterval(stop_iteration * Δt / 2)) return simulation end @@ -99,7 +99,7 @@ function setup_xy_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir= dir = dir, force = true, field_slicer = nothing, prefix = @sprintf("forced_free_slip_xy_Nx%d_Δt%.1e", Nx, Δt), - time_interval=stop_iteration * Δt / 2) + schedule = TimeInterval(stop_iteration * Δt / 2)) return simulation end diff --git a/verification/convergence_tests/ConvergenceTests/OneDimensionalCosineAdvectionDiffusion.jl b/verification/convergence_tests/ConvergenceTests/OneDimensionalCosineAdvectionDiffusion.jl index 90d2c5138b..7a482bee3a 100644 --- a/verification/convergence_tests/ConvergenceTests/OneDimensionalCosineAdvectionDiffusion.jl +++ b/verification/convergence_tests/ConvergenceTests/OneDimensionalCosineAdvectionDiffusion.jl @@ -71,7 +71,7 @@ function run_test(; Nx, Δt, stop_iteration, U = 1, κ = 1e-4, u = (x, y, z) -> c(y, x, z, 0, U, κ), c = (x, y, z) -> c(y, x, z, 0, U, κ)) - simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, iteration_interval=stop_iteration) + simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, schedule=IterationInterval(stop_iteration)) @info "Running 1D in y cosine advection diffusion test for u and cy with Ny = $Nx and Δt = $Δt ($(typeof(advection)))..." run!(simulation) diff --git a/verification/convergence_tests/ConvergenceTests/TwoDimensionalDiffusion.jl b/verification/convergence_tests/ConvergenceTests/TwoDimensionalDiffusion.jl index 26e9fe4fca..39f73ee903 100644 --- a/verification/convergence_tests/ConvergenceTests/TwoDimensionalDiffusion.jl +++ b/verification/convergence_tests/ConvergenceTests/TwoDimensionalDiffusion.jl @@ -47,7 +47,7 @@ function setup_simulation(; Nx, Δt, stop_iteration, architecture=CPU(), dir=DAT simulation.output_writers[:fields] = JLD2OutputWriter(model, model.tracers; dir = dir, force = true, field_slicer = nothing, prefix = @sprintf("%s_%s_diffusion_Nx%d_Δt%.1e", "$(topo[1])", "$(topo[2])", Nx, Δt), - time_interval = stop_iteration * Δt / 10) + schedule = TimeInterval(stop_iteration * Δt / 10)) end return simulation diff --git a/verification/lid_driven_cavity/lid_driven_cavity.jl b/verification/lid_driven_cavity/lid_driven_cavity.jl index 2a1654a900..c6a0c7febe 100644 --- a/verification/lid_driven_cavity/lid_driven_cavity.jl +++ b/verification/lid_driven_cavity/lid_driven_cavity.jl @@ -49,7 +49,7 @@ function simulate_lid_driven_cavity(; Re, N, end_time) output_attributes = Dict("ζ" => Dict("longname" => "vorticity", "units" => "1/s")) field_output_writer = - NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", time_interval=0.1, + NetCDFOutputWriter(model, fields, filename="lid_driven_cavity_Re$Re.nc", schedule=TimeInterval(0.1), global_attributes=global_attributes, output_attributes=output_attributes, dimensions=dims) diff --git a/verification/stommel_gyre/stommel_gyre_advection.jl b/verification/stommel_gyre/stommel_gyre_advection.jl index 021f54572f..54c8c091ad 100644 --- a/verification/stommel_gyre/stommel_gyre_advection.jl +++ b/verification/stommel_gyre/stommel_gyre_advection.jl @@ -59,7 +59,7 @@ function setup_simulation(N, T, CFL, ϕₐ, advection_scheme; u, v) output_attributes = Dict("c" => Dict("longname" => "passive tracer")) simulation.output_writers[:fields] = - NetCDFOutputWriter(model, fields, filename=filename, time_interval=0.01, + NetCDFOutputWriter(model, fields, filename=filename, schedule=TimeInterval(0.01), global_attributes=global_attributes, output_attributes=output_attributes) return simulation diff --git a/verification/stratified_couette_flow/stratified_couette_flow.jl b/verification/stratified_couette_flow/stratified_couette_flow.jl index d1afffe147..e4cdad81b8 100644 --- a/verification/stratified_couette_flow/stratified_couette_flow.jl +++ b/verification/stratified_couette_flow/stratified_couette_flow.jl @@ -185,7 +185,7 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1, field_writer = JLD2OutputWriter(model, fields, dir=base_dir, prefix=prefix * "_fields", - init=init_save_parameters_and_bcs, time_interval=10, + init=init_save_parameters_and_bcs, schedule=TimeInterval(10), force=true, verbose=true) ##### @@ -209,7 +209,7 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1, profile_writer = JLD2OutputWriter(model, profiles, dir=base_dir, prefix=prefix * "_profiles", - init=init_save_parameters_and_bcs, time_interval=1, + init=init_save_parameters_and_bcs, schedule=TimeInterval(1), force=true, verbose=true) ##### @@ -225,7 +225,7 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1, statistics_writer = JLD2OutputWriter(model, statistics, dir=base_dir, prefix=prefix * "_statistics", - init=init_save_parameters_and_bcs, time_interval=1/2, + init=init_save_parameters_and_bcs, schedule=TimeInterval(1/2), force=true, verbose=true) ##### diff --git a/verification/thermal_bubble/thermal_bubble.jl b/verification/thermal_bubble/thermal_bubble.jl index 845bac0b38..76ac76d841 100644 --- a/verification/thermal_bubble/thermal_bubble.jl +++ b/verification/thermal_bubble/thermal_bubble.jl @@ -49,7 +49,7 @@ function setup_simulation(N, advection_scheme) global_attributes = Dict("N" => N, "advection_scheme" => string(typeof(advection_scheme))) simulation.output_writers[:fields] = - NetCDFOutputWriter(model, fields, filename=filename, time_interval=100, global_attributes=global_attributes) + NetCDFOutputWriter(model, fields, filename=filename, schedule=TimeInterval(100), global_attributes=global_attributes) return simulation end