{ "cells": [ { "cell_type": "markdown", "id": "129e6481", "metadata": {}, "source": [ "# Robin MH: a Robin multiplier is not physical Darcy flux\n", "\n", "We write the Robin weak equations directly in UFL. Our independently manufactured data are $p=\\sin(2\\pi x)\\sin(2\\pi y)$, $K=I$, $f=8\\pi^2p$, and homogeneous exterior pressure. Follow **meshes → local forms → global forms → assemble → solve → physical fields**.\n", "\n", "The modified multiplier is $\\lambda_T=(q-p\\sigma)\\cdot n_T$, where $q=-\\nabla p$ and $\\sigma=\\nu(x-x_0)/2$. Consequently this multiplier cannot be plotted or tested as $q\\cdot n_T$. The method follows [Barrenechea, Gomes and Paredes (2024)](https://doi.org/10.1137/22M1542556)." ] }, { "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", "# Download checksum-verified support files; extraction does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/e4b740c6b2ca1ddd8b998294688f5e0c56f9053698db6a3951993ec778114fbc/robin_mh_workflow-companion.zip\"\n", "COMPANION_SHA256 = \"e4b740c6b2ca1ddd8b998294688f5e0c56f9053698db6a3951993ec778114fbc\"\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/robin_mh_workflow.ipynb\", directory=ROOT)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "bb8ff1f1", "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 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)" ] }, { "cell_type": "markdown", "id": "74b9d60f", "metadata": {}, "source": [ "## 1. Define macro, local and interface spaces\n", "\n", "The hierarchy owns mesh association. The interface binding owns geometric incidence and numbering; variational signs are declared separately." ] }, { "cell_type": "code", "execution_count": null, "id": "ed5da6d3", "metadata": {}, "outputs": [], "source": [ "macro = TriangleMesh.unit_square(2)\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": "e48343f9", "metadata": {}, "source": [ "## 2. Declare the local Robin equations\n", "\n", "$$\n", "(\\nabla p_T,\\nabla v)_T+\\langle\\sigma\\cdot n_T\\,p_T,v\\rangle_{\\partial T}\n", "+\\langle\\lambda_T,v\\rangle_{\\partial T}=(f,v)_T.\n", "$$\n", "\n", "Choose $x_0=(0,0)$ and $\\nu=1/4$. On this unit-square problem, $R=\\sqrt{2}/2$ gives the admissible sufficient bound $\\nu\\le 1/(4R^2)=1/2$. The Robin operator has no constant kernel: it retains no pressure mean. The global test remains the independent negative pressure trace pairing." ] }, { "cell_type": "code", "execution_count": null, "id": "ebec8e1e", "metadata": {}, "outputs": [], "source": [ "def local_equations(local: LocalContext):\n", " \"\"\"Declare Robin diffusion and its independent pressure trace test.\"\"\"\n", " space = local.native_space(degree=3)\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\": 12})\n", " forcing = 8 * np.pi**2 * ufl.sin(2 * np.pi * x[0]) * ufl.sin(2 * np.pi * x[1])\n", " ds = ufl.Measure(\"ds\", domain=space.mesh, metadata={\"quadrature_degree\": 12})\n", " normal = ufl.FacetNormal(space.mesh)\n", " sigma = 0.25 * x / 2\n", " local.field(\"pressure\", space)\n", " return local.equations(\n", " a=ufl.inner(ufl.grad(p), ufl.grad(v)) * dx + ufl.dot(sigma, normal) * p * v * ds,\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", " )" ] }, { "cell_type": "markdown", "id": "0cb70e89", "metadata": {}, "source": [ "## 3. Declare global pressure continuity and boundary moments\n", "\n", "$$\n", "-\\sum_T\\langle p_T,\\mu_T\\rangle_{\\partial T}=-\\langle p_D,\\mu\\rangle_{\\partial\\Omega}.\n", "$$\n", "\n", "The trace maps carry canonical incidence. The negative sign is the mathematical global-test convention. Zero retained coordinates reflect the invertible local Robin operator." ] }, { "cell_type": "code", "execution_count": null, "id": "f2d6e986", "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=0)" ] }, { "cell_type": "markdown", "id": "b6032099", "metadata": {}, "source": [ "## 4. Assemble and solve\n", "\n", "Assembly solves the independent local equations and sums their interface contributions. Solving reconstructs named local fields in their executed coefficient bases." ] }, { "cell_type": "code", "execution_count": null, "id": "d1415266", "metadata": {}, "outputs": [], "source": [ "with threadpool_limits(1):\n", " system = assemble(problem)\n", " coefficients = solve(system)\n", "pressure_fields = coefficients.field(\"pressure\")" ] }, { "cell_type": "markdown", "id": "6d8accfc", "metadata": {}, "source": [ "## 5. Integrate physical errors\n", "\n", "Pressure and its broken gradient are integrated separately with positive order-12 quadrature. No algebraic residual is used as a substitute for these errors." ] }, { "cell_type": "code", "execution_count": null, "id": "63dce12d", "metadata": {}, "outputs": [], "source": [ "from pymhm.fem.scalar.operators import triangle_quadrature\n", "bary, weights = triangle_quadrature(12)\n", "pressure_error, gradient_error = 0.0, 0.0\n", "for fine, field in zip(local_meshes, pressure_fields, strict=True):\n", " points = np.einsum(\"qi,tia->tqa\", bary, fine.points[fine.cells])\n", " flat = points.reshape(-1, 2)\n", " ep = (field.evaluate(flat) - exact_pressure(flat)).reshape(points.shape[:2])\n", " eg = (field.gradient(flat) - exact_gradient(flat)).reshape(*points.shape[:2], 2)\n", " pressure_error += np.einsum(\"t,q,tq->\", fine.areas, weights, ep**2)\n", " gradient_error += np.einsum(\"t,q,tq->\", fine.areas, weights, np.sum(eg**2, axis=-1))\n", "print({\"original_equation_residual\": coefficients.raw_residual,\n", " \"pressure_L2\": float(np.sqrt(pressure_error)),\n", " \"broken_gradient_L2\": float(np.sqrt(gradient_error))})\n", "assert coefficients.raw_residual < 1e-10" ] }, { "cell_type": "markdown", "id": "02962a7f", "metadata": {}, "source": [ "## 6. Plot analytical, numerical and error fields\n", "\n", "Independent local panels retain one-sided values. Every panel shows the actual macro partition." ] }, { "cell_type": "code", "execution_count": null, "id": "a31233db", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 2, figsize=(8.5, 7), 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", " flux_magnitude = np.linalg.norm(-field.gradient(points), axis=1)\n", " panels.append((fine, exact, values, flux_magnitude))\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", "flux_max = max(np.max(flux) for _, _, _, flux in panels)\n", "titles = (\"Analytical pressure\", \"Hybrid pressure\", \"Pressure error\", \"Hybrid Darcy flux magnitude\")\n", "for index, (ax, title) in enumerate(zip(axes.flat, titles, strict=True)):\n", " for fine, exact, values, flux in panels:\n", " value = (exact, values, values - exact, flux)[index]\n", " limit = (common, common, error_max, flux_max)[index]\n", " artist = ax.tripcolor(*fine.points.T, fine.cells, facecolors=value,\n", " vmin=0.0 if index == 3 else -limit,\n", " vmax=limit, cmap=\"viridis\" if index == 3 else \"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": "c4ad6236", "metadata": {}, "source": [ "## Independently refined classical primal Galerkin baseline\n", "\n", "The classical baseline uses one globally conforming P2 space on independent 16×16, 32×32 and 64×64 triangular grids, with the same physical coefficient, source and boundary values:\n", "\n", "$$\n", "(\\nabla p_h^{\\mathrm{CG}},\\nabla v_h)=(f,v_h),\\qquad\n", "p_h^{\\mathrm{CG}}\\big\\vert_{\\partial\\Omega}=p_D.\n", "$$\n", "\n", "The shared volume operator integrates this already stated form. The explicit essential lifting $A_{ff}p_f=f_f-A_{fb}p_D$ selects its classical boundary values. Check this reference's own errors before comparing fields. The multiscale solve above retains its simpler macro mesh and separate local spaces; it does not use this refined global grid. An exactly representable polynomial has roundoff-level classical errors on all three grids." ] }, { "cell_type": "code", "execution_count": null, "id": "efbb139c", "metadata": {}, "outputs": [], "source": [ "from threadpoolctl import threadpool_limits\n", "from pymhm.fem.scalar.triangle import scalar_operators, nodal_space, tabulate\n", "from pymhm.fem.scalar.operators import triangle_quadrature\n", "from pymhm.linalg.linear import solve_linear\n", "\n", "classical_rows = []\n", "bary, weights = triangle_quadrature(12)\n", "for resolution in (16, 32, 64):\n", " reference_mesh = TriangleMesh.unit_square(resolution)\n", " A, _, forcing = scalar_operators(reference_mesh, 2, diffusion=1.0, source=source, order=12)\n", " dofs, nodes = nodal_space(reference_mesh, 2)\n", " boundary = np.flatnonzero(np.any(np.isclose(nodes, 0.0) | np.isclose(nodes, 1.0), axis=1))\n", " free = np.setdiff1d(np.arange(len(nodes)), boundary)\n", " reference_pressure = np.zeros(len(nodes))\n", " reference_pressure[boundary] = exact_pressure(nodes[boundary])\n", " rhs = forcing[free] - A[free][:, boundary] @ reference_pressure[boundary]\n", " with threadpool_limits(1):\n", " reference_pressure[free] = solve_linear(A[free][:, free], rhs)\n", " _, _, basis, gradient, _ = tabulate(reference_mesh, 2, bary)\n", " points = np.einsum(\"qi,tia->tqa\", bary, reference_mesh.points[reference_mesh.cells])\n", " difference = reference_pressure[dofs] @ basis.T - exact_pressure(points.reshape(-1, 2)).reshape(points.shape[:2])\n", " error = float(np.sqrt(reference_mesh.areas @ (difference**2 @ weights)))\n", " classical_rows.append({\"resolution\": resolution, \"elements\": len(reference_mesh.cells),\n", " \"pressure_L2\": error})\n", "for row in classical_rows:\n", " print(row)\n", "# Polynomial patches can be exactly represented, in which case all levels are\n", "# at roundoff. Otherwise verify the classical reference's own refinement.\n", "assert (max(row[\"pressure_L2\"] for row in classical_rows) < 1e-10\n", " or all(b[\"pressure_L2\"] < a[\"pressure_L2\"] for a, b in zip(classical_rows, classical_rows[1:])))\n", "centers = reference_mesh.points[reference_mesh.cells].mean(axis=1)\n", "_, _, center_basis, center_gradient, _ = tabulate(reference_mesh, 2, np.full((1, 3), 1/3))\n", "reference_values = (reference_pressure[dofs] @ center_basis.T)[:, 0]\n", "reference_flux = -np.einsum(\"ti,tqia->tqa\", reference_pressure[dofs], center_gradient)[:, 0]\n", "fig, axes = plt.subplots(1, 2, figsize=(9, 4), layout=\"constrained\")\n", "for ax, values, title in zip(axes, (reference_values, np.linalg.norm(reference_flux, axis=1)),\n", " (\"Fine classical P2 pressure\", \"Fine classical Darcy flux magnitude\"), strict=True):\n", " artist = ax.tripcolor(*reference_mesh.points.T, reference_mesh.cells, facecolors=values, cmap=\"viridis\")\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": "416bbace", "metadata": {}, "source": [ "## 7. Keep the Robin and local refinement terms distinct\n", "\n", "The complete smooth acquisition uses the identity coefficient on the unit square with the paper's independently differentiated data $p=\\sin(6\\pi x)\\sin(14\\pi y)$ and $f=232\\pi^2p$, P3 local pressure, two local subdivisions, P1 Robin traces and $\\nu=1/4$. L-shaped **macroelements** partition the square; the physical domain has no reentrant corner.\n", "\n", "Corollary 2.7 of [Barrenechea, Gomes and Paredes (2024)](https://doi.org/10.1137/22M1542556) separates the local Galerkin and skeletal energy contributions:\n", "\n", "$$\n", "\\lVert K^{1/2}\\nabla(p-p_{Hh})\\rVert\n", "\\lesssim H^{\\ell+1}\\lvert K\\nabla p+p\\sigma\\rvert_{\\ell+1}\n", "+h^k\\lvert p\\rvert_{k+1}.\n", "$$\n", "\n", "Thus $\\ell=1$, $k=3$ and a fixed local subdivision ratio give second-order physical flux. Pressure order three is a measurement for these smooth data. The Robin multiplier is not the physical Darcy flux, so the error below uses $q=-K\\nabla p$. The record states its assembly/error quadrature, equation checks and execution sources. It is read as attributed numerical data; this cell does not execute the campaign again." ] }, { "cell_type": "markdown", "id": "ee70cae4", "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": "4ade1e6c", "metadata": {}, "source": [ "import json\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/mh/comparison.json'\n", "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n", "refinement_rows = [refinement_row for refinement_row in refinement_data['smooth'] if refinement_row['mesh'] == 'L-polygons' and refinement_row['trace_degree'] == 1]\n", "refinement_error_keys = {'pressure L2': 'pressure_error_l2', 'physical Darcy flux L2': 'flux_error_l2'}\n", "refinement_targets = {'physical Darcy flux L2': 2.0}\n", "refinement_rows = sorted(refinement_rows, key=lambda refinement_row: refinement_row['macro_diameter'], reverse=True)\n", "refinement_sizes = np.asarray([refinement_row['macro_diameter'] 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, 'macro_diameter', refinement_error_keys, refinement_targets, method='robin-mh', spaces='P3/r2 locals; P1 Robin trace; nu=1/4', rate_provenance='Barrenechea, Gomes and Paredes (2024), Corollary 2.7: energy order two under Assumption A1, matching regularly refined local partitions and the stated pressure/modified-flux regularity; pressure order three observed separately.', refinement='macro diameter', 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": "c74dff83", "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." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }