{ "cells": [ { "cell_type": "markdown", "id": "4a42675c", "metadata": {}, "source": [ "# Darcy on an original SPE10 layer: a small macro mesh and complete material data\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", "Build the local/global API explicitly for the complete layer 36 of SPE10 Model 2: all $60\\times220$ material pixels are retained, without cropping, smoothing or contrast clipping. MHM uses a $6\\times11$ macro grid. The classical reference uses a much finer single conforming grid.\n", "\n", "The geometry, layer and Darcy boundary data follow [Paredes, Valentin and Versieux (2024)](https://doi.org/10.1016/j.cam.2023.115415). Their local refinement is not specified; ours is explicit, so this tutorial does not claim a matched discrete reproduction.\n", "\n", "Install `pymhm[notebooks,visualization]` and the compatible native DOLFINx/UFL backend, then open this notebook with `jupyter lab` in a writable directory. Physical coefficients, local/global equations and conforming reference forms are declared below. Importable scalar helpers handle physical evaluation, norms, plots and archives. Fine reference solves may need several minutes and substantial memory.\n", "\n", "The primal hybrid construction follows [Harder, Paredes and Valentin (2013)](https://doi.org/10.1016/j.jcp.2013.03.019); the material layer and face-based comparison above follow [Paredes, Valentin and Versieux (2024)](https://doi.org/10.1016/j.cam.2023.115415).\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": "8f13eeb3", "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/1d2ec8b0d25291cfc5e4f3f4165491bc39ff665c3c0d7594dd8b2216d1eeed13/darcy_spe10_layer-companion.zip\"\n", "COMPANION_SHA256 = \"1d2ec8b0d25291cfc5e4f3f4165491bc39ff665c3c0d7594dd8b2216d1eeed13\"\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_spe10_layer.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", " dirichlet_solve,\n", " evaluate_bound_pressure,\n", " evaluate_qk,\n", " fine_quad_panel,\n", " broken_quad_panel,\n", " fine_line_profile,\n", " grid_quadrature,\n", " load_spe10_layer,\n", " native_scalar_space,\n", " physical_errors,\n", " plot_convergence,\n", " plot_field_panels,\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", "from pymhm.materials.cartesian import CartesianCellField\n", "from matplotlib.collections import LineCollection\n" ] }, { "cell_type": "markdown", "id": "aa2843cc", "metadata": {}, "source": [ "## 1. Read the unchanged data and verify provenance\n", "\n", "The compact archive preserves the original OPM horizontal entries $K_{xx},K_{yy}$. The acquisition helper verifies SHA256 and returns `CartesianCellField`: it selects a declared side at an interface, without averaging. Provenance identifies the project, revision, source URL and axis convention. Coordinates are feet, permeability is mD, and the prescribed pressure difference is dimensionless." ] }, { "cell_type": "code", "execution_count": null, "id": "20fd5705", "metadata": {}, "outputs": [], "source": [ "permeability, provenance = load_spe10_layer(ROOT, number=36)\n", "print(\n", " {\n", " \"source\": provenance[\"repository\"],\n", " \"revision\": provenance[\"revision\"],\n", " \"layer\": provenance[\"selected_layer\"],\n", " \"spacing_ft\": permeability.spacing,\n", " \"shape\": permeability.values.shape,\n", " \"spectral_contrast\": float(\n", " permeability.values[..., 0, 0].max() / permeability.values[..., 0, 0].min()\n", " ),\n", " }\n", ")\n", "macro = CartesianMacroMesh(6, 11, (0.0, 1200.0, 0.0, 2200.0))\n", "pressure_boundary = lambda points: 1 - points[:, 1] / 2200.0\n", "source = 0.0\n", "# K in mD, coordinates in ft; pressure difference is dimensionless in this benchmark.\n", "# Darcy flux therefore has the corresponding benchmark units mD/ft.\n", "local_degree, local_refinement, trace_segments = 1, 80, 16\n", "skeleton = SkeletonSpace(\n", " macro, tuple(FaceSpace.uniform(1, trace_segments, continuous=True) for _ in macro.faces)\n", ")\n" ] }, { "cell_type": "markdown", "id": "411cbdd8", "metadata": {}, "source": [ "## 2. Same material, operator and boundary data\n", "\n", "$$\n", "\\begin{aligned}\n", "\\boldsymbol q&=-K\\nabla p,\\\\\n", "\\nabla\\cdot\\boldsymbol q&=0\n", "\\quad\\text{in }(0,1200)\\times(0,2200),\\\\\n", "p(x,0)&=1,\\qquad p(x,2200)=0,\\\\\n", "\\boldsymbol q\\cdot\\boldsymbol n&=0\n", "\\quad\\text{on }x=0\\text{ and }x=1200.\n", "\\end{aligned}\n", "$$\n", "\n", "`pressure_boundary` supplies values only to the horizontal Dirichlet faces. Vertical faces have homogeneous physical no-flow data. Every macroelement measures $200\\times200$ ft; material pixels measure $20\\times10$ ft. A local refinement divisible by 20 resolves their interfaces.\n", "\n", "Our explicitly integrated local grids align with every material pixel boundary. Gauss integration inside these elements therefore integrates each original constant pixel value. Changing the discretization never changes $K$. Our local refinement is 80 in each direction, while the macro grid still has only 66 elements. For this fitted grid, a tensor-valued native DG0 coefficient represents each unchanged pixel exactly. The optional alternative operator shown later also supports exact pixel intersections for nonaligned grids." ] }, { "cell_type": "markdown", "id": "0ae1ef38", "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", "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 plotted raw flux $-K\\nabla p_h$ is not an $H(\\mathrm{div})$ reconstruction and is not guaranteed to conserve on every fine cell.\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": "d16da2f8", "metadata": {}, "outputs": [], "source": [ "def ufl_coefficient_source(domain: Any, material: CartesianCellField) -> tuple[Any, Any]:\n", " \"\"\"Represent unchanged pixel tensors in native DG0 on a pixel-fitted grid.\"\"\"\n", " import basix.ufl\n", " from dolfinx import fem\n", "\n", " tensor_space = fem.functionspace(\n", " domain, basix.ufl.element(\"DG\", \"quadrilateral\", 0, shape=(2, 2))\n", " )\n", " coefficient = fem.Function(tensor_space)\n", " # DG0 interpolation points are cell interiors, so each fitted cell retains its incident pixel.\n", " coefficient.interpolate(lambda x: material(x[:2].T).reshape(-1, 4).T)\n", " force = fem.Constant(domain, np.float64(0.0))\n", " return coefficient, 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", " ratios = np.asarray(self.data.spacing) / fine.spacing\n", " assert np.allclose(ratios, np.round(ratios)) and np.all(ratios >= 1), (\n", " \"DG0 must represent pixels exactly: local cells must fit material boundaries\"\n", " )\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": "e7436676", "metadata": {}, "source": [ "## 4. Declare global terms and solve MHM\n", "\n", "The sidewalls fix the physical normal-flux multiplier to zero. The horizontal pressure boundary contributes the negative right-hand side for $C=-B^T$. One constant mode is retained per macroelement. Prescribed pressure removes the global gauge.\n", "\n", "The visible provider and global declaration define the entire formulation. `assemble` and `system.solve()` own the generic condensation, global solution and reconstruction." ] }, { "cell_type": "markdown", "id": "bcac92a9", "metadata": {}, "source": [ "### Evaluate pressure in its executed coefficient basis\n", "\n", "Each named pressure view stores its executed native basis and the corresponding portable Qk coefficient order. `evaluate_bound_pressure` selects a macrocell and evaluates that view's independent coefficients and physical gradients. The solution coefficients represent pressure; the material tensor $K$ is applied separately to obtain the physical Darcy flux $-K\\nabla p_h$." ] }, { "cell_type": "code", "execution_count": null, "id": "57cc5c4e", "metadata": {}, "outputs": [], "source": [ "no_flow = {int(face): 0.0 for face in macro.boundary_faces if abs(macro.normals[face, 0]) > 0.5}\n", "boundary, fixed = boundary_data(skeleton, pressure_boundary, no_flow, order=5)\n", "provider = DarcyLocalProvider(macro, skeleton, permeability, local_degree, local_refinement, 5)\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", "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", " \"macro_cells\": len(macro.cells),\n", " \"local_refinement\": local_refinement,\n", " \"trace_degree\": 1,\n", " \"trace_segments\": trace_segments,\n", " \"global_unknowns\": system.matrix.shape[0],\n", " \"fixed_trace_coefficients\": len(fixed),\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": "4714de6c", "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": "2dc2f174", "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, local_degree, permeability=permeability, source=source, 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": "d3489a65", "metadata": {}, "source": [ "## 5. Conforming Galerkin references and their own uncertainty\n", "\n", "There is no analytical solution here. Assemble classical conforming $Q_3$ on three pixel-aligned grids: $60\\times220$, $120\\times440$ and $240\\times880$. Impose $p=1$ and $p=0$ strongly on horizontal boundaries. The sidewalls use the natural zero-flux condition with the same physical Darcy sign.\n", "\n", "The fine field is a **numerical reference**, not an exact solution. Measure differences between successive references in pressure, vector flux and flux energy. High contrast and intersecting material interfaces do not inherit smooth-problem rates.\n", "\n", "The common integration grid resolves pixels, local element boundaries and reference boundaries. Norms evaluate each executed coefficient basis directly, rather than interpolating archived images. This is an independent classical assembly sharing basis kernels, not execution of an external paper code." ] }, { "cell_type": "code", "execution_count": null, "id": "39bfa59d", "metadata": {}, "outputs": [], "source": [ "reference_rows, reference_evaluators, reference_meshes, reference_coefficients = [], [], [], []\n", "for nx, ny in ((60, 220), (120, 440), (240, 880)):\n", " print(\n", " f\"Assembling conforming Q3 reference {nx} x {ny}; {(3 * nx + 1) * (3 * ny + 1)} coefficients\",\n", " flush=True,\n", " )\n", " fine = CartesianMacroMesh(nx, ny, macro.bounds)\n", " domain, space, mapping = native_scalar_space(fine, 3)\n", " p, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n", " K, f = ufl_coefficient_source(domain, permeability)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 7})\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, 3)\n", " exterior = np.flatnonzero(np.isclose(nodes[:, 1], 0.0) | np.isclose(nodes[:, 1], 2200.0))\n", " coefficients = dirichlet_solve(\n", " a_cg, load_cg, mapping[exterior], pressure_boundary(nodes[exterior])\n", " )[mapping]\n", " reference_evaluators.append(partial(evaluate_qk, fine, 3, coefficients))\n", " reference_meshes.append(fine)\n", " reference_coefficients.append(coefficients)\n", " reference_rows.append({\"shape\": (nx, ny), \"unknowns\": len(nodes)})\n", "reference_rows\n" ] }, { "cell_type": "code", "execution_count": null, "id": "27fb6c82", "metadata": {}, "outputs": [], "source": [ "error_points, error_weights = grid_quadrature(macro.bounds, (480, 880), order=4)\n", "reference_refinement = [\n", " physical_errors(a, b, permeability, error_points, error_weights)\n", " for a, b in zip(reference_evaluators[:-1], reference_evaluators[1:])\n", "]\n", "reference = reference_evaluators[-1]\n", "comparison = physical_errors(mhm_evaluator, reference, permeability, error_points, error_weights)\n", "print(\"Reference refinement:\", reference_refinement)\n", "print(\"MHM versus finest conforming reference:\", comparison)\n", "# Refinement evidence is a field measurement, not only a linear residual.\n", "assert reference_refinement[-1][\"pressure_L2\"] < reference_refinement[0][\"pressure_L2\"]\n", "assert reference_refinement[-1][\"flux_L2\"] < reference_refinement[0][\"flux_L2\"]\n" ] }, { "cell_type": "markdown", "id": "1f48614c", "metadata": {}, "source": [ "### Plot the baseline refinement evidence\n", "\n", "The plot uses successive **reference differences**, not errors against an exact solution. The observed slope measures how those increments decrease under simultaneous reference refinement; it does not establish the asymptotic error rate of this heterogeneous problem. Pressure and physical flux can settle at different speeds." ] }, { "cell_type": "code", "execution_count": null, "id": "cade0c04", "metadata": {}, "outputs": [], "source": [ "figure = plot_convergence(\n", " [20.0, 10.0],\n", " {\n", " \"reference pressure L2 increment\": [r[\"pressure_L2\"] for r in reference_refinement],\n", " \"reference Darcy flux L2 increment\": [r[\"flux_L2\"] for r in reference_refinement],\n", " },\n", ")\n", "for axis in figure.axes:\n", " axis.set_xlabel(\"Reference x-cell width [ft]\")\n", "plt.show()\n", "print(\n", " \"Reference uncertainty/MHM difference ratios:\",\n", " {\n", " field: reference_refinement[-1][field] / comparison[field]\n", " for field in (\"pressure_L2\", \"flux_L2\", \"flux_energy\")\n", " },\n", ")\n" ] }, { "cell_type": "markdown", "id": "812f5a39", "metadata": {}, "source": [ "## 6. Fields with the actual macro mesh\n", "\n", "Relate the permeability channels to the physical flux and compare pressure under identical geometry and boundary data. Each panel retains independent macrocell values. The flux fields are raw $-K\\nabla p_h$, not $H(\\mathrm{div})$ reconstructions or averaged multipliers.\n", "\n", "Permeability pixels are displayed exactly with nearest-neighbor rendering. The flux-magnitude display uses a logarithmic color scale to show channels across the large contrast; the physical field norms above retain the untransformed vector flux.\n", "\n", "Each flux panel samples the actual fine elements independently. The tensor is evaluated at the interior of its pixel-fitted element and applied to that element's own gradient at all four corners. Separate vertices preserve both material-interface and finite-element limits, instead of querying a global coefficient at a shared pixel vertex." ] }, { "cell_type": "code", "execution_count": null, "id": "499115ba", "metadata": {}, "outputs": [], "source": [ "fig, axis = plt.subplots(figsize=(6, 8), layout=\"constrained\")\n", "artist = axis.imshow(\n", " np.log10(permeability.values[..., 0, 0]).T,\n", " origin=\"lower\",\n", " extent=macro.bounds,\n", " interpolation=\"nearest\",\n", " aspect=\"equal\",\n", ")\n", "axis.add_collection(LineCollection(macro.points[macro.faces], colors=\"0.15\", linewidths=0.7))\n", "axis.set(title=\"SPE10 layer 36: log10 Kxx [mD]\", xlabel=\"x [ft]\", ylabel=\"y [ft]\")\n", "fig.colorbar(artist, ax=axis, shrink=0.9)\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1793ecf1", "metadata": {}, "outputs": [], "source": [ "mhm_pressure = broken_quad_panel(system, solution, local_degree, permeability, \"pressure\")\n", "reference_pressure = fine_quad_panel(\n", " reference_meshes[-1], 3, reference_coefficients[-1], permeability, \"pressure\"\n", ")\n", "pressure_difference = (*mhm_pressure[:2], mhm_pressure[2] - reference(mhm_pressure[0])[0])\n", "mhm_flux = broken_quad_panel(system, solution, local_degree, permeability, \"flux_magnitude\")\n", "reference_flux = fine_quad_panel(\n", " reference_meshes[-1], 3, reference_coefficients[-1], permeability, \"flux_magnitude\"\n", ")\n", "# Display log10(magnitude expressed in mD/ft); physical norms remain untransformed.\n", "mhm_flux = (*mhm_flux[:2], np.log10(np.maximum(mhm_flux[2], 1e-12)))\n", "reference_flux = (*reference_flux[:2], np.log10(np.maximum(reference_flux[2], 1e-12)))\n", "panels = {\n", " \"MHM pressure\": mhm_pressure,\n", " \"Fine Galerkin Q3 pressure\": reference_pressure,\n", " \"Pressure difference MHM − reference\": pressure_difference,\n", " \"log10 MHM Darcy flux magnitude [mD/ft]\": mhm_flux,\n", " \"log10 reference Darcy flux magnitude [mD/ft]\": reference_flux,\n", "}\n", "plot_field_panels(macro, panels, figsize=(15, 12))\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "a3232ef8", "metadata": {}, "source": [ "## 7. Macro conservation and profiles with both interface limits\n", "\n", "The local constant mode gives $\\int_{\\partial T}s_{TF}\\lambda=\\int_T f=0$. Integrate the oriented skeletal flux in every macroelement to verify this **macro** balance. This does not prove fine-cell conservation of the raw gradient flux.\n", "\n", "Evaluate the $x=199$ ft profile separately on the eleven macro intervals. Preserve both interface limits, and mark macroface crossings with vertical lines.\n", "The profile uses each explicitly selected fine element and its constant incident tensor. It retains both endpoint limits with disconnected line segments, so neither material jumps nor macroface traces are averaged.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "90011ab2", "metadata": {}, "outputs": [], "source": [ "macro_balance = []\n", "for cell in range(len(macro.cells)):\n", " total = 0.0\n", " for side, face in enumerate(macro.cell_faces[cell]):\n", " t, w = skeleton.faces[face].quadrature(5)\n", " normal_flux = skeleton.faces[face].evaluate(t) @ solution.trace[skeleton.dofs(int(face))]\n", " total += macro.signs[cell, side] * macro.lengths[face] * (w @ normal_flux)\n", " macro_balance.append(total)\n", "print(\"Maximum integrated macro imbalance:\", float(np.max(np.abs(macro_balance))))\n", "fig, axes = plt.subplots(1, 2, figsize=(12, 4), constrained_layout=True)\n", "for row in range(11):\n", " field = pressure_fields[row * 6]\n", " y, p, qy = fine_line_profile(\n", " field.mesh, local_degree, field.portable_coefficients, permeability, 199.0\n", " )\n", " axes[0].plot(y, p, color=\"C0\", label=\"MHM\" if row == 0 else None)\n", " axes[1].plot(y, qy, color=\"C0\", label=\"MHM\" if row == 0 else None)\n", "y, pref, qref = fine_line_profile(\n", " reference_meshes[-1], 3, reference_coefficients[-1], permeability, 199.0\n", ")\n", "axes[0].plot(y, pref, color=\"C1\", label=\"Fine Galerkin\")\n", "axes[1].plot(y, qref, color=\"C1\", label=\"Fine Galerkin\")\n", "for ax, title in zip(axes, (\"Pressure at x=199 ft\", \"Darcy flux y at x=199 ft\")):\n", " for y_face in range(200, 2200, 200):\n", " ax.axvline(y_face, color=\"0.7\", lw=0.7)\n", " ax.set_xlabel(\"y [ft]\")\n", " ax.set_title(title)\n", " ax.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "2f53b0b2", "metadata": {}, "source": [ "## 8. What this reference comparison supports\n", "\n", "Successive reference differences quantify baseline resolution. Compare them with MHM–reference differences: when baseline uncertainty matters for a field norm, refine the reference further. Pressure often settles before physical flux under this contrast; similar pressure plots do not certify flux accuracy.\n", "\n", "For the stated discretizations, the finest reference increment is about 5.4% in physical flux, compared with a 15.6% MHM–reference difference. The flux comparison remains limited by reference resolution: this is a difference against a finite approximation, not the true flux error. Pressure is more settled, with a 0.021% reference increment compared with a 1.05% MHM–reference difference. Further reference refinement is required for more precise flux accuracy.\n", "\n", "The method uses 66 macroelements and the declared spaces: local $Q_1$ pressure and a continuous piecewise $P_1$ trace on each macroface. Increase `local_refinement` or `trace_segments` explicitly for further studies, preserving the original material. Do not alter the approximation family solely to fit an image.\n", "\n", "Primary data source: [OPM/opm-data, SPE10 Model 2](https://github.com/OPM/opm-data/tree/eaa2261683a97027e057c2bc49612ad1c86390b3/spe10model2). Acquisition prints the exact revision and SHA256. The opening attribution distinguishes the literature problem from a matched reproduction. This notebook executes PyMHM and a separately assembled classical baseline, not an external article implementation." ] }, { "cell_type": "markdown", "id": "references-darcy_spe10_layer", "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", "- Diego Paredes, Frédéric Valentin, and Henrique M. Versieux (2024). *Revisiting the robustness of the multiscale hybrid-mixed method: The face-based strategy*, Journal of Computational and Applied Mathematics 436, 115415. [DOI: 10.1016/j.cam.2023.115415](https://doi.org/10.1016/j.cam.2023.115415)." ] } ], "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": "ce5b45a67c55e6a92a468fb36879a283c60c58ceb2bd8437607fdf1c18612f68" } }, "nbformat": 4, "nbformat_minor": 5 }