High-Order Pore-Level Multiscale Method (hPLMM)

• High-Order Pore-Level Multiscale Method (hPLMM) is a computational framework that embeds fine-scale pore physics into coarse-scale models through enriched mortar and interface representations.
• It leverages localized basis constructions and adaptive enrichment strategies to enhance simulation accuracy, reducing GMRES iterations by up to five times compared to low-order methods.
• hPLMM extends to inertial pore-flow and multi-continuum problems, offering a scalable approach for simulating complex mechanics in fractured and heterogeneous media.

Searching arXiv for the specified papers and closely related hPLMM material. High-Order Pore-Level Multiscale Method (hPLMM) denotes a multiscale computational framework in which fine-scale pore or microstructural physics are encoded into localized basis functions, reduced spaces, or coarse operators so that large, ill-conditioned systems on geometrically and materially complex media can be solved with substantially reduced cost while retaining pore-sensitive accuracy. In the most explicit formulation available in the supplied literature, hPLMM is presented as a two-level preconditioner for linear elasticity on arbitrary porous and fractured solids, where “high-order” refers to enriched interface representations through mortar spaces rather than rigid interface traces (Khan et al., 24 Sep 2025). The same literature also uses hPLMM in a broader methodological sense: as a natural high-order upgrade path for pore-scale Oseen multiscale methods in non-periodic porous media, and as a close analogue of higher-order homogenization, localized orthogonal decomposition, and generalized multiscale finite element constructions for heterogeneous multi-continuum flow (Muljadi, 2016, Dong et al., 7 Apr 2026).

1. Terminology and conceptual scope

In the elasticity formulation, hPLMM generalizes the earlier low-order Pore-Level Multiscale Method (PLMM) by replacing one displacement degree of freedom per interface component with multiple mortar nodes and associated mortar basis functions. The defining distinction is that PLMM treats interfaces as effectively rigid, whereas hPLMM allows non-uniform interface displacements and therefore captures local bending and twisting moments under challenging loading conditions (Khan et al., 24 Sep 2025). In this usage, “high-order” is an interface-space statement: the approximation space on subdomain contacts is enriched so that the coarse space represents more than rigid-body transmission.

The supplied literature also attaches “high-order” to several other, compatible mechanisms. In the Oseen-based pore-flow setting, the low-order Crouzeix–Raviart MsFEM is explicitly described as a foundation on which a high-order pore-level multiscale method can be built by increasing local finite-element order, enforcing higher moments on coarse edges, adding bubble-type enrichments, and performing adaptive enrichment in regions of strong inertia (Muljadi, 2016). In the scalar elliptic setting, high-order refers to coarse polynomial degree $p$ and to a localized multiscale space achieving energy-norm convergence of order $H^{p+2}$ under regularity assumptions on the right-hand side $f$, without restrictive assumptions on the coefficient, the domain, or the exact solution (Dong et al., 2022).

A common misconception is that hPLMM must mean higher polynomial degree alone. The supplied papers do not support that reduction. They instead indicate several distinct but related meanings: enriched interface traces, higher-order edge moments, second-order asymptotic reconstructions, and higher-order localized coarse spaces. This suggests that hPLMM is best understood as a class of pore-aware high-fidelity multiscale constructions whose defining feature is not a single discretization choice but the systematic elevation of coarse-scale representation quality.

2. Governing formulations and model classes

The elasticity version of hPLMM targets the small-strain linear-elastic deformation problem posed on the solid phase $\Omega_s\subset\mathbb{R}^D$, with

$$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$

together with Dirichlet and Neumann boundary conditions on $\Gamma^D$ and $\Gamma^N$. The baseline material law is isotropic with Lamé parameters $\lambda,\mu$, but spatially varying stiffness and fractures are allowed, including stiffness contrast up to $10^6$ in the reported experiments (Khan et al., 24 Sep 2025).

In pore-scale inertial flow, the relevant microscopic model is the steady Oseen approximation in the pore space $\Omega_{\text{pore}}$,

$H^{p+2}$0

with no slip on grain boundaries, non-periodic grain patterns, and Reynolds numbers ranging from creeping-flow regimes to local values of approximately $H^{p+2}$1 in the heterogeneous-$H^{p+2}$2 tests. The method is formulated specifically to avoid any periodicity assumption and to operate on randomly placed grains and highly non-periodic microstructures (Muljadi, 2016).

