Convergence lab

How to read it

The errors are relative L2 norms over the whole domain, inclusion included: ‖f_h − f‖ / ‖f‖ with ‖g‖² = ∫ g² dA, integrated element by element with the same 7-point, degree-5 quadrature rule the assembly uses, the finite element velocity (nodal plus bubble) and stress evaluated at the quadrature points against the exact fields there, with the analytical branch continued smoothly through each element. In the sliver between the polygonal interface and the true circle, a separate polar Gauss rule subtracts the continued matrix-side contribution and adds the correct inclusion-side contribution. This accounts for the stress jump without counting the sliver twice. The order column is the observed convergence rate: when halving the element size h divides the error by 2^p, the method is of order p.

For the quadratic velocity of MILAMIN’s 7-node Crouzeix–Raviart element, theory promises third-order velocities and second-order stresses in these norms. The velocity measures ≈2 here, and that is the first lesson of this page: the circle exists in the mesh only as straight element edges, and the polygonal approximation of the curved interface is itself an O(h²) error that caps the velocity order. The control experiment proves it: on a square domain, where straight edges are exact, the same element measures its full O(h³·⁰) velocity order – with constant-pressure and rotational (non-trivial pressure) manufactured solutions alike (scripts/convtest.ts in the repository). Recovering cubic order on curved boundaries is what curved (isoparametric) element edges are for; MILAMIN's meshes from Triangle have straight edges, so its isoparametric mapping is affine, and like most practical codes it trades that last order for refinement.

The stress tells the second lesson: it converges at first order, not second, and the table says why. The maximum shear stress jumps by a factor of nearly the viscosity ratio across the interface, and in the sliver between the polygon and the true circle the mesh has the wrong material. The stress error there is O(1) regardless of refinement; the sliver's area shrinks as O(h²), so its contribution to the L2 error falls only as O(h), and from the third level on it is more than 90% of the squared stress error. Away from the interface the stress converges at its natural second order. A velocity is continuous across the interface, which is why it does not suffer the same way: the price of a viscosity jump on a polygonal interface is paid exactly where the error map lights up.

The error map tells the spatial story: the solution is most exact far from the inclusion and the error concentrates at the interface, where the polygonal approximation of the circle, the viscosity jump and the pressure jump all live. The meshes of this study are refined uniformly, the interface only through the number of polygon vertices, so the map is a fair picture of where refinement would pay.

Two things are being checked at once here. Verification: the discrete solution approaches the exact one, so the assembly, solver and boundary conditions are implemented correctly – the same analytical solutions guard this codebase in its verification suite, 46 checks, run before every release. And discretization economics: a level in the table halves the element size h, which in two dimensions quadruples the number of elements and unknowns; what that buys is in the order columns, four times smaller velocity errors and two times smaller stress errors per level.

Reference: Dabrowski, M., Krotkiewski, M. & Schmid, D. W. (2008), MILAMIN: MATLAB-based finite element method solver for large problems, Geochemistry, Geophysics, Geosystems 9, Q04030, doi:10.1029/2007GC001719, section 4.