{ "cells": [ { "cell_type": "markdown", "id": "fbeb3592", "metadata": {}, "source": [ "# Flux reconstruction, error indicators and adaptive decisions\n", "\n", "This lesson starts with user-written UFL local/global equations. It then defines a physical coefficient record, reconstructs an H(div) flux, evaluates an estimator and chooses cells for refinement. Numerical algorithms remain in the package; no physical-model constructor selects the equations.\n", "\n", "The primary case is $p=\\sin(2\\pi x)\\sin(2\\pi y)$, $K=I$, $f=8\\pi^2p$, with homogeneous Dirichlet pressure on the unit square. The formulation follows [Harder, Paredes and Valentin (2013)](https://doi.org/10.1016/j.jcp.2013.03.019). Flux recovery and the four-term estimator follow [Barrenechea, Martins, Pereira and Valentin (2026)](https://doi.org/10.1137/24M1673073).\n", "\n", "Install the native UFL backend as explained in the installation guide. The workflow is **meshes → local forms → global forms → assemble → solve → recovery → indicators → marking**." ] }, { "cell_type": "code", "execution_count": null, "id": "companion-bootstrap", "metadata": { "tags": [ "pymhm-companion-bootstrap" ] }, "outputs": [], "source": [ "from pathlib import Path\n", "import os\n", "import sys\n", "from pymhm.io.workspace import workspace_from_archive\n", "\n", "# Preserve this exact notebook as the numerical source checked by its records.\n", "NOTEBOOK_RELATIVE = Path(\"notebooks/darcy/reconstruction_and_indicators.ipynb\")\n", "source_declaration = os.environ.get(\"PYMHM_NOTEBOOK_SOURCE\")\n", "source_candidates = (\n", " (Path(source_declaration),) if source_declaration else\n", " (Path.cwd() / NOTEBOOK_RELATIVE.name, Path.cwd() / NOTEBOOK_RELATIVE)\n", ")\n", "NOTEBOOK_SOURCE = next((path for path in source_candidates if path.is_file()), None)\n", "if NOTEBOOK_SOURCE is None:\n", " raise FileNotFoundError(\n", " \"Keep this notebook in the launch directory, or set PYMHM_NOTEBOOK_SOURCE \"\n", " \"to its actual file before starting the kernel.\"\n", " )\n", "notebook_source_bytes = NOTEBOOK_SOURCE.read_bytes()\n", "\n", "# Download checksum-verified support files; extraction does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/4ad6aef8aa4059979d69bc5584e87cc7f20954afc1980aa30902d1d1839087aa/reconstruction_and_indicators-companion.zip\"\n", "COMPANION_SHA256 = \"4ad6aef8aa4059979d69bc5584e87cc7f20954afc1980aa30902d1d1839087aa\"\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", "# Prepare only this notebook's declared inputs, without executing its equations.\n", "from scripts.notebook_reproduction import notebook_workspace\n", "\n", "ROOT = notebook_workspace(\"darcy/reconstruction_and_indicators.ipynb\", directory=ROOT)\n", "\n", "# Stage actual source bytes, including their independent provenance fingerprint.\n", "workspace_notebook = ROOT / NOTEBOOK_RELATIVE\n", "workspace_notebook.parent.mkdir(parents=True, exist_ok=True)\n", "workspace_notebook.write_bytes(notebook_source_bytes)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "76e5d6b2", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import ufl\n", "from pymhm import (\n", " Equation, FaceSpace, LocalContext, MeshHierarchy, SkeletonSpace, TriangleMesh,\n", " assemble, bind_interface, bind_problem, columns, solve,\n", ")\n", "from pymhm.postprocessing.solutions import DarcySolution\n", "from pymhm.recovery.moments import reconstruct_darcy_moments\n", "from pymhm.estimators.darcy import estimate_darcy_error\n", "from pymhm.adaptivity.darcy import mark_dorfler\n", "from threadpoolctl import threadpool_limits\n", "\n", "\n", "def exact_pressure(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Analytical pressure; it vanishes on every exterior edge.\"\"\"\n", " return np.sin(2 * np.pi * points[:, 0]) * np.sin(2 * np.pi * points[:, 1])\n", "\n", "\n", "def exact_gradient(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Independently differentiated physical pressure gradient.\"\"\"\n", " x, y = 2 * np.pi * points.T\n", " return 2 * np.pi * np.column_stack((np.cos(x) * np.sin(y), np.sin(x) * np.cos(y)))\n", "\n", "\n", "def source(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Negative Laplacian computed independently of the discrete operator.\"\"\"\n", " return 8 * np.pi**2 * exact_pressure(points)\n", "RUN_STUDY = os.environ.get(\"PYMHM_NOTEBOOK_STUDY\", \"0\") == \"1\"\n" ] }, { "cell_type": "markdown", "id": "cc81c629", "metadata": {}, "source": [ "## 1. Choose meshes and approximation spaces\n", "\n", "For reconstruction degree $m=2$ and trace degree $\\ell=1$, choose local degree $k=3$. This satisfies the two-dimensional restriction $k\\ge\\ell+2$ and $\\ell\\le m\\le k$. The estimator additionally needs convex macrotriangles, a globally conforming fine partition, identity diffusion and homogeneous Dirichlet data. These are mathematical restrictions, not conventions inferred by the software." ] }, { "cell_type": "code", "execution_count": null, "id": "c1045b06", "metadata": {}, "outputs": [], "source": [ "macro = TriangleMesh.unit_square(8)\n", "local_meshes = tuple(macro.submesh(cell, 2) for cell in range(len(macro.cells)))\n", "hierarchy = MeshHierarchy(macro, local_meshes)\n", "skeleton = SkeletonSpace(macro, tuple(FaceSpace.uniform(1) for _ in macro.faces))\n", "interface = bind_interface(skeleton, convention=\"normal\")" ] }, { "cell_type": "markdown", "id": "543be992", "metadata": {}, "source": [ "## 2. Declare local equations\n", "\n", "$$\n", "\\begin{aligned}\n", "(\\nabla p_T,\\nabla v)_T+\\langle\\lambda_T,v\\rangle_{\\partial T}&=(f,v)_T,\\\\\n", "-\\langle p_T,\\mu_T\\rangle_{\\partial T}&=g_T(\\mu_T).\n", "\\end{aligned}\n", "$$\n", "\n", "The local volume form has a constant kernel. Its physical volume moment fixes the local complement and retains one mean per macrocell. The interface binding supplies normal incidence and numbering; the minus sign in the global test pairing is part of the formulation." ] }, { "cell_type": "code", "execution_count": null, "id": "a13d2886", "metadata": {}, "outputs": [], "source": [ "def local_equations(local: LocalContext, degree: int = 3, forcing=None):\n", " \"\"\"Declare primal diffusion, independent trace tests and physical mean.\"\"\"\n", " space = local.native_space(degree=degree)\n", " p, v = ufl.TrialFunction(space.space), ufl.TestFunction(space.space)\n", " x = ufl.SpatialCoordinate(space.mesh)\n", " dx = ufl.Measure(\"dx\", domain=space.mesh, metadata={\"quadrature_degree\": 20})\n", " forcing = (8 * np.pi**2 * ufl.sin(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])\n", " if forcing is None else forcing(space.mesh))\n", " local.field(\"pressure\", space)\n", " return local.equations(\n", " a=ufl.inner(ufl.grad(p), ufl.grad(v)) * dx,\n", " L=forcing * v * dx,\n", " b=local.trace_pairings(lambda phi, ds: phi * v * ds),\n", " c=local.trace_pairings(lambda phi, ds: -phi * p * ds, axis=\"rows\"),\n", " kernel=np.ones((space.size, 1)), moments=columns(v * dx),\n", " )" ] }, { "cell_type": "markdown", "id": "7f732dae", "metadata": {}, "source": [ "## 3. Declare global boundary moments\n", "\n", "$$\n", "-\\sum_T\\langle p_T,\\mu_T\\rangle_{\\partial T}\n", "=-\\langle p_D,\\mu\\rangle_{\\partial\\Omega},\\qquad p_D=0.\n", "$$\n", "\n", "The Dirichlet load is therefore zero. Its sign is kept explicit so that changing the boundary data does not change the formulation accidentally." ] }, { "cell_type": "code", "execution_count": null, "id": "14661ea7", "metadata": {}, "outputs": [], "source": [ "def global_equation(global_context):\n", " \"\"\"Supply the weak Dirichlet moments in the bound interface layout.\"\"\"\n", " boundary, fixed = global_context.boundary_data(0.0, order=8)\n", " assert not fixed\n", " return Equation(0, global_context.trace_load(-boundary))\n", "\n", "\n", "problem = bind_problem(hierarchy, interface, local_equations,\n", " global_equation=global_equation, retained=1)" ] }, { "cell_type": "markdown", "id": "e9109ea7", "metadata": {}, "source": [ "## 4. Assemble, solve and name the physical result\n", "\n", "Named fields preserve the native executed coefficient basis. `portable_coefficients` converts through its stored mapping; reshaping native vectors would not be equivalent. `DarcySolution` below is a data record for established reconstruction and estimator operations. Constructing it performs no PDE solve.\n", "\n", "The `raw_residual` printed below is the original **global compatibility** diagnostic for trace and retained equations. It is distinct from every original local volume row and from a physical L2 error. The companion native original-row observer separately checks the same eight-subdivision P2/P0 and P3/P1 equations; their maximum relative local defects are 3.32e-15 and 1.04e-14, with full uncondensed defects relative to the physical load 5.67e-15 and 1.76e-14. The localized control also passes: 4.20e-15 and 6.12e-15. All use the unchanged 1e-10 equation gate.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c83d0324", "metadata": {}, "outputs": [], "source": [ "with threadpool_limits(1):\n", " system = assemble(problem)\n", " coefficients = solve(system)\n", "pressure_fields = coefficients.field(\"pressure\")\n", "solution = DarcySolution(\n", " skeleton=skeleton, local_meshes=local_meshes,\n", " pressure=tuple(field.portable_coefficients for field in pressure_fields),\n", " flux=(), hybrid=coefficients, formulation=\"primal\", permeability=1.0,\n", " source=source, degree=3, quadrature_order=10,\n", ")\n", "print({\"original_equation_residual\": coefficients.raw_residual,\n", " \"pressure_L2\": solution.l2_error(exact_pressure, order=12)})" ] }, { "cell_type": "markdown", "id": "2ee81265", "metadata": {}, "source": [ "## 5. Reconstruct a flux and distinguish its conservation tests\n", "\n", "The canonical RT moment reconstruction preserves boundary normal moments from the skeleton, averages interior normal moments and preserves raw interior vector moments. It is not the energy-minimizing RT0 equilibration. Its divergence balance is against continuous macro-local tests, not individual discontinuous fine-cell constants.\n", "\n", "$$\n", "\\begin{aligned}\n", "\\langle q_h\\cdot n,\\phi\\rangle_E&=\\text{declared normal moments},\\\\\n", "(q_h,\\psi)_t&=(-\\nabla p_h,\\psi)_t.\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "e5949f07", "metadata": {}, "outputs": [], "source": [ "recovered = reconstruct_darcy_moments(solution, degree=2, quadrature_order=10)\n", "normal_defect = max(np.max(abs(row)) for row in recovered.normal_flux_residuals())\n", "continuous_defect = max(np.max(abs(row)) for row in recovered.continuous_moment_residuals())\n", "fine_defect = max(np.max(abs(row)) for row in recovered.fine_conservation_residuals())\n", "assert normal_defect < 1e-10\n", "assert continuous_defect < 1e-10\n", "print({\"normal_moment_defect\": normal_defect,\n", " \"continuous_test_balance_defect\": continuous_defect,\n", " \"fine_cell_balance_defect_not_imposed\": fine_defect,\n", " \"recovered_flux_L2\": recovered.flux_l2_error(lambda x: -exact_gradient(x), order=12)})" ] }, { "cell_type": "markdown", "id": "a6c3b6b7", "metadata": {}, "source": [ "## 6. Evaluate the four-term indicator\n", "\n", "$$\n", "\\eta^2=\\sum_T\\big[(\\eta_{1,T}+\\eta_{3,T}+\\eta_{\\mathrm{osc},T})^2\n", "+\\eta_{2,T}^2\\big].\n", "$$\n", "\n", "The contributions measure flux defect, potential nonconformity, divergence-projection defect and source oscillation. The recovered potential is a separate continuous field. For this unit-diffusion case, compare the estimator with the independently integrated broken energy error; numerical quadrature does not provide an interval-certified bound." ] }, { "cell_type": "code", "execution_count": null, "id": "13198af0", "metadata": {}, "outputs": [], "source": [ "estimate = estimate_darcy_error(solution, homogeneous_dirichlet=True,\n", " degree=2, quadrature_order=10)\n", "energy_error = estimate.energy_error(exact_gradient, order=12)\n", "assert max(estimate.equilibrium_defect) < 1e-10\n", "assert estimate.total >= energy_error\n", "print({\"energy_error\": energy_error, \"estimator\": estimate.total,\n", " \"effectivity\": estimate.total / energy_error})" ] }, { "cell_type": "markdown", "id": "b3bae16b", "metadata": {}, "source": [ "## 7. Plot physical fields with the actual macro mesh\n", "\n", "Each field panel has its own color scale. Numerical values and errors are sampled independently in each local mesh, preserving interface sides." ] }, { "cell_type": "code", "execution_count": null, "id": "10207c43", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), layout=\"constrained\")\n", "all_exact, all_values, panels = [], [], []\n", "for fine, field in zip(local_meshes, pressure_fields, strict=True):\n", " points = fine.points[fine.cells].mean(axis=1)\n", " exact, values = exact_pressure(points), field.evaluate(points)\n", " panels.append((fine, exact, values))\n", " all_exact.extend(exact)\n", " all_values.extend(values)\n", "common = max(np.max(np.abs(all_exact)), np.max(np.abs(all_values)))\n", "error_max = max(np.max(abs(values - exact)) for _, exact, values in panels)\n", "for index, (ax, title) in enumerate(zip(axes, (\"Analytical pressure\", \"MHM pressure\", \"Pressure error\"), strict=True)):\n", " for fine, exact, values in panels:\n", " value = (exact, values, values - exact)[index]\n", " limit = error_max if index == 2 else common\n", " artist = ax.tripcolor(*fine.points.T, fine.cells, facecolors=value,\n", " vmin=-limit, vmax=limit, cmap=\"coolwarm\")\n", " for edge in macro.faces:\n", " ax.plot(*macro.points[edge].T, color=\"black\", lw=0.5, alpha=0.7)\n", " ax.set(title=title, xlabel=\"x\", ylabel=\"y\", aspect=\"equal\")\n", " fig.colorbar(artist, ax=ax)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "90de3acb", "metadata": {}, "source": [ "## 8. Mark cells; choose the discretization control to refine\n", "\n", "Dörfler marking selects a minimal sorted collection carrying the prescribed fraction of the squared indicator:\n", "\n", "$$\n", "\\sum_{T\\in\\mathcal M}\\eta_T^2\\ge\\theta\\sum_T\\eta_T^2.\n", "$$\n", "\n", "Marking is distinct from refining. Before changing a mesh, choose whether the dominant error comes from the macro partition, local approximation, or trace resolution. Refinement ancestry must preserve material and boundary markers. Rebuild the equations on the new spaces and compare error against work; a general adaptive sequence has no predetermined uniform-grid rate." ] }, { "cell_type": "code", "execution_count": null, "id": "b2d90cce", "metadata": {}, "outputs": [], "source": [ "marked = mark_dorfler(estimate.local_squared, theta=0.5)\n", "assert estimate.local_squared[marked].sum() >= 0.5 * estimate.local_squared.sum()\n", "fig, ax = plt.subplots(figsize=(5, 4), layout=\"constrained\")\n", "artist = ax.tripcolor(*macro.points.T, macro.cells, facecolors=estimate.local_squared,\n", " edgecolors=\"black\", cmap=\"viridis\")\n", "centers = macro.points[macro.cells[marked]].mean(axis=1)\n", "ax.scatter(*centers.T, marker=\"x\", color=\"white\", s=60, label=\"Dörfler marked\")\n", "ax.set(xlabel=\"x\", ylabel=\"y\", title=\"Squared macro indicators\", aspect=\"equal\")\n", "ax.legend()\n", "fig.colorbar(artist, ax=ax)\n", "plt.show()\n", "print({\"marked_macro_cells\": np.flatnonzero(marked).tolist(), \"marked_fraction\":\n", " float(estimate.local_squared[marked].sum() / estimate.local_squared.sum())})" ] }, { "cell_type": "markdown", "id": "af727e74", "metadata": {}, "source": [ "## 9. Refine a localized physical solution, rebuild and compare work\n", "\n", "Uniform refinement first establishes the expected smooth-case order. Adaptivity\n", "then uses a different smooth pressure with a localized peak,\n", "\n", "$$\n", "p(x,y)=x(1-x)y(1-y)\\exp\\!\\left[-80\\big((x-0.3)^2+(y-0.4)^2\\big)\\right].\n", "$$\n", "\n", "It vanishes on the entire boundary. The following NumPy source is differentiated\n", "independently from the UFL source. Both define the same unit-diffusion operator.\n", "This manufactured adaptive test uses the estimator of Barrenechea et al. (2026);\n", "it is an analytical extension, rather than their heterogeneous-material example.\n", "The error is compared against the exact gradient and against uniform meshes at\n", "similar work. Adaptive mesh diameters are nonuniform: the horizontal coordinate\n", "is the total number of retained and trace unknowns, not a fictitious uniform $H$." ] }, { "cell_type": "code", "execution_count": null, "id": "a4e7c426", "metadata": {}, "outputs": [], "source": [ "from functools import partial\n", "from time import perf_counter\n", "from pymhm.meshes.refinement import refine_triangles\n", "from examples.tutorial_recovery_acquisition import recovery_diagnostics, recovery_display_samples, read_recovery_record\n", "\n", "\n", "def solve_declared(macro, ell, source_values, forcing=None):\n", " \"\"\"Assemble the already declared forms on caller-chosen meshes and spaces.\"\"\"\n", " meshes = tuple(macro.submesh(i, 2) for i in range(len(macro.cells)))\n", " trace = SkeletonSpace(macro, tuple(FaceSpace.uniform(ell) for _ in macro.faces))\n", " problem = bind_problem(\n", " MeshHierarchy(macro, meshes), bind_interface(trace, convention=\"normal\"),\n", " partial(local_equations, degree=ell + 2, forcing=forcing),\n", " global_equation=global_equation, retained=1,\n", " )\n", " with threadpool_limits(1):\n", " coefficients = solve(assemble(problem))\n", " result = DarcySolution(\n", " skeleton=trace, local_meshes=meshes,\n", " pressure=tuple(field.portable_coefficients for field in coefficients.field(\"pressure\")),\n", " flux=(), hybrid=coefficients, formulation=\"primal\", permeability=1.0,\n", " source=source_values, degree=ell + 2, quadrature_order=16,\n", " )\n", " return result\n", "\n", "\n", "def localized_pressure(points):\n", " \"\"\"A boundary-compatible smooth peak at (0.3, 0.4).\"\"\"\n", " x, y = points.T\n", " return x * (1-x) * y * (1-y) * np.exp(-80*((x-.3)**2 + (y-.4)**2))\n", "\n", "\n", "def localized_gradient(points):\n", " \"\"\"Independent first derivatives of the exact localized pressure.\"\"\"\n", " x, y = points.T\n", " g = np.exp(-80*((x-.3)**2 + (y-.4)**2))\n", " X, Y = x*(1-x), y*(1-y)\n", " return g[:, None] * np.column_stack((Y*(1-2*x-160*(x-.3)*X), X*(1-2*y-160*(y-.4)*Y)))\n", "\n", "\n", "def localized_source(points):\n", " \"\"\"Independent negative Laplacian, rather than a discrete manufactured load.\"\"\"\n", " x, y = points.T\n", " g = np.exp(-80*((x-.3)**2 + (y-.4)**2))\n", " X, Y = x*(1-x), y*(1-y)\n", " dxx = Y * (-2-320*(x-.3)*(1-2*x) + X*(25600*(x-.3)**2-160))\n", " dyy = X * (-2-320*(y-.4)*(1-2*y) + Y*(25600*(y-.4)**2-160))\n", " return -(dxx+dyy)*g\n", "\n", "\n", "def localized_forcing(domain):\n", " \"\"\"Differentiate the analytical UFL expression independently of NumPy.\"\"\"\n", " x = ufl.SpatialCoordinate(domain)\n", " p = x[0]*(1-x[0])*x[1]*(1-x[1])*ufl.exp(-80*((x[0]-.3)**2+(x[1]-.4)**2))\n", " return -ufl.div(ufl.grad(p))\n", "\n", "\n", "if RUN_STUDY:\n", " adaptive_rows = []\n", " adaptive_mesh = TriangleMesh.unit_square(4)\n", " for step in range(8):\n", " with threadpool_limits(1):\n", " adaptive_solution = solve_declared(adaptive_mesh, 1, localized_source, localized_forcing)\n", " row, adaptive_estimate = recovery_diagnostics(\n", " adaptive_solution, localized_pressure, localized_gradient,\n", " )\n", " row[\"step\"] = step\n", " row[\"global_unknowns\"] = row[\"trace_dofs\"] + row[\"macro_cells\"]\n", " if adaptive_rows:\n", " assert row[\"energy_error\"] < adaptive_rows[-1][\"energy_error\"]\n", " adaptive_rows.append(row)\n", " print({key: row[key] for key in (\"step\", \"macro_cells\", \"global_unknowns\", \"energy_error\", \"effectivity\")})\n", " if step < 7:\n", " marked = mark_dorfler(adaptive_estimate.local_squared, theta=0.5)\n", " assert adaptive_estimate.local_squared[marked].sum() >= 0.5*adaptive_estimate.local_squared.sum()\n", " adaptive_mesh = refine_triangles(adaptive_mesh, marked).mesh\n", "else:\n", " adaptive_data = read_recovery_record(\n", " ROOT/\"examples/results/tutorial-methods/adaptivity-current.json\",\n", " notebook=ROOT/\"notebooks/darcy/reconstruction_and_indicators.ipynb\",\n", " )\n", " adaptive_rows = adaptive_data[\"adaptive\"]\n", " adaptive_mesh = TriangleMesh(adaptive_data[\"final_macro_points\"], adaptive_data[\"final_macro_cells\"])\n", " for row in adaptive_rows:\n", " print({key: row[key] for key in (\"step\", \"macro_cells\", \"global_unknowns\", \"energy_error\", \"effectivity\")})\n" ] }, { "cell_type": "markdown", "id": "0f5d94e5", "metadata": {}, "source": [ "## 10. Recover a conforming potential and a separately equilibrated RT0 flux\n", "\n", "The Oswald potential averages incident finite-element nodal traces in the global continuous P3 space, then imposes the declared homogeneous boundary data. Its global coefficient vector supplies a single value at every shared node. This is a different field from the original broken pressure.\n", "\n", "For strict fine-cell balance, instead solve the weighted RT0 minimization with the same physical source and a **degree-zero aligned** skeletal trace. Reassemble those newly declared spaces first; do not pass incompatible P1 trace data to the RT0 operator. The minimization is an available reconstruction, not an assertion of equivalence to the published RT2 moment operator." ] }, { "cell_type": "code", "execution_count": null, "id": "d3f4a41f", "metadata": {}, "outputs": [], "source": [ "from pymhm.estimators.darcy import recover_potential\n", "from pymhm.recovery.equilibrated import equilibrate_flux\n", "from pymhm.fem.scalar.triangle import nodal_space\n", "\n", "potential = recover_potential(solution, homogeneous_dirichlet=True)\n", "_, potential_points = nodal_space(potential.mesh, potential.degree)\n", "boundary_nodes = np.any(np.isclose(potential_points, 0.0) | np.isclose(potential_points, 1.0), axis=1)\n", "assert np.max(abs(potential.values[boundary_nodes])) == 0.0\n", "print({\"conforming_potential_L2\": potential.l2_error(exact_pressure, order=12),\n", " \"potential_boundary_defect\": float(np.max(abs(potential.values[boundary_nodes])))})\n", "\n", "rt0_skeleton = SkeletonSpace(macro)\n", "rt0_interface = bind_interface(rt0_skeleton, convention=\"normal\")\n", "rt0_problem = bind_problem(hierarchy, rt0_interface, local_equations,\n", " global_equation=global_equation, retained=1)\n", "with threadpool_limits(1):\n", " rt0_system = assemble(rt0_problem)\n", " rt0_coefficients = solve(rt0_system)\n", "rt0_solution = DarcySolution(\n", " skeleton=rt0_skeleton, local_meshes=local_meshes,\n", " pressure=tuple(field.portable_coefficients for field in rt0_coefficients.field(\"pressure\")),\n", " flux=(), hybrid=rt0_coefficients, formulation=\"primal\", permeability=1.0,\n", " source=source, degree=3, quadrature_order=10,\n", ")\n", "equilibrated = equilibrate_flux(rt0_solution)\n", "fine_balance = max(np.max(abs(row)) for row in equilibrated.conservation_residuals())\n", "assert fine_balance < 1e-10\n", "print({\"equilibrated_fine_cell_balance\": fine_balance,\n", " \"equilibrated_Darcy_flux_L2\": equilibrated.l2_error(lambda x: -exact_gradient(x), order=12)})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 11. Establish the asymptotic regime before interpreting recovery\n", "\n", "Barrenechea et al. (2026), Theorem 4.7 and Remark 5.1, give\n", "\n", "$$\n", "\\begin{aligned}\n", "\\lVert q-q_h^{\\mathrm{RT}}\\rVert_{0,\\Omega}&\\le C(h^k+H^{\\ell+1}+h^{\\ell+1}),\\\\\n", "\\lVert f-\\Pi_{\\Omega,m}\\operatorname{div}q_h^{\\mathrm{RT}}\\rVert_{0,\\Omega}&\\le Ch^{m+1}.\n", "\\end{aligned}\n", "$$\n", "\n", "Here $h=H/2$, $k=\\ell+2$, and $m=2$. Thus the recovered physical flux has\n", "order one for $\\ell=0$ and two for $\\ell=1$, while its **projected** divergence\n", "has order three. The raw divergence is a distinct quantity and is reported\n", "separately. Local degree, RT degree and mesh ratios remain fixed throughout.\n", "The estimator's reliability bound and stable effectivity are checked at every\n", "level; no theorem requires a universal gain in recovered pressure error.\n", "\n", "The acquisition helper only measures already defined fields. It does not choose\n", "the operator or spaces. Set `PYMHM_NOTEBOOK_STUDY=1` to execute every study solve;\n", "otherwise the same current numerical record supplies the convergence figures." ] }, { "cell_type": "code", "metadata": {}, "source": [ "from examples.tutorial_recovery_acquisition import write_recovery_record, recovery_provenance\n", "from examples.tutorial_convergence import plot_method_series, asymptotic_summary\n", "\n", "study_path = ROOT / \"examples/results/tutorial-methods/recovery-current.json\"\n", "RUN_STUDY = os.environ.get(\"PYMHM_NOTEBOOK_STUDY\", \"0\") == \"1\"\n", "source_paths = [ROOT/\"notebooks/darcy/reconstruction_and_indicators.ipynb\", ROOT/\"examples/tutorial_recovery_acquisition.py\", ROOT/\"pixi.lock\"]\n", "if RUN_STUDY:\n", " executed_sources = recovery_provenance(ROOT, source_paths)\n", " uniform_rows = []\n", " for ell in (0, 1):\n", " for n in (4, 8, 16, 32):\n", " started = perf_counter()\n", " with threadpool_limits(1):\n", " state = solve_declared(TriangleMesh.unit_square(n), ell, source)\n", " row, state_estimate = recovery_diagnostics(state, exact_pressure, exact_gradient)\n", " row.update(macro_resolution=n, wall_seconds=perf_counter()-started)\n", " uniform_rows.append(row)\n", " print({key: row[key] for key in (\"trace_degree\", \"macro_resolution\", \"energy_error\", \"recovered_flux_l2\", \"projected_divergence_l2\", \"effectivity\", \"wall_seconds\")})\n", " write_recovery_record(\n", " study_path, root=ROOT,\n", " source_paths=source_paths, source_manifest=executed_sources,\n", " data={\"case\": \"sin(2*pi*x)sin(2*pi*y), identity diffusion, zero Dirichlet\",\n", " \"scope\": \"Current native UFL analytical study of the published section 6.1 data\",\n", " \"reference\": \"10.1137/24M1673073\", \"local_subdivisions\": 2,\n", " \"assembly_quadrature_degree\": 20, \"uniform\": uniform_rows},\n", " )\n", "else:\n", " uniform_rows = read_recovery_record(study_path, notebook=ROOT/\"notebooks/darcy/reconstruction_and_indicators.ipynb\")[\"uniform\"]\n", "\n", "for ell in (0, 1):\n", " rows = [row for row in uniform_rows if row[\"trace_degree\"] == ell]\n", " h = np.array([row[\"macro_diameter\"] for row in rows])\n", " errors = {\"recovered Darcy flux L2\": np.array([row[\"recovered_flux_l2\"] for row in rows]),\n", " \"broken energy error\": np.array([row[\"energy_error\"] for row in rows]),\n", " \"projected divergence L2\": np.array([row[\"projected_divergence_l2\"] for row in rows])}\n", " series = {\"method\": f\"RT2 recovery, trace degree {ell}\", \"sizes\": h.tolist(),\n", " \"errors\": {key: values.tolist() for key, values in errors.items()},\n", " \"rates\": {key: (np.log(values[:-1]/values[1:])/np.log(h[:-1]/h[1:])).tolist() for key, values in errors.items()},\n", " \"expected\": {\"recovered Darcy flux L2\": ell+1, \"broken energy error\": ell+1, \"projected divergence L2\": 3},\n", " \"refinement\": \"Macro diameter H, with h=H/2\",\n", " \"spaces\": f\"P{ell+2} local pressure; P{ell} normal trace; RT2 reconstruction\",\n", " \"rate_provenance\": \"Barrenechea et al. (2026), Theorem 4.7 and Remark 5.1\"}\n", " print(asymptotic_summary(series))\n", " plot_method_series(series)\n", " plt.show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 12. Compare a classical refined conforming Galerkin reference\n", "\n", "The baseline is assembled independently on a global continuous P3 space.\n", "The local MHM equations are not reused as its assembled matrix. Three global\n", "meshes verify its own refinement error against the analytical pressure and\n", "flux; the analytical solution remains the exact comparator. Each physical norm\n", "is integrated with two quadrature orders. This separates reference accuracy\n", "from the approximation quality of the coarse-interface method." ] }, { "cell_type": "code", "metadata": {}, "source": [ "from examples.introduction.scalar import native_scalar_space, dirichlet_solve\n", "from examples.introduction.transport import scalar_error_norms\n", "from pymhm.core.equations import compile_form\n", "from pymhm.fem.scalar.triangle import nodal_space\n", "\n", "classical_rows = []\n", "for n in (16, 32, 64):\n", " fine = TriangleMesh.unit_square(n)\n", " domain, V, mapping = native_scalar_space(fine, 3)\n", " p, v = ufl.TrialFunction(V), ufl.TestFunction(V)\n", " x = ufl.SpatialCoordinate(domain)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 20})\n", " f = 8*np.pi**2*ufl.sin(2*np.pi*x[0])*ufl.sin(2*np.pi*x[1])\n", " matrix = compile_form(ufl.inner(ufl.grad(p), ufl.grad(v))*dx)\n", " load = compile_form(f*v*dx)\n", " _, nodes = nodal_space(fine, 3)\n", " exterior = np.flatnonzero(np.any(np.isclose(nodes, 0)|np.isclose(nodes, 1), axis=1))\n", " with threadpool_limits(1):\n", " values = dirichlet_solve(matrix, load, mapping[exterior], np.zeros(len(exterior)))[mapping]\n", " low = scalar_error_norms((fine,), (values,), 3, exact_pressure, exact_gradient, 1.0, order=12)\n", " high = scalar_error_norms((fine,), (values,), 3, exact_pressure, exact_gradient, 1.0, order=16)\n", " assert abs(low[\"flux_l2\"]-high[\"flux_l2\"]) < 1e-7*high[\"flux_l2\"]\n", " classical_rows.append({\"n\": n, \"unknowns\": len(nodes), **high})\n", " print(classical_rows[-1])" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 13. Inspect the finest recovered fields and the adaptive error/work curve\n", "\n", "The analytical pressure, original MHM pressure, signed pressure error, analytical\n", "Darcy-flux magnitude, recovered RT2 Darcy-flux magnitude and vector-flux error\n", "are sampled separately on each fine triangle. The actual macro partition is\n", "shown in every panel. Convergence uses integrated errors, not these samples.\n", "The final adaptive mesh and local indicators show where the localized peak\n", "requires extra resolution. The uniform comparator solves the same localized\n", "operator and data at a fixed polynomial degree." ] }, { "cell_type": "code", "metadata": {}, "source": [ "from examples.tutorial_recovery_plots import plot_recovery_fields, plot_adaptive_study, plot_estimator_study\n", "\n", "with threadpool_limits(1):\n", " display_solution = solve_declared(TriangleMesh.unit_square(16), 1, source)\n", " display_row, display_estimate = recovery_diagnostics(display_solution, exact_pressure, exact_gradient)\n", " display_samples = recovery_display_samples(display_solution, display_estimate, exact_pressure, exact_gradient)\n", "plot_recovery_fields(display_samples)\n", "plt.show()\n", "plot_estimator_study(uniform_rows)\n", "plt.show()\n", "\n", "if RUN_STUDY:\n", " localized_uniform = []\n", " for n in (4, 8, 16, 32):\n", " with threadpool_limits(1):\n", " state = solve_declared(TriangleMesh.unit_square(n), 1, localized_source, localized_forcing)\n", " row, _ = recovery_diagnostics(state, localized_pressure, localized_gradient)\n", " row[\"global_unknowns\"] = row[\"trace_dofs\"] + row[\"macro_cells\"]\n", " localized_uniform.append(row)\n", " final_indicators = adaptive_estimate.local_squared\n", "else:\n", " localized_uniform = adaptive_data[\"uniform\"]\n", " final_indicators = np.asarray(adaptive_data[\"final_indicators\"])\n", "plot_adaptive_study(adaptive_rows, localized_uniform, adaptive_mesh, final_indicators)\n", "plt.show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 14. Use another multiscale problem as a local operator\n", "\n", "A hierarchy with outer and child macro partitions replaces each outer local solve by\n", "a global problem on four child macrotriangles. Leaf pressure remains P3 on two local edge subdivisions and\n", "child flux traces remain P1. Two independent P1 segments on each parent face\n", "span exactly the two child boundary traces. The automatic restriction preserves\n", "normal signs and segment coordinates; it is not an approximate projection.\n", "The child problem leaves exterior boundary values to its parent.\n", "\n", "The constant child pressure mode is zero on trace coordinates and one on each\n", "child retained mean. Its physical moment uses child areas. These moments retain the physical volume\n", "average of each nested problem's domain. The finest named pressure fields are\n", "obtained with `recursive=True`; no leaf is solved a second time during recovery.\n", "Nested condensation does not introduce a new convergence theorem. When its\n", "trace spaces span the same leaf discretization, it must reproduce the flattened\n", "MHM solution and inherit that discretization's physical error rate." ] }, { "cell_type": "code", "metadata": {}, "source": [ "def recursive_provider(local: LocalContext):\n", " \"\"\"Declare a child global problem inside this parent local macrotriangle.\"\"\"\n", " child_macro = local.mesh\n", " leaf_meshes = tuple(child_macro.submesh(i, 2) for i in range(len(child_macro.cells)))\n", " child_trace = SkeletonSpace(child_macro, tuple(FaceSpace.uniform(1) for _ in child_macro.faces))\n", " child = bind_problem(\n", " MeshHierarchy(child_macro, leaf_meshes), bind_interface(child_trace, convention=\"normal\"),\n", " local_equations, retained=1,\n", " )\n", " kernel = np.r_[np.zeros(child_trace.size), np.ones(len(child_macro.cells))][:, None]\n", " moments = np.r_[np.zeros(child_trace.size), child_macro.areas][:, None]\n", " return local.nested(child, kernel=kernel, moments=moments)\n", "\n", "\n", "if RUN_STUDY:\n", " recursive_rows = []\n", " for n in (2, 4, 8, 16):\n", " outer = TriangleMesh.unit_square(n)\n", " children = tuple(outer.submesh(i, 2) for i in range(len(outer.cells)))\n", " outer_trace = SkeletonSpace(outer, tuple(FaceSpace.uniform(1, 2) for _ in outer.faces))\n", " recursive_problem = bind_problem(\n", " MeshHierarchy(outer, children), bind_interface(outer_trace, convention=\"normal\"),\n", " recursive_provider, global_equation=global_equation, retained=1,\n", " )\n", " with threadpool_limits(1):\n", " recursive_coefficients = solve(assemble(recursive_problem))\n", " leaves = recursive_coefficients.field(\"pressure\", recursive=True)\n", " leaf_meshes = tuple(field.mesh for field in leaves)\n", " low = scalar_error_norms(leaf_meshes, leaves, 3, exact_pressure, exact_gradient, 1.0, order=12)\n", " high = scalar_error_norms(leaf_meshes, leaves, 3, exact_pressure, exact_gradient, 1.0, order=16)\n", " assert abs(low[\"flux_l2\"]-high[\"flux_l2\"]) < 1e-7*high[\"flux_l2\"]\n", " assert recursive_coefficients.raw_residual < 1e-10\n", " assert max(child.raw_residual for child in recursive_coefficients.children) < 1e-10\n", " # Flatten the identical child macro mesh, including its actual diagonals.\n", " points = np.concatenate([child.points for child in children])\n", " _, first, inverse = np.unique(points, axis=0, return_index=True, return_inverse=True)\n", " cells, offset = [], 0\n", " for child in children:\n", " cells.extend(inverse[offset+child.cells])\n", " offset += len(child.points)\n", " flat_mesh = TriangleMesh(points[first], np.array(cells))\n", " with threadpool_limits(1):\n", " flat = solve_declared(flat_mesh, 1, source)\n", " flat_fields = flat.hybrid.field(\"pressure\")\n", " discrepancy = max(float(np.max(abs(a.portable_coefficients-b.portable_coefficients)))\n", " for a, b in zip(leaves, flat_fields, strict=True))\n", " assert discrepancy < 1e-9\n", " row = {\"outer_resolution\": n, \"child_diameter\": np.sqrt(2)/(2*n),\n", " \"outer_trace_dofs\": outer_trace.size, \"leaf_macro_cells\": len(flat_mesh.cells),\n", " \"physical_errors\": high, \"flattened_coefficient_difference\": discrepancy,\n", " \"original_equation_residual\": recursive_coefficients.raw_residual}\n", " recursive_rows.append(row)\n", " print(row)\n", "\n", "else:\n", " recursive_rows = read_recovery_record(\n", " ROOT/\"examples/results/tutorial-methods/recursive-current.json\",\n", " notebook=ROOT/\"notebooks/darcy/reconstruction_and_indicators.ipynb\",\n", " )[\"recursive\"]\n", " for row in recursive_rows:\n", " print(row)\n", "\n", "h = np.array([row[\"child_diameter\"] for row in recursive_rows])\n", "errors = np.array([row[\"physical_errors\"][\"flux_l2\"] for row in recursive_rows])\n", "series = {\"method\": \"Recursive primal MHM\", \"sizes\": h.tolist(),\n", " \"errors\": {\"Darcy flux L2\": errors.tolist()},\n", " \"rates\": {\"Darcy flux L2\": (np.log(errors[:-1]/errors[1:])/np.log(h[:-1]/h[1:])).tolist()},\n", " \"expected\": {\"Darcy flux L2\": 2}, \"refinement\": \"Child macro diameter, fixed outer/child/local mesh ratios\",\n", " \"spaces\": \"P3 leaf pressure, P1 child traces, two P1 segments per parent face\",\n", " \"rate_provenance\": \"Exact equality to flattened primal MHM; Harder et al. (2013)\"}\n", "print(asymptotic_summary(series))\n", "plot_method_series(series)\n", "plt.show()" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "74fc045c", "metadata": {}, "source": [ "## References and extensions\n", "\n", "- Harder, Paredes and Valentin (2013), *A family of Multiscale Hybrid-Mixed finite element methods for the Darcy equation with rough coefficients*, [DOI](https://doi.org/10.1016/j.jcp.2013.03.019).\n", "- Barrenechea, Martins, Pereira and Valentin (2026), [DOI](https://doi.org/10.1137/24M1673073).\n", "\n", "For a declared SPD coefficient, use the energy-normalized estimator and supply certified ellipticity bounds when required. RT0 energy-minimizing equilibration instead imposes fine-cell averages and needs degree-zero trace segments aligned with fine boundary faces. Flow and elasticity indicators have their own hypotheses. The corresponding tutorials explain these contracts rather than silently applying the scalar indicator to other equations." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }