{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MH²M with multiscale permeability: pressure, conormal and trace\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", "The **three-field Multiscale Method (MH²M)** uses local pressure, local boundary conormal and a shared global pressure trace. We declare all three equations explicitly and let `LocalEquations` and `MultiscaleProblem` eliminate the local fields.\n", "\n", "Our two-dimensional P1/P0/P1 family follows [de Barros, Madureira and Valentin (2026, version 3)](https://arxiv.org/abs/2404.16978v3), equations (7), (28)–(29) and section 6.2. The oscillatory material and meshes below define an original introductory case, not a reproduction of the article's heterogeneous figures.\n", "\n", "All problem data, boundary conditions, provider, classical baseline, reconstruction, norms and plots are defined in notebook cells. PyMHM supplies only generic finite-element and algebraic operations. Open the downloaded notebook in a writable directory:\n", "\n", "```bash\n", "jupyter lab mh2m_multiscale.ipynb\n", "```\n", "\n", "Execution is serial; process workers need importable provider 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": "mh2m_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/35d2e7fcf53362a350d2f63714e035ebaac3c9388aeef7d0605955acc0154ea9/mh2m_multiscale-companion.zip\"\n", "COMPANION_SHA256 = \"35d2e7fcf53362a350d2f63714e035ebaac3c9388aeef7d0605955acc0154ea9\"\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/mh2m_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 import solve\n", "from pymhm.fem.traces.pressure_2d import PressureTraceSpace\n" ], "id": "mh2m_multiscale-01" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. The physical problem and the conormal sign\n", "\n", "Use the same scalar Darcy problem as the MsHHO introduction:\n", "\n", "$$\n", "\\begin{aligned}\n", " -\\nabla\\cdot(a_\\varepsilon\\nabla p)&=1 &&\\text{in }(0,1)^2,\\\\\n", " p&=0 &&\\text{on its boundary},\\\\\n", " a_\\varepsilon(x,y)&=2+\\sin(2\\pi x/\\varepsilon)\n", " \\sin(2\\pi y/\\varepsilon),\n", " &\\varepsilon&=1/8,\\\\\n", " q&=-a_\\varepsilon\\nabla p,\n", " &\\eta_K&=a_\\varepsilon\\nabla p\\cdot n_K=-q\\cdot n_K.\n", "\\end{aligned}\n", "$$\n", "\n", "The article's **conormal** $\\eta_K$ has the opposite sign to physical outward flux. Each macrocell has its private copy on a shared face; pressure trace $\\rho$ is global. The coefficient has ellipticity bounds 1 and 3.\n", "\n", "Our 32 macrotriangles do not resolve eight oscillations per direction. Each local mesh has 16 subdivisions per macroedge, keeping the material inside assembly and physical flux norms. This is fine local resolution on a simple macro mesh, not replacement of the coefficient by its average.\n" ], "id": "mh2m_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": "mh2m_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": "mh2m_multiscale-04" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Choose compatible spaces before assembling operators\n", "\n", "Use $k=0$ in the two-dimensional family of section 6.2:\n", "\n", "$$\n", "\\begin{aligned}\n", " \\Gamma_{H_\\Gamma}&=P_1^{\\mathrm{C0}}(\\text{pressure segments}),\\\\\n", " \\Lambda_{H_\\Lambda}(K)&=P_0^{\\mathrm{broken}}(\\text{conormal segments}),\\\\\n", " V_h(K)&=P_1^{\\mathrm{C0}}(\\text{local triangles}),\\\\\n", " H_\\Gamma&\\sim H/4,\\\\\n", " H_\\Lambda&\\sim H/8,\\\\\n", " h&\\sim H/16.\n", "\\end{aligned}\n", "$$\n", "\n", "Each pressure segment contains **two** conormal segments (M2), and each conormal segment contains **two** local boundary edges (M1). The uniform partitions match between incident macrocells. These mesh conditions accompany the chosen degrees; matrix dimensions alone do not establish compatibility.\n", "\n", "`PressureTraceSpace` shares pressure values at macrovertices and numbers interior face nodes. `SkeletonSpace` supplies the broken conormal basis, with a private copy per macrocell. Only `gamma` coordinates are global.\n" ], "id": "mh2m_multiscale-05" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "local_degree = 1\n", "gamma_segments, conormal_segments = 4, 8\n", "gamma = PressureTraceSpace.uniform(macro, degree=1, segments=gamma_segments)\n", "conormal = SkeletonSpace(\n", " macro,\n", " tuple(FaceSpace.uniform(0, conormal_segments) for _ in macro.faces),\n", ")\n", "assert conormal_segments == 2 * gamma_segments\n", "assert local_refinement == 2 * conormal_segments\n", "pressure_interface = bind_interface(gamma, convention=\"value\")\n", "conormal_interface = bind_interface(conormal, convention=\"value\")\n", "print(\n", " {\n", " \"global_pressure_trace_coordinates\": gamma.size,\n", " \"local_conormal_unknowns_per_macrocell\": 3 * conormal_segments,\n", " \"local_P1_nodes_per_macrocell\": (local_refinement + 1) * (local_refinement + 2) // 2,\n", " }\n", ")\n" ], "id": "mh2m_multiscale-06" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Write the three variational equations\n", "\n", "For each macrocell $K$,\n", "\n", "$$\n", "\\begin{aligned}\n", " (a_\\varepsilon\\nabla p_K,\\nabla v)_K\n", " -\\langle\\eta_K,v\\rangle_{\\partial K}&=(1,v)_K,\\\\\n", " -\\langle\\mu,p_K\\rangle_{\\partial K}\n", " +\\langle\\mu,\\rho\\rangle_{\\partial K}&=0,\\\\\n", " \\sum_K\\langle\\eta_K,\\xi\\rangle_{\\partial K}&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "The second equation matches **trace moments** against P0; it does not impose pointwise equality. The third equation balances conormal against continuous pressure trace tests.\n", "\n", "With $A_{ij}=a_K(\\phi_j,\\phi_i)$, $B_{ij}=\\int_{\\partial K}\\phi_i\\mu_j$ and $P_{j\\ell}=\\int_{\\partial K}\\mu_j\\psi_\\ell$,\n", "\n", "$$\n", "\\begin{aligned}\n", " \\begin{bmatrix}A&-B\\\\-B^T&0\\end{bmatrix}\n", " \\begin{bmatrix}p_K\\\\\\eta_K\\end{bmatrix}\n", " +\\begin{bmatrix}0\\\\P\\end{bmatrix}\\rho_K\n", " &=\\begin{bmatrix}f_K\\\\0\\end{bmatrix},\\\\\n", " \\sum_K P^T\\eta_K&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "Write the boundary pairing $B$ as `phi * v * ds` in UFL, selecting the private conormal space as its interface. `local.interface_pairing(conormal_interface)` integrates $P$ on the union of the pressure and conormal partitions and accumulates shared pressure vertices. Both are unsigned local boundary pairings: their mathematical signs are the explicit **$-B$ and $-B^T$** in the saddle. The context owns face ordering and coefficient transport.\n" ], "id": "mh2m_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", "`local.native_space` binds the declared fine mesh and equispaced element. The local saddle keeps its native pressure coordinates; `compile_form` assembles the expressions and releases native assembly resources. The same weak form can return portable nodal matrices for the independent classical baseline. The following section declares the two boundary spaces and their pairings.\n", "\n" ], "id": "mh2m_multiscale-08" }, { "cell_type": "code", "metadata": {}, "source": [ "import basix.ufl\n", "import ufl\n", "\n", "\n", "def user_volume_forms(\n", " fine: TriangleMesh,\n", " degree: int,\n", " *,\n", " binding: Any = None,\n", " portable: bool = True,\n", ") -> tuple[Any, Any, np.ndarray]:\n", " \"\"\"Assemble the user weak form in native or declared portable nodal coordinates.\"\"\"\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", " if portable:\n", " return A[order][:, order].tocsc(), M[order][:, order].tocsc(), F[order]\n", " return A, M, F\n" ], "execution_count": null, "outputs": [], "id": "mh2m_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": "mh2m_multiscale-10" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "def local_three_field_equations(local: LocalContext) -> LocalEquations:\n", " \"\"\"Declare native pressure/conormal equations and the shared pressure balance.\"\"\"\n", " fine = local.mesh\n", " element = basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " local_degree,\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " )\n", " pressure_space = local.native_space(element)\n", " stiffness, mass, force = user_volume_forms(\n", " fine,\n", " local_degree,\n", " binding=pressure_space,\n", " portable=False,\n", " )\n", " v = ufl.TestFunction(pressure_space.space)\n", " boundary_forms = local.trace_pairings(\n", " lambda phi, ds: phi * v * ds,\n", " interface=conormal_interface,\n", " )\n", " B = assemble_pairing(boundary_forms.forms)\n", " P = local.interface_pairing(conormal_interface, order=4)\n", "\n", " # These block dimensions express the three-field formulation itself.\n", " npressure, nconormal = stiffness.shape[0], B.shape[1]\n", " local_saddle = sparse.bmat([[stiffness, -B], [-B.T, None]], format=\"csc\")\n", " coupling = np.vstack((np.zeros((npressure, P.shape[1])), P))\n", " pressure_reconstruction = np.hstack((np.eye(npressure), np.zeros((npressure, nconormal))))\n", " local.field(\"pressure\", pressure_space, reconstruction=pressure_reconstruction)\n", " return local.equations(\n", " a=local_saddle,\n", " L=np.r_[force, np.zeros(nconormal)],\n", " b=coupling,\n", " c=coupling.T,\n", " metadata={\n", " \"mesh\": fine,\n", " \"A\": stiffness,\n", " \"F\": force,\n", " \"B\": B,\n", " \"P\": P,\n", " \"pressure_count\": npressure,\n", " \"volume_moments\": np.asarray(mass.sum(axis=1)).ravel(),\n", " },\n", " )\n" ], "id": "mh2m_multiscale-11" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Declare the shared trace equation and recover pressure\n", "\n", "The stiffness $A_K$ alone has a constant kernel. For the compatible spaces above, the complete local saddle is invertible: boundary pairings fix the constant and trace moments. We neither invent a kernel for that saddle nor pin an arbitrary pressure node.\n", "\n", "The local unknown of `LocalEquations` is $(p_K,\\eta_K)$. The coupling `b` is $[0;P]$, and the global balance `c` is its transpose. Generic assembly eliminates both local fields, leaving $\\rho$. Dirichlet fixes its exterior nodes, so no global pressure gauge remains.\n", "\n", "The article's Neumann maps use a **boundary-mean-zero** decomposition. The full saddle below retains the complete pressure and does not select any nodal or volume-mean complement. It expresses the same three-field equations directly. The field record below associates the pressure coefficients with their executed P1 basis; it contains no formulation or solver.\n" ], "id": "mh2m_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", " pressure_interface,\n", " local_three_field_equations,\n", " global_equation=Equation(0, 0),\n", " retained=0,\n", " fixed=lambda global_context: global_context.fix_faces(macro.boundary_faces, 0),\n", ")\n", "with threadpool_limits(1):\n", " system = assemble(problem)\n", " solution = solve(system)\n", "pressure_fields = solution.field(\"pressure\")\n", "physical_fields = tuple(\n", " ScalarField(field.mesh, local_degree, field.portable_coefficients) for field in pressure_fields\n", ")\n", "global_free_unknowns = gamma.size - len(problem.fixed)\n", "conormal_fields = tuple(\n", " mixed[data[\"pressure_count\"] :]\n", " for data, mixed in zip(system.local_metadata, solution.fields, strict=True)\n", ")\n", "print({\"global_pressure_trace_free\": global_free_unknowns, \"global_residual\": solution.residual})\n" ], "id": "mh2m_multiscale-14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5. Check injectivity, trace moments and macro conservation\n", "\n", "The constant test $v=1$ gives\n", "\n", "$$\n", " -\\int_{\\partial K}\\eta_K=\\int_K1=\\lvert K\\rvert.\n", "$$\n", "\n", "Thus physical outward flux $q\\cdot n_K=-\\eta_K$ conserves on each macrocell. Check this balance, $B^Tp_K=P\\rho_K$, the constant stiffness kernel and the original pressure equation. A rank check on $B$ confirms that no conormal mode is invisible on the first cell.\n", "\n", "These are execution diagnostics. Uniform stability depends on the compatible spaces and M1/M2 conditions, not merely a small algebraic residual. Raw $-a_\\varepsilon\\nabla p_h$ need not equal $-\\eta_K$ pointwise on a boundary and is not claimed to be an $H(\\mathrm{div})$ field on the fine mesh.\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": "mh2m_multiscale-15" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "diagnostics = {\n", " \"constant_kernel_relative\": 0.0,\n", " \"trace_moment_max\": 0.0,\n", " \"macro_conservation_max\": 0.0,\n", " \"original_local_relative_residual\": 0.0,\n", "}\n", "basis_digests = []\n", "archive = {\n", " \"macro_points\": macro.points,\n", " \"macro_cells\": macro.cells,\n", " \"pressure_trace\": solution.trace,\n", " \"gamma_nodes\": gamma.nodes,\n", "}\n", "for cell, (data, field, eta) in enumerate(\n", " zip(system.local_metadata, physical_fields, conormal_fields, strict=True)\n", "):\n", " A, F, B, P = data[\"A\"], data[\"F\"], data[\"B\"], data[\"P\"]\n", " executed_pressure = pressure_fields[cell].coefficients\n", " rho = 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[\"trace_moment_max\"] = max(\n", " diagnostics[\"trace_moment_max\"],\n", " float(np.max(abs(B.T @ executed_pressure - P @ rho))),\n", " )\n", " diagnostics[\"macro_conservation_max\"] = max(\n", " diagnostics[\"macro_conservation_max\"],\n", " float(abs(B.sum(axis=0) @ eta + F.sum())),\n", " )\n", " defect = A @ executed_pressure - B @ eta - F\n", " scale = max(np.linalg.norm(F), np.linalg.norm(A @ executed_pressure), np.linalg.norm(B @ eta))\n", " diagnostics[\"original_local_relative_residual\"] = max(\n", " diagnostics[\"original_local_relative_residual\"],\n", " float(np.linalg.norm(defect) / scale),\n", " )\n", " basis_digests.append(executed_array_digest(B, P))\n", " descriptor = pressure_fields[cell].definition.descriptor\n", " binding = solution.trace_bindings[cell]\n", " archive_cell(\n", " archive,\n", " cell,\n", " field,\n", " executed=pressure_fields[cell],\n", " trace_binding=binding,\n", " conormal=eta,\n", " boundary_pairing=B,\n", " pressure_trace_pairing=P,\n", " )\n", "assert np.linalg.matrix_rank(system.local_metadata[0][\"B\"]) == 3 * conormal_segments\n", "assert max(diagnostics.values()) < 1e-10\n", "print(diagnostics)\n" ], "id": "mh2m_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": "mh2m_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": "mh2m_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": "mh2m_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": "mh2m_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=\"MH²M\",\n", ")\n", "plt.show()\n" ], "id": "mh2m_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": "mh2m_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": "mh2m_multiscale-27" }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "record = {\n", " \"method\": \"MH²M\",\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", " \"mh2m_multiscale\",\n", " record,\n", " archive,\n", " {\n", " \"physical-fields\": figure,\n", " \"permeability\": fig_material,\n", " \"reference-refinement\": refinement_figure,\n", " },\n", " )\n", ")\n" ], "id": "mh2m_multiscale-28" }, { "cell_type": "markdown", "id": "f0e9e696", "metadata": {}, "source": [ "## 8. Verify a compatible smooth three-field family\n", "\n", "The heterogeneous teaching problem above is distinct from the smooth rate qualification. On the unit square with identity permeability, use the independently differentiated polynomial data\n", "\n", "$$\n", "\\begin{aligned}\n", "p(x,y)&=x(x-1)y(y-1),\\\\\n", "f(x,y)&=-2\\bigl[x(x-1)+y(y-1)\\bigr].\n", "\\end{aligned}\n", "$$\n", "\n", "Exterior pressure is zero. The three-field declarations remain the equations written above. Choose continuous Gamma P2 on one macroface segment, Lambda P1 on **two segments**, and local P2 pressure with **four subdivisions**. Two Lambda segments per Gamma segment and two fine edges per Lambda segment satisfy M2 and M1 with fixed ratio two. The macro sequence is $n=2,4,8,16,32$.\n", "\n", "The compatible estimate of [de Barros, Madureira and Valentin (2026)](https://arxiv.org/abs/2404.16978v3) gives second-order broken gradient for this smooth family. Third-order pressure is observed. The paper's section-8.1 experimental P2/P1/P2 family with one Lambda segment and two local subdivisions is discussed under Remark 16; it is a different discretization, and is not substituted for the M1/M2-compatible family below.\n", "\n", "Each current acquisition checks the original local equations, Lambda pressure pairing and free Gamma equations, and integrates the physical fields independently at quadrature orders 10 and 12." ] }, { "cell_type": "markdown", "id": "0b0189dd", "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": "828b3604", "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/mh2m-compatible-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', 'broken gradient L2': 'gradient_l2'}\n", "refinement_targets = {'broken gradient 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='mh2m', spaces='Gamma P2, one segment; Lambda P1, two segments; local P2/r4; M1=M2=2', rate_provenance='de Barros, Madureira and Valentin (2026), version 3, Theorem 19: energy order two under the two-dimensional M1/M2 Fortin conditions, regular aligned partitions and the stated pressure/source-lift/volume-flux regularity; pressure order three 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": "18ecb70e", "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-mh2m_multiscale", "metadata": {}, "source": [ "## References\n", "\n", "- Franklin de Barros, Alexandre L. Madureira, and Frédéric Valentin (2026). *A three-field Multiscale Method*. arXiv preprint, version 3, 5 August 2026; first submitted 25 April 2024. [arXiv: 2404.16978v3](https://arxiv.org/abs/2404.16978v3)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12" }, "pymhm": { "problem": "introduction", "methods": [ "MH2M" ], "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": "484d3227933298856a85dc24309b6aedf198358ea3f1faa23d50238dbb6ff2b8" } }, "nbformat": 4, "nbformat_minor": 5 }