Radial finite elements and the per-degree operator
After the spectral reduction, each spherical-harmonic degree \(j\) is an independent 1D radial boundary-value problem. VILMA2 discretizes it with finite elements in radius, assembling the symmetric saddle-point operator that is the Hessian of the energy functional.
Mesh and bases
The mesh spans the whole sphere \(\langle 0, a\rangle\) — the inviscid core is meshed too — with piecewise-constant \(\rho_k,\mu_k\) per element. Each element’s size follows the target for its own depth (default 5 km above 70 km, 10 km to 670 km, 40 km below, after VEGA; namelist mesh_dr, mesh_depth, mesh_dr_scale). Layer interfaces are always nodes, and within a layer the elements share \(\int \mathrm{d}r/h\) equally, so a thick layer is meshed like a stack of thin ones. M3–L70–V01 uses 218 nodes (217 elements), the 11-layer PREM 220 nodes (219 elements).
The table below shows how the solution converges with the mesh. It scales the whole spacing by mesh_dr_scale and compares the elastic Love numbers of M3–L70–V01 with the normal-mode table (see Love numbers) and the displacement of the low-viscosity-zone benchmark (see LVZ). It is produced by make diag_mesh_spacing and bin/diag_mesh_spacing.x love|lvz3d.
mesh_dr_scale |
elements | max error \(h/\ell/k\), \(j \le 48\) [%] | max error \(\ell/k\), \(j \le 256\) [%] | LVZ \(U\) at 200 yr [m] |
|---|---|---|---|---|
| 2 | 110 | 0.05 / 0.40 / 0.76 | 29 / 46 | −1.2079 |
| 1 (default) | 217 | 0.06 / 0.15 / 0.25 | 7.6 / 11.7 | −1.2180 |
| 0.5 | 434 | 0.06 / 0.08 / 0.11 | 1.95 / 2.97 | −1.2186 |
| 0.25 | 867 | 0.06 / 0.07 / 0.08 | 0.54 / 0.79 | — |
The high-degree errors of \(\ell\) and \(k\) fall by about a factor of four each time the spacing is halved. They are therefore a second-order discretization error of the mesh, not a limit of the reference. The vertical number \(h\) and the low degrees are converged at the default spacing.
Two finite-element bases are used (Martinec 2000, §6):
- P1 “tent” functions \(\psi_k(r)\) for the nodal fields: \(U,V,F = \sum_k [\,\cdot\,]_k\,\psi_k\) (continuous, one value per node).
- P0 (piecewise-constant) functions \(\xi_k(r)\) for the pressure: \(\Pi = \sum_k \Pi_k\,\xi_k\) (one value per element).
This P1(U,V,F) / P0(\(\Pi\)) pairing is the inf–sup-stable choice for the incompressible saddle-point problem: the pressure enforces \(\operatorname{div}\mathbf{u}=0\) weakly through its coupling block.
Element bilinear forms
The element matrices are built from the closed-form radial integrals of Martinec’s Appendix C (the \(I_1\ldots I_7\), \(K_1\ldots K_3\) families), implemented and unit-tested in vilma_radial_integrals:
| Block | Equation | Integrals | Couples |
|---|---|---|---|
| Pressure / incompressibility \(\delta E_{\text{press}}\) | eq. 82 | \(K_1,K_2\) | \(\Pi \leftrightarrow (U' + 2U/r - JV/r)\) |
| Shear stiffness \(\delta E_{\text{shear}}\) | eq. 80 | \(I_1,I_3,I_6\) (\(\times\mu_k\)) | \(U,V\) |
| Self-gravity \(\delta E_{\text{grav}}\) | eq. 81 | \(I_2,I_4,I_5,I_7\), \(\tfrac{1}{4\pi G}\) | \(U,V,F\) |
| Rigid-mode penalty \(\delta E_{\text{uniq}}\) | eq. 83 | \(K_3\) | \(U,V\) at \(j=1\) |
The singular \(I_7\) (the \(1/r\) term carrying \(R_k\)) is skipped in the innermost element where \(R_1 = 0\), avoiding a \(0\cdot\infty\); the \(r^2\) weighting handles regularity at the centre so no explicit centre boundary condition is needed.
A symmetric operator — and the bug that hid in it
Laid out node-interleaved as \([U\ V\ F\mid\Pi]\), the per-degree system is
\[ \begin{bmatrix} A & B^{\mathsf T}\\ B & 0\end{bmatrix} \begin{bmatrix} d\\ \Pi\end{bmatrix} = \begin{bmatrix} f\\ 0\end{bmatrix}, \]
with \(A\) the shear+gravity stiffness and \(B\) the incompressibility coupling. The whole operator is symmetric — it is the second variation (Hessian) of the energy functional \(E\), hence self-transpose by construction. In particular the self-gravity \(U\!\leftrightarrow\! F\) coupling is a transpose pair: the potential-gradient body force \(\int\rho_0(\mathrm{d}F/\mathrm{d}r)\,\delta U\,r^2\) and the Poisson source \(\int\rho_0 U(\mathrm{d}\,\delta F/\mathrm{d}r)\,r^2\) discretize to \(I^2_{\beta\alpha}\) and \(I^2_{\alpha\beta}\) respectively.
The discretized self-gravity block written out in Martinec (2000) eq. (81) appears to contain a typo: it gives the \(U\)–\(F\) coupling as \(I^2_{\alpha\beta}\) where the symmetric partner \(I^2_{\beta\alpha}\) is required. Because the full operator is the Hessian of the energy functional \(E\), it must be symmetric; the corresponding term in the continuous form (eq. 65) makes this clear on term-by-term inspection. Using eq. (81) as written breaks the operator symmetry and softens the elastic low-degree Love numbers markedly — \(h_e(2)\) was \(\sim 47\%\) below the benchmark in our own tests, while every analytic limit and the fluid limit still passed. Anyone implementing this block should take eq. (65) as the arbiter, use \(I^2_{\beta\alpha}\) for the \(U\)–\(F\) coupling, and assert \(\lVert A - A^{\mathsf T}\rVert = 0\) on the assembled matrix.
Boundary and regularity conditions
- Surface — the exterior potential match (\(\varphi_1\propto r^{-(j+1)}\) outside) adds a \((j+1)F(a)\) term on the \(F\!-\!F\) diagonal at the surface node; the load enters as the right-hand side \(-a^2\sigma\,g_0(a)\) on \(U(a)\) and \(-a^2\sigma\) on \(F(a)\) (eq. 84).
- Centre \(r=0\) — no explicit BC (meshed through, \(r^2\) weighting + \(R_1=0\)).
- Core–mantle boundary — none imposed; free-slip emerges from \(\mu=0\) in the meshed core.
- Density-jump interfaces — natural conditions of the weak form.
Degree 1: the geocenter, kept sparse
Degree 1 carries the rigid-translation null space, removed by \(E_{\text{uniq}}\) (eq. 83). That penalty is a rank-1 term whose coefficient is \(\sim 10^{16}\times\) the band entries — i.e. it is already a hard constraint \(w^{\mathsf T}d = 0\) (zero volume-integrated displacement, a centre-of-figure-like gauge) in all but name. Imposing it as a dense penalty would fill the band; instead it borders the band with one KKT row/column:
\[ \begin{bmatrix} A_{\text{band}} & w\\ w^{\mathsf T} & 0\end{bmatrix} \begin{bmatrix} d\\ \lambda\end{bmatrix} = \begin{bmatrix} f\\ 0\end{bmatrix}. \]
The solve runs at \(\text{ndof}+1\) internally but keeps the physical \(\text{ndof}\) interface, so the time stepper needs no degree-1 special case. The same KKT system also yields the null mode itself (zero physical right-hand side, border value 1). Because a rigid translation carries no strain, the reference frame is a post-solve choice: adding a multiple of the null mode that makes the exterior degree-1 potential vanish, \(F_1(a) = 0\), puts the state in the CM frame (Blewitt 2003) (deg1_frame = "cm", the default). The geocenter is thus the physical degree-1 signal, not a hard-coded \(\{h_1,l_1,k_1\}\). See Sea-level equation.
The toroidal operator
When the toroidal field is carried, \(W\) has its own per-degree operator rather than a fifth interleaved field: \(W\) couples to nothing else on the left-hand side, so the spheroidal band is untouched. It is the \(W\) part of the shear energy (eq. 80), assembled from the same \(I_1,I_3,I_6\) integrals. It is tridiagonal and symmetric, with fluid-core nodes pinned to \(W = 0\) and, at degree 1, one KKT border per solid shell removing the rigid rotation. See Toroidal flow.
vilma_radial_fe assembles build_dense_operator (the reference dense operator, with the dense \(E_{\text{uniq}}\) penalty for cross-checking) and the banded operator used in production; vilma_radial_integrals holds the Appendix-C element integrals. Love numbers are extracted by loading_love — see Love numbers. Validated by test_assembly, test_love.