Tutorial: 2D Poisson with volumential + pytential

This notebook is a practical tutorial for new volumential users. We solve a smooth Poisson problem on a nontrivial starfish domain using:

• a volume potential on box-quadrature points,

• harmonic continuation of source values outside the physical domain,

• a harmonic boundary correction to enforce Dirichlet data.

Because we use a manufactured exact solution, every stage can be validated with concrete error metrics.

Roadmap

What this tutorial covers:

1. Problem setup: geometry, exact solution, and forcing term.

2. One-time infrastructure: OpenCL context and near-field table cache.

3. Reusable solver: a single run_poisson_case(...) driver.

4. Single run diagnostics: inspect geometry, source extension behavior, and error maps.

5. Uniform co-refinement: refine volume + boundary + QBX together and measure convergence.

6. Result interpretation: summarize convergence and solver health.

7. Boundary-focused AMR: refine only the volume boxes cut by the boundary, and compare error per source point and where the error sits against uniform refinement, at the probes and at the volume nodes.

Tip: this notebook is designed so you can copy run_poisson_case(...) into your own project and swap in your geometry and boundary data.

from __future__ import annotations

from functools import partial
import numpy as np
from scipy.spatial import cKDTree

import pyopencl as cl
import pyopencl.array

from arraycontext import flatten, unflatten
from meshmode.dof_array import DOFArray
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import InterpolatoryQuadratureSimplexGroupFactory
from meshmode.mesh.generation import make_curve_mesh, starfish

from boxtree.array_context import PyOpenCLArrayContext as BoxtreePyOpenCLArrayContext
from boxtree.traversal import FMMTraversalBuilder

from pytools.obj_array import new_1d as obj_array_1d

from pytential.array_context import PyOpenCLArrayContext
from pytential.qbx import QBXLayerPotentialSource
from pytential.target import PointsTarget

import volumential.meshgen as mg
from volumential.expansion_wrangler_fpnd import (
    FPNDExpansionWrangler,
    FPNDTreeIndependentDataForWrangler,
)
from volumential.function_extension import compute_harmonic_extension
from volumential.nearfield_potential_table import DuffyBuildConfig
from volumential.table_manager import NearFieldInteractionTableManager
from volumential.tree_interactive_build import build_particle_tree_from_box_tree
from volumential.volume_fmm import drive_volume_fmm, interpolate_volume_potential

from sumpy.expansion import DefaultExpansionFactory
from sumpy.kernel import LaplaceKernel

1) Problem Definition

We build a manufactured-solution test case so the true error is known everywhere.

• Domain boundary: starfish(t) from meshmode.

• Exact solution: smooth oscillatory modes plus a Gaussian bump.

• Forcing: f = -laplacian(u_exact).

• Evaluation set: fixed interior probe points inside the physical starfish domain (probe_scale=1.00 by default) so boundary-region error is reflected in reported metrics.

from matplotlib.path import Path

dim = 2
bbox_a, bbox_b = -1.30, 1.30
q_order = 4

table_filename = 'nft_poisson2d_starfish_demo.sqlite'
table_root_extent = bbox_b - bbox_a

bdry_order = 4
gmres_tolerance = 1.0e-8


def u_exact(x, y):
    phase1 = 3.0 * x + 2.0 * y
    phase2 = 4.0 * x - 1.5 * y
    phase3 = 2.5 * x - 3.5 * y

    smooth_wave = 0.55 * np.sin(phase1) + 0.35 * np.cos(phase2)
    modulated_wave = 0.25 * np.exp(0.4 * x - 0.3 * y) * np.sin(phase3)

    bump_r2 = (x - 0.25) ** 2 + (y + 0.15) ** 2
    bump = 0.30 * np.exp(-5.0 * bump_r2)

    return smooth_wave + modulated_wave + bump


def laplacian_u_exact(x, y):
    phase1 = 3.0 * x + 2.0 * y
    phase2 = 4.0 * x - 1.5 * y

    term1 = -0.55 * (3.0**2 + 2.0**2) * np.sin(phase1)
    term2 = -0.35 * (4.0**2 + 1.5**2) * np.cos(phase2)

    a, b = 0.4, -0.3
    c, d = 2.5, -3.5
    phase3 = c * x + d * y
    envelope = np.exp(a * x + b * y)
    term3 = 0.25 * envelope * (
        (a * a + b * b - c * c - d * d) * np.sin(phase3)
        + 2.0 * (a * c + b * d) * np.cos(phase3)
    )

    sigma = 5.0
    bump_r2 = (x - 0.25) ** 2 + (y + 0.15) ** 2
    term4 = 0.30 * (4.0 * sigma * sigma * bump_r2 - 4.0 * sigma) * np.exp(-sigma * bump_r2)

    return term1 + term2 + term3 + term4


def rhs_f(x, y):
    return -laplacian_u_exact(x, y)


def make_starfish_path(scale=1.0, npts=4096):
    t = np.linspace(0.0, 1.0, npts, endpoint=False)
    curve = starfish(t)
    vertices = np.ascontiguousarray(np.vstack([scale * curve[0], scale * curve[1]]).T)
    closed_vertices = np.vstack([vertices, vertices[0]])
    return Path(closed_vertices), closed_vertices


starfish_path, starfish_vertices = make_starfish_path(scale=1.0, npts=4096)
probe_scale = 1.00
probe_path, _ = make_starfish_path(scale=probe_scale, npts=4096)

probe_axis = np.linspace(-1.15, 1.15, 65)
probe_xx, probe_yy = np.meshgrid(probe_axis, probe_axis, indexing='xy')
probe_candidates = np.ascontiguousarray(np.vstack([probe_xx.ravel(), probe_yy.ravel()]).T)
probe_mask = probe_path.contains_points(probe_candidates)
probe_points_host = np.ascontiguousarray(probe_candidates[probe_mask])

print(f'fixed interior probe points (starfish scale={probe_scale:.2f}):', probe_points_host.shape[0])

2) Infrastructure Setup

Create OpenCL/QBX contexts and build or load the near-field interaction table once.

This step is intentionally separate so repeated parameter studies can reuse the same cached table.

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
actx = PyOpenCLArrayContext(queue)

probe_target = PointsTarget(actx.freeze(actx.from_numpy(np.ascontiguousarray(probe_points_host.T))))

tm = NearFieldInteractionTableManager(
    table_filename,
    root_extent=table_root_extent,
    queue=queue,
)
build_config = DuffyBuildConfig(
    radial_rule='tanh-sinh-fast',
    regular_quad_order=8,
    radial_quad_order=21,
)

nftable, _ = tm.get_table(
    dim,
    'Laplace',
    q_order,
    force_recompute=False,
    queue=queue,
    build_config=build_config,
)

