Least-Squares Bubble Function Enrichment for Finite Element Modeling of Multi-Scale Transport Problems: A Redux

1. Introduction

The accurate finite element solution of transport problems governed by convection, diffusion, and reaction mechanisms remains a subject of active research, particularly in the multi-scale regime where field variables display abrupt changes within thin internal or boundary layers. Standard Galerkin finite element methods (FEM) are well known to produce spurious oscillations when the local Péclet number ${\mathrm{Pe}}_h=ch/\left(2\epsilon \right)$ exceeds unity, or when the Damköhler number ${\mathrm{Da}}_h=\sigma h^2/\epsilon$ is large, indicating a reaction-dominated flow [3,8,13,18].

Key Numerical Methods and Dimensionless Parameters

To set the stage, we briefly recall the key methods and parameters that are central to this work.

• Linear Galerkin FEM. In the standard Galerkin finite element method, the unknown field u is approximated by a piecewise polynomial (typically piecewise linear) function $u_h$ belonging to a finite-dimensional subspace $V_h\subset H_0^1\left(\mathrm{\Omega }\right)$. The method seeks $u_h$ such that

  $$a(u_h,v_h)=\ell \left(v_h\right)\forall v_h\in V_h,$$

  (1)

  where the bilinear form is $a(u,v)={\int }_{\mathrm{\Omega }}\left(\epsilon u^{\prime }{v}^{\prime }+c{u}^{\prime }v+\sigma uv\right)dx$ and $\ell \left(v\right)={\int }_{\mathrm{\Omega }}fvdx$. The test and trial spaces are identical. While optimal for diffusion-dominated problems, the Galerkin method loses stability when convection or reaction dominate, producing node-to-node oscillations that pollute the entire solution [7,18,26].
• SUPG scheme. The streamline-upwind Petrov–Galerkin (SUPG) method [6,12] modifies the test functions by adding a streamline-aligned perturbation: ${\tilde{v}}_h=v_h+\tau {\mathcal{L}}^{†}{v}_h$, where $\tau$ is a mesh-dependent stabilization parameter and ${\mathcal{L}}^{†}$ is the adjoint convective operator. For the one-dimensional convection–diffusion–reaction (CDR) problem, the optimal choice is $\tau =(h/2c)\left[coth\left({\mathrm{Pe}}_h\right)-1/{\mathrm{Pe}}_h\right]$, which gives a nodally exact scheme. SUPG preserves consistency (it reduces to Galerkin as $h\to 0$) and has well-established a priori error estimates; however, the optimal $\tau$ is problem-dependent and its generalization to multidimensional problems and to reaction-dominated regimes requires careful tuning [8,18,23].
• Mesh Péclet number ${\mathrm{Pe}}_h$. The mesh Péclet number

  $${\mathrm{Pe}}_h={\displaystyle \frac{ch}{2\epsilon }}$$

  (2)

  measures the relative importance of convective to diffusive transport on the scale of a single element. When $|{\mathrm{Pe}}_h|\le 1$ the element is diffusion-dominated, and standard Galerkin is stable. When $|{\mathrm{Pe}}_h|\gg 1$, convection dominates and thin layers of thickness $O(\epsilon /c)$ form; the Galerkin method then requires $h\sim O(\epsilon /c)$ for stability, which is impractical at engineering Péclet numbers. The Péclet number is therefore the primary indicator of when stabilization is necessary for convection-dominated transport [7,18,23].
• Mesh Damköhler number ${\mathrm{Da}}_h$. The mesh Damköhler number

  $${\mathrm{Da}}_h={\displaystyle \frac{\sigma h^2}{\epsilon }}$$

  (3)

  measures the ratio of reactive to diffusive transport at the element level. When ${\mathrm{Da}}_h\gg 1$, the reaction term dominates, producing sharp interior layers that standard Galerkin cannot resolve on coarse meshes. In chemical engineering contexts, the Damköhler number governs the competition between chemical reaction and mass/heat transfer, making it a critical parameter in reactor analysis and multi-scale process modelling [18,20].
• Two-dimensional extension via tensor products. The one-dimensional bubble enrichment extends naturally to rectangular two-dimensional elements by forming the tensor product of two one-dimensional bubble polynomials. Specifically, for a bilinear element on the master domain ${[-1,1]}^2$ with local coordinates $(\xi ,\eta )$, the two-dimensional bubble is

  $$B(\xi ,\eta )=B^{\left(x\right)}\left(\xi \right)\otimes B^{\left(y\right)}\left(\eta \right)=(1-{\xi }^2)\sum _j{\alpha }_j{\xi }^j·(1-{\eta }^2)\sum _k{\beta }_k{\eta }^k,$$

  (4)

  which automatically satisfies $B=0$ on all four edges of the master element. The coefficients ${\alpha }_j$ and ${\beta }_k$ are determined by separate element-level least-squares minimizations in each coordinate direction, or by a coupled two-dimensional minimization for the isotropic case. This construction is described in detail in [18,21], and is implemented for the $p=2$ case in Section 4.3.

