OceananigansLagrangianFilter.jl Library

OnlineFilterConfig

The configuration of an online Lagrangian filter, made with the keyword constructor below.

OnlineFilterConfig(; grid::AbstractGrid,
                        output_filename::String = "online_filtered_output.jld2",
                        var_names_to_filter::Tuple{Vararg{String}},
                        velocity_names::Tuple{Vararg{String}},
                        N::Union{Int, Nothing} = nothing,
                        freq_c::Union{Real, Nothing} = nothing,
                        filter_params::Union{NamedTuple, Nothing} = nothing,
                        compute_maps::Bool = true,
                        regrid_to_mean::Bool = true,
                        compute_mean_velocities::Bool = true,
                        lagrangian::Bool = true,
                        label::String = "",
                        relaxation_mask_func::Union{Function, Nothing} = nothing,
                        relaxation_params::Union{NamedTuple, Nothing} = nothing,
                        frequency_scale_func::Union{Function, Nothing} = nothing,
                        frequency_scale_params::Union{NamedTuple, Nothing} = nothing,
                        frequency_scale_dependencies::Tuple{Vararg{Symbol}} = ()

                        )

Constructs a configuration object for online Lagrangian filtering of Oceananigans data. This function validates the time specifications and filter parameters before creating the OnlineFilterConfig object.

Keyword arguments

• grid: (required) The grid of the model.
• output_filename: The filename for the output of the online filtered data. Default: "online_filtered_output.jld2".
• var_names_to_filter: (required) A Tuple of Strings specifying the names of the tracer variables to be filtered.
• velocity_names: (required) A Tuple of Strings specifying the names of the model's velocity fields to be used for advection.
• N, freq_c: Parameters for a Butterworth filter. N is the order of the filter, and freq_c is the cutoff frequency. These are used to automatically generate filter_params if not provided. Must be specified together if filter_params is not given.
• filter_params: A NamedTuple containing the coefficients for a custom filter. Only filterparams OR N and `freqc` should be given.
• compute_maps: A Bool indicating whether to solve for and output the maps (displacements from the mean position). Default: true.
• regrid_to_mean: A Bool indicating whether the filtered fields should be interpolated to the mean position after the simulation (i.e. calculate the generalised Lagrangian mean), using regrid_to_mean_position!. This requires the maps, so compute_maps is set to true if needed. Only possible with normalised filter coefficients on a RectilinearGrid or LatitudeLongitudeGrid, and set to false otherwise. On a LatitudeLongitudeGrid it assumes the grid is locally Cartesian over the trajectory excursions (see Regridding to the mean position). Default: true.
• compute_mean_velocities: A Bool indicating whether to compute and output the mean velocities (from the maps). Default: true.
• lagrangian: true for the Lagrangian filter, or false for the Eulerian filter, which filters at fixed points: the filtered variables mustn't be advected (set the model's tracer_advection with filtered_vars_advection), and the outputs are named _Eulerian_filtered. The maps and mean velocities can still be computed, but not regridded. Default: true.
• label: A String label for the variables that will be created to pass to the model. For use when multiple filter configurations are to be run at the same time. Default: "".
• relaxation_mask_func: A Function defining the mask for relaxing the filtered variables towards their local equilibrium (see Boundary relaxation), e.g. near the boundaries of the domain. Should be 1 for full relaxation, and 0 for no relaxation. Arguments should be the non-flat spatial coordinates and relaxation_params. Relaxation is only included if this is given. Default nothing.
• relaxation_params: A NamedTuple of parameters for relaxation_mask_func, which must include the relaxation timescale. Default nothing.
• frequency_scale_func: A Function giving the frequency scale $M$, which locally multiplies the filter's frequencies (e.g. the cutoff frequency), so that the filter can vary in space, and in time through frequency_scale_dependencies. Must be positive everywhere throughout the run (with frequency_scale_dependencies, it is only checked for the initial fields), and GPU-compatible if running on a GPU. Arguments should be the non-flat spatial coordinates, the values of the frequency_scale_dependencies, and frequency_scale_params. The filter is the same everywhere if this is nothing. See Varying cutoff frequency. Default nothing.
• frequency_scale_params: A NamedTuple of parameters for frequency_scale_func. Default nothing.
• frequency_scale_dependencies: A Tuple of Symbols naming fields at cell centres whose values are passed to frequency_scale_func, e.g. (:zeta,). They can't be the velocities: compute what the frequency scale needs (e.g. the vorticity) as a field at cell centres instead. Offline, they must be in var_names_to_filter or extra_auxiliary_var_names. Online, they can be any field of the model, e.g. one of its auxiliary_fields. Default ().

Example:

using OceananigansLagrangianFilter
using Oceananigans.Units

Nx = 50
Nz = 20
L = 10kilometers 
H = 100meters 

grid = RectilinearGrid(CPU(),size = (Nx, Nz), 
                       x = (-L/2, L/2),
                       z = (-H, 0),
                       topology = (Periodic, Flat, Bounded))

filter_config = OnlineFilterConfig( grid = grid,
                                    output_filename = "test_filter.jld2",
                                    var_names_to_filter = ("b","T"), 
                                    velocity_names = ("u","w"),
                                    N = 2,
                                    freq_c = 1e-4/2)

