`FastIsostasy.ForwardPlan`

— Type`ForwardPlan`

Allias for in-place precomputed plans from FFTW or CUFFT. Used to compute forward FFT.

`FastIsostasy.InversePlan`

— Type`InversePlan`

Allias for in-place precomputed plans from FFTW or CUFFT. Used to compute inverse FFT.

`FastIsostasy.ComputationDomain`

— Type```
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)`

`FastIsostasy.CurrentState`

— Type`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.

`FastIsostasy.DenseOutputs`

— Type`DenseOutputs()`

Return a struct containing the fields of viscous displacement, viscous displacement rate, elastic displacement, geoid displacement, sea level and the computation time resulting from solving a `FastIsoProblem`

.

`FastIsostasy.FastIsoProblem`

— Type```
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}`

.

`FastIsostasy.FastIsoTools`

— Type`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 geoid displacement, plans for FFTs, interpolators of the load and the viscosity over time and preallocated arrays.

`FastIsostasy.InversionConfig`

— Type`InversionConfig`

Struct containing configuration parameters for a [`InversionProblem`

].

Need to choose regularization factor α ∈ (0,1], When you have enough observation data α=1: no regularization

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

`FastIsostasy.InversionData`

— Type`InversionData`

Struct containing data (either observational or output of a golden standard model) for a [`InversionProblem`

].

`FastIsostasy.InversionProblem`

— Type`InversionProblem`

Struct containing variables and configs for the inversion of Solid-Earth parameter fields. For now, only viscosity can be inverted but future versions will support lithosphere rigidity. For now, the unscented Kalman inversion is the only method available but ensemble Kalman inversion will be available in future.

`FastIsostasy.LayeredEarth`

— Type`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`

.

`FastIsostasy.OceanSurfaceChange`

— Type`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)`

`FastIsostasy.PhysicalConstants`

