Time stepping and the field driver

The viscoelastic evolution is an explicit march in time built directly on the elastic operator. This page covers how a step is taken, how the memory state is carried, and how the field-level driver feeds the sea-level equation without corrupting that state.

One step

The \(\omega=1\) Maxwell scheme (see Rheology) keeps the left-hand operator fixed and updates only the right-hand side. A step at degree \(j\) is:

Form the dissipative forcing \(-\int\boldsymbol{\tau}^{V,i-1}:\delta\boldsymbol{\varepsilon}\) from the previous memory stress (a two-point radial Gauss quadrature of the spectral double-dot over the four spheroidal tensor components).
Add the surface-load right-hand side for the current load.
Back-substitute through the once-factored operator (no re-factorization).
Update the memory stress \(\boldsymbol{\tau}^{V,i} = (1-M)\boldsymbol{\tau}^{V,i-1} - 2\mu M\,\boldsymbol{\varepsilon}^{i}\) from the new strain.

Elastic layers (\(M\to 0\)) freeze; fluid layers (\(\mu=0\)) carry no memory. The 1D validation test_relax shows a held load relaxing monotonically from the elastic Love number to the fluid limit \(-(2j+1)/3\), with \(t_{\text{relax}}\propto\eta\) and convergence under step refinement.

This explicit update is conditionally stable — the step is capped by the shortest Maxwell time, so a short dynamical \(\Delta t\) is sub-cycled within the coupling interval (the fe_timestep controller sizes it; scheme = "trap" selects the 2nd-order adaptive variant).

The viscoelastic field driver

The SLE needs a field load whose every coefficient \((\ell,m)\) has its own independent load history — the ice and ocean evolve differently everywhere. The ve_response driver provides this while remaining safe to call from inside an iterative solver, using the affine structure of the response at fixed time:

\[ u_{\ell m} \;=\; g_u(\ell)\,\sigma_{\ell m} \;+\; \text{drift}_{\ell m}, \]

where \(g_u(\ell)\) is the instantaneous (elastic) gain and \(\text{drift}_{\ell m}\) comes from the frozen past-relaxation memory. The step is bracketed:

begin_step freezes the drift — one memory-forcing solve per coefficient, done once at the start of the step;
apply() is then a cheap affine evaluation the SLE fixed point may call repeatedly;
commit_step advances the Maxwell memory once, with the converged load.

Each complex \((\ell,m)\) history is stored as two real histories (real/imaginary), since the operator and \(M=\mu\Delta t/\eta\) are real. The per-element Maxwell kernel is shared verbatim with the 1D stepper — test_ve_response confirms the field driver’s first step equals elastic_response exactly, and that a held degree-2 field reproduces the 1D ve_degree history to \(\sim 10^{-13}\).

Memory storage and restart

The per-\((\ell,m)\) memory and drift are stored in degree-grouped slot order (an \(\ell m \leftrightarrow k\) map), so the inner loop is contiguous and per-degree, allowing operator reuse without scattering the (large) memory arrays. The only prognostic state carried across a restart is the viscoelastic memory — the Maxwell stress \(\boldsymbol{\tau}\) — plus the integrator state and the model time. Reference fields (equilibrium bedrock, reference ice/topo) are static and validated on read. test_restart confirms a checkpoint-restart reproduces an uninterrupted run bit-for-bit. A lower-resolution restart is interpolated up to the model grid (degree-grouped spectral block copy), so a run can be spun up cheaply at coarse resolution and continued at fine resolution.

Coupling cadence

For host coupling (CLIMBER-X), the model sub-cycles its short dynamical step \(\Delta t\) within a longer coupling interval: ice thickness is held across the interval, the SLE is solved each \(\Delta t\), and the bedrock \(z_{\text{bed}} = z_{\text{bed,eq}} - \text{rsl}\) is returned at the coupling cadence. The reference state and its datum are chosen by i_eq (e.g. a present-day RTopo reference, or an equilibrium configuration); a non-equilibrium start (a held LGM load) is brought into balance by a spin-up (solid_earth_spinup) that relaxes the held load to bed-stationary equilibrium (rate below equil_rate_tol) — capped by equil_time_max — with an optional cheap 1-D pre-equilibration (pre_spinup_1d) before switching the transient to the 3-D path.

Where this lives in the code

fe_viscoelastic (1D stepper + shared kernel), fe_response (ve_response field driver, begin_step/commit_step), fe_io (netCDF restart), and fe_coupling (host init/update contract). The cost of begin_step is the main performance driver and is addressed in Solver.