Data API

Our data is principally four dimensional: $x$, $y$, time and facies. Nevertheless, in our naming we drop the facies dimension, so a full data set is reffered to as a DataVolume and a section is reffered to as a DataSlice.

Principally, we'd like to be able to do post-processing independent from knowing the input. In more complex scenarios (like variable subsidence) this poses a bit of a challenge on the data being stored.

Data Hierarchy

---
title: Concept hierarchy
---
flowchart LR
    Bundle --> Pack --> Data --> DataHeader
    Bundle --> Header --> Axes
    Pack --> Header

Most functions for obtaining reduced or transformed data are defined on the level of a Pack.

Header

We store the different axes in a Vector.

@kwdef struct Axes
    x::Vector{Length}
    y::Vector{Length}
    t::Vector{Time}
end

The Header struct contains an index of all data sets contained in a bundle.

@kwdef struct DataHeader
    kind::Symbol
    slice::Slice2
    write_interval::Int
end

@kwdef struct Header
    tag::String
    axes::Axes

    Δt::Time
    time_steps::Int
    grid_size::NTuple{2,Int}
    n_facies::Int

    initial_topography::Matrix{Amount}
    sea_level::Vector{Length}
    subsidence_rate::Rate
    data_sets::Dict{Symbol,DataHeader}
    attributes::Dict{String,Any} = Dict()
end

data_sets(header::Header) = header.data_sets
data_volumes(header::Header) = [s for (s, h) in header.data_sets if h.kind == :volume]
data_slices(header::Header) = [s for (s, h) in header.data_sets if h.kind == :slice]
data_columns(header::Header) = [s for (s, h) in header.data_sets if h.kind == :column]

data_sets(bundle::AbstractBundle) = bundle |> header |> data_sets
data_volumes(bundle::AbstractBundle) = bundle |> header |> data_volumes
data_slices(bundle::AbstractBundle) = bundle |> header |> data_slices
data_columns(bundle::AbstractBundle) = bundle |> header |> data_columns

Data

The principle task for the Data struct is to hold the simulation data as written for each time step (modulo write interval). We also store the slice and write_interval so that the dimensions of a Data object are understood with respect to the full run.

In this definition, D is the dimensionality of the data, 1 for column, 2 for slice and 3 for volume. Because many members also have a facies dimension, F should always be D+1.

@kwdef struct Data{F,D}
    slice::Slice2
    write_interval::Int
    # Julia doesn't allow to say Array{Amount,D+1} here
    disintegration::Array{Amount,F}
    production::Array{Amount,F}
    deposition::Array{Amount,F}
    bathymetry::Array{Amount,D}
    active_layer::Union{Array{Amount,F}, Nothing} = nothing
    stratigraphic_column::Ref{Union{Array{Amount,F}, Nothing}} = nothing
end

const DataVolume = Data{4,3}
const DataSlice = Data{3,2}
const DataColumn = Data{2,1}

disintegration(v::Data) = v.disintegration
production(v::Data) = v.production
deposition(v::Data) = v.deposition
active_layer(v::Data) = v.active_layer
bathymetry(v::Data) = v.bathymetry

Note, that the stratigraphic_column field is there for caching computations. We make this field mutable by using a Ref. We can further slice data using the index operator.

Base.getindex(v::Data{F,D}, args...) where {F,D} =
    let k = count_ints(args...)
        Data{F - k,D - k}(
            reduce_slice(v.slice, args...),
            v.write_interval,
            v.disintegration[:, args..., :],
            v.production[:, args..., :],
            v.deposition[:, args..., :],
            v.bathymetry[args..., :],
            v.active_layer == nothing ? nothing : v.active_layer[:, args..., :],
            nothing)  # stratigraphic_column: reset so it is recomputed for the slice
    end

Pack

A Pack is the combination of a Header and Data field.

struct Pack{F, D}
    header::Header
    data::Data{F, D}