For overlapping multi-continuum media, the higher-order homogenization line considers microscopic pressure fields $H^{p+2}$3, $H^{p+2}$4, coupled by strong exchange terms of order $H^{p+2}$5, with micro-periodic coefficients $H^{p+2}$6, $H^{p+2}$7, and $H^{p+2}$8. The resulting homogenized equations contain effective storage coefficients $H^{p+2}$9, effective permeability tensors $f$0, effective transfer coefficients $f$1, and cross-gradient tensors $f$2, so the macroscopic model is not merely dual-porosity diffusion but an enhanced multi-continuum system with additional coupling terms (Dong et al., 7 Apr 2026).

Related formulations in the supplied corpus broaden the admissible physics further. The conforming primal-dual mixed method treats multiscale Darcy flow on two macro-regions with weak interface laws for pressure and flux balance, while the GMsFEM formulations address multicontinuum poroelasticity and slightly compressible single-phase flow with local spectral coarse spaces and residual-driven online enrichment (Morales, 2017, Tyrylgin et al., 2019, Fu et al., 2022). Taken together, these formulations show that hPLMM is not restricted to one PDE family; rather, it is a multiscale design pattern spanning elasticity, inertial pore flow, Darcy-type exchange, poroelastic coupling, and compressible flow.

3. Local basis construction and interface coupling

The elasticity hPLMM begins from a decomposition of the solid domain into non-overlapping grain grids $f$3, shared contact interfaces $f$4, and thin overlapping contact grids $f$5. On each grain grid, local shape functions and a correction function are computed by solving elasticity problems with zero-stress conditions on void boundaries and unknown interface displacements represented through mortars. The fine-scale displacement on each grain grid is approximated by a superposition of one correction function and multiple shape functions weighted by coarse unknowns associated with mortar nodes and local directions (Khan et al., 24 Sep 2025).

The mortar representation is the key coupling mechanism. Instead of a constant interface trace, the interface displacement functions are expanded in mortar basis functions $f$6, where $f$7 indexes mortar nodes placed quasi-uniformly on each interface and $f$8 indexes local displacement directions. The global coarse problem is obtained by enforcing mortar-weighted traction equilibrium and kinematic continuity across every contact. This produces a coarse linear system whose solution reconstructs a first-pass fine-scale approximation and, in the algebraic realization, defines the prolongation matrix used by the two-level preconditioner (Khan et al., 24 Sep 2025).

In the Oseen MsFEM lineage, local basis construction is instead edge-centered and nonconforming. For each coarse edge $f$9, one solves a local Oseen problem on the two-element patch $\Omega_s\subset\mathbb{R}^D$0, with edge-average constraints

$\Omega_s\subset\mathbb{R}^D$1

The resulting Crouzeix–Raviart-type basis functions satisfy only weak continuity,

$\Omega_s\subset\mathbb{R}^D$2

so coarse edges may cut through grains and irregular interfaces without boundary-fitted meshing or oversampling (Muljadi, 2016).

Other supplied frameworks show alternative coupling principles that are structurally relevant to hPLMM. The primal-dual mixed formulation uses asymmetric regularity—$\Omega_s\subset\mathbb{R}^D$3 velocity and $\Omega_s\subset\mathbb{R}^D$4 pressure in one subregion, $\Omega_s\subset\mathbb{R}^D$5 pressure and $\Omega_s\subset\mathbb{R}^D$6-type velocity structure in another—and enforces interface laws through bilinear forms rather than strong continuity constraints (Morales, 2017). The GMsFEM variants construct local snapshot spaces on coarse neighborhoods $\Omega_s\subset\mathbb{R}^D$7, solve spectral problems to select dominant modes, and then assemble global multiscale spaces; in compressible flow, residual-driven online basis functions are computed from local residual problems and added adaptively when the offline space is insufficient, especially for singular sources (Tyrylgin et al., 2019, Fu et al., 2022).

4. High-order mechanisms

The most explicit high-order mechanism in hPLMM is the mortar interface space. Gaussian mortars are built from normalized Gaussian weights centered at mortar nodes, with spread controlled by a parameter $\Omega_s\subset\mathbb{R}^D$8; algebraic mortars are obtained by solving tangential interface elasticity problems. In both cases, the mortar functions form partition-of-unity-type interface spaces and permit multiple displacement modes per interface. When the number of mortar nodes per interface is $\Omega_s\subset\mathbb{R}^D$9, hPLMM reduces to low-order PLMM; increasing $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$0 enriches the interface representation and allows non-rigid interfacial deformation (Khan et al., 24 Sep 2025).

