{ "cells": [ { "cell_type": "markdown", "id": "b0faffd8", "metadata": {}, "source": [ "# Multiscale Darcy: a two-level primal MHM workflow\n", "\n", "The original MHM is hybrid because it couples broken local pressure fields\n", "through a skeleton multiplier. Its global problem is mixed: the unknowns\n", "are normal Darcy-flux coefficients and one pressure mean per macroelement.\n", "The local finite element method here is primal. Local Neumann responses\n", "contain the full permeability field; no homogenized permeability is\n", "substituted. See [Harder, Paredes and Valentin (2013)](https://doi.org/10.1016/j.jcp.2013.03.019).\n", "\n", "This tutorial uses **two-level MHM**. The first level discretizes the\n", "skeletal normal flux and retains one pressure mean per macroelement.\n", "The second level approximates the local Neumann responses in fine-mesh\n", "$Q_2$ spaces. With exact local responses, the skeletal discretization\n", "defines the one-level method; see\n", "[Araya et al. (2013)](https://doi.org/10.1137/120888223).\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", "This notebook builds primal MHM step by step: bind the macro and local meshes, declare the spaces, write local UFL equations and the global balance, then assemble and solve. A small macro mesh controls the global problem; independent local fine meshes resolve the material oscillations. The baseline is classical conforming Galerkin assembled on several much finer global meshes.\n", "\n", "Install `pymhm[notebooks,visualization]` and the compatible native DOLFINx/UFL backend, then open this notebook with `jupyter lab` in a writable directory. Run the cells in order. Physical data and equations are declared here; downloaded companion helpers provide evaluation, norms, plots and archives.\n", "\n", "The manufactured solution distinguishes the finite-reference error from the multiscale error. This is an analytical control devised for this tutorial, rather than a matched paper reproduction.\n", "\n", "The local Neumann responses, skeletal flux unknowns and retained cell constants follow [Harder, Paredes and Valentin (2013)](https://doi.org/10.1016/j.jcp.2013.03.019). The oscillatory manufactured data and refinement sequence are defined for this tutorial.\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" ] }, { "cell_type": "code", "execution_count": null, "id": "e1b1f91c", "metadata": {}, "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/3473ec3bc17a7c7dc8fdf16c6b2dd1057bdbbd6ddbddc6e6cf4b2c9b83b2f027/darcy_multiscale_convergence-companion.zip\"\n", "COMPANION_SHA256 = \"3473ec3bc17a7c7dc8fdf16c6b2dd1057bdbbd6ddbddc6e6cf4b2c9b83b2f027\"\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/darcy_multiscale_convergence.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 dataclasses import dataclass\n", "from functools import partial\n", "from typing import Any\n", "from pymhm import Equation, LocalEquations, assemble, columns\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", " Array,\n", " dirichlet_solve,\n", " evaluate_bound_pressure,\n", " evaluate_qk,\n", " grid_quadrature,\n", " native_scalar_space,\n", " observed_rates,\n", " physical_errors,\n", " plot_convergence,\n", " plot_field_panels,\n", " rectangle_panel,\n", ")\n", "import ufl\n", "from pymhm import ExecutionConfig, CartesianMacroMesh\n", "from pymhm.fem.scalar.quadrilateral import qk_space, quadrilateral_operators\n", "from pymhm.fem.scalar.operators import boundary_data\n" ] }, { "cell_type": "markdown", "id": "8d906eb6", "metadata": {}, "source": [ "## 1. Physical problem and an independently derived source\n", "\n", "On $\\Omega=(0,1)^2$, Darcy's law and conservation are\n", "\n", "$$\n", "\\boldsymbol q=-K\\nabla p,\\qquad\\nabla\\cdot\\boldsymbol q=f,\n", "\\qquad p=0\\quad\\text{on }\\partial\\Omega.\n", "$$\n", "\n", "Choose $\\varepsilon=1/8$, $a=1.5$ and\n", "\n", "$$\n", "\\begin{aligned}\n", "K(x,y)&=\\mathrm{diag}\\!\\left(e^{a\\sin(\\omega x)},e^{a\\cos(\\omega y)}\\right),\n", "\\qquad\\omega=2\\pi/\\varepsilon,\\\\\n", "p_\\star(x,y)&=\\sin(\\pi x)\\sin(\\pi y).\n", "\\end{aligned}\n", "$$\n", "\n", "The uniform spectral contrast is $e^{2a}\\simeq20.1$. Eight oscillations fit in each direction, even when only four macroelements fit in that direction. Differentiate **before** discretizing:\n", "\n", "$$\n", "f=\\pi^2(K_{xx}+K_{yy})p_\\star\n", " -(\\partial_xK_{xx})\\partial_xp_\\star\n", " -(\\partial_yK_{yy})\\partial_yp_\\star.\n", "$$\n", "\n", "The data class below implements these expressions directly. Its source never applies an assembled matrix to the exact coefficients. The coefficient is smooth, but its derivatives depend on $\\varepsilon$; our convergence measurements keep $\\varepsilon$ fixed." ] }, { "cell_type": "code", "execution_count": null, "id": "44b39512", "metadata": {}, "outputs": [], "source": [ "@dataclass(frozen=True)\n", "class OscillatoryDarcyData:\n", " \"\"\"Smooth anisotropic diffusion with declared wavelength and contrast.\n", "\n", " ``Kxx=exp(a*sin(2*pi*x/epsilon))`` and\n", " ``Kyy=exp(a*cos(2*pi*y/epsilon))``; cross entries vanish. Thus the\n", " uniform spectral contrast is ``exp(2*a)``, independently of the mesh.\n", " \"\"\"\n", "\n", " epsilon: float = 0.125\n", " amplitude: float = 1.5\n", "\n", " def permeability(self, points: Array) -> Array:\n", " \"\"\"Return the physical SPD tensor at XY points, with no averaging.\"\"\"\n", " phase = 2 * np.pi * points / self.epsilon\n", " tensors = np.zeros((len(points), 2, 2))\n", " tensors[:, 0, 0] = np.exp(self.amplitude * np.sin(phase[:, 0]))\n", " tensors[:, 1, 1] = np.exp(self.amplitude * np.cos(phase[:, 1]))\n", " return tensors\n", "\n", " def pressure(self, points: Array) -> Array:\n", " \"\"\"Return ``sin(pi*x)*sin(pi*y)`` on the unit square.\"\"\"\n", " return np.sin(np.pi * points[:, 0]) * np.sin(np.pi * points[:, 1])\n", "\n", " def gradient(self, points: Array) -> Array:\n", " \"\"\"Return the independently differentiated exact pressure gradient.\"\"\"\n", " x, y = np.pi * points.T\n", " return np.pi * np.column_stack((np.cos(x) * np.sin(y), np.sin(x) * np.cos(y)))\n", "\n", " def flux(self, points: Array) -> Array:\n", " \"\"\"Return exact physical Darcy flux ``-K grad(p)``.\"\"\"\n", " return -np.einsum(\"qab,qb->qa\", self.permeability(points), self.gradient(points))\n", "\n", " def source(self, points: Array) -> Array:\n", " \"\"\"Return independently differentiated ``-div(K grad(p))``.\n", "\n", " This expression uses analytic derivatives of both K and p, never the\n", " assembled operator or a discretely manufactured right-hand side.\n", " \"\"\"\n", " omega = 2 * np.pi / self.epsilon\n", " k = self.permeability(points)\n", " derivative_x = self.amplitude * omega * np.cos(omega * points[:, 0]) * k[:, 0, 0]\n", " derivative_y = -self.amplitude * omega * np.sin(omega * points[:, 1]) * k[:, 1, 1]\n", " gradient = self.gradient(points)\n", " return (\n", " np.pi**2 * (k[:, 0, 0] + k[:, 1, 1]) * self.pressure(points)\n", " - derivative_x * gradient[:, 0]\n", " - derivative_y * gradient[:, 1]\n", " )\n" ] }, { "cell_type": "markdown", "id": "choose-scales", "metadata": {}, "source": [ "## 2. Choose the macro mesh, local meshes and trace space\n", "\n", "Choose the discretization of the global and local problems:\n", "\n", "- **Macro mesh:** `macro` partitions the unit square into macroelements. Here $H$ denotes the side length of each square.\n", "- **Local meshes and spaces:** `local_refinement=r` divides each macroelement into $r\\times r$ fine squares, with side length $h=H/r$. `local_degree=k` sets the local pressure space $Q_k$.\n", "- **Global trace space:** the trace degree $\\ell$ and `trace_segments` select the normal-flux approximation on each macroface. Segmentation refines this global space.\n", "\n", "These choices configure the two-level discretization. The coefficient period $\\varepsilon=1/8$ is a physical parameter of the permeability and remains fixed throughout the refinement study.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "39766b49", "metadata": {}, "outputs": [], "source": [ "data = OscillatoryDarcyData(epsilon=1 / 8, amplitude=1.5)\n", "exact = lambda points: (data.pressure(points), data.gradient(points))\n", "macro = CartesianMacroMesh(4, 4)\n", "# H=1/4; h_local=H/r=1/64; degree k=2; trace degree ell=1.\n", "local_degree, local_refinement, trace_segments = 2, 16, 2\n", "# The oscillatory, nonpolynomial material requires independently qualified integration.\n", "assembly_order = 26 # Gauss points per axis; UFL quadrature degree 51.\n", "skeleton = SkeletonSpace(\n", " macro, tuple(FaceSpace.uniform(1, trace_segments, continuous=True) for _ in macro.faces)\n", ")\n", "print(\n", " {\n", " \"macro_cells\": len(macro.cells),\n", " \"H\": 1 / 4,\n", " \"local_h\": 1 / 64,\n", " \"epsilon\": data.epsilon,\n", " \"spectral_contrast\": float(np.exp(2 * data.amplitude)),\n", " }\n", ")\n" ] }, { "cell_type": "markdown", "id": "77647269", "metadata": {}, "source": [ "## 3. Translate the weak formulation into local and global objects\n", "\n", "On each macroelement $T$, pressure uses a continuous $Q_k$ fine space. The skeletal multiplier $\\lambda$ is the physical normal Darcy flux in the unique orientation $\\boldsymbol n_F$ of face $F$. The normal interface binding derives $s_{TF}=\\boldsymbol n_T\\cdot\\boldsymbol n_F$ from the macroface topology. Write the trace integrals in UFL; `LocalContext` supplies their basis supports and coefficient maps.\n", "\n", "$$\n", "\\begin{aligned}\n", "a_T(p_T,v)+\\langle s_{TF}\\lambda,v\\rangle_{\\partial T}&=(f,v)_T,\\\\\n", "a_T(p_T,v)&=\\int_T K\\nabla p_T\\cdot\\nabla v.\n", "\\end{aligned}\n", "$$\n", "\n", "| Mathematical term | Code declaration |\n", "| --- | --- |\n", "| Local energy and source | `a=inner(K*grad(p),grad(v))*dx`, `L=f*v*dx` in UFL |\n", "| Oriented normal-flux pairing | `local.trace_pairings(lambda phi, ds: phi*v*ds)` |\n", "| Constant kernel and physical average | `kernel`, `moments` in `LocalEquations` |\n", "| Weak pressure continuity | `local.trace_pairings(lambda phi, ds: -phi*p*ds, axis=\"rows\")` |\n", "| Macroface coordinates | `bind_interface(skeleton, convention=\"normal\")` supplies local/global maps |\n", "\n", "The local Neumann kernel is constant. The pressure decomposition is\n", "\n", "$$\n", "p_T=c_T+T_T\\lambda+\\widehat T_Tf,\n", "\\qquad\\int_TT_T\\lambda=\\int_T\\widehat T_Tf=0.\n", "$$\n", "\n", "The source response and each interface response are solved on the\n", "mean-zero complement. The constant $c_T$ is a global physical unknown.\n", "The retained compatibility equation supplies the missing constant test,\n", "\n", "$$\n", "\\int_{\\partial T}s_{TF}\\lambda\\,ds=\\int_Tf\\,dx.\n", "$$\n", "\n", "Together with the pressure-coupling moments, this is the original MHM\n", "global problem. It is not imposed by setting all local pressure means\n", "to zero.\n", "\n", "Interior faces have zero pressure-jump moments; Dirichlet faces have prescribed pressure moments. With the declared sign $C=-B^T$, the additional global right-hand side is\n", "\n", "$$\n", "-\\sum_T\\langle s_{TF}\\mu,p_T\\rangle_{\\partial T}\n", "=-\\langle\\mu,p_D\\rangle_{\\Gamma_D}.\n", "$$\n", "\n", "`Equation(0, ...)` stores this additional global term. One constant mode per macroelement gives the compatibility equation imposing **macro conservation**. The Darcy flux is evaluated as $q_h=-K\\nabla p_h$.\n", "\n", "All local solves use 26 Gauss points per coordinate direction, declared in UFL\n", "as `quadrature_degree=51`. This integrates the oscillatory coefficient and\n", "manufactured source on the coarsest local mesh as well as the finer meshes.\n", "Assembly quadrature is an independent numerical choice: increasing it changes\n", "neither the local polynomial degree nor the skeleton space. Field errors use a\n", "separate, fixed integration grid, so changes in the discretization choices\n", "can be compared consistently.\n", "\n", "The provider below declares the executable UFL weak forms and each coupling sign. `local.native_space` owns the native geometry and coefficient convention, while `local.trace_pairings` binds both interface forms independently. The provider supplies no face DOF indices or nodal permutation. The generic `assemble` operation owns elimination, assembly and reconstruction; no PDE-specific problem constructor is used." ] }, { "cell_type": "code", "execution_count": null, "id": "c7aa44aa", "metadata": {}, "outputs": [], "source": [ "def ufl_coefficient_source(domain: Any, physical_data: OscillatoryDarcyData) -> tuple[Any, Any]:\n", " \"\"\"Declare the same smooth K and independently differentiated manufactured f in UFL.\"\"\"\n", " x = ufl.SpatialCoordinate(domain)\n", " a, omega = physical_data.amplitude, 2 * np.pi / physical_data.epsilon\n", " kxx = ufl.exp(a * ufl.sin(omega * x[0]))\n", " kyy = ufl.exp(a * ufl.cos(omega * x[1]))\n", " tensor = ufl.as_matrix(((kxx, 0.0), (0.0, kyy)))\n", " exact_pressure = ufl.sin(np.pi * x[0]) * ufl.sin(np.pi * x[1])\n", " pressure_x = np.pi * ufl.cos(np.pi * x[0]) * ufl.sin(np.pi * x[1])\n", " pressure_y = np.pi * ufl.sin(np.pi * x[0]) * ufl.cos(np.pi * x[1])\n", " kxx_x = a * omega * ufl.cos(omega * x[0]) * kxx\n", " kyy_y = -a * omega * ufl.sin(omega * x[1]) * kyy\n", " force = np.pi**2 * (kxx + kyy) * exact_pressure - kxx_x * pressure_x - kyy_y * pressure_y\n", " return tensor, force\n", "\n", "\n", "@dataclass(frozen=True)\n", "class DarcyLocalProvider:\n", " \"\"\"Declare primal Qk UFL energy, physical flux coupling and pressure moments.\"\"\"\n", "\n", " macro: CartesianMacroMesh\n", " skeleton: SkeletonSpace\n", " data: object\n", " degree: int\n", " refinement: int\n", " quadrature_order: int\n", "\n", " def __call__(self, local: LocalContext) -> LocalEquations:\n", " \"\"\"Declare K grad(p)·grad(v), f v, oriented normal flux and the constant average.\"\"\"\n", " cell, fine = local.cell, local.mesh\n", " import basix\n", " import basix.ufl\n", "\n", " element = basix.ufl.element(\n", " \"Lagrange\",\n", " \"quadrilateral\",\n", " self.degree,\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " )\n", " binding = local.native_space(element)\n", " domain, space = binding.mesh, binding.space\n", " mapping = binding.mapping\n", " p, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n", " K, f = ufl_coefficient_source(domain, self.data)\n", " dx = ufl.Measure(\n", " \"dx\", domain=domain, metadata={\"quadrature_degree\": 2 * self.quadrature_order - 1}\n", " )\n", " # These executed UFL forms are the local weak formulation, directly.\n", " a = ufl.inner(K * ufl.grad(p), ufl.grad(v)) * dx\n", " load = f * v * dx\n", " local.field(\"pressure\", binding)\n", " b = local.trace_pairings(lambda phi, ds: phi * v * ds)\n", " c = local.trace_pairings(lambda phi, ds: -phi * p * ds, axis=\"rows\")\n", " area = float(fine.areas.sum())\n", " return local.equations(\n", " a=a,\n", " L=load,\n", " b=b,\n", " c=c,\n", " kernel=np.ones((len(mapping), 1)),\n", " moments=columns((v / area) * dx),\n", " metadata=(fine, mapping),\n", " )\n" ] }, { "cell_type": "markdown", "id": "b608ce79", "metadata": {}, "source": [ "## 4. Construct, assemble and solve the global problem\n", "\n", "`retained=1` retains one pressure-average coordinate per macroelement; `bind_problem` derives its global coordinate layout. The global equation imposes pressure through boundary moments, without fixing local pressure nodes. Dirichlet data remove the global constant-pressure ambiguity, so no additional physical gauge is required.\n", "\n", "The local problems compute source responses and responses to trace basis functions. Their coefficients are recovered after solving for the macroface and retained coordinates. They do not create a large conforming global mesh." ] }, { "cell_type": "markdown", "id": "3ef6eca6", "metadata": {}, "source": [ "### Evaluate pressure in its executed coefficient basis\n", "\n", "Before reconstructing fields, define how the x-fastest Qk coefficients represent pressure and physical gradients. The second function selects a macrocell and evaluates its independent local coefficients. This is field evaluation, not another solve; it preserves the distinction between broken pressure and physical Darcy flux." ] }, { "cell_type": "code", "execution_count": null, "id": "132619bc", "metadata": {}, "outputs": [], "source": [ "provider = DarcyLocalProvider(macro, skeleton, data, local_degree, local_refinement, assembly_order)\n", "boundary, fixed = boundary_data(skeleton, data.pressure, order=7)\n", "hierarchy = MeshHierarchy(\n", " macro, tuple(macro.submesh(cell, local_refinement) for cell in range(len(macro.cells)))\n", ")\n", "interface = bind_interface(skeleton, convention=\"normal\")\n", "problem = bind_problem(\n", " hierarchy,\n", " interface,\n", " provider,\n", " global_equation=lambda global_problem: Equation(0, global_problem.trace_load(-boundary)),\n", " retained=1,\n", " fixed=fixed,\n", ")\n", "system = assemble(problem, execution=ExecutionConfig(\"serial\", native_threads=1))\n", "solution = system.solve()\n", "pressure_fields = solution.field(\"pressure\")\n", "local_meshes = tuple(field.mesh for field in pressure_fields)\n", "mhm_fields = tuple(field.portable_coefficients for field in pressure_fields)\n", "mhm_evaluator = partial(evaluate_bound_pressure, macro, pressure_fields)\n", "print(\n", " {\n", " \"global_unknowns\": system.matrix.shape[0],\n", " \"largest_local_unknowns\": max(len(v) for v in solution.fields),\n", " \"original_equations_relative_residual\": solution.raw_residual,\n", " }\n", ")\n", "# Named fields carry their mesh and executed basis; no index map is needed to evaluate.\n", "first_point = macro.points[macro.cells[0]].mean(axis=0, keepdims=True)\n", "print(\"First macrocell pressure at its center:\", pressure_fields[0].evaluate(first_point))\n" ] }, { "cell_type": "markdown", "id": "7b0dc99c", "metadata": {}, "source": [ "### Optional numerical equivalence after the UFL definition\n", "\n", "The main formulation above is the executed UFL weak form. Once its operator and coordinate conventions are explicit, `quadrilateral_operators` provides an alternative assembled representation. The following control verifies equivalence on one local mesh; it does not select the global/local problem or hide its forms. The classical baseline below remains an explicit UFL assembly." ] }, { "cell_type": "code", "execution_count": null, "id": "6f23b2e4", "metadata": {}, "outputs": [], "source": [ "declared = problem.local_provider(0)\n", "fine, mapping = declared.metadata\n", "ufl_a = compile_form(declared.a)\n", "ufl_load = compile_form(declared.L, (len(mapping),))\n", "ready_a, _, ready_load = quadrilateral_operators(\n", " fine,\n", " local_degree,\n", " permeability=data.permeability,\n", " source=data.source,\n", " order=provider.quadrature_order,\n", ")\n", "operator_defect = np.max(np.abs((ufl_a[mapping][:, mapping] - ready_a).data), initial=0.0)\n", "load_defect = float(np.max(np.abs(ufl_load[mapping] - ready_load), initial=0.0))\n", "print({\"UFL_vs_ready_operator_max\": float(operator_defect), \"UFL_vs_ready_source_max\": load_defect})\n", "assert operator_defect < 1e-10 * max(1.0, float(np.max(np.abs(ready_a.data))))\n", "assert load_defect < 1e-10 * max(1.0, float(np.max(np.abs(ready_load))))\n" ] }, { "cell_type": "markdown", "id": "a4bf75da", "metadata": {}, "source": [ "## 5. Classical primal Galerkin on independently refined global meshes\n", "\n", "Now assemble **one global continuous matrix**, without an MHM skeleton or local condensation. The operator, permeability, source and boundary values are identical. Strong Dirichlet values eliminate boundary nodes; `dirichlet_solve`, defined below, performs only that generic algebraic elimination.\n", "\n", "$$\n", "A^{\\rm CG}_{ij}=\\int_\\Omega K\\nabla\\phi_j\\cdot\\nabla\\phi_i,\n", "\\qquad F^{\\rm CG}_i=\\int_\\Omega f\\phi_i.\n", "$$\n", "\n", "Evaluate pressure and raw gradients directly from the executed nodal basis. Integrate on a common grid resolving every compared finite-element mesh, with physical Gauss weights. The field measurements are pressure $L^2$, vector Darcy flux $L^2$, and flux energy\n", "\n", "$$\n", "\\left(\\int_\\Omega\\boldsymbol e_q^T K^{-1}\\boldsymbol e_q\\right)^{1/2}.\n", "$$\n", "\n", "The baseline is assembled independently as a different classical discretization, while sharing PyMHM basis and integration kernels. It is not an independent external-code comparison. Successive references are compared, and the analytical solution separately verifies their errors." ] }, { "cell_type": "code", "execution_count": null, "id": "ab98496c", "metadata": {}, "outputs": [], "source": [ "reference_rows, reference_evaluators = [], []\n", "error_points, error_weights = grid_quadrature(macro.bounds, (128, 128), order=7)\n", "for n in (32, 64, 128):\n", " fine = CartesianMacroMesh(n, n)\n", " domain, space, mapping = native_scalar_space(fine, 2)\n", " p, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n", " K, f = ufl_coefficient_source(domain, data)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 2 * assembly_order - 1})\n", " a_cg = compile_form(ufl.inner(K * ufl.grad(p), ufl.grad(v)) * dx)\n", " load_cg = compile_form(f * v * dx)\n", " _, nodes = qk_space(fine, 2)\n", " exterior = np.flatnonzero(np.any(np.isclose(nodes, 0) | np.isclose(nodes, 1), axis=1))\n", " coefficients = dirichlet_solve(\n", " a_cg, load_cg, mapping[exterior], data.pressure(nodes[exterior])\n", " )[mapping]\n", " evaluator = partial(evaluate_qk, fine, 2, coefficients)\n", " reference_evaluators.append(evaluator)\n", " errors = physical_errors(evaluator, exact, data.permeability, error_points, error_weights)\n", " reference_rows.append({\"n\": n, \"unknowns\": len(nodes), **errors})\n", "reference_rows\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3733ae08", "metadata": {}, "outputs": [], "source": [ "reference_refinement = [\n", " physical_errors(first, second, data.permeability, error_points, error_weights)\n", " for first, second in zip(reference_evaluators[:-1], reference_evaluators[1:])\n", "]\n", "print(\"Successive reference differences:\", reference_refinement)\n", "reference = reference_evaluators[-1]\n", "comparison = physical_errors(\n", " mhm_evaluator, reference, data.permeability, error_points, error_weights\n", ")\n", "exact_errors = physical_errors(mhm_evaluator, exact, data.permeability, error_points, error_weights)\n", "print(\"MHM versus fine CG:\", comparison)\n", "print(\"MHM versus exact solution:\", exact_errors)\n", "# The exact solution keeps this conclusion independent of reference uncertainty.\n", "assert reference_rows[-1][\"flux_L2\"] < reference_rows[0][\"flux_L2\"]\n", "assert reference_rows[-1][\"pressure_L2\"] < reference_rows[0][\"pressure_L2\"]\n" ] }, { "cell_type": "markdown", "id": "262d63b7", "metadata": {}, "source": [ "## 6. Plot the coefficient, pressure and physical flux\n", "\n", "The overlay is the actual macro mesh used by MHM, on every analytical, numerical and reference panel. Duplicate interface coordinates belong to separate local triangulations, preserving both one-sided values. Neither the broken pressure nor its raw flux is smoothed across macrofaces.\n", "\n", "The permeability tensor has two distinct diagonal components: $K_{xx}$ controls transport in the x direction and $K_{yy}$ in the y direction. Plot both with the same physical color limits to compare their oscillations. The field panels then place the analytical pressure and Darcy flux beside MHM and the refined classical reference." ] }, { "cell_type": "code", "execution_count": null, "id": "ec0d74f2", "metadata": {}, "outputs": [], "source": [ "permeability_panels = {\n", " \"Permeability $K_{xx}$\": rectangle_panel(\n", " macro, exact, quantity=\"permeability_xx\", permeability=data.permeability\n", " ),\n", " \"Permeability $K_{yy}$\": rectangle_panel(\n", " macro, exact, quantity=\"permeability_yy\", permeability=data.permeability\n", " ),\n", "}\n", "permeability_bounds = (np.exp(-data.amplitude), np.exp(data.amplitude))\n", "plot_field_panels(\n", " macro,\n", " permeability_panels,\n", " figsize=(12, 4.8),\n", " color_limits={label: permeability_bounds for label in permeability_panels},\n", ")\n", "plt.show()\n", "\n", "panels = {\n", " \"Analytical pressure\": rectangle_panel(macro, exact),\n", " \"MHM pressure\": rectangle_panel(macro, mhm_evaluator),\n", " \"Fine CG pressure\": rectangle_panel(macro, reference),\n", " \"Analytical Darcy flux magnitude\": rectangle_panel(\n", " macro, exact, quantity=\"flux_magnitude\", permeability=data.permeability\n", " ),\n", " \"MHM Darcy flux magnitude\": rectangle_panel(\n", " macro, mhm_evaluator, quantity=\"flux_magnitude\", permeability=data.permeability\n", " ),\n", " \"Fine CG Darcy flux magnitude\": rectangle_panel(\n", " macro, reference, quantity=\"flux_magnitude\", permeability=data.permeability\n", " ),\n", "}\n", "plot_field_panels(macro, panels, figsize=(15, 9))\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "61158abc", "metadata": {}, "source": [ "## 7. Separate macro size, local size and macroface resolution\n", "\n", "Configure the macro partition, local meshes and skeleton space separately to study their effects on the error. Keep $k=2$ and $\\ell=1$ and run:\n", "\n", "1. $H=1/2,1/4,1/8$ with fixed $h=1/64$ and two segments per face;\n", "2. $h=1/8,1/16,1/32,1/64$ with fixed $H=1/4$ and two segments per face;\n", "3. one, two and four segments per face with fixed $H=1/4$ and $h=1/64$.\n", "\n", "Each configuration uses the same qualified assembly quadrature (26 Gauss\n", "points per axis), including the coarsest local mesh. The error integration\n", "remains fixed independently of these refinements.\n", "\n", "With $s$ segments per macroface, each segment has length $H/s$. In the first experiment, changing $H$ therefore also changes the physical trace-segment length, even though the segment count stays fixed. In the second experiment, the macro mesh and trace space stay fixed; in the third, both volume meshes stay fixed.\n", "\n", "The driver below repeats the same visible local/global declarations. Successive observed rates use $\\log(e_i/e_{i+1})/\\log(\\eta_i/\\eta_{i+1})$, where $\\eta$ is $H$, $h$, or segment length. A plateau can reflect an error contribution from a discretization choice held fixed. Section 8 studies the asymptotic rates for its stated smooth problem and approximation spaces." ] }, { "cell_type": "code", "execution_count": null, "id": "782346c5", "metadata": {}, "outputs": [], "source": [ "def run_discretization(\n", " n: int, refinement: int, segments: int, *, quadrature_order: int = assembly_order\n", ") -> dict:\n", " \"\"\"Repeat the declared local/global equations for one explicit resolution.\"\"\"\n", " grid = CartesianMacroMesh(n, n)\n", " trace = SkeletonSpace(\n", " grid, tuple(FaceSpace.uniform(1, segments, continuous=True) for _ in grid.faces)\n", " )\n", " local = DarcyLocalProvider(grid, trace, data, 2, refinement, quadrature_order)\n", " boundary, fixed = boundary_data(trace, data.pressure, order=7)\n", " hierarchy = MeshHierarchy(\n", " grid, tuple(grid.submesh(cell, refinement) for cell in range(len(grid.cells)))\n", " )\n", " declared = bind_problem(\n", " hierarchy,\n", " bind_interface(trace, convention=\"normal\"),\n", " local,\n", " global_equation=lambda global_problem: Equation(0, global_problem.trace_load(-boundary)),\n", " retained=1,\n", " fixed=fixed,\n", " )\n", " assembled = assemble(declared, execution=ExecutionConfig(\"serial\", native_threads=1))\n", " resolved = assembled.solve()\n", " evaluate = partial(evaluate_bound_pressure, grid, resolved.field(\"pressure\"))\n", " return {\n", " \"H\": 1 / n,\n", " \"h\": 1 / (n * refinement),\n", " \"segments\": segments,\n", " \"assembly_gauss_points_per_axis\": quadrature_order,\n", " \"trace_segment_length\": 1 / (n * segments),\n", " \"global_unknowns\": assembled.matrix.shape[0],\n", " \"residual\": resolved.raw_residual,\n", " **physical_errors(evaluate, exact, data.permeability, error_points, error_weights),\n", " }\n", "\n", "\n", "macro_rows = [run_discretization(n, 64 // n, 2) for n in (2, 4, 8)]\n", "local_rows = [run_discretization(4, r, 2) for r in (2, 4, 8, 16)]\n", "trace_rows = [run_discretization(4, 16, s) for s in (1, 2, 4)]\n", "for name, rows, variable in (\n", " (\"H\", macro_rows, \"H\"),\n", " (\"h\", local_rows, \"h\"),\n", " (\"macroface segment\", trace_rows, \"trace_segment_length\"),\n", "):\n", " scales = [row[variable] for row in rows]\n", " print(name, rows)\n", " print(\"Observed pressure rates:\", observed_rates(scales, [r[\"pressure_L2\"] for r in rows]))\n", " print(\"Observed flux rates:\", observed_rates(scales, [r[\"flux_L2\"] for r in rows]))\n", " figure = plot_convergence(\n", " scales,\n", " {\n", " \"pressure L2\": [r[\"pressure_L2\"] for r in rows],\n", " \"Darcy flux L2\": [r[\"flux_L2\"] for r in rows],\n", " },\n", " )\n", " for axis in figure.axes:\n", " axis.set_xlabel(name)\n", " plt.show()\n" ] }, { "cell_type": "markdown", "id": "0f647c4a", "metadata": {}, "source": [ "## 8. Qualify the smooth asymptotic rate separately\n", "\n", "The oscillatory coefficient above resolves distinct macro, local and trace scales. Its short refinement sequences measure that chosen material; a smooth energy theorem is not transferred to them without its regularity and scale assumptions.\n", "\n", "For an independent smooth check, take the unit square, $K=I$ and\n", "\n", "$$\n", "p=\\cos(2\\pi x)\\cos(2\\pi y),\\qquad\n", "f=8\\pi^2p,\\qquad q=-\\nabla p.\n", "$$\n", "\n", "The exterior normal flux is zero and the physical pressure has zero integral. The local spaces are P1 pressure with P0 normal-flux traces, preserving their compatible constant kernel and integral gauge. The macro sequence is $n=4,8,16,32,64$. The one-level energy order one and pressure order two follow Theorem 4.1 of [Araya et al. (2013)](https://doi.org/10.1137/120888223), with the stated solution and dual regularity. This two-level P1 study also contains local finite element error; the measured orders qualify this sequence, rather than proving the theorem for arbitrary local spaces.\n", "\n", "The record identifies the executed physical problem, spaces and source provenance. This analytical qualification uses its declared mesh family; it is not labeled a matching reconstruction of the paper's original meshes. The cell below recomputes the orders from the attributed full field norms." ] }, { "cell_type": "markdown", "id": "3ee4afdb", "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": "76714b90", "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/darcy_2013_comparison.json'\n", "refinement_data = json.loads(refinement_record.read_text(encoding='utf-8'))\n", "refinement_rows = refinement_data['results']\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='primal-mhm', spaces='P1 local pressure; P0 normal-flux trace', rate_provenance='Harder, Paredes and Valentin (2013): declared finite P1/P0 study; energy order one and pressure order two observed. Araya et al. (2013), Theorem 4.1 analyzes exact local responses; finite locals add local approximation error.', 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": "fbfb9e0c", "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-darcy_multiscale_convergence", "metadata": {}, "source": [ "## References\n", "\n", "- Christopher Harder, Diego Paredes, and Frédéric Valentin (2013). *A family of Multiscale Hybrid-Mixed finite element methods for the Darcy equation with rough coefficients*, Journal of Computational Physics 245, 107–130. [DOI: 10.1016/j.jcp.2013.03.019](https://doi.org/10.1016/j.jcp.2013.03.019).\n", "\n", "- Rodolfo Araya, Christopher Harder, Diego Paredes, and Frédéric Valentin (2013). *Multiscale Hybrid-Mixed Method*, SIAM Journal on Numerical Analysis 51(6), 3505–3531. [DOI: 10.1137/120888223](https://doi.org/10.1137/120888223).\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" }, "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": "8ad37e168a5a221747703cb5cd194ee144a3fdf65b0a8983c157f1be936b55e9" } }, "nbformat": 4, "nbformat_minor": 5 }