Externalised public API

plot_transect(sim::Simulation, vars;
    ii = sim.domain.mx:sim.domain.nx,
    jj = sim.domain.my,
    analytic_cylinder_solution = false)

Plot the results of a FastIsostasy simulation along a transect. The transect is defined by the indices ii and jj, which specify the range of grid points to plot. The vars argument specifies which variables to plot (e.g., :u, :h, etc.). If analytic_cylinder_solution is set to true, the function will also plot the analytical solution for a cylindrical load.

plot_load(domain, H_ice)

Plot the ice load as a function of position. The domain argument specifies the spatial domain of the simulation, and H_ice is a matrix representing the ice thickness at each grid point.

plot_earth(sim::Simulation)
plot_earth(domain, solidearth)

Plot the solid-Earth parameters used in the simulation.

plot_out_at_time(sim, vars, t)

Plot the output of a FastIsostasy simulation at a specific time t. The vars argument specifies which variables to plot.

plot_out_over_time(sim, var, t_vec, copts)

Plot the output of a FastIsostasy simulation for a specific variable var over a range of times specified in t_vec. The copts argument can be used to customize the plot (e.g., colors, line styles, etc.).

plot_computation_time(sim)

Plot the simulation year as a function of computation time. This can be useful for assessing the performance of the simulation.

AbstractInversion

Supertype of the inversion problems, IceLoadInversion and ParameterInversion. Both share the same machinery — decode θ into a Simulation template, run the forward model to every observation time, and score the result with loss — and differ only in intent and defaults.

IceLoadInversion(sim, encoding, observations; regularizations = (),
    diffmode = TangentMode(), lossmodel = DefaultLoss())

Infer an unknown ice load (and, jointly, solid-Earth parameters) from surface observations. encoding maps the parameter vector θ onto the ice-thickness snapshots and any other model inputs via reconstruct! — e.g. Test1Encoding, which writes time-varying Vialov domes plus a viscosity field.

Arguments and keywords are shared with ParameterInversion:

• sim: the Simulation template. It is reset to its initial condition on every loss evaluation, so one problem can be evaluated repeatedly.
• encoding: an AbstractEncoding, or nothing for full-field control (see ParameterInversion).
• observations: a collection of Observations. All must currently carry the same observable type — a mixed collection is abstractly typed and trips Enzyme.
• regularizations: a collection of AbstractRegularizations (e.g. SurfaceSmoothnessReg) added to the misfit.
• diffmode: TangentMode (forward, low-dimensional θ) or AdjointMode (reverse, cost independent of length(θ)).
• lossmodel: an AbstractLoss defining the misfit norm.

See also loss, gradient!, solve!.

ParameterInversion(sim, encoding, observations; regularizations = (),
    diffmode = TangentMode(), lossmodel = DefaultLoss())

Calibrate solid-Earth parameters (mantle viscosity, densities, …) against observed deformation, with the ice load held known and fixed. Same fields and keywords as IceLoadInversion; the two differ only in intent.

The encoding decides the parameterization:

• an AbstractEncoding such as Test2Encoding reduces the unknowns to a handful of numbers (a background viscosity plus Gaussian anomalies and densities), which suits TangentMode;
• nothing selects full-field control: θ is the flattened log10 effective-viscosity field (column-major, so reshape(θ, size(domain.X)) is the map) and the physical field is 10^θ. There is then one unknown per grid cell, which is only affordable with AdjointMode, whose gradient cost is independent of length(θ).

See also loss, gradient!, solve!.

loss(prob, θ) -> scalar

Decode θ into prob.sim, run the forward model to every observation time, and return misfit(prob.lossmodel, preds, prob.observations) plus the regularization penalties. This is the AD-free objective the diff engines differentiate.

gradient!(g, prob, θ)

In-place gradient of loss(prob, θ) w.r.t. θ. Dispatches on prob.diffmode (TangentMode → forward, AdjointMode → checkpointed reverse), implemented by the Enzyme extension. The 3-arg dispatcher lives in core (so it always resolves); the mode-specific 4-arg method falls back to an informative error here and is overridden by the extension once loaded.

loss_and_gradient!(g, prob, θ) -> loss_value

In-place gradient of loss(prob, θ) w.r.t. θ, written into g, returning the loss value alongside it. Under TangentMode the primal falls out of the same forward pass that computes each directional tangent (Enzyme.ForwardWithPrimal), so this is one fewer forward run per optimizer iteration than calling loss and gradient! separately. Implemented by the Enzyme extension; see gradient! for the dispatcher/fallback split.

solve!(prob, optimizer; kwargs...)

Run the inversion with optimizer. Implemented by the Optim extension.

AbstractDiffMode

Supertype of the differentiation modes, TangentMode (forward) and AdjointMode (reverse). A diff mode carries policy — which Enzyme mode to use, the tangent batch size / checkpoint schedule, and whether a low-dimensional encoding is required — without referencing Enzyme, so the core package stays AD-free; the extensions map it onto the actual engine.