In the Oseen setting, the supplied literature identifies several natural routes to a high-order pore-level method. These include higher-order local finite elements such as Q2–Q1 or Q2–Q2 in the basis solves, higher-order nonconforming constraints that enforce edge moments

$$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$1

bubble-type enrichments supported inside coarse elements, and adaptive enrichment by additional local problems in regions with strong local Reynolds number or recirculation (Muljadi, 2016). The significance is that the nonconforming coarse space can be enriched without abandoning the core advantages of weak continuity and grain-crossing coarse meshes.

Higher-order homogenization supplies another mechanism: systematic asymptotic reconstruction. In the HOMS formulation, the microscopic solution is expanded as

$$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$2

where first-order cell functions $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$3 and $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$4 capture gradient and transfer effects, while second-order cell functions $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$5, $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$6, $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$7, $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$8, and $$\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{f}, \qquad \boldsymbol{\sigma}(\boldsymbol{u})=\mathbf{C}:\nabla^s\boldsymbol{u},$$9 encode transient, curvature, directional coupling, and refined transfer corrections. The resulting second-order reconstruction $\Gamma^D$0 is the high-order multiscale field used to recover pore-scale oscillations from the homogenized solution (Dong et al., 7 Apr 2026).

A different high-order mechanism appears in the localized orthogonal decomposition framework. There, the multiscale basis functions are $\Gamma^D$1, where $\Gamma^D$2 is a localized corrector operator and $\Gamma^D$3 is a modified projector that stabilizes localization by treating the lowest-order modes differently from the higher ones. The method constructs coarse problem-adapted ansatz spaces by solving auxiliary problems on local subdomains and is specifically designed to cure the error deterioration observed when local subdomains are not large enough (Dong et al., 2022). This is directly relevant to hPLMM because pore-level local solves are expensive; any stabilization of localization radius materially affects feasibility.

5. Accuracy, convergence, and computational behavior

For the elasticity-preconditioning realization, the reported behavior is dominated by two observations. First, the first-pass solution produced by the hPLMM coarse space is already substantially better than that of low-order PLMM, especially under shear loading, where rigid-interface traces are inadequate. Second, when used inside GMRES, hPLMM combined with the contact–grain smoother reduces iteration counts by factors up to approximately $\Gamma^D$4 relative to PLMM, GDSW, and cAMG, with the largest improvement occurring when the number of mortar nodes per interface increases from $\Gamma^D$5 to $\Gamma^D$6 (Khan et al., 24 Sep 2025). The same study reports robustness on disk packs, sandstone variants, a highly heterogeneous DARCY domain, and 3D bone microstructure, including cases with stiffness contrast $\Gamma^D$7, fractures, and complex void–solid topology.

The Oseen MsFEM results establish monotone convergence in non-periodic porous configurations. In Poiseuille flow without grains and $\Gamma^D$8, the $\Gamma^D$9 error decreases from approximately $\Gamma^N$0 at $\Gamma^N$1 to approximately $\Gamma^N$2 at $\Gamma^N$3. In pattern (a) with $\Gamma^N$4 grains, $\Gamma^N$5, and $\Gamma^N$6, the $\Gamma^N$7 error decreases from approximately $\Gamma^N$8 to approximately $\Gamma^N$9. In pattern (b) with $\lambda,\mu$0 small grains, $\lambda,\mu$1, and $\lambda,\mu$2, the $\lambda,\mu$3 error decreases from approximately $\lambda,\mu$4 to approximately $\lambda,\mu$5. For a heterogeneous $\lambda,\mu$6 reaching local $\lambda,\mu$7, the $\lambda,\mu$8 error decreases from approximately $\lambda,\mu$9 to approximately $10^6$0 (Muljadi, 2016). These data are used in the supplied text to argue that weak continuity across coarse edges is sufficient even at large $10^6$1 and in worst-case grain-edge intersections.

