FastIsostasy.ForwardPlanType
ForwardPlan

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

FastIsostasy.InversePlanType
InversePlan

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

FastIsostasy.ComputationDomainType
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.CurrentStateType
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.DenseOutputsType
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.FastIsoProblemType
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.FastIsoToolsType
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.InversionConfigType
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.InversionDataType
InversionData

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

FastIsostasy.InversionProblemType
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.LayeredEarthType
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.OceanSurfaceChangeType
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.PhysicalConstantsType
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.ReferenceEarthModelType
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.analytic_solutionMethod
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_greenintegrandMethod
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.extract_inversionMethod
extract_inversion()

Extract results of parameter inversion from the priors and ukiobj that resulted from solve!(paraminv::InversionProblem).

FastIsostasy.gauss_distrMethod
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_effective_viscosityMethod
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_elasticgreenMethod
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_greenintegrand_coeffsMethod
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_loadgreenMethod
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_quad_coeffsMethod
get_quad_coeffs(T, n)

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

FastIsostasy.get_rMethod
get_r(x::T, y::T) where {T<:Real}

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

FastIsostasy.get_rigidityMethod
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.kernelpromoteMethod
kernelpromote(X, arraykernel)

Promote X to the kernel (Array or CuArray) specified by arraykernel.

FastIsostasy.latlon2stereoMethod
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_datasetMethod
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"
FastIsostasy.load_premMethod
load_prem()

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

FastIsostasy.matrifyMethod
matrify(x, Nx, Ny)

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

FastIsostasy.quadrature1DMethod
quadrature1D(f, n, x1, x2)

Compute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.

FastIsostasy.quadrature2DMethod
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_cpuMethod
reinit_structs_cpu(Omega, p)

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

FastIsostasy.savefipMethod
savefip(filename, fip; T = Float32)

Save the output of fip::FastIsoProblem as NetCDF file under filename.

FastIsostasy.scalefactorMethod
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.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.stereo2latlonMethod
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_scalingMethod
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_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_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.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.DataModule
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.


Advanced
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).

Warning

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).

Warning

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.