SINDBAD

julia

Simulation

The Simulation module provides the core functionality for running the SINDBAD Terrestrial Ecosystem Model (TEM). It includes utilities for preparing model-ready objects, managing spinup processes and running models.

Purpose

This module integrates various components and utilities required to execute the SINDBAD TEM, including precomputations, spinup, and time loop simulations. It supports parallel execution and efficient handling of large datasets.

Dependencies

Related (SINDBAD ecosystem)

TimeSamplers: Temporal helpers used in time-loop workflows.
OmniTools: Shared utilities used during simulation and output handling.

External (third-party)

ComponentArrays: Hierarchical state/parameter containers.
ProgressMeter: Progress bars for long-running simulations.
ThreadPools: Threaded parallelism helpers.

Internal (within Sindbad)

Sindbad.DataLoaders
Sindbad.Setup
Sindbad.Types
SindbadTEM

Included Files

utilsSimulation.jl: Core helpers for forcing slices, output handling, progress/logging, and simulation orchestration.
deriveSpinupForcing.jl: Derive spinup forcing data for steady-state initialization.
prepTEMOut.jl: Prepare output containers/structures for efficient writing.
prepTEM.jl: Prepare model-ready inputs and configuration for runs.
runTEMForLocation.jl: Run the TEM for a single location (optional spinup + main loop).
runTEMInSpace.jl: Run the TEM across spatial grids (parallel execution).
utilsLazy.jl: Helpers of lazy runs on YAXArrays cubes: the wrapper that hides cube types, the cache and block sizing, and the zarr store helpers.
prepTEMLazy.jl: Build the lazy run of the TEM on YAXArrays cubes with xmap (prepTEMLazy).
runTEMLazy.jl: Run the lazy run block by block into a zarr store (runTEMLazy), and the TEM of each grid in it (TEMLazy).
spinupTEM.jl: Spinup routines to initialize the model to steady state.
spinupSequence.jl: Sequential spinup loops for iterative refinement.

Notes

The package is designed to be modular and extensible, allowing users to customize and extend its functionality for specific use cases.
It integrates tightly with the SINDBAD framework, leveraging shared types and utilities from Setup.
Some spinup solvers are enabled via an optional extension:
- NLsolve → SindbadNLsolveExt (see ext/SindbadNLsolveExt/).

source

Functions

TEMLazy

Sindbad.Simulation.TEMLazy Function

julia

TEMLazy(map_cubes; loc_land::NamedTuple, tem::NamedTuple, selected_models, forcing_vars::AbstractArray)

Arguments:

map_cubes: collection/tuple of all input and output cubes from mapCube
loc_land: initial SINDBAD land with all fields and subfields
tem: a nested NT with necessary information of helpers, models, and spinup needed to run SINDBAD TEM and models
selected_models: a tuple of all models selected in the given model structure
forcing_vars: forcing variables
forcing_sources: the index of the input cube of each forcing variable, see shareForcingCubes. Variables with the same source share one cube.
parameter_to_index: the index of each parameter per location in the models, or nothing when there are none. The values of the location then set the parameters of selected_models before the run.
parameter_time: the getTimeVaryingInfo of the parameters, or nothing when no parameter varies in time. The models of the location are then a TimeVaryingModels, which the time loop updates.

source

Code

julia

function TEMLazy(map_cubes...;selected_models, forcing_vars, forcing_sources=Tuple(eachindex(forcing_vars)), loc_land::NamedTuple, output_vars, tem::NamedTuple, clean_data, restart_vars=(), parameter_to_index=nothing, parameter_time=nothing)
    # Make NaN check here instead of AllNaN filters


    outputs, inputs, restarts, parameter_inputs = unpackGridCubes(map_cubes; output_vars, forcing_vars, restart_vars, n_inputs=maximum(forcing_sources))
    selected_models = setGridParameters(selected_models, parameter_to_index, parameter_inputs)
    selected_models = setTimeVaryingModels(selected_models, parameter_inputs, parameter_time)
    # the same rule as filterNanGrids: a grid with an invalid forcing value at any
    # time step is not run, and its output and restart stay NaN
    if any(in_forcing -> any(is_invalid_number, in_forcing), inputs)
        foreach(xout -> fill!(xout, NaN), outputs)
        foreach(xout -> fill!(xout, NaN), restarts)
        return nothing
    end

    # ? apply clean_data fields to input data points
    _num_type, _forcing_data_info = clean_data

    # each variable is cleaned with its own units and bounds, also when it shares its
    # input cube with another variable
    inputs = map(enumerate(forcing_sources)) do (i, i_source)
        _data_info = _forcing_data_info[i]
        map(data_point -> cleanData(data_point, _data_info, _num_type), inputs[i_source])
    end
    loc_forcing = (; Pair.(forcing_vars, inputs)...)
    # the restart cubes are filled after the spinup step with `save_restart`
    store_restart = land -> fillGridRestart(restarts, land, restart_vars)
    # every pixel gets its own copy, as the land has mutable arrays and the pixels can
    # run on different threads
    coreTEMLazy!(outputs, selected_models, loc_forcing, deepcopy(loc_land), tem, store_restart)
    return nothing
end

function TEMLazyForGrid(map_cubes...; tem_kwargs::Base.RefValue{Any})
    TEMLazy(map_cubes...; tem_kwargs[]...)
    return nothing
end

accumulateSpinupForcing

Sindbad.Simulation.accumulateSpinupForcing Function

julia

accumulateSpinupForcing(spin_sequences, forcing, tem_helpers, spinup_forcing)

Derive the spinup forcing for each sequence in spin_sequences and collect them into a NamedTuple keyed by the forcing name, deriving each distinct name only once.

Arguments:

spin_sequences: a tuple of spinup sequences still to be processed
forcing: a forcing NT that contains the forcing time series set for a location
tem_helpers: helper NT with necessary objects for model run and type consistencies
spinup_forcing: the NamedTuple of forcings accumulated so far

Notes:

The sequences are consumed recursively with Base.tail rather than with a for loop, so that each sequence is handled by its own method specialisation. The forcing name is read from the sequence's type parameter, which lets the "already derived" check resolve at compile time.

source

Code

julia

function accumulateSpinupForcing end

accumulateSpinupForcing(::Tuple{}, _, _, spinup_forcing) = spinup_forcing

accumulateSpinupForcing(spin_sequences::Tuple{SpinupStepWithAggregator{:none},Vararg}, forcing, tem_helpers, spinup_forcing) = accumulateSpinupForcing(Base.tail(spin_sequences), forcing, tem_helpers, spinup_forcing)

function accumulateSpinupForcing(spin_sequences::Tuple{SpinupStepWithAggregator{F},Vararg}, forcing, tem_helpers, spinup_forcing) where {F}
    if F ∉ keys(spinup_forcing)
        seq_forc = getSpinupForcing(forcing, first(spin_sequences), tem_helpers.vals.forcing_types)
        spinup_forcing = set_namedtuple_field(spinup_forcing, (F, seq_forc))
    end
    return accumulateSpinupForcing(Base.tail(spin_sequences), forcing, tem_helpers, spinup_forcing)
end

function accumulateSpinupForcing(spin_sequences::Tuple{SpinupStepWithAggregator{F},Vararg}, forcing, tem_helpers, spinup_forcing) where {F}
    if F ∉ keys(spinup_forcing)
        seq_forc = getSpinupForcing(forcing, first(spin_sequences), tem_helpers.vals.forcing_types)
        spinup_forcing = set_namedtuple_field(spinup_forcing, (F, seq_forc))
    end
    return accumulateSpinupForcing(Base.tail(spin_sequences), forcing, tem_helpers, spinup_forcing)
end

coreTEM

Sindbad.Simulation.coreTEM Function

julia

coreTEM(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_land, tem_info, spinup_mode)

Runs the SINDBAD Terrestrial Ecosystem Model (TEM) for a single location, with or without spinup, based on the specified spinup_mode.

Arguments:

selected_models: A tuple of all models selected in the given model structure.
loc_forcing: A forcing NamedTuple containing the time series of environmental drivers for a single location.
loc_spinup: A NamedTuple with the spinup sequence of the location and the forcing derived from it, used to initialize the model to a steady state (only used if spinup is enabled).
loc_forcing_t: A forcing NamedTuple for a single location and a single time step.
loc_land: Initial SINDBAD land NamedTuple with all fields and subfields.
tem_info: A helper NamedTuple containing necessary objects for model execution and type consistencies.
spinup_mode: A type that determines whether spinup is included or excluded

Returns:

land_time_series: A vector of SINDBAD land states for each time step after the model simulation.

source

Code

julia

function coreTEM end

function coreTEM(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_land, tem_info, spinup_mode)

    land_prec = precomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers)

    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, spinup_mode)

    land_time_series = timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_spin, tem_info, tem_info.run.debug_model)
    return land_time_series
end

function coreTEM(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_land, tem_info, spinup_mode)

    land_prec = precomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers)

    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, spinup_mode)

    land_time_series = timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_spin, tem_info, tem_info.run.debug_model)
    return land_time_series
end

function coreTEM(selected_models, loc_forcing, loc_spinup, loc_forcing_t, land_time_series, loc_land, tem_info, spinup_mode)
    land_prec = precomputeTEM(selected_models, loc_forcing_t, loc_land, tem_info.model_helpers)

    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, spinup_mode)

    timeLoopTEM(selected_models, loc_forcing, loc_forcing_t, land_time_series, land_spin, tem_info, tem_info.run.debug_model)
    return nothing
end

coreTEM!

Sindbad.Simulation.coreTEM! Function

julia

coreTEM!(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_output, loc_land, tem_info)
coreTEM!(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_output, loc_restart, loc_land, tem_info)

Executes the core SINDBAD Terrestrial Ecosystem Model (TEM) for a single location, including precomputations, spinup, and the main time loop.

Arguments:

selected_models: A tuple of all models selected in the given model structure, or a TimeVaryingModels when parameters vary in time. Precompute and spinup then use the parameters of the first time step, and timeLoopTEM! updates them.
loc_forcing: A forcing NamedTuple containing the time series of environmental drivers for a single location.
loc_spinup: A NamedTuple with the spinup sequence of the location and the forcing derived from it, used to initialize the model to a steady state (only used if spinup is enabled).
loc_forcing_t: A forcing NamedTuple for a single location and a single time step.
loc_output: An output array or view for storing the model outputs for a single location.
loc_restart: The restart arrays of the location, where the pools after the spinup step with save_restart are stored. With nothing, as in the 7-argument form used by optimization and hybrid runs, nothing is stored.
loc_land: Initial SINDBAD land NamedTuple with all fields and subfields.
tem_info: A helper NamedTuple containing necessary objects for model execution and type consistencies.

Details

Executes the main TEM simulation logic with the provided parameters and data. Handles both regular simulation and spinup modes based on the spinup_mode flag.

Extended help

Precomputations:
- The function runs precomputeTEM to prepare the land state for the simulation.

source

Code

julia

function coreTEM!(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_output, loc_land, tem_info)
    return coreTEM!(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_output, nothing, loc_land, tem_info)
end

function coreTEM!(selected_models, loc_forcing, loc_spinup, loc_forcing_t, loc_output, loc_restart, loc_land, tem_info)
    # update the loc_forcing with the actual location
    loc_forcing_t = getForcingForTimeStep(loc_forcing, loc_forcing_t, 1, tem_info.vals.forcing_types)
    # precompute and spinup use the parameters of the first time step when they vary
    # in time, and the time loop updates them
    loc_models = initialModels(selected_models)
    # run precompute
    land_prec = precomputeTEM(loc_models, loc_forcing_t, loc_land, tem_info.model_helpers)
    # run spinup, which stores the pools for the restart file after the spinup step with
    # `save_restart`
    land_spin = spinupTEM(loc_models, loc_spinup, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM, loc_restart)

    timeLoopTEM!(selected_models, loc_forcing, loc_forcing_t, loc_output, land_spin, tem_info.vals.forcing_types, tem_info.model_helpers, tem_info.vals.output_vars, tem_info.n_timesteps, tem_info.run.debug_model)
    return nothing
end

dropVariablesFromStore

Sindbad.Simulation.dropVariablesFromStore Function

julia

dropVariablesFromStore(store_path, names)

Deletes the variables names from the zarr store at store_path.

Notes:

A zarr store keeps the values of each dimension, e.g., lat or d_cEco, as an array of its own. A dimension that no variable uses after the deletion is deleted too, so the store has no leftover dimensions.
A lazy run uses it to delete the restart variables from all_variables.zarr once they are saved to the restart file.

source

Code

julia

function dropVariablesFromStore(store_path, names)
    foreach(name -> rm(joinpath(store_path, string(name)); recursive=true, force=true), names)
    array_dims = Dict{String,Vector{String}}()
    for name ∈ readdir(store_path)
        zattrs_file = joinpath(store_path, name, ".zattrs")
        isfile(joinpath(store_path, name, ".zarray")) && isfile(zattrs_file) || continue
        array_dims[name] = String.(parsefile(zattrs_file; dicttype=Dict{String,Any})["_ARRAY_DIMENSIONS"])
    end
    is_coordinate(name) = array_dims[name] == [name]
    used = union(String[], (dims for (name, dims) ∈ array_dims if !is_coordinate(name))...)
    for name ∈ keys(array_dims)
        is_coordinate(name) && name ∉ used && rm(joinpath(store_path, name); recursive=true)
    end
    return nothing
end

getAllSpinupForcing

Sindbad.Simulation.getAllSpinupForcing Function

julia

getAllSpinupForcing(forcing, spin_seq, tem_helpers)

prepare the spinup forcing all forcing setups in the spinup sequence

Arguments:

forcing: a forcing NT that contains the forcing time series set for ALL locations
spin_seq: a sequence of information to carry out spinup at different steps with information on models to use, forcing, stopping critera, etc.
tem_helpers: helper NT with necessary objects for model run and type consistencies

Returns:

spinup_forcing::NamedTuple: A NamedTuple containing aggregated forcing data for each spinup sequence

Examples

julia

julia> using Sindbad

julia> # Prepare spinup forcing for all sequences
julia> # spinup_forcing = getAllSpinupForcing(forcing, spin_sequences, tem_helpers)

source

Code

julia

function getAllSpinupForcing(forcing, spin_sequences::Tuple{Vararg{SpinupStepWithAggregator}}, tem_helpers)
    return accumulateSpinupForcing(spin_sequences, forcing, tem_helpers, (;))
end

getForcingForTimeStep

Sindbad.Simulation.getForcingForTimeStep Function

julia

getForcingForTimeStep(forcing, loc_forcing_t, ts, Val{forc_with_type})

Get forcing values for a specific time step based on the forcing type.

Arguments:

forcing: a forcing NT that contains the forcing time series set for ALL locations
loc_forcing_t: a forcing NT for a single timestep to be reused in every time step
ts: time step to get the forcing for
forc_with_type: Value type parameter specifying the forcing type

source

getLocData

Sindbad.Simulation.getLocData Function

julia

getLocData(forcing, output_array, loc_ind)

Arguments:

forcing: a forcing NT that contains the forcing time series set for ALL locations
output_array: an output array/view for ALL locations
loc_ind: a tuple with the spatial indices of the data for a given location

source

julia

getLocData(forcing, output_array, loc_ind)

Arguments:

output_array: an output array/view for ALL locations
loc_ind: a tuple with the spatial indices of the data for a given location

source

julia

getLocData(forcing, output_array, loc_ind)

Arguments:

forcing: a forcing NT that contains the forcing time series set for ALL locations
loc_ind: a tuple with the spatial indices of the data for a given location

source

Code

julia

function getLocData(forcing::NamedTuple, output_array::AbstractArray, loc_ind)
    loc_forcing = getLocData(forcing, loc_ind)
    loc_output = getLocData(output_array, loc_ind)
    return loc_forcing, loc_output
end

function getLocData(output_array::AbstractArray, loc_ind)
    loc_output = map(output_array) do a
        view_at_trailing_indices(a, loc_ind)
    end
    return loc_output
end

function getLocData(forcing::NamedTuple, loc_ind)
    loc_forcing = map(forcing) do a
        # a variable with fewer dimensions than the domain's spatial index tuple has none of
        # the domain's spatial dimensions at all (e.g. a globally-uniform forcing such as
        # atmospheric CO2 with no lat/lon axis) -- it is constant across the domain, so it is
        # passed through unindexed rather than attempting a spatial view it cannot support
        if ndims(a) < length(loc_ind)
            Array(a)
        else
            view_at_trailing_indices(a, loc_ind) |> Array
        end
    end
    return loc_forcing
end

getLoopRanges

Sindbad.Simulation.getLoopRanges Function

julia

getLoopRanges(info, xmap_cubes; grids_per_thread=8)

Returns the cache, the block of grids and the run filter of a lazy run, as max_cache, loop_ranges and run_filter of runTEMLazy.

Arguments:

info: a SINDBAD NT with all information needed for setup and execution of an experiment
xmap_cubes: the xmap field returned by prepTEMLazy
grids_per_thread: the number of valid grids per thread that an "auto" block aims for

Returns:

the cache in bytes
the number of steps of each loop dimension of the xmap in a block, or nothing to let DiskArrayEngine choose the block
a ValidGridsFilter that skips the blocks without a valid grid, or nothing when no forcing without time marks the valid grids

Notes:

exe_rules.yax_block_grids sets the block, with the sizes in the order of data_dimension.space:
- "auto" or null: whole chunks of the input cubes with time. The step in each spatial dimension is the least common multiple of their chunk sizes, so no such chunk is read for more than one block. The block grows by whole steps, keeping its sides near equal, while its fullest block has at most grids_per_thread valid grids per thread and it fits the memory, see getMaxGridsPerBlock. With PmapParallelization, these are the threads of one worker. The valid grids are counted per block, so a block over land gets the intended work and one over ocean has less, instead of every block getting the average.
- [n_lat, n_lon]: a fixed block of grids, e.g. for benchmarks.
- "optimizer": the block that DiskArrayEngine optimizes from its model of read time. That model does not count the invalid grids, nor that a chunk of a forcing with all time steps is decompressed whole for every block that reads it.
With every setting, a block without a valid grid is skipped and not read.
A block is cut at the end of each dimension. A spatial subset that does not start at a chunk boundary shifts the chunks against the blocks, so some chunks are then read by two blocks.
The output chunks follow the block in space.

source

Code

julia

function getLoopRanges(info, xmap_cubes; grids_per_thread=8)
    setting = something(get(info.experiment.exe_rules, :yax_block_grids, "auto"), "auto")
    data_dimension = info.experiment.data_settings.forcing.data_dimension
    space_dims = Symbol.(data_dimension.space)
    time_name = Symbol(data_dimension.time)
    in_cubes = xmap_cubes.in_cubes
    op = first(values(xmap_cubes.dataset.cubes)).data.op
    # the loop dimensions of the xmap are those of all its inputs, in this order,
    # followed by those that only the outputs have
    input_names = DD.name.(DD.combinedims(in_cubes...; val=true, type=false, msg=nothing))
    length(input_names) <= length(op.windowsize) || error("the $(length(input_names)) dimensions of the input cubes are more than the $(length(op.windowsize)) loop dimensions of the lazy run")
    i_space = map(d -> findfirst(==(d), input_names), Tuple(space_dims))
    any(isnothing, i_space) && error("the spatial dimensions $(space_dims) are not all among the dimensions $(input_names) of the input cubes")
    space_sizes = map(i -> op.windowsize[i], i_space)
    expected_sizes = map(d -> length(DD.dims(first(c for c ∈ in_cubes if DD.hasdim(c, d)), d)), Tuple(space_dims))
    space_sizes == expected_sizes || error("the loop sizes $(space_sizes) of the spatial dimensions $(space_dims) do not match their lengths $(expected_sizes)")
    valid_grids = getValidGrids(in_cubes[1:xmap_cubes.n_forcing], space_dims, time_name)
    valid_sums = getValidGridSums(valid_grids)
    run_filter = isnothing(valid_sums) ? nothing : ValidGridsFilter(valid_sums, i_space)
    setting == "optimizer" && return getYaxMaxCache(info, xmap_cubes), nothing, run_filter
    bytes_per_grid = sum(cube -> getBytesPerGrid(cube, space_dims), (in_cubes..., values(xmap_cubes.dataset.cubes)...))
    block = if setting == "auto"
        steps = map(d -> getChunkStep(in_cubes, d, time_name), Tuple(space_dims))
        max_valid = grids_per_thread * getThreadsPerBlock(info.helpers.run.parallelization)
        max_grids = getMaxGridsPerBlock(info, bytes_per_grid)
        fits = b -> prod(b) <= max_grids && maximum(getValidCountsPerBlock(valid_sums, b, space_sizes)) <= max_valid
        growBlock(steps, space_sizes, fits)
    elseif setting isa AbstractVector && length(setting) == length(space_dims)
        map((n, s) -> clamp(Int(n), 1, s), Tuple(setting), space_sizes)
    else
        error("experiment.exe_rules.yax_block_grids must be \"auto\", \"optimizer\", or the grids of a block in each of $(space_dims), got `$(setting)`")
    end
    loop_steps = ntuple(length(op.windowsize)) do i
        j = findfirst(==(i), i_space)
        isnothing(j) ? op.windowsize[i] : block[j]
    end
    n_grids = prod(block)
    n_blocks = prod(cld.(space_sizes, block))
    valid_counts = getValidCountsPerBlock(valid_sums, block, space_sizes)
    n_skipped = count(iszero, valid_counts)
    parallelization = info.helpers.run.parallelization
    run_threads = "$(getProcsPerBlock(parallelization)) x $(getThreadsPerBlock(parallelization))"
    print_info(getLoopRanges, @__FILE__, @__LINE__, "$(n_blocks) blocks of $(join(block, "x")) grids in $(join(space_dims, "x")), up to $(maximum(valid_counts)) valid grids per block for $(run_threads) threads, $(round(n_grids * bytes_per_grid / 1e9; digits=2)) GB per block. $(n_skipped) blocks without a valid grid are skipped.", n_m=4)
    return n_grids * bytes_per_grid, loop_steps, run_filter