print('near-field table ready:', table_filename)

3) Solver Module

run_poisson_case(...) executes one full Poisson solve and returns diagnostics.

Pipeline per run:

1. Build boundary discretization + QBX objects.

2. Build (uniform or boundary-refined) volume mesh and quadrature points.

3. Evaluate interior source values from f and harmonically continue exterior source values.

4. Compute volume potential with volumential FMM.

5. Solve harmonic boundary correction for g - v|_boundary.

6. Report manufactured-solution errors and GMRES states.

This is the function to reuse when you run your own parameter sweeps.

def build_boundary_qbx(nelements_bdry, qbx_order):
    bdry_mesh = make_curve_mesh(
        starfish,
        np.linspace(0.0, 1.0, nelements_bdry + 1),
        bdry_order,
    )

    density_discr = Discretization(
        actx,
        bdry_mesh,
        InterpolatoryQuadratureSimplexGroupFactory(bdry_order),
    )

    qbx = QBXLayerPotentialSource(
        density_discr,
        fine_order=4 * bdry_order,
        qbx_order=qbx_order,
        fmm_order=False,
    )

    bdry_nodes = actx.thaw(density_discr.nodes())

    flat_bdry_nodes = flatten(bdry_nodes, actx, leaf_class=DOFArray)
    bdry_x_host = actx.to_numpy(flat_bdry_nodes[0])
    bdry_y_host = actx.to_numpy(flat_bdry_nodes[1])
    bdry_points_host = np.ascontiguousarray(np.vstack([bdry_x_host, bdry_y_host]))

    g_bdry_flat_host = np.ascontiguousarray(u_exact(bdry_x_host, bdry_y_host))
    g_bdry_flat = actx.from_numpy(g_bdry_flat_host)
    g_bdry = unflatten(bdry_nodes[0], g_bdry_flat, actx)

    return qbx, density_discr, bdry_nodes, g_bdry, bdry_points_host


def refine_volume_mesh_near_boundary(
    volume_mesh,
    boundary_path,
    boundary_vertices,
    *,
    max_level=10,
    max_passes=32,
    near_boundary_factor=0.20,
):
    # `max_level` is a tree level, with the root box at level 0: leaves at
    # `max_level` are not refined further.
    boundary_vertices = np.ascontiguousarray(boundary_vertices[:-1])
    boundary_kdtree = cKDTree(boundary_vertices)

    n_refined_total = 0
    n_passes = 0

    for _ in range(max_passes):
        leaf_boxes = volume_mesh._leaf_boxes()
        leaf_levels = volume_mesh._leaf_levels()

        if len(leaf_boxes) == 0:
            break

        eligible = leaf_levels < max_level
        if not np.any(eligible):
            break

        centers = volume_mesh._leaf_centers()
        sizes = volume_mesh._leaf_side_lengths()
        half = 0.5 * sizes

        x = centers[:, 0]
        y = centers[:, 1]

        corners = np.ascontiguousarray(
            np.vstack(
                [
                    np.column_stack([x - half, y - half]),
                    np.column_stack([x + half, y - half]),
                    np.column_stack([x + half, y + half]),
                    np.column_stack([x - half, y + half]),
                ]
            )
        )

        edge_midpoints = np.ascontiguousarray(
            np.vstack(
                [
                    np.column_stack([x, y - half]),
                    np.column_stack([x + half, y]),
                    np.column_stack([x, y + half]),
                    np.column_stack([x - half, y]),
                ]
            )
        )

        stencil_points = np.ascontiguousarray(
            np.vstack([corners, edge_midpoints, centers])
        )
        stencil_inside = boundary_path.contains_points(stencil_points).reshape(9, -1).T
        stencil_mixed = np.any(stencil_inside, axis=1) & (~np.all(stencil_inside, axis=1))

        dist_to_boundary, _ = boundary_kdtree.query(centers)
        near_boundary = dist_to_boundary <= (
            near_boundary_factor * np.sqrt(2.0) * sizes + 1.0e-14
        )

        refine_leaf_mask = eligible & (stencil_mixed | near_boundary)
        if not np.any(refine_leaf_mask):
            break

        refine_flags = np.zeros(volume_mesh.boxtree.nboxes, dtype=bool)
        refine_flags[leaf_boxes[refine_leaf_mask]] = True
        coarsen_flags = np.zeros_like(refine_flags)

        n_refined_total += int(np.count_nonzero(refine_leaf_mask))
        n_passes += 1

        volume_mesh.boxtree.refine_and_coarsen(
            refine_flags=refine_flags,
            coarsen_flags=coarsen_flags,
            error_on_ignored_flags=False,
        )

    leaf_levels = volume_mesh._leaf_levels()
    max_leaf_level = int(np.max(leaf_levels)) if len(leaf_levels) else -1

    return {
        'passes': int(n_passes),
        'n_refined_total': int(n_refined_total),
        'n_active_cells': int(volume_mesh.n_active_cells()),
        'max_leaf_level': max_leaf_level,
    }


# `run_poisson_case(node_errors=True)` also measures the error at the interior
# volume nodes, leaving out those closer than `node_boundary_clearance` to the
# boundary. The distance is taken to the starfish sampled finely enough for it
# to be accurate to about 1e-4.
node_boundary_clearance = 1.0e-3
starfish_fine_kdtree = cKDTree(
    np.ascontiguousarray(np.array(starfish(np.linspace(0.0, 1.0, 2**16, endpoint=False))).T)
)