— Type`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).

`FastIsostasy.ReferenceEarthModel`

— Type`ReferenceEarthModel`

Return a struct with vectors containing the:

- radius (distance from Earth center),
- depth (distance from Earth surface),
- density,
- P-wave velocities,
- S-wave velocities,

which are typically used to characterize the properties of a spherically symmetrical solid Earth.

`FastIsostasy.ReferenceState`

— Type`ReferenceState`

Return a struct containing the reference state.

`FastIsostasy.SolverOptions`

— Type`Options`

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

.

`FastIsostasy.analytic_solution`

— Method`analytic_solution(r, t, c, p, H0, R0, domains; n_quad_support)`

Return the analytic solution of the bedrock displacement resulting from a cylindrical ice load with radius `R0`

and height `H0`

on a flat Earth represented by an elastic plate overlaying a viscous half space. Parameters are provided in `c, p`

. The points at which the solution is computed are specified by the distance `r`

to the center of the domain. The time at which the solution is computed is specified by `t`

.

`FastIsostasy.apply_bc!`

— Method`apply_bc!(u::M, bcm::M, nbc::T)`

A generic function to update `u`

such that the boundary conditions are respected on average, which is the only way to impose them for Fourier collocation.

`FastIsostasy.build_greenintegrand`

— Method```
build_greenintegrand(distance::Vector{T},
greenintegrand_coeffs::Vector{T}) where {T<:AbstractFloat}
```

Compute the integrands of the Green's function resulting from a load at a given `distance`

and based on provided `greenintegrand_coeffs`

. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

`FastIsostasy.corner_bc!`

— Method`corner_bc!(u, Nx, Ny, offset)`

Shift `u`

s.t. its mean value at the corner points is `offset`

.

`FastIsostasy.dist2angulardist`

— Method`dist2angulardist(r::Real)`

Convert Euclidean to angular distance along great circle.

`FastIsostasy.dudt_isostasy!`

— Method`dudt_isostasy!()`

Update the displacement rate `dudt`

of the viscous response.

`FastIsostasy.edge_bc!`

— Method`edge_bc!(u, Nx, Ny)`

Shift `u`

s.t. its mean value at on the edge is 0. Same as in bueler-fast-2007.

`FastIsostasy.extract_inversion`

— Method`extract_inversion()`

Extract results of parameter inversion from the `priors`

and `ukiobj`

that resulted from `solve!(paraminv::InversionProblem)`

.

`FastIsostasy.gauss_distr`

— Method```
gauss_distr(X::KernelMatrix{T}, Y::KernelMatrix{T},
mu::Vector{<:Real}, sigma::Matrix{<:Real})
```

Compute `Z = f(X,Y)`

with `f`

a Gaussian function parametrized by mean `mu`

and covariance `sigma`

.

`FastIsostasy.get_differential_fourier`

— Method`get_differential_fourier(W, N2)`

Compute the matrices representing the differential operators in the fourier space.

`FastIsostasy.get_effective_viscosity`

— Method```
get_effective_viscosity(
layer_viscosities::Vector{KernelMatrix{T}},
layers_thickness::Vector{T},
Omega::ComputationDomain{T, M},
) where {T<:AbstractFloat}
```

Compute equivalent viscosity for multilayer model by recursively applying the formula for a halfspace and a channel from Lingle and Clark (1975).

`FastIsostasy.get_elasticgreen`

— Method`get_elasticgreen(Omega, quad_support, quad_coeffs)`

Integrate load response over field by using 2D quadrature with specified support points and associated coefficients.

`FastIsostasy.get_geoidgreen`

— Method`get_geoidgreen(fip::FastIsoProblem)`

Return the Green's function used to compute the geoid anomaly as in ^{[Coulon2021]}.

`FastIsostasy.get_greenintegrand_coeffs`

— Method`get_greenintegrand_coeffs(T)`

Return the load response coefficients with type `T`

. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

`FastIsostasy.get_loadgreen`

— Method```
get_loadgreen(r::T, rm::Vector{T}, greenintegrand_coeffs::Vector{T},
interp_greenintegrand_::Interpolations.Extrapolation) where {T<:AbstractFloat}
```

Compute the integrands of the Green's function resulting from a load at a given `distance`

and based on provided `greenintegrand_coeffs`

. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

`FastIsostasy.get_normalized_lin_transform`

— Method`get_normalized_lin_transform(x1, x2)`

Return parameters of linear function mapping `x1, x2`

onto `-1, 1`

.

`FastIsostasy.get_quad_coeffs`

— Method`get_quad_coeffs(T, n)`

Return support points and associated coefficients with specified Type for Gauss-Legendre quadrature.

`FastIsostasy.get_r`

— Method`get_r(x::T, y::T) where {T<:Real}`

Get euclidean distance of point (x, y) to origin.

`FastIsostasy.get_rigidity`

— Method`get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}`

Compute rigidity `D`

based on thickness `t`

, Young modulus `E`

and Poisson ration `nu`

.

`FastIsostasy.init`

— Method`init(fip)`

Initialize an `ode::CoupledODEs`

, aimed to be used in `step!`

.

`FastIsostasy.kernelpromote`

— Method`kernelpromote(X, arraykernel)`

Promote X to the kernel (`Array`

or `CuArray`

) specified by `arraykernel`

.

`FastIsostasy.latlon2stereo`

— Method`latlon2stereo(lat, lon, lat0, lon0)`

Compute stereographic projection (x,y) for a given latitude `lat`

longitude `lon`

, reference latitude `lat0`

and reference longitude `lon0`

. Optionally one can provide `lat::KernelMatrix`

and `lon::KernelMatrix`

if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).

`FastIsostasy.load_dataset`

— Method`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:

- "ICE6G
*D": ice loading history from ICE6G*D. - "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"

`FastIsostasy.load_prem`

— Method`load_prem()`

Load Preliminary Reference Earth Model (PREM) from Dzewonski and Anderson (1981).

`FastIsostasy.lon360tolon180`

— Method`lon360tolon180(lon, X)`

Convert longitude and field from `lon=0:360`

to `lon=-180:180`

.

`FastIsostasy.m_per_sec2mm_per_yr`

— Method`m_per_sec2mm_per_yr(dudt::Real)`

Convert displacement rate `dudt`

from $ m \, s^{-1} $ to $ mm \, \mathrm{yr}^{-1} $.

`FastIsostasy.matrify`

— Method`matrify(x, Nx, Ny)`

Generate a vector of constant matrices from a vector of constants.

`FastIsostasy.meshgrid`

— Method`meshgrid(x, y)`

Return a 2D meshgrid spanned by `x, y`

.

`FastIsostasy.normalized_lin_transform`

— Method`normalized_lin_transform(y, m, p)`

Apply normalized linear transformation with slope `m`

and bias `p`

on `y`

.

`FastIsostasy.quadrature1D`

— Method`quadrature1D(f, n, x1, x2)`

Compute 1D Gauss-Legendre quadrature of `f`

between `x1`

and `x2`

based on `n`

support points.

