{ "cells": [ { "cell_type": "markdown", "id": "ca17c16d", "metadata": {}, "source": [ "# Mixed-local MHM: flux, pressure and boundary pressure\n", "\n", "We write every local block and the global boundary equation, then assemble, solve and evaluate the fields. No named Darcy constructor chooses the problem. The source is independently derived from $p=1+x^2+y^2$ and $q=-2(x,y)$: $f=\\nabla\\cdot q=-4$. Exterior pressure equals the analytical value.\n", "\n", "The example uses portable Basix RT0/P0 local kernels. RT0 contains this affine radial flux; its P0 pressure is an approximation and should not be declared exact. The general mixed-local construction follows [Durán, Devloo, Gomes and Valentin (2019)](https://doi.org/10.1016/j.cma.2019.05.013)." ] }, { "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/8f0c39ef7008398e2b53e90bc1b2d17075cbc14c8029634ed151fb89f11eca38/mixed_mhm_workflow-companion.zip\"\n", "COMPANION_SHA256 = \"8f0c39ef7008398e2b53e90bc1b2d17075cbc14c8029634ed151fb89f11eca38\"\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/mixed_mhm_workflow.ipynb\", directory=ROOT)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e8e39b39", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import sparse\n", "from pymhm import Equation, LocalContext, MeshHierarchy, SkeletonSpace, TriangleMesh\n", "from pymhm import assemble, bind_interface, bind_problem, solve\n", "from pymhm.fem.scalar.operators import rt0_operators, rt0_evaluate, face_integration, triangle_quadrature\n", "\n", "\n", "def exact_pressure(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Analytical nonhomogeneous boundary pressure.\"\"\"\n", " return 1 + np.sum(points**2, axis=1)\n", "\n", "\n", "def exact_flux(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Physical Darcy flux for identity permeability.\"\"\"\n", " return -2 * points\n", "\n", "\n", "def source(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Independently differentiated flux divergence.\"\"\"\n", " return np.full(len(points), -4.0)" ] }, { "cell_type": "markdown", "id": "7e934095", "metadata": {}, "source": [ "## 1. Define macro and local meshes\n", "\n", "The P0 normal-flux skeleton controls the global approximation. Each macrotriangle has its own red-refined RT0/P0 local mesh. Numbering and face incidence come from the mesh and shared pairing owner." ] }, { "cell_type": "code", "execution_count": null, "id": "d242efe2", "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)\n", "interface = bind_interface(skeleton, convention=\"normal\")" ] }, { "cell_type": "markdown", "id": "e2074074", "metadata": {}, "source": [ "## 2. Translate the mixed weak problem into blocks\n", "\n", "$$\n", "\\begin{aligned}\n", "(K^{-1}q,r)_T-(p,\\nabla\\cdot r)_T+\\langle\\eta,r\\cdot n_T\\rangle_{\\partial T}&=0,\\\\\n", "-(\\nabla\\cdot q,w)_T&=-(f,w)_T,\\\\\n", "\\langle q\\cdot n_T-\\lambda_T,\\xi\\rangle_{\\partial T}&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "The flux, pressure and private boundary-pressure coefficients form $u_T=(q_T,p_T,\\eta_T)$. The first two volume operators are $M$ and $D$; $N$ pairs boundary pressure with fine normal-flux coefficients. $F$ transports the skeletal normal-flux density to those integrated boundary coefficients.\n", "\n", "The common Basix/FEM operations integrate these forms in their declared coefficient basis; the block matrix and signs below are the user's formulation. The final trace map already includes canonical normal incidence, so `coordinates=\"global\"` explicitly avoids applying it again." ] }, { "cell_type": "code", "execution_count": null, "id": "f6017bd8", "metadata": {}, "outputs": [], "source": [ "def local_equations(local: LocalContext):\n", " \"\"\"Declare the three mixed equations and their independent global test.\"\"\"\n", " fine = local.mesh\n", " M, D, force = rt0_operators(fine, 1.0, source, order=6)\n", " _, F = face_integration(macro, local.cell, fine, skeleton)\n", " nq, npres, nb = len(fine.faces), len(fine.cells), len(fine.boundary_faces)\n", " N = sparse.coo_matrix((np.ones(nb), (fine.boundary_faces, np.arange(nb))),\n", " shape=(nq, nb)).tocsc()\n", " zero = sparse.csc_matrix((npres, nb))\n", " a = sparse.bmat([[M, -D.T, N], [-D, None, zero],\n", " [N.T, zero.T, None]], format=\"csc\")\n", " b = np.zeros((nq + npres + nb, F.shape[1]))\n", " b[nq + npres:] = -F\n", " L = np.r_[np.zeros(nq), -force, np.zeros(nb)]\n", " constant_shift = np.r_[np.zeros(nq), np.ones(npres + nb)][:, None]\n", " pressure_integral = np.r_[np.zeros(nq), fine.areas, np.zeros(nb)][:, None]\n", " return local.equations(a=a, L=L, b=b, c=-b.T,\n", " kernel=constant_shift, moments=pressure_integral,\n", " coordinates=\"global\", metadata=(fine, nq, npres))" ] }, { "cell_type": "markdown", "id": "973393b5", "metadata": {}, "source": [ "## 3. Define the global boundary equation\n", "\n", "The joint local kernel shifts pressure and auxiliary boundary pressure by the same constant. Its physical moment integrates only the volume pressure. The global pressure test uses $C=-B^T$; this mixed coefficient convention gives the **positive** Dirichlet boundary load. The primal pressure convention has a different load sign." ] }, { "cell_type": "code", "execution_count": null, "id": "a5a3cb40", "metadata": {}, "outputs": [], "source": [ "def global_equation(global_context):\n", " \"\"\"Declare weak exterior pressure in the mixed trace-test convention.\"\"\"\n", " boundary, fixed = global_context.boundary_data(exact_pressure, order=6)\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": "1a57bc8a", "metadata": {}, "source": [ "## 4. Assemble, solve and split the declared fields\n", "\n", "The local block layout was explicitly declared, so the result can be interpreted without guessing coefficient locations. Auxiliary boundary pressure is not included in the physical pressure norm." ] }, { "cell_type": "code", "execution_count": null, "id": "c20e163f", "metadata": {}, "outputs": [], "source": [ "system = assemble(problem)\n", "coefficients = solve(system)\n", "fields = []\n", "for metadata, values in zip(system.local_metadata, coefficients.fields, strict=True):\n", " fine, nq, npres = metadata\n", " fields.append((fine, values[:nq], values[nq:nq + npres]))\n", "assert coefficients.raw_residual < 1e-10\n", "print({\"original_equation_residual\": coefficients.raw_residual})" ] }, { "cell_type": "markdown", "id": "4843daba", "metadata": {}, "source": [ "## 5. Verify physical flux and fine-cell conservation\n", "\n", "Integrate each field on its own fine cells. The RT0 physical flux is evaluated with the shared Piola/basis owner. Fine-cell conservation is imposed by the pressure/divergence equation; it is stronger than merely checking each macroelement's total outflow." ] }, { "cell_type": "code", "execution_count": null, "id": "22ce1f42", "metadata": {}, "outputs": [], "source": [ "bary, weights = triangle_quadrature(8)\n", "pressure_squared, flux_squared, balance = 0.0, 0.0, []\n", "for fine, q, p in fields:\n", " points = np.einsum(\"qi,tia->tqa\", bary, fine.points[fine.cells])\n", " exact_p = exact_pressure(points.reshape(-1, 2)).reshape(points.shape[:2])\n", " exact_q = exact_flux(points.reshape(-1, 2)).reshape((*points.shape[:2], 2))\n", " q_values = rt0_evaluate(fine, q, bary)\n", " pressure_squared += fine.areas @ ((p[:, None] - exact_p)**2 @ weights)\n", " flux_squared += fine.areas @ (np.sum((q_values - exact_q)**2, axis=2) @ weights)\n", " _, divergence, force = rt0_operators(fine, 1.0, source, order=8)\n", " balance.append(np.max(abs(divergence @ q - force)))\n", "assert np.sqrt(flux_squared) < 1e-10\n", "assert max(balance) < 1e-10\n", "print({\"pressure_L2\": float(np.sqrt(pressure_squared)),\n", " \"physical_Darcy_flux_L2\": float(np.sqrt(flux_squared)),\n", " \"fine_cell_balance_max\": float(max(balance))})" ] }, { "cell_type": "markdown", "id": "428016dc", "metadata": {}, "source": [ "## 6. Plot analytical pressure, mixed pressure and the physical flux\n", "\n", "All panels include the actual macro mesh. Piecewise-constant pressure is displayed independently on each fine cell; no artificial interpolation makes it look continuous." ] }, { "cell_type": "code", "execution_count": null, "id": "c75a46ec", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), layout=\"constrained\")\n", "for index, title in enumerate((\"Analytical pressure\", \"RT0/P0 pressure\", \"Physical Darcy flux magnitude\")):\n", " for fine, q, p in fields:\n", " centers = fine.points[fine.cells].mean(axis=1)\n", " values = (exact_pressure(centers), p,\n", " np.linalg.norm(rt0_evaluate(fine, q, np.array([[1/3, 1/3, 1/3]]))[:, 0], axis=1))[index]\n", " artist = axes[index].tripcolor(*fine.points.T, fine.cells, facecolors=values,\n", " vmin=0 if index == 2 else 1,\n", " vmax=2*np.sqrt(2) if index == 2 else 3)\n", " for edge in macro.faces:\n", " axes[index].plot(*macro.points[edge].T, color=\"black\", lw=0.5)\n", " axes[index].set(title=title, xlabel=\"x\", ylabel=\"y\", aspect=\"equal\")\n", " fig.colorbar(artist, ax=axes[index])\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": "e022bbff", "metadata": {}, "source": [ "## 7. Verify a smooth mixed refinement\n", "\n", "The polynomial patch above verifies the flux, pressure convention and fine-cell conservation. Its exactly representable flux cannot give a meaningful convergence rate. We therefore use a separate smooth unit-square sequence:\n", "\n", "$$\n", "p=\\sin(2\\pi x)\\sin(2\\pi y),\\qquad\n", "f=8\\pi^2p,\\qquad q=-\\nabla p.\n", "$$\n", "\n", "Exterior pressure is zero. The three mixed saddle equations remain the ones written above; their local spaces are now RT1/P1 with P1 normal-flux traces and two fine subdivisions per macro edge. Mathematical RT1 corresponds to Basix degree two. Macro resolutions are $n=1,2,4,8,16,32$.\n", "\n", "For this smooth convex-square problem and fixed $h=H/2$, theorem 4, equations (53)–(55), of [Durán, Devloo, Gomes and Valentin (2019)](https://doi.org/10.1016/j.cma.2019.05.013) gives order two for physical flux and for the un-enriched P1 pressure. The enhanced pressure estimate requires the corresponding interior enrichment and does not apply to this pair.\n", "\n", "The recorded acquisition integrates both physical fields at independent quadrature orders 10 and 12, and checks original local/global equations, fine-cell divergence moments and oriented normal-flux moments. These controls differ from the error against the analytical fields. RT1/BDM and enriched families replace the compatible local operators; the declared global workflow is unchanged." ] }, { "cell_type": "markdown", "id": "1a0a8fd4", "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": "5b0fdb69", "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/mixed-mhm-current.json'\n", "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n", "refinement_rows = refinement_data['rows']\n", "refinement_error_keys = {'pressure L2': 'pressure_l2', 'physical Darcy flux L2': 'flux_l2'}\n", "refinement_targets = {'pressure L2': 2.0, '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='mixed-mhm', spaces='RT1/P1 locals; P1 normal-flux trace; two local edge subdivisions', rate_provenance='Durán, Devloo, Gomes and Valentin (2019), Theorem 4, equations (53)–(55): smooth compatible un-enriched RT1/P1 spaces on aligned shape-regular partitions with bounded commuting projections; pressure L2 estimates also require dual regularity.', 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": "66d8e4b3", "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" } }, "nbformat": 4, "nbformat_minor": 5 }