def run_poisson_case(
    *,
    vol_nlevels,
    nelements_bdry,
    qbx_order,
    fmm_order,
    adaptive_boundary_max_level=None,
    adaptive_max_passes=32,
    adaptive_near_boundary_factor=0.20,
    node_errors=False,
    return_fields=False,
    return_mesh=False,
):
    qbx, density_discr, bdry_nodes, g_bdry, bdry_points_host = build_boundary_qbx(
        nelements_bdry=nelements_bdry,
        qbx_order=qbx_order,
    )

    volume_mesh = mg.MeshGen2D(
        degree=q_order,
        nlevels=vol_nlevels,
        a=bbox_a,
        b=bbox_b,
        queue=queue,
    )

    if adaptive_boundary_max_level is not None and adaptive_boundary_max_level > vol_nlevels:
        # `adaptive_boundary_max_level` counts levels the way `vol_nlevels`
        # does: MeshGen2D(nlevels=n) has its leaves at tree level n - 1, so the
        # boundary leaves are refined down to tree level
        # adaptive_boundary_max_level - 1, the leaf size of a uniform mesh
        # with vol_nlevels=adaptive_boundary_max_level.
        mesh_info = refine_volume_mesh_near_boundary(
            volume_mesh,
            starfish_path,
            starfish_vertices,
            max_level=adaptive_boundary_max_level - 1,
            max_passes=adaptive_max_passes,
            near_boundary_factor=adaptive_near_boundary_factor,
        )
    else:
        leaf_levels = volume_mesh._leaf_levels()
        mesh_info = {
            'passes': 0,
            'n_refined_total': 0,
            'n_active_cells': int(volume_mesh.n_active_cells()),
            'max_leaf_level': int(np.max(leaf_levels)) if len(leaf_levels) else -1,
        }

    source_points_host = np.ascontiguousarray(volume_mesh.get_q_points())
    source_weights_host = np.ascontiguousarray(volume_mesh.get_q_weights())

    rhs_full_host = rhs_f(source_points_host[:, 0], source_points_host[:, 1])
    interior_source_mask = starfish_path.contains_points(source_points_host)
    exterior_source_mask = ~interior_source_mask

    source_vals_host = np.ascontiguousarray(rhs_full_host.copy())

    rhs_bdry_flat_host = np.ascontiguousarray(rhs_f(bdry_points_host[0], bdry_points_host[1]))
    rhs_bdry = unflatten(
        bdry_nodes[0],
        actx.from_numpy(rhs_bdry_flat_host),
        actx,
    )

    gmres_rhs_ext_state = 'not-needed'
    if np.any(exterior_source_mask):
        exterior_source_points_host = np.ascontiguousarray(source_points_host[exterior_source_mask])
        exterior_source_target = PointsTarget(
            actx.freeze(actx.from_numpy(np.ascontiguousarray(exterior_source_points_host.T)))
        )

        rhs_exterior, dbg_rhs_ext = compute_harmonic_extension(
            queue,
            exterior_source_target,
            qbx,
            density_discr,
            rhs_bdry,
            loc_sign=+1,
            representation_mode='auto',
            target_association_tolerance=0.05,
            gmres_tolerance=gmres_tolerance,
            actx=actx,
        )
        source_vals_host[exterior_source_mask] = actx.to_numpy(rhs_exterior)
        gmres_rhs_ext_state = dbg_rhs_ext['gmres_result'].state

    all_targets_host = np.ascontiguousarray(
        np.hstack([probe_points_host.T, bdry_points_host])
    )

    target_points = obj_array_1d(
        [
            cl.array.to_device(queue, np.ascontiguousarray(all_targets_host[iaxis]))
            for iaxis in range(dim)
        ]
    )

    source_vals = cl.array.to_device(queue, np.ascontiguousarray(source_vals_host))
    source_weights = cl.array.to_device(queue, source_weights_host)

    bt_actx = BoxtreePyOpenCLArrayContext(queue)
    source_tree = build_particle_tree_from_box_tree(
        bt_actx,
        volume_mesh.boxtree,
        source_points_host,
    )

    trav_builder = FMMTraversalBuilder(bt_actx)
    source_trav, _ = trav_builder(bt_actx, source_tree)

    knl = LaplaceKernel(dim)
    expn_factory = DefaultExpansionFactory()
    local_expn_class = expn_factory.get_local_expansion_class(knl)
    mpole_expn_class = expn_factory.get_multipole_expansion_class(knl)

    tree_indep = FPNDTreeIndependentDataForWrangler(
        ctx,
        partial(mpole_expn_class, knl),
        partial(local_expn_class, knl),
        [knl],
        exclude_self=False,
    )

    wrangler = FPNDExpansionWrangler(
        tree_indep=tree_indep,
        traversal=source_trav,
        near_field_table=nftable,
        dtype=np.float64,
        fmm_level_to_order=lambda kernel, kernel_args, tree, lev: fmm_order,
        quad_order=q_order,
        queue=queue,
    )

    (vf_source,) = drive_volume_fmm(
        source_trav,
        wrangler,
        source_vals * source_weights,
        source_vals,
    )

    vf_all = interpolate_volume_potential(
        target_points,
        source_trav,
        wrangler,
        vf_source,
    )
    vf_all_host = vf_all.get() if hasattr(vf_all, 'get') else np.asarray(vf_all)
    n_probe = probe_points_host.shape[0]

    vf_probe_host = vf_all_host[:n_probe]
    vf_bdry_host = vf_all_host[n_probe:]

    g_bdry_flat_host = np.ascontiguousarray(u_exact(bdry_points_host[0], bdry_points_host[1]))
    corr_bdry_flat = actx.from_numpy(np.ascontiguousarray(g_bdry_flat_host - vf_bdry_host))
    corr_bdry = unflatten(bdry_nodes[0], corr_bdry_flat, actx)

    corr_probe, dbg_corr = compute_harmonic_extension(
        queue,
        probe_target,
        qbx,
        density_discr,
        corr_bdry,
        loc_sign=-1,
        representation_mode='auto',
        target_association_tolerance=0.05,
        gmres_tolerance=gmres_tolerance,
        actx=actx,
    )
    corr_probe_host = actx.to_numpy(corr_probe)

    u_num_host = vf_probe_host + corr_probe_host
    u_ref_host = u_exact(probe_points_host[:, 0], probe_points_host[:, 1])

    abs_err = np.abs(u_num_host - u_ref_host)
    rel_l2 = np.linalg.norm(abs_err) / np.linalg.norm(u_ref_host)

    if node_errors:
        # The solution at the interior volume nodes, where the volume potential
        # comes straight from the FMM; at the probes it is interpolated from
        # the nodes of the leaf that holds the probe. The correction is
        # evaluated directly at both, from the same boundary density.
        node_dist, _ = starfish_fine_kdtree.query(source_points_host)
        node_mask = interior_source_mask & (node_dist > node_boundary_clearance)
        node_points_host = np.ascontiguousarray(source_points_host[node_mask])
        node_target = PointsTarget(
            actx.freeze(actx.from_numpy(np.ascontiguousarray(node_points_host.T)))
        )
        corr_node = dbg_corr['eval_ext_f'](node_target)
        vf_source_host = vf_source.get() if hasattr(vf_source, 'get') else np.asarray(vf_source)
        u_node_ref_host = u_exact(node_points_host[:, 0], node_points_host[:, 1])
        node_abs_err = np.abs(
            vf_source_host[node_mask] + actx.to_numpy(corr_node) - u_node_ref_host
        )
        # relative L2 error with the mesh's own quadrature weights
        node_weights_host = source_weights_host[node_mask]
        node_rel_l2 = np.sqrt(
            np.sum(node_weights_host * node_abs_err**2)
            / np.sum(node_weights_host * u_node_ref_host**2)
        )

    result = {
        'vol_nlevels': vol_nlevels,
        'nelements_bdry': nelements_bdry,
        'qbx_order': qbx_order,
        'fmm_order': fmm_order,
        # leaf side of the mesh before any adaptive refinement
        'h': (bbox_b - bbox_a) / (2 ** (vol_nlevels - 1)),
        # side of the smallest leaf after adaptive refinement
        'h_min': (bbox_b - bbox_a) / (2 ** int(mesh_info['max_leaf_level'])),
        'adaptive_boundary_max_level': adaptive_boundary_max_level,
        'adaptive_near_boundary_factor': float(adaptive_near_boundary_factor),
        'n_sources': int(source_points_host.shape[0]),
        'n_cells': int(mesh_info['n_active_cells']),
        'mesh_max_leaf_level': int(mesh_info['max_leaf_level']),
        'mesh_refine_passes': int(mesh_info['passes']),
        'mesh_refined_leaf_boxes': int(mesh_info['n_refined_total']),
        'n_source_interior': int(np.count_nonzero(interior_source_mask)),
        'n_source_exterior': int(np.count_nonzero(exterior_source_mask)),
        'n_probe': int(probe_points_host.shape[0]),
        'rel_l2_err': float(rel_l2),
        'max_abs_err': float(abs_err.max()),
        'p95_abs_err': float(np.percentile(abs_err, 95)),
        'gmres_rhs_ext_state': gmres_rhs_ext_state,
        'gmres_corr_state': dbg_corr['gmres_result'].state,
    }

    if node_errors:
        result.update(
            {
                'n_node': int(node_points_host.shape[0]),
                'node_rel_l2_err': float(node_rel_l2),
            }
        )
        if return_fields:
            result.update(
                {
                    'node_points_host': node_points_host,
                    'node_abs_err': node_abs_err,
                }
            )

    if return_fields:
        result.update(
            {
                'source_points_host': source_points_host,
                'source_vals_host': source_vals_host,
                'u_num_host': u_num_host,
                'u_ref_host': u_ref_host,
                'abs_err': abs_err,
                'bdry_points_host': bdry_points_host,
            }
        )

    if return_mesh:
        result.update(
            {
                'mesh_leaf_centers_host': np.ascontiguousarray(volume_mesh._leaf_centers()),
                'mesh_leaf_side_lengths_host': np.ascontiguousarray(volume_mesh._leaf_side_lengths()),
            }
        )

    return result