`FastIsostasy.quadrature2D`

— Method`quadrature2D(f, x, w, x1, x2, y1, y2)`

Return the integration of `f`

over [`x1, x2`

] x [`y1, y2`

] with `x, w`

some pre-computed support points and coefficients of the Gauss-Legendre quadrature.

`FastIsostasy.reinit_structs_cpu`

— Method`reinit_structs_cpu(Omega, p)`

Reinitialize `Omega::ComputationDomain`

and `p::LayeredEarth`

on the CPU, mostly for post-processing purposes.

`FastIsostasy.samesize_conv_indices`

— Method`samesize_conv_indices(N, M)`

Get the start and end indices required for a `samesize_conv`

`FastIsostasy.savefip`

— Method`savefip(filename, fip; T = Float32)`

Save the output of `fip::FastIsoProblem`

as NetCDF file under `filename`

.

`FastIsostasy.scalefactor`

— Method```
scalefactor(lat::T, lon::T, lat0::T, lon0::T) where {T<:Real}
scalefactor(lat::M, lon::M, lat0::T, lon0::T) where {T<:Real, M<:KernelMatrix{T}}
```

Compute scaling factor of stereographic projection for a given `(lat, lon)`

and reference `(lat0, lon0)`

. Angles must be provided in radians. Reference: ^{[Snyder1987]}, p. 157, eq. (21-4).

`FastIsostasy.seconds2years`

— Method`seconds2years(t::Real)`

Convert input time `t`

from seconds to years.

`FastIsostasy.simple_euler!`

— Method`simple_euler!()`

Update the state `u`

by performing an explicit Euler integration of its derivative `dudt`

over a time step `dt`

.

`FastIsostasy.solve!`

— Method`solve!(fip)`

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

.

`FastIsostasy.solve!`

— Method`solve!(paraminv::InversionProblem)`

Return `priors`

and `ukiobj`

that allow to extract the results of the parameter inversion as initialized in `paraminv`

.

`FastIsostasy.step!`

— Method`step!(fip)`

Step `fip::FastIsoProblem`

over `tspan`

and based on `ode::CoupledODEs`

, typically obtained by `init`

.

`FastIsostasy.stereo2latlon`

— Method`stereo2latlon(x, y, lat0, lon0)`

Compute the inverse stereographic projection `(lat, lon)`

based on Cartesian coordinates `(x,y)`

and for a given reference latitude `lat0`

and reference longitude `lon0`

. Optionally one can provide `x::KernelMatrix`

and `y::KernelMatrix`

if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!

Convert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).

`FastIsostasy.surfacechange_residual!`

— Method```
surfacechange_residual!(osc, z, delta_V)
surfacechange_residual!(Vresidual, z, z_k, A_itp, delta_V)
```

Update the residual of the piecewise linear approximation, used to solve the sea-level/ocean-surface nonlinearity.

`FastIsostasy.three_layer_scaling`

— Method```
three_layer_scaling(Omega::ComputationDomain, kappa::T, visc_ratio::T,
channel_thickness::T)
```

Return the viscosity scaling for a three-layer model and based on a the wave number `kappa`

, the `visc_ratio`

and the `channel_thickness`

. Reference: Bueler et al. 2007, below equation 15.

`FastIsostasy.update_V_af!`

— Method`update_V_af!(fip::FastIsoProblem)`

Update the volume contribution from ice above floatation as in goelzer-brief-2020 (eq. 13). Note: we do not use (eq. 1) as it is only a special case of (eq. 13) that does not allow a correct representation of external sea-level forcings.

`FastIsostasy.update_V_den!`

— Method`update_V_den!(fip::FastIsoProblem)`

Update the volume contribution associated with the density difference between meltwater and sea water, as in goelzer-brief-2020 (eq. 10).

`FastIsostasy.update_V_pov!`

— Method`update_V_pov!(fip::FastIsoProblem)`

