API reference

ComputationDomain
ComputationDomain(W, n)
ComputationDomain(Wx, Wy, Nx, Ny)

Return a struct containing all information related to geometry of the domain and potentially used parallelism. To initialize one with 2*W and 2^n grid cells:

Omega = ComputationDomain(W, n)

If a rectangular domain is needed, run:

Omega = ComputationDomain(Wx, Wy, Nx, Ny)

PhysicalConstants

Return a struct containing important physical constants. Comes with default values that can however be changed by the user, for instance by running:

c = PhysicalConstants(rho_ice = 0.93)   # (kg/m^3)

All constants are given in SI units (kilogram, meter, second).

LayeredEarth(Omega; layer_boundaries, layer_viscosities)

Return a struct containing all information related to the lateral variability of solid-Earth parameters. To initialize with values other than default, run:

Omega = ComputationDomain(3000e3, 7)
lb = [100e3, 300e3]
lv = [1e19, 1e21]
p = LayeredEarth(Omega, layer_boundaries = lb, layer_viscosities = lv)

which initializes a lithosphere of thickness $T_1 = 100 \mathrm{km}$, a viscous channel between $T_1$and $T_2 = 200 \mathrm{km}$and a viscous halfspace starting at $T_2$. This represents a homogenous case. For heterogeneous ones, simply make lb::Vector{Matrix}, lv::Vector{Matrix} such that the vector elements represent the lateral variability of each layer on the grid of Omega::ComputationDomain.

ReferenceState

Return a struct containing the reference state.

CurrentState

Return a mutable struct containing the geostate which will be updated over the simulation. The geostate contains all the states of the [FastIsoProblem] to be solved.

FastIsoTools(Omega, c, p)

Return a struct containing pre-computed tools to perform forward-stepping of the model. This includes the Green's functions for the computation of the lithosphere and the SSH perturbation, plans for FFTs, interpolators of the load and the viscosity over time and preallocated arrays.

Options

Return a struct containing the options relative to solving a FastIsoProblem.

OceanSurfaceChange(; z0 = 0.0)

Return a mutable struct OceanSurfaceChange containing:

• z_k: the GMSL at current time step k.
• A_k: the ocean surface at current time step k.
• z: a vector of GMSL values used as knots for interpolation.
• A: a vector of ocean surface values used as knots for interpolation.
• A_itp: an interpolator of ocean surface over depth. Bias-free for present-day.
• A_pd: the present-day ocean surface.
• res: residual of the nonlinear equation solved numerically.

An osc::OceanSurfaceChange can be used as function to update osc.z_k and osc.A_k based on osc.A_itp and an input delta_V by running:

osc(delta_V)

FastIsoProblem(Omega, c, p, t_out)
FastIsoProblem(Omega, c, p, t_out, Hice)
FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice)

Return a struct containing all the other structs needed for the forward integration of the model over Omega::ComputationDomain with parameters c::PhysicalConstants and p::LayeredEarth. The outputs are stored at t_out::Vector{<:AbstractFloat}.

solve!(fip)

Solve the isostatic adjustment problem defined in fip::FastIsoProblem.

init(fip)

Initialize an ode::CoupledODEs, aimed to be used in step!.

step!(fip)

Step fip::FastIsoProblem over tspan and based on ode::CoupledODEs, typically obtained by init.

update_diagnostics!(dudt, u, fip, t)

Update all the diagnotisc variables, i.e. all fields of fip.now apart from the displacement, which requires an integrator.

lv_elva!(dudt, u, fip, t)

Update the displacement rate dudt of the viscous response according to LV-ELVA.

update_elasticresponse!(fip::FastIsoProblem)

Update the elastic response by convoluting the Green's function with the load anom. To use coefficients differing from ^{[Farrell1972]}, see FastIsoTools.

InversionConfig

Struct containing configuration parameters for a [InversionProblem].

Fields

• method::Any: Inversion method to use.
• paramspriors::NamedTuple: Prior information about the parameters to invert.
• N_iter::Int: Number of iterations for the inversion.
• α_reg::Real: Regularization factor. When you have enough observation data α=1 (no regularization)
• update_freq::Int: Update frequency for the inversion.

1 : approximate posterior cov matrix with an uninformative prior. 0 : weighted average between posterior cov matrix with an uninformative prior and prior.

• n_samples::Int: Number of samples for the inversion.
• scale_obscov::Real: Scaling factor for the observational covariance matrix.

InversionData

Struct containing the inversion data.

Fields

• t::Vector{T}: Time vector.
• nt::Int: Number of time steps.
• X::Vector{M}: Ground truth input (forcing).
• Y::Vector{M}: Ground truth response.
• mask::BitMatrix: Region of interest.
• countmask::Int: count(mask) = number of cells used for inversion.

InversionProblem

Struct containing variables and configs for the inversion of Solid-Earth parameter fields. InversionProblem needs to be initialized using inversion_problem. For now, the unscented Kalman inversion is the only method available.

Fields

• fip::FastIsoProblem: FastIsoProblem object.
• config::InversionConfig: Configuration for the inversion.
• data::InversionData: Data for the inversion.
• reduction::R: Parameter reduction method.
• priors::PD: Prior distribution.
• ukiobj::EKP: Unscented Kalman inversion object.
• error::V: Error vector.
• out::Vector{V}: Output vector.
• G_ens::M: Ensemble of the covariance matrix.

inversion_problem(fip, config, data, reduction, priors; save_stride_iter::Int = 1)

Generate an inversion problem for the given fip::FastIsoProblem object.

ParameterReduction

Abstract type for parameter reduction methods. Any subtype must implement the reconstruct!(fip, theta) method, which assigns the reconstructed parameter values to fip::FastIsoProblem.

print_inversion_evolution(paraminv, n, ϕ_n, reduction)

Print the inversion evolution.

extract_inversion(paraminv, n)

Extract the inversion results to compare them with the ground truth.

reconstruct!(fip, params, reduction)

Reconstruct the parameter values from reduction and update fip accordingly.

extract_output(fip, reduction, data)

Extract the output of the forward run for the inversion.

load_dataset(name) → (dims), field, interpolator

Return the dims::Tuple{Vararg{Vector}}, the field<:Array and the interpolator corresponding to a data set defined by a unique name::String. For instance:

(lon180, lat, t), Hice, Hice_itp = load_dataset("ICE6G_D")

Following options are available for parameter fields:

• "ICE6GD": ice loading history from ICE6GD.
• "Wiens2022": viscosity field from (Wiens et al. 2022)
• "Lithothickness_Pan2022": lithospheric thickness field from (Pan et al. 2022)
• "Viscosity_Pan2022": viscosity field from (Pan et al. 2022)

Following options are available for model results:

• "Spada2011"
• "LatychevGaussian"
• "LatychevICE6G"

reinit_structs_cpu(Omega, p)

Reinitialize Omega::ComputationDomain and p::LayeredEarth on the CPU, mostly for post-processing purposes.

write_out!(now, out, k)

Write results in output vectors.