end

getOutDims

Sindbad.Simulation.getOutDims Function

julia

getOutDims(info, forcing_helpers[, ::OutputArray | ::OutputMArray | ::OutputSizedArray | ::OutputYAXArray])

Retrieves the dimensions for SINDBAD output based on the specified array backend.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
forcing_helpers: A NamedTuple with information on forcing sizes and dimensions.
::OutputArray: (Optional) A type dispatch for using a base Array as the array backend.
::OutputMArray: (Optional) A type dispatch for using MArray as the array backend.
::OutputSizedArray: (Optional) A type dispatch for using SizedArray as the array backend.
::OutputYAXArray: (Optional) A type dispatch for using YAXArray as the array backend.

Returns:

A vector of output dimensions, where each dimension is represented as a tuple of Dim objects.

Notes:

For OutputArray, OutputMArray, and OutputSizedArray, all dimensions are included.
For OutputYAXArray, spatial dimensions are excluded, and metadata is added for each variable.

Examples

julia

julia> using Sindbad

julia> # Get output dimensions with default array type
julia> # outdims = getOutDims(info, forcing_helpers)

julia> # Get output dimensions with YAXArray
julia> # outdims = getOutDims(info, forcing_helpers, OutputYAXArray())

source

Code

julia

function getOutDims end

