{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7ccaee94",
   "metadata": {},
   "source": [
    "# MHM-USFEM: difficult local reaction–diffusion problems\n",
    "\n",
    "\n",
    "The main workflow is **meshes → spaces → local equations → global balance → assemble → solve → fields and errors**. `MeshHierarchy` associates macro and local meshes; `bind_interface` and `LocalContext` own supported numbering, geometric orientation and native coordinate conversions. The mathematical forms, physical trace meaning, local modes and gauge remain explicit in the cells below. Fully manual/custom spaces use the same numerical owners; see the [custom-interface notebook](https://github.com/ipes-lncc/pymhm/blob/main/notebooks/foundations/operators/custom_interface.ipynb).\n",
    "\n",
    "We explicitly construct `Equation`, `LocalEquations` and `MultiscaleProblem`, then connect each code block to the physical formulation. **MHM-USFEM** is the name used here for scalar **MHM-UNUSUAL** of [Martins, Pacazuca and Valentin (2025), Eqs. (14)–(15), §4.1](https://doi.org/10.55592/cilamce2025.v5i.14270). This is a scalar reaction–diffusion example; the separate Stokes–Brinkman tutorial treats velocity and pressure. Advection is zero here, and the negative residual stabilization is different from SUPG.\n",
    "\n",
    "Open the downloaded notebook with `jupyter lab` after installing `pymhm[notebooks,visualization]` and the compatible native DOLFINx/UFL backend. Exact data, sources, local and global forms, boundary conditions, independent references, norms and plots are all defined below. PyMHM supplies generic UFL compilation, oriented trace integration, condensation and checked linear algebra.\n",
    "\n",
    "We compare an intentionally underresolved P1/P0 control with a refined P1/P0 family that partitions each macroface and satisfies the paper's red-refinement condition. Both use the same physical equation and the published stabilization parameter. Reduced nodal oscillation can coexist with a larger integrated error; we report both. Our declared SW–NE connectivity, strong exterior Dirichlet enforcement and local meshes do not constitute a reproduction of the paper's historical figures.\n",
    "\n",
    "\n",
    "The PyMHM distribution contains only the library. The first cell explicitly\n",
    "downloads a checksum-verified companion archive and prepares the declared data\n",
    "using its local Python helpers. You can inspect these support files in `ROOT`.\n",
    "Complete studies and historical replay retain their existing opt-in flags.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ee88c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import os\n",
    "import sys\n",
    "from pymhm.io.workspace import workspace_from_archive\n",
    "\n",
    "# Download verified support files; this operation does not execute them.\n",
    "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/e482bd3b768df30dca25cc72b5e147b051f60b59f6d8b7c0d87bc976b931b487/mhm_usfem_rad-companion.zip\"\n",
    "COMPANION_SHA256 = \"e482bd3b768df30dca25cc72b5e147b051f60b59f6d8b7c0d87bc976b931b487\"\n",
    "WORKSPACE = Path(\n",
    "    os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd() / \".pymhm-companions\" / COMPANION_SHA256)\n",
    ")\n",
    "ROOT = workspace_from_archive(COMPANION_URL, sha256=COMPANION_SHA256, directory=WORKSPACE)\n",
    "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n",
    "if str(ROOT) not in sys.path:\n",
    "    sys.path.insert(0, str(ROOT))\n",
    "\n",
    "# Explicitly prepare declared data with the downloaded Python helpers.\n",
    "from scripts.notebook_reproduction import notebook_workspace\n",
    "\n",
    "ROOT = notebook_workspace(\"introduction/mhm_usfem_rad.ipynb\", directory=ROOT)\n",
    "root = ROOT\n",
    "print(\"Workspace:\", ROOT)\n",
    "\n",
    "\n",
    "import json\n",
    "import numpy as np\n",
    "from pymhm import TriangleMesh, FaceSpace, SkeletonSpace\n",
    "from pymhm.core.equations import Equation, LocalEquations\n",
    "from pymhm.core.multiscale import assemble\n",
    "\n",
    "REPORTS = ROOT / \"build/introduction\"\n",
    "REPORTS.mkdir(parents=True, exist_ok=True)\n",
    "np.set_printoptions(precision=5, suppress=True)\n",
    "from pymhm.fem.scalar.triangle import nodal_space\n",
    "\n",
    "from dataclasses import dataclass, replace\n",
    "from typing import Any\n",
    "from threadpoolctl import threadpool_limits\n",
    "\n",
    "threadpool_limits(limits=1)\n",
    "from pymhm import MeshHierarchy, LocalContext, bind_interface, bind_problem\n",
    "\n",
    "from functools import partial\n",
    "from examples.introduction.transport import (\n",
    "    scalar_error_norms,\n",
    "    scalar_reference,\n",
    "    record_rad_case,\n",
    "    plot_rad_primary,\n",
    "    plot_rad_severe,\n",
    "    plot_rad_refined,\n",
    "    plot_rad_profiles,\n",
    ")\n",
    "from examples.introduction.vector import execution_provenance\n",
    "\n",
    "measure_case = partial(record_rad_case, reports=REPORTS, root=ROOT)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rad-runtime-policy",
   "metadata": {},
   "source": [
    "### Runtime resources (independent of the formulation)\n",
    "\n",
    "The equations, approximation spaces, quadrature and error checks below do not depend on the worker count. We schedule independent local solves with `ExecutionConfig`; the coordinator adds their global contributions in a deterministic order, including contributions to shared faces. The conforming reference solves and postprocessing remain serial.\n",
    "\n",
    "The default `PYMHM_RAD_EXECUTION=auto` runs small cases serially. It selects processes for the resolved family when there are at least 512 macroelements, eight local subdivisions and four face segments. This choice follows a representative timing control for this tutorial, rather than a general performance claim. The default upper bound is 32 workers, limited by this process's CPU affinity. Set `PYMHM_RAD_WORKERS` to override that upper bound, or `PYMHM_RAD_EXECUTION=serial`, `thread` or `process` to explicitly select a backend. Passing an `ExecutionConfig` directly to `solve_rad_case` bypasses the automatic choice.\n",
    "\n",
    "Each worker uses one native BLAS/OpenMP thread. A bounded pipeline holds at most twice the active worker count in flight, and each solve caps the active worker count at its number of macroelements. These options change execution resources only.\n",
    "\n",
    "A notebook-owned temporary directory holds the local UFL compiler cache. For every solve configuration, we assemble the first two macrotriangles serially before starting workers: these triangles cover the two orientations in this mesh. This warms the volume and face kernels for the chosen diffusion, stabilization and face partition. The optional process worker module below is generated in this same owned directory from the notebook's functions. The cache remains alive throughout the study and is removed after the final result is saved. The 7200-second compiler timeout is not a solver tolerance.\n",
    "\n",
    "After the final cleanup, rerun both the runtime-resource and worker-export cells before solving again. If the study is interrupted, let PyMHM finish releasing its workers, then call `RAD_JIT_CACHE.cleanup()` or restart the kernel before restarting the setup.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rad-runtime-configuration",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tempfile import TemporaryDirectory\n",
    "from pymhm.execution.cpu import ExecutionConfig\n",
    "\n",
    "try:\n",
    "    RAD_AVAILABLE_CPUS = len(os.sched_getaffinity(0))\n",
    "except (AttributeError, OSError):\n",
    "    RAD_AVAILABLE_CPUS = os.cpu_count() or 1\n",
    "\n",
    "RAD_POLICY = os.environ.get(\"PYMHM_RAD_EXECUTION\", \"auto\")\n",
    "if RAD_POLICY not in {\"auto\", \"serial\", \"thread\", \"process\"}:\n",
    "    raise ValueError(\"PYMHM_RAD_EXECUTION must be 'auto', 'serial', 'thread', or 'process'\")\n",
    "RAD_WORKERS = int(os.environ.get(\"PYMHM_RAD_WORKERS\", str(min(32, RAD_AVAILABLE_CPUS))))\n",
    "RAD_EXECUTION = ExecutionConfig(\n",
    "    backend=\"process\" if RAD_POLICY == \"auto\" else RAD_POLICY,\n",
    "    workers=RAD_WORKERS,\n",
    "    native_threads=1,\n",
    "    batch_size=2 * RAD_WORKERS,\n",
    "    pipeline=True,\n",
    ")\n",
    "RAD_JIT_CACHE = TemporaryDirectory(prefix=\"pymhm-rad-jit-\")\n",
    "RAD_JIT_OPTIONS = {\"cache_dir\": Path(RAD_JIT_CACHE.name), \"timeout\": 7200}\n",
    "print(\n",
    "    \"Local execution policy:\", RAD_POLICY,\n",
    "    \"workers:\", RAD_EXECUTION.effective_workers,\n",
    "    \"native threads per worker:\", RAD_EXECUTION.native_threads,\n",
    "    \"pending bound:\", RAD_EXECUTION.effective_batch_size,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d652fb1",
   "metadata": {},
   "source": [
    "## 1. Identify the equation, boundary data and difficult length scale\n",
    "\n",
    "On the unit square, reaction and source are one and diffusion is a positive constant:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "-\\epsilon\\Delta u+u&=1, & q&=-\\epsilon\\nabla u,\\\\\n",
    "u&=0\\quad\\text{on }x=0,1, & q\\cdot n&=0\\quad\\text{on }y=0,1.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The exact solution has two layers with width of order $\\sqrt\\epsilon$:\n",
    "\n",
    "$$\n",
    "u(x,y)=1-\\frac{e^{-x/\\sqrt\\epsilon}+e^{-(1-x)/\\sqrt\\epsilon}}{1+e^{-1/\\sqrt\\epsilon}}.\n",
    "$$\n",
    "\n",
    "Differentiating gives $-\\epsilon u''+u=1$, independently of the discrete operator. At $\\epsilon=10^{-5}$, a coarse uniform local mesh cannot resolve the layer. We measure nodal overshoot above the physical upper bound one, nodal error and integrated errors. The coarse control keeps the same spaces as diffusion decreases; a separate resolution family tests whether the physical fields improve when both the local and trace approximation spaces are refined.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "692777ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass(frozen=True)\n",
    "class ReactionLayer:\n",
    "    \"\"\"Equation (22): ``-epsilon*Delta u+u=1`` with vertical zero Dirichlet data.\n",
    "\n",
    "    Horizontal boundaries carry zero outward physical flux ``-epsilon*grad(u).n``.\n",
    "    ``epsilon`` is positive and the layer width is its square root. Evaluation\n",
    "    avoids positive exponential arguments, including in very thin layers.\n",
    "    \"\"\"\n",
    "\n",
    "    epsilon: float\n",
    "\n",
    "    def __post_init__(self) -> None:\n",
    "        \"\"\"Require a finite positive physical diffusion coefficient.\"\"\"\n",
    "        if not np.isfinite(self.epsilon) or self.epsilon <= 0:\n",
    "            raise ValueError(\"epsilon must be finite and positive\")\n",
    "\n",
    "    def value(self, points: np.ndarray) -> np.ndarray:\n",
    "        \"\"\"Evaluate the exact scalar on physical unit-square coordinates.\"\"\"\n",
    "        scale = 1 / np.sqrt(self.epsilon)\n",
    "        x = points[:, 0]\n",
    "        return 1 - (np.exp(-scale * x) + np.exp(-scale * (1 - x))) / (1 + np.exp(-scale))\n",
    "\n",
    "    def gradient(self, points: np.ndarray) -> np.ndarray:\n",
    "        \"\"\"Evaluate the independently differentiated physical gradient.\"\"\"\n",
    "        scale = 1 / np.sqrt(self.epsilon)\n",
    "        x = points[:, 0]\n",
    "        derivative = scale * (np.exp(-scale * x) - np.exp(-scale * (1 - x))) / (1 + np.exp(-scale))\n",
    "        return np.column_stack((derivative, np.zeros(len(points))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22078d7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "EPSILON = 1e-3\n",
    "# This scientific family fixes scalar P1 and multiplier P0.\n",
    "LOCAL_SUBDIVISIONS = 2\n",
    "LOCAL_QUADRATURE_DEGREE, OPERATOR_CHECK_ORDER, ERROR_ORDER = 8, 12, 24\n",
    "truth = ReactionLayer(EPSILON)\n",
    "\n",
    "# The independently differentiated equation is -epsilon*u'' + u = 1.\n",
    "x = np.linspace(0, 1, 101)\n",
    "points = np.column_stack((x, np.full_like(x, 0.37)))\n",
    "scale = 1 / np.sqrt(EPSILON)\n",
    "second = -(scale**2) * (np.exp(-scale * x) + np.exp(-scale * (1 - x))) / (1 + np.exp(-scale))\n",
    "np.testing.assert_allclose(-EPSILON * second + truth.value(points), 1, atol=2e-15)\n",
    "np.testing.assert_allclose(truth.gradient(points)[:, 1], 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98657b19",
   "metadata": {},
   "source": [
    "## 2. Write the local equations and skeletal balance\n",
    "\n",
    "Let $B_K$ integrate the globally oriented flux basis against the trace of a local test function. In the matrix equations below, $u_K$ is the local scalar coefficient vector and $\\lambda$ is the global trace coefficient vector; their basis evaluations give the physical scalar and flux multiplier. The API expresses\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "A_Ku_K+B_K\\lambda&=F_K,\\\\\n",
    "\\sum_K C_Ku_K&=0, & C_K&=-B_K^\\mathsf T.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Summing the $C_Ku_K$ equations imposes continuity of scalar trace moments on internal faces. The normal interface binding supplies the two opposite incident orientations; `b` and `c` are declared independently with UFL boundary forms. On active flux faces, $\\lambda=q\\cdot n_F$ uses a single global face normal; with zero advection there is no Robin correction to the physical diffusive flux.\n",
    "\n",
    "The complete strong reaction–diffusion operator is $Lw=w-\\epsilon\\Delta w$. USFEM changes the local bilinear and linear forms:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_K(u,v)&=(\\epsilon\\nabla u,\\nabla v)_K+(u,v)_K\n",
    " -\\sum_{T\\subset K}\\tau_T(Lu,Lv)_T,\\\\\n",
    "\\ell_K(v)&=(1,v)_K-\\sum_{T\\subset K}\\tau_T(1,Lv)_T,\\\\\n",
    "\\tau_T&=\\frac{m h_T^2}{\\max\\{m h_T^2,2\\epsilon\\}+2\\epsilon},\\qquad m=\\tfrac13.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Locals are P1, hence $\\Delta v_h=0$ inside each fine triangle and $Lv_h=v_h$. This justifies the published $m=1/3$ for this space. Below we declare the Galerkin form directly in UFL, then add the **negative** residual form and its matching source. `assemble_form` executes native assembly; the native space binding owns the coefficient convention, while context trace pairings supply the declared UFL interface integrals.\n",
    "\n",
    "Reaction makes the local operator invertible: constants are not a kernel, and no retained coarse amplitudes are required. On vertical exterior faces, zero scalar data are imposed strongly by removing their nodal coordinates and exterior flux columns. Horizontal exterior faces prescribe zero normal flux. Vertical Dirichlet faces are absent from flux coupling: their zero dummy trace coordinates do not prescribe the physical flux. The global balance is tested on the free internal trace coordinates.\n",
    "\n",
    "For P1/P0 in two dimensions, Theorem 1 requires one red refinement of a simplicial local mesh that matches the macroface partition. With one P0 segment per macroface, two local subdivisions meet this condition. With four P0 segments per face, we declare a matching uniform four-subdivision mesh and take its red refinement: eight local subdivisions. This condition concerns compatibility of the trace and local spaces; a positive local matrix alone does not establish global stability.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f68970c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import basix.ufl\n",
    "import dolfinx\n",
    "import ufl\n",
    "\n",
    "\n",
    "def rad_forms(domain: Any, *, epsilon: float, stabilized: bool) -> tuple[Any, Any, Any]:\n",
    "    \"\"\"Translate the local reaction–diffusion bilinear form and load into UFL.\"\"\"\n",
    "    V = dolfinx.fem.functionspace(domain, basix.ufl.element(\"Lagrange\", \"triangle\", 1))\n",
    "    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)\n",
    "    dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": LOCAL_QUADRATURE_DEGREE})\n",
    "    a = (epsilon * ufl.inner(ufl.grad(u), ufl.grad(v)) + u * v) * dx\n",
    "    L = v * dx\n",
    "    if stabilized:\n",
    "        Lu, Lv = u - epsilon * ufl.div(ufl.grad(u)), v - epsilon * ufl.div(ufl.grad(v))\n",
    "        h = ufl.CellDiameter(domain)\n",
    "        tau = (h**2 / 3) / (ufl.max_value(h**2 / 3, 2 * epsilon) + 2 * epsilon)\n",
    "        a -= tau * Lu * Lv * dx\n",
    "        L -= tau * Lv * dx\n",
    "    return V, a, L\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c665a8c",
   "metadata": {},
   "source": [
    "The native mesh and coefficient conversions come from the shared backend binding. The provider below writes the volume and boundary forms, chooses its physical boundary conditions and registers the scalar field. The explicit restriction to free scalar coordinates imposes the declared vertical Dirichlet data; it is part of the formulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "940ae1a0",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import basix\n",
    "\n",
    "\n",
    "from pymhm.backends.spaces import create_native_mesh\n",
    "from pymhm.backends.forms import assemble_form, assemble_pairing\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf104b08",
   "metadata": {},
   "source": [
    "`LocalEquations` now receives the executed forms through the shared native-space binding. `LocalContext` owns the mesh and interface maps. The following provider applies the vertical Dirichlet nodes, keeps the horizontal physical-flux convention, and exposes the two oriented trace blocks directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acca5234",
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_rad(\n",
    "    local: LocalContext,\n",
    "    *,\n",
    "    macro: TriangleMesh,\n",
    "    skeleton: SkeletonSpace,\n",
    "    epsilon: float,\n",
    "    stabilized: bool,\n",
    "    local_subdivisions: int = LOCAL_SUBDIVISIONS,\n",
    ") -> LocalEquations:\n",
    "    \"\"\"Declare and execute the local scalar UFL form, then apply explicit mixed data.\"\"\"\n",
    "    fine = local.mesh\n",
    "    domain = create_native_mesh(fine)\n",
    "    V, a, L = rad_forms(domain, epsilon=epsilon, stabilized=stabilized)\n",
    "    # Generic native assembly and coordinate permutation, independent of the PDE.\n",
    "    binding = local.native_space(V)\n",
    "    order = binding.mapping\n",
    "    native_A = assemble_form(a, jit_options=RAD_JIT_OPTIONS)\n",
    "    native_F = assemble_form(L, jit_options=RAD_JIT_OPTIONS)\n",
    "    A, F = native_A[order][:, order], native_F[order]\n",
    "    _, nodes = nodal_space(fine, 1)\n",
    "    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)\n",
    "    pair_b = local.trace_pairings(lambda phi, ds: phi * v * ds)\n",
    "    pair_c = local.trace_pairings(lambda phi, ds: -phi * u * ds, axis=\"rows\")\n",
    "    B = assemble_pairing(\n",
    "        pair_b.forms, axis=\"columns\", jit_options=RAD_JIT_OPTIONS\n",
    "    )[order]\n",
    "    C = assemble_pairing(\n",
    "        pair_c.forms, axis=\"rows\", jit_options=RAD_JIT_OPTIONS\n",
    "    )[:, order]\n",
    "    fixed_nodes = np.flatnonzero(\n",
    "        np.isclose(nodes[:, 0], 0, atol=1e-12) | np.isclose(nodes[:, 0], 1, atol=1e-12)\n",
    "    )\n",
    "    free = np.setdiff1d(np.arange(len(nodes)), fixed_nodes)\n",
    "    local.field(\"scalar\", binding, reconstruction=binding.to_native(np.eye(len(nodes))[:, free]))\n",
    "    B, C = B[free], C[:, free]\n",
    "    return local.equations(\n",
    "        a=A[free][:, free],\n",
    "        L=F[free],\n",
    "        b=B,\n",
    "        c=C,\n",
    "        metadata={\n",
    "            \"mesh\": fine,\n",
    "            \"free\": free,\n",
    "            \"nodes\": len(nodes),\n",
    "            \"local_subdivisions\": local_subdivisions,\n",
    "        },\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "rad-spawn-runtime-explanation",
   "metadata": {},
   "source": [
    "### Optional runtime setup: importable functions for spawned workers\n",
    "\n",
    "The mathematical definitions above are the complete local formulation. The next cell concerns Python process execution only: a spawned worker must import its callable from a module, whereas a notebook function lives in an interactive session.\n",
    "\n",
    "We use `inspect.getsource` to export **the exact functions just defined above**, their imports and the declared constants to a temporary module. There is no second maintained copy of the operator or stabilization. PyMHM still owns process scheduling, local elimination, ordered assembly and native resource cleanup. You can select serial or thread execution with the same mathematical functions; those modes use the original notebook functions directly.\n",
    "\n",
    "The module name identifies the complete generated source payload. If you edit a function or constant and rerun this cell, workers import the new implementation rather than a previously imported notebook callback.\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "rad-spawn-runtime-module",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import hashlib\n",
    "import importlib\n",
    "import inspect\n",
    "\n",
    "RAD_PROVIDER_SOURCE = \"\\n\\n\".join(\n",
    "    inspect.getsource(function).rstrip() for function in (rad_forms, local_rad)\n",
    ")\n",
    "RAD_WORKER_IMPORTS = \"\"\"from pathlib import Path\n",
    "from typing import Any\n",
    "import numpy as np\n",
    "import basix.ufl\n",
    "import dolfinx\n",
    "import ufl\n",
    "from pymhm import TriangleMesh, SkeletonSpace, LocalContext\n",
    "from pymhm.core.equations import LocalEquations\n",
    "from pymhm.fem.scalar.triangle import nodal_space\n",
    "from pymhm.backends.spaces import create_native_mesh\n",
    "from pymhm.backends.forms import assemble_form, assemble_pairing\n",
    "\"\"\"\n",
    "RAD_WORKER_CONSTANTS = (\n",
    "    f\"LOCAL_SUBDIVISIONS = {LOCAL_SUBDIVISIONS!r}\\n\"\n",
    "    f\"LOCAL_QUADRATURE_DEGREE = {LOCAL_QUADRATURE_DEGREE!r}\\n\"\n",
    "    f\"RAD_JIT_OPTIONS = {{'cache_dir': Path({RAD_JIT_CACHE.name!r}), 'timeout': 7200}}\\n\"\n",
    ")\n",
    "RAD_WORKER_PAYLOAD = (\n",
    "    RAD_WORKER_IMPORTS + \"\\n\" + RAD_WORKER_CONSTANTS + \"\\n\" + RAD_PROVIDER_SOURCE + \"\\n\"\n",
    ")\n",
    "RAD_WORKER_MODULE_NAME = (\n",
    "    \"_pymhm_rad_workers_\" + hashlib.sha256(RAD_WORKER_PAYLOAD.encode()).hexdigest()\n",
    ")\n",
    "RAD_WORKER_FILE = Path(RAD_JIT_CACHE.name) / (RAD_WORKER_MODULE_NAME + \".py\")\n",
    "RAD_WORKER_FILE.write_text(RAD_WORKER_PAYLOAD, encoding=\"utf-8\")\n",
    "sys.path.insert(0, RAD_JIT_CACHE.name)\n",
    "importlib.invalidate_caches()\n",
    "RAD_WORKER_MODULE = importlib.import_module(RAD_WORKER_MODULE_NAME)\n",
    "RAD_PROVIDER_SOURCE_SHA256 = hashlib.sha256(RAD_PROVIDER_SOURCE.encode()).hexdigest()\n",
    "print(\"Spawned workers import the notebook's exact provider source:\", RAD_PROVIDER_SOURCE_SHA256)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7e8f9cf",
   "metadata": {},
   "source": [
    "The primary operator above is the executed UFL form. As a secondary convenience check, `scalar_operators` can assemble its unstabilized part. We compare its restricted matrix and load against our declared Galerkin equation on one representative macrocell; it is not used to define the tutorial's MHM operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c9ed264",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymhm.fem.scalar.triangle import scalar_operators\n",
    "\n",
    "probe_macro = TriangleMesh.unit_square(2)\n",
    "probe_skeleton = SkeletonSpace(probe_macro, tuple(FaceSpace.uniform(0) for _ in probe_macro.faces))\n",
    "probe_problem = bind_problem(\n",
    "    MeshHierarchy(\n",
    "        probe_macro,\n",
    "        tuple(\n",
    "            probe_macro.submesh(cell, LOCAL_SUBDIVISIONS) for cell in range(len(probe_macro.cells))\n",
    "        ),\n",
    "    ),\n",
    "    bind_interface(probe_skeleton, convention=\"normal\"),\n",
    "    lambda local: local_rad(\n",
    "        local, macro=probe_macro, skeleton=probe_skeleton, epsilon=EPSILON, stabilized=False\n",
    "    ),\n",
    ")\n",
    "probe_equation = probe_problem.local_provider(0)\n",
    "ready_A, _, ready_F = scalar_operators(\n",
    "    probe_equation.metadata[\"mesh\"],\n",
    "    1,\n",
    "    diffusion=EPSILON,\n",
    "    reaction=1.0,\n",
    "    source=1.0,\n",
    "    order=OPERATOR_CHECK_ORDER,\n",
    ")\n",
    "free = probe_equation.metadata[\"free\"]\n",
    "np.testing.assert_allclose(probe_equation.a.toarray(), ready_A[free][:, free].toarray(), atol=1e-14)\n",
    "np.testing.assert_allclose(probe_equation.L, ready_F[free], atol=1e-14)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46af816",
   "metadata": {},
   "source": [
    "## 3. Declare and solve the global problem\n",
    "\n",
    "`Equation(0,0)` adds zero to the sum of the local balances. Homogeneous boundary data already enter the local restrictions and fixed exterior trace coefficients. `assemble` eliminates each independent local operator; the global solve then reconstructs the local coefficients. No call selects a PDE from a model name."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05a5308d",
   "metadata": {},
   "source": [
    "We measure the broken scalar error and physical flux $-\\epsilon\\nabla u_h$ against the analytical fields. The following quadrature utility also records nodal extrema, since a small integral error alone can miss oscillations in a thin layer.\n",
    "Norm integration, archive replay and field/profile displays are importable from `examples/introduction/transport.py` and `examples/introduction/vector.py`. The local forms and the complete global solve remain above; the conforming P2 forms below are separately declared.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33246e0c",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Physical scalar/gradient/flux norms are integrated by scalar_error_norms (imported above).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11700a98",
   "metadata": {},
   "source": [
    "### Record the numerical coordinates for reproducibility\n",
    "\n",
    "The next utility saves the actual local coordinate basis, coefficients, orientation maps and reconstructed fields as each solve finishes. It contains no physical operator. Replaying the saved basis is checked with one and two BLAS threads; this positive-reaction example has no local nullspace modes. You can read the global construction immediately below without studying the archive format first.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e18e867c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_rad_case(\n",
    "    n,\n",
    "    epsilon,\n",
    "    stabilized,\n",
    "    *,\n",
    "    face_subdivisions=1,\n",
    "    local_subdivisions=LOCAL_SUBDIVISIONS,\n",
    "    execution: ExecutionConfig | None = None,\n",
    "):\n",
    "    \"\"\"Declare the mesh, multiplier space and scalar global equation for one configuration.\"\"\"\n",
    "    macro = TriangleMesh.unit_square(n)\n",
    "    skeleton = SkeletonSpace(\n",
    "        macro, tuple(FaceSpace.uniform(0, face_subdivisions) for _ in macro.faces)\n",
    "    )\n",
    "    if execution is None:\n",
    "        execution = RAD_EXECUTION\n",
    "        resolved_parallel_case = (\n",
    "            len(macro.cells) >= 512 and local_subdivisions >= 8 and face_subdivisions >= 4\n",
    "        )\n",
    "        if RAD_POLICY == \"auto\" and not resolved_parallel_case:\n",
    "            execution = replace(execution, backend=\"serial\")\n",
    "    local_workers = min(execution.effective_workers, len(macro.cells))\n",
    "    case_execution = replace(\n",
    "        execution,\n",
    "        workers=local_workers,\n",
    "        batch_size=min(execution.effective_batch_size, 2 * local_workers),\n",
    "    )\n",
    "    local_function = (\n",
    "        RAD_WORKER_MODULE.local_rad if execution.backend == \"process\" else local_rad\n",
    "    )\n",
    "    provider = partial(\n",
    "        local_function,\n",
    "        macro=macro,\n",
    "        skeleton=skeleton,\n",
    "        epsilon=epsilon,\n",
    "        stabilized=stabilized,\n",
    "        local_subdivisions=local_subdivisions,\n",
    "    )\n",
    "    hierarchy = MeshHierarchy(\n",
    "        macro, tuple(macro.submesh(c, local_subdivisions) for c in range(len(macro.cells)))\n",
    "    )\n",
    "    # Vertical scalar data are strong in local_rad; horizontal physical flux is zero.\n",
    "    fixed = {int(d): 0.0 for f in macro.boundary_faces for d in skeleton.dofs(int(f))}\n",
    "    problem = bind_problem(\n",
    "        hierarchy,\n",
    "        bind_interface(skeleton, convention=\"normal\"),\n",
    "        provider,\n",
    "        global_equation=Equation(0, 0),\n",
    "        retained=0,\n",
    "        fixed=fixed,\n",
    "    )\n",
    "    # Compile both triangle orientations on the coordinator before starting workers.\n",
    "    # This covers this coefficient/stabilization/face-partition configuration.\n",
    "    for cell in range(min(2, len(macro.cells))):\n",
    "        problem.local_provider(cell)\n",
    "    system = assemble(problem, execution=case_execution)\n",
    "    return macro, skeleton, system, system.solve()\n",
    "\n",
    "\n",
    "METHODS = ((\"MHM-Galerkin\", False), (\"MHM-USFEM\", True))\n",
    "cases, rows, state_archives = {}, [], []\n",
    "for n in (2, 4, 8, 16):\n",
    "    for method, stabilized in METHODS:\n",
    "        case, row, archive = measure_case(\n",
    "            *solve_rad_case(n, EPSILON, stabilized),\n",
    "            truth,\n",
    "            order=ERROR_ORDER,\n",
    "            name=f\"rad-n{n}-{method.lower().replace(' ', '-')}\",\n",
    "            named=True,\n",
    "            method=method,\n",
    "            n=n,\n",
    "            H=1 / n,\n",
    "        )\n",
    "        cases[n, method] = case\n",
    "        rows.append(row)\n",
    "        state_archives.append(archive)\n",
    "rows\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9260da5",
   "metadata": {},
   "source": [
    "## 4. Independently assemble and refine a classical conforming reference\n",
    "\n",
    "The CG2 reference uses the same diffusion, unit reaction/source and mixed boundary data, with a globally conforming field and no skeleton or stabilization. Its real UFL form is independently assembled by DOLFINx below. We impose zero scalar values strongly on the vertical boundary and leave the horizontal zero-flux condition natural.\n",
    "\n",
    "The exact solution lets us verify true reference errors on several fine meshes. For the severe $\\epsilon=10^{-5}$ challenge, additional meshes have 256, 512 and 1024 divisions in the layer-normal direction and 16 in the invariant direction. This declared anisotropy resolves the physical scale without changing the operator or boundary data. A numerical reference remains a numerical field, even when an exact field is also available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fcb1bf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def classical_rad_forms(space, epsilon):\n",
    "    \"\"\"Declare the separate conforming P2 operator, source and exact field.\"\"\"\n",
    "    u, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n",
    "    domain = space.ufl_domain()\n",
    "    dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 20})\n",
    "    a = (epsilon * ufl.inner(ufl.grad(u), ufl.grad(v)) + u * v) * dx\n",
    "    x = ufl.SpatialCoordinate(domain)\n",
    "    exact = 1 - (ufl.exp(-x[0] / np.sqrt(epsilon)) + ufl.exp(-(1 - x[0]) / np.sqrt(epsilon))) / (\n",
    "        1 + np.exp(-1 / np.sqrt(epsilon))\n",
    "    )\n",
    "    return a, v * dx, exact, dx\n",
    "\n",
    "\n",
    "reference_rows, references = [], {}\n",
    "REFERENCE_GRIDS = [(EPSILON, n, n) for n in (32, 64, 128)] + [\n",
    "    (1e-5, n, 16) for n in (256, 512, 1024)\n",
    "]\n",
    "for epsilon, nx, ny in REFERENCE_GRIDS:\n",
    "    case, row, archive = scalar_reference(\n",
    "        epsilon, nx, ny, classical_rad_forms, ReactionLayer(epsilon), reports=REPORTS, root=ROOT\n",
    "    )\n",
    "    references[epsilon, nx] = case\n",
    "    reference_rows.append(row)\n",
    "    state_archives.append(archive)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7350bffe",
   "metadata": {},
   "source": [
    "## 5. Inspect fields and observed convergence rates\n",
    "\n",
    "Every spatial panel highlights the **actual MHM macro mesh**, including analytical and reference panels. Flux means the raw physical field $-\\epsilon\\nabla u_h$; the skeletal solve does not make it an H(div) reconstruction.\n",
    "\n",
    "The sequence keeps P1 locals, P0 faces and two local subdivisions per macro edge, varying only the macro mesh. Reported rates are measured from this layer problem; smooth-problem orders are not assumed. Reference errors are measured separately against the exact solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d46b1919",
   "metadata": {},
   "source": [
    "For two successive resolutions, the measured slope is $r=\\log(e_1/e_2)/\\log(H_1/H_2)$. These small plotting utilities compute that quantity directly and keep error and rate axes separate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fb6ce7b",
   "metadata": {},
   "source": [
    "The display utility evaluates each local polynomial independently, including its one-sided boundary values. Every panel overlays the actual macro mesh and has its own color scale. It changes only the display sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44d44d80",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_rad_primary(cases, rows, reference_rows, references, truth, reports=REPORTS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "279a1c96",
   "metadata": {},
   "source": [
    "## 6. Preserve a deliberately underresolved control\n",
    "\n",
    "Fix the macro mesh at $n=8$, with one P0 segment per macroface and two local subdivisions. Both methods use the same scalar PDE and approximation spaces. Decreasing $\\epsilon$ makes the physical layer thinner without increasing resolution. This control demonstrates the approximation limits of these spaces.\n",
    "\n",
    "For constant diffusion and P1 locals, the elementwise Laplacian vanishes. On the uniform fine mesh, write $\\eta=1-\\tau$. The stabilized matrix and source reduce exactly to\n",
    "\n",
    "$$\n",
    "A_{\\mathrm{US}}=\\epsilon K+\\eta M,\n",
    "\\qquad F_{\\mathrm{US}}=\\eta F.\n",
    "$$\n",
    "\n",
    "Here $K$ and $M$ are the ordinary P1 stiffness and mass matrices. In the reaction-dominated branch, dividing by $\\eta$ exposes the discrete diffusion scale\n",
    "\n",
    "$$\n",
    "\\epsilon_{\\mathrm{effective}}\n",
    "=\\frac{\\epsilon}{1-\\tau}\n",
    "=\\epsilon+\\frac{h_T^2}{6}.\n",
    "$$\n",
    "\n",
    "This identity explains why the underresolved stabilized profile can have a broader transition. It does not change the physical diffusion or justify changing the method parameter. For this control, $h_T=\\sqrt{2}/16\\approx0.0884$, whereas the physical layer width at $\\epsilon=10^{-5}$ is $\\sqrt{\\epsilon}\\approx0.00316$.\n",
    "\n",
    "The tables keep overshoot, nodal error, scalar L2 and physical-flux L2 visible. USFEM need not improve every measure or preserve nodal bounds. The profiles retain independent one-sided values at every macroface intersection.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49601836",
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep, challenging = [], {}\n",
    "macro = TriangleMesh.unit_square(8)\n",
    "for epsilon in dict.fromkeys((1e-2, EPSILON, 1e-5)):\n",
    "    for method, stabilized in METHODS:\n",
    "        case, row, archive = measure_case(\n",
    "            *solve_rad_case(8, epsilon, stabilized),\n",
    "            ReactionLayer(epsilon),\n",
    "            order=40,\n",
    "            name=f\"rad-eps{epsilon:g}-{method.lower()}\",\n",
    "            epsilon=epsilon,\n",
    "            method=method,\n",
    "        )\n",
    "        row[\"overshoot\"] = max(0.0, row[\"nodal_maximum\"] - 1)\n",
    "        challenging[epsilon, method] = case[2:]\n",
    "        sweep.append(row)\n",
    "        state_archives.append(archive)\n",
    "sweep\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c11b064",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_rad_severe(macro, challenging, references, sweep, ReactionLayer, reports=REPORTS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3c89bac",
   "metadata": {},
   "source": [
    "## 7. Refine the local and trace spaces together\n",
    "\n",
    "Keep P1 local polynomials and P0 on **each face segment**, but use four equal segments per macroface and eight local subdivisions. The eight-subdivision mesh is one red refinement of the matching four-subdivision mesh. We vary macro divisions through $n=4,8,16$ for $\\epsilon=10^{-5}$ and also compute $\\epsilon=10^{-3}$ at $n=16$ for the profile comparison.\n",
    "\n",
    "The finest case has 512 macrotriangles and four P0 segments per face, as reported for the paper's analytical comparison. The paper states $\\bar H=H/4$ and $h=\\bar H$, but does not specify the historical matching triangulation and connectivity sufficiently to recover the local meshes. Our uniform red-refined realization is declared here. Moreover, the paper's continuous hybrid formulation imposes Dirichlet trace moments weakly; this tutorial imposes the exterior data strongly on local nodes. These are material differences when interpreting the discrete results. We compare the same analytical PDE and an admissible discretization; we do not claim a matched reproduction of Fig. 2.\n",
    "\n",
    "The global equation is still the sum of local scalar-trace balances. Only the geometry and the number of face segments change below; the UFL forms, stabilization sign, source and boundary convention remain those already defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acc11e75",
   "metadata": {},
   "outputs": [],
   "source": [
    "resolved_cases, resolved_rows = {}, []\n",
    "FACE_SUBDIVISIONS, RESOLVED_LOCAL_SUBDIVISIONS = 4, 8\n",
    "RESOLVED_LEVELS = {1e-5: (4, 8, 16)}\n",
    "RESOLVED_LEVELS.setdefault(EPSILON, (16,))\n",
    "for epsilon, resolutions in RESOLVED_LEVELS.items():\n",
    "    for n in resolutions:\n",
    "        for method, stabilized in METHODS:\n",
    "            case, row, archive = measure_case(\n",
    "                *solve_rad_case(\n",
    "                    n,\n",
    "                    epsilon,\n",
    "                    stabilized,\n",
    "                    face_subdivisions=FACE_SUBDIVISIONS,\n",
    "                    local_subdivisions=RESOLVED_LOCAL_SUBDIVISIONS,\n",
    "                ),\n",
    "                ReactionLayer(epsilon),\n",
    "                order=ERROR_ORDER,\n",
    "                name=f\"rad-resolved-eps{epsilon:g}-n{n}-{method.lower()}\",\n",
    "                epsilon=epsilon,\n",
    "                method=method,\n",
    "                n=n,\n",
    "                H=1 / n,\n",
    "                face_subdivisions=FACE_SUBDIVISIONS,\n",
    "                local_subdivisions=RESOLVED_LOCAL_SUBDIVISIONS,\n",
    "            )\n",
    "            row[\"overshoot\"] = max(0.0, row[\"nodal_maximum\"] - 1)\n",
    "            resolved_cases[epsilon, n, method] = case\n",
    "            resolved_rows.append(row)\n",
    "            state_archives.append(archive)\n",
    "            print(\"Refined P1/P0\", epsilon, n, method, row)\n",
    "resolved_rows\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc5642d0",
   "metadata": {},
   "source": [
    "### Inspect the refined physical fields and measured rates\n",
    "\n",
    "The following severe-layer panels use the finest declared P1/P0 family. They show both methods, their scalar errors and the independently refined CG2 baseline. Every panel overlays the actual 512-triangle macro mesh. The convergence figure reports observed scalar and physical-flux errors and slopes; no smooth-problem rate is transferred to this singularly perturbed case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7612f014",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_rad_refined(resolved_cases, resolved_rows, references, ReactionLayer, reports=REPORTS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6bc689e",
   "metadata": {},
   "source": [
    "### Verify physical integration at both resolution extremes\n",
    "\n",
    "The severe coarse control uses Duffy rules with 40 and 48 Gaussian points per coordinate. Both finest refined methods use 24 and 32. These comparisons verify the reported scalar and raw physical-flux norms; an algebraic residual is recorded separately. Each solved state also archives its executed local basis and passes replay with one and two BLAS threads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "707fa33d",
   "metadata": {},
   "outputs": [],
   "source": [
    "quadrature_controls: list[dict[str, Any]] = []\n",
    "for family, orders in ((\"underresolved\", (40, 48)), (\"refined\", (24, 32))):\n",
    "    for method in (\"MHM-Galerkin\", \"MHM-USFEM\"):\n",
    "        if family == \"underresolved\":\n",
    "            meshes, fields = challenging[1e-5, method]\n",
    "            measured = next(\n",
    "                row for row in sweep if row[\"epsilon\"] == 1e-5 and row[\"method\"] == method\n",
    "            )\n",
    "        else:\n",
    "            meshes, fields = resolved_cases[1e-5, 16, method][2:]\n",
    "            measured = next(\n",
    "                row\n",
    "                for row in resolved_rows\n",
    "                if row[\"epsilon\"] == 1e-5 and row[\"n\"] == 16 and row[\"method\"] == method\n",
    "            )\n",
    "        higher = scalar_error_norms(\n",
    "            meshes,\n",
    "            fields,\n",
    "            1,\n",
    "            ReactionLayer(1e-5).value,\n",
    "            ReactionLayer(1e-5).gradient,\n",
    "            1e-5,\n",
    "            order=orders[1],\n",
    "        )\n",
    "        norms = (\"scalar_l2\", \"flux_l2\")\n",
    "        np.testing.assert_allclose(\n",
    "            [higher[norm] for norm in norms],\n",
    "            [measured[norm] for norm in norms],\n",
    "            rtol=1e-8,\n",
    "        )\n",
    "        quadrature_controls.append(\n",
    "            dict(\n",
    "                family=family,\n",
    "                method=method,\n",
    "                epsilon=1e-5,\n",
    "                orders=list(orders),\n",
    "                norms={norm: [measured[norm], higher[norm]] for norm in norms},\n",
    "            )\n",
    "        )\n",
    "quadrature_controls\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e9db587",
   "metadata": {},
   "source": [
    "## 8. Compare coarse and refined one-sided profiles\n",
    "\n",
    "The top row preserves the $n=8$, unsplit-face, two-local-subdivision control. The bottom row uses $n=16$, four P0 segments per face and eight local subdivisions. Each panel includes the exact field and its corresponding fine CG2 reference. These meshes remain much smaller in the global multiscale solve than the reference mesh.\n",
    "\n",
    "Horizontal profiles retain a separate segment for each incident macrocell. Vertical markers identify actual macroface crossings; neighboring endpoint values are evaluated independently instead of being averaged. Profiles need not agree pointwise across a macroface because the global equations enforce the declared trace moments.\n",
    "\n",
    "For the coarse P1/P0 control at $\\epsilon=10^{-5}$, MHM-USFEM has smaller overshoot and maximum nodal error than MHM-Galerkin, but a larger scalar $L^2$ error. At $\\epsilon=10^{-3}$, its maximum nodal error is larger.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a33a189",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "provenance = execution_provenance(\n",
    "    \"notebooks/introduction/mhm_usfem_rad.ipynb\",\n",
    "    root=ROOT,\n",
    "    execution={\n",
    "        \"policy\": RAD_POLICY,\n",
    "        \"parallel_backend\": RAD_EXECUTION.backend,\n",
    "        \"automatic_process_configuration\": \"at least 512 macroelements, 8 local subdivisions and 4 face segments\",\n",
    "        \"workers\": RAD_EXECUTION.effective_workers,\n",
    "        \"native_threads\": RAD_EXECUTION.native_threads,\n",
    "        \"pipeline\": RAD_EXECUTION.pipeline,\n",
    "        \"pending_bound\": RAD_EXECUTION.effective_batch_size,\n",
    "        \"workers_capped_to_macroelement_count\": True,\n",
    "        \"provider_source_sha256\": RAD_PROVIDER_SOURCE_SHA256,\n",
    "        \"jit_cache\": \"notebook-owned temporary directory; serial two-orientation warm-up\",\n",
    "    },\n",
    "    literature_comparison=\"same analytical PDE; declared SW-NE topology and strong exterior Dirichlet data; historical matching local triangulation/connectivity is unresolved; no matched-figure reproduction\",\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00a63d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_rad_profiles(\n",
    "    macro,\n",
    "    challenging,\n",
    "    resolved_cases,\n",
    "    references,\n",
    "    ReactionLayer,\n",
    "    epsilon_primary=EPSILON,\n",
    "    underresolved_subdivisions=LOCAL_SUBDIVISIONS,\n",
    "    face_subdivisions=FACE_SUBDIVISIONS,\n",
    "    refined_subdivisions=RESOLVED_LOCAL_SUBDIVISIONS,\n",
    "    reports=REPORTS,\n",
    ")\n",
    "_ = (REPORTS / \"mhm-usfem-rad.json\").write_text(\n",
    "    json.dumps(\n",
    "        {\n",
    "            \"problem\": \"CILAMCE 2025 section 4.1 analytical reaction layer\",\n",
    "            \"source_url\": \"https://doi.org/10.55592/cilamce2025.v5i.14270\",\n",
    "            \"epsilon\": EPSILON,\n",
    "            \"local_degree\": 1,\n",
    "            \"trace_degree\": 0,\n",
    "            \"local_subdivisions\": LOCAL_SUBDIVISIONS,\n",
    "            \"local_ufl_quadrature_degree\": LOCAL_QUADRATURE_DEGREE,\n",
    "            \"operator_check_duffy_order\": OPERATOR_CHECK_ORDER,\n",
    "            \"reference_backend\": f\"DOLFINx {dolfinx.__version__}; independent CG2 UFL\",\n",
    "            \"boundary\": \"strong u=0 vertical, zero physical flux horizontal\",\n",
    "            \"convergence\": rows,\n",
    "            \"reference_refinement\": reference_rows,\n",
    "            \"provenance\": provenance,\n",
    "            \"state_archives\": state_archives,\n",
    "            \"reaction_sweep\": sweep,\n",
    "            \"underresolved_configuration\": {\n",
    "                \"macro_divisions\": 8,\n",
    "                \"face_subdivisions\": 1,\n",
    "                \"local_subdivisions\": LOCAL_SUBDIVISIONS,\n",
    "                \"role\": \"fixed-space resolution control; stabilization has no nodal-bound guarantee\",\n",
    "            },\n",
    "            \"refined_convergence\": resolved_rows,\n",
    "            \"refined_configuration\": {\n",
    "                \"face_subdivisions\": 4,\n",
    "                \"local_subdivisions\": 8,\n",
    "                \"compatibility\": f\"P1/P0; face subdivisions={FACE_SUBDIVISIONS}, local subdivisions={RESOLVED_LOCAL_SUBDIVISIONS}\",\n",
    "                \"boundary_variant\": \"strong exterior nodal Dirichlet data; continuous paper hybridization is weak\",\n",
    "                \"literature_status\": \"same PDE and declared admissible realization, not matched historical Fig2\",\n",
    "                \"unresolved_paper_input\": \"matching local triangulation and macro connectivity for Fig2\",\n",
    "            },\n",
    "            \"quadrature_controls\": quadrature_controls,\n",
    "        },\n",
    "        indent=2,\n",
    "    )\n",
    "    + \"\\n\"\n",
    ")\n",
    "\n",
    "# All local assembly is complete; release the notebook-owned compiler cache.\n",
    "sys.path.remove(RAD_JIT_CACHE.name)\n",
    "sys.modules.pop(RAD_WORKER_MODULE_NAME, None)\n",
    "RAD_JIT_CACHE.cleanup()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "767c87ad",
   "metadata": {},
   "source": [
    "## Qualify smooth convergence separately from the sharp layer\n",
    "\n",
    "An underresolved singular perturbation is not an asymptotic smooth-problem test. Keep the reaction-diffusion operator, but take $\\epsilon=\\sigma=1$, $f=1$, zero exterior scalar values on the two vertical boundaries, and homogeneous normal data on the horizontal boundaries:\n",
    "\n",
    "$$\n",
    "u(x,y)=1-\\frac{e^{-x}+e^{-(1-x)}}{1+e^{-1}}.\n",
    "$$\n",
    "\n",
    "The trace degree is $\\ell=0$. Local P1 scalar with one red refinement satisfies the admissible $k=\\ell+1$ refined alternative of theorem 1 in [Martins, Pacazuca and Valentin (2025)](https://doi.org/10.55592/cilamce2025.v5i.14270). Preserve the method's negative residual pairing. The macro sequence is $n=2,4,8,16,32$, with two local subdivisions at every level.\n",
    "\n",
    "The smooth broken-gradient reference order is one, as observed in the paper’s Figure 1; order-two scalar error is also observed. Theorem 1, estimate (21), includes an additional fine-cell Péclet factor and does not give a separate scalar L2 estimate. This qualified sequence does not claim a diffusion-uniform rate for the sharp layers considered above. The attributed record preserves its formulation, error integration and source provenance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "637bd223",
   "metadata": {},
   "source": [
    "### Recompute every measured order\n",
    "\n",
    "For successive physical errors $E_{i-1},E_i$ at the actual refinement sizes, compute\n",
    "\n",
    "$$\n",
    "r_i=\\frac{\\log(E_{i-1}/E_i)}{\\log(h_{i-1}/h_i)}.\n",
    "$$\n",
    "\n",
    "The first code lines read one attributed numerical record. They preserve every level and display the complete error/order table. The plotting helper only measures and plots these errors: it constructs no local or global PDE. Targets below are declared from the stated hypotheses, rather than fitted from the data. The final four measured levels give three consecutive orders and a fitted terminal slope. For each declared target $q$, a nearly constant $E/h^q$ provides a second view of the asymptotic regime.\n"
   ]
  },
  {
   "cell_type": "code",
   "id": "17c0f3e4",
   "metadata": {},
   "source": [
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Markdown, display\n",
    "from examples.tutorial_convergence import refinement_series, asymptotic_summary, plot_method_series\n",
    "refinement_record = ROOT / 'examples/results/unusual/analytical.json'\n",
    "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n",
    "refinement_rows = [dict(H=2 ** 0.5 / refinement_row['n'], **refinement_row) for refinement_row in refinement_data['rows'] if refinement_row['case'] == 'smooth' and refinement_row['method'] == 'unusual']\n",
    "refinement_error_keys = {'scalar L2': 'scalar_l2', 'broken gradient L2': 'gradient_l2'}\n",
    "refinement_targets = {'broken gradient L2': 1.0}\n",
    "refinement_rows = sorted(refinement_rows, key=lambda refinement_row: refinement_row['H'], reverse=True)\n",
    "refinement_sizes = np.asarray([refinement_row['H'] for refinement_row in refinement_rows])\n",
    "refinement_field_errors = {refinement_field: np.asarray([refinement_row[refinement_key] for refinement_row in refinement_rows]) for refinement_field, refinement_key in refinement_error_keys.items()}\n",
    "refinement_orders = {refinement_field: np.log(refinement_values[:-1] / refinement_values[1:]) / np.log(refinement_sizes[:-1] / refinement_sizes[1:]) for refinement_field, refinement_values in refinement_field_errors.items()}\n",
    "refinement_header = ['Refinement size'] + [refinement_column for refinement_field in refinement_error_keys for refinement_column in (refinement_field, 'Order')]\n",
    "refinement_lines = [' | '.join(refinement_header), ' | '.join(['---'] * len(refinement_header))]\n",
    "for refinement_index, refinement_size in enumerate(refinement_sizes):\n",
    "    refinement_values = [f'{refinement_size:.6g}']\n",
    "    for refinement_field in refinement_error_keys:\n",
    "        refinement_values.extend([f'{refinement_field_errors[refinement_field][refinement_index]:.6e}', '—' if refinement_index == 0 else f'{refinement_orders[refinement_field][refinement_index - 1]:.3f}'])\n",
    "    refinement_lines.append(' | '.join(refinement_values))\n",
    "display(Markdown('\\n'.join(refinement_lines)))\n",
    "refinement_series_data = refinement_series(refinement_record, refinement_rows, 'H', refinement_error_keys, refinement_targets, method='mhm-usfem', spaces='P1/r2 locals; P0 trace; epsilon=sigma=1; negative residual pairing', rate_provenance='Martins, Pacazuca and Valentin (2025), Theorem 1: admissible reaction-diffusion energy bound with local Peclet factor; Figure 1 smooth numerical orders', root=ROOT)\n",
    "print({'record': refinement_series_data['record'], 'sha256': refinement_series_data['sha256']})\n",
    "for refinement_field, refinement_summary in asymptotic_summary(refinement_series_data).items():\n",
    "    print(refinement_field, {'last_three_orders': [round(refinement_order, 3) for refinement_order in refinement_summary['orders']], 'fitted_order': round(refinement_summary['fitted_order'], 3), 'target': refinement_summary['target'], 'normalized_amplitude_ratio': refinement_summary.get('amplitude_ratio')})\n",
    "plot_method_series(refinement_series_data)\n",
    "plt.show()\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "ad2fd3d1",
   "metadata": {},
   "source": [
    "The upper panel preserves all coarse and fine measurements. Dashed lines show declared target powers anchored at the finest measured error. The shaded interval always contains the final four levels; it does not select points to improve a fitted slope. Read the consecutive orders together with the target-normalized errors and the stated quadrature/local-equation controls. A slope alone does not establish the hypotheses of an error estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "references-mhm_usfem_rad",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "- Larissa Martins, Juan F. Pacazuca, and Frédéric Valentin (2025). *A Multiscale Hybrid-Mixed Method with Local Stabilization*. Proceedings of the Ibero-Latin American Congress on Computational Methods in Engineering, CILAMCE 2025, volume 5, article 14270; published online 18 March 2026. [DOI: 10.55592/cilamce2025.v5i.14270](https://doi.org/10.55592/cilamce2025.v5i.14270).\n",
    "\n",
    "- L. P. Franca and Frédéric Valentin (2000). *On an improved unusual stabilized finite element method for the advective–reactive–diffusive equation*. Computer Methods in Applied Mechanics and Engineering 190(13–14), 1785–1800. [DOI: 10.1016/S0045-7825(00)00190-0](https://doi.org/10.1016/S0045-7825(00)00190-0)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (Pixi introduction)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "pymhm_editorial_revision": {
   "date": "2026-10-06",
   "scope": "Markdown citations and bibliography only",
   "preserved": "Code cells, execution counts, execution timestamps and numerical outputs",
   "preserved_code_cells_sha256": "cc96b267e57bc0b205b1c871e58f3f320a1ad634a3fe521a4e41782c67e6f7b6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}