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 vilma_timestep controller sizes the sub-steps according to the memory scheme (see Rheology):
"fe"(default): the 1st-order update has no embedded error signal, so \(\Delta t\) is set a priori by the stability ceiling \(M = \mu\Delta t/\eta \le\)cfl(default 1; it must lie in (0, 2), the update’s stability range), taken over the fastest Maxwell rate present. The interval is split into equal sub-steps of at most that length.dt_maxlowers the sub-step further.dt_minraises it, also above the stability ceiling, so it should stay at its default of 0 (none) for this scheme. A guard rolls a sub-step back and halves it if the memory norm becomes non-finite. It is a safety net and normally never fires."trap": each candidate step is taken once at \(\Delta t\) and once as two \(\Delta t/2\) halves. Their difference estimates the local error of the memory, which is compared withatol\(+\)rtol\(\times\) the memory norm (defaults \(10^{-3}\) and \(10^{-4}\)). A rejected step is rolled back and retried smaller. An accepted one keeps the fine result. The next step size follows from the error estimate with the factorsafety(0.9). Its change is limited to betweenshrink_min(0.2) andgrow_max(5), and its value to the rangedt_min–dt_max. A step atdt_minis accepted even above tolerance. The first trial step isdt_init(default 0: the whole interval). The step size then carries over from one interval to the next.
The dt_* parameters are given in years. The coupling interval itself is not a parameter: the host, or the forcing time axis of the standalone driver, passes it with each update. The ice load is interpolated linearly between the interval endpoints. scheme is "fe" or "trap" (see Rheology).
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_stepfreezes 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;- the memory advance then moves the Maxwell memory forward once, with the converged load.
The memory advance has two entry points. commit_step takes a given load and is used by the response-level tests. The sea-level solve instead uses prepare_endpoint, advance_endpoint and finalize_step, which let the load and the memory converge together. With scheme = "fe" both do one explicit update. With scheme = "trap" the update needs the end-of-step strain, so the step endpoint is iterated: the drift is re-solved against the trial \(\boldsymbol{\tau}^{V,i}\) until the surface drift stops changing (relative change below \(10^{-6}\)). In the sea-level solve the water load is re-converged in each of these passes, up to sle_max_mem_iter passes (default 20).
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.
A restart file is written without rounding. It carries the following state:
| State | Variables | Carried when |
|---|---|---|
| Maxwell memory stress, four spheroidal channels or six with the toroidal field | tau_* |
viscoelastic response |
| start-of-step load of the trapezoidal scheme, and its flag | sigma_n_*, sigma_primed |
viscoelastic response |
| model time | time |
always |
| next trial step of the step controller | dt_try_s |
always |
| polar motion and the memory of both degree-2 channels | rot_m_*, rot_load_mem, rot_tidal_mem |
rotation on |
| host-resolved terms of the last update (the start of the next ramp) | w_prev, I_g_prev, C_prev, F_prev |
host-resolved terms |
| surface fields of the last update: ice, relative sea level and its two parts, the bed rate, bed, ocean function, barystatic scalars | h_ice, rsl, dz_bed, dz_ss, dz_bed_rate, z_bed, C_ocean, bsl, … |
always |
The memory stress is by far the largest part. The surface fields follow from the memory and the load, but they are restored as well: the ice is the start of the next load ramp, and the relative sea level seeds the next sea-level solve. The elastic and the null response carry no memory, so their restart holds the other rows only. The full variable table is input/vilma-variables.md.
The reference fields are static: the equilibrium bed and ice and, with host-resolved terms, the host reference. They are written once and compared with the initialized model on read. A mismatch stops the run. A file with host-resolved terms cannot be read into a run without them, nor the reverse.
A restart also records the Earth its memory belongs to: checksums of the radial mesh (mesh_checksum), of the radial viscosity (visc_radial_checksum) and of the 3-D field (visc_3d_checksum). A restart from another mesh is refused. A restart from another viscosity is refused too, unless restart_other_visc = .true. (&ctl): its memory then seeds the run (e.g. a 1-D spin-up for a 3-D run) rather than continuing it. Other parameters (the time scheme, the sea-level settings) are not checked.
A run resumes at the restart’s time. The driver starts the forcing window at that time, and a time_init that selects another start stops the run, so a restart cannot replay the forcing by mistake. restart_time_reset = .true. (&ctl) instead takes the restart’s state as the state at the start of the window, whatever time the file carries.
test_restart asserts that a restarted run reproduces an uninterrupted one bit for bit, with rotation on, for a radially symmetric Earth and for a 3-D field with toroidal channels. test_restart_3d asserts the same after the 1-D pre-spin-up, and test_restart_host covers host-resolved terms.
A lower-resolution restart can be read into a higher-resolution run. The memory of the shared low degrees is copied (a contiguous block in the degree-grouped order) and the higher degrees start at zero. The surface fields, the trapezoidal load and the host-resolved terms live on the file’s grid and are not transferred. The caller re-derives the surface fields with a zero-length update. A run can thus be spun up at coarse resolution and continued at fine resolution, but that continuation is not bit for bit. A higher-resolution file cannot be read into a lower-resolution run.
Coupling cadence
For host coupling (CLIMBER-X), the model sub-cycles its short dynamical step \(\Delta t\) within a longer coupling interval. The ice thickness ramps linearly from the previous update’s value to the new one 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 is the relaxed state. The memory is zero there, and relative sea level is measured against it. A host passes the reference at initialization. The standalone driver chooses it with i_eq (&ctl):
i_eq |
Reference bed and ice |
|---|---|
| 0 | the start slice of the forcing |
| 1 (default) | a present-day reference (z_bed_ref_file, h_ice_ref_file), so that rsl is close to zero at the present day |
| 2 | an equilibrium state read from z_bed_eq_file and h_ice_eq_file |
| 3 | the present-day reference bed plus an rsl field from rsl_restart_file |
Spin-up
A run that starts away from the reference (a held LGM load, say) starts with zero memory, as if the load had just been applied. The optional spin-up (solid_earth_spinup) removes this. It holds the start-slice ice and relaxes the model toward isostatic equilibrium, with the reference kept as the datum. The transient then starts with viscous memory. The spin-up has two phases, both off by default:
- Full-model phase (
equil_time_max\(> 0\), in years). The held load is relaxed in passes of growing length: 200 yr first, then 1.6 times longer each pass, at most 60 passes. After each pass the mean bed velocity is the global mean of the absolute bed change over the pass, divided by the pass length. The phase ends when this falls belowequil_rate_tol(default \(10^{-3}\) m/yr).equil_time_maxcaps the total relaxation time. If the cap is reached first, the run continues and prints a warning with the remaining rate. - 1-D pre-phase (
pre_spinup_1d = .true.). The same relaxation first runs on a laterally uniform field, the lateral mean of \(\log_{10}\eta\) of the 3-D field at each radius. It uses the same convergence test and has no time cap. The 3-D field is then restored and the memory is kept. This phase is cheap, because no element is on the 3-D path. It is independent ofequil_time_max.
The rate of each pass is measured over exactly that pass’s memory advance. The criterion is a global mean. It bounds the mean bed velocity, not the velocity at every point, and it is not a measure of how far the relaxation has gone. The default equil_rate_tol stops the spin-up when the bed moves slower than \(10^{-3}\) m/yr. In a test with a held polar cap relaxed over 500 kyr, this is reached at about 83 % of the full relaxation. Full relaxation over such spans is not affordable in realistic cases, so the spin-up gives a practical reference, not a fully relaxed one. For a reproducible reference, fix the length of the spin-up with equil_time_max and set equil_rate_tol = 0, as the reference of the VILMA1 versus VILMA2 deglaciation does. Elements with very high viscosity do not equilibrate within any practical cap, which is what the lid rule addresses. test_spinup asserts that the spin-up relaxes the held load well beyond a single cold step, and that the state after it is near-stationary. The spin-up is not available with solver = "v1".
A zero-length update (dt = 0) advances neither time nor memory. It solves the sea-level equation under the ice it is handed, against the current memory, and refreshes the surface fields, so relative sea level, the ocean function and h_ice describe the same ice. It is idempotent: a second call with the same ice changes nothing. The spin-up begins with one, and a run without spin-up (a restarted run included) begins with one, to seed the start-slice ice. A negative dt stops the run.
vilma_viscoelastic (1D stepper + shared kernel), vilma_response (ve_response field driver, begin_step/commit_step), vilma_io (netCDF restart), and vilma_coupling (host init/update contract, spin-up), vilma_timestep (the sub-step controller) and vilma_sle (the load–memory iteration). The cost of the drift solve in begin_step is addressed in Solver.