function getOutDims(info, forcing_helpers)
    outdims = getOutDims(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims
end

function getOutDims(info, forcing_helpers)
    outdims = getOutDims(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims
end

function getOutDims(info, forcing_helpers, ::Union{OutputArray, OutputMArray, OutputSizedArray})
    outdims_pairs = getOutDimsPairs(info.output, forcing_helpers)
    outdims = map(outdims_pairs) do dim_pairs
        od = []
        for _dim in dim_pairs
            push!(od, YAXArrays.Dim{first(_dim)}(last(_dim)))
        end
        Tuple(od)
    end
    return outdims
end

function getOutDims(info, forcing_helpers, ::OutputYAXArray)
    outdims_pairs = getOutDimsPairs(info.output, forcing_helpers);
    space_dims = Symbol.(info.experiment.data_settings.forcing.data_dimension.space)
    var_dims = map(outdims_pairs) do dim_pairs
        od = []
        for _dim in dim_pairs
            if first(_dim) ∉ space_dims
                push!(od, DD.rebuild(DD.Dimensions.name2dim(first(_dim)),(last(_dim))))
            end
        end
        Tuple(od)
    end
    v_index = 1
    outdims = map(info.output.variables) do vname_full
        vname = string(last(vname_full))
        _properties = collectMetadata(info, vname_full)
        vdims = var_dims[v_index]
        # outformat = info.settings.experiment.model_output.format
        outformat = info.output.format # ? test and see if this works, which it should!
        backend = outformat == "nc" ? :netcdf : :zarr
        out_dim = YAXArrays.XOutput(vdims...)#=;
        properties = _properties,
        path=joinpath(info.output.dirs.data, "$(vname).$(outformat)"),
        backend=backend,
        overwrite=true)=#
        v_index += 1
        out_dim
    end
    return outdims
end

function getOutDimsArrays end

function getOutDimsArrays(info, forcing_helpers)
    outdims, outarray = getOutDimsArrays(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers)
    outdims, outarray = getOutDimsArrays(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers, oarr::OutputArray)
    outdims = getOutDims(info, forcing_helpers, oarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers, omarr::OutputMArray)
    outdims = getOutDims(info, forcing_helpers, omarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    marray = MArray{Tuple{size(outarray)...},eltype(outarray)}(undef)
    return outdims, marray
end

function getOutDimsArrays(info, forcing_helpers, osarr::OutputSizedArray)
    outdims = getOutDims(info, forcing_helpers, osarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    sized_array = SizedArray{Tuple{size(outarray)...},eltype(outarray)}(undef)
    return outdims, sized_array
end

function getOutDimsArrays(info, forcing_helpers, oayax::OutputYAXArray)
    outdims = getOutDims(info, forcing_helpers, oayax)
    outarray = nothing
    return outdims, outarray
end

function getOutDimsPairs(tem_output, forcing_helpers; dthres=1, with_time=true)
    forcing_axes = forcing_helpers.axes
    dim_loops = first.(forcing_axes)
    axes_dims_pairs = []
    if !isnothing(forcing_helpers.dimensions.permute)
        dim_perms = Symbol.(forcing_helpers.dimensions.permute)
        if dim_loops !== dim_perms
            for ix in eachindex(dim_perms)
                dp_i = dim_perms[ix]
                dl_ind = findall(x -> x == dp_i, dim_loops)[1]
                f_a = forcing_axes[dl_ind]
                ax_dim = Pair(first(f_a), last(f_a))
                push!(axes_dims_pairs, ax_dim)
            end
        end
    else
        axes_dims_pairs = map(x -> Pair(first(x), last(x)), forcing_axes)
    end
    opInd = 1
    outdims_pairs = map(tem_output.variables) do vname_full
        depth_info = tem_output.depth_info[opInd]
        depth_size = first(depth_info)
        depth_name = last(depth_info)
        od = []
        with_time && push!(od, axes_dims_pairs[1])
        if depth_size > dthres
            if depth_size == 1
                depth_name = "idx"
            end
            push!(od, Pair(Symbol(depth_name), (1:depth_size)))
        end
        foreach(axes_dims_pairs[2:end]) do f_d
            push!(od, f_d)
        end
        opInd += 1
        Tuple(od)
    end
    return outdims_pairs
end

getOutDimsArrays

Sindbad.Simulation.getOutDimsArrays Function

julia

getOutDimsArrays(info, forcing_helpers[, ::OutputArray | ::OutputMArray | ::OutputSizedArray | ::OutputYAXArray])

Retrieves the dimensions and corresponding data for SINDBAD output based on the specified array backend.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
forcing_helpers: A NamedTuple with information on forcing sizes and dimensions.
::OutputArray: (Optional) A type dispatch for using a base Array as the array backend.
::OutputMArray: (Optional) A type dispatch for using MArray as the array backend.
::OutputSizedArray: (Optional) A type dispatch for using SizedArray as the array backend.
::OutputYAXArray: (Optional) A type dispatch for using YAXArray as the array backend.

Returns:

A tuple (outdims, outarray):
- outdims: A vector of output dimensions, where each dimension is represented as a tuple of Dim objects.
- outarray: The corresponding data array, initialized based on the specified array backend.

Notes:

For OutputArray, OutputMArray, and OutputSizedArray, the data array is initialized with the appropriate backend type.
For OutputYAXArray, the data array is set to nothing, as YAXArray handles data differently.

Examples:

Using a base Array:

julia

outdims, outarray = getOutDimsArrays(info, forcing_helpers, OutputArray())

Using MArray:

julia

outdims, outarray = getOutDimsArrays(info, forcing_helpers, OutputMArray())

Using SizedArray:

julia

outdims, outarray = getOutDimsArrays(info, forcing_helpers, OutputSizedArray())

Using YAXArray:

julia

outdims, outarray = getOutDimsArrays(info, forcing_helpers, OutputYAXArray())

Default behavior:

julia

outdims, outarray = getOutDimsArrays(info, forcing_helpers)

source

Code

julia

function getOutDimsArrays end

function getOutDimsArrays(info, forcing_helpers)
    outdims, outarray = getOutDimsArrays(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers)
    outdims, outarray = getOutDimsArrays(info, forcing_helpers, info.helpers.run.output_array_type)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers, oarr::OutputArray)
    outdims = getOutDims(info, forcing_helpers, oarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    return outdims, outarray
end

function getOutDimsArrays(info, forcing_helpers, omarr::OutputMArray)
    outdims = getOutDims(info, forcing_helpers, omarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    marray = MArray{Tuple{size(outarray)...},eltype(outarray)}(undef)
    return outdims, marray
end

function getOutDimsArrays(info, forcing_helpers, osarr::OutputSizedArray)
    outdims = getOutDims(info, forcing_helpers, osarr)
    outarray = getNumericArrays(info, forcing_helpers.sizes)
    sized_array = SizedArray{Tuple{size(outarray)...},eltype(outarray)}(undef)
    return outdims, sized_array
end

function getOutDimsArrays(info, forcing_helpers, oayax::OutputYAXArray)
    outdims = getOutDims(info, forcing_helpers, oayax)
    outarray = nothing
    return outdims, outarray
end

getRestartDimsArrays

Sindbad.Simulation.getRestartDimsArrays Function

julia

getRestartDimsArrays(info, forcing_helpers, loc_land)

Prepares the variables, dimensions and arrays of the restart file, which holds the main pools after the spinup step with save_restart.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
forcing_helpers: A NamedTuple with information on forcing sizes and dimensions.
loc_land: The land of one location after a model run of one time step.

Returns:

A NamedTuple with:
- variables: the field and subfield pairs from getRestartVars.
- dims_pairs: the dimension pairs of each variable without the time dimension. A layer dimension is only included for variables with more than one layer.
- dims: the Dim objects built from dims_pairs, used for saving.
- data: one array per variable with the layout (1, depth, space...).

Notes:

The arrays keep a leading time dimension of size 1 so that getLocData and setOutputForTimeStep! work on them the same way as on the output arrays.

source

Code

julia

function getRestartDimsArrays(info, forcing_helpers, loc_land)
    restart_vars = getRestartVars(info, loc_land)
    # the dimension names follow the output settings, but the sizes always come from
    # land so that every layer of a variable is stored
    depth_info = map(restart_vars, getLandVarDepthInfo(info, loc_land, restart_vars)) do var_pair, d_info
        length(getproperty(getproperty(loc_land, first(var_pair)), last(var_pair))), last(d_info)
    end
    restart_info = (; variables=restart_vars, depth_info=depth_info)
    dims_pairs = getOutDimsPairs(restart_info, forcing_helpers; with_time=false)
    restart_dims = map(dims_pairs) do dim_pairs
        Tuple(YAXArrays.Dim{first(_dim)}(last(_dim)) for _dim in dim_pairs)
    end
    num_type = info.helpers.numbers.num_type
    space_sizes = Base.tail(values(forcing_helpers.sizes))
    restart_array = map(depth_info) do d_info
        fill(num_type(NaN), 1, first(d_info), space_sizes...)
    end
    return (; variables=restart_vars, dims_pairs=dims_pairs, dims=restart_dims, data=[restart_array...])
end

getSpatialInfo

Sindbad.Simulation.getSpatialInfo Function

julia

getSpatialInfo(forcing_helpers)
getSpatialInfo(forcing, filter_nan_pixels)

get the information of the indices of the data to run the model for. The second variant additionally filters out the grids with invalid forcing, see filterNanGrids.

Arguments:

forcing: a forcing NT that contains the forcing time series set for ALL locations

source

Code

julia

function getSpatialInfo end

function getSpatialInfo(forcing_helpers)
    @debug "     getting the space locations to run the model loop"
    forcing_sizes = forcing_helpers.sizes
    loopvars = collect(keys(forcing_sizes))
    additionaldims = setdiff(loopvars, [Symbol(forcing_helpers.dimensions.time)])
    spacesize = values(forcing_sizes[additionaldims])
    loc_space_maps = vec(collect(Iterators.product(Base.OneTo.(spacesize)...)))
    loc_space_maps = map(loc_space_maps) do loc_names
        map(zip(loc_names, additionaldims)) do (loc_index, lv)
            lv => loc_index
        end
    end
    loc_space_maps = Tuple(loc_space_maps)
    space_ind = Tuple([Tuple(last.(loc_space_map)) for loc_space_map ∈ loc_space_maps])
    return space_ind, loc_space_maps
end

function getSpatialInfo(forcing_helpers)
    @debug "     getting the space locations to run the model loop"
    forcing_sizes = forcing_helpers.sizes
    loopvars = collect(keys(forcing_sizes))
    additionaldims = setdiff(loopvars, [Symbol(forcing_helpers.dimensions.time)])
    spacesize = values(forcing_sizes[additionaldims])
    loc_space_maps = vec(collect(Iterators.product(Base.OneTo.(spacesize)...)))
    loc_space_maps = map(loc_space_maps) do loc_names
        map(zip(loc_names, additionaldims)) do (loc_index, lv)
            lv => loc_index
        end
    end
    loc_space_maps = Tuple(loc_space_maps)
    space_ind = Tuple([Tuple(last.(loc_space_map)) for loc_space_map ∈ loc_space_maps])
    return space_ind, loc_space_maps
end

function getSpatialInfo(forcing, filter_nan_pixels)
    space_ind, loc_space_maps = getSpatialInfo(forcing.helpers)
    loc_space_maps = filterNanGrids(forcing, loc_space_maps, filter_nan_pixels)
    space_ind = Tuple([Tuple(last.(loc_space_map)) for loc_space_map ∈ loc_space_maps])
    return space_ind
end

getSpinupSequence

Sindbad.Simulation.getSpinupSequence Function

julia

getSpinupSequence(method, spinup_config, loc_forcing, helpers_dates, aggregator_cache)

Build the spinup sequence of one location for the given sequence method.

Arguments:

method: a SINDBAD SpinupSequence that determines how the steps are chosen
spinup_config: the spinup settings, with the sequence method and its options
loc_forcing: a forcing NT for a single location, from which a method may read variables
helpers_dates: a NT with date-related helpers of the experiment
aggregator_cache: a Dict of already built temporal aggregators, shared across locations

Returns:

a tuple of SpinupStepWithAggregator, one per step of the sequence

spinup sequence method:

SpinupSequence

Abstract type for methods that build the spinup sequence of a location in SINDBAD

Available methods/subtypes:

SequenceDefault: spin the selected spinup models up on the mean seasonal cycle, bracketed by a run over all years, and use this when the experiment settings give no sequence
SequenceList: use the explicit list of spinup sequences given in the experiment settings
SequenceWithAge: build the spinup sequence from the disturbance year of the location so that the years since disturbance are replayed
SequenceWithSpinupInput: build the spinup sequence with spinup input from previously saved cEco and TWS pools through forcing json

Extended help

Notes:

The sequence is built here rather than while reading the settings because a method may need the forcing of a location, which only exists once the data is loaded.
To add a method, define a type that subtypes SpinupSequence in SimulationTypes.jl and a getSpinupSequence method for it here. The method name in the settings is the type name in snake case, so SequenceWithAge is written as sequence_with_age.

source

Code

julia

function getSpinupSequenceSteps(sequence)
    steps = [0]
    for (seq_index, spin_seq) ∈ enumerate(sequence)
        append!(steps, fill(seq_index, spin_seq.n_repeat * spinupStoreSteps(spin_seq)))
    end
    return steps
end

getYaxMaxCache

Sindbad.Simulation.getYaxMaxCache Function

julia

getYaxMaxCache(info, in_cubes, out_cubes=(); n_forcing=length(in_cubes),
    grids_per_thread=4, memory_fraction=0.66)
getYaxMaxCache(info, xmap_cubes)

Returns the cache size of a cube run, compute_to_zarr or mapCube, from the yax_max_cache setting. A number in the setting is used as it is, and with null the cache is computed.

Arguments:

info: a SINDBAD NT with all information needed for setup and execution of an experiment
in_cubes: the input cubes of the run, starting with the forcing
out_cubes: the output cubes of the run, empty when they are not built yet
xmap_cubes: the xmap field returned by prepTEMLazy, which holds both
n_forcing: the number of forcing cubes at the start of in_cubes
grids_per_thread: the number of valid grids per thread in a block
memory_fraction: the largest share of the available memory for the run

Notes:

A cube run reads, runs and writes one block of grids at a time, and the threads only share the grids of that block. A grid with invalid forcing returns at once, so a block with fewer valid grids than threads leaves the other threads idle. The computed cache fits grids_per_thread valid grids per thread.
The bytes of a grid are summed over the input and output cubes. A cube without a spatial dimension is read once for all grids of a block and is not counted.
The share of valid grids is counted from the forcing without time, since reading the forcing with time costs as much as the run.
The cache only holds the input and output buffers of a block. Reading or writing a block can copy each buffer once more, and every thread holds the data of its own grid, so the run is taken to need twice the cache and one grid per thread. This is kept within memory_fraction of the memory that is available when the cache is computed, which excludes what the session already holds.
With PmapParallelization, the threads are those of one worker, and every worker holds its own block, see getCacheMemoryLimit.

source

Code

julia

function getYaxMaxCache(info, in_cubes, out_cubes=(); n_forcing=length(in_cubes), grids_per_thread=4, memory_fraction=0.66)
    max_cache = get(info.experiment.exe_rules, :yax_max_cache, nothing)
    max_cache isa Real && return max_cache
    isnothing(max_cache) || error("experiment.exe_rules.yax_max_cache must be a number of bytes, or null to compute it, got `$(max_cache)`")
    data_dimension = info.experiment.data_settings.forcing.data_dimension
    space_dims = Symbol.(data_dimension.space)
    time_name = Symbol(data_dimension.time)
    bytes_per_grid = sum(cube -> getBytesPerGrid(cube, space_dims), (in_cubes..., out_cubes...))
    valid_fraction = getValidGridFraction(in_cubes[1:n_forcing], space_dims, time_name)
    parallelization = info.helpers.run.parallelization
    n_threads = getThreadsPerBlock(parallelization)
    n_procs = getProcsPerBlock(parallelization)
    n_grids = ceil(Int, grids_per_thread * n_threads / valid_fraction)
    cache = n_grids * bytes_per_grid
    print_info(getYaxMaxCache, @__FILE__, @__LINE__, "computed cache of $(round(cache / 1e9; digits=2)) GB: $(n_grids) grids per block at $(round(bytes_per_grid / 1e6; digits=2)) MB each, for $(grids_per_thread) valid grids on each of $(n_threads) threads with $(round(100 * valid_fraction; digits=1))% of grids valid", n_m=4)
    memory_limit, available_memory = getCacheMemoryLimit(bytes_per_grid, parallelization; memory_fraction=memory_fraction)
    if cache > memory_limit
        @warn "The computed cache of $(round(cache / 1e9; digits=2)) GB needs up to $(round(((n_procs == 1 ? 2 : 2 * n_procs + 1) * cache + n_procs * n_threads * bytes_per_grid) / 1e9; digits=2)) GB for the run, more than $(round(Int, 100 * memory_fraction))% of the $(round(available_memory / 1e9; digits=2)) GB of available memory. It is set to $(round(memory_limit / 1e9; digits=2)) GB, so some threads stay idle."
        cache = memory_limit
    end
    return cache
end

prepSequenceDefault

Sindbad.Simulation.prepSequenceDefault Function

julia

prepSequenceDefault(; nrepeat_base=200, forcing_msc="day_MSC")

Build the steps of the default spinup sequence: the selected spinup models are run over the mean seasonal cycle, bracketed on either side by a single run over all years.

Arguments:

nrepeat_base=200 [default]: the repeats of the mean seasonal cycle
forcing_msc: the name of the mean seasonal cycle spinup forcing set

Outputs

a vector of spinup sequence steps, to be given types by getSpinupSequenceWithTypes

source

Code

julia

function prepSequenceDefault(; nrepeat_base=200, forcing_msc="day_MSC")
    return [
        Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => "all_years", "n_repeat" => 1),
        Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => forcing_msc, "n_repeat" => nrepeat_base),
        Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => "all_years", "n_repeat" => 1),
    ]
end

prepSequenceWithAge

Sindbad.Simulation.prepSequenceWithAge Function

julia

prepSequenceWithAge(year_disturbance; nrepeat_base=200, year_start=1979, forcing_msc="day_MSC")

Build the steps of a spinup sequence that replays the age of a location, i.e. the years between its disturbance and the start of the experiment.

Arguments:

year_disturbance: the year in which the location was disturbed, as an integer
nrepeat_base=200 [default]: the repeats of the mean seasonal cycle before the disturbance
year_start: 1979 [default] start year of the experiment, as an integer
forcing_msc: the name of the mean seasonal cycle spinup forcing set

Outputs

a vector of spinup sequence steps, to be given types by getSpinupSequenceWithTypes

source

Code

julia

function prepSequenceWithAge(year_disturbance; nrepeat_base=200, year_start = 1979, forcing_msc = "day_MSC")
    nrepeat_age = nrepeatYearsAge(year_disturbance; year_start)
    sequence = [
        Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => "all_years", "n_repeat" => 1),
        Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => forcing_msc, "n_repeat" => nrepeat_base),
        Dict("spinup_mode" => "eta_scale_H", "forcing" => nothing, "models" => nothing, "n_repeat" => 1, "save_restart" => true),
    ]
    if nrepeat_age >= 0
        sequence = [
            Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => "all_years", "n_repeat" => 1),
            Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => forcing_msc, "n_repeat" => nrepeat_base),
            Dict("spinup_mode" => "clearance_A0", "forcing" => nothing, "models" => nothing, "n_repeat" => 1),
            Dict("spinup_mode" => "eta_scale_H", "forcing" => nothing, "models" => nothing, "n_repeat" => 1, "save_restart" => true),
            Dict("spinup_mode" => "loop_TEM", "models" => "spinup", "forcing" => forcing_msc, "n_repeat" => nrepeat_age),
        ]
    end
    return sequence
end

prepSpinupOut

Sindbad.Simulation.prepSpinupOut Function

julia

prepSpinupOut(info, run_helpers)