Update the volume contribution to the ocean (from isostatic adjustement in ocean regions), which corresponds to the "potential ocean volume" in goelzer-brief-2020 (eq. 14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcings.

`FastIsostasy.update_bsl!`

— Method`update_bsl!(fip::FastIsoProblem)`

Update the sea-level contribution of melting above floatation and density correction. Note that this differs from goelzer-brief-2020 (eq. 12) because the ocean surface is not assumed to be constant. Furthermore, the contribution to ocean volume from the bedrock uplift is not included here since the volume displaced on site is arguably blanaced by the depression of the peripherial forebulge.

`FastIsostasy.update_diagnostics!`

— Method`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.

`FastIsostasy.update_elasticresponse!`

— Method`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.

`FastIsostasy.update_geoid!`

— Method`update_geoid!(fip::FastIsoProblem)`

Update the geoid by convoluting the Green's function with the load anom.

`FastIsostasy.update_loadcolumns!`

— Method`update_loadcolumns!(fip::FastIsoProblem, u::KernelMatrix{T}, H_ice::KernelMatrix{T})`

Update the load columns of a `::CurrentState`

.

`FastIsostasy.update_seasurfaceheight!`

— Method`update_seasurfaceheight!(fip::FastIsoProblem)`

Update the sea-level by adding the various contributions as in coulon-contrasting-2021. Here, the constant term is used to impose a zero geoid perturbation in the far field rather thank for mass conservation and is embedded in convolution operation.

`FastIsostasy.where_significant`

— Method`where_significant()`

Find points of parameter field that can be inverted. We here assume that

`FastIsostasy.write_out!`

— Method`write_out!(fip::FastIsoProblem)`

Write results in output vectors if the load is updated internally. If the load is updated externally, the user is responsible for writing results.

`FastIsostasy.years2seconds`

— Method`years2seconds(t::Real)`

Convert input time `t`

from years to seconds.

`FastIsostasy.Data`

— Module`Module Data`

The module Data is created in the module where `@init_parallel_stencil`

is called from. It provides the following types:

`Data.Number`

The type of numbers used by @zeros, @ones, @rand and @fill and in all array types of module `Data`

(selected with argument `numbertype`

of `@init_parallel_stencil`

).

`Data.Array{ndims}`

Expands to `Data.Array{numbertype, ndims}`

, where `numbertype`

is the datatype selected with `@init_parallel_stencil`

and the datatype `Data.Array`

is chosen to be compatible with the package for parallelization selected with `@init_parallel_stencil`

(Array for Threads, CUDA.CuArray or CUDA.CuDeviceArray for CUDA and AMDGPU.ROCArray or AMDGPU.ROCDeviceArray for AMDGPU; `@parallel`

and `@parallel_indices`

convert CUDA.CuArray and AMDGPU.ROCArray automatically to CUDA.CuDeviceArray and AMDGPU.ROCDeviceArray in kernels when required).

`Data.CellArray{ndims}`

Expands to `Data.CellArray{numbertype, ndims}`

, where `numbertype`

is the datatype selected with `@init_parallel_stencil`

and the datatype `Data.CellArray`

is chosen to be compatible with the package for parallelization selected with `@init_parallel_stencil`

(CPUCellArray for Threads, CuCellArray or CuDeviceCellArray for CUDA and ROCCellArray or ROCDeviceCellArray for AMDGPU; `@parallel`

and `@parallel_indices`

convert CellArray automatically to DeviceCellArray when required).

`Data.Cell{S}`

Expands to `Union{StaticArrays.SArray{S, numbertype}, StaticArrays.FieldArray{S, numbertype}}`

, where `numbertype`

is the datatype selected with `@init_parallel_stencil`

.

`Data.DeviceArray{ndims}`

Expands to `Data.DeviceArray{numbertype, ndims}`

, where `numbertype`

is the datatype selected with `@init_parallel_stencil`

and the datatype `Data.DeviceArray`

is chosen to be compatible with the package for parallelization selected with `@init_parallel_stencil`

(Array for Threads, CUDA.CuDeviceArray for CUDA AMDGPU.ROCDeviceArray for AMDGPU).

This datatype is not intended for explicit manual usage. `@parallel`

and `@parallel_indices`

convert CUDA.CuArray and AMDGPU.ROCArray automatically to CUDA.CuDeviceArray and AMDGPU.ROCDeviceArray in kernels when required.

`Data.DeviceCellArray{ndims}`

Expands to `Data.DeviceCellArray{numbertype, ndims}`

, where `numbertype`

is the datatype selected with `@init_parallel_stencil`

and the datatype `Data.DeviceCellArray`

is chosen to be compatible with the package for parallelization selected with `@init_parallel_stencil`

(CPUCellArray for Threads, CuDeviceCellArray for CUDA and ROCDeviceCellArray for AMDGPU).

This datatype is not intended for explicit manual usage. `@parallel`

and `@parallel_indices`

convert CUDA.CuArray and AMDGPU.ROCArray automatically to CUDA.CuDeviceArray and AMDGPU.ROCDeviceArray in kernels when required.