4) Single-Resolution Walkthrough

Run one representative configuration and inspect:

• geometry and sampling locations,

• interior/exterior source split,

• where harmonic continuation changes source values,

• numerical solution and interior error distribution.

demo_result = run_poisson_case(
    vol_nlevels=4,
    nelements_bdry=192,
    qbx_order=4,
    fmm_order=14,
    return_fields=True,
)

n_src_total = demo_result['n_sources']
n_src_ext = demo_result['n_source_exterior']
ext_ratio = (n_src_ext / n_src_total) if n_src_total else 0.0

print('demo configuration:')
print(
    '  vol_nlevels=', demo_result['vol_nlevels'],
    'nelements_bdry=', demo_result['nelements_bdry'],
    'qbx_order=', demo_result['qbx_order'],
    'fmm_order=', demo_result['fmm_order'],
)
print('  n_sources=', demo_result['n_sources'], 'n_cells=', demo_result['n_cells'], 'n_probe=', demo_result['n_probe'])
print('  mesh max leaf level=', demo_result['mesh_max_leaf_level'], 'refine passes=', demo_result['mesh_refine_passes'])
print('  source split (in, out)=', demo_result['n_source_interior'], demo_result['n_source_exterior'])
print(f'  exterior-source fraction={ext_ratio:.3f}')
print('  gmres states (rhs-ext, corr)=', demo_result['gmres_rhs_ext_state'], demo_result['gmres_corr_state'])
print('  rel_l2_err=', demo_result['rel_l2_err'])
print('  max_abs_err=', demo_result['max_abs_err'])
print('  p95_abs_err=', demo_result['p95_abs_err'])

import matplotlib.pyplot as plt

source_points = demo_result['source_points_host']
source_vals = demo_result['source_vals_host']
u_num = demo_result['u_num_host']
u_ref = demo_result['u_ref_host']
abs_err = demo_result['abs_err']

interior_src_mask = starfish_path.contains_points(source_points)
exterior_src_mask = ~interior_src_mask
rhs_true_on_sources = rhs_f(source_points[:, 0], source_points[:, 1])
rhs_ext_delta = source_vals - rhs_true_on_sources

_, probe_vertices_vis = make_starfish_path(scale=probe_scale, npts=4096)

fig, ax = plt.subplots(2, 3, figsize=(16, 9.5), constrained_layout=True)

ax00 = ax[0, 0]
ax00.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.6, label='domain boundary')
ax00.plot(probe_vertices_vis[:, 0], probe_vertices_vis[:, 1], color='tab:green', lw=1.1, ls='--', label=f'probe contour ({probe_scale:.2f}x)')
ax00.scatter(probe_points_host[:, 0], probe_points_host[:, 1], s=5, color='tab:blue', alpha=0.30, label='probe points')
ax00.set_title('Domain and Probe Set')
ax00.set_aspect('equal')
ax00.legend(loc='upper right', fontsize=8)

ax01 = ax[0, 1]
ax01.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.0)
ax01.scatter(source_points[interior_src_mask, 0], source_points[interior_src_mask, 1], s=5, alpha=0.35, color='tab:orange', label='interior sources')
ax01.scatter(source_points[exterior_src_mask, 0], source_points[exterior_src_mask, 1], s=5, alpha=0.35, color='tab:red', label='exterior sources (extended)')
ax01.set_title('Volume Source Point Split')
ax01.set_aspect('equal')
ax01.legend(loc='upper right', fontsize=8)

ax02 = ax[0, 2]
sc_rhs = ax02.scatter(source_points[:, 0], source_points[:, 1], c=rhs_ext_delta, s=8, cmap='coolwarm')
ax02.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.0)
ax02.set_title(r'RHS continuation delta: $f_{used} - f_{analytic}$')
ax02.set_aspect('equal')
cb_rhs = fig.colorbar(sc_rhs, ax=ax02, shrink=0.84)
cb_rhs.set_label('delta')

ax10 = ax[1, 0]
sc_u = ax10.scatter(probe_points_host[:, 0], probe_points_host[:, 1], c=u_num, s=10, cmap='viridis')
ax10.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.0)
ax10.set_title('Numerical Solution on Probe Set')
ax10.set_aspect('equal')
cb_u = fig.colorbar(sc_u, ax=ax10, shrink=0.84)
cb_u.set_label('u_num')