# output
┌ Info: Advection for Lagrangian filtering will be performed using full model velocities u, v, and w.
└         Maps (displacements from the mean position) will be computed corresponding to velocities: ("u", "w").
[ Info: Mean velocities corresponding to ("u", "w") will be computed.
[ Info: Variables to be filtered: ("b", "T"). Ensure these are valid tracer or auxiliary field names in the simulation.
[ Info: Setting filter parameters to use Butterworth order 2, cutoff frequency 5.0e-5
OnlineFilterConfig(50×1×20 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [-5000.0, 5000.0) regularly spaced with Δx=200.0
├── Flat y
└── Bounded  z ∈ [-100.0, 0.0]     regularly spaced with Δz=5.0, "test_filter.jld2", ("b", "T"), ("u", "w"), (a1 = 1.421067568548072e-20, b1 = -7.071067811865475e-5, c1 = 3.535533905932738e-5, d1 = -3.535533905932738e-5, N_coeffs = 1, outputs = :combined), true, true, true, true, "", nothing, nothing, nothing, nothing, ())

OfflineFilterConfig

The configuration of an offline Lagrangian filter, run with run_offline_Lagrangian_filter, made with the keyword constructor below.

OfflineFilterConfig(; original_data_filename::String,
                    var_names_to_filter::Tuple{Vararg{String}},
                    velocity_names::Tuple{Vararg{String}},
                    extra_auxiliary_var_names::Tuple{Vararg{String}} = (),
                    T_start::Union{Real,Nothing} = nothing,
                    T_end::Union{Real,Nothing} = nothing,
                    T::Union{Real,Nothing} = nothing,
                    architecture::AbstractArchitecture = CPU(),
                    T_out::Union{Real,Nothing} = nothing,
                    N::Union{Int, Nothing} = nothing,
                    freq_c::Union{Real, Nothing} = nothing,
                    filter_params::Union{NamedTuple, Nothing} = nothing,
                    Δt::Union{Real,Nothing} = nothing,
                    compute_maps::Bool = true,
                    regrid_to_mean::Bool = true,
                    forward_output_filename::String = "forward_output.jld2",
                    backward_output_filename::String = "backward_output.jld2",
                    output_filename::String = "filtered_output.jld2",
                    compute_mean_velocities::Bool = true,
                    delete_intermediate_files::Bool = true,
                    compute_Eulerian_filter::Bool = false,
                    output_netcdf::Bool = false,
                    output_original_data::Bool = true,
                    lagrangian::Bool = true,
                    advection::AbstractAdvectionScheme = WENO(),
                    grid::Union{AbstractGrid, Nothing} = nothing,
                    label::String = "",
                    relaxation_mask_func::Union{Function, Nothing} = nothing,
                    relaxation_params::Union{NamedTuple, Nothing} = nothing,
                    frequency_scale_func::Union{Function, Nothing} = nothing,
                    frequency_scale_params::Union{NamedTuple, Nothing} = nothing,
                    frequency_scale_dependencies::Tuple{Vararg{Symbol}} = ())

Constructs a configuration object for offline Lagrangian filtering of Oceananigans data. This function validates the input data file, time specifications, and filter parameters before creating the OfflineFilterConfig object.

Keyword arguments

• original_data_filename: (required) The path to the JLD2 or NetCDF file containing the original Oceananigans output data.
• var_names_to_filter: (required) A Tuple of Strings specifying the names of the tracer variables to be filtered.
• velocity_names: (required) A Tuple of Strings specifying the names of the velocity fields in the data file to be used for advection.
• extra_auxiliary_var_names: A Tuple of Strings naming extra variables in the data file to read in as auxiliary fields of the filter model, without filtering them, e.g. for use in frequency_scale_func. They are interpolated in time like the filtered variables, but are not reversed in sign on the backward pass. Default: ().
• T_start: Start time for the filter. Must be within the time range of the data. If not given, defaults to either T_end - T (if they are given), or the start time of the original data.
• T_end: End time for the filter. Must be within the time range of the data. If not given, defaults to either Tstart + T (after Tstart given or computed, if T is given), or the end time of the original data.
• T: Duration of the filtering. If not given, defaults to Tend - Tstart (after Tstart and Tend are given or computed).
• architecture: The architecture (CPU or GPU) to be used for the filtering computation. Default: CPU().
• T_out: The output time step for the filtered data. If nothing, it defaults to the time step of the original data.
• N, freq_c: Parameters for a Butterworth squared filter. N is the order of the filter, and freq_c is the cutoff frequency. These are used to automatically generate filter_params if not provided. Must be specified together if filter_params is not given.
• filter_params: A NamedTuple containing the coefficients for a custom filter. Only filterparams OR N and `freqc` should be given.
• Δt: The time step for the internal Lagrangian filter simulation. If nothing, it defaults to T_out / 10, but this may not be appropriate.
• compute_maps: A Bool indicating whether to solve for and output the maps (displacements from the mean position). Default: true.
• regrid_to_mean: A Bool indicating whether to interpolate the filtered fields to the mean position (i.e. calculate the generalised Lagrangian mean). This requires the maps, so compute_maps is set to true if needed. Only possible for a Lagrangian filter with normalised filter coefficients on a RectilinearGrid or LatitudeLongitudeGrid, and set to false otherwise. On a LatitudeLongitudeGrid it assumes the grid is locally Cartesian over the trajectory excursions (see Regridding to the mean position). Default: true.
• forward_output_filename: The filename for the output of the forward filter pass. Default: "forward_output.jld2".
• backward_output_filename: The filename for the output of the backward filter pass. Default: "backward_output.jld2".
• output_filename: The filename for the final combined and mapped output. Default: "filtered_output.jld2".
• compute_mean_velocities: A Bool indicating whether to compute and output the mean velocities (from the maps). Default: true.
• delete_intermediate_files: A Bool indicating whether to delete forward_output.jld2 and backward_output.jld2 after the final combined file is created. Default: true.
• compute_Eulerian_filter: A Bool indicating whether to also compute an Eulerian-mean-based filter for comparison. Default: false.
• output_netcdf: A Bool indicating whether to also convert the final JLD2 output file to a NetCDF file. Default: false.
• output_original_data: A Bool indicating whether to include the original data in the final output file for comparison. Default: true.
• lagrangian: true for the Lagrangian filter, or false for the Eulerian filter, which filters at fixed points: the filtered variables aren't advected, and the outputs are named _Eulerian_filtered. The maps and mean velocities can still be computed, but not regridded. Default: true.
• advection: The advection scheme to use for the Lagrangian filter simulation (not used if lagrangian = false). Default: WENO(). Using lower-order schemes may be a source of error.
• grid: The grid for the simulation. If nothing, the grid is inferred from the original_data_filename (preferred option)
• label: A String label for the variables that will be created to pass to the model. For use when multiple filter configurations are to be run at the same time. Default: "".
• relaxation_mask_func: A Function defining the mask for relaxing the filtered variables towards their local equilibrium (see Boundary relaxation), e.g. near the boundaries of the domain. Should be 1 for full relaxation, and 0 for no relaxation. Arguments should be the non-flat spatial coordinates and relaxation_params. Relaxation is only included if this is given. Default nothing.
• relaxation_params: A NamedTuple of parameters for relaxation_mask_func, which must include the relaxation timescale. Default nothing.
• frequency_scale_func: A Function giving the frequency scale $M$, which locally multiplies the filter's frequencies (e.g. the cutoff frequency), so that the filter can vary in space, and in time through frequency_scale_dependencies. Must be positive everywhere throughout the run (with frequency_scale_dependencies, it is only checked for the initial fields), and GPU-compatible if running on a GPU. Arguments should be the non-flat spatial coordinates, the values of the frequency_scale_dependencies, and frequency_scale_params. The filter is the same everywhere if this is nothing. See Varying cutoff frequency. Default nothing.
• frequency_scale_params: A NamedTuple of parameters for frequency_scale_func. Default nothing.
• frequency_scale_dependencies: A Tuple of Symbols naming fields at cell centres whose values are passed to frequency_scale_func, e.g. (:zeta,). They can't be the velocities: compute what the frequency scale needs (e.g. the vorticity) as a field at cell centres instead. Offline, they must be in var_names_to_filter or extra_auxiliary_var_names. Online, they can be any field of the model, e.g. one of its auxiliary_fields. Default ().

Example:

using OceananigansLagrangianFilter
using Oceananigans.Units
path_to_sim = "../test/data/reference_sim.jld2"
filter_config = OfflineFilterConfig(original_data_filename=path_to_sim, 
                                    output_filename = "output_file.jld2", 
                                    var_names_to_filter = ("b",), 
                                    velocity_names = ("u","w"), 
                                    architecture = CPU(), 
                                    Δt = 20minutes, 
                                    T_out = 1hour, 
                                    N = 2, 
                                    freq_c = 1e-4/2, 
                                    compute_mean_velocities = true, 
                                    output_netcdf = true, 
                                    delete_intermediate_files = true, 
                                    compute_Eulerian_filter = true) 

# output
┌ Info: Advection for Lagrangian filtering will be performed using only velocities ("u", "w") -
│ any other velocity components will be zero by default. Maps (displacements from the mean position) will
└ be computed corresponding to velocities: ("u", "w").
[ Info: Mean velocities corresponding to ("u", "w") will be computed.
[ Info: Filter interval will be from T_start=0.0 to T_end=86400.0, duration T=86400.0
[ Info: Setting filter parameters to use Butterworth squared, order 2, cutoff frequency 5.0e-5
OfflineFilterConfig("../test/data/reference_sim.jld2", ("b",), ("u", "w"), (), 0.0, 86400.0, 86400.0, CPU(), 3600.0, (a1 = 1.767766952966369e-5, b1 = 1.767766952966369e-5, c1 = 3.535533905932738e-5, d1 = 3.535533905932738e-5, N_coeffs = 1, outputs = :combined, sine_parity = :even), 1200.0, true, true, "forward_output.jld2", "backward_output.jld2", "output_file.jld2", true, true, true, true, true, true, WENO{3, Float64, Nothing}(order=5)
├── buffer_scheme: WENO{2, Float64, Nothing}(order=3)
│   └── buffer_scheme: Centered(order=2)
└── advecting_velocity_scheme: Centered(order=2), 10×1×10 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [-5000.0, 5000.0) regularly spaced with Δx=1000.0
├── Flat y
└── Bounded  z ∈ [-100.0, 0.0]     regularly spaced with Δz=10.0, "", nothing, nothing, nothing, nothing, ())

run_offline_Lagrangian_filter(config::OfflineFilterConfig; batch_size::Union{Int, Nothing} = nothing)

Runs an offline Lagrangian filter on an Oceananigans FieldTimeSeries dataset as configured by config.

This function performs a series of steps to filter the data:

1. Run forward simulation: A LagrangianFilter model is created and run forward in time to compute the first half of the filter contributions. Input data is read directly from the original source file via a BufferedDataReader — no intermediate file is written.
2. Run backward simulation: A second BufferedDataReader is created for the backward pass and the simulation is run again. Velocity reversal is handled on-the-fly during interpolation.
3. Combine results: The forward and backward simulation outputs are summed to produce the final filtered data.
4. Post-processing: Optional post-processing steps are performed, including regridding the data to the mean position, computing a comparative Eulerian filter, and converting the output file to NetCDF.

The filter can also be run in batches of batch_size terms (each term of the filter, e.g. each frequency of a spectrum filter, has a cosine and a sine tracer for each filtered variable and map), one after the other. Each batch runs steps 1-3 with a smaller model, writing to its own file (<output_filename>_batch<n>.jld2), which is merged into output_filename as soon as the batch is done and then deleted (unless delete_intermediate_files = false). Outputs of separate terms (outputs = :separate) are numbered as in the whole filter, and combined outputs are the sums of those of the batches. The output is then the same as without batches (combined outputs to the Float32 precision they are saved with), and the post-processing in step 4 is done once on the merged output. Batches of the same size have the same model types, so only the first is compiled (and the last, if it is smaller). Batches help when a filter has many terms, as the compilation time grows faster than the number of terms, and a smaller model needs less memory.

Arguments

• config: An instance of OfflineFilterConfig that specifies all parameters and file paths for the filtering process.
• batch_size: The number of filter terms in each batch, or nothing (default) to run all the terms together.

AbstractDataSource

Supertype for objects that provide frame-by-frame read access to simulation output, without requiring an intermediate copy of the data to be written to disk.

Subtypes must implement:

• stored_times(source) → Vector{Float64} of ascending physical times in the Tstart..Tend window
• read_frame!(cpu_field::Field, source, varname::String, idx::Int) → fills cpu_field with frame idx

BufferedDataReader

Reads simulation output frame-by-frame from disk into a two-slot GPU buffer, interpolating in time on demand. Data is read directly from the original file, with no intermediate file.

Two buffer slots each hold a full set of GPU Fields (one per variable). As the simulation advances, the slot holding the older frame is recycled: the swap is just an integer flip (lo_slot = 3 - lo_slot) with zero allocation.

Forward direction : physical time = Tstart + simt (buffer scans file forward) Backward direction : physical time = Tend - simt (buffer scans file backward; velocity signs are negated during interpolation)

JLD2DataSource

Wraps a JLD2 file of Oceananigans output, providing indexed frame reads without writing any intermediate files.

Fields

• filename : path to the JLD2 file
• var_names : tracer/auxiliary variable names to read
• vel_names : velocity variable names to read
• stored_times : physical times for the Tstart..Tend window (ascending)
• stored_iterations : corresponding JLD2 iteration keys

JLD2DataSource(filename, var_names, vel_names; T_start=nothing, T_end=nothing)

Open filename and build a source covering the time window [T_start, T_end].

NetCDFDataSource

Wraps a NetCDF file of Oceananigans (or compatible) output, providing indexed frame reads without writing any intermediate files.

Fields

• filename : path to the NetCDF file
• var_names : tracer/auxiliary variable names to read
• vel_names : velocity variable names to read
• stored_times : physical times for the Tstart..Tend window (ascending)
• stored_indices : corresponding indices into the NetCDF time dimension
• time_name : name of the time coordinate in the file (default "time")

NetCDFDataSource(filename, var_names, vel_names;
                 T_start=nothing, T_end=nothing, time_name="time")

Open filename and build a source covering the time window [T_start, T_end].

advance_buffer!(reader, sim_t)

Ensure the two buffer slots bracket sim_t. Reuses whichever frame is still valid (forward or backward slide) so at most one new frame is read per call.

create_buffered_reader(config; direction=:forward, locations=Dict(), time_name="time")

Build a BufferedDataReader for the data file named in config, covering the time window [config.T_start, config.T_end].

Keyword arguments

• direction : :forward (default) or :backward
• locations : Dict{String, NTuple{3, DataType}} mapping variable names to their spatial location, e.g. Dict("u" => (Face,Center,Center)). Only required for NetCDF sources (JLD2 sources infer location from file metadata). Otherwise the location is read from the variable's attributes in the file, as written by Oceananigans' NetCDFWriter.
• time_name : name of the time coordinate in a NetCDF file (default "time")

create_data_source(config; time_name="time")

The data source (a JLD2DataSource, or a NetCDFDataSource for a file ending in .nc or .netcdf) for the original data file named in config, covering the time window [config.T_start, config.T_end], for the variables to filter, config.extra_auxiliary_var_names and the velocities. time_name is the name of the time coordinate in a NetCDF file.

interpolate_to_model!(model, reader, sim_t)

Linearly interpolate the two buffered frames at sim_t and write the result directly into model.velocities and model.auxiliary_fields, including the halo cells shared by the model and saved data grids (see shared_halo_regions).

Velocities are negated for backward-direction filtering.

update_input_data!(model, reader::BufferedDataReader)

Callback that reads directly from the source file via a two-frame GPU buffer. Calls advance_buffer! to slide the window if needed, then linearly interpolates into the model's velocity and auxiliary fields. Designed to be used at callsite UpdateStateCallsite(), so that the fields are updated at each substep of multi-stage time steppers.

compute_Eulerian_filter!(config::AbstractConfig)

Computes the Eulerian filter for specified variables and writes the results to a combined output file. At each grid point, the unfiltered fields are filtered in time with the filter's weight function. Each term of the weight function is a decaying complex exponential, so the filtered fields at all the output times are built up as running sums: forwards in time over the earlier times, and offline also backwards over the later times. Each field is read once in each direction, so the cost grows linearly with the number of saved times. Offline, the sums are kept in a temporary file next to the output file until they are written. The results are saved to the same JLD2 file in a new group with the _Eulerian_filtered suffix, at the output times. The running sums are computed on the filter's architecture (the offline architecture, or that of the online model grid), so on a GPU if the filter runs on one.

Offline, the unfiltered fields are read from the original data file, as the Lagrangian filter reads them: at the saved times between T_start and T_end, and linearly interpolated between them, so the result doesn't depend on the output interval T_out. Online, the unfiltered fields are read from the output file, so they must be added to the outputs (see the online example), and saved often enough to resolve the motions being filtered out; otherwise, run the filter with lagrangian = false.

This method serves as a benchmark for comparison with the main Lagrangian filter.

If the filter outputs its terms separately (outputs = :separate), each term is Eulerian filtered with its own weight function, and named like the Lagrangian outputs with _Eulerian_filtered in the stem, e.g. b_Eulerian_filtered_C1_scaled and b_Eulerian_filtered_S1_scaled.

The ends of the data are treated as by the Lagrangian filter, so that the two can be compared: each field is taken to be constant before the start (see initialise_filtered_vars_from_model), and offline, the sums over later times start at the end from the sums over earlier times, as the backward pass of the offline filter starts from the end of the forward pass. Between saved times each field is taken to vary linearly, and the weight function is integrated exactly, so a field that is constant in time is unchanged by a normalised filter.

With a frequency scale $M$ that only depends on the coordinates (no frequency_scale_dependencies), the weight function at each point is $M G(M \tau)$, with $M$ evaluated where each variable is saved (e.g. at the u points for u). A frequency scale with frequency_scale_dependencies varies in time, which isn't supported here: for the Eulerian filter with that frequency scale, run the filter with lagrangian = false instead.

Arguments

• config: An instance of AbstractConfig containing the file path, variable names, and filter parameters.

compute_time_shift!(config::AbstractConfig)

Computes the time shift for a forward filter (either online or forward-only offline) based on its coefficients and writes the shifted time series to the output file.

The time shift is computed as the time delay introduced by the filter's transfer function. This new time series is stored in a new group called timeseries/t_shifted within the output JLD2 file.

Arguments

• config: A configuration object of type AbstractFilterConfig which contains the output_filename and filter_params (filter coefficients).

create_filtered_vars(config::AbstractConfig)

Creates a Tuple of Symbols representing the names of the filtered tracer variables.

• For a single-exponential filter (N_coeffs = 0.5), the function generates names with a _C1 suffix.
• For a multi-coefficient filter (N_coeffs > 0.5), it generates pairs of names for each coefficient, suffixed with _C# and _S#, where # is the coefficient index.

If compute_maps or compute_mean_velocities is enabled in the configuration, additional symbols are created for the spatial mapping variables corresponding to each velocity component, prefixed with xi_ and suffixed with the corresponding coefficient names.

Arguments

• config: An instance of AbstractConfig containing the names of the variables to filter, the filter parameters, and the compute_maps and compute_mean_velocities booleans.

Returns

A Tuple of Symbols representing the names of the filtered variables to be used as tracers in the simulation.

create_forcing(filtered_vars::Tuple{Vararg{Symbol}}, config::AbstractConfig)

Create the forcing for the filtered tracers named in filtered_vars (from create_filtered_vars): the right-hand sides of the filter equations for each filtered variable, and for the map variables if they are needed, plus the boundary relaxation if config.relaxation_mask_func is given. The filter's frequencies are multiplied by the frequency scale if config.frequency_scale_func is given.

Returns a NamedTuple with one DiscreteForcing for each filtered tracer.

create_original_vars(config::AbstractConfig)

Creates a NamedTuple to serve as auxiliary fields for the original variables in a simulation. The fields are initialised from the input data to ensure that they are in the correct location. - the value of the data is unimportant.

Arguments

• config: An instance of AbstractConfig containing the names of the variables and the simulation grid.

Returns

A NamedTuple where each key is a Symbol of a variable name to be filtered, or of an extra auxiliary variable (config.extra_auxiliary_var_names), and each value is a field for that variable.

create_output_fields(model::AbstractModel, config::AbstractConfig; direction::Symbol = :forward)

Reconstructs the final output fields from the model's tracers and auxiliary fields. This function performs the following steps:

1. Reconstructs filtered variables: For each variable to be filtered, it sums the contributions from the individual filter coefficients (_C and _S tracers) using the coefficients from filter_params.
2. Reconstructs spatial mapping fields: If spatial mapping is enabled, the function also reconstructs the xi_ fields that represent the filtered position.
3. Builds mean velocities: If compute_mean_velocities is true, the function reconstructs the mean velocity fields using the xi_ fields.
4. Includes original data: The original data is added to the output dictionary for comparison and analysis if config.output_original_data is true.

By default (outputs = :combined in filter_params), each filtered quantity is output as a single field. With outputs = :separate, each (cosine or sine) term of the filter is output separately, named after its tracer with a _scaled suffix, e.g. b_C1_scaled, xi_u_C1_scaled, and u_C1_scaled for the mean velocities. With sine_parity = :odd, the sine terms are odd in time, so their backward-pass contributions are negated. See filter_output_options.

For the offline filter, the outputs of each pass are defined as that pass's contribution to the filtered fields, so that the forward and backward outputs are summed (see sum_forward_backward_contributions!). The backward pass is run with the velocities negated, so its mean velocity outputs are negated (direction = :backward). With a frequency scale, the mean velocities are the velocity of the mean position (see Varying cutoff frequency).

Arguments

• model: An instance of an AbstractModel containing the tracer and auxiliary fields.
• config: An instance of AbstractConfig with the names of the variables, velocity components, and filter parameters.

Keyword arguments

• direction: :forward (default, also used for the online filter) or :backward, the direction of the offline filter pass that the outputs are for.

Returns

A Dict where keys are the names of the output fields (e.g., var_name_Lagrangian_filtered, xi_vel_name, var_name) and values are the corresponding reconstructed Fields.

filtered_vars_advection(config::AbstractConfig; advection = WENO())

The tracer advection schemes for the filtered variables of the online filter, to merge into the model's tracer_advection: advection for the Lagrangian filter, or nothing (no advection) for the Eulerian filter (lagrangian = false). Tracers missing from a NamedTuple tracer_advection are advected with Centered(), so list all the other tracers too, e.g.

tracer_advection = merge((; b = WENO(), T = WENO()), filtered_vars_advection(filter_config))

get_frequency_response(; freq::AbstractArray, filter_params::NamedTuple, direction::Symbol = :both, term = nothing)

Calculates the frequency response $\hat{G}(\omega)$ of the filter, i.e. the Fourier transform of the weight function (see get_weight_function), so that filtering $e^{i\omega t}$ gives $\hat{G}(\omega) e^{i\omega t}$.

Keyword arguments

• freq: A vector of frequencies $\omega$ (in radians per unit time).
• filter_params: A NamedTuple containing the filter coefficients (a, b, c, d) and the number of coefficient pairs (N_coeffs).
• direction: :both (default) for the offline filter, :forward for the online filter (or the forward pass of the offline filter), or :backward for the backward pass of the offline filter.
• term: nothing (default) for the whole filter, or a single term, labelled as in the names of the separate outputs (see create_output_fields), e.g. "C1" (for b_C1_scaled) or "S2".

Returns

• A vector Ghat of the frequency response at each frequency in freq. This is real for the offline filter (direction = :both) if the weight function is even in time (sine_parity = :even, the default), and complex otherwise.

get_weight_function(; t::AbstractArray, tref::Real, filter_params::NamedTuple, direction::Symbol = :both, term = nothing)

Computes the weight function of the filter, which determines how much each point in the timeseries t contributes to the filtered value at a reference time tref. For the offline filter this is

$G(t) = \sum_i e^{-c_i |\tau|} \left[a_i \cos(d_i \tau) + b_i \sin(d_i |\tau|)\right], \quad \tau = t_{ref} - t,$

or with sine terms $b_i \sin(d_i \tau)$ if filter_params.sine_parity = :odd.

Keyword arguments

• t: A collection of time points in the timeseries.
• tref: The reference time at which the filter is being evaluated.
• filter_params: A NamedTuple containing the coefficients (a, b, c, d) and the number of coefficient pairs (N_coeffs).
• direction: :both (default) for the offline filter, :forward to keep only times before tref (as for the online filter, or the forward pass of the offline filter), or :backward to keep only times after tref (the backward pass of the offline filter).
• term: nothing (default) for the whole filter, or a single term, labelled as in the names of the separate outputs (see create_output_fields), e.g. "C1" (for b_C1_scaled) or "S2".

Returns

• A vector of weights G, with the same dimensions as t, representing the value of the weight function at each time point relative to tref.

initialise_filtered_vars_from_data(model, reader::BufferedDataReader, config)

Initialises the model's tracer fields, which represent the components of the filtered variables, to the (scaled) original data at sim_t = 0 (i.e. the physical start/end of the filter interval for forward/backward runs). This improves the "spin-up" of the filter simulation by providing a good starting point. The data is read into the model's velocities and auxiliary fields, as at each time step, and the filtered variables are then initialised from the model with initialise_filtered_vars_from_model.

• For a single-exponential filter (N_coeffs = 0.5), only the _C1 tracer exists and is initialised.
• For a multi-coefficient filter (N_coeffs > 0.5), both the _C and _S fields for each coefficient are initialised.

If maps are being computed (compute_maps or compute_mean_velocities), the maps are initialised from the velocities in the same way.

initialise_filtered_vars_from_model(model::AbstractModel,config::AbstractConfig)

Initializes the model's filtered tracer fields using the actual tracer fields that are assumed to have been already set. This improves the "spin-up" of the filter simulation by providing a good starting point.

The initialization formula depends on the number of filter coefficients (N_coeffs):

• For a single-exponential filter (N_coeffs = 0.5), only the _C1 tracer exists and is initialized.
• For a multi-coefficient filter (N_coeffs > 0.5), both the _C and _S tracers for each coefficient are initialized.

Both the filtered variables and the maps are initialised.

Arguments

• model: The AbstractModel whose tracers are to be initialized.
• config: An instance of AbstractConfig with the filter parameters.

jld2_to_netcdf(jld2_filename::String, nc_filename::String)

Converts a JLD2 output file generated by an Oceananigans simulation into a standard NetCDF file. This function is useful for post-processing and for sharing data with other tools that expect the NetCDF format.

The conversion process involves the following steps:

1. Read JLD2 data: Opens the input JLD2 file and reads the grid, time, and all timeseries variables.
2. Create NetCDF file: Creates a new NetCDF file with a .nc extension.
3. Define dimensions: Defines NetCDF dimensions based on the grid sizes and staggered locations (e.g., x_caa for cell centers, x_faa for cell faces). Only interior points are included, not halo cells.
4. Define grid variables: Writes the grid coordinates, the cell spacings at each location (e.g., dz_aaf, or dx_fca on a LatitudeLongitudeGrid, where it varies with latitude) and metadata (e.g., Lx, Ny) as variables to the NetCDF file.
5. Write timeseries data: Iterates through each variable in the JLD2 file's timeseries, determines its location on the grid, and writes the interior of the data (dropping any halo cells) to a new variable in the NetCDF file.
6. Add metadata: Adds attributes to each variable, including boundary conditions and units, for better documentation.

Arguments

• jld2_filename: A String specifying the path to the input JLD2 file.
• nc_filename: A String specifying the path for the output NetCDF file.

regrid_to_mean_position!(config::AbstractConfig; extra_vars_to_regrid = (), architecture = filter_architecture(config),
                         maxiter = 50)

Regrid the Lagrangian filtered fields to the mean position, and save each as a new variable with the suffix _at_mean in the output file.

Each filtered field $f^*$ is stored at the regular grid points, which label the particles, but belongs at their mean position $\Xi(\boldsymbol{x}) = \boldsymbol{x} + \boldsymbol{\xi}(\boldsymbol{x})$. For each grid point $\boldsymbol{y}$, the label $\boldsymbol{x} = \Xi^{-1}(\boldsymbol{y})$ is found by solving $\boldsymbol{x} + \boldsymbol{\xi}(\boldsymbol{x}) = \boldsymbol{y}$ with Newton's method, with the maps $\boldsymbol{\xi}$ interpolated linearly between cell centres, and the regridded field is $f^*(\boldsymbol{x})$. This runs on the CPU or GPU, on a RectilinearGrid or LatitudeLongitudeGrid. On a LatitudeLongitudeGrid the maps are converted from metres to degrees, which assumes the grid is locally Cartesian over the trajectory excursions (see Regridding to the mean position).

The regridded field is NaN at grid points that are not the mean position of any particle in the fluid: where the label $\boldsymbol{x}$ lies outside the domain, or in an immersed solid (in a solid cell, or below the bottom height). This happens, for example, near a boundary where the particles' mean positions are pulled away from it. It is also NaN at grid points in an immersed solid (including those below the bottom height in a partial cell), and where Newton's method does not converge. Beyond bounded boundaries, the maps and fields are extrapolated linearly, so nothing is assumed about the flow at the boundary. A warning is given if the map folds (so that some grid points are the mean position of more than one label).

Arguments

• config: The OfflineFilterConfig or OnlineFilterConfig used for the filtering.
• extra_vars_to_regrid: Additional filtered variables in the output file to regrid.
• architecture: The architecture to regrid on. Default: the offline filter's architecture, or the architecture of the online model's grid.
• maxiter: The maximum number of Newton iterations. Default: 50.

set_offline_BW2_filter_params(; N::Int=1, freq_c::Real=1)

Calculates the coefficients for a filter that has a frequency response given by a Butterworth filter with order N and cutoff frequency freq_c, squared.

Uses N exponentials and N/2 sets of coefficients (a,b,c,d). N should therefore be even, since exponentials come in pairs to ensure a real-valued filter.

However, the special case N=1 is allowed, which gives a single (real) exponential filter.

Frequency response: Ghat(omega) = 1 / (1 + (omega / freq_c)^(2*N)) Real filter shape: G(t) = sum_{i=1}^{N/2} exp(-c_i*abs(t))*(a_i*cos(d_i * abs(t)) + b_i*sin(d_i * abs(t)))

This function supports two types of filters:

• A single exponential filter when N=1. This is a special case that generates two coefficients instead of 4. The unidirectional filter is a single exponential, and N_coeffs = 0.5. Only a1 and c1 are returned.
• A Butterworth squared filter for N>1. This generates N/2 sets of coefficients (a,b,c,d), representing a filter of order N. The coefficients are computed based on the filter's order and cutoff frequency.

Arguments

• N: The order parameter for the filter. N=1 for a single exponential. For N>1, the filter's order is N. Must be a non-negative even integer.
• freq_c: The cutoff frequency of the filter. Must be a real number.

Returns

• A NamedTuple containing the filter coefficients, N_coeffs (the number of coefficient pairs), and the output options outputs = :combined and sine_parity = :even (the defaults, see filter_output_options).

set_offline_spectrum_filter_params(; alpha::Real, freqs::AbstractVector, normalisation::Symbol = :spectral)

Coefficients for an exponentially-windowed spectral filter, which extracts the signal in a narrow band around each frequency in freqs. The window is w(t) = exp(-alpha*|t|), giving the in-phase (C) and quadrature (S) weight functions at each frequency omega_n:

C(t) = A_n exp(-alpha*|t|) cos(omega_n*t)
S(t) = A_n exp(-alpha*|t|) sin(omega_n*t)

so that a_n = b_n = A_n, c_n = alpha and d_n = omega_n. The S weight functions are odd in t (sine_parity = :odd), and each term is output separately (outputs = :separate).

alpha sets the frequency resolution: the window lasts about 1/alpha, so frequencies closer than about alpha are not resolved.

Two normalisations are available:

• :spectral (default): A_n = sqrt(alpha), so that the window has unit energy, ∫ (A_n w)^2 dt = 1. Then C^2 + S^2 estimates the (two-sided) power spectral density of the signal at omega_n, normalised so that the variance is ∫ S(omega) domega / 2π.
• :unit_gain: A_n = alpha (alpha^2 + 4omega_n^2) / (2alpha^2 + 4omega_n^2), so that C passes a signal at exactly omega_n through unchanged: C is the band-passed signal, and S is its quadrature (a quarter period later). Then sqrt(C^2 + S^2) is its envelope and atan(S, C) its instantaneous phase (omega_n*t plus the phase of the signal). S has gain 4omega_n^2 / (2alpha^2 + 4omega_n^2), so the envelope and phase are accurate to about alpha^2 / (2omega_n^2). At omega_n = 0 this reduces to alpha/2, the single exponential of set_offline_BW2_filter_params(N = 1).

The S terms can be left out of the outputs by setting their coefficients b_n to zero. See set_online_spectrum_filter_params for the online filter.

Arguments

• alpha: Decay rate of the exponential window. Must be positive.
• freqs: Frequencies (radians per unit time) at which to extract the signal.
• normalisation: :spectral (default) or :unit_gain.

Returns

• A NamedTuple of coefficients, N_coeffs, outputs = :separate and sine_parity = :odd.

set_online_BW_filter_params(; N::Int=1, freq_c::Real=1)

Calculates the coefficients for a filter that has a frequency response given by a Butterworth filter with order N and cutoff frequency freq_c. Note that the frequency response is not squared, like in the offline forward-backward filter, and the frequency response is not real-valued, implying a nonlinear phase shift.

Uses N exponentials and N/2 sets of coefficients (a,b,c,d). N should therefore be even, since exponentials come in pairs to ensure a real-valued filter.

However, the special case N=1 is allowed, which gives a single (real) exponential filter.

Frequency response: abs(Ghat(omega)) = 1 / sqrt(1 + (omega / freq_c)^(2*N)) Real filter shape: G(t) = sum_{i=1}^{N/2} exp(-c_i*t)* (a_i*cos(d_i * t) + b_i*sin(d_i * t)) for t>=0, and 0 for t<0.

This function supports two types of filters:

• A single exponential filter when N=1. This is a special case that generates two coefficients instead of 4. The unidirectional filter is a single exponential, and N_coeffs = 0.5. Only a1 and c1 are returned.
• A Butterworth filter for N>1. This generates N/2 sets of coefficients (a,b,c,d), representing a filter of order N. The coefficients are computed based on the filter's order and cutoff frequency.

Arguments

• N: The order parameter for the filter. N=1 for a single exponential. For N>1, the filter's order is N. Must be a non-negative even integer.
• freq_c: The cutoff frequency of the filter. Must be a real number.

Returns

• A NamedTuple containing the filter coefficients, N_coeffs (the number of coefficient pairs), and the output option outputs = :combined (the default, see filter_output_options).

set_online_spectrum_filter_params(; alpha::Real, freqs::AbstractVector, normalisation::Symbol = :spectral)

Coefficients for an exponentially-windowed spectral filter for the online filter: the causal version of set_offline_spectrum_filter_params, which only uses the past. The window is w(t) = exp(-alpha*t) for t > 0 (the time before the present), giving the in-phase (C) and quadrature (S) weight functions at each frequency omega_n:

C(t) = A_n exp(-alpha*t) cos(omega_n*t)
S(t) = A_n exp(-alpha*t) sin(omega_n*t)

so that a_n = b_n = A_n, c_n = alpha and d_n = omega_n, and each term is output separately (outputs = :separate).

alpha sets the frequency resolution: the window lasts about 1/alpha, so frequencies closer than about alpha are not resolved. Because the window only uses the past, changes in the amplitude of the signal show up in the envelope sqrt(C^2 + S^2) after a delay of about 1/alpha (so t - 1/alpha is the time for the envelope or power, rather than the time computed by compute_time_shift!).

Two normalisations are available:

• :spectral (default): A_n = sqrt(2alpha), so that the window has unit energy, ∫ (A_n w)^2 dt = 1. Then, as for the offline filter, C^2 + S^2 estimates the (two-sided) power spectral density of the signal at omega_n, normalised so that the variance is ∫ S(omega) domega / 2π.
• :unit_gain: A_n = alpha sqrt(alpha^2 + 4omega_n^2) / sqrt(alpha^2 + omega_n^2), so that C passes a signal at exactly omega_n with its amplitude preserved, and S gives its quadrature with gain omega_n / sqrt(alpha^2 + omega_n^2). Then sqrt(C^2 + S^2) is its envelope and atan(S, C) its instantaneous phase (omega_n*t plus the phase of the signal). Unlike the offline filter, the one-sided window shifts the phases slightly: C lags the signal and S leads its quadrature, each by about alpha / (2omega_n) radians, so the envelope and phase have a ripple of relative size about alpha / (2omega_n) at twice the frequency (for a ripple below 1%, use alpha ≲ 0.02 omega_n). At omega_n = 0 this reduces to alpha, the single exponential of set_online_BW_filter_params(N = 1).

The S terms can be left out of the outputs by setting their coefficients b_n to zero.

Arguments

• alpha: Decay rate of the exponential window. Must be positive.
• freqs: Frequencies (radians per unit time) at which to extract the signal.
• normalisation: :spectral (default) or :unit_gain.

Returns

• A NamedTuple of coefficients, N_coeffs and outputs = :separate.

sum_forward_backward_contributions!(config::AbstractConfig)

Combines the output from the forward and backward filter simulations into a single output file. This function performs the final step of the offline filter algorithm by summing the contributions from each pass.

The function performs the following steps:

1. Initializes the combined file: A new JLD2 file is created to store the final output.
2. Copies metadata and unfiltered data: The file structure, metadata, and the original, unfiltered data are copied from the forward output file.
3. Sums filtered contributions: For each filtered variable, the data from the backward output file is loaded as a FieldTimeSeries. The data is then interpolated to match the time steps of the forward simulation, and the two datasets are summed and written to the combined output file.

The outputs of each pass are that pass's contribution to the filtered fields (see create_output_fields), so all filtered fields, including the mean velocities, are summed.

Arguments

• config: An instance of AbstractConfig containing the file paths and variable names.
• extra_filtered_var_names::Tuple{Vararg{String}}=(): Optional tuple of additional filtered variable names that have been calculated by the filter and also need to be combined. The outputs of each pass should be defined as its contribution, so that they are summed (e.g. negating quantities that change sign under time reversal in the backward pass, as create_output_fields does for the mean velocities).
• extra_original_data_names::Tuple{Vararg{String}}=(): Optional tuple of additional original names that have been output and should be copied to the combined output file.

zero_closure_for_filtered_vars(config::AbstractConfig)

The diffusivities for the filtered variables of the online filter, all zero, to merge into the κ of the model's closure, e.g.

closure = ScalarDiffusivity(ν = 1e-6, κ = merge((; b = 1e-6), zero_closure_for_filtered_vars(filter_config)))

The filtered variables are then not diffused. Includes the maps if they are solved for.

Arguments

• config: An instance of AbstractConfig with the filter parameters.