{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MsHHO with multiscale permeability: from the problem to the API\n", "\n", "\n", "The main workflow is **meshes → spaces → local equations → global balance → assemble → solve → fields and errors**. `MeshHierarchy` associates macro and local meshes; `bind_interface` and `LocalContext` own supported numbering, geometric orientation and native coordinate conversions. The mathematical forms, physical trace meaning, local modes and gauge remain explicit in the cells below. Fully manual/custom spaces use the same numerical owners; see the [custom-interface notebook](https://github.com/ipes-lncc/pymhm/blob/main/notebooks/foundations/operators/custom_interface.ipynb).\n", "\n", "We build a simple **Multiscale Hybrid High-Order** Darcy problem and explain every local and global block. All case definitions, equations, boundary conditions, reference solves, norms and plots are visible in this notebook. The weak operator is written in UFL. PyMHM provides generic compilation, trace pairings, constrained reconstruction and linear algebra.\n", "\n", "The method family follows [Chaumont-Frelet, Ern, Lemaire and Valentin (2022)](https://doi.org/10.1051/m2an/2021082), sections 4–5. Their local spaces are defined by exactly solved PDEs. We use a finite local Galerkin realization on declared fine meshes. This original introductory problem is not a reproduction of an article figure or table.\n", "\n", "Run with the installed package in a writable directory and the compatible native DOLFINx/UFL backend:\n", "\n", "```bash\n", "jupyter lab mshho_multiscale.ipynb\n", "```\n", "\n", "Prerequisites: a weak Poisson formulation, basic NumPy and Jupyter. Cell-defined providers are used in serial execution; process workers require importable callables.\n", "\n", "Field evaluation, norms, plots and executed-array archives use the importable\n", "[supporting scalar helpers](https://github.com/ipes-lncc/pymhm/blob/main/examples/introduction/scalar.py).\n", "The physical data and local/global variational equations remain explicit below.\n", "\n", "\n", "The PyMHM distribution contains only the library. The first cell explicitly\n", "downloads a checksum-verified companion archive and prepares the declared data\n", "using its local Python helpers. You can inspect these support files in `ROOT`.\n", "Complete studies and historical replay retain their existing opt-in flags.\n" ], "id": "mshho_multiscale-00" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "from pathlib import Path\n", "import os\n", "import sys\n", "from pymhm.io.workspace import workspace_from_archive\n", "\n", "# Download verified support files; this operation does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/c6b3e4e96bcbd504442e8fb9d7fc8f9474603b5609d123d79225f62d97167e2c/mshho_multiscale-companion.zip\"\n", "COMPANION_SHA256 = \"c6b3e4e96bcbd504442e8fb9d7fc8f9474603b5609d123d79225f62d97167e2c\"\n", "WORKSPACE = Path(\n", " os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd() / \".pymhm-companions\" / COMPANION_SHA256)\n", ")\n", "ROOT = workspace_from_archive(COMPANION_URL, sha256=COMPANION_SHA256, directory=WORKSPACE)\n", "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n", "if str(ROOT) not in sys.path:\n", " sys.path.insert(0, str(ROOT))\n", "\n", "# Explicitly prepare declared data with the downloaded Python helpers.\n", "from scripts.notebook_reproduction import notebook_workspace\n", "\n", "ROOT = notebook_workspace(\"introduction/mshho_multiscale.ipynb\", directory=ROOT)\n", "root = ROOT\n", "print(\"Workspace:\", ROOT)\n", "\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from typing import Any\n", "from scipy import sparse\n", "from threadpoolctl import threadpool_limits\n", "from pymhm import Equation, LocalEquations, assemble\n", "from pymhm import MeshHierarchy, LocalContext, bind_interface, bind_problem\n", "from pymhm.core.equations import compile_form\n", "from pymhm.fem.traces.interval import FaceSpace, SkeletonSpace\n", "from examples.introduction.scalar import (\n", " plot_material,\n", " plot_reference_differences,\n", " plot_scalar_comparison,\n", " executed_array_digest,\n", " ScalarField,\n", " archive_cell,\n", " compare_scalar_fields,\n", " dirichlet_solve,\n", " native_scalar_space,\n", " save_scalar_report,\n", ")\n", "from pymhm.fem.scalar.triangle import nodal_space\n", "from pymhm.meshes.triangle import TriangleMesh\n", "from pymhm.backends.forms import assemble_pairing\n", "\n", "plt.rcParams.update({\"figure.dpi\": 110, \"font.size\": 10})\n", "from pymhm.core.moments import energy_reconstruction\n" ], "id": "mshho_multiscale-01" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Define the physical problem, material period and meshes\n", "\n", "On $\\Omega=(0,1)^2$, solve\n", "\n", "$$\n", "\\begin{aligned}\n", " -\\nabla\\cdot(a_\\varepsilon\\nabla p)&=1 &&\\text{in }\\Omega,\\\\\n", " p&=0 &&\\text{on }\\partial\\Omega,\\\\\n", " a_\\varepsilon(x,y)&=2+\n", " \\sin(2\\pi x/\\varepsilon)\\sin(2\\pi y/\\varepsilon),\n", " &\\varepsilon&=1/8,\\\\\n", " q&=-a_\\varepsilon\\nabla p.\n", "\\end{aligned}\n", "$$\n", "\n", "The bounds $1\\le a_\\varepsilon\\le3$ give uniform ellipticity. The coefficient period $\\varepsilon=1/8$ is physical data: eight material oscillations fit in each direction. The macro and local meshes are discretization choices. Our **macro mesh** has four squares per direction, split into 32 triangles; each macrocell contains material variation. Its independent local Galerkin problem uses 16 subdivisions per macroedge. The coarse unknowns are cell and face moments. Cell moments are eliminated locally, leaving a global system of face moments; the fine local fields resolve the coefficient variation.\n", "\n", "`TriangleMesh` defines geometry. Below, we write and execute the form $\\int_K a_\\varepsilon\\nabla u\\cdot\\nabla v$ and load $\\int_K v$ in UFL before introducing optional assembly conveniences.\n" ], "id": "mshho_multiscale-02" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "epsilon = 1 / 8\n", "\n", "\n", "def permeability(points: np.ndarray) -> np.ndarray:\n", " \"\"\"Evaluate aε I through its positive scalar factor, with ε=1/8.\"\"\"\n", " phase = 2 * np.pi * points / epsilon\n", " return 2 + np.sin(phase[:, 0]) * np.sin(phase[:, 1])\n", "\n", "\n", "macro = TriangleMesh.unit_square(4)\n", "source = 1.0\n", "local_refinement = 16\n", "assembly_order = 8\n", "print(\n", " {\n", " \"macro_triangles\": len(macro.cells),\n", " \"epsilon\": epsilon,\n", " \"H_max\": float(macro.lengths.max()),\n", " \"h_max\": float(macro.lengths.max() / local_refinement),\n", " \"ellipticity_bounds\": [1.0, 3.0],\n", " \"source\": source,\n", " }\n", ")\n" ], "id": "mshho_multiscale-03" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "fig_material = plot_material(macro, permeability, resolution=64, limits=(1, 3))\n", "plt.show()\n" ], "id": "mshho_multiscale-04" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. The unknowns are pressure moments\n", "\n", "Choose $m=k_F=1$: three P1 volume moments per macrotriangle and two P1 moments per macroface. We store **integrals** rather than polynomial projection coefficients. The two coordinate systems are related by invertible polynomial mass matrices:\n", "\n", "$$\n", "\\begin{aligned}\n", " z_{K,i}&=\\int_K p_h r_{K,i},\n", " &r_K&=(1,\\widehat x,\\widehat y),\\\\\n", " z_{F,j}&=\\int_F p_h\\psi_{F,j},\n", " &\\psi_F&=(1,2t-1).\n", "\\end{aligned}\n", "$$\n", "\n", "Here $\\widehat x,\\widehat y$ are Cartesian coordinates scaled to the macrocell bounding box. They span P1 on each triangle. Both incident macrocells share $z_F$; their complete fine nodal fields remain independent. Exterior face moments are zero by Dirichlet data.\n", "\n", "`SkeletonSpace` numbers these face coordinates. Its name does not determine their physical meaning: **these are pressure moments**, not flux coefficients. Our local finite space is continuous P2; the moment constraints must be independent in that space.\n" ], "id": "mshho_multiscale-05" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "local_degree = 2\n", "skeleton = SkeletonSpace(macro, tuple(FaceSpace.uniform(1) for _ in macro.faces))\n", "fixed = {int(dof): 0.0 for face in macro.boundary_faces for dof in skeleton.dofs(int(face))}\n", "global_free_unknowns = skeleton.size - len(fixed)\n", "print(\n", " {\n", " \"face_moments_total\": skeleton.size,\n", " \"face_moments_free\": global_free_unknowns,\n", " \"cell_moments_per_macrocell\": 3,\n", " \"local_polynomial_degree\": local_degree,\n", " }\n", ")\n" ], "id": "mshho_multiscale-06" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Construct the local basis by constrained energy minimization\n", "\n", "Let $A_K$ be the fine stiffness and $C_K$ the matrix of volume and face functionals. The reconstruction $R_K$ gives the minimum-energy field with prescribed moments:\n", "\n", "$$\n", "\\begin{aligned}\n", " \\begin{bmatrix}A_K&C_K\\\\C_K^T&0\\end{bmatrix}\n", " \\begin{bmatrix}R_K\\\\M_K\\end{bmatrix}\n", " &=\\begin{bmatrix}0\\\\I\\end{bmatrix},\\\\\n", " C_K^TR_K&=I,\n", " &E_K&=R_K^TA_KR_K.\n", "\\end{aligned}\n", "$$\n", "\n", "The Neumann stiffness has the constant kernel. Independent moments fix it; the stiffness must also be coercive on the moment-zero space. A small residual alone does not prove these properties.\n", "\n", "The P1 cell tests are represented exactly at P2 nodes, so `mass @ cell_tests` integrates their volume moments. The face functionals are written directly as `phi * v * ds` with UFL. The value interface binding owns their geometry and numbering: **pressure** moments have no normal-incidence sign. The volume coefficients use the native space binding's executed nodal convention.\n", "\n", "`energy_reconstruction` solves the generic $A/C$ saddle above. It receives no physical-model or method selector. The constant source belongs to P1: in integral-moment coordinates the projected-source functional is simply $z_{K,0}$, hence load $(1,0,0)$ on cell coordinates and zero on face coordinates.\n" ], "id": "mshho_multiscale-07" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write the local differential operator in UFL\n", "\n", "Define the trial pressure `p` and test pressure `v` on the local Pk space. The UFL expressions below are the executable weak form, mass pairing and source. Changing this expression defines a new local operator; the multiscale contracts do not select a physical model.\n", "\n", "$$\n", "a_K(p,v)=\\int_K a_\\varepsilon\\nabla p\\cdot\\nabla v,\n", "\\qquad \\ell_K(v)=\\int_K v.\n", "$$\n", "\n", "The small geometry adapter keeps the same fine triangles and equispaced nodal basis. It maps native coefficients by their physical coordinates; it contains no PDE. `compile_form` assembles the expressions and releases native assembly resources. The interface pairing is introduced separately, with its declared normal convention.\n", "\n" ], "id": "mshho_multiscale-08" }, { "cell_type": "code", "metadata": {}, "source": [ "import basix.ufl\n", "import ufl\n", "from pymhm.backends.forms import assemble_pairing\n", "\n", "\n", "def user_volume_forms(\n", " fine: TriangleMesh, degree: int, *, binding: Any = None\n", ") -> tuple[Any, Any, np.ndarray]:\n", " \"\"\"Execute the user-written weak operator, mass and source in the stated basis.\"\"\"\n", " if binding is None:\n", " _, V, order = native_scalar_space(fine, degree)\n", " else:\n", " V, order = binding.space, binding.mapping\n", " p, v = ufl.TrialFunction(V), ufl.TestFunction(V)\n", " x = ufl.SpatialCoordinate(V.mesh)\n", " a_epsilon = 2 + ufl.sin(2 * np.pi * x[0] / epsilon) * ufl.sin(2 * np.pi * x[1] / epsilon)\n", " dx = ufl.Measure(\"dx\", domain=V.mesh, metadata={\"quadrature_degree\": 2 * assembly_order})\n", " a = a_epsilon * ufl.inner(ufl.grad(p), ufl.grad(v)) * dx\n", " mass = p * v * dx\n", " L = source * v * dx\n", " A, M, F = compile_form(a), compile_form(mass), compile_form(L)\n", " return A[order][:, order].tocsc(), M[order][:, order].tocsc(), F[order]\n" ], "execution_count": null, "outputs": [], "id": "mshho_multiscale-09" }, { "cell_type": "code", "metadata": {}, "source": [ "first_fine = macro.submesh(0, local_refinement)\n", "A_ufl, M_ufl, F_ufl = user_volume_forms(first_fine, local_degree)\n", "print(\n", " {\n", " \"operator\": \"user-written UFL\",\n", " \"local_nodal_unknowns\": A_ufl.shape[0],\n", " \"constant_kernel_relative\": float(\n", " np.linalg.norm(A_ufl @ np.ones(A_ufl.shape[0])) / sparse.linalg.norm(A_ufl)\n", " ),\n", " }\n", ")\n", "form_equivalence = {} # Filled by the optional convenience check at the end.\n" ], "execution_count": null, "outputs": [], "id": "mshho_multiscale-10" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "def local_moment_equations(local: LocalContext) -> LocalEquations:\n", " \"\"\"Declare MsHHO P1 cell/face moments on a local conforming P2 space.\"\"\"\n", " cell, fine = local.cell, local.mesh\n", " binding = local.native_space(\n", " basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " local_degree,\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " )\n", " )\n", " stiffness, mass, force = user_volume_forms(fine, local_degree, binding=binding)\n", " _, nodes = nodal_space(fine, local_degree)\n", " vertices = macro.points[macro.cells[cell]]\n", " scaled_xy = 2 * (nodes - vertices.min(axis=0)) / np.ptp(vertices, axis=0) - 1\n", " cell_tests = np.column_stack((np.ones(len(nodes)), scaled_xy))\n", " volume_moments = mass @ cell_tests\n", " v = ufl.TestFunction(binding.space)\n", " face_functionals = local.trace_pairings(lambda phi, ds: phi * v * ds)\n", " face_moments = assemble_pairing(face_functionals.forms, axis=\"columns\")[binding.mapping]\n", " moments = np.column_stack((volume_moments, face_moments))\n", " reconstruction, energy = energy_reconstruction(stiffness, moments)\n", "\n", " # The declared source is constant. Since the first cell moment is\n", " # integral_K p, its load coordinate is source; all other coordinates vanish.\n", " moment_load = np.r_[source, np.zeros(moments.shape[1] - 1)]\n", " ncell = volume_moments.shape[1]\n", " return local.equations(\n", " a=energy[:ncell, :ncell],\n", " L=moment_load[:ncell],\n", " b=energy[:ncell, ncell:],\n", " c=energy[ncell:, :ncell],\n", " coordinates=\"global\",\n", " d=energy[ncell:, ncell:],\n", " g=moment_load[ncell:],\n", " metadata={\n", " \"mesh\": fine,\n", " \"R\": reconstruction,\n", " \"C\": moments,\n", " \"A\": stiffness,\n", " \"F\": force,\n", " \"energy\": energy,\n", " },\n", " )\n" ], "id": "mshho_multiscale-11" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. The global problem contains face moments\n", "\n", "Partition $E_K$ into cell ($c$) and face ($F$) blocks. The projected-source formulation becomes\n", "\n", "$$\n", "\\begin{aligned}\n", " E_{cc}z_K+E_{cF}z_{\\partial K}&=(1,0,0)^T,\\\\\n", " \\sum_K\\big(E_{Fc}z_K+E_{FF}z_{\\partial K}\\big)&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "These are exactly `LocalEquations`: `a,L,b` define the first equation; `c,d,g` contribute the second. `assemble` eliminates the three cell moments and adds shared face contributions.\n", "\n", "`Equation(0, zeros)` supplies no additional global term. `retained=0` retains no extra modes: the fine stiffness kernel has already been fixed in the constrained reconstruction, and cell moments are ordinary eliminated local unknowns. Prescribed exterior pressure moments remove the global pressure gauge.\n", "\n", "After solving, `solution.fields` contains **cell moments**. Multiplying the complete local moment vector by the executed $R_K$ recovers physical nodal pressure. The small field record below only associates coefficients with their declared mesh and Pk basis; it is not a solver.\n" ], "id": "mshho_multiscale-12" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "hierarchy = MeshHierarchy(\n", " macro, tuple(macro.submesh(cell, local_refinement) for cell in range(len(macro.cells)))\n", ")\n", "problem = bind_problem(\n", " hierarchy,\n", " bind_interface(skeleton, convention=\"value\"),\n", " local_moment_equations,\n", " global_equation=Equation(0, 0),\n", " retained=0,\n", " fixed=fixed,\n", ")\n", "with threadpool_limits(1):\n", " system = assemble(problem)\n", " solution = system.solve()\n", "\n", "# solution.fields contains cell moments, not physical nodal pressures.\n", "physical_fields = tuple(\n", " ScalarField(\n", " data[\"mesh\"], local_degree, data[\"R\"] @ np.r_[cell_moments, solution.local_trace(cell)]\n", " )\n", " for cell, (data, cell_moments) in enumerate(\n", " zip(system.local_metadata, solution.fields, strict=True)\n", " )\n", ")\n", "print(\n", " {\n", " \"global_free_unknowns\": global_free_unknowns,\n", " \"global_residual\": solution.residual,\n", " \"local_nodal_unknowns_per_macrocell\": len(physical_fields[0].values),\n", " }\n", ")\n" ], "id": "mshho_multiscale-14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5. Check the kernel and the executed moment basis\n", "\n", "Check $A_K1=0$, $C_K^TR_K=I$, and the moments of the recovered pressure. Full pressure may jump across a macroface; its P1 face moments agree.\n", "\n", "The finite local Galerkin realization does not automatically inherit the ideal article space's $H(\\mathrm{div})$ flux property. We report the raw physical flux $-a_\\varepsilon\\nabla p_h$ inside each fine element and do not assert fine-cell conservation.\n", "\n", "Archive each executed $R_K$ with its SHA256 digest and the complete moment vector. Field replay uses that same matrix; persisted moment coefficients must not be applied to a numerically changed basis.\n", "\n", "The physical reconstruction reads `solution.local_trace(cell)`: the bound space supplies the executed local face coefficients, including any declared basis changes. The value convention used here has no outward-normal sign. The method's reconstruction matrix remains an explicit part of its mathematical definition.\n" ], "id": "mshho_multiscale-15" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "diagnostics = {\n", " \"constant_kernel_relative\": 0.0,\n", " \"moment_identity_max\": 0.0,\n", " \"recovered_moment_max\": 0.0,\n", "}\n", "basis_digests = []\n", "archive = {\"macro_points\": macro.points, \"macro_cells\": macro.cells, \"face_moments\": solution.trace}\n", "for cell, (data, field, cell_moments) in enumerate(\n", " zip(system.local_metadata, physical_fields, solution.fields, strict=True)\n", "):\n", " A, C, R = data[\"A\"], data[\"C\"], data[\"R\"]\n", " target = np.r_[cell_moments, solution.local_trace(cell)]\n", " diagnostics[\"constant_kernel_relative\"] = max(\n", " diagnostics[\"constant_kernel_relative\"],\n", " float(np.linalg.norm(A @ np.ones(A.shape[0])) / np.linalg.norm(A.data)),\n", " )\n", " diagnostics[\"moment_identity_max\"] = max(\n", " diagnostics[\"moment_identity_max\"],\n", " float(np.max(abs(C.T @ R - np.eye(C.shape[1])))),\n", " )\n", " diagnostics[\"recovered_moment_max\"] = max(\n", " diagnostics[\"recovered_moment_max\"],\n", " float(np.max(abs(C.T @ field.values - target))),\n", " )\n", " basis_digests.append(executed_array_digest(R))\n", " archive_cell(archive, cell, field, reconstruction=R, moments=C, moment_coordinates=target)\n", "assert max(diagnostics.values()) < 1e-10\n", "print(diagnostics)\n" ], "id": "mshho_multiscale-16" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6. A classical primal Galerkin baseline, with its own refinement check\n", "\n", "On a fine **globally conforming** mesh, the baseline solves the same operator, material, source and boundary conditions:\n", "\n", "$$\n", "\\begin{aligned}\n", " p_{\\mathrm{CG}}&\\in V_{\\mathrm{CG}}\\subset H_0^1(\\Omega),\\\\\n", " (a_\\varepsilon\\nabla p_{\\mathrm{CG}},\\nabla v)_\\Omega\n", " &=(1,v)_\\Omega,\\qquad v\\in V_{\\mathrm{CG}}.\n", "\\end{aligned}\n", "$$\n", "\n", "Use P2 on 32, 64 and 128 squares per direction, each split into two triangles. Assembly and boundary identification are visible below; `solve_linear` owns the free algebraic solve. This classical global assembly is separate from multiscale condensation. Both use DOLFINx/UFL with the same weak form; their global assembly procedures and boundary spaces are separate.\n", "\n", "The finest classical field is a **numerical reference**, not an exact solution. First compare 32→64 and 64→128 in physical norms. The plotted difference rate is a reference-resolution indicator; two consecutive differences do not constitute a rigorous reference-error bound.\n" ], "id": "mshho_multiscale-17" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "references = []\n", "reference_dimensions = []\n", "for resolution in (32, 64, 128):\n", " fine_global = TriangleMesh.unit_square(resolution)\n", " A_ref, _, F_ref = user_volume_forms(fine_global, 2)\n", " _, xy = nodal_space(fine_global, 2)\n", " boundary = np.flatnonzero(np.any(np.isclose(xy, 0) | np.isclose(xy, 1), axis=1))\n", " p_ref = dirichlet_solve(A_ref, F_ref, boundary, np.zeros(len(boundary)))\n", " references.append(ScalarField(fine_global, 2, p_ref))\n", " reference_dimensions.append(\n", " {\n", " \"squares_per_direction\": resolution,\n", " \"triangles\": len(fine_global.cells),\n", " \"free_unknowns\": len(p_ref) - len(boundary),\n", " }\n", " )\n", " print(reference_dimensions[-1])\n", "reference = references[-1]\n" ], "id": "mshho_multiscale-19" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "reference_refinement = [\n", " compare_scalar_fields((coarser,), finer, permeability, order=8)\n", " for coarser, finer in zip(references[:-1], references[1:])\n", "]\n", "for interval, norms in zip((\"32 → 64\", \"64 → 128\"), reference_refinement):\n", " print(interval, norms)\n", "\n", "# These rates concern differences between successive reference meshes.\n", "# They are indicators, not exact-solution error rates.\n", "refinement_figure = plot_reference_differences(\n", " [1 / 64, 1 / 128],\n", " reference_refinement,\n", ")\n", "\n", "plt.show()\n", "\n", "# A reference must resolve its own physical fields before serving as a baseline.\n", "assert reference_refinement[-1][\"pressure_relative\"] < reference_refinement[0][\"pressure_relative\"]\n", "assert reference_refinement[-1][\"flux_relative\"] < reference_refinement[0][\"flux_relative\"]\n" ], "id": "mshho_multiscale-21" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 7. Compare pressure and physical flux on a common fine partition\n", "\n", "The reference mesh resolves every local mesh interface in this nested, uniformly refined case. Integrate on its triangles and use its field norms as denominators:\n", "\n", "$$\n", "\\begin{aligned}\n", " e_p&=\\frac{\\lVert p_h-p_{\\mathrm{ref}}\\rVert_{L^2(\\Omega)}}\n", " {\\lVert p_{\\mathrm{ref}}\\rVert_{L^2(\\Omega)}},\\\\\n", " e_q&=\\frac{\\lVert q_h-q_{\\mathrm{ref}}\\rVert_{L^2(\\Omega)}}\n", " {\\lVert q_{\\mathrm{ref}}\\rVert_{L^2(\\Omega)}},\\\\\n", " e_E&=\\frac{\\lVert a_\\varepsilon^{-1/2}(q_h-q_{\\mathrm{ref}})\\rVert_{L^2(\\Omega)}}\n", " {\\lVert a_\\varepsilon^{-1/2}q_{\\mathrm{ref}}\\rVert_{L^2(\\Omega)}}.\n", "\\end{aligned}\n", "$$\n", "\n", "The comparison code keeps the coefficient inside flux and energy norms. Quadrature orders 8 and 10 independently check norm integration. At volume quadrature points, evaluate the original coefficients directly; do not compare images or interpolated color maps.\n", "\n", "For plotting, sample every fine triangle separately. Duplicate interface points retain separate incident values. Pressure and flux-magnitude panels share scales with their reference panels; each difference has its own colorbar. All panels show the actual macro mesh.\n" ], "id": "mshho_multiscale-22" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "errors_q8 = compare_scalar_fields(physical_fields, reference, permeability, order=8)\n", "errors_q10 = compare_scalar_fields(physical_fields, reference, permeability, order=10)\n", "quadrature_difference = max(\n", " abs(errors_q8[key] - errors_q10[key]) / max(abs(errors_q10[key]), np.finfo(float).tiny)\n", " for key in (\"pressure_L2\", \"flux_L2\", \"flux_energy\")\n", ")\n", "print(\n", " {\"physical_relative_errors\": errors_q10, \"relative_norm_change_q8_q10\": quadrature_difference}\n", ")\n", "assert quadrature_difference < 1e-5\n", "\n", "figure = plot_scalar_comparison(\n", " macro,\n", " physical_fields,\n", " reference,\n", " permeability,\n", " numerical_label=\"MsHHO\",\n", ")\n", "plt.show()\n" ], "id": "mshho_multiscale-24" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Optional: use a prepared operator for this weak form\n", "\n", "For the scalar diffusion form written above, `scalar_operators` provides a\n", "portable Basix assembly convenience. Defining a new physical operator does not\n", "require this function: the main computation already assembles the UFL expression\n", "written by the user. The following check verifies agreement in the same local\n", "finite-element space and basis before substituting this convenience.\n" ], "id": "mshho_multiscale-26" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "from pymhm.fem.scalar.triangle import scalar_operators\n", "\n", "A_prepared, M_prepared, F_prepared = scalar_operators(\n", " first_fine,\n", " local_degree,\n", " diffusion=permeability,\n", " source=source,\n", " order=assembly_order,\n", ")\n", "form_equivalence = {\n", " \"stiffness_relative\": float(\n", " sparse.linalg.norm(A_ufl - A_prepared) / sparse.linalg.norm(A_prepared)\n", " ),\n", " \"mass_relative\": float(sparse.linalg.norm(M_ufl - M_prepared) / sparse.linalg.norm(M_prepared)),\n", " \"source_relative\": float(np.linalg.norm(F_ufl - F_prepared) / np.linalg.norm(F_prepared)),\n", "}\n", "assert max(form_equivalence.values()) < 1e-10\n", "print(form_equivalence)\n" ], "id": "mshho_multiscale-27" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "record = {\n", " \"method\": \"MsHHO\",\n", " \"physical_problem\": f\"-div(aε grad(p))={source:g}, p|boundary=0\",\n", " \"source\": float(source),\n", " \"epsilon\": epsilon,\n", " \"macro_triangles\": len(macro.cells),\n", " \"local_refinement\": local_refinement,\n", " \"local_degree\": local_degree,\n", " \"assembly_order\": assembly_order,\n", " \"error_orders\": [8, 10],\n", " \"global_free_unknowns\": global_free_unknowns,\n", " \"reference_dimensions\": reference_dimensions,\n", " \"reference_refinement\": reference_refinement,\n", " \"physical_errors\": errors_q10,\n", " \"quadrature_relative_difference\": quadrature_difference,\n", " \"global_residual\": solution.residual,\n", " \"diagnostics\": diagnostics,\n", " \"basis_digests\": basis_digests,\n", " \"form_equivalence\": form_equivalence,\n", "}\n", "print(\n", " save_scalar_report(\n", " ROOT,\n", " \"mshho_multiscale\",\n", " record,\n", " archive,\n", " {\n", " \"physical-fields\": figure,\n", " \"permeability\": fig_material,\n", " \"reference-refinement\": refinement_figure,\n", " },\n", " )\n", ")\n" ], "id": "mshho_multiscale-28" }, { "cell_type": "markdown", "id": "ef98bf8b", "metadata": {}, "source": [ "## 8. Separate smooth source and energy convergence\n", "\n", "The oscillatory teaching problem above is not replaced by a smooth one. To qualify the asymptotic rate independently, use the identity coefficient on the unit square and\n", "\n", "$$\n", "p=\\sin(\\pi x)\\sin(\\pi y),\\qquad\n", "f=2\\pi^2p,\\qquad q=-\\nabla p.\n", "$$\n", "\n", "Exterior pressure is zero. Keep $m=\\ell=0$ cell/face moments, the projected source, and P3 local energy reconstruction with two subdivisions. The macro sequence is $n=1,2,4,8,16,32$.\n", "\n", "Theorem 6.3 of [Chaumont-Frelet, Ern, Lemaire and Valentin (2022)](https://doi.org/10.1051/m2an/2021082) separates the $H^{m+2}$ source term and $H^{\\ell+1}$ flux term, under the stated source/flux regularity assumptions. The energy target here is order one. Pressure order two is a measured observation. The ideal method uses exact local reconstructions; this executed finite Galerkin realization independently controls its reconstruction and moment equations rather than identifying P3/r2 with an exact local solve.\n", "\n", "The current norm-only record retains order-10/order-12 independent physical error integration, original cell/free-face moment equations and reconstruction moment constraints. Their algebraic moment residuals are distinct from the L2 pressure/flux errors plotted below." ] }, { "cell_type": "markdown", "id": "dbdf3dad", "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": "64e90f7a", "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/mshho-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 = {'physical Darcy flux L2': 1.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='mshho', spaces='P0 cell and face moments; P3/r2 energy reconstruction; projected source', rate_provenance='Chaumont-Frelet et al. (2022), Theorem 6.3: first-order energy for exact local PDE reconstructions, source in H^1 and A grad(p) in H^1 on regular matching macrocell submeshes; finite P3/r2 reconstruction adds local approximation error; pressure order two observed separately.', 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": "ce058b07", "metadata": {}, "source": [ "The upper panel preserves all coarse and fine measurements. Dashed lines show declared target powers anchored at the finest measured error. The shaded interval always contains the final four levels; it does not select points to improve a fitted slope. Read the consecutive orders together with the target-normalized errors and the stated quadrature/local-equation controls. A slope alone does not establish the hypotheses of an error estimate." ] }, { "cell_type": "markdown", "id": "references-mshho_multiscale", "metadata": {}, "source": [ "## References\n", "\n", "- Théophile Chaumont-Frelet, Alexandre Ern, Simon Lemaire, and Frédéric Valentin (2022). *Bridging the Multiscale Hybrid-Mixed and Multiscale Hybrid High-Order Methods*, ESAIM: M2AN 56, 261–285. [DOI: 10.1051/m2an/2021082](https://doi.org/10.1051/m2an/2021082)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12" }, "pymhm": { "problem": "introduction", "methods": [ "MsHHO" ], "source_outputs": "empty" }, "pymhm_editorial_revision": { "date": "2026-10-06", "scope": "Markdown citations and bibliography only", "preserved": "Code cells, execution counts, execution timestamps and numerical outputs", "preserved_code_cells_sha256": "b933d084b0b56067c96ab6cf9443f234d4c1d1f1418e7bf2c41fcc46ea26b89a" } }, "nbformat": 4, "nbformat_minor": 5 }