{ "cells": [ { "cell_type": "markdown", "id": "014234ef", "metadata": {}, "source": [ "# Stokes–Brinkman boundary layers: MHM and MHM-USFEM\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 the local mixed equations and the global skeletal equation explicitly. The classical baseline is an independently assembled, globally conforming **Taylor–Hood P2/P1** solution on three fine meshes. We compare MHM with Taylor–Hood locals against MHM-USFEM with stabilized P2/P2 locals. Both MHM configurations use the **same local velocity mesh, velocity degree and macroface space**. We also assess the published single-element USFEM degree family independently of that Taylor–Hood comparison.\n", "\n", "The analytical problem is the boundary-layer example in [Araya, Harder, Poza and Valentin (2017), §3.1.2](https://doi.org/10.1016/j.cma.2017.05.027); its [2016 author preprint](https://www.ci2ma.udec.cl/pdf/pre-publicaciones2/2016/pp16-15.pdf) gives the equations and exact fields. This notebook uses a declared SW–NE triangulation. We distinguish the paper's single-element local degree family from the locally refined Taylor–Hood/USFEM comparison and from a resolved subface control. Geometry, coefficient, source and boundary data are the same throughout.\n", "\n", "Open the downloaded notebook with `jupyter lab` after installing `pymhm[notebooks,visualization]` and the compatible native DOLFINx/UFL backend. Every physical field, source, variational form, boundary condition, pressure gauge and measurement is defined in the cells below. PyMHM assembles the executed UFL forms, integrates oriented traces, condenses local equations and solves the global system.\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": "9b89c794", "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/4b921d8289bdf4527c20b34c7a2bf234b6cc7abf8db1dc6f9dd8a6ec428cd7e4/stokes_brinkman_boundary_layer-companion.zip\"\n", "COMPANION_SHA256 = \"4b921d8289bdf4527c20b34c7a2bf234b6cc7abf8db1dc6f9dd8a6ec428cd7e4\"\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/stokes_brinkman_boundary_layer.ipynb\", directory=ROOT)\n", "root = ROOT\n", "print(\"Workspace:\", ROOT)\n", "\n", "import os\n", "\n", "import json\n", "import numpy as np\n", "from pymhm import TriangleMesh, FaceSpace, SkeletonSpace\n", "from pymhm.core.equations import Equation, LocalEquations\n", "from pymhm.core.multiscale import assemble\n", "from pymhm.fem.scalar.operators import triangle_quadrature, boundary_data\n", "from pymhm.fem.scalar.triangle import nodal_space, tabulate\n", "\n", "REPORTS = ROOT / \"build/introduction\"\n", "REPORTS.mkdir(parents=True, exist_ok=True)\n", "np.set_printoptions(precision=6, suppress=True)\n", "\n", "from dataclasses import dataclass\n", "from typing import Any\n", "from threadpoolctl import threadpool_limits\n", "\n", "threadpool_limits(limits=1)\n", "\n", "from pymhm import MeshHierarchy, LocalContext, bind_interface, bind_problem\n", "from pymhm.backends.spaces import bind_space\n", "\n", "from functools import partial\n", "from examples.introduction.vector import (\n", " flow_error_norms,\n", " record_brinkman_case,\n", " brinkman_reference,\n", " execution_provenance,\n", " plot_brinkman_convergence,\n", " plot_brinkman_family_rates,\n", " plot_brinkman_family_fields,\n", " plot_brinkman_enriched_fields,\n", " plot_brinkman_profiles,\n", ")\n", "\n", "measure_case = partial(record_brinkman_case, reports=REPORTS, root=ROOT)\n" ] }, { "cell_type": "markdown", "id": "brinkman-owned-jit-explanation", "metadata": {}, "source": [ "### Runtime setup: an owned UFL compiler cache\n", "\n", "The local volume, pressure-mean and face forms use a notebook-owned temporary compiler cache. Its directory controls kernel storage only: operators, quadrature, approximation spaces, pressure gauge and solver checks are unchanged. We allow 7200 seconds for compilation; this timeout does not change any numerical tolerance.\n", "\n", "Keep this cache alive while running the study. The final cell releases it after assembly, reference solves, error checks and plots are complete. Before another solve, rerun the cache-setup cell. After an interruption, call `BRINKMAN_JIT_CACHE.cleanup()` or restart the kernel before rerunning the setup. The local/global assembly in this tutorial remains serial.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "brinkman-owned-jit-cache", "metadata": {}, "outputs": [], "source": [ "from tempfile import TemporaryDirectory\n", "\n", "BRINKMAN_JIT_CACHE = TemporaryDirectory(prefix=\"pymhm-brinkman-jit-\")\n", "BRINKMAN_JIT_OPTIONS = {\"cache_dir\": Path(BRINKMAN_JIT_CACHE.name), \"timeout\": 7200}\n" ] }, { "cell_type": "markdown", "id": "4ba2e144", "metadata": {}, "source": [ "## 1. State the operator and derive the source independently\n", "\n", "On the unit square, prescribe velocity on the whole boundary and fix the volume mean of pressure:\n", "\n", "$$\n", "\\begin{aligned}\n", "-\\nu\\Delta\\boldsymbol u+\\gamma\\boldsymbol u+\\nabla p&=\\boldsymbol f,\\\\\n", "\\nabla\\cdot\\boldsymbol u&=0,\\\\\n", "\\boldsymbol u\\big\\vert_{\\partial\\Omega}&=\\boldsymbol u_*\\big\\vert_{\\partial\\Omega},\n", "&\\int_\\Omega p&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "The published example sets $\\nu=10^{-2}$ and $\\gamma=1$. Define\n", "\n", "$$\n", "\\begin{aligned}\n", "g(s)&=\\frac{e^{(s-1)/\\nu}-e^{-1/\\nu}}{1-e^{-1/\\nu}},\\\\\n", "\\boldsymbol u_*(x,y)&=(y-g(y),\\ x-g(x)),\n", "&p_*(x,y)&=x-y.\n", "\\end{aligned}\n", "$$\n", "\n", "The two velocity components depend on the opposite coordinate, so divergence is exactly zero. Since $g''(s)=e^{(s-1)/\\nu}/[\\nu^2(1-e^{-1/\\nu})]$, the independent source is\n", "\n", "$$\n", "\\boldsymbol f(x,y)=\\nu\\big(g''(y),g''(x)\\big)\n", "+\\gamma\\boldsymbol u_*(x,y)+(1,-1).\n", "$$\n", "\n", "The PDE coefficient is **$\\nu$**, whereas the exponential layer width is **$\\nu$** too; this is a manufactured example, not a homogeneous-forcing layer of width $\\sqrt{\\nu/\\gamma}$. In particular, replacing $-\\nu\\Delta$ by $-\\nu^2\\Delta$ would define another PDE. The stable exponential expression avoids overflow." ] }, { "cell_type": "code", "execution_count": null, "id": "420d943d", "metadata": {}, "outputs": [], "source": [ "@dataclass(frozen=True)\n", "class BrinkmanLayer:\n", " \"\"\"The unit-square exact fields of Araya et al. (2017), section 3.1.2.\n", "\n", " The operator is ``-nu*Delta u+gamma*u+grad(p)`` and not ``-nu**2*Delta u``.\n", " The exact exponential layer has width ``nu``. Velocity is prescribed on\n", " the whole exterior; pressure ``x-y`` has zero volume mean. This example\n", " uses vector-Laplacian pseudostress, not symmetric Cauchy stress.\n", " \"\"\"\n", "\n", " viscosity: float = 1e-2\n", " drag: float = 1.0\n", "\n", " def __post_init__(self) -> None:\n", " \"\"\"Require positive viscosity and positive constant resistance for this tutorial.\"\"\"\n", " if not np.isfinite(self.viscosity) or self.viscosity <= 0:\n", " raise ValueError(\"viscosity must be finite and positive\")\n", " if not np.isfinite(self.drag) or self.drag <= 0:\n", " raise ValueError(\n", " \"this tutorial requires finite positive drag; zero drag needs a local kernel\"\n", " )\n", "\n", " def layer(self, coordinate: np.ndarray) -> np.ndarray:\n", " \"\"\"Evaluate ``(exp((s-1)/nu)-exp(-1/nu))/(1-exp(-1/nu))`` stably.\"\"\"\n", " return (np.exp((coordinate - 1) / self.viscosity) - np.exp(-1 / self.viscosity)) / (\n", " -np.expm1(-1 / self.viscosity)\n", " )\n", "\n", " def layer_derivative(self, coordinate: np.ndarray, order: int) -> np.ndarray:\n", " \"\"\"Differentiate the exponential once or twice in physical coordinates.\"\"\"\n", " if order not in (1, 2):\n", " raise ValueError(\"layer derivative order must be one or two\")\n", " return np.exp((coordinate - 1) / self.viscosity) / (\n", " self.viscosity**order * -np.expm1(-1 / self.viscosity)\n", " )\n", "\n", " def velocity(self, points: np.ndarray) -> np.ndarray:\n", " \"\"\"Return exact velocity ``(y-g(y), x-g(x))``.\"\"\"\n", " reversed_points = points[:, ::-1]\n", " return reversed_points - self.layer(reversed_points)\n", "\n", " def pressure(self, points: np.ndarray) -> np.ndarray:\n", " \"\"\"Return the mean-zero exact pressure ``x-y``.\"\"\"\n", " return points[:, 0] - points[:, 1]\n", "\n", " def gradient(self, points: np.ndarray) -> np.ndarray:\n", " \"\"\"Return velocity Jacobians with axes ``(point, component, derivative)``.\"\"\"\n", " result = np.zeros((len(points), 2, 2))\n", " result[:, 0, 1] = 1 - self.layer_derivative(points[:, 1], 1)\n", " result[:, 1, 0] = 1 - self.layer_derivative(points[:, 0], 1)\n", " return result\n", "\n", " def source(self, points: np.ndarray) -> np.ndarray:\n", " \"\"\"Apply the declared operator analytically, including ``grad(p)=(1,-1)``.\"\"\"\n", " return (\n", " self.viscosity * self.layer_derivative(points[:, ::-1], 2)\n", " + self.drag * self.velocity(points)\n", " + np.array([1.0, -1.0])\n", " )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f15a60ab", "metadata": {}, "outputs": [], "source": [ "NU, GAMMA = 1e-2, 1.0\n", "truth = BrinkmanLayer(NU, GAMMA)\n", "# The main family fixes P2 velocity, P1/P2 pressure and P0 multipliers.\n", "LOCAL_SUBDIVISIONS = 4\n", "LOCAL_QUADRATURE_DEGREE, ERROR_ORDER = 28, 20\n", "INVERSE_M = 1 / 100\n", "points = np.array([[0.21, 0.34], [0.72, 0.81], [0.97, 0.98]])\n", "np.testing.assert_allclose(np.trace(truth.gradient(points), axis1=1, axis2=2), 0.0)\n", "np.testing.assert_allclose(\n", " truth.source(points)\n", " - GAMMA * truth.velocity(points)\n", " - NU * truth.layer_derivative(points[:, ::-1], 2),\n", " np.tile([1.0, -1.0], (len(points), 1)),\n", " atol=1e-12,\n", ")\n", "# Assumption (M), 2D: this Taylor--Hood submesh has interior vertices.\n", "probe_local_mesh = TriangleMesh.unit_square(1).submesh(0, LOCAL_SUBDIVISIONS)\n", "boundary_vertices = np.unique(probe_local_mesh.faces[probe_local_mesh.boundary_faces])\n", "local_interior_vertices = len(probe_local_mesh.points) - len(boundary_vertices)\n", "assert local_interior_vertices > 0\n", "print(\"Taylor-Hood local interior vertices:\", local_interior_vertices)\n" ] }, { "cell_type": "markdown", "id": "bb0638c1", "metadata": {}, "source": [ "The NumPy expressions above provide physical data and independent error evaluation. The same explicit analytical data are written below as UFL expressions so that both local equations and the independent classical assembly integrate the actual exponential source." ] }, { "cell_type": "code", "execution_count": null, "id": "ac9d2329", "metadata": {}, "outputs": [], "source": [ "import basix.ufl\n", "import dolfinx\n", "import ufl\n", "\n", "\n", "def brinkman_ufl_data(domain: Any) -> tuple[Any, Any, Any]:\n", " \"\"\"Write the exact velocity, pressure and independently derived source as UFL data.\"\"\"\n", " x = ufl.SpatialCoordinate(domain)\n", " denominator = -np.expm1(-1 / NU)\n", "\n", " def g(s: Any) -> Any:\n", " \"\"\"Evaluate the analytical exponential layer in UFL coordinates.\"\"\"\n", " return (ufl.exp((s - 1) / NU) - np.exp(-1 / NU)) / denominator\n", "\n", " uexact = ufl.as_vector((x[1] - g(x[1]), x[0] - g(x[0])))\n", " pexact = x[0] - x[1]\n", " f = (\n", " ufl.as_vector(\n", " (\n", " ufl.exp((x[1] - 1) / NU) / (NU * denominator),\n", " ufl.exp((x[0] - 1) / NU) / (NU * denominator),\n", " )\n", " )\n", " + GAMMA * uexact\n", " + ufl.as_vector((1.0, -1.0))\n", " )\n", " return uexact, pexact, f\n" ] }, { "cell_type": "markdown", "id": "08ceb47a", "metadata": {}, "source": [ "## 2. Translate the mathematical mixed form directly into UFL\n", "\n", "Use the **vector Laplacian** form and its pseudostress, not the symmetric-strain Stokes operator:\n", "\n", "$$\n", "\\begin{aligned}\n", "a_K((\\boldsymbol u,p),(\\boldsymbol v,q))\n", "&=(\\nu\\nabla\\boldsymbol u,\\nabla\\boldsymbol v)_K\n", " +(\\gamma\\boldsymbol u,\\boldsymbol v)_K\\\\\n", "&\\quad-(p,\\nabla\\cdot\\boldsymbol v)_K-(q,\\nabla\\cdot\\boldsymbol u)_K.\n", "\\end{aligned}\n", "$$\n", "\n", "The pressure test is negated to obtain a symmetric saddle matrix. For USFEM, the strong residual contains **$+\\nabla p$**:\n", "\n", "$$\n", "\\begin{aligned}\n", "R(\\boldsymbol u,p)&=-\\nu\\Delta\\boldsymbol u+\\gamma\\boldsymbol u+\\nabla p,\\\\\n", "a_K^{\\rm US}&=a_K-\\sum_{T\\subset K}\\tau_T(R(\\boldsymbol u,p),R(\\boldsymbol v,q))_T,\\\\\n", "\\ell_K^{\\rm US}&=(\\boldsymbol f,\\boldsymbol v)_K\n", " -\\sum_{T\\subset K}\\tau_T(\\boldsymbol f,R(\\boldsymbol v,q))_T,\\\\\n", "\\tau_T&=\\frac{h_T^2}{\\max\\{\\gamma h_T^2,4\\nu/m\\}+4\\nu/m}.\n", "\\end{aligned}\n", "$$\n", "\n", "These are the symmetric-test version of Eqs. (41)–(42) of the article. With constant scalar drag, the 2017 minimum-eigenvalue and 2025 maximum-eigenvalue conventions coincide. Here $h_T$ is the actual fine-triangle diameter.\n", "\n", "The comparison pair is P2/P2 with an unsplit P0 trace. It satisfies the sufficient condition $k-\\ell\\geq d$ of [Araya et al. (2025), Theorem 4.1](https://doi.org/10.1137/24M1649368). Four local subdivisions provide interior vertices for the Taylor–Hood comparison, as required by Assumption (M) in Section 4.2. The single-element degree family is stated separately below. We declare the conservative inverse parameter $m=1/100$, then verify\n", "\n", "$$\n", "m h_T^2\\lVert\\Delta v_h\\rVert_T^2\\leq\\lVert\\nabla v_h\\rVert_T^2\n", "$$\n", "\n", "on the actual P2 right triangles. This is an operator bound, not a value inferred from convergence. The Taylor–Hood P2/P1 submesh has three interior vertices on each macrotriangle and satisfies the stated local-mesh hypothesis; neither configuration introduces artificial pressure regularization." ] }, { "cell_type": "code", "execution_count": null, "id": "1563013d", "metadata": {}, "outputs": [], "source": [ "from pymhm.fem.inequalities import laplacian_inverse_bound\n", "\n", "# The bound is an operation on polynomial derivatives, independent of the PDE.\n", "# Basix supplies the actual degree-dependent values; PyMHM owns the eigenproblem.\n", "inverse_parameters, inverse_controls = {}, []\n", "probe = TriangleMesh.unit_square(1)\n", "bary, weights = triangle_quadrature(8)\n", "for degree in (2, 3, 4):\n", " _, _, _, gradient, hessian = tabulate(probe, degree, bary)\n", " h = probe.lengths[probe.cell_faces].max(axis=1)\n", " selected = laplacian_inverse_bound(gradient, hessian, weights, h)\n", " np.testing.assert_allclose(selected, selected[0], rtol=1e-12, atol=1e-15)\n", " inverse_parameters[degree] = float(selected[0])\n", " stiffness = np.einsum(\"q,tqia,tqja->tij\", weights, gradient, gradient)\n", " laplacian = np.trace(hessian, axis1=-2, axis2=-1)\n", " residual_gram = np.einsum(\"q,tqi,tqj->tij\", weights, laplacian, laplacian)\n", " margins = np.linalg.eigvalsh(\n", " selected[:, None, None] * h[:, None, None] ** 2 * residual_gram - stiffness\n", " )[:, -1]\n", " assert margins.max() < 1e-10 * max(1.0, np.abs(stiffness).max())\n", " inverse_controls.append(\n", " dict(degree=degree, m=float(selected[0]), maximum_signed_margin=float(margins.max()))\n", " )\n", "# The comparison's declared conservative P2 parameter also obeys the bound.\n", "_, _, _, gradient, hessian = tabulate(probe, 2, bary)\n", "stiffness = np.einsum(\"q,tqia,tqja->tij\", weights, gradient, gradient)\n", "laplacian = np.trace(hessian, axis1=-2, axis2=-1)\n", "residual_gram = np.einsum(\"q,tqi,tqj->tij\", weights, laplacian, laplacian)\n", "bound_margins = np.linalg.eigvalsh(INVERSE_M * h[:, None, None] ** 2 * residual_gram - stiffness)[\n", " :, -1\n", "]\n", "assert bound_margins.max() < 1e-10\n", "inverse_controls\n" ] }, { "cell_type": "code", "execution_count": null, "id": "82b21132", "metadata": {}, "outputs": [], "source": [ "def brinkman_forms(\n", " domain: Any,\n", " *,\n", " stabilized: bool,\n", " velocity_degree: int = 2,\n", " inverse_m: float = INVERSE_M,\n", " quadrature_degree: int = LOCAL_QUADRATURE_DEGREE,\n", ") -> tuple[Any, Any, Any, Any]:\n", " \"\"\"Declare the executed vector-Laplacian mixed form and its USFEM residual.\"\"\"\n", " pressure_degree = velocity_degree if stabilized else velocity_degree - 1\n", " element = basix.ufl.mixed_element(\n", " [\n", " basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " velocity_degree,\n", " shape=(2,),\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " ),\n", " basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " pressure_degree,\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " ),\n", " ]\n", " )\n", " W = dolfinx.fem.functionspace(domain, element)\n", " u, p = ufl.TrialFunctions(W)\n", " v, q = ufl.TestFunctions(W)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": quadrature_degree})\n", " _, _, f = brinkman_ufl_data(domain)\n", " a = (\n", " NU * ufl.inner(ufl.grad(u), ufl.grad(v))\n", " + GAMMA * ufl.inner(u, v)\n", " - p * ufl.div(v)\n", " - q * ufl.div(u)\n", " ) * dx\n", " L = ufl.inner(f, v) * dx\n", " if stabilized:\n", " R = -NU * ufl.div(ufl.grad(u)) + GAMMA * u + ufl.grad(p)\n", " Rtest = -NU * ufl.div(ufl.grad(v)) + GAMMA * v + ufl.grad(q)\n", " h = ufl.CellDiameter(domain)\n", " tau = h**2 / (ufl.max_value(GAMMA * h**2, 4 * NU / inverse_m) + 4 * NU / inverse_m)\n", " a -= tau * ufl.inner(R, Rtest) * dx\n", " L -= tau * ufl.inner(f, Rtest) * dx\n", " return W, a, L, q * dx\n" ] }, { "cell_type": "markdown", "id": "3b132919", "metadata": {}, "source": [ "The native mesh and coefficient conversions come from the shared backend binding. The provider below writes the mixed velocity–pressure volume forms and both oriented interface pairings, and registers velocity and pressure fields. Exterior velocity data enter the global trace equation; the physical pressure gauge is imposed after assembly." ] }, { "cell_type": "code", "execution_count": null, "id": "221f81c7", "metadata": { "jupyter": { "source_hidden": true } }, "outputs": [], "source": [ "import basix\n", "\n", "\n", "from pymhm.backends.spaces import create_native_mesh, coefficient_map\n", "from pymhm.backends.forms import assemble_form, assemble_pairing\n", "\n", "\n", "def mixed_nodal_blocks(\n", " W: Any,\n", " a: Any,\n", " L: Any,\n", " mean_form: Any,\n", " fine: TriangleMesh,\n", " pressure_degree: int,\n", " velocity_degree: int = 2,\n", ") -> tuple[Any, np.ndarray, np.ndarray, int]:\n", " \"\"\"Compile supplied vector/scalar forms in declared canonical nodal coordinates.\"\"\"\n", " binding = bind_space(fine, W)\n", " order = np.r_[\n", " coefficient_map(binding, component=0),\n", " coefficient_map(binding, component=1),\n", " ]\n", " native_A = assemble_form(a, jit_options=BRINKMAN_JIT_OPTIONS)\n", " native_F = assemble_form(L, jit_options=BRINKMAN_JIT_OPTIONS)\n", " A, F = native_A[order][:, order], native_F[order]\n", " pressure_weights = assemble_form(mean_form, jit_options=BRINKMAN_JIT_OPTIONS)[order]\n", " return A, F, pressure_weights, len(nodal_space(fine, velocity_degree)[1])\n" ] }, { "cell_type": "markdown", "id": "b3763b92", "metadata": {}, "source": [ "## 3. State trace orientation, boundary data and the physical gauge\n", "\n", "With the global normal $n_F$, the multiplier is the **negative physical pseudotraction in that normal direction**:\n", "\n", "$$\n", "\\lambda_F=(-\\nu\\nabla\\boldsymbol u+pI)n_F.\n", "$$\n", "\n", "`trace_coupling` multiplies that globally oriented basis by each incident macrocell's outward-normal sign and integrates it against velocity. Thus local equations are $A_Kw_K+B_K\\lambda=F_K$, where $w_K$ stores coefficients of the mixed physical fields $(\\boldsymbol u_K,p_K)$ in the executed local basis. Choosing $C_K=-B_K^\\mathsf T$ makes the global equation\n", "\n", "$$\n", "-\\sum_K B_K^\\mathsf T w_K=-\\langle\\boldsymbol u_*,\\mu\\rangle_{\\partial\\Omega}.\n", "$$\n", "\n", "The right-hand side is a velocity moment, not prescribed traction. All exterior velocity data are imposed weakly by this equation. `Equation(0,-boundary)` makes its sign explicit. Positive drag removes the velocity-translation kernel, so there are **no retained kernel amplitudes** here.\n", "\n", "Whole-boundary velocity data leave a global pressure constant undetermined. `mean_constraint` lifts the explicit local pressure weights into a global row. We constrain the **physical reconstructed pressure integral**; pinning a skeletal coefficient would generally choose a different gauge." ] }, { "cell_type": "markdown", "id": "960d530b", "metadata": {}, "source": [ "The provider below binds the mixed velocity–pressure space to the local mesh, writes both UFL interface pairings independently, and registers named velocity and pressure fields. A shared native coefficient adapter supplies the executed mixed representation used by the following norms and archives. The pressure mean form remains explicit for the global physical gauge." ] }, { "cell_type": "code", "execution_count": null, "id": "aed5cce1", "metadata": {}, "outputs": [], "source": [ "def local_brinkman(\n", " local: LocalContext,\n", " *,\n", " macro: TriangleMesh,\n", " skeleton: SkeletonSpace,\n", " stabilized: bool,\n", " subdivisions: int = LOCAL_SUBDIVISIONS,\n", " velocity_degree: int = 2,\n", " inverse_m: float = INVERSE_M,\n", " quadrature_degree: int = LOCAL_QUADRATURE_DEGREE,\n", ") -> LocalEquations:\n", " \"\"\"Connect the declared mixed UFL equations to their oriented skeletal coordinates.\"\"\"\n", " fine = local.mesh\n", " pressure_degree = velocity_degree if stabilized else velocity_degree - 1\n", " W, a, L, mean_form = brinkman_forms(\n", " create_native_mesh(fine),\n", " stabilized=stabilized,\n", " velocity_degree=velocity_degree,\n", " inverse_m=inverse_m,\n", " quadrature_degree=quadrature_degree,\n", " )\n", " A, F, pweights, nv = mixed_nodal_blocks(\n", " W, a, L, mean_form, fine, pressure_degree, velocity_degree\n", " )\n", " binding = local.native_space(W)\n", " u, p = ufl.TrialFunctions(W)\n", " v, q = ufl.TestFunctions(W)\n", " pair_b = local.trace_pairings(lambda phi, ds: ufl.inner(phi, v) * ds)\n", " pair_c = local.trace_pairings(lambda phi, ds: -ufl.inner(phi, u) * ds, axis=\"rows\")\n", " order = np.r_[coefficient_map(binding, component=0), coefficient_map(binding, component=1)]\n", " B = assemble_pairing(\n", " pair_b.forms, axis=\"columns\", jit_options=BRINKMAN_JIT_OPTIONS\n", " )[order]\n", " C = assemble_pairing(\n", " pair_c.forms, axis=\"rows\", jit_options=BRINKMAN_JIT_OPTIONS\n", " )[:, order]\n", " # Preserve the declared (canonical velocity, canonical pressure) record layout.\n", " identity = np.eye(len(F))\n", " native_layout = binding.to_native(identity[: 2 * nv], component=0) + binding.to_native(\n", " identity[2 * nv :], component=1\n", " )\n", " local.field(\"velocity\", binding, component=0, reconstruction=native_layout)\n", " local.field(\"pressure\", binding, component=1, reconstruction=native_layout)\n", " return local.equations(\n", " a=A,\n", " L=F,\n", " b=B,\n", " c=C,\n", " metadata={\n", " \"mesh\": fine,\n", " \"velocity_nodes\": nv,\n", " \"velocity_degree\": velocity_degree,\n", " \"inverse_m\": inverse_m if stabilized else None,\n", " \"pressure_degree\": pressure_degree,\n", " \"pressure_weights\": pweights,\n", " },\n", " )\n" ] }, { "cell_type": "markdown", "id": "be15288e", "metadata": {}, "source": [ "Before the refinement loop, define the measurements: velocity and pressure $L^2$ errors, velocity-gradient error, broken pseudostress error, the fine-cell divergence norm and the separate macro mass defect. These quadrature operations evaluate fields; they do not define a local PDE." ] }, { "cell_type": "markdown", "id": "37bf2d0f", "metadata": {}, "source": [ "The measurements use the general Basix kernel with the declared nodal order and affine geometry. The adapter below requests values and the first derivatives required by the norms; it does not select or assemble an operator.\n", "Physical norms, executed-state archives and one-sided field/profile displays are importable from `examples/introduction/vector.py`. The local forms and the complete global solve remain in this notebook, together with separately declared conforming reference forms.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2c389259", "metadata": { "jupyter": { "source_hidden": true } }, "outputs": [], "source": [ "# Physical field norms, macro mass and pressure integral use flow_error_norms (imported above).\n" ] }, { "cell_type": "markdown", "id": "03d6d16f", "metadata": {}, "source": [ "### Record the numerical coordinates for reproducibility\n", "\n", "The next utility saves the actual local coordinate basis, coefficients, orientation maps and reconstructed fields as each solve finishes. It contains no physical operator. Replaying the saved basis is checked with one and two BLAS threads; this positive-reaction or positive-drag example has no local nullspace modes. You can read the global construction immediately below without studying the archive format first." ] }, { "cell_type": "code", "execution_count": null, "id": "3504a5bd", "metadata": {}, "outputs": [], "source": [ "def solve_brinkman_case(\n", " n,\n", " stabilized,\n", " *,\n", " subdivisions=LOCAL_SUBDIVISIONS,\n", " degree=2,\n", " trace_degree=0,\n", " segments=1,\n", " inverse_m=INVERSE_M,\n", "):\n", " \"\"\"Declare the mesh, trace space, weak velocity data and physical pressure gauge.\"\"\"\n", " macro = TriangleMesh.unit_square(n)\n", " skeleton = SkeletonSpace(\n", " macro, tuple(FaceSpace.uniform(trace_degree, segments) for _ in macro.faces), components=2\n", " )\n", " boundary, fixed = boundary_data(skeleton, truth.velocity, {}, order=32)\n", " assert not fixed # All exterior velocity data enter the global trace load.\n", " provider = lambda local: local_brinkman(\n", " local,\n", " macro=macro,\n", " skeleton=skeleton,\n", " stabilized=stabilized,\n", " subdivisions=subdivisions,\n", " velocity_degree=degree,\n", " inverse_m=inverse_m,\n", " )\n", " hierarchy = MeshHierarchy(\n", " macro, tuple(macro.submesh(c, subdivisions) for c in range(len(macro.cells)))\n", " )\n", " problem = bind_problem(\n", " hierarchy,\n", " bind_interface(skeleton, convention=\"normal\"),\n", " provider,\n", " global_equation=Equation(0, -boundary),\n", " retained=0,\n", " )\n", " system = assemble(problem)\n", " gauge = system.mean_constraint(\n", " [data[\"pressure_weights\"] for data in system.local_metadata], 0.0\n", " )\n", " solution = system.solve(constraints=[gauge])\n", " return macro, skeleton, system, solution, gauge\n", "\n", "\n", "METHODS = ((\"MHM Taylor-Hood\", False), (\"MHM-USFEM\", True))\n", "cases, rows, state_archives = {}, [], []\n", "for n in (4, 8, 16):\n", " for method, stabilized in METHODS:\n", " macro, skeleton, system, solution, gauge = solve_brinkman_case(n, stabilized)\n", " case, row, archive = measure_case(\n", " macro,\n", " skeleton,\n", " system,\n", " solution,\n", " truth,\n", " velocity_degree=2,\n", " pressure_degree=2 if stabilized else 1,\n", " order=ERROR_ORDER,\n", " name=f\"brinkman-n{n}-{method.lower().replace(' ', '-')}\",\n", " named=True,\n", " gauge=gauge if n == 4 else None,\n", " method=method,\n", " n=n,\n", " )\n", " cases[n, method] = case\n", " rows.append(row)\n", " state_archives.append(archive)\n", " print(\n", " method,\n", " n,\n", " {key: row[key] for key in (\"velocity_l2\", \"pressure_l2\", \"macro_mass_defect\")},\n", " )\n" ] }, { "cell_type": "markdown", "id": "d6017d6d", "metadata": {}, "source": [ "### Assess the published single-element degree family\n", "\n", "[Araya et al. (2017), Section 3.1 and Figures 6–8](https://doi.org/10.1016/j.cma.2017.05.027) uses one local triangle per macrotriangle and the equal-order family\n", "\n", "$$\n", "\\begin{aligned}\n", "\\ell&=0,1,2, & k&=\\ell+2,\\\\\n", "\\Lambda_H\\vert_F&=[P_\\ell(F)]^2,\n", "& (\\boldsymbol u_h,p_h)\\vert_K&\\in[P_k(K)]^2\\times P_k(K).\n", "\\end{aligned}\n", "$$\n", "\n", "We select these spaces explicitly, keeping the same analytical problem and global pressure gauge. The `velocity_degree` argument changes the UFL element, its executed field representation and interface integration together. There are no fine local submeshes in this experiment. The mesh parameter is the actual macrotriangle diameter $H=\\sqrt2/n$.\n", "\n", "The paper reports eventual rates $\\ell+2$ for velocity $L^2$ and $\\ell+1$ for pressure $L^2$ and broken velocity-gradient error, with loss of rates while the layer is unresolved. These are literature observations for this example; the measured slopes below determine which regime our resolutions reach. The inverse coefficient is computed from each executed polynomial space, without using the exact solution or fitted error. The historical numerical inverse constants and connectivity are not supplied by the article, so this is a comparison of its PDE and degree family rather than a literal reproduction of its plotted values.\n", "\n", "The default introductory profile uses $n=8,16,32$ for each of the three degrees.\n", "It measures the same PDE and published approximation family on three resolutions;\n", "these levels can remain preasymptotic, so they do not establish the eventual\n", "literature rates. The main Taylor–Hood/USFEM study at $n=4,8,16$, all three refined\n", "classical references, and the enriched face/local-space control are unchanged.\n", "\n", "Set `PYMHM_FULL_STUDY=1` before execution to acquire the original full\n", "qualification sequence: $n=8,16,32,64,128$ for $\\ell=0$ and\n", "$n=8,16,32,64$ for $\\ell=1,2$. The selected profile and actual levels are recorded\n", "in the numerical provenance; either profile computes its own fields and slopes.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cc12d53b", "metadata": {}, "outputs": [], "source": [ "# Set PYMHM_FULL_STUDY=1 to extend each family to its complete qualification grids.\n", "FULL_STUDY = os.environ.get(\"PYMHM_FULL_STUDY\", \"0\") == \"1\"\n", "PUBLISHED_FULL_LEVELS = {0: (8, 16, 32, 64, 128), 1: (8, 16, 32, 64), 2: (8, 16, 32, 64)}\n", "PUBLISHED_LEVELS = PUBLISHED_FULL_LEVELS if FULL_STUDY else {ell: (8, 16, 32) for ell in (0, 1, 2)}\n", "STUDY_PROFILE = \"full qualification\" if FULL_STUDY else \"introductory three-mesh refinement profile\"\n", "print(\"Single-element family profile:\", STUDY_PROFILE, PUBLISHED_LEVELS)\n", "published_cases, published_rows = {}, []\n", "for ell, levels in PUBLISHED_LEVELS.items():\n", " degree = ell + 2\n", " for n in levels:\n", " macro, skeleton, system, solution, _ = solve_brinkman_case(\n", " n,\n", " True,\n", " subdivisions=1,\n", " degree=degree,\n", " trace_degree=ell,\n", " inverse_m=inverse_parameters[degree],\n", " )\n", " case, row, archive = measure_case(\n", " macro,\n", " skeleton,\n", " system,\n", " solution,\n", " truth,\n", " velocity_degree=degree,\n", " pressure_degree=degree,\n", " order=ERROR_ORDER,\n", " name=f\"brinkman-single-ell{ell}-n{n}\",\n", " method=\"MHM-USFEM single element\",\n", " ell=ell,\n", " local_degree=degree,\n", " n=n,\n", " local_subdivisions=1,\n", " inverse_m=inverse_parameters[degree],\n", " )\n", " published_cases[ell, n] = case\n", " published_rows.append(row)\n", " state_archives.append(archive)\n", " print(\n", " \"single-element USFEM\",\n", " ell,\n", " degree,\n", " n,\n", " {key: row[key] for key in (\"velocity_l2\", \"pressure_l2\", \"velocity_h1_seminorm\")},\n", " )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c23166f3", "metadata": {}, "outputs": [], "source": [ "assembly_quadrature_controls = []\n", "probe_macro = TriangleMesh.unit_square(PUBLISHED_LEVELS[0][0])\n", "for ell in (0, 1, 2):\n", " degree = ell + 2\n", " probe_skeleton = SkeletonSpace(\n", " probe_macro, tuple(FaceSpace.uniform(ell) for _ in probe_macro.faces), components=2\n", " )\n", " candidates = np.flatnonzero(\n", " np.any(np.isclose(probe_macro.points[probe_macro.cells, 0], 1.0), axis=1)\n", " )\n", " cell_index = int(candidates[0]) # A macrotriangle touching the actual layer x=1.\n", " hierarchy = MeshHierarchy(\n", " probe_macro, tuple(probe_macro.submesh(c, 1) for c in range(len(probe_macro.cells)))\n", " )\n", " probe_problem = bind_problem(\n", " hierarchy, bind_interface(probe_skeleton, convention=\"normal\"), lambda local: None\n", " )\n", " operators = [\n", " local_brinkman(\n", " probe_problem.local_context(cell_index),\n", " macro=probe_macro,\n", " skeleton=probe_skeleton,\n", " stabilized=True,\n", " subdivisions=1,\n", " velocity_degree=degree,\n", " inverse_m=inverse_parameters[degree],\n", " quadrature_degree=q,\n", " )\n", " for q in (LOCAL_QUADRATURE_DEGREE, max(40, LOCAL_QUADRATURE_DEGREE + 12))\n", " ]\n", " standard, higher = operators\n", " np.testing.assert_allclose(standard.a.toarray(), higher.a.toarray(), rtol=1e-11, atol=1e-13)\n", " np.testing.assert_allclose(standard.L, higher.L, rtol=1e-10, atol=1e-13)\n", " assembly_quadrature_controls.append(\n", " dict(\n", " ell=ell,\n", " local_degree=degree,\n", " cell=cell_index,\n", " quadrature_degrees=[LOCAL_QUADRATURE_DEGREE, max(40, LOCAL_QUADRATURE_DEGREE + 12)],\n", " matrix_maximum_difference=float(np.max(np.abs((standard.a - higher.a).toarray()))),\n", " load_maximum_difference=float(np.max(np.abs(standard.L - higher.L))),\n", " )\n", " )\n", "assembly_quadrature_controls\n" ] }, { "cell_type": "markdown", "id": "ae0edc1c", "metadata": {}, "source": [ "### Resolve the layer while retaining a simple macro mesh\n", "\n", "The unsplit P0 sequence and the degree family assess their measured layer regime separately. A small algebraic residual and macro mass defect do not imply a resolved velocity gradient. To separate this limitation from the local mixed solver, we add a **distinct resolution control**: eight local subdivisions per macro edge and eight independent P0 subfaces per original macroface, on macro grids $n=4,8,16$. The actual mesh scales are $\\mathcal H=\\sqrt{2}/n$ for macrotriangles and $\\bar H=h=\\mathcal H/8$ for their aligned subface and local partitions.\n", "\n", "Both mixed methods retain P2 velocity and the same refined local mesh; pressure is P1 for Taylor–Hood and P2 for USFEM. Trace integration resolves the declared subfaces and includes the same orientation convention. This enriched skeletal discretization is reported separately from the unsplit polynomial experiment. We assess its physical errors, pressure gauge and macro balance through the same assembly; we do not transfer the paper's unsplit-face convergence theorem to this control.\n", "\n", "At $n=16$, the global geometry has only 512 macrotriangles, versus 32768 triangles in the classical reference. Fine resolution is confined to independent local problems and subface coefficients. Increasing only local resolution while leaving the trace fixed would retain the original face approximation error.\n", "\n", "The initial enriched mesh remains relatively coarse compared with the layer width. Three mesh resolutions provide two observed slopes without assuming that the first pair is asymptotic." ] }, { "cell_type": "code", "execution_count": null, "id": "f7256793", "metadata": {}, "outputs": [], "source": [ "control_rows, enriched_cases = [], {}\n", "TRACE_SUBFACES, CONTROL_LOCAL_SUBDIVISIONS = 8, 8\n", "for n in (4, 8, 16):\n", " for method, stabilized in METHODS:\n", " macro, skeleton, system, solution, _ = solve_brinkman_case(\n", " n, stabilized, subdivisions=CONTROL_LOCAL_SUBDIVISIONS, segments=TRACE_SUBFACES\n", " )\n", " case, row, archive = measure_case(\n", " macro,\n", " skeleton,\n", " system,\n", " solution,\n", " truth,\n", " velocity_degree=2,\n", " pressure_degree=2 if stabilized else 1,\n", " order=ERROR_ORDER,\n", " name=f\"brinkman-enriched-n{n}-{method.lower().replace(' ', '-')}\",\n", " method=method,\n", " n=n,\n", " trace_subfaces=TRACE_SUBFACES,\n", " local_subdivisions=CONTROL_LOCAL_SUBDIVISIONS,\n", " )\n", " enriched_cases[n, method] = case\n", " control_rows.append(row)\n", " state_archives.append(archive)\n", " print(\n", " \"enriched\",\n", " method,\n", " n,\n", " {key: row[key] for key in (\"velocity_l2\", \"pressure_l2\", \"velocity_h1_seminorm\")},\n", " )\n" ] }, { "cell_type": "markdown", "id": "4dfba502", "metadata": {}, "source": [ "## 4. Independently assemble a classical Taylor–Hood baseline\n", "\n", "The fine reference is one conforming P2/P1 mesh, with no skeleton, local condensation or residual stabilization. DOLFINx/UFL assembles the same operator with the independently derived source. Velocity is prescribed **strongly** on all boundary nodes; pressure has the same zero-volume-mean gauge. Strong and weak boundary enforcement give different discrete spaces, so identical coefficients are not expected.\n", "\n", "The native pressure test $q$ supplies the volume-moment row. The reference solve eliminates the fixed velocity coordinates and adds one exact mean constraint. It uses the package's checked linear solver without changing its tolerances. We measure true errors in three fine meshes because an exact solution is available; a fine numerical solution remains a numerical baseline." ] }, { "cell_type": "code", "execution_count": null, "id": "43fa3d8c", "metadata": {}, "outputs": [], "source": [ "def classical_brinkman_forms(space):\n", " \"\"\"Declare the separate conforming P2/P1 energy, source and pressure moment.\"\"\"\n", " u, p = ufl.TrialFunctions(space)\n", " v, q = ufl.TestFunctions(space)\n", " domain = space.ufl_domain()\n", " uexact, pexact, f = brinkman_ufl_data(domain)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 20})\n", " a = (\n", " NU * ufl.inner(ufl.grad(u), ufl.grad(v))\n", " + GAMMA * ufl.inner(u, v)\n", " - p * ufl.div(v)\n", " - q * ufl.div(u)\n", " ) * dx\n", " return a, ufl.inner(f, v) * dx, q * dx, uexact, pexact, dx\n", "\n", "\n", "references, reference_rows = {}, []\n", "for n in (32, 64, 128):\n", " case, row, archive = brinkman_reference(\n", " n, classical_brinkman_forms, truth, reports=REPORTS, root=ROOT\n", " )\n", " references[n] = case\n", " reference_rows.append(row)\n", " state_archives.append(archive)\n" ] }, { "cell_type": "markdown", "id": "1811ca6f", "metadata": {}, "source": [ "## 5. Measure convergence, without assuming the asymptotic regime\n", "\n", "$H=\\sqrt2/n$ denotes the actual maximum macrotriangle diameter. The Taylor–Hood/USFEM comparison has four local subdivisions per macro edge; the published degree-family experiment has one local triangle per macrotriangle. Their trace spaces and resulting measured errors are reported independently. The study keeps local degree, trace degree, subdivision factor, operator and boundary data fixed while halving $H$.\n", "\n", "A layer may be unresolved at the first levels. Plotting the observed rates makes this visible. The published smooth-space orders cannot be asserted before the exponential layer and the local approximation error are adequately resolved. The smallest singular value checks uniqueness of a particular gauged matrix; it does not establish a uniform inf-sup bound for an arbitrary mesh family.\n", "\n", "The macro mass defect is $\\max_K\\left|\\int_K\\nabla\\cdot\\boldsymbol u_h\\right|$. It is distinct from the reported fine-cell divergence $L^2$ norm: macro conservation does not make the raw velocity pointwise divergence-free." ] }, { "cell_type": "markdown", "id": "3aef63f0", "metadata": {}, "source": [ "For two successive resolutions, the measured slope is $r=\\log(e_1/e_2)/\\log(H_1/H_2)$. These small plotting utilities compute that quantity directly and keep error and rate axes separate." ] }, { "cell_type": "code", "execution_count": null, "id": "d65b0485", "metadata": {}, "outputs": [], "source": [ "plot_brinkman_convergence(rows, reference_rows, control_rows, reports=REPORTS)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c43ab7c0", "metadata": {}, "outputs": [], "source": [ "published_rates = plot_brinkman_family_rates(published_rows, reports=REPORTS)\n" ] }, { "cell_type": "markdown", "id": "7140781d", "metadata": {}, "source": [ "## 6. Inspect velocity, pressure, error and one-sided profiles\n", "\n", "Every spatial panel highlights the actual **macro mesh**, including analytical and classical-reference panels. The fields are sampled from their original polynomials. Pressure remains independently reconstructed in each macrocell; coincident endpoints are not averaged. Velocity magnitude labels refer to velocity, not Darcy flux.\n", "\n", "We show the enriched $n=16$ MHM fields and a detailed profile in the layer near $x=1$. The source has nonzero pressure gradient, so pressure accuracy is a separate measurement rather than a zero-pressure special case." ] }, { "cell_type": "markdown", "id": "c14dbe0d", "metadata": {}, "source": [ "The display utility evaluates each local polynomial independently, including its one-sided boundary values. Every panel overlays the actual macro mesh and has its own color scale. It changes only the display sampling." ] }, { "cell_type": "markdown", "id": "9b9dde21", "metadata": {}, "source": [ "The additional single-element figure evaluates the $\\ell=2$, P4/P4 degree-family solution on its own macro mesh. Its coefficients and pressure mean come from that experiment, with no subface enrichment. The refined Taylor–Hood solution remains the independently assembled numerical baseline.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cb592c79", "metadata": {}, "outputs": [], "source": [ "plot_brinkman_family_fields(published_cases, PUBLISHED_LEVELS, references, truth, reports=REPORTS)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "38bc9ad2", "metadata": {}, "outputs": [], "source": [ "plot_brinkman_enriched_fields(enriched_cases, references, truth, reports=REPORTS)\n" ] }, { "cell_type": "markdown", "id": "34a568a9", "metadata": {}, "source": [ "The classical pressure reference uses the same physical zero-mean gauge. The additional comparison below shows that field and the separate pressure errors of both multiscale methods." ] }, { "cell_type": "markdown", "id": "c4707152", "metadata": {}, "source": [ "Horizontal profiles retain a separate segment for each incident macrocell. Vertical markers identify macroface crossings; neighboring endpoint values are evaluated independently instead of being averaged." ] }, { "cell_type": "code", "execution_count": null, "id": "bb8cdcf2", "metadata": {}, "outputs": [], "source": [ "plot_brinkman_profiles(enriched_cases, truth, reports=REPORTS)\n" ] }, { "cell_type": "markdown", "id": "8dda511e", "metadata": {}, "source": [ "## 7. Check quadrature and record the precise method provenance\n", "\n", "Native local assembly uses UFL quadrature degree 28, and physical errors use 20 Gauss points in each Duffy coordinate. Reintegrating the finest MHM solution with 28 points checks that the layer norms are not quadrature artifacts. The baseline errors use a separately specified UFL degree-20 rule on the fine conforming meshes.\n", "\n", "The raw pseudostress $\\nu\\nabla\\boldsymbol u_h-p_hI$ is a broken derived field. Its $L^2$ error is reported, but we do not label it H(div) or infer fine-cell mass conservation. The singular-value and residual checks supplement the physical convergence and macro-balance measurements." ] }, { "cell_type": "code", "execution_count": null, "id": "e41d1c63", "metadata": { "jupyter": { "source_hidden": true } }, "outputs": [], "source": [ "provenance = execution_provenance(\n", " \"notebooks/introduction/stokes_brinkman_boundary_layer.ipynb\",\n", " root=ROOT,\n", " native_jit_cache=\"notebook-owned temporary directory; removed after final result save\",\n", " study_profile=STUDY_PROFILE,\n", " single_element_levels=PUBLISHED_LEVELS,\n", " full_qualification_levels=PUBLISHED_FULL_LEVELS,\n", " literature_comparison=\"same analytical PDE; declared tutorial discretization, no matched-figure reproduction\",\n", ")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "76667dc4", "metadata": {}, "outputs": [], "source": [ "quadrature_checks = {}\n", "for method in (\"MHM Taylor-Hood\", \"MHM-USFEM\"):\n", " macro, meshes, u, p, pd = cases[16, method]\n", " higher = flow_error_norms(meshes, u, p, 2, pd, truth, order=28)\n", " measured = next(row for row in rows if row[\"n\"] == 16 and row[\"method\"] == method)\n", " for key in (\"velocity_l2\", \"pressure_l2\", \"velocity_h1_seminorm\", \"pseudostress_l2\"):\n", " np.testing.assert_allclose(higher[key], measured[key], rtol=1e-8, atol=1e-12)\n", " quadrature_checks[method] = higher\n", "for method in (\"MHM Taylor-Hood\", \"MHM-USFEM\"):\n", " macro, meshes, u, p, pd = enriched_cases[16, method]\n", " higher = flow_error_norms(meshes, u, p, 2, pd, truth, order=28)\n", " measured = next(r for r in control_rows if r[\"n\"] == 16 and r[\"method\"] == method)\n", " for key in (\"velocity_l2\", \"pressure_l2\", \"velocity_h1_seminorm\", \"pseudostress_l2\"):\n", " np.testing.assert_allclose(higher[key], measured[key], rtol=1e-8, atol=1e-12)\n", " quadrature_checks[method + \" enriched\"] = higher\n", "for ell in (0, 1, 2):\n", " selected = [row for row in published_rows if row[\"ell\"] == ell]\n", " finest = selected[-1]\n", " macro, meshes, u, p, degree = published_cases[ell, finest[\"n\"]]\n", " higher = flow_error_norms(meshes, u, p, degree, degree, truth, order=28)\n", " for key in (\"velocity_l2\", \"pressure_l2\", \"velocity_h1_seminorm\", \"pseudostress_l2\"):\n", " np.testing.assert_allclose(higher[key], finest[key], rtol=1e-8, atol=1e-12)\n", " quadrature_checks[f\"USFEM single-element ell={ell}\"] = higher\n", "(REPORTS / \"stokes-brinkman-boundary-layer.json\").write_text(\n", " json.dumps(\n", " {\n", " \"problem\": \"Araya et al. 2017 section 3.1.2, independently derived source\",\n", " \"source_url\": \"https://doi.org/10.1016/j.cma.2017.05.027\",\n", " \"viscosity\": NU,\n", " \"drag\": GAMMA,\n", " \"velocity_degree\": 2,\n", " \"local_subdivisions\": LOCAL_SUBDIVISIONS,\n", " \"pressure_degrees\": {\"MHM Taylor-Hood\": 1, \"MHM-USFEM\": 2},\n", " \"trace_degree\": 0,\n", " \"inverse_m\": INVERSE_M,\n", " \"inverse_bound_margin\": float(bound_margins.max()),\n", " \"boundary\": \"full velocity Dirichlet, weak MHM / strong conforming reference\",\n", " \"pressure_gauge\": \"zero physical volume mean\",\n", " \"reference_backend\": f\"DOLFINx {dolfinx.__version__}; independent P2/P1 UFL\",\n", " \"local_ufl_quadrature_degree\": LOCAL_QUADRATURE_DEGREE,\n", " \"error_duffy_order\": ERROR_ORDER,\n", " \"convergence\": rows,\n", " \"reference_refinement\": reference_rows,\n", " \"quadrature_checks\": quadrature_checks,\n", " \"enriched_control\": control_rows,\n", " \"single_element_family\": published_rows,\n", " \"single_element_rates\": published_rates,\n", " \"inverse_controls\": inverse_controls,\n", " \"assembly_quadrature_controls\": assembly_quadrature_controls,\n", " \"taylor_hood_local_interior_vertices\": local_interior_vertices,\n", " \"mesh_size_convention\": \"actual maximum macrotriangle diameter H=sqrt(2)/n\",\n", " \"theory_url\": \"https://doi.org/10.1137/24M1649368\",\n", " \"state_archives\": state_archives,\n", " \"provenance\": provenance,\n", " },\n", " indent=2,\n", " )\n", ")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "brinkman-owned-jit-cleanup", "metadata": {}, "outputs": [], "source": [ "# Local assembly and all postprocessing are complete.\n", "BRINKMAN_JIT_CACHE.cleanup()\n" ] }, { "cell_type": "markdown", "id": "08430e85", "metadata": {}, "source": [ "The resulting rates describe these declared two-dimensional spaces and resolutions. The exact PDE and local residual are taken from [Araya et al. (2017)](https://doi.org/10.1016/j.cma.2017.05.027); the tutorial uses the explicitly declared SW–NE mesh and distinguishes the single-element P2/P2, P3/P3 and P4/P4 degree family from the refined P2 comparison. Reported pseudostress is a broken raw derived field. The independently assembled classical comparison uses [DOLFINx](https://docs.fenicsproject.org/dolfinx/) through UFL and the project's checked sparse linear solver.\n", "\n", "This configuration requires positive drag. In the pure Stokes limit, the vector-Laplacian local operator has two retained velocity translations in two dimensions; their basis moments and coarse amplitudes must be declared in the general problem API. That limit needs its own kernel configuration." ] }, { "cell_type": "markdown", "id": "a783dc03", "metadata": {}, "source": [ "## 8. Compare the two smooth approximation families\n", "\n", "The boundary-layer study above measures a difficult, distinct regime. To assess asymptotic approximation rates, use the smooth polynomial Stokes data with $\\nu=1$, $\\gamma=0$ and zero-volume-mean pressure:\n", "\n", "$$\n", "\\begin{aligned}\n", " \\psi(x,y)&=128x^2(1-x)^2y^2(1-y)^2,\\\\\n", " u&=(\\partial_y\\psi,-\\partial_x\\psi),\\\\\n", " p&=150(x-\\tfrac12)(y-\\tfrac12),\\qquad\n", " f=-\\Delta u+\\nabla p.\n", "\\end{aligned}\n", "$$\n", "\n", "These are the analytical data evaluated in the attributed PyMHM records. Section 3.1.1 of Araya et al. (2017) uses a negative streamfunction prefactor, so its velocity and force differ. This comparison assesses smooth convergence with the explicitly displayed data.\n", "\n", "The same displayed volume form loses its reaction term. In this pure Stokes case, its two local null modes are $Z_j=(e_j,0)$, $j=1,2$. Particular local responses have zero velocity moments $\\int_T u_j$; their two translation amplitudes remain global unknowns (`retained=2`). The total physical velocity is not assigned zero mean. The pressure gauge remains the single physical condition $\\int_\\Omega p=0$.\n", "\n", "The two rate studies deliberately retain their different approximation spaces:\n", "\n", "- **MHM-USFEM:** local stabilized P3/P3 on one triangle, P1 vector traces and a crisscross macro grid, $H=1/n$.\n", "- **MHM with Taylor–Hood locals:** stable P2/P1 on 16 fine triangles per macrotriangle (four subdivisions per edge), P1 vector traces and the diagonal macro grid, $H=\\sqrt2/n$.\n", "\n", "The Taylor–Hood local mesh satisfies the interior-vertex Assumption (M). Its P2/P1 locals and P1 trace have $k-\\ell=1$, outside the sufficient $k-\\ell\\ge2$ condition of [Araya et al. (2025), Theorem 4.3](https://doi.org/10.1137/24M1649368). The concrete saddle operator and trace lifting are checked independently; these finite checks do not extend that theorem. The USFEM P3/P3–P1 family meets Theorem 4.2, which estimates the broken velocity–pressure product norm rather than an additional velocity L2 order. Order-three velocity L2 and order-two pressure L2 are observed reference orders for these smooth studies. [Araya et al. (2017), Section 2.2](https://doi.org/10.1016/j.cma.2017.05.027) states the general local stability and trace-lifting requirements.\n", "\n", "Read the attributed acquisition records below. The graphs retain every measured level and show all consecutive orders; their highlighted window always uses the final four levels. Reading a numerical record does not recompute its PDE or change its recorded source identity. The independent local/global equation checks, pressure gauge and quadrature controls are recorded separately from the physical field errors.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e54a141a", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from IPython.display import SVG, display\n", "\n", "from examples.tutorial_convergence import refinement_series, plot_method_series, asymptotic_summary\n", "\n", "usfem_path = ROOT / \"examples/results/stokes-adaptive/polynomial-uniform-nu1-g0-l1.json\"\n", "usfem_record = json.loads(usfem_path.read_text())\n", "usfem_rows = [dict(row, H=2 / np.sqrt(row[\"macro_cells\"])) for row in usfem_record[\"rows\"]]\n", "usfem_series = refinement_series(\n", " usfem_path, usfem_rows, \"H\",\n", " {\"velocity L2\": \"velocity_l2\", \"pressure L2\": \"pressure_l2\"},\n", " {\"velocity L2\": 3.0, \"pressure L2\": 2.0}, method=\"stokes-brinkman\",\n", " spaces=\"Stabilized P3/P3 on one triangle; P1 vector traces; crisscross macro mesh\",\n", " rate_provenance=\"Araya et al. (2017): observed velocity L2 order; Araya et al. (2025), Theorem 4.2: product-norm bound\",\n", " refinement=\"macro diameter\", root=ROOT,\n", ")\n", "\n", "hood_path = ROOT / \"examples/results/tutorial-methods/stokes-taylor-hood-asymptotic.json\"\n", "hood_record = json.loads(hood_path.read_text())\n", "hood_series = refinement_series(\n", " hood_path, hood_record[\"rows\"], \"H\",\n", " {\"velocity L2\": \"velocity_l2\", \"pressure L2\": \"pressure_l2\"},\n", " {\"velocity L2\": 3.0, \"pressure L2\": 2.0}, method=\"stokes-galerkin\",\n", " spaces=\"Taylor-Hood P2/P1, r4 local mesh; unsplit vector P1 traces; diagonal macro mesh\",\n", " rate_provenance=\"Observed smooth refined Taylor-Hood family; k-ell=1 is outside the sufficient gap of Theorem 4.3 (2025)\",\n", " refinement=\"macro diameter\", root=ROOT,\n", ")\n", "for series in (usfem_series, hood_series):\n", " print(series[\"method\"], \"input\", series[\"record\"], \"sha256\", series[\"sha256\"])\n", " print(asymptotic_summary(series))\n", " figure = plot_method_series(series)\n", " convergence_path = REPORTS / (series[\"method\"] + \"-smooth-convergence.svg\")\n", " figure.savefig(convergence_path, bbox_inches=\"tight\",\n", " metadata={\"Creator\": \"IPES Research Group\", \"Rights\": \"CC BY 4.0\"})\n", " display(SVG(filename=str(convergence_path)))\n", " plt.close(figure)\n" ] }, { "cell_type": "markdown", "id": "references-stokes_brinkman_boundary_layer", "metadata": {}, "source": [ "## References\n", "\n", "- Rodolfo Araya, Christopher Harder, Abner H. Poza, and Frédéric Valentin (2017). *Multiscale hybrid-mixed method for the Stokes and Brinkman equations—The method*, Computer Methods in Applied Mechanics and Engineering 324, 29–53. [DOI: 10.1016/j.cma.2017.05.027](https://doi.org/10.1016/j.cma.2017.05.027). Author preprint (2016): [CI²MA Preprint 2016-15](https://www.ci2ma.udec.cl/pdf/pre-publicaciones2/2016/pp16-15.pdf).\n", "\n", "- Rodolfo Araya, Christopher Harder, Abner H. Poza, and Frédéric Valentin (2025). *Multiscale Hybrid-Mixed Methods for the Stokes and Brinkman Equations—A Priori Analysis*, SIAM Journal on Numerical Analysis 63(2), 588–618. [DOI: 10.1137/24M1649368](https://doi.org/10.1137/24M1649368)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (Pixi introduction)", "language": "python", "name": "python3" }, "language_info": { "name": "python" }, "pymhm_editorial_revision": { "date": "2026-10-06", "preserved": "Code cells, execution counts, execution timestamps and numerical outputs", "preserved_code_cells_sha256": "3ed2c097fdfe4175e70462f98499ef3a6884cca3a7d22383bf61eee263ade661", "scope": "Markdown citations and bibliography only" } }, "nbformat": 4, "nbformat_minor": 5 }