Rule of thumb: forward mode costs one forward run per parameter, so it suits encoded (low-dimensional) θ; reverse mode costs one checkpointed sweep regardless of length(θ), so it is what makes full-field inversion feasible.

TangentMode(; batch = 8)

Forward-mode (tangent) differentiation. Cost scales with the number of parameters, so it requires a low number of parameter to estimate (which can be achieved via a low-dimensional encoding). batch is the number of tangent directions propagated together (chunking).

AdjointMode(; checkpoint_every = 10)

Reverse-mode (adjoint) differentiation with periodic checkpointing. Cost is roughly independent of the parameter dimension, so full 2-D fields can be inverted without an encoding. checkpoint_every is the number of save/output intervals between stored state snapshots.

AbstractObservable

Supertype of the observable tags — lightweight, field-free types naming what is measured, which an Observation pairs with where and when.

Available tags: VerticalUpliftObservable (u + ue), VerticalUpliftRateObservable (dudt) and RelativeSeaLevelObservable ((z_ss − z_ss_ref) − (u + ue)).

Total vertical bedrock displacement (viscous + elastic), u + ue.

Vertical bedrock uplift rate, dudt (viscous rate; the elastic contribution is piecewise-constant between sparse diagnostic updates and is neglected here).

Relative sea level: sea-surface change minus bedrock uplift, (z_ss − z_ss_ref) − (u + ue).

Observation(tag, points, times, data; σ = 1)

A set of measurements of tag::AbstractObservable at grid points (Vector{CartesianIndex{2}}) and times. data is ordered points-fastest, then times: [p1@t1, p2@t1, …, p1@t2, …], i.e. length npoints * ntimes. σ is a scalar or a per-entry vector of the same length.

SimulatedObservable(tag, points, times, sim)

A virtual station recording tag::AbstractObservable at grid points (Vector{CartesianIndex{2}}) and times during a forward run. data is flat, points-fastest then times (Observation's convention); k is the cursor index into times of the next pending recording. Construct after sim exists (its linear indices are backend-promoted onto sim's field array family), then attach with attach_simobs!(sim, tag, points, times).

attach_simobs!(sim, tag, points, times) -> SimulatedObservable

Build a SimulatedObservable from sim's current field layout and append it to sim.simobs; run!/step! then record it automatically at each of times.

AbstractEncoding{T}

Supertype of the encodings, which map a low-dimensional parameter vector θ onto the high-dimensional model inputs (effective viscosity, densities, ice-thickness snapshots) through reconstruct!. Encoding the unknowns is what makes forward-mode (TangentMode) inversion affordable.

Concrete encodings: Test1Encoding (joint Vialov ice + bimodal viscosity), Test2Encoding (4-Gaussian viscosity + densities). EOFEncoding / AutoEncoding / VariationalAutoEncoding are stubs for trained decoders.

Passing encoding = nothing to an inversion instead selects full-field control (one unknown per grid cell); see ParameterInversion.

An encoding must be Enzyme-legal: indexed reads of θ and in-place broadcasts into sim fields, no closures over untracked state.

nparams(encoding) -> Int

The number of parameters the encoding expects, i.e. the required length(θ).

reconstruct!(sim, params, reduction)

Reconstruct the parameter values from reduction and update sim accordingly.

Test1Encoding(knot_times, radii, visc_amps; scale = ones)

Encoding for Test 1. Fixed config: knot_times (K ice-interpolation times), radii (the 3 fixed Vialov radii Lᵢ), visc_amps (the 2 fixed log10 anomaly amplitudes, e.g. (-1, +1) decades).

θ layout (length 3K + 13): [Hc₁(1:K), Hc₂(1:K), Hc₃(1:K), x₁,y₁, x₂,y₂, x₃,y₃, log10η_bg, μ₁ₓ,μ₁ᵧ,σ₁, μ₂ₓ,μ₂ᵧ,σ₂].

scale (length 3K + 13, default all-ones) rescales each parameter before it is used: the physical value is θ[i] * scale[i]. This lets the optimization variable θ be dimensionless and O(1) even though the physical parameters span many orders of magnitude (metres of thickness, metres of position, decades of viscosity) — essential for L-BFGS conditioning. With the default ones, θ is the physical parameter vector directly. The mapping is a plain broadcast, so it stays Enzyme-legal.

Test2Encoding(; scale = ones(19))

Encoding for Test 2. θ layout (length 19): [log10η_bg, (μₓ,μᵧ,σ,amp)×4, ρ_uppermantle, ρ_litho]. Ice thickness is not touched (assumed known). Concrete Float64 (the inversion-run precision decision, §1 Precision row).

scale (length 19, default all-ones) rescales each parameter before use: the physical value is θ[i] * scale[i]. As with Test1Encoding, this lets the optimization variable θ be dimensionless and O(1) even though the physical parameters span decades of viscosity, metres of anomaly position/width and thousands of kg/m³ of density — essential for L-BFGS conditioning. The mapping is a plain broadcast, so it stays Enzyme-legal.

