Skip to content

chc.galerkin

Stability tier: experimental.

1D Galerkin finite elements with a tridiagonal (Thomas / progonka) solve.

Weak form of -u'' = f on [0,1] with u(0)=u(1)=0 and a piecewise-linear hat basis: the 1D analogue of the bilinear (Q1) elements of poisson_2d. Hat stiffness assembles to the tridiagonal stencil (1/h)[-1, 2, -1], solved by the Thomas sweep (the "progonka" kernel of Marchuk-Agoshkov projection-grid methods).

That operator is symmetric positive definite, where testing with the trial space is optimal by Cea's lemma. The second half of the module is the case where that fails: adding advection, -eps u'' + s u' = 0, makes the stencil antisymmetric and the matrix non-SPD, and above the cell Peclet number Pe = s*h/(2*eps) = 1 the same scheme oscillates. convection_diffusion_1d carries the Petrov-Galerkin test space that fixes it; the threshold, the oscillation criterion and the optimal parameter are derived in validation/convection_diffusion.mac and proved in proofs/convection_diffusion.v.

ConvectionDiffusionCurve dataclass

Evidence for the cell-Peclet threshold and for what a Petrov-Galerkin test space buys.

thomas_solve(sub, diag, sup, rhs)

Solve a tridiagonal system (Thomas / progonka). sub[0] and sup[-1] are unused.

poisson_1d(f, n)

FEM solution of -u'' = f on [0,1], u(0)=u(1)=0, on n uniform elements.

Returns the node coordinates (length n+1) and the nodal solution (with the zero boundaries).

poisson_2d(f, n)

FEM solution of -Δu = f on [0,1]^2, u=0 on the boundary, bilinear (Q1) elements.

The 2D analogue of poisson_1d: a bilinear tensor-product basis assembled element-by-element on an n x n grid; f is vectorised over the grids. Returns the 1D node coordinates and the (n+1, n+1) nodal solution.

optimal_upwind(peclet)

The SUPG parameter coth(Pe) - 1/Pe that makes 1-D convection-diffusion nodally exact.

Derived in validation/convection_diffusion.mac STEP 4: demanding that the stabilised amplification factor equal the exact exp(2*Pe) forces the effective cell Peclet to be tanh(Pe), and solving Pe/(1 + alpha*Pe) = tanh(Pe) gives this. Christie-Griffiths- Mitchell-Zienkiewicz (1976); Brooks-Hughes (1982).

coth(Pe) - 1/Pe is a difference of two quantities that both blow up as Pe -> 0, so the naive form loses every significant digit there. The series Pe/3 - Pe^3/45 takes over below the crossover, where it is already accurate to well past double precision.

convection_diffusion_exact(x, eps, s)

Exact solution of -eps u'' + s u' = 0 on [0,1] with u(0)=0, u(1)=1.

Written about the layer at x = 1 (validation/convection_diffusion.mac STEP 1d) rather than as (exp(s x/eps) - 1)/(exp(s/eps) - 1): the two are equal, but the latter overflows once s/eps passes ~700 while this form stays bounded, and the convection-dominated regime is the whole point of the problem.

convection_diffusion_1d(eps, s, n, *, alpha=0.0)

Petrov-Galerkin FEM for -eps u'' + s u' = 0 on [0,1], u(0)=0, u(1)=1.

Unlike poisson_1d, the operator here is not symmetric: the advection term s integral(u' phi_i) assembles antisymmetrically to +-s/2, so the stiffness matrix is not SPD and testing with the trial space is no longer the natural choice.

alpha selects the test space w_i = phi_i + alpha*(h/2)*phi_i'. For piecewise-linear elements the added piece contributes exactly an artificial diffusion alpha*s*h/2 (the phi_i'' term vanishes inside every element), so the assembled scheme is the central one with eps -> eps*(1 + alpha*Pe). That is the Petrov-Galerkin matrix, not a shortcut around it.

  • alpha = 0 -- Bubnov-Galerkin (test space = trial space), central differencing. Oscillates once the cell Peclet Pe = s*h/(2*eps) exceeds 1.
  • alpha = 1 -- full upwind. Never oscillates, but its amplification factor is exactly 1 + 2*Pe, the first two terms of exp(2*Pe), so it is only first-order accurate.
  • alpha = optimal_upwind(Pe) -- nodally exact.

Returns the node coordinates (length n+1) and the nodal solution including both boundaries.

convection_diffusion_certificate(*, eps=0.01, s=1.0, n_coarse=20, n_fine=200)

Exhibit the cell-Peclet instability, its threshold, and the Petrov-Galerkin cure.

Derived in validation/convection_diffusion.mac, machine-checked in proofs/convection_diffusion.v. Three schemes on the same mesh, one of which must fail:

  • Bubnov-Galerkin above Pe = 1 oscillates, because its amplification factor (1+Pe)/(1-Pe) turns negative there while the exact exp(2*Pe) never does;
  • full upwind never oscillates but its amplification is exactly 1 + 2*Pe, the first two terms of exp(2*Pe), so it is first-order;
  • the optimal alpha matches the exact factor and is nodally exact.

supg_error is reported both absolutely and in units of the working machine epsilon. The absolute number is meaningless without the dtype -- it is 6e-9 under float32 and 1e-16 under float64 -- and a threshold tuned to one of those would silently pass on the other.