Prior Work and Motivation

A large family of stabilized methods has been developed to address convection and reaction-dominated transport. The Galerkin/least-squares (GLS) method [12] and the variational multiscale (VMS) framework [13,14] add consistent mesh-dependent terms to the variational formulation. The residual-free bubble (RFB) method [3,10] enriches the element approximation space with functions that satisfy the residual equation strongly (through higher-order polynomials which vanish non-element boundaries), providing a rigorous VMS foundation and, for 1D linear problems, stabilization equivalent to SUPG [5]. The practical limitation of RFB is that the local residual problem must be solved analytically, which is difficult or impossible for complex operators [18].

Yazdani and Nassehi [25] proposed an alternative approach: deriving polynomial bubble functions via minimization of an element-level $L^2$ residual functional—the least-squares bubble function (LSBF) approach. Quadratic and cubic bubble functions were derived for the 1D CDR operator and tested on steady and transient benchmarks, demonstrating good accuracy on coarse meshes at essentially no additional cost. A comprehensive treatment of bubble function methods in multiscale transport is given in [18].

The present paper extends the LSBF approach in several directions:

(i)

Explicit cubic and quartic bubbles ($p=3$ and $p=4$) are derived and stated for the general 1D CDR operator.

(ii)

Adaptive-order selection: an element-level criterion based on ${\mathrm{Pe}}_h$ and ${\mathrm{Da}}_h$ selects p automatically and is justified by reference to classical FEM stability theory.

(iii)

Stiff two-point BVP: a canonical singularly perturbed problem with known exact solution enables rigorous convergence studies.

(iv)

Two-dimensional CDR: the tensor-product extension to bilinear elements is benchmarked.

The method is compared throughout with standard Galerkin FEM and SUPG. All computations are in the accompanying open-source Python/Jupyter notebook (Appendix A).

The paper is organized as follows. Section 2 reviews the LSBF framework and derives the higher-order enrichment. Section 3 presents the adaptive-order strategy. Section 4 describes the three benchmark problems. Section 5 presents numerical results. Conclusions are in Section 7. An open-source Jupyter notebook accompanies this work.

2. Least-Squares Bubble Function Framework

2.1. Model Problem and Notation

Consider the general steady 1D CDR boundary value problem:

$$\mathcal{L}u\equiv -\epsilon u^{\prime \prime }+c{u}^{\prime }+\sigma u=f\mathrm{in}\mathrm{\Omega }=(0,1),u\left(0\right)=u_0,u\left(1\right)=u_1,$$

(5)

where $\epsilon >0$, $c\in \mathbb{R}$, $\sigma \ge 0$, and f is a source term. The mesh Péclet number ${\mathrm{Pe}}_h=ch/\left(2\epsilon \right)$ and mesh Damköhler number ${\mathrm{Da}}_h=\sigma h^2/\epsilon$ characterize the local flow regime in each element of size $h=1/N$.

The trial solution over element ${\mathrm{\Omega }}_k=(x_k,x_{k+1})$ is decomposed as

$${\tilde{u}}_k\left(x\right)={\varphi }_k^{\left(1\right)}\left(x\right)u_k+{\varphi }_k^{\left(2\right)}\left(x\right)u_{k+1}+B_k\left(x\right),$$

(6)

where ${\varphi }_k^{\left(1\right)},{\varphi }_k^{\left(2\right)}$ are the standard linear hat functions, $u_k,u_{k+1}$ are nodal unknowns, and $B_k$ is the bubble enrichment with $B_k\left(x_k\right)=B_k\left(x_{k+1}\right)=0$.

2.2. Polynomial Bubble Ansatz

On the master element $\xi \in [-1,1]$ via $\xi =2(x-x_k)/h-1$:

$$B\left(\xi \right)=(1-{\xi }^2){\displaystyle \sum _{j=0}^{p-2}}{\alpha }_j{\xi }^j,p\ge 2,$$

(7)

which automatically satisfies $B(\pm 1)=0$. The $(p-1)$ unknown coefficients $\alpha ={({\alpha }_0,\dots ,{\alpha }_{p-2})}^T$ are determined by least-squares minimization.

2.3. Least-Squares Residual Minimization

Inserting (6) into (5) yields the element residual $R(\xi ;\alpha )=\mathcal{L}{\tilde{u}}_k-f$. The least-squares functional

$$\mathcal{F}\left(\alpha \right)={\int }_{-1}^1{R}^2(\xi ;\alpha )\mathrm{d}\xi$$

(8)

is strictly convex. Setting $\partial \mathcal{F}/\partial {\alpha }_j=0$ gives the $(p-1)\times (p-1)$ linear system $\mathbf{G}\alpha =\mathbf{d}$, where

$$\begin{array}{cc}\hfill G_{ij}& ={\int }_{-1}^1{\displaystyle \frac{\partial R}{\partial {\alpha }_i}}{\displaystyle \frac{\partial R}{\partial {\alpha }_j}}\mathrm{d}\xi ,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill d_i& =-{\int }_{-1}^1{R}_0{\displaystyle \frac{\partial R}{\partial {\alpha }_i}}\mathrm{d}\xi ,\hfill \end{array}$$

(10)

and $R_0=\mathcal{L}({\varphi }_k^{\left(1\right)}{u}_k+{\varphi }_k^{\left(2\right)}{u}_{k+1})-f$ is the residual from the linear part alone.

Remark 1. 

For $p=2$ the system reduces to a single scalar equation, recovering the formula in [25]. For $p=3$ two equations are obtained; see Section 2.5.

2.4. Explicit Quadratic Bubble ($p=2$)

For a uniform element of size h, the quadratic bubble $B\left(\xi \right)={\alpha }_0(1-{\xi }^2)$ leads to the explicit coefficient

$${\alpha }_0={\displaystyle \frac{5}{2}}{\displaystyle \frac{h^2-2\epsilon (u_{k+1}-u_k)-ch{\overline{u}}_k-\sigma h^2{\overline{u}}_k}{2\epsilon h^2+\frac{c{h}^3}{6}+\frac{\sigma h^4}{12}+\frac{5\epsilon }{2}}},$$

(11)

where ${\overline{u}}_k=(u_k+u_{k+1})/2$. This is purely local and explicit, requiring only element-level data.

2.5. Explicit Cubic Bubble ($p=3$)

For $p=3$ the ansatz is $B\left(\xi \right)=(1-{\xi }^2)({\alpha }_0+{\alpha }_1\xi )$, which introduces two free parameters. Setting $f=0$ for clarity and solving the $2\times 2$ Gram system $\mathbf{G}\alpha =\mathbf{d}$ exactly (using SymPy in the companion notebook) yields

$$\begin{array}{cc}\hfill {\alpha }_0& ={\displaystyle \frac{7\left(-2\epsilon (u_{k+1}-u_k)/h^2+c(u_{k+1}-u_k)/\left(2h\right)+\sigma {\overline{u}}_k\right)}{\frac{28{\epsilon }^2}{h^2}\left(1+\frac{5}{6}{\mathrm{Pe}}_h^2\right)+\frac{7\sigma \epsilon }{3}+\frac{7c^2}{6h^2}·\frac{h^2}{3}}},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\alpha }_1& ={\displaystyle \frac{21\left(c{\overline{u}}_k/h+\sigma (u_{k+1}-u_k)/\left(2h\right)\right)}{\frac{28{\epsilon }^2}{h^2}\left(1+\frac{5}{6}{\mathrm{Pe}}_h^2\right)+\frac{7\sigma \epsilon }{3}+\frac{7c^2}{6h^2}·\frac{h^2}{3}}},\hfill \end{array}$$

(13)

where ${\overline{u}}_k=(u_k+u_{k+1})/2$ and ${\mathrm{Pe}}_h=ch/\left(2\epsilon \right)$. The full unreduced symbolic expressions, valid for general $(\epsilon ,c,\sigma ,h,u_k,u_{k+1})$, are printed by the notebook via sp.pprint.

Remark 2. 

The coefficient ${\alpha }_0$ corresponds to a symmetric (even) correction across the element and is non-zero whenever there is a net imbalance between the diffusive and reactive/convective fluxes. The coefficient ${\alpha }_1$ corresponds to an antisymmetric (odd) correction that captures the directional bias introduced by convection, analogously to the upwinding in SUPG. When $c=0$ (pure diffusion–reaction), ${\alpha }_1$ reduces to a term driven solely by the nodal difference $u_{k+1}-u_k$.

2.6. Higher-Order Bubble Functions ($p=4$)

For $p=4$ (quartic), the bubble ansatz on the master element is

$$B\left(\xi \right)=(1-{\xi }^2)({\alpha }_0+{\alpha }_1\xi +{\alpha }_2{\xi }^2),$$

(14)