end

bathymetry(v::Pack) = bathymetry(v.data)
production(v::Pack) = production(v.data)
deposition(v::Pack) = deposition(v.data)
active_layer(v::Pack) = active_layer(v.data)
disintegration(v::Pack) = disintegration(v.data)

Base.getindex(p::Pack, args...) = Pack(p.header, getindex(p.data, args...))

Bundle

A Bundle combines multiple data sets. It should support listing all data sets, and obtaining a header and obtaining a data set by key-value association.

Slices

In the Input struct the user can specify a dictionary of OutputSpec, specifying how much and at which interval to write output. Typically, you'd want a full topographic output with lower time resolution, and choose a transect with full time resolution. For example, on the ALCAP model:

const INPUT = ALCAP.Input(
    box = Box{Coast}(grid_size = (300, 150), phys_scale = 50.0u"m"),
    time = TimeProperties(Δt = 50.0u"yr", steps = 20_000),
    output = Dict(
        :topography => OutputSpec(write_interval = 200),
        :profile => OutputSpec(slice = (:, 75))),

    ...)  # add more options

Saving the full output of this simulation would take several hundreds of gigabytes, not gargantuan, but a bit unwieldy if you want to save many simulation runs. With this output specification, we cut down on this significantly.

@kwdef struct Input <: AbstractInput
    output = Dict(:full => OutputSpec((:,:), 1))
end

The default is to write all output, which is fine for smaller runs. The slice argument of OutputSpec will accept three different forms:

• (:, :) (default) output the full area of the model.
• (<n>, :) or (:, <n>) output a slice of the model, either with a fixed $x$ or a fixed $y$ coordinate. In our examples we always have the $x$ axis orthogonal to the shoreline, so slicing with a fixed $y$ (the second form) is what we use.
• (<m>, <n>), output a column of the model. If you have a very precise experiment workflow, this could be of use. You'll have to specify each column as a separate output. Most of the time though, we can extract columns from slice data in a post-processing stage, so if all your columns have the same $y$ coordinate, taking a slice is the preferred option.

If an output with a write_interval that does not divide the total number of timesteps is specified, the output array will contain div(time.steps, write_interval) values and the value for the remainder of the steps at the end of the model will not be recorded (since this will not cover an equal interval of time as the rest of the values).

Slices refer to the way the full DataVolume is sliced into smaller DataSlice and DataColumn objects. If a slice is (:, :) we are dealing with the full volume. If one of the members of the pair is an integer, we have a data slice and with two integers we are dealing with a column. We have some helper functions to deal with these distinctions:

count_ints(::Int, args...) = 1 + count_ints(args...)
count_ints(_, args...) = count_ints(args...)
count_ints() = 0

reduce_slice(s::Tuple{Colon,Colon}, x, y) = (x, y)
reduce_slice(s::Tuple{Int,Colon}, y::Int) = (s[1], y)
reduce_slice(s::Tuple{Colon,Int}, x::Int) = (x, s[2])

The slice information is heterogenous, so it is stored as a string in the HDF5 file. We need to be able to parse it back.

function parse_slice(s::AbstractString)
    if s == ":"
        return (:)
    end

    elements = split(s, ":")
    if length(elements) == 1
        return parse(Int, s)
    end

    a, b = elements
    return parse(Int, a):parse(Int, b)
end

parse_multi_slice(s::AbstractString) = Slice2(parse_slice.(split(s, ",")))

data_kind(::Int, ::Int) = :column
data_kind(::Int, _) = :slice
data_kind(_, ::Int) = :slice
data_kind(_, _) = :volume
data_kind(spec::OutputSpec) = data_kind(spec.slice...)

Module

module Storage

using Unitful
using .Iterators: repeated

import ...CarboKitten: time_axis, box_axes, OutputSpec, AbstractOutput, AbstractInput, AbstractState
import ..Abstract: AbstractBundle, header