The higher-order homogenization results provide both asymptotic and computational accuracy statements. Under the stated regularity and symmetry assumptions, the HOMS reconstruction satisfies an $10^6$2 error bound in both $10^6$3 and $10^6$4. Numerically, in a 2D porous example with $10^6$5, HOMS achieves relative $10^6$6 errors of approximately $10^6$7 at $10^6$8 and relative $10^6$9 seminorm errors of approximately $\Omega_{\text{pore}}$0–$\Omega_{\text{pore}}$1, while reducing runtime from approximately $\Omega_{\text{pore}}$2 s for the full fine-scale FEM to approximately $\Omega_{\text{pore}}$3 s for the cell-plus-homogenized HOMS workflow (Dong et al., 7 Apr 2026).

The supporting frameworks exhibit analogous convergence patterns. The conforming primal-dual mixed method proves well-posedness through a Babuška–Brezzi argument and reports $\Omega_{\text{pore}}$4 convergence for velocity in $\Omega_{\text{pore}}$5 and $\Omega_{\text{pore}}$6, $\Omega_{\text{pore}}$7 for $\Omega_{\text{pore}}$8 in $\Omega_{\text{pore}}$9, and $H^{p+2}$00 for $H^{p+2}$01 in $H^{p+2}$02 with $H^{p+2}$03 in $H^{p+2}$04 for the RT(0)–P1–P0 discretization (Morales, 2017). The high-order localized elliptic method achieves energy-norm errors of order $H^{p+2}$05 for sufficiently smooth $H^{p+2}$06 and, unlike the earlier Maier localization, shows stagnation rather than deterioration when $H^{p+2}$07 decreases at fixed localization depth $H^{p+2}$08 (Dong et al., 2022). In compressible-flow GMsFEM, residual-driven online basis functions are reported as particularly important for singular sources, with the numerical tables showing sharp reductions in both $H^{p+2}$09 and $H^{p+2}$10 errors once even a small number of online functions is added (Fu et al., 2022).

6. Applications, limitations, and directions of development

The elasticity hPLMM is presented for applications that require repeated solution of large, ill-conditioned linear systems on arbitrary porous or fractured structures. The reported application areas include risk analysis of subsurface $H^{p+2}$11 storage and optimization of porous materials for batteries, prosthetics, and aircraft (Khan et al., 24 Sep 2025). In these settings, the physical interpretation of the coarse space matters: each coarse degree of freedom corresponds to a local interface motion, and each local basis function encodes how that motion propagates through nearby grains or structural members.

In flow problems, the supplied literature positions hPLMM-type constructions as tools for both direct simulation and upscaling. The Oseen MsFEM is described as naturally suited for estimating effective permeability $H^{p+2}$12, inertial coefficients such as Forchheimer $H^{p+2}$13, and variability induced by non-periodic microstructure, without relying on periodic representative elementary volumes (Muljadi, 2016). The multi-continuum homogenization framework similarly converts unit-cell solutions into effective storage, permeability, transfer, and cross-gradient couplings (Dong et al., 7 Apr 2026). The GMsFEM and primal-dual mixed papers then show how such upscaled structures can be solved efficiently on large domains with explicit interface exchange, fracture coupling, or poroelastic feedback (Morales, 2017, Tyrylgin et al., 2019, Fu et al., 2022).

Several limitations are explicit. In the penalized Oseen formulation, penalty parameters must be tuned relative to $H^{p+2}$14 and $H^{p+2}$15 to avoid excessive stiffness or numerical locking, and robust linear solvers or preconditioners are required for large-scale extensions (Muljadi, 2016). In the elasticity hPLMM, increasing the number of mortar nodes raises coarse-space quality but also enlarges the coarse solve; the supplied study notes that gains beyond about $H^{p+2}$16 suggest the need for multi-level coarsening (Khan et al., 24 Sep 2025). In the HOMS setting, the method relies on periodic microstructure and asymptotic scale separation, whereas the Oseen and hPLMM elasticity formulations do not (Dong et al., 7 Apr 2026, Muljadi, 2016).

A further point of interpretation concerns scope. The supplied literature does not present a single universal hPLMM for all porous-media PDEs. What it does show, very consistently, is a reusable architecture: local subproblems posed on pore-aware neighborhoods or cells; enriched coarse representations that preserve interface, inertial, or transfer physics; and algebraic or variational couplings designed to remain accurate on non-periodic, high-contrast, or multi-continuum media. A plausible implication is that hPLMM is emerging less as a single named algorithm than as a technically coherent multiscale paradigm linking pore-level fidelity to scalable coarse computation.