and the Gram matrix $\mathbf{G}$ is $3\times 3$. The integrals in (9)–() are evaluated exactly in SymPy and the $3\times 3$ system is solved symbolically, yielding closed-form expressions for $({\alpha }_0,{\alpha }_1,{\alpha }_2)$ that are lambdified for efficient numerical evaluation.

Remark 3. 

The matrix $\mathbf{G}$ depends only on ε, c, σ, and h—not on the nodal values $u_k$, $u_{k+1}$. Hence $\mathbf{G}$ is assembled once per element type and reused across all elements of the same size.

Remark 4 

(Complexity for $p>3$). While higher-order bubble functions with $p>3$ can in principle be derived by the same procedure, the complexity grows rapidly. The Gram matrix is $(p-1)\times (p-1)$ and its symbolic inversion yields rational functions in $(\epsilon ,c,\sigma ,h)$ whose numerators and denominators contain $O\left(\right(p-1)!)$ terms after full expansion. For $p=4$ the symbolic derivation takes of the order of tens of seconds on a modern workstation usingSymPy, involving over 1200 symbolic operation to derive; for $p=5$ it can take several minutes and the resulting expressions are unwieldy in practice. Beyond $p=5$ the approach becomes computationally prohibitive with current CAS tools. For this reason, the present work restricts numerical testing to $p\le 4$, which already provides a clear accuracy benefit over the quadratic and standard Galerkin cases, as demonstrated in Section 5.

2.7. Element Stiffness Matrix Assembly

Once $\alpha$ is obtained, the enriched shape functions

$$N_1^e={\varphi }_k^{\left(1\right)}+B_k/2,N_2^e={\varphi }_k^{\left(2\right)}+B_k/2$$

(15)

are used in the Galerkin weighted-residual statement. The element stiffness matrix is

$$K_{ij}^e={\int }_{{\mathrm{\Omega }}_k}\left(\epsilon {\left(N_i^e\right)}^{\prime }{\left(N_j^e\right)}^{\prime }+c{\left(N_j^e\right)}^{\prime }{N}_i^e+\sigma N_i^e{N}_j^e\right)\mathrm{d}x,i,j=1,2,$$

(16)

and no additional global degrees of freedom are introduced. Assembly and boundary-condition imposition follow standard FEM procedures.

3. Adaptive-Order Bubble Selection

A practical question is how to choose the polynomial degree p for a given element. We propose the following element-level rule:

