{ "cells": [ { "cell_type": "markdown", "id": "b148f0aa", "metadata": {}, "source": [ "# Residual Petrov–Galerkin MHM: write the jump equation and enrichment\n", "\n", "We declare the primal local equations, add a user-written global residual form, and lift its boundary residual through the same local solver owner. The physical pressure is $p=\\sin(2\\pi x)\\sin(2\\pi y)$, with identity diffusion and independently derived source $8\\pi^2p$. Exterior pressure is zero.\n", "\n", "The construction follows [Fernando, Martins, Pereira and Valentin (2023)](https://doi.org/10.1007/s40314-023-02304-y). The trace multiplier of the base field differs from the enriched conservative normal flux. Neither volume gradient is automatically H(div)-conforming." ] }, { "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/87b99ef4019e62e7a343eae837aa8b0b526eba208077835a1e333f7340011e21/pgmhm_workflow-companion.zip\"\n", "COMPANION_SHA256 = \"87b99ef4019e62e7a343eae837aa8b0b526eba208077835a1e333f7340011e21\"\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/pgmhm_workflow.ipynb\", directory=ROOT)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d91fb504", "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)\n", "from scipy import sparse\n", "from pymhm import with_global_equation\n", "from pymhm.fem.scalar.triangle import scalar_operators, trace_coupling\n", "from pymhm.fem.traces.jump import face_jump_form\n", "from pymhm.postprocessing.nodal import nodal_field\n", "from pymhm.postprocessing.fields import DiscreteField\n", "from pymhm.core.reconstruction import reconstruct_local\n", "from pymhm.postprocessing.solutions import PGMHMSolution" ] }, { "cell_type": "markdown", "id": "f2f85d16", "metadata": {}, "source": [ "## 1. Define spaces\n", "\n", "Choose trace degree one and local degree three. In two dimensions this meets the energy-estimate restriction $k\\ge\\ell+d$ and $\\ell\\ge1$. Select $\\alpha=0.1$ and certified diffusion lower bound one." ] }, { "cell_type": "code", "execution_count": null, "id": "b46cf3f8", "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\")\n", "degree, alpha, lower = 3, 0.1, 1.0" ] }, { "cell_type": "markdown", "id": "4a9dc068", "metadata": {}, "source": [ "## 2. Define the primal volume and independent trace forms\n", "\n", "$$\n", "A_T(p,v)=(\\nabla p,\\nabla v)_T,\\qquad\n", "B_T(\\lambda,v)=\\langle\\lambda,v\\rangle_{\\partial T},\\qquad\n", "C_T(p,\\mu)=-\\langle p,\\mu\\rangle_{\\partial T}.\n", "$$\n", "\n", "For reference, the volume integrand is `ufl.inner(ufl.grad(p), ufl.grad(v)) * dx`. Here the shared Basix `scalar_operators` operation integrates this declared volume form in the portable nodal basis used by the common jump operator. It supplies element tabulation and quadrature; it does not select the multiscale method. The block signs, kernel, physical moment and global equation below are explicit. `trace_coupling` already applies canonical incidence, hence `coordinates=\"global\"`." ] }, { "cell_type": "code", "execution_count": null, "id": "b99ec598", "metadata": {}, "outputs": [], "source": [ "def local_equations(local: LocalContext):\n", " \"\"\"Declare volume diffusion, physical mean and signed scalar trace forms.\"\"\"\n", " fine = local.mesh\n", " A, mass, force = scalar_operators(fine, degree, diffusion=1.0, source=source, order=12)\n", " B = trace_coupling(macro, local.cell, fine, skeleton, degree)\n", " constant = np.ones((A.shape[0], 1))\n", " field = nodal_field(\"pressure\", fine, degree)\n", " local.field(\"pressure\", evaluator=field.evaluator,\n", " gradient_evaluator=field.gradient_evaluator, basis_id=field.basis_id)\n", " return local.equations(a=A, L=force, b=B, c=-B.T,\n", " kernel=constant, moments=mass @ constant,\n", " coordinates=\"global\", metadata=(fine,))\n", "\n", "\n", "def global_equation(global_context):\n", " \"\"\"Homogeneous weak Dirichlet pressure in the independent trace test.\"\"\"\n", " boundary, _ = global_context.boundary_data(0.0, order=8)\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)\n", "with threadpool_limits(1):\n", " base_system = assemble(problem)" ] }, { "cell_type": "markdown", "id": "4a90494d", "metadata": {}, "source": [ "## 3. Write the global residual jump equation\n", "\n", "Let $J_F z+j_F^f$ be the one-sided pressure jump reconstructed from reduced coordinates $z$, and set $\\tau_F=\\alpha/(2H_F)$. The additional global terms are\n", "\n", "$$\n", "\\begin{aligned}\n", "D_{\\mathrm{jump}}&=\\sum_F J_F^T W_F\\tau_FJ_F,\\\\\n", "g_{\\mathrm{jump}}&=\\sum_FJ_F^TW_F\\tau_F(p_D-j_F^f).\n", "\\end{aligned}\n", "$$\n", "\n", "`face_jump_form` supplies executed response maps, common quadrature and one-sided trace evaluation. The multiplication and global accumulation below implement the two displayed forms. Coincident global entries are summed by sparse COO assembly." ] }, { "cell_type": "code", "execution_count": null, "id": "c99f9e17", "metadata": {}, "outputs": [], "source": [ "penalties = tuple(face_jump_form(base_system, skeleton, face, degree, 12,\n", " alpha, lower, 0.0)\n", " for face in range(len(macro.faces)))\n", "rows, columns_, entries = [], [], []\n", "load = np.zeros_like(base_system.rhs)\n", "for data in penalties:\n", " weighted = data.weights * data.coefficient\n", " block = data.jump.T @ (weighted[:, None] * data.jump)\n", " forcing = data.jump.T @ (weighted * (data.prescribed - data.source_jump))\n", " rows.extend(np.repeat(data.indices, len(data.indices)))\n", " columns_.extend(np.tile(data.indices, len(data.indices)))\n", " entries.extend(block.ravel())\n", " np.add.at(load, data.indices, forcing)\n", "D_jump = sparse.coo_matrix((entries, (rows, columns_)),\n", " shape=base_system.matrix.shape).tocsc()\n", "system = with_global_equation(base_system, Equation(D_jump, load))\n", "with threadpool_limits(1):\n", " coefficients = solve(system)" ] }, { "cell_type": "markdown", "id": "13c6f8e7", "metadata": {}, "source": [ "## 4. Define and solve the enrichment problem\n", "\n", "The normal correction is $\\delta\\lambda_F=-\\tau_F([p]-p_D)$. Its local boundary load is assembled from the declared one-sided pairing, and the correction has zero physical volume mean:\n", "\n", "$$\n", "A_T\\delta p_T=-B_T\\delta\\lambda_T,\\qquad \\int_T\\delta p_T=0.\n", "$$\n", "\n", "The same `reconstruct_local` owner solves this additional local load. No factorization or basis algorithm is copied into the notebook. The final `PGMHMSolution` is a physical data record for norms and conservation; it does not construct or choose the problem." ] }, { "cell_type": "code", "execution_count": null, "id": "632a2c1d", "metadata": {}, "outputs": [], "source": [ "coordinates = np.r_[coefficients.trace, *coefficients.coarse]\n", "extra = [np.zeros(len(response.problem.load)) for response in system.responses]\n", "for data in penalties:\n", " jump = data.jump @ coordinates[data.indices] + data.source_jump - data.prescribed\n", " normal_correction = -data.coefficient * jump\n", " for cell, sign, evaluation in data.sides:\n", " extra[cell] += sign * (evaluation.T @ (data.weights * normal_correction))\n", "enriched = []\n", "for response, base, boundary_load in zip(system.responses, coefficients.fields, extra, strict=True):\n", " local_problem = response.problem.with_load(-boundary_load)\n", " correction = reconstruct_local(local_problem, np.zeros(local_problem.coupling.shape[1]),\n", " np.zeros(local_problem.coarse_basis.shape[1]))\n", " enriched.append(base + correction)\n", "solution = PGMHMSolution(skeleton, local_meshes, coefficients.fields, tuple(enriched),\n", " coefficients, system, 1.0, source, degree, 12,\n", " alpha, lower, penalties, tuple(extra))\n", "pressure_fields = tuple(DiscreteField(nodal_field(\"pressure\", fine, degree), values)\n", " for fine, values in zip(local_meshes, enriched, strict=True))\n", "print({\"pressure_L2\": solution.l2_error(exact_pressure, order=12, enriched=True),\n", " \"Darcy_flux_L2\": solution.flux_l2_error(lambda x: -exact_gradient(x), order=12, enriched=True),\n", " \"enriched_macro_balance\": float(np.max(abs(solution.conservation_residuals())))})\n", "assert np.max(abs(solution.conservation_residuals())) < 1e-10" ] }, { "cell_type": "markdown", "id": "b8290585", "metadata": {}, "source": [ "## 5. Plot the enriched physical pressure\n", "\n", "The macro mesh is shown on every analytical, numerical and error panel." ] }, { "cell_type": "code", "execution_count": null, "id": "4e784874", "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": "f8d6cbe4", "metadata": {}, "source": [ "## 6. Reach the smooth enriched-flux regime\n", "\n", "The smooth unit-square sequence uses the same declared residual-jump equations and enrichment, with $p=\\sin(2\\pi x)\\sin(2\\pi y)$, independently derived $f=8\\pi^2p$, identity permeability and homogeneous exterior pressure. Its spaces are P3 local pressure on one triangle, P1 traces and $\\alpha=0.1$, on all macro resolutions $n=2,4,8,16,32,64$.\n", "\n", "The energy estimates in section 5 of [Fernando et al. (2023)](https://doi.org/10.1007/s40314-023-02304-y) separate $h^k$ and $H^{\\ell+1}$, with $\\ell\\geq1$ and $k\\geq\\ell+d$. Here $\\ell=1$, $k=3$, $d=2$; the second-order term governs the enriched physical flux. Pressure order three is observed. The reconstructed vector field is the enriched volume gradient, and does not become a globally H(div)-conforming flux merely because macro balance holds.\n", "\n", "The acquisition checks the original enriched local nodal equations, the global jump equation and the separate macro balance. Independent order-10 and order-12 quadratures control the full pressure and physical vector-flux errors. Read all six levels below, including their coarse behavior." ] }, { "cell_type": "markdown", "id": "6693a804", "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": "18927ec0", "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/tutorial-methods/pgmhm-current.json'\n", "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n", "refinement_rows = refinement_data['rows']\n", "refinement_error_keys = {'enriched pressure L2': 'pressure_l2', 'enriched Darcy flux L2': 'flux_l2'}\n", "refinement_targets = {'enriched 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='pgmhm', spaces='P3/r1 locals; P1 trace; residual enrichment; alpha=0.1', rate_provenance='Fernando et al. (2023), section 5: the energy estimate requires ell>=1, k>=ell+d, the stated regularity and sufficiently small alpha for uniform stability; energy order two and pressure order three are measured for the declared finite P3/P1 family with alpha=0.1.', 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": "9b2e9aa4", "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 }