ax11 = ax[1, 1]
sc_e = ax11.scatter(probe_points_host[:, 0], probe_points_host[:, 1], c=np.log10(abs_err + 1.0e-14), s=10, cmap='magma')
ax11.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.0)
ax11.set_title(r'Probe error map: $\log_{10}(|u_{num}-u_{exact}|)$')
ax11.set_aspect('equal')
cb_e = fig.colorbar(sc_e, ax=ax11, shrink=0.84)
cb_e.set_label('log10 abs err')

ax12 = ax[1, 2]
ax12.hist(abs_err, bins=40, color='tab:purple', alpha=0.85)
ax12.axvline(np.percentile(abs_err, 50), color='tab:blue', ls='--', lw=1.5, label='p50')
ax12.axvline(np.percentile(abs_err, 95), color='tab:red', ls='--', lw=1.5, label='p95')
ax12.set_title('Probe Error Distribution')
ax12.set_xlabel('absolute error')
ax12.set_ylabel('count')
ax12.legend(loc='upper right', fontsize=8)

fig.suptitle(
    'Single-Resolution Diagnostics'
    + f"\nrel L2={demo_result['rel_l2_err']:.3e}, max={demo_result['max_abs_err']:.3e}, p95={demo_result['p95_abs_err']:.3e}",
    fontsize=12,
)

plt.show()

5) Uniform Co-Refinement Study

To estimate asymptotic behavior, refine together:

• volume mesh level (vol_nlevels),

• boundary element count (nelements_bdry),

• QBX order (qbx_order),

• FMM order (fmm_order, mildly increased).

This keeps one subsystem from becoming the dominant error floor too early.

The loop starts at co_refinement_start_level and stops at co_refinement_max_level, or earlier once the relative L2 error is below co_refinement_plateau_target (1e-10) and one more level gains less than co_refinement_plateau_log10_delta (0.1 decade). At the default settings the error stops improving at about 5e-9 (level 9), above that target, so the loop runs to level 10; this cell takes most of the notebook’s run time.

co_refinement_start_level = 3
co_refinement_max_level = 10
co_refinement_plateau_target = 1.0e-10
co_refinement_min_runs = 4
co_refinement_plateau_log10_delta = 0.10  # <=0.10 decade improvement (~<=26%) treated as plateau


def co_refinement_cfg_for_level(level):
    step = level - co_refinement_start_level
    return {
        'vol_nlevels': level,
        'nelements_bdry': int(96 * (2**step)),
        'qbx_order': min(3 + step, 9),
        'fmm_order': 12 + 2 * step,
    }


co_results = []

print(f'Co-refinement (harmonic correction) on fixed interior probe set (starfish scaled by {probe_scale:.2f}):')
print(f'- Start level={co_refinement_start_level}, max level={co_refinement_max_level}')
print(f'- Plateau target: rel L2 <= {co_refinement_plateau_target:.1e} with < {co_refinement_plateau_log10_delta:.2f} decade gain')
print(' level   n_bdry  qbx fmm      h      n_src   rel_l2_err   max_abs_err    p95_abs   rate  gmres')

prev_err = None
plateau_reached = False
plateau_reason = ''

for level in range(co_refinement_start_level, co_refinement_max_level + 1):
    cfg = co_refinement_cfg_for_level(level)
    row = run_poisson_case(**cfg, return_fields=False)
    co_results.append(row)

    err = row['rel_l2_err']
    rate = np.log2(prev_err / err) if (prev_err is not None and err > 0.0) else np.nan
    rate_str = f'{rate:5.2f}' if np.isfinite(rate) else '  -  '

    gmres_ok = (row['gmres_rhs_ext_state'] == 'success' and row['gmres_corr_state'] == 'success')

    print(
        f"{row['vol_nlevels']:6d}"
        f" {row['nelements_bdry']:8d}"
        f" {row['qbx_order']:4d}"
        f" {row['fmm_order']:3d}"
        f" {row['h']:8.5f}"
        f" {row['n_sources']:8d}"
        f" {row['rel_l2_err']:12.4e}"
        f" {row['max_abs_err']:12.4e}"
        f" {row['p95_abs_err']:10.3e}"
        f" {rate_str:>5s}"
        f" {'ok' if gmres_ok else 'fail'}"
    )

    log10_gain = np.log10(prev_err) - np.log10(err) if (prev_err is not None and prev_err > 0.0 and err > 0.0) else np.nan

    if (
        len(co_results) >= co_refinement_min_runs
        and err <= co_refinement_plateau_target
        and np.isfinite(log10_gain)
        and log10_gain < co_refinement_plateau_log10_delta
    ):
        plateau_reached = True
        plateau_reason = (
            f"level {level}: rel L2={err:.3e}, decade gain={log10_gain:.3f} < {co_refinement_plateau_log10_delta:.3f}"
        )
        break

    prev_err = err

if plateau_reached:
    print(f'Plateau reached: {plateau_reason}')
elif co_results and co_results[-1]['rel_l2_err'] <= co_refinement_plateau_target:
    print(
        f"Reached rel L2 <= {co_refinement_plateau_target:.1e} at level {co_results[-1]['vol_nlevels']}, "
        'but plateau criterion was not met before max level.'
    )
else:
    print(
        f"Did not reach rel L2 <= {co_refinement_plateau_target:.1e} by level {co_refinement_max_level}."
    )

hs = np.array([row['h'] for row in co_results])
errs = np.array([row['rel_l2_err'] for row in co_results])
max_errs = np.array([row['max_abs_err'] for row in co_results])
nsrc = np.array([row['n_sources'] for row in co_results])
levels = np.array([row['vol_nlevels'] for row in co_results])

if len(errs) >= 2 and np.all(errs > 0.0):
    slope, intercept = np.polyfit(np.log(hs), np.log(errs), 1)
    print(f'Harmonic-correction log-log slope (error vs h): {slope:.3f}')
else:
    slope = np.nan
    print('Not enough positive error data for slope estimate.')

fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), constrained_layout=True)

ax0 = axes[0]
ax0.loglog(hs, errs, 'o-', lw=2.0, ms=8, label='rel L2 error')
ax0.loglog(hs, max_errs, 's--', lw=1.6, ms=6, label='max abs error')
ax0.axhline(co_refinement_plateau_target, color='tab:red', ls='--', lw=1.0, alpha=0.65, label='1e-10 target')
if np.isfinite(slope):
    fit = np.exp(intercept) * hs**slope
    ax0.loglog(hs, fit, ':', lw=1.6, label=f'fit slope={slope:.2f}')
ax0.set_xlabel('h (leaf side of the volume mesh)')
ax0.set_ylabel('error')
ax0.set_title('Error vs h (uniform co-refinement)')
ax0.grid(True, which='both', ls=':', alpha=0.5)
ax0.legend()

ax1 = axes[1]
ax1.loglog(nsrc, errs, 'o-', lw=2.0, ms=8, color='tab:blue')
for x, y, lev in zip(nsrc, errs, levels):
    ax1.annotate(f'l{lev}', (x, y), textcoords='offset points', xytext=(4, 4), fontsize=8)
ax1.set_xlabel('number of volume source points')
ax1.set_ylabel('relative L2 error')
ax1.set_title('Accuracy-cost trend')
ax1.grid(True, which='both', ls=':', alpha=0.5)

plt.show()

6) Interpreting Convergence Results

Use the diagnostics below to answer practical questions:

• Do errors decrease consistently under co-refinement?

• Is GMRES robust across runs?

• What is the observed error-vs-resolution slope?

• Is the final accuracy good enough for your target workflow?

best = min(co_results, key=lambda row: row['rel_l2_err'])
last = co_results[-1]
first = co_results[0]

overall_reduction = first['rel_l2_err'] / last['rel_l2_err'] if last['rel_l2_err'] > 0 else np.inf

print('Summary observations:')
print(f"- Best co-refined run: level={best['vol_nlevels']} with rel L2={best['rel_l2_err']:.4e}")
print(f"- Relative L2 reduction (first -> last): {overall_reduction:.2f}x")
print(f"- Final relative L2 error: {last['rel_l2_err']:.4e}")
print(f"- Final max abs error: {last['max_abs_err']:.4e}")

all_gmres_ok = all(
    row['gmres_rhs_ext_state'] == 'success' and row['gmres_corr_state'] == 'success'
    for row in co_results
)
print(f"- GMRES successful on all co-refinement runs: {all_gmres_ok}")

if plateau_reached:
    print(f"- Plateau stop condition reached: {plateau_reason}")
else:
    print('- Plateau stop condition not reached before max level cap.')

if np.isfinite(slope):
    if slope > 3.0:
        slope_note = 'strong high-order trend for this smooth manufactured case'
    elif slope > 1.0:
        slope_note = 'clear convergence trend'
    elif slope > 0.0:
        slope_note = 'convergent but moderate asymptotic rate'
    else:
        slope_note = 'non-convergent trend (investigate resolution balance)'
    print(f"- Fitted slope interpretation: {slope:.3f} -> {slope_note}")

7) Boundary-Focused AMR Study

The harmonic continuation in run_poisson_case (step 3 of its pipeline) makes the source continuous across the boundary but not smooth there: its normal derivative jumps. On a leaf that the boundary cuts, the polynomial of degree q_order - 1 in each variable that represents the source on the leaf cannot follow that kink; on every other leaf the source is smooth. Refining only the cut leaves goes after that one weakness, and this section measures what it buys.

AMR marker used here:

• a leaf is cut when its 9-point stencil (corners, edge midpoints, center) has points both inside and outside the starfish;

• each pass splits the cut leaves that are still coarser than the target, and the tree is kept 2:1 balanced, so some uncut neighbors are split as well;

• adaptive_boundary_max_level counts levels the way vol_nlevels does: amr-l5-to7 starts from the vol_nlevels=5 mesh and refines the cut leaves down to the leaf size of a uniform vol_nlevels=7 mesh;

• adaptive_near_boundary_factor=0.0 turns off the distance-based marker, so only cut leaves are marked.

The cases are uniform meshes with vol_nlevels 5, 6 and 7, and the vol_nlevels=5 mesh with its cut leaves refined to levels 6, 7 and 8. All six use the same boundary discretization (768 elements, QBX order 6, FMM order 18), so only the volume mesh changes from one case to the next.

The error is measured in two places. At the probes, the full interior probe set (probe_scale=1.00), the volume potential is interpolated from the nodes of the leaf that holds the probe (interpolate_volume_potential). At the interior volume nodes (node_errors=True) it comes straight from the FMM, and node_rel_l2 is the relative L2 error there with the mesh’s quadrature weights. The boundary correction is evaluated directly at both. Nodes closer than node_boundary_clearance (1e-3) to the boundary are left out, at most 2% of them: a few nodes of the amr-l5-to8 mesh lie within 2e-5 of the curve, where evaluating the correction gives errors of order one.

Besides the table of mesh sizes and errors, the cell prints:

• error per source point: every case against uniform-l5, and amr-l5-to7 against uniform-l6, which has about as many source points, at the probes and at the nodes;

• where the error sits: the probes and nodes are split into those inside leaves that amr-l5-to7 refined, a band along the boundary, and the rest; for every case the cell prints the RMS error in both groups and the band’s share of the squared probe error.

The first figure plots error against source points for both kinds of refinement, at the probes and at the nodes, and the RMS probe error of the two groups; the second shows the meshes and probe error maps of uniform-l5, amr-l5-to7 and uniform-l6 on one color scale.

amr_common = {
    'nelements_bdry': 768,
    'qbx_order': 6,
    'fmm_order': 18,
    'adaptive_near_boundary_factor': 0.00,
}

amr_compare_plan = [
    {'label': 'uniform-l5', 'vol_nlevels': 5, 'adaptive_boundary_max_level': None},
    {'label': 'amr-l5-to6', 'vol_nlevels': 5, 'adaptive_boundary_max_level': 6},
    {'label': 'amr-l5-to7', 'vol_nlevels': 5, 'adaptive_boundary_max_level': 7},
    {'label': 'amr-l5-to8', 'vol_nlevels': 5, 'adaptive_boundary_max_level': 8},
    {'label': 'uniform-l6', 'vol_nlevels': 6, 'adaptive_boundary_max_level': None},
    {'label': 'uniform-l7', 'vol_nlevels': 7, 'adaptive_boundary_max_level': None},
]

amr_results = []

print('Boundary refinement against uniform refinement (manufactured-solution error):')
print('rel_l2_err, max_abs_err and p95_abs are taken at the probes, node_rel_l2 at the interior volume nodes.')
print('       case  n_cells    n_src    h_min passes   rel_l2_err   max_abs_err    p95_abs  node_rel_l2  gmres')

for cfg in amr_compare_plan:
    try:
        row = run_poisson_case(
            **amr_common,
            vol_nlevels=cfg['vol_nlevels'],
            adaptive_boundary_max_level=cfg['adaptive_boundary_max_level'],
            node_errors=True,
            return_fields=True,
            return_mesh=True,
        )
    except Exception as exc:
        row = {
            'label': cfg['label'],
            'status': 'fail',
            'error': repr(exc),
        }
        amr_results.append(row)
        print(f"{cfg['label']:>11s}  FAIL  {row['error']}")
        continue

    row['label'] = cfg['label']
    row['status'] = 'ok'
    amr_results.append(row)

    gmres_ok = (row['gmres_rhs_ext_state'] == 'success' and row['gmres_corr_state'] == 'success')

    print(
        f"{row['label']:>11s}"
        f" {row['n_cells']:8d}"
        f" {row['n_sources']:8d}"
        f" {row['h_min']:8.5f}"
        f" {row['mesh_refine_passes']:6d}"
        f" {row['rel_l2_err']:12.4e}"
        f" {row['max_abs_err']:12.4e}"
        f" {row['p95_abs_err']:10.3e}"
        f" {row['node_rel_l2_err']:12.4e}"
        f" {'ok' if gmres_ok else 'fail'}"
    )