Collects the stored spinup of every location into arrays for the spinup file.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
run_helpers: the run helpers of an eager run, after runTEM!

Returns:

nothing when store_spinup is false. Otherwise a NamedTuple with:
- variables: the field and subfield pairs of the stored pools, e.g. (:pools, :cEco). The step index has nothing as its pair.
- names: the variable names in the file, the pool names and spinup_sequence_step
- dims: the Dim objects of each variable
- data: one array per variable
- sequences: the spinup sequence of every location
- space_ind: the spatial indices of every location, in the order of sequences

Notes:

Each pool has the dimensions (spinup_step, layer, space...). The layer and spatial dimensions are those of the same pool in the restart file, so the layer dimension is dropped for a pool with one layer.
The locations can have sequences of different lengths, e.g. with sequence_with_age. The spinup_step dimension has the length of the longest one, and the rows after the end of a shorter sequence are NaN.
spinup_sequence_step has the dimensions (spinup_step, space...). It is the index of the step of the spinup sequence that wrote each row, 0 for the state before the spinup, and -1 for the padding rows.

source

Code

julia

function prepSpinupOut end

prepSpinupOut(info, run_helpers) = prepSpinupOut(info, run_helpers, info.helpers.run.store_spinup)

prepSpinupOut(_, _, ::DoNotStoreSpinup) = nothing

function prepSpinupOut(info, run_helpers, ::DoStoreSpinup)
    restart = run_helpers.restart
    space_land = run_helpers.space_land
    space_ind = run_helpers.space_ind
    sequences = [loc_spinup.sequence for loc_spinup ∈ run_helpers.space_spinup]
    sequence_steps = map(getSpinupSequenceSteps, sequences)
    n_max = maximum(length, sequence_steps)
    step_dim = YAXArrays.Dim{:spinup_step}(1:n_max)
    num_type = info.helpers.numbers.num_type
    variables = []
    names = Symbol[]
    dims = []
    data = []
    space_dims = ()
    for pool_name ∈ (:cEco, :TWS)
        store_name = Symbol("spinup_", pool_name)
        loc_store = getproperty(first(space_land).states, store_name)
        isnothing(loc_store) && continue
        var_index = findfirst(==((:pools, pool_name)), restart.variables)
        isnothing(var_index) && continue
        n_layers = size(loc_store, 2)
        space_sizes = size(restart.data[var_index])[3:end]
        pool_data = fill(num_type(NaN), n_max, n_layers, space_sizes...)
        for (loc_land, loc_ind) ∈ zip(space_land, space_ind)
            loc_store = getproperty(loc_land.states, store_name)
            view_at_trailing_indices(pool_data, loc_ind)[1:size(loc_store, 1), :] .= loc_store
        end
        pool_dims = restart.dims[var_index]
        if n_layers == 1
            pool_data = dropdims(pool_data; dims=2)
        end
        space_dims = Tuple(pool_dims)[end-length(space_sizes)+1:end]
        push!(variables, (:pools, pool_name))
        push!(names, pool_name)
        push!(dims, (step_dim, pool_dims...))
        push!(data, pool_data)
    end
    isempty(data) && return nothing
    space_sizes = size(last(data))[end-length(space_dims)+1:end]
    step_data = fill(-1, n_max, space_sizes...)
    for (loc_steps, loc_ind) ∈ zip(sequence_steps, space_ind)
        view_at_trailing_indices(step_data, loc_ind)[1:length(loc_steps)] .= loc_steps
    end
    push!(variables, nothing)
    push!(names, :spinup_sequence_step)
    push!(dims, (step_dim, space_dims...))
    push!(data, step_data)
    return (; variables, names, dims, data, sequences, space_ind)
end

function prepSpinupOut(info, run_helpers, ::DoStoreSpinup)
    restart = run_helpers.restart
    space_land = run_helpers.space_land
    space_ind = run_helpers.space_ind
    sequences = [loc_spinup.sequence for loc_spinup ∈ run_helpers.space_spinup]
    sequence_steps = map(getSpinupSequenceSteps, sequences)
    n_max = maximum(length, sequence_steps)
    step_dim = YAXArrays.Dim{:spinup_step}(1:n_max)
    num_type = info.helpers.numbers.num_type
    variables = []
    names = Symbol[]
    dims = []
    data = []
    space_dims = ()
    for pool_name ∈ (:cEco, :TWS)
        store_name = Symbol("spinup_", pool_name)
        loc_store = getproperty(first(space_land).states, store_name)
        isnothing(loc_store) && continue
        var_index = findfirst(==((:pools, pool_name)), restart.variables)
        isnothing(var_index) && continue
        n_layers = size(loc_store, 2)
        space_sizes = size(restart.data[var_index])[3:end]
        pool_data = fill(num_type(NaN), n_max, n_layers, space_sizes...)
        for (loc_land, loc_ind) ∈ zip(space_land, space_ind)
            loc_store = getproperty(loc_land.states, store_name)
            view_at_trailing_indices(pool_data, loc_ind)[1:size(loc_store, 1), :] .= loc_store
        end
        pool_dims = restart.dims[var_index]
        if n_layers == 1
            pool_data = dropdims(pool_data; dims=2)
        end
        space_dims = Tuple(pool_dims)[end-length(space_sizes)+1:end]
        push!(variables, (:pools, pool_name))
        push!(names, pool_name)
        push!(dims, (step_dim, pool_dims...))
        push!(data, pool_data)
    end
    isempty(data) && return nothing
    space_sizes = size(last(data))[end-length(space_dims)+1:end]
    step_data = fill(-1, n_max, space_sizes...)
    for (loc_steps, loc_ind) ∈ zip(sequence_steps, space_ind)
        view_at_trailing_indices(step_data, loc_ind)[1:length(loc_steps)] .= loc_steps
    end
    push!(variables, nothing)
    push!(names, :spinup_sequence_step)
    push!(dims, (step_dim, space_dims...))
    push!(data, step_data)
    return (; variables, names, dims, data, sequences, space_ind)
end

prepTEM

Sindbad.Simulation.prepTEM Function

julia

prepTEM(forcing::NamedTuple, info::NamedTuple)
prepTEM(selected_models, forcing::NamedTuple, info::NamedTuple)
prepTEM(selected_models, forcing::NamedTuple, observations::NamedTuple, info::NamedTuple)

Prepares the SINDBAD Terrestrial Ecosystem Model (TEM) for execution by setting up the necessary inputs, outputs, and configurations with different variants for different experimental setups.

Arguments:

selected_models: A tuple of all models selected in the given model structure.
forcing::NamedTuple: A forcing NamedTuple containing the time series of environmental drivers for all locations.
observations::NamedTuple: A NamedTuple containing observational data for model validation.
info::NamedTuple: A nested NamedTuple containing necessary information, including:
- Helpers for running the model.
- Model configurations.
- Spinup settings.

Returns:

run_helpers: A NamedTuple containing preallocated data and configurations required to run the TEM, including:
- Spatial forcing data.
- Per-location spinup sequences and the forcing derived from them.
- Output arrays.
- Land variables.
- Temporal and spatial indices.
- Model and helper configurations.

Notes:

The function dynamically prepares the required data structures based on the specified PreAllocputType in info.helpers.run.land_output_type.
It handles spatial and temporal data preparation, including filtering NaN pixels, initializing land variables, and setting up forcing and output arrays.
This function is a key step in preparing the SINDBAD TEM for execution.

Examples:

Preparing TEM with observations:

julia

selected_models = (model1, model2)
forcing = (data = ..., variables = ...)
observations = (data = ..., variables = ...)
info = (helpers = ..., models = ..., spinup = ...)
run_helpers = prepTEM(selected_models, forcing, observations, info)

Preparing TEM without observations:

julia

selected_models = (model1, model2)
forcing = (data = ..., variables = ...)
info = (helpers = ..., models = ..., spinup = ...)
run_helpers = prepTEM(selected_models, forcing, info)

source

Code

julia

function prepTEM end

function prepTEM(forcing::NamedTuple, info::NamedTuple)
    selected_models = info.models.forward
    return prepTEM(selected_models, forcing, info)
end

function prepTEM(forcing::NamedTuple, info::NamedTuple)
    selected_models = info.models.forward
    return prepTEM(selected_models, forcing, info)
end

function prepTEM(selected_models, forcing::NamedTuple, info::NamedTuple)
    print_info(prepTEM, @__FILE__, @__LINE__, "preparing to run terrestrial ecosystem model (TEM)", n_f=1)
    output = prepTEMOut(info, forcing.helpers)
    print_info(prepTEM, @__FILE__, @__LINE__, "  preparing helpers for running model experiment", n_f=4)
    run_helpers = helpPrepTEM(selected_models, info, forcing, output, info.helpers.run.land_output_type)
    print_info_separator()

    return run_helpers
end

function prepTEM(selected_models, forcing::NamedTuple, observations::NamedTuple, info::NamedTuple)
    print_info(prepTEM, @__FILE__, @__LINE__, "preparing to run terrestrial ecosystem model (TEM)", n_f=1)
    output = prepTEMOut(info, forcing.helpers)
    run_helpers = helpPrepTEM(selected_models, info, forcing, observations, output, info.helpers.run.land_output_type)
    print_info_separator()

    return run_helpers
end

prepTEMLazy

Sindbad.Simulation.prepTEMLazy Function

julia

prepTEMLazy(selected_models, forcing::NamedTuple, info::NamedTuple)

Builds the lazy run of the TEM over the forcing cubes with xmap, without running it. runTEMLazy runs it.

Arguments:

selected_models: a tuple of all models selected in the given model structure
forcing: a forcing NT that contains the forcing time series set for ALL locations
info: a SINDBAD NT with all information needed for setup and execution of an experiment

Returns:

a NamedTuple with the output Dataset, the restart Dataset with the pools after spinup, the restart_vars written to it, and the spinup_mode of the step that saved the restart

source

Code

julia