using ...Algorithms.StratigraphicColumn: stratigraphic_column!

export Header, Data, DataVolume, DataSlice, DataColumn, Pack, Axes, DataHeader
export data_kind, data_sets, data_volumes, data_slices, data_columns
export sediment_thickness, water_depth, surface_heights, stratigraphic_column

const Length = typeof(1.0u"m")
const Time = typeof(1.0u"Myr")
const Slice2 = NTuple{2,Union{Int,Colon,UnitRange{Int}}}
const Amount = typeof(1.0u"m")
const Sediment = typeof(1.0u"m")
const Rate = typeof(1.0u"m/Myr")

<<slice-helpers>>
<<data-header>>
<<data-data>>
<<data-pack>>

"""
    stratigraphic_column(data)

Given a data set, compute the stratigraphic column. Result is memoised in
`data.stratigraphic_column` so repeated calls are free.
"""
function stratigraphic_column(data::Data{F, D}) where {F, D}
    if data.stratigraphic_column[] === nothing
        net_deposition = data.deposition .- data.disintegration
        for c in eachslice(net_deposition, dims=(1:D...,))
            stratigraphic_column!(c)
        end
        data.stratigraphic_column[] = net_deposition
    end
    return data.stratigraphic_column[]
end

stratigraphic_column(p::Pack{F, D}) where {F, D} = stratigraphic_column(p.data)

"""
    water_depth(header, data)

Compute the water depth function for the given data set.
"""
function water_depth(header::Header, data::Data{F, D}) where {F, D}
    sl = reshape(header.sea_level[1:data.write_interval:end], (repeated(1, D-1)..., :))
    return sl .- data.bathymetry
end

water_depth(p::Pack{F, D}) where {F, D} = water_depth(p.header, p.data)

"""
    sediment_thickness(data)

Compute the sediment thickness at each moment in the run by taking the cumulative
sum of the net deposition (deposition - disintegration) at each moment.
"""
function sediment_thickness(data::Data{F, D}) where {F, D}
    net_deposition = dropdims(sum(data.deposition .- data.disintegration, dims=1), dims=1)
    for c in eachslice(net_deposition, dims=(1:D-1...,))
        for i in 2:length(c)
            c[i] += c[i-1]
        end
    end
    return net_deposition
end

sediment_thickness(p::Pack{F, D}) where {F, D} = sediment_thickness(p.data)

"""
    surface_heights(header, data)

Compute the sediment surface height at every `(spatial..., time)` cell,
accounting for subsidence and net deposition. Returns an array of shape
`(spatial..., n_t+1)` where the first time entry is the initial topography
minus total subsidence and subsequent entries accumulate the preserved
sediment column.

Works with `DataColumn`, `DataSlice`, and `DataVolume`.
"""
function surface_heights(header::Header, data::Data{F, D}) where {F, D}
    total_subsidence = (header.axes.t[end] - header.axes.t[1]) * header.subsidence_rate
    initial_topography = header.initial_topography[data.slice...]
    sc = stratigraphic_column(data)
    # Sum over the facies dimension (dim 1), yielding (spatial..., n_t)
    sc_sum = dropdims(sum(sc, dims=1), dims=1)
    # Cumulative sediment accumulation along the time axis (last dim)
    accumulated = cumsum(sc_sum, dims=ndims(sc_sum))
    n_t = size(sc_sum, ndims(sc_sum))
    h0 = initial_topography .- total_subsidence
    # Build result array: shape (spatial..., n_t+1)
    sz = (size(sc_sum)[1:end-1]..., n_t + 1)
    h = Array{eltype(h0)}(undef, sz...)
    selectdim(h, ndims(h), 1) .= h0
    selectdim(h, ndims(h), 2:n_t+1) .= h0 .+ accumulated
    return h
end

surface_heights(p::Pack{F, D}) where {F, D} = surface_heights(p.header, p.data)

end