successful = {row['label']: row for row in amr_results if row.get('status') == 'ok'}
ok_rows = [row for row in amr_results if row.get('status') == 'ok']


def smaller_of(label_a, err_a, label_b, err_b):
    if err_a <= err_b:
        return f'{label_a} {err_b / err_a:.2f}x smaller'
    return f'{label_b} {err_a / err_b:.2f}x smaller'


# Error per source point: every refined case against the mesh it started
# from, then boundary refinement against uniform refinement at about the
# same number of source points; at the probes and at the nodes.
print()
print('Error per source point (rel L2 at the probes, and at the nodes):')
if 'uniform-l5' in successful:
    base_row = successful['uniform-l5']
    for row in ok_rows:
        if row['label'] == 'uniform-l5':
            continue
        err_gain = base_row['rel_l2_err'] / row['rel_l2_err']
        node_err_gain = base_row['node_rel_l2_err'] / row['node_rel_l2_err']
        src_ratio = row['n_sources'] / base_row['n_sources']
        print(
            f"  {row['label']:>11s} vs uniform-l5: {err_gain:6.2f}x smaller at the probes,"
            f" {node_err_gain:6.2f}x at the nodes, for {src_ratio:5.2f}x the source points"
        )
if 'amr-l5-to7' in successful and 'uniform-l6' in successful:
    amr_row = successful['amr-l5-to7']
    uni_row = successful['uniform-l6']
    print(
        f"  at about equal cost ({amr_row['n_sources']} and {uni_row['n_sources']} source points):"
        f" at the probes {smaller_of('amr-l5-to7', amr_row['rel_l2_err'], 'uniform-l6', uni_row['rel_l2_err'])},"
        f" at the nodes {smaller_of('amr-l5-to7', amr_row['node_rel_l2_err'], 'uniform-l6', uni_row['node_rel_l2_err'])}"
    )

# Where the error sits: split the probes into those inside leaves that
# amr-l5-to7 refined (a band along the boundary) and the rest, and look at
# every case's error on the two groups.
in_band = None
if 'amr-l5-to7' in successful:
    band_row = successful['amr-l5-to7']
    leaf_sides = band_row['mesh_leaf_side_lengths_host']
    refined_leaf = leaf_sides < 0.75 * band_row['h']
    refined_centers = band_row['mesh_leaf_centers_host'][refined_leaf]
    refined_half = 0.5 * leaf_sides[refined_leaf]

    def in_refined_leaves(points, chunk_size=1024):
        # compare the points with the refined leaves a chunk at a time, so that
        # the tens of thousands of nodes of the finest meshes do not need a
        # points-by-leaves array at once
        inside = np.empty(len(points), dtype=bool)
        for start in range(0, len(points), chunk_size):
            chunk = points[start:start + chunk_size]
            offset = np.abs(chunk[:, None, :] - refined_centers[None, :, :]).max(axis=2)
            inside[start:start + chunk_size] = np.any(offset <= refined_half[None, :] + 1.0e-12, axis=1)
        return inside

    in_band = in_refined_leaves(probe_points_host)

    print()
    print(
        f'Where the error sits: {np.count_nonzero(in_band)} of {in_band.size} probes'
        ' lie in leaves that amr-l5-to7 refined (the band), the rest in leaves it left alone.'
    )
    print('RMS error at the probes, the band share of their squared error, and RMS error at the nodes:')
    print('       case  rms_err_band  rms_err_rest  band_share  node_rms_band  node_rms_rest')
    for row in ok_rows:
        err = row['abs_err']
        row['rms_err_band'] = float(np.sqrt(np.mean(err[in_band] ** 2)))
        row['rms_err_rest'] = float(np.sqrt(np.mean(err[~in_band] ** 2)))
        row['band_share'] = float(np.sum(err[in_band] ** 2) / np.sum(err**2))
        node_err = row['node_abs_err']
        node_in_band = in_refined_leaves(row['node_points_host'])
        row['node_rms_err_band'] = float(np.sqrt(np.mean(node_err[node_in_band] ** 2)))
        row['node_rms_err_rest'] = float(np.sqrt(np.mean(node_err[~node_in_band] ** 2)))
        print(
            f"{row['label']:>11s}"
            f" {row['rms_err_band']:13.4e}"
            f" {row['rms_err_rest']:13.4e}"
            f" {row['band_share']:11.3f}"
            f" {row['node_rms_err_band']:14.4e}"
            f" {row['node_rms_err_rest']:14.4e}"
        )

if ok_rows:
    fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), constrained_layout=True)

    ax0 = axes[0]
    families = [
        ('uniform refinement', ['uniform-l5', 'uniform-l6', 'uniform-l7'], 'o', 'tab:gray'),
        ('boundary refinement from l5', ['uniform-l5', 'amr-l5-to6', 'amr-l5-to7', 'amr-l5-to8'], 's', 'tab:orange'),
    ]
    for family_label, family, marker, color in families:
        rows = [successful[label] for label in family if label in successful]
        ax0.loglog(
            [row['n_sources'] for row in rows],
            [row['rel_l2_err'] for row in rows],
            marker + '-', lw=2.0, ms=7, color=color, label=f'{family_label}, probes',
        )
        ax0.loglog(
            [row['n_sources'] for row in rows],
            [row['node_rel_l2_err'] for row in rows],
            marker + '--', lw=1.4, ms=5, color=color, mfc='none', label=f'{family_label}, nodes',
        )
    for row in ok_rows:
        ax0.annotate(row['label'], (row['n_sources'], row['rel_l2_err']), textcoords='offset points', xytext=(4, 4), fontsize=8)
    ax0.set_xlabel('number of volume source points')
    ax0.set_ylabel('relative L2 error')
    ax0.set_title('Error against source points')
    ax0.grid(True, which='both', ls=':', alpha=0.5)
    ax0.legend(loc='lower left', fontsize=8)

    ax1 = axes[1]
    if in_band is not None:
        xpos = np.arange(len(ok_rows))
        ax1.bar(xpos - 0.2, [row['rms_err_band'] for row in ok_rows], width=0.4, color='tab:orange', label='probes in the refined band')
        ax1.bar(xpos + 0.2, [row['rms_err_rest'] for row in ok_rows], width=0.4, color='tab:blue', label='all other probes')
        ax1.set_xticks(xpos)
        ax1.set_xticklabels([row['label'] for row in ok_rows], rotation=30, ha='right')
        ax1.set_yscale('log')
        ax1.legend(loc='upper right', fontsize=8)
    ax1.set_ylabel('RMS absolute error')
    ax1.set_title('Error in the refined band and elsewhere')
    ax1.grid(True, axis='y', which='both', ls=':', alpha=0.5)

    plt.show()