function prepTEMLazy(selected_models, forcing::NamedTuple, info::NamedTuple)
    # forcing/input information, one cube per source so that no source is read twice
    forcing_cubes, forcing_dims, forcing_sources = shareForcingCubes(info, forcing)
    incubes = map(hideCubeType, forcing_cubes);
    indims = forcing_dims;
    n_forcing_cubes = length(incubes)
    # a spatial parameter has no core dimensions, so each location gets one value, and a
    # spatiotemporal one has time, so each location gets its values over time like a
    # forcing. The values set the parameters of the models of the location in
    # TEMLazy.
    parameter_cubes, parameter_dims, parameter_to_index, parameter_time = getParameterCubes(selected_models, forcing)
    incubes = (incubes..., parameter_cubes...)
    indims = (indims..., parameter_dims...)
    # information for running model
    run_helpers = prepTEM(forcing, info);
    _num_type = Val{info.helpers.numbers.num_type}()
    _forcing_data_info = getForcingDataInfo(info, forcing.variables)
    # every location starts from this land, so the land keeps one type in all models
    loc_land = getLandTemplate(selected_models, forcing, info, run_helpers.tem_info, run_helpers.loc_land, _forcing_data_info)
    checkOutputLayers(info.output.variables, info.output.depth_info, loc_land)
    #output = XOutput.(getproperty.(run_helpers.output_dims, :axisdesc))
    output = run_helpers.output_dims
    # the template land has every pool
    restart = getRestartDimsArrays(info, forcing.helpers, loc_land)
    restart_vars = restart.variables
    restart_output = getRestartOutDims(info, restart.dims_pairs)
    all_output = (output..., restart_output...)
    # DiskArrayEngine needs all output cubes to have the same dimensions, so every cube
    # is padded with singleton versions of the dimensions of the other cubes.
    alloutdims = union(getproperty.(all_output, :outaxes)...)
    new_output = map(all_output) do out
        singleton_dims = setdiff(alloutdims, out.outaxes)
        newdims = (out.outaxes..., DD.rebuild.(singleton_dims, (1:1,))...) 
        YAXArrays.XOutput(newdims, out.destroyaxes, out.outtype, out.properties)
    end
    tem_kwargs = (selected_models=selected_models,
        forcing_vars=forcing.variables,
        forcing_sources=forcing_sources,
        output_vars=run_helpers.output_vars,
        loc_land=loc_land,
        tem=run_helpers.tem_info,
        clean_data=(; _num_type, _forcing_data_info),
        restart_vars=restart_vars,
        parameter_to_index=parameter_to_index,
        parameter_time=parameter_time)
    outcubes = xmap(TEMLazyForGrid,
        (incubes .⊘ indims)...;
    # the arguments become part of the type of every function of the xmap run, so they
    # are passed behind a Ref{Any} to keep those types small and fast to compile
    function_kwargs = (tem_kwargs=Ref{Any}(tem_kwargs),),
        output = new_output,
        # xmap runs the pixels of a block on one thread unless this is set. The number of
        # threads of the Julia session decides how many are used.
        allow_threads = true,
        #outdims=run_helpers.output_dims,
        #max_cache=info.experiment.exe_rules.yax_max_cache,
        #ispar=false,
        )
    # the cubes are kept in a Vector{Any}, as code for the tuple of all their types is
    # slow to compile
    outcubes = outcubes isa YAXArray ? Any[outcubes] : collect(Any, outcubes)
    n_out = length(output)
    varnames = Tuple(getUniqueVarNames(info.output.variables))
    restart_names = Tuple(last.(restart_vars))
    # besides the padding, xmap adds a singleton for every core dimension of the forcing
    # that the output does not have, e.g. the soil layers of the texture forcing
    space_dims = Symbol.(info.experiment.data_settings.forcing.data_dimension.space)
    padded_names = Vector{Symbol}[setdiff(DD.name.(DD.dims(o_cube)), (DD.name.(out.outaxes)..., space_dims...)) for (out, o_cube) ∈ zip(collect(Any, all_output), outcubes)]
    dropped_cubes = Any[dropPaddedDims(o_cube, padded) for (o_cube, padded) ∈ zip(outcubes, padded_names)]
    output_ds = datasetFromCubes(varnames, dropped_cubes[1:n_out])
    restart_ds = datasetFromCubes(restart_names, dropped_cubes[(n_out+1):end])
    # a pool can also be an output variable, so the restart cubes get their own keys
    restart_keys = map(r -> Symbol("restart_", r), restart_names)
    xmap_keys = (varnames..., restart_keys...)
    xmap_cubes = (; dataset=datasetFromCubes(xmap_keys, outcubes),
        padded_names=Dict(zip(xmap_keys, padded_names)),
        output_names=varnames, restart_names=restart_names, restart_keys=restart_keys,
        in_cubes=incubes, n_forcing=n_forcing_cubes)
    return (; output=output_ds, restart=restart_ds, restart_vars=restart_vars, spinup_mode=getSaveRestartSpinupMode(info), xmap=xmap_cubes)
end

prepTEMOut

Sindbad.Simulation.prepTEMOut Function

julia

prepTEMOut(info::NamedTuple, forcing_helpers::NamedTuple)

Prepares the output NamedTuple required for running the Terrestrial Ecosystem Model (TEM) in SINDBAD.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
forcing_helpers: A NamedTuple with information on forcing sizes and dimensions.

Returns:

A NamedTuple output_tuple containing:

land_init: The initial land state from info.land_init.
dims: A vector of output dimensions, where each dimension is represented as a tuple of Dim objects.
data: A vector of numeric arrays initialized for output variables.
variables: A list of output variable names.
Additional fields for optimization output if optimization is enabled.

Notes:

The function initializes the output dimensions and data arrays based on the specified array backend (info.helpers.run.output_array_type).
If optimization is enabled (info.helpers.run.run_optimization), additional fields for optimized parameters are added to the output.
The function uses helper functions like getOutDimsArrays and setupOptiOutput to prepare the output.

Examples:

Basic usage:

julia

output_tuple = prepTEMOut(info, forcing_helpers)

Accessing output fields:

julia

dims = output_tuple.dims
data = output_tuple.data
variables = output_tuple.variables

source

Code

julia

function prepTEMOut(info::NamedTuple, forcing_helpers::NamedTuple)
    print_info(prepTEMOut, @__FILE__, @__LINE__, "preparing output helpers for Terrestrial Ecosystem Model (TEM)", n_f=4)
    land = info.helpers.land_init
    output_tuple = (;)
    output_tuple = set_namedtuple_field(output_tuple, (:land_init, land))
    @debug "     prepTEMOut: getting out variables, dimension and arrays..."
    outdims, outarray = getOutDimsArrays(info, forcing_helpers, info.helpers.run.output_array_type)
    output_tuple = set_namedtuple_field(output_tuple, (:dims, outdims))
    output_tuple = set_namedtuple_field(output_tuple, (:data, outarray))
    output_tuple = set_namedtuple_field(output_tuple, (:variables, info.output.variables))

    output_tuple = setupOptiOutput(info, output_tuple, info.helpers.run.run_optimization)
    @debug "\n----------------------------------------------\n"
    return output_tuple
end

runSpinupSequences

Sindbad.Simulation.runSpinupSequences Function

julia

runSpinupSequences(spin_seqs, selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart=nothing)

Runs every spinup sequence of spin_seqs in order, threading land through the sequences.

Arguments:

spin_seqs: a tuple of spinup sequences, each with the forcing name, number of timesteps, number of repeats, and the spinup mode
selected_models: a tuple of all models selected in the given model structure
spinup_forcings: a forcing NT with one entry per distinct spinup forcing name
loc_forcing_t: a forcing NT for a single location and a single time step
land: SINDBAD NT input to the spinup of TEM during which subfield(s) of pools are overwritten
tem_info: helper NT with necessary objects for model run and type consistencies
loc_restart: where the pools are stored for the restart file after the step whose save_restart is true, or nothing

Returns:

the updated land

Notes:

The sequences are consumed recursively with Base.tail rather than with a for loop. spin_seqs is a heterogeneous tuple, so a for loop indexes it with a runtime state, which puts the whole tuple and every extracted element on the heap once per iteration. Recursion specialises one method per sequence and removes those allocations.

source

Code

julia

function runSpinupSequences end

runSpinupSequences(::Tuple{}, selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart=nothing) = land

function runSpinupSequences(spin_seqs::Tuple, selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart=nothing)
    spin_seq = first(spin_seqs)
    n_repeat = spin_seq.n_repeat
    sel_forcing = sequenceForcing(spinup_forcings, spin_seq)
    spinup_models = getSpinupModels(selected_models, spin_seq.model_indices)
    land = spinupSequence(spinup_models, sel_forcing, loc_forcing_t, land, tem_info, spin_seq.n_timesteps, n_repeat, spin_seq.spinup_mode)
    spin_seq.save_restart && setRestart!(loc_restart, land, tem_info)
    return runSpinupSequences(Base.tail(spin_seqs), selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart)
end

function runSpinupSequences(spin_seqs::Tuple, selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart=nothing)
    spin_seq = first(spin_seqs)
    n_repeat = spin_seq.n_repeat
    sel_forcing = sequenceForcing(spinup_forcings, spin_seq)
    spinup_models = getSpinupModels(selected_models, spin_seq.model_indices)
    land = spinupSequence(spinup_models, sel_forcing, loc_forcing_t, land, tem_info, spin_seq.n_timesteps, n_repeat, spin_seq.spinup_mode)
    spin_seq.save_restart && setRestart!(loc_restart, land, tem_info)
    return runSpinupSequences(Base.tail(spin_seqs), selected_models, spinup_forcings, loc_forcing_t, land, tem_info, loc_restart)
end

runTEM

Sindbad.Simulation.runTEM Function

julia

runTEM(forcing::NamedTuple, info::NamedTuple)
runTEM(selected_models, forcing::NamedTuple, loc_spinup, loc_forcing_t, loc_land::NamedTuple, tem_info::NamedTuple)
runTEM(selected_models, loc_forcing::NamedTuple, loc_spinup, loc_forcing_t, land_time_series, loc_land::NamedTuple, tem_info::NamedTuple)

Runs the SINDBAD Terrestrial Ecosystem Model (TEM) for a single location, with or without spinup, based on the provided configurations. The two main variants are the ones with and without the preallocated land time series. The shorthand version with two input arguments calls the one without preallocated land time series.

Arguments:

selected_models: A tuple of all models selected in the given model structure.
forcing::NamedTuple: A forcing NamedTuple containing the time series of environmental drivers for all locations.
loc_spinup: A NamedTuple with the spinup sequence of the location and the forcing derived from it, used to initialize the model to a steady state.
loc_forcing_t: A forcing NamedTuple for a single location and a single time step.
loc_land::NamedTuple: Initial SINDBAD land NamedTuple with all fields and subfields.
tem_info::NamedTuple: A nested NamedTuple containing necessary information, including:
- Model helpers.
- Spinup configurations.
- Debugging options.
- Output configurations.

Returns:

LandWrapper: A wrapper containing the time series of SINDBAD land states for each time step after the model simulation.

Notes:

The function internally calls coreTEM to handle the main simulation logic.
If spinup is enabled (DoSpinupTEM), the function runs the spinup process before the main simulation.
If spinup is disabled (DoNotSpinupTEM), the function directly runs the main simulation.
The function prepares the necessary inputs and configurations using prepTEM before executing the simulation.

Examples

julia

julia> using Sindbad

julia> # Run TEM with spinup (shorthand version)
julia> # land_time_series = runTEM(forcing, info)

julia> # Run TEM with spinup (detailed version)
julia> # land_time_series = runTEM(selected_models, forcing, loc_spinup, loc_forcing_t, loc_land, tem_info)

julia> # Run TEM without spinup
julia> # land_time_series = runTEM(selected_models, forcing, nothing, loc_forcing_t, loc_land, tem_info)

source

Code

julia

