Numerical Accuracy¶
Each solution method in gFlex has a distinct accuracy profile.
Analytical solutions (sas and sas_ng)¶
The superposition-of-analytical-solutions methods are exact for a plate
of constant elastic thickness under the no_outside_loads boundary condition
(zero deflection at infinity). In one dimension the Green’s function is the
exponential sinusoid of Eq. (3) in Wickert (2016); in two dimensions it is the
zeroth-order Kelvin function \(\mathrm{kei}\), evaluated via
scipy.special.kei(). Floating-point errors in these special-function
evaluations are at machine precision and negligible in practice.
The only meaningful source of error is a domain that is too small: if the
load produces non-trivial deflection at the domain boundary, the
no_outside_loads assumption is violated and the solution is physically
wrong — not because of numerics, but because the boundary condition does not
match the situation. Use gflex.flexural_wavelengths() to estimate the
flexural parameter \(\alpha\) and ensure the domain extends several
\(\alpha\) beyond the loaded region. Unlike the finite difference solver,
there is no grid-spacing error; output sampling density affects only the
display, not the solution.
FFT spectral solver¶
The FFT solver is spectrally accurate for uniform elastic thickness: for a
smooth load field the solution error decays faster than any finite power of the
grid spacing, limited in practice only by floating-point arithmetic. For
periodic boundary conditions (periodic on all sides) the solution is exact
to machine precision. For all other boundary conditions the load is
zero-padded by \(4\alpha\) on each side before the transform, which
approximates the no_outside_loads condition; the padding introduces a small
error near the domain edges that decreases as the pad width increases relative
to the load’s flexural footprint. For typical geoscience applications the
\(4\alpha\) default padding is more than sufficient.
Because the FFT assumes a single, spatially uniform flexural rigidity \(D\), it cannot represent variable-\(T_e\) problems; those require the finite difference solver.
Finite difference solver (fd)¶
Comparison of numerical (FD) and analytical (SAS) solutions in one dimension
(a) and two dimensions (c), and their differences (b, d), for a
100 km central line load / circular load. The 1–2 m offset in (b) is due
primarily to the no_outside_loads BC of the analytical solution versus the
zero_displacement_zero_slope BC of the FD solution; the cross-shaped residual in
(d) reflects boundary effects along the longer diagonal boundaries.
Reproduced from Wickert (2016), Fig. 3;
CC BY 3.0.¶
2-D finite-difference solver: second-order convergence¶
The 2-D finite-difference solver (method = 'fd') is second-order accurate in space: halving the grid spacing
\(\Delta x\) reduces the numerical error by a factor of approximately
four (\(\mathcal{O}(\Delta x^2)\)).
This is verified two ways in gFlex:
Method of Manufactured Solutions (MMS) — a formal unit test (
tests.test_2D_FD.test_2d_fd_convergence_order()) that runs with every CI build. A cosine deflection field is chosen as the exact solution; the corresponding load is derived analytically and fed to the solver. Three grid refinements (N = 30, 60, 120 cells per side) on a periodic domain confirm a convergence rate > 1.8 at all refinement levels.Kelvin-function benchmark — a grid-convergence study comparing the FD solver to the analytical Kelvin-function solution for an infinite plate under a point load. Results are summarised in the table below.
Kelvin-function benchmark setup¶
Young’s modulus \(E\) |
65 GPa |
Elastic thickness \(T_e\) |
10 km |
Poisson’s ratio \(\nu\) |
0.25 |
Mantle density \(\rho_m\) |
3300 kg m⁻³ |
Infill density \(\rho_\text{fill}\) |
0 kg m⁻³ (air) |
\(g\) |
9.8 m s⁻² |
Flexural parameter \(\alpha\) |
≈ 21 km |
Domain |
600 km × 600 km, point load at centre |
Boundary conditions |
|
Grid spacings tested |
\(\Delta x\) = 10 000, 5 000, 2 500, 1 250 m |
Convergence orders (coarse → fine)¶
The table gives the local convergence order between successive grid refinements at several distances from the load, expressed as multiples of the flexural parameter α.
\(r/\alpha\) |
\(\Delta x\) 10→5 km |
\(\Delta x\) 5→2.5 km |
\(\Delta x\) 2.5→1.25 km |
|---|---|---|---|
0.97 |
2.20 |
2.07 |
2.01 |
1.46 |
2.68 |
2.18 |
2.04 |
1.95 |
1.07 |
1.88 |
1.97 |
2.92 |
−0.07 |
1.77 |
1.95 |
At the finest spacing (\(\Delta x = 1250\) m, \(\Delta x/\alpha \approx 0.06\)) relative errors are 0.02–0.1 % everywhere except very close to the singularity of the point load.
The sub-second-order rates at \(r/\alpha\) = 1.95 and 2.92 on the coarsest grids reflect pre-asymptotic effects near the forebulge zero-crossing, where the solution changes sign and the absolute error temporarily dominates the relative error. All rates converge to ≈ 2 at finer resolution, confirming asymptotic second-order behaviour.
Practical guidance¶
For most geoscience applications, \(\Delta x / \alpha \leq 0.1\) (ten
cells per flexural parameter) is sufficient to keep errors below 1 %. At
\(\Delta x / \alpha \approx 0.5\) (two cells per flexural parameter)
errors can be several percent, particularly near the load centre and the
forebulge. Use gflex.flexural_wavelengths() to estimate \(\alpha\)
before choosing a grid spacing.
1-D zero_displacement_zero_slope boundary condition: MMS verification¶
The zero_displacement_zero_slope (clamped) boundary condition enforces
\(w = 0\) and \(dw/dx = 0\) at the plate edge. The manufactured
solution
satisfies both conditions at \(\xi = 0\) and \(\xi = 1\) exactly. Its fourth derivative is constant, so the manufactured load
includes a spatially varying elastic-foundation term — making this a nontrivial test of the full governing equation, not just the biharmonic operator alone. The error metric is the \(L^\infty\) relative error:
Physical parameters: \(T_e = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km.
Results¶
\(n_x\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
1.84 × 10⁻³ |
51 |
12.0 |
4.55 × 10⁻⁴ |
101 |
6.0 |
1.14 × 10⁻⁴ |
201 |
3.0 |
2.85 × 10⁻⁵ |
401 |
1.5 |
7.12 × 10⁻⁶ |
801 |
0.75 |
1.78 × 10⁻⁶ |
Convergence slope (finest three points): \(\mathcal{O}(\Delta x^{2.00})\).
The test function
tests.test_bc_mms.TestClampedBC1D.test_convergence_order()
runs with every CI build.
1-D and 2-D zero_moment_zero_shear boundary condition: MMS verification¶
The manufactured solution
satisfies all four free-end boundary conditions (\(w''=w'''=0\)) at both ends exactly. Its fourth derivative is
giving the manufactured load
The 2-D extension uses the separable solution \(w = -W_0\,f(\xi)\,f(\eta)\) with \(f(t) = t^4(1-t)^4\), which satisfies the free-end condition on all four sides.
Physical parameters: \(T_e = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km, \(W_0 = 25600\) m (giving \(|w_\mathrm{exact}|_\mathrm{max} = 100\) m).
Results (1-D)¶
\(n_x\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
1.72 × 10⁻² |
51 |
12.0 |
4.40 × 10⁻³ |
101 |
6.0 |
1.11 × 10⁻³ |
201 |
3.0 |
2.77 × 10⁻⁴ |
401 |
1.5 |
6.93 × 10⁻⁵ |
801 |
0.75 |
1.73 × 10⁻⁵ |
Convergence slope (finest three points): \(\mathcal{O}(\Delta x^{2.00})\).
Results (2-D)¶
\(n_x = n_y\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
1.76 × 10⁻² |
51 |
12.0 |
4.53 × 10⁻³ |
101 |
6.0 |
1.15 × 10⁻³ |
201 |
3.0 |
2.85 × 10⁻⁴ |
Convergence slope (finest two points): \(\mathcal{O}(\Delta x^{2.01})\).
The test functions
tests.test_bc_mms.TestFreeEndBC1D.test_convergence_order() and
tests.test_bc_mms.TestFreeEndBC2D.test_convergence_order() run with
every CI build.
2-D zero_displacement_zero_slope boundary condition: MMS verification¶
The exact solution
satisfies all four clamped boundary conditions (\(w = 0\), \(\partial w/\partial n = 0\)) exactly. Because \(g''''(t) = 24\) (constant), the manufactured load
includes a spatially varying elastic-foundation term. The error metric is the \(L^\infty\) relative error.
Physical parameters: \(T_e = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km, \(W_0 = 1600\) m (max \(|w_\mathrm{exact}| = 6.25\) m).
Results¶
\(n_x = n_y\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
1.61 × 10⁻³ |
51 |
12.0 |
4.31 × 10⁻⁴ |
101 |
6.0 |
1.10 × 10⁻⁴ |
201 |
3.0 |
2.75 × 10⁻⁵ |
Convergence slope (finest two points): \(\mathcal{O}(\Delta x^{1.99})\).
The test function
tests.test_bc_mms.TestClampedBC2D.test_convergence_order()
runs with every CI build.
zero_slope_zero_shear (mirror) boundary condition: MMS verification¶
The zero_slope_zero_shear (mirror) boundary condition enforces
\(dw/dx = 0\) and \(d^3w/dx^3 = 0\) at the plate edge. A
cosine satisfies both conditions identically, giving the manufactured
solution
and the manufactured load
The 2-D extension uses the separable solution \(w = -W_0\cos(\pi\xi)\cos(\pi\eta)\), for which the load is
Physical parameters: \(T_e = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km, \(W_0 = 1\) m.
Results (1-D)¶
\(n_x\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
9.50 × 10⁻⁶ |
51 |
12.0 |
2.38 × 10⁻⁶ |
101 |
6.0 |
5.94 × 10⁻⁷ |
201 |
3.0 |
1.49 × 10⁻⁷ |
401 |
1.5 |
3.81 × 10⁻⁸ |
Convergence slopes (coarse → fine): \(\mathcal{O}(\Delta x^{2.06})\), \(\mathcal{O}(\Delta x^{2.03})\), \(\mathcal{O}(\Delta x^{2.01})\), \(\mathcal{O}(\Delta x^{1.97})\).
Results (2-D)¶
\(n_x = n_y\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
3.76 × 10⁻⁵ |
51 |
12.0 |
9.40 × 10⁻⁶ |
101 |
6.0 |
2.35 × 10⁻⁶ |
201 |
3.0 |
5.88 × 10⁻⁷ |
Convergence slopes (coarse → fine): \(\mathcal{O}(\Delta x^{2.06})\), \(\mathcal{O}(\Delta x^{2.03})\), \(\mathcal{O}(\Delta x^{2.01})\).
The test functions tests.test_bc_mms.TestMirrorBC1D.test_convergence_order()
and tests.test_bc_mms.TestMirrorBC2D.test_convergence_order() run with
every CI build and assert slopes > 1.8 at all refinement levels.
1-D and 2-D zero_displacement_zero_moment boundary condition: MMS verification¶
The zero_displacement_zero_moment (pinned / simply-supported) boundary
condition enforces \(w = 0\) (zero deflection) and \(d^2w/dx^2 = 0\)
(zero bending moment) at the plate edge. A sine satisfies both conditions
exactly, giving the manufactured solution
and the manufactured load
The 2-D extension uses the separable solution \(w = -W_0\sin(\pi\xi)\sin(\pi\eta)\), for which the load is
Physical parameters: \(T_e = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.81\) m s⁻², \(L = 600\) km, \(W_0 = 1600\) m.
Results (1-D)¶
\(n_x\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
9.49 × 10⁻⁶ |
51 |
12.0 |
2.37 × 10⁻⁶ |
101 |
6.0 |
5.94 × 10⁻⁷ |
201 |
3.0 |
1.49 × 10⁻⁷ |
401 |
1.5 |
3.81 × 10⁻⁸ |
Convergence slopes (coarse → fine): \(\mathcal{O}(\Delta x^{2.00})\), \(\mathcal{O}(\Delta x^{2.00})\), \(\mathcal{O}(\Delta x^{2.00})\), \(\mathcal{O}(\Delta x^{1.96})\).
Results (2-D)¶
\(n_x = n_y\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
3.75 × 10⁻⁵ |
51 |
12.0 |
9.39 × 10⁻⁶ |
101 |
6.0 |
2.35 × 10⁻⁶ |
201 |
3.0 |
5.87 × 10⁻⁷ |
Convergence slopes (coarse → fine): \(\mathcal{O}(\Delta x^{2.00})\), \(\mathcal{O}(\Delta x^{2.00})\), \(\mathcal{O}(\Delta x^{2.00})\).
The test functions tests.test_bc_mms.TestPinnedBC1D.test_convergence_order()
and tests.test_bc_mms.TestPinnedBC2D.test_convergence_order() run with
every CI build.
Variable elastic thickness: 1-D MMS verification¶
The finite-difference solver supports spatially variable elastic thickness \(T_e(x)\). To verify convergence, a linearly varying rigidity
is constructed by setting
so that \(D(\xi) = E T_e^3/[12(1-\nu^2)]\) ranges from \(D_0\) at \(\xi = 0\) to \(1.5\,D_0\) at \(\xi = 1\). The exact solution is the clamped shape function
Expanding the governing equation \(\partial^2/\partial x^2\!\left[D\,\partial^2 w/\partial x^2\right] + \Delta\rho g\,w = -q_s\) for linear \(D\) (so \(D'' = 0\)) and noting that \(w'''' = 24/L^4\) (constant) gives the manufactured load
Physical parameters: \(T_{e,\mathrm{ref}} = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km, \(W_0 = 1\) m, clamped BCs on both ends.
Results¶
\(n_x\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
2.07 × 10⁻³ |
51 |
12.0 |
5.10 × 10⁻⁴ |
101 |
6.0 |
1.27 × 10⁻⁴ |
201 |
3.0 |
3.17 × 10⁻⁵ |
401 |
1.5 |
7.92 × 10⁻⁶ |
801 |
0.75 |
1.98 × 10⁻⁶ |
Convergence slopes (finest four points): \(\mathcal{O}(\Delta x^{2.01})\), \(\mathcal{O}(\Delta x^{2.01})\), \(\mathcal{O}(\Delta x^{2.00})\).
The test function tests.test_bc_mms.TestVariableTe1D.test_convergence_order()
runs with every CI build.
Variable elastic thickness: 2-D MMS verification¶
The 2-D variable-\(T_e\) test uses a rigidity that varies in both directions:
constructed by setting
The exact solution is the same separable clamped shape as the constant-\(T_e\) 2-D case:
Because \(D\) varies in both \(x\) and \(y\), the governing equation (van Wees & Cloetingh 1994) contains cross-derivative terms that are absent in the 1-D case. For bilinear \(D\) — where \(\partial^2 D/\partial x^2 = \partial^2 D/\partial y^2 = 0\) — the equation reduces to
where the last three terms arise from spatial variation of \(D\) and are absent when \(D\) is uniform. The final term uses coefficient \(2(1-\nu)\) rather than \(2\), reflecting the full van Wees & Cloetingh form (see Theory and Numerics). A non-zero cross-derivative term exercises the off-diagonal entries of the variable-\(T_e\) finite-difference stencil and provides a stronger verification than the 1-D strip test alone. The manufactured load \(q_s\) is derived analytically from the equation above with \(w = w_\mathrm{exact}\).
Physical parameters: \(T_{e,\mathrm{ref}} = 30\) km, \(E = 65\) GPa, \(\nu = 0.25\), \(\rho_m = 3300\) kg m⁻³, \(\rho_\mathrm{fill} = 0\), \(g = 9.8\) m s⁻², \(L = 600\) km, clamped BCs on all four sides.
Results¶
\(n_x = n_y\) |
\(\Delta x\) [km] |
\(L^\infty\) error |
|---|---|---|
26 |
24.0 |
1.86 × 10⁻³ |
51 |
12.0 |
4.98 × 10⁻⁴ |
101 |
6.0 |
1.26 × 10⁻⁴ |
201 |
3.0 |
3.16 × 10⁻⁵ |
Convergence slopes (coarse → fine): \(\mathcal{O}(\Delta x^{1.90})\), \(\mathcal{O}(\Delta x^{1.98})\), \(\mathcal{O}(\Delta x^{2.00})\).
The test function tests.test_bc_mms.TestVariableTe2D.test_convergence_order()
runs with every CI build.