EOFEncoding

Empirical-orthogonal-function (linear) decoder: field = mean + Modes * θ. Stub — the decoder application must be implemented Enzyme-legally.

AutoEncoding

Autoencoder decoder. Stub.

VariationalAutoEncoding

Variational-autoencoder decoder. Stub.

AbstractLoss

Supertype of the loss models, which customize only the misfit term of loss — the norm applied to the (prediction, data, σ) triples — via misfit(lossmodel, preds, observations). The default is DefaultLoss.

A loss model is stored as a field of the inversion problem (like diffmode), so any arrays it carries (e.g. a noise-covariance factor) are picked up automatically by the AD engines. Implement a new one by subtyping and adding a misfit method.

DefaultLoss()

The weighted-L2 misfit ½·Σₒ‖(predₒ − dataₒ)/σₒ‖², summed over observations.

misfit(lossmodel, preds, observations) -> scalar

The data-misfit term of loss: the norm of the residual between the predictions preds (one vector per observation, already gathered by the forward run) and the observed data. Extension point of AbstractLoss — subtype it and add a method here to plug in a custom norm.

AbstractRegularization

Supertype of the regularizers and soft bounds added to the data misfit in loss. Each implements penalty(reg, sim, θ), evaluated with sim already holding the decoded parameters, so field-based penalties can read the model state directly.

Concrete types: TikhonovReg (the one configurable regularizer, with the L2Reg and SurfaceSmoothnessReg shorthands) and DecodedBounds (soft two-sided hinge bounds).

TikhonovReg(target::AbstractRegTarget, order::AbstractRegOrder = Order0();
    λ = 1, weights = nothing)

One configurable Tikhonov-style regularizer; see the module docstring above for target/order. weights, if given, must match the target's shape/length and only applies under Order0().

L2Reg(λ)

Thin constructor over TikhonovReg: magnitude penalty on the full raw parameter vector, λ · ‖θ‖².

SurfaceSmoothnessReg(λ)

Thin constructor over TikhonovReg: gradient penalty on the implied ice surface, λ · Σₖ ‖∇(Hₖ + b_ref)‖².

DecodedBounds(quantity, lo, hi, λ)

Soft two-sided hinge penalty keeping a decoded quantity within [lo, hi]: λ · Σ [relu(lo − d)² + relu(d − hi)²], where d = decoded(quantity, sim) is either a field or a scalar. quantity is a BoundedQuantity (Log10Viscosity, UpperMantleDensity, LithoDensity).

AbstractRegTarget

Supertype of the TikhonovReg targets, i.e. what is penalised: ThetaTarget (the raw parameter vector or a subset), FieldTarget (a decoded field or scalar) or SurfaceTarget (the implied ice surface at every snapshot).

ThetaTarget(idx = Colon())

Regularization target: the raw parameter vector θ (encoded space), or a subset of it selected by idx (anything that indexes a vector — Colon(), UnitRange, Vector{Int}). Only meaningful with Order0() — θ has no spatial structure to take a gradient of.

FieldTarget(quantity::BoundedQuantity)

Regularization target: a decoded field or scalar, decoded(quantity, sim). quantity is one of the BoundedQuantity selectors shared with DecodedBounds (Log10Viscosity, UpperMantleDensity, LithoDensity).

SurfaceTarget()

Regularization target: the implied ice surface s = Hₖ + b_ref at every ice snapshot Hₖ (IceLoadInversion), with a fixed reference bed b_ref = sim.ref.z_b. Order1() gives the original SurfaceSmoothnessReg penalty.

AbstractRegOrder

Supertype of the TikhonovReg orders, i.e. how the target is penalised: Order0 penalises magnitude, Order1 penalises the spatial gradient (smoothness).

Magnitude penalty Σ wᵢ xᵢ² (optionally weighted).

Gradient penalty Σ ‖∇x‖² via the FD stencils (field/surface targets only).

BoundedQuantity

Supertype of the decoded-quantity selectors, which name a physical quantity to penalise without closing over state (Enzyme-friendly). Used by FieldTarget and DecodedBounds.

Available: Log10Viscosity, UpperMantleDensity, LithoDensity.

Log10Viscosity()

Selects log10.(sim.solidearth.effective_viscosity) — a field.

UpperMantleDensity()

Selects sim.solidearth.rho_uppermantle — a scalar.

LithoDensity()

Selects sim.solidearth.rho_litho — a scalar.

StateSnapshot(sim)

Allocate a buffer holding a full copy of sim's mutable state (its CurrentState, its barystatic-sea-level object, and the clock). Reuse it across many snapshot!/restore! calls.

snapshot!(buf::StateSnapshot, sim) -> buf

Copy sim's current mutable state into the preallocated buf.

restore!(sim, buf::StateSnapshot) -> sim

Write the state stored in buf back into sim, so a forward run resumed from here reproduces the trajectory captured at snapshot! time.