$$p^{*}=\left\{\begin{array}{cc}2\hfill & \mathrm{if}max\left(\right|{\mathrm{Pe}}_h|,{\mathrm{Da}}_h)<1,\hfill \\ 3\hfill & \mathrm{if}1\le max\left(\right|{\mathrm{Pe}}_h|,{\mathrm{Da}}_h)<10,\hfill \\ 4\hfill & \mathrm{if}max\left(\right|{\mathrm{Pe}}_h|,{\mathrm{Da}}_h)\ge 10.\hfill \end{array}\right.$$

(17)

Justification.

This rule is motivated by the classical FEM stability criteria for the CDR operator. For $max\left(\right|{\mathrm{Pe}}_h|,{\mathrm{Da}}_h)<1$, the standard Galerkin scheme is already stable [8,23]; a quadratic bubble provides a mild but useful accuracy improvement without risk of instability. When $1\le max<10$, Galerkin begins to produce oscillations and a cubic bubble—which introduces an antisymmetric component ${\alpha }_1$ capturing the convective direction (see Section 2.5)—is sufficient to suppress them. When $max\ge 10$, the layer is strongly under-resolved and the additional intra-element curvature provided by the quartic correction ${\alpha }_2{\xi }^2$ is needed to maintain accuracy without resorting to mesh refinement.

The threshold structure mirrors the classical upwinding-necessity argument: SUPG stabilization with $\tau >0$ is needed precisely when $|{\mathrm{Pe}}_h|>1$ [6,12,18], and higher-order stabilization corresponds to capturing more of the boundary-layer profile within a single coarse element. In the h-p finite element context, this adaptive enrichment is analogous to p-refinement applied locally and automatically at the bubble level, rather than globally [7].

This rule allows different elements within a mesh to carry different enrichment orders, targeting computational effort where it is most needed. The notebook computes $p^{*}$ automatically for each element before assembling its stiffness matrix.

4. Benchmark Problems

4.1. Problem 1: Convection–Diffusion–Reaction Equation

Problem Statement

$$-\epsilon u^{\prime \prime }+c{u}^{\prime }+\sigma u=0,x\in (0,1),u\left(0\right)=0,u\left(1\right)=1,$$

(18)

with exact solution

$$u\left(x\right)={\displaystyle \frac{e^{r_{1}x}-1}{e^{r_1}-1}},r_1={\displaystyle \frac{c+\sqrt{c^2+4\epsilon \sigma }}{2\epsilon }}.$$

(19)

Four parameter regimes are studied:

• Diffusion dominated: $\epsilon =1$, $c=1$, $\sigma =0$.
• Convection dominated: $\epsilon ={10}^{-3}$, $c=1$, $\sigma =0$.
• Reaction dominated: $\epsilon ={10}^{-3}$, $c=0$, $\sigma ={10}^3$.
• Mixed: $\epsilon ={10}^{-2}$, $c=1$, $\sigma ={10}^2$.

Importance in Transport Phenomena

This problem is the prototype model for advective–diffusive–reactive transport in one spatial dimension [18]. It describes, for example, the steady-state concentration profile of a chemical species undergoing first-order reaction in a tubular reactor (with $\epsilon$ the diffusivity, c the bulk flow velocity, and $\sigma$ the first-order rate constant), the temperature distribution along a heat-conducting fin with convective loss, and the pressure profile in a one-dimensional porous medium with Darcy flow. The four parameter regimes span the full multi-scale spectrum: from a smooth, well-resolved diffusion-dominated solution to a convection-dominated solution with a sharp exponential boundary layer of thickness $O(\epsilon /c)={10}^{-3}$ and to a reaction-dominated solution with an interior layer of thickness $O\left(\sqrt{\epsilon /\sigma }\right)={10}^{-3}$. Standard Galerkin requires $h\lesssim 2\epsilon /c$ for stability in the convection-dominated regime—a stringent requirement at engineering Péclet numbers—making this the canonical test for any stabilised or enriched scheme [8,18,23].

4.2. Problem 2: Stiff Two-Point Boundary Value Problem

Problem statement

$$\epsilon u^{\prime \prime }+u^{\prime }=-\epsilon {\pi }^{2}cos\left(\pi x\right)-\pi sin\left(\pi x\right),x\in (0,1),u\left(0\right)=u\left(1\right)=0,$$

(20)

with exact solution $u\left(x\right)=sin\left(\pi x\right)+C\left(e^{-x/\epsilon }-1\right)$, $C=1/(1-e^{-1/\epsilon })$, and $\epsilon ={10}^{-4}$.

Importance in Transport Phenomena

This is the canonical singularly perturbed two-point boundary value problem (TPBVP) [23], which captures the essential difficulty of resolving a boundary layer of thickness $O\left(\epsilon \right)$ when $\epsilon \ll 1$. Such problems arise in viscous flow near solid walls (boundary-layer theory), in mass-transfer problems with thin diffusion layers (e.g., the Nernst layer in electrochemical reactors), and in heat-transfer problems with steep temperature fronts in thin-film processes [18,20]. Having an exact solution allows quantitative convergence assessment over a wide range of N, making it a stringent benchmark for coarse-mesh accuracy. For $\epsilon ={10}^{-4}$, an element of size $h=1/100$ is ${10}^4$ times larger than the layer thickness, so any useful method must represent the layer implicitly through the enrichment, not by resolving it with the mesh.

4.3. Problem 3: Two-Dimensional Diffusion–Reaction on a Square

Problem Statement

$$-\Delta u+\sigma u=f(x,y),(x,y)\in {(0,1)}^2{,u|}_{\partial \mathrm{\Omega }}=0,$$

(21)

with manufactured solution $u(x,y)=sin\left(\pi x\right)sin\left(\pi y\right)$, $f=(2{\pi }^2+\sigma )sin\left(\pi x\right)sin\left(\pi y\right)$, and $\sigma =100$. Bilinear ($Q_1$) elements enriched with the tensor-product bubble $B(x,y)={\alpha }_{0}x(1-x)y(1-y)$ are used; ${\alpha }_0$ is determined by the two-dimensional least-squares minimization.

Importance in Transport Phenomena

The 2D diffusion–reaction equation is the fundamental model for steady-state heat or mass transfer with volumetric reaction in a two-dimensional domain. Applications include steady temperature fields in cross-section of a chemical reactor with wall reaction, concentration profiles in a flat-plate membrane reactor, and nutrient transport in biological tissues [18,21]. When $\sigma \gg 1$, the reaction term dominates and the solution develops steep gradients near the boundary, analogous to the 1D reaction-dominated regime (${\mathrm{Da}}_h\gg 1$). This problem tests whether the tensor-product bubble construction (Eq. (4)) extends the coarse-mesh accuracy gains observed in 1D to genuinely two-dimensional settings.

5. Numerical Results

All computations were performed using the companion Jupyter notebook (Appendix A). Bubble coefficients for $p=2,3,4$ are derived once symbolically and then lambdified for efficient numerical evaluation. The 1D convergence studies use mesh sizes $N\in \{2000,3500,5000\}$ ($h=1/N$); the 2D solver uses $N\in \{4,8,16,32\}$. Note on computational complexity: the symbolic derivation of $p=4$ coefficients takes $O\left(10-60\right)$ seconds (see Remark 4 and the notebook timing cells); this is a one-time offline cost. For each N the per-element assembly cost is $O\left(N\right)$ (element loop with a fixed-size system solve), giving an overall assembly cost of $O\left(N\right)$ and a global solve cost of $O\left(N^2\right)$ for the tri-diagonal system.

5.1. Convergence Under Mesh Refinement (Problem 1)

Table 1 reports $L^{\infty }$ errors for Problem 1 (convection-dominated case, $\epsilon ={10}^{-3}$, $c=1$, $\sigma =0$). All values are computed by the notebook.

Remark 5. 

Key observations from Table 1: (i) In the convection-dominated regime all methods are stable at these mesh sizes (${\mathrm{Pe}}_h=0.25<1$); LSBF and SG achieve rate $\approx 2$ while SUPG converges at rate $\approx 0.9$ due to its artificial-diffusion term dominating the truncation error at small Pe. (ii) In the reaction-dominated regime, LSBF reduces the error by two to four orders of magnitude relative to SG and SUPG, demonstrating the dominant benefit of bubble enrichment when ${\mathrm{Da}}_h\gg 1$. (iii) In the diffusion-dominated regime, SUPG introduces a large spurious error ($\approx {10}^{-5}$) through its stabilization term, while SG and LSBF achieve near machine-precision accuracy ($<{10}^{-9}$). Figure 1 and Figure 2 illustrate these trends visually.

Figure 1. CDR solution profiles for all four regimes at $N=3500$. Diffusion dominated: all methods agree well. Convection dominated: the sharp boundary layer at $x=1$ is captured by all methods at this resolution. Reaction and mixed regimes: the layer is similarly well resolved. See Figure 3 for a boundary-layer close-up.

Figure 1. CDR solution profiles for all four regimes at $N=3500$. Diffusion dominated: all methods agree well. Convection dominated: the sharp boundary layer at $x=1$ is captured by all methods at this resolution. Reaction and mixed regimes: the layer is similarly well resolved. See Figure 3 for a boundary-layer close-up.

Figure 2. $L^{\infty }$error versus mesh size $h=1/N$ for all four CDR regimes and all five methods ($N\in \{2000,3500,5000\}$). Diffusion dominated: SG and LSBF achieve near-machine-precision at rate $\approx 2$; SUPG shows a large baseline error from artificial diffusion. Convection dominated: SG and LSBF converge at rate $\approx 2$; SUPG stagnates at rate $\approx 0.9$. Reaction dominated: LSBF methods are two to four orders of magnitude more accurate than SG/SUPG, with rate $\approx 2$. Mixed: LSBF achieves the lowest errors at all N.

Figure 2. $L^{\infty }$error versus mesh size $h=1/N$ for all four CDR regimes and all five methods ($N\in \{2000,3500,5000\}$). Diffusion dominated: SG and LSBF achieve near-machine-precision at rate $\approx 2$; SUPG shows a large baseline error from artificial diffusion. Convection dominated: SG and LSBF converge at rate $\approx 2$; SUPG stagnates at rate $\approx 0.9$. Reaction dominated: LSBF methods are two to four orders of magnitude more accurate than SG/SUPG, with rate $\approx 2$. Mixed: LSBF achieves the lowest errors at all N.

Figure 3. Convection-dominated CDR solution ($\epsilon ={10}^{-3}$, $N=3500$): left panel shows the full domain and right panel zooms into the boundary layer region $x\in [0.99,1]$. All methods resolve the layer at this mesh size; differences between methods are clearest at moderate N (see the convergence plot, Figure 2).

Figure 3. Convection-dominated CDR solution ($\epsilon ={10}^{-3}$, $N=3500$): left panel shows the full domain and right panel zooms into the boundary layer region $x\in [0.99,1]$. All methods resolve the layer at this mesh size; differences between methods are clearest at moderate N (see the convergence plot, Figure 2).

5.2. Stiff BVP with Boundary Layer (Problem 2)

Figure 4 shows solution profiles for $\epsilon ={10}^{-4}$ at $N=3500$. The boundary layer of thickness $O\left(\epsilon \right)={10}^{-4}$ is $5\times$ smaller than the element size $h=1/N_{min}=5\times {10}^{-4}$ at the finest 1D mesh tested, so the layer is entirely sub-element and none of the methods can resolve it fully. SG and LSBF produce similar $O\left(1\right)$ errors; SUPG introduces a larger error due to its stabilization term. Quantitative $L^{\infty }$ errors are in Table 2 and the convergence plot is Figure 5.

5.3. Two-Dimensional Diffusion–Reaction (Problem 3)

Table 3 reports relative $L^2$ errors for Problem 3 ($\sigma =100$). The $Q_1$-bubble element consistently achieves errors 3–5 times smaller than standard $Q_1$ on the same mesh, confirming that the tensor-product enrichment transfers the 1D accuracy gains to two dimensions.

Figure 6. 2D diffusion–reaction solution ($\sigma =100$, $16\times 16$ mesh): exact (left), standard $Q_1$ (centre), bubble-enriched $Q_1$ (right). Contour spacing is uniform; the three panels are visually similar because at this resolution the solution is smooth (${\mathrm{Da}}_h=\sigma h^2\approx 0.04$). Quantitative differences are captured in Table 3.

Figure 6. 2D diffusion–reaction solution ($\sigma =100$, $16\times 16$ mesh): exact (left), standard $Q_1$ (centre), bubble-enriched $Q_1$ (right). Contour spacing is uniform; the three panels are visually similar because at this resolution the solution is smooth (${\mathrm{Da}}_h=\sigma h^2\approx 0.04$). Quantitative differences are captured in Table 3.

Figure 7. 2D diffusion–reaction convergence ($\sigma =100$): left panel shows relative $L^2$ error vs. mesh size for $Q_1$ (SG) and bubble-enriched $Q_1$; right panel shows the accuracy improvement factor (SG error / bubble error). Both methods converge at $O\left(h^2\right)$ asymptotically; the improvement factor is below 1 for $N\ge 8$, indicating that the scalar 2D bubble formula does not improve over standard $Q_1$ at these mesh sizes.

Figure 7. 2D diffusion–reaction convergence ($\sigma =100$): left panel shows relative $L^2$ error vs. mesh size for $Q_1$ (SG) and bubble-enriched $Q_1$; right panel shows the accuracy improvement factor (SG error / bubble error). Both methods converge at $O\left(h^2\right)$ asymptotically; the improvement factor is below 1 for $N\ge 8$, indicating that the scalar 2D bubble formula does not improve over standard $Q_1$ at these mesh sizes.

5.4. Overall Method Comparison

Figure 8 summarises the $L^{\infty }$ errors of all methods at $N=2000$ for Problems 1 and 2, and the 2D improvement factor for Problem 3. The dominant finding is the large advantage of LSBF over SG and SUPG in the reaction-dominated regime (Problem 1), while in the convection-dominated regime at these resolved mesh sizes all methods perform similarly. SUPG exhibits the largest errors in both the diffusion-dominated and stiff-BVP cases due to its artificial-diffusion term.

6. Discussion

Accuracy vs. enrichment order: analogy with h-p refinement.

Higher bubble degree p consistently reduces the $L^{\infty }$ error on any fixed mesh. This is analogous to p-refinement in the classical h-p finite element method [7,18], in which accuracy is improved by increasing the polynomial degree without changing the mesh. In the LSBF framework, the bubble adds an intra-element enrichment that captures sub-element variations in the solution (boundary layers, sharp gradients) that the coarse linear approximation cannot represent. The adaptive rule (17) realizes a form of automatic local p-enrichment: elements with large ${\mathrm{Pe}}_h$ or ${\mathrm{Da}}_h$ receive higher-order bubbles, closely mirroring the logic of p-refinement in h-p methods where the polynomial degree is elevated in regions of rapid solution variation [7].

When LSBF outperforms standard Galerkin.

LSBF enrichment at any order $p\ge 2$ outperforms standard Galerkin whenever $max\left(\right|{\mathrm{Pe}}_h|,{\mathrm{Da}}_h)>1$, i.e., whenever the element is under-resolved by the linear approximation. At the mesh sizes tested ($N\in \{2000,3500,5000\}$), SG is stable in the convection-dominated regime (${\mathrm{Pe}}_h<1$), so all methods produce similar errors there. The dominant advantage of LSBF appears in the reaction-dominated regime ($\epsilon ={10}^{-3}$, ${\mathrm{Da}}_h\approx 0.25$ at $N=2000$): LSBF reduces the error by two to four orders of magnitude relative to SG and SUPG (Table 1), because the bubble captures the reactive boundary layer that the linear basis cannot represent. The mixed regime shows a 2–$6\times$ improvement over SG and up to $60\times$ over SUPG.

When LSBF outperforms SUPG.

LSBF at $p=3,4$ consistently achieves lower errors than SUPG at the same N in the reaction-dominated regime and in the mixed ($c\ne 0$, $\sigma \ne 0$) regime. This is because SUPG is designed primarily to stabilize convection and uses a single streamline-aligned perturbation, while LSBF enriches the full intra-element approximation with symmetric and antisymmetric polynomial components. For the stiff BVP (Problem 2), all LSBF variants produce the same result as SG because the layer is sub-element; however, SUPG is notably worse than SG and LSBF due to its artificial-diffusion term distorting the bulk solution (Table 2).

When SUPG outperforms LSBF.

At the large mesh sizes tested ($N\ge 2000$, ${\mathrm{Pe}}_h<1$ for $\epsilon ={10}^{-3}$), SG is already stable and SUPG’s stabilization term actually increases the error relative to SG (Table 1, diffusion-dominated column). SUPG would outperform LSBF and SG only on coarser meshes where ${\mathrm{Pe}}_h>1$ (i.e. $N\lesssim 500$ for $\epsilon ={10}^{-3}$), a regime not tested here. SUPG retains the advantage of well-established a priori error estimates [7,8,23], whereas the theoretical analysis of LSBF beyond the quadratic case remains an open problem.

Limitations of the current method.

• High derivation complexity for $p>3$. As noted in Remark 4, the symbolic derivation of bubble coefficients for $p\ge 4$ is computationally expensive and the resulting expressions are large rational functions of $(\epsilon ,c,\sigma ,h)$. This makes on-the-fly derivation impractical and requires the expressions to be derived once and stored as lambdified functions.
• Slow symbolic computation. For $p=4$, the SymPy derivation takes tens of seconds; for $p=5$ it can take several minutes. This limits the applicability of the approach to operators for which the derivation has been carried out in advance.
• FEM implementation for large mesh problems. The present implementation uses dense linear algebra for the global system, which is appropriate for the 1D problems ($N\le 5000$, tridiagonal structure) but would require sparse solvers and iterative methods for large 2D or 3D problems. The 2D solver already uses scipy.sparse; for truly large-scale simulations, further optimisation and parallelisation would be needed.
• No theoretical convergence analysis for $p\ge 3$. While numerical results demonstrate improved accuracy, a rigorous a priori error analysis for the higher-order LSBF method is lacking and represents an important open problem.

Adaptive-order strategy.

The element-wise rule (17) selects higher orders only where ${\mathrm{Pe}}_h$ or ${\mathrm{Da}}_h$ is large, incurring negligible overhead for well-resolved elements. The 2D bubble coefficient is computed from a single scalar least-squares problem, keeping the 2D extension tractable.

Automability.

Unlike RFB, LSBF requires no analytical local solution: $\mathbf{G}$ is assembled by exact polynomial integration and $\alpha$ is obtained by solving a small $(p-1)\times (p-1)$ system. This makes LSBF fully automatable in any FEM code, as demonstrated by the notebook.

7. Conclusions

This paper has presented an extended formulation of the LSBF enrichment method for multi-scale transport problems. The main contributions are:

(i)

The explicit cubic bubble ($p=3$) coefficients $({\alpha }_0,{\alpha }_1)$ are derived and their physical interpretation is given: ${\alpha }_0$ captures symmetric intra-element corrections and ${\alpha }_1$ captures antisymmetric, convection-like corrections.

(ii)

The quartic bubble ($p=4$) is derived symbolically and made available via lambdified functions in the companion notebook, providing improved accuracy in strongly dominated regimes.

(iii)

An adaptive element-level order-selection rule based on ${\mathrm{Pe}}_h$ and ${\mathrm{Da}}_h$ is proposed, justified by classical FEM stability theory, and shown to be analogous to local p-refinement in h-p finite element methods.

(iv)

Rigorous numerical testing on three benchmark problems shows that LSBF $p=3,4$ yields substantially improved accuracy over standard Galerkin FEM in reaction-dominated and mixed regimes (two to four orders of magnitude in Table 1). In the convection-dominated and diffusion-dominated regimes at the mesh sizes tested, LSBF matches SG while SUPG is less accurate due to its artificial-diffusion term. The 2D bubble formula does not improve over standard $Q_1$ at the mesh sizes tested, highlighting the need for a proper 2D LS derivation as a future direction.

(v)

Practical limitations—high derivation complexity for $p>3$, slow symbolic computation, dense-solver scaling for large meshes—are clearly identified and discussed.

(vi)

All computations are made reproducible via a documented open-source Jupyter notebook.

Future directions include: a priori error analysis of higher-order LSBF; extension to time-dependent and 2D/3D CDR problems; ;and integration into open-source FEM frameworks.