function runTEMOne(selected_models, loc_forcing, land_init, tem, spinup_sequence)
    loc_forcing_t = getForcingForTimeStep(loc_forcing, loc_forcing, 1, tem.vals.forcing_types)
    loc_land = buildFullLand(selected_models, loc_forcing_t, land_init, tem.model_helpers)
    loc_land = addSpinupStore(loc_land, spinup_sequence, tem.run.store_spinup)
    checkSpinupInput(loc_land, loc_forcing_t, spinup_sequence)
    return loc_forcing_t, loc_land
end

runTEM!

Sindbad.Simulation.runTEM! Function

julia

runTEM!(selected_models, forcing::NamedTuple, info::NamedTuple)
runTEM!(forcing::NamedTuple, info::NamedTuple)
runTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_land, tem_info::NamedTuple)
runTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_restart, space_land, tem_info::NamedTuple)

Runs the SINDBAD Terrestrial Ecosystem Model (TEM) for all locations and time steps using preallocated arrays as the model data backend. This function supports multiple configurations for efficient execution.

Arguments:

For the first variant:

selected_models: A tuple of all models selected in the given model structure.
forcing::NamedTuple: A forcing NamedTuple containing the time series of environmental drivers for all locations.
info::NamedTuple: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.

For the second variant:

forcing::NamedTuple: A forcing NamedTuple containing the time series of environmental drivers for all locations.
info::NamedTuple: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.

For the third variant:

space_selected_models: A vector with the tuple or long tuple of the selected models of each location, from prepTEM.
space_forcing: A collection of forcing NamedTuples for multiple locations, replicated to avoid data races during parallel execution.
space_spinup: A collection of spinup NamedTuples for multiple locations, each with the sequence of that location and the forcing derived from it, replicated to avoid data races during parallel execution.
loc_forcing_t: A forcing NamedTuple for a single location and a single time step.
space_output: A collection of output arrays/views for multiple locations, replicated to avoid data races during parallel execution.
space_land: A collection of initial SINDBAD land NamedTuples for multiple locations, ensuring that the model states for one location do not overwrite those of another.
tem_info::NamedTuple: A helper NamedTuple containing necessary objects for model execution and type consistencies.

For the fourth variant:

Same as the third variant, plus space_restart, the per-location views of the restart arrays from prepTEM, where the pools after the spinup step with save_restart are stored. The third variant passes nothing for every location, so the optimization and cost runs store no restart data.

The third and fourth variants take the keyword show_progress, which shows a progress bar over the locations when it is true. It is false by default, so that the many runs of an optimization stay quiet.

Returns:

output_array: A preallocated array containing the simulation results for all locations and time steps.

Notes:

Precomputations:
- The function uses prepTEM to prepare the necessary inputs and configurations for the simulation.
Parallel Execution:
- The function uses parallelizeTEM! to distribute the simulation across multiple locations using the specified parallelization backend (Threads.@threads or qbmap).
Core Execution:
- For each location, the function calls coreTEM! to execute the TEM simulation, including spinup (if enabled) and the main time loop.
Data Safety:
- The function ensures data safety by replicating forcing, output, and land data for each location, avoiding data races during parallel execution.

Examples

julia

julia> using Sindbad

julia> # Run TEM with preallocated arrays (shorthand)
julia> # output_array = runTEM!(forcing, info)

julia> # Run TEM with preallocated arrays (detailed)
julia> # output_array = runTEM!(selected_models, forcing, info)

julia> # Run TEM with precomputed helpers
julia> # runTEM!(selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_land, tem_info)

source

Code

julia

function runTEM! end

function runTEM!(selected_models, forcing::NamedTuple, info::NamedTuple)
    run_helpers = prepTEM(selected_models, forcing, info)
    runTEM!(run_helpers.space_selected_models, run_helpers.space_forcing, run_helpers.space_spinup, run_helpers.loc_forcing_t, run_helpers.space_output, run_helpers.space_land, run_helpers.tem_info)
    return run_helpers.output_array
end

function runTEM!(selected_models, forcing::NamedTuple, info::NamedTuple)
    run_helpers = prepTEM(selected_models, forcing, info)
    runTEM!(run_helpers.space_selected_models, run_helpers.space_forcing, run_helpers.space_spinup, run_helpers.loc_forcing_t, run_helpers.space_output, run_helpers.space_land, run_helpers.tem_info)
    return run_helpers.output_array
end

function runTEM!(forcing::NamedTuple, info::NamedTuple)
    run_helpers = prepTEM(forcing, info)
    runTEM!(run_helpers.space_selected_models, run_helpers.space_forcing, run_helpers.space_spinup, run_helpers.loc_forcing_t, run_helpers.space_output, run_helpers.space_land, run_helpers.tem_info)
    return run_helpers.output_array
end

function runTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_land, tem_info::NamedTuple; show_progress=false)
    # one `nothing` per location, so no pools are stored for the restart file
    no_restart = Fill(nothing, length(space_forcing))
    return runTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, no_restart, space_land, tem_info; show_progress=show_progress)
end

function runTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_restart, space_land, tem_info::NamedTuple; show_progress=false)
    parallelizeTEM!(space_selected_models, space_forcing, space_spinup, loc_forcing_t, space_output, space_restart, space_land, tem_info, tem_info.run.parallelization; show_progress=show_progress)
    return nothing
end

runTEMLazy

Sindbad.Simulation.runTEMLazy Function

julia

runTEMLazy(xmap_cubes, store_path; max_cache=5e8, loop_ranges=nothing,
    run_filter=nothing, variable_properties=nothing,
    global_properties=Dict{String,Any}(), save_kwargs=(;),
    parallelization=ThreadsParallelization())

Runs the lazy run of prepTEMLazy and writes all its output and restart cubes in one pass to the zarr store at store_path.

Arguments:

xmap_cubes: the xmap field returned by prepTEMLazy
store_path: the path of the zarr store that holds the computed cubes
max_cache: the cache size for the block that DiskArrayEngine optimizes
loop_ranges: the steps of each loop dimension in a block, from getLoopRanges, or nothing to let DiskArrayEngine choose the block
run_filter: a ValidGridsFilter from getLoopRanges that skips the blocks without a valid grid, or nothing to run every block
variable_properties: a Dict of the attributes of each output variable by name, or nothing
global_properties: the global attributes of the store
save_kwargs: the keyword arguments of savedataset, e.g., the compressor
parallelization: how the blocks run, see getBlockRunner

Returns:

a NamedTuple with the output and restart Datasets, read from the store

Notes:

All cubes of prepTEMLazy come from one xmap. Saving them one by one runs the model once per cube, while runXmap runs it once for all of them.
The cubes keep their padded dimensions until they are computed, since the run can only fuse the direct outputs of the xmap. The store has them removed.
With the attributes, the store can be the final output file. The restart cubes stay in it until they are saved to their own file, see dropVariablesFromStore.

source

Code

julia

function runTEMLazy(xmap_cubes, store_path; max_cache=5e8, loop_ranges=nothing, run_filter=nothing, variable_properties=nothing, global_properties=Dict{String,Any}(), save_kwargs=(;), parallelization=ThreadsParallelization())
    dataset = xmap_cubes.dataset
    if !isnothing(variable_properties)
        cubes = Any[addPropsToCube(cube, get(variable_properties, name, nothing)) for (name, cube) ∈ dataset.cubes]
        dataset = datasetFromCubes(collect(keys(dataset.cubes)), cubes)
    end
    computed = runXmap(dataset, store_path, xmap_cubes.padded_names; max_cache=max_cache, loop_ranges=loop_ranges, run_filter=run_filter, properties=global_properties, save_kwargs=save_kwargs, parallelization=parallelization)
    output_ds = datasetFromCubes(xmap_cubes.output_names, Any[computed.cubes[name] for name ∈ xmap_cubes.output_names])
    restart_ds = datasetFromCubes(xmap_cubes.restart_names, Any[computed.cubes[key] for key ∈ xmap_cubes.restart_keys])
    return (; output=output_ds, restart=restart_ds)
end

setOutputForTimeStep!

Sindbad.Simulation.setOutputForTimeStep! Function

julia

setOutputForTimeStep!(outputs, land, ts, Val{output_vars})

Arguments:

outputs: vector of model output vectors
land: a core SINDBAD NT that contains all variables for a given time step that is overwritten at every timestep
ts: time step
::Val{output_vars}: a dispatch for vals of the output variables to generate the function

source

setupOptiOutput

Sindbad.Simulation.setupOptiOutput Function

julia

setupOptiOutput(info::NamedTuple, output::NamedTuple[, ::DoRunOptimization | ::DoNotRunOptimization])

Creates the output fields needed for the optimization experiment.

Arguments:

info: A SINDBAD NamedTuple containing all information needed for setup and execution of an experiment.
output: A base output NamedTuple to which optimization-specific fields will be added.
::DoRunOptimization: (Optional) A type dispatch indicating that optimization is enabled. Adds fields for optimized parameters.
::DoNotRunOptimization: (Optional) A type dispatch indicating that optimization is not enabled. Returns the base output unchanged.

Returns:

A NamedTuple containing the base output fields, with additional fields for optimization if enabled.

Notes:

When optimization is enabled, the function:
- Adds a parameter_dim field to the output, which includes the parameter dimension and metadata.
- Creates an OutDims object for storing optimized parameters, with the appropriate backend and file path.
When optimization is not enabled, the function simply returns the input output NamedTuple unchanged.

Examples:

With optimization enabled:

julia

output = setupOptiOutput(info, output, DoRunOptimization())

Without optimization:

julia

output = setupOptiOutput(info, output, DoNotRunOptimization())

source

Code

julia

function setupOptiOutput end

function setupOptiOutput(info::NamedTuple, output::NamedTuple, ::DoRunOptimization)
    @debug "     setupOptiOutput: getting parameter output for optimization..."
    params = info.optimization.parameter_table.name_full    
    paramaxis = YAXArrays.Dim{:parameter}(params)
    outformat = info.output.format
    backend = outformat == "nc" ? :netcdf : :zarr
    od = YAXArrays.XOutput(paramaxis)
        #path=joinpath(info.output.dirs.optimization,
        #    "optimized_parameters.$(outformat)"),
        #backend=backend,
        #overwrite=true)
    # list of parameter
    output = set_namedtuple_field(output, (:parameter_dim, od))
    return output
end

function setupOptiOutput(info::NamedTuple, output::NamedTuple, ::DoRunOptimization)
    @debug "     setupOptiOutput: getting parameter output for optimization..."
    params = info.optimization.parameter_table.name_full    
    paramaxis = YAXArrays.Dim{:parameter}(params)
    outformat = info.output.format
    backend = outformat == "nc" ? :netcdf : :zarr
    od = YAXArrays.XOutput(paramaxis)
        #path=joinpath(info.output.dirs.optimization,
        #    "optimized_parameters.$(outformat)"),
        #backend=backend,
        #overwrite=true)
    # list of parameter
    output = set_namedtuple_field(output, (:parameter_dim, od))
    return output
end