map_labels = [label for label in ['uniform-l5', 'amr-l5-to7', 'uniform-l6'] if label in successful]
if map_labels:
    from matplotlib.collections import LineCollection

    def box_segments_from_centers_sizes(centers, sizes):
        cx = np.asarray(centers[:, 0], dtype=np.float64)
        cy = np.asarray(centers[:, 1], dtype=np.float64)
        h = 0.5 * np.asarray(sizes, dtype=np.float64)
        x0 = cx - h
        x1 = cx + h
        y0 = cy - h
        y1 = cy + h

        segs = np.empty((4 * cx.size, 2, 2), dtype=np.float64)
        segs[0::4, 0, 0], segs[0::4, 0, 1] = x0, y0
        segs[0::4, 1, 0], segs[0::4, 1, 1] = x1, y0
        segs[1::4, 0, 0], segs[1::4, 0, 1] = x1, y0
        segs[1::4, 1, 0], segs[1::4, 1, 1] = x1, y1
        segs[2::4, 0, 0], segs[2::4, 0, 1] = x1, y1
        segs[2::4, 1, 0], segs[2::4, 1, 1] = x0, y1
        segs[3::4, 0, 0], segs[3::4, 0, 1] = x0, y1
        segs[3::4, 1, 0], segs[3::4, 1, 1] = x0, y0
        return segs

    map_rows = [successful[label] for label in map_labels]
    log_errs = [np.log10(row['abs_err'] + 1.0e-14) for row in map_rows]
    vmin = min(float(np.min(v)) for v in log_errs)
    vmax = max(float(np.max(v)) for v in log_errs)

    fig, axes = plt.subplots(2, len(map_rows), figsize=(5.2 * len(map_rows), 10.0), constrained_layout=True, squeeze=False)

    for icol, (row, log_err) in enumerate(zip(map_rows, log_errs)):
        ax = axes[0, icol]
        segs = box_segments_from_centers_sizes(
            row['mesh_leaf_centers_host'],
            row['mesh_leaf_side_lengths_host'],
        )
        lc = LineCollection(segs, colors='tab:blue', linewidths=0.22, alpha=0.45)
        ax.add_collection(lc)
        ax.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.3)
        ax.set_xlim(bbox_a, bbox_b)
        ax.set_ylim(bbox_a, bbox_b)
        ax.set_aspect('equal')
        ax.set_title(f"{row['label']}\n{row['n_cells']} cells, {row['n_sources']} source points")

        ax = axes[1, icol]
        sc = ax.scatter(probe_points_host[:, 0], probe_points_host[:, 1], c=log_err, s=10, cmap='magma', vmin=vmin, vmax=vmax)
        ax.plot(starfish_vertices[:, 0], starfish_vertices[:, 1], color='black', lw=1.0)
        ax.set_xlim(bbox_a, bbox_b)
        ax.set_ylim(bbox_a, bbox_b)
        ax.set_aspect('equal')
        ax.set_title(f"{row['label']}: rel L2={row['rel_l2_err']:.2e}")

    cb = fig.colorbar(sc, ax=axes[1, :], shrink=0.84)
    cb.set_label(r'$\log_{10}|u_{num}-u_{exact}|$')
    fig.suptitle('Volume meshes (top) and probe error maps on one color scale (bottom)', fontsize=12)
    plt.show()

What the comparison shows

A run at the default settings on a CPU OpenCL device printed the numbers below; other devices can differ in the last digits. The band holds 1127 of the 2491 probes.

case	cells	source points	rel L2, probes	RMS in the band, probes	RMS elsewhere, probes	rel L2, nodes	RMS elsewhere, nodes
uniform-l5	256	4096	4.82e-05	3.48e-05	1.66e-05	2.80e-05	6.02e-06
amr-l5-to6	454	7264	2.38e-05	9.88e-06	1.52e-05	2.25e-06	8.44e-07
amr-l5-to7	994	15904	2.06e-05	1.47e-06	1.52e-05	5.04e-07	2.36e-07
amr-l5-to8	2254	36064	1.86e-05	1.35e-06	1.37e-05	1.44e-07	1.16e-07
uniform-l6	1024	16384	3.23e-06	2.29e-06	1.19e-06	2.25e-06	8.41e-07
uniform-l7	4096	65536	6.15e-07	4.52e-07	1.98e-07	4.77e-07	1.93e-07

• At the probes. One pass of boundary refinement makes the error 2.0x smaller for 1.8x the source points; further passes add little: 2.3x for 3.9x the source points after two passes, 2.6x for 8.8x after three. On uniform-l5 the band carries 78% of the squared error, and two passes make the band’s RMS error about 24x smaller, but the error elsewhere barely moves (RMS 1.66e-05 on uniform-l5, 1.37e-05 after three passes), so after two passes it is 99% of the squared error. uniform-l6, with about as many source points as amr-l5-to7 (16384 and 15904), has a 6.4x smaller error.

• At the nodes. The error is set by how fine the cut leaves are: amr-l5-to6 and uniform-l6 both reach 2.25e-06, and amr-l5-to7 and uniform-l7 reach 5.0e-07 and 4.8e-07, though the boundary-refined meshes have only 44% and 24% of the source points. At about equal cost amr-l5-to7 is 4.5x more accurate than uniform-l6, and in the error-against-source-points plot every boundary-refined mesh lies below the uniform ones. The error falls even at the nodes of the leaves that boundary refinement leaves at the vol_nlevels=5 size: from 6.0e-06 to 2.4e-07 in two passes.

• What is left at the probes. The probes outside the band lie in those same vol_nlevels=5 leaves. After two passes their error is 64x that of the nodes around them (RMS 1.52e-05 and 2.36e-07), so it comes from interpolating the potential from those nodes to the probes, not from the solution at the nodes. The boundary discretization limits neither measure: with the same one, uniform-l7 reaches 6.2e-07 at the probes and amr-l5-to8 1.4e-07 at the nodes.

So on this problem boundary refinement buys the accuracy of the solution at the nodes: refining only the cut leaves does about what refining every leaf to the same size does, with a fraction of the source points. It does not buy the interpolation to points inside the leaves it leaves alone, and that is why uniform refinement looks better at the probes. Where the solution is needed at such points, the leaves that hold them have to be fine enough for interpolating the potential there as well.