function setupOptiOutput(info::NamedTuple, output::NamedTuple, ::DoNotRunOptimization)
    return output
end

spinup

Sindbad.Simulation.spinup Function

julia

spinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, spinup_mode::SpinupMode)

Runs the spinup process for the SINDBAD Terrestrial Ecosystem Model (TEM) to initialize the model to a steady state. The spinup process updates the state variables (e.g., pools) using various spinup methods.

Arguments:

spinup_models: A tuple of a subset of all models in the given model structure that are selected for spinup.
spinup_forcing: A forcing NamedTuple containing the time series of environmental drivers for the spinup process.
loc_forcing_t: A forcing NamedTuple for a single location and a single time step.
land: A SINDBAD NamedTuple containing all variables for a given time step, which is overwritten at every timestep.
tem_info: A helper NamedTuple containing necessary objects for model execution and type consistencies.
n_timesteps: The number of timesteps for the spinup process.
spinup_mode::SpinupMode: A type dispatch that determines the spinup method to be used.

Returns:

land: The updated SINDBAD NamedTuple containing the final state of the model after the spinup process.

SpinupMode

Abstract type for model spinup modes in SINDBAD

Available methods/subtypes:

ClearanceA0: cut the vegetation carbon pools down to c_remain, as EtaScaleA0H does without the scaling by ηH
EtaScaleA0H: scale carbon pools using diagnostic scalars for ηH and c_remain
EtaScaleA0HCWD: scale carbon pools of CWD (cLitSlow) using ηH and set vegetation pools to c_remain
EtaScaleAH: scale carbon pools using diagnostic scalars for ηH and ηA
EtaScaleAHCWD: scale carbon pools of CWD (cLitSlow) using ηH and scale vegetation pools by ηA
EtaScaleH: scale only the heterotrophic carbon pools using diagnostic scalars for ηH
LoopTEM: run the TEM over the spinup forcing with the models of the step
NlsolveFixedpointTrustregionCEco: use a fixed-point nonlinear solver with trust region for carbon pools (cEco)
NlsolveFixedpointTrustregionCEcoTWS: use a fixed-point nonlinear solver with trust region for both cEco and TWS
NlsolveFixedpointTrustregionTWS: use a fixed-point nonlinearsolver with trust region for Total Water Storage (TWS)
ODEAutoTsit5Rodas5: use the AutoVern7(Rodas5) method from DifferentialEquations.jl for solving ODEs
ODEDP5: use the DP5 method from DifferentialEquations.jl for solving ODEs
ODETsit5: use the Tsit5 method from DifferentialEquations.jl for solving ODEs
SSPDynamicSSTsit5: use the SteadyState solver with DynamicSS and Tsit5 methods
SSPSSRootfind: use the SteadyState solver with SSRootfind method
SpinupInput: replace the pools with the forcing variables of the same name, e.g., the main pools of a restart file, and fill their component pools
Spinup_TWS: Spinup spinup_mode for Total Water Storage (TWS)
Spinup_cEco: Spinup spinup_mode for cEco
Spinup_cEco_TWS: Spinup spinup_mode for cEco and TWS

Extended help

Notes:

The spinup process can use different methods depending on the spinup_mode, including fixed-point solvers, ODE solvers, and steady-state solvers.
The function dynamically selects the appropriate spinup method based on the spinup_mode dispatch type.
For ODE-based methods, the function uses DifferentialEquations.jl to solve the spinup equations.
For steady-state solvers, the function uses methods like DynamicSS or SSRootfind to find equilibrium states.

Examples:

Running spinup with selected models:

julia

land = spinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, LoopTEM())

Running spinup with ODE solver (Tsit5):

julia

land = spinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, ODETsit5())

Running spinup with fixed-point solver for cEco and TWS:

julia

land = spinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, NlsolveFixedpointTrustregionCEcoTWS())

Running spinup with steady-state solver (SSRootfind):

julia

land = spinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, SSPSSRootfind())

source

Code

julia

function spinupTEMLazy end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoSpinupTEM, store_restart=nothing)
    loc_spinup = getLocSpinup(loc_forcing, tem_info)
    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM, store_restart)
    return land_spin
end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoSpinupTEM, store_restart=nothing)
    loc_spinup = getLocSpinup(loc_forcing, tem_info)
    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM, store_restart)
    return land_spin
end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoNotSpinupTEM, store_restart=nothing)
    setRestart!(store_restart, land_prec, tem_info)
    return land_prec
end

spinupTEM

Sindbad.Simulation.spinupTEM Function

julia

spinupTEM(selected_models, loc_spinup, loc_forcing_t, land, tem_info, spinup_mode, loc_restart=nothing)

The main spinup function that handles the spinup method based on inputs from spinup.json. Either the spinup is loaded or/and run using spinup functions for different spinup methods.

Arguments:

selected_models: a tuple of all models selected in the given model structure
loc_spinup: a NT with the spinup sequence of the location and the forcing derived from it
loc_forcing_t: a forcing NT for a single location and a single time step
land: SINDBAD NT input to the spinup of TEM during which subfield(s) of pools are overwritten
tem_info: helper NT with necessary objects for model run and type consistencies
loc_restart: where the pools are stored for the restart file, after the step of the sequence whose save_restart is true. With nothing, the default, nothing is stored. Without spinup, the pools of the input land are stored.
spinup_mode: A type dispatch that determines whether spinup is included or excluded:
- ::DoSpinupTEM: Runs the spinup process before the main simulation. Set spinup_TEM to true in the flag section of experiment_json.
- ::DoNotSpinupTEM: Skips the spinup process and directly runs the main simulation. Set spinup_TEM to false in the flag section of experiment_json.

Examples

julia

julia> using Sindbad

julia> # Run spinup with DoSpinupTEM
julia> # land_spinup = spinupTEM(selected_models, forcing, loc_forcing_t, land, tem_info, DoSpinupTEM())

julia> # Skip spinup with DoNotSpinupTEM
julia> # land = spinupTEM(selected_models, forcing, loc_forcing_t, land, tem_info, DoNotSpinupTEM())

Notes:

When DoSpinupTEM is used:
- The function runs the spinup sequences
When DoNotSpinupTEM is used:
- The function skips the spinup process returns the land as is`

source

Code

julia

function spinupTEMLazy end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoSpinupTEM, store_restart=nothing)
    loc_spinup = getLocSpinup(loc_forcing, tem_info)
    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM, store_restart)
    return land_spin
end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoSpinupTEM, store_restart=nothing)
    loc_spinup = getLocSpinup(loc_forcing, tem_info)
    land_spin = spinupTEM(selected_models, loc_spinup, loc_forcing_t, land_prec, tem_info, tem_info.run.spinup_TEM, store_restart)
    return land_spin
end

function spinupTEMLazy(selected_models, loc_forcing, loc_forcing_t, land_prec, tem_info, ::DoNotSpinupTEM, store_restart=nothing)
    setRestart!(store_restart, land_prec, tem_info)
    return land_prec
end

timeLoopTEMSpinup

Sindbad.Simulation.timeLoopTEMSpinup Function

julia

timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps)

do/run the time loop of the spinup models to update the pool. Note that, in this function, the time series is not stored and the land/land is overwritten with every iteration. Only the state at the end is returned

Arguments:

spinup_models: a tuple of a subset of all models in the given model structure that is selected for spinup
spinup_forcing: a selected/sliced/computed forcing time series for running the spinup sequence for a location
loc_forcing_t: a forcing NT for a single location and a single time step
land: SINDBAD NT input to the spinup of TEM during which subfield(s) of pools are overwritten
tem_info: helper NT with necessary objects for model run and type consistencies
n_timesteps: number of time steps

source

julia

timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, store_spinup_mode)

Runs the time loop of the spinup models like the method without store_spinup_mode, and with DoStoreSpinup writes cEco and TWS to the spinup store after every time step.

Notes:

Only spinup with LoopTEM calls this method. The solver based spinup modes call the method without store_spinup_mode, so their trial runs are not stored.

source

Code

julia

function timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps)
    for ts ∈ 1:n_timesteps
        f_ts = getForcingForTimeStep(spinup_forcing, loc_forcing_t, ts, tem_info.vals.forcing_types)
        land = computeTEM(spinup_models, f_ts, land, tem_info.model_helpers)
    end
    return land
end

function timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, ::DoNotStoreSpinup)
    return timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps)
end

function timeLoopTEMSpinup(spinup_models, spinup_forcing, loc_forcing_t, land, tem_info, n_timesteps, store_spinup_mode::DoStoreSpinup)
    for ts ∈ 1:n_timesteps
        f_ts = getForcingForTimeStep(spinup_forcing, loc_forcing_t, ts, tem_info.vals.forcing_types)
        land = computeTEM(spinup_models, f_ts, land, tem_info.model_helpers)
        storeSpinupStep!(land, store_spinup_mode)
    end
    return land
end

TEM

Modules

Extensions

Functions

TEMLazy

accumulateSpinupForcing

coreTEM

coreTEM!

dropVariablesFromStore

getAllSpinupForcing

getForcingForTimeStep

getLocData

getLoopRanges

getOutDims

getOutDimsArrays

getRestartDimsArrays

getSpatialInfo

getSpinupSequence

getYaxMaxCache

prepSequenceDefault

prepSequenceWithAge

prepSpinupOut

prepTEM

prepTEMLazy

prepTEMOut

runSpinupSequences

runTEM

runTEM!

runTEMLazy

setOutputForTimeStep!

setupOptiOutput

spinup

spinupTEM

timeLoopTEMSpinup

Functions ​

TEMLazy ​

accumulateSpinupForcing ​

coreTEM ​

coreTEM! ​

dropVariablesFromStore ​

getAllSpinupForcing ​

getForcingForTimeStep ​

getLocData ​

getLoopRanges ​

getOutDims ​

getOutDimsArrays ​

getRestartDimsArrays ​

getSpatialInfo ​

getSpinupSequence ​

getYaxMaxCache ​

prepSequenceDefault ​

prepSequenceWithAge ​

prepSpinupOut ​

prepTEM ​

prepTEMLazy ​

prepTEMOut ​

runSpinupSequences ​

runTEM ​

runTEM! ​

runTEMLazy ​

setOutputForTimeStep! ​

setupOptiOutput ​

spinup ​

spinupTEM ​

timeLoopTEMSpinup ​

Functions

TEMLazy

accumulateSpinupForcing

coreTEM

coreTEM!

dropVariablesFromStore

getAllSpinupForcing

getForcingForTimeStep

getLocData

getLoopRanges

getOutDims

getOutDimsArrays

getRestartDimsArrays

getSpatialInfo

getSpinupSequence

getYaxMaxCache

prepSequenceDefault

prepSequenceWithAge

prepSpinupOut

prepTEM

prepTEMLazy

prepTEMOut

runSpinupSequences

runTEM

runTEM!

runTEMLazy

setOutputForTimeStep!

setupOptiOutput

spinup

spinupTEM

timeLoopTEMSpinup