experimental/adrh: boundary kinds (prescribed inflow flux, weak Dirichlet IIPG/NIPG/SIPG) against Proteus

experimental/adrh/README.md

@@ -29,8 +29,24 @@ against 0.39 s for the Python path's kernels alone with its quadrature arrays al
   coupled system 2.1x, Poisson 1.3x; with all-true masks the unrolled kernel costs the same as the runtime loops, so the gain
   is all sparsity.
 * `ad_regression.py` – steady advection-diffusion with the manufactured sine solution
-  (`AnalyticalSolutions.LinearADR_Sine`, c = 0) on regular tets: the engine, a plain Proteus model (`ad_sine_p.py`,
-  `ad_sine_n.py`) and the exact solution. Engine and Proteus nodal solutions agree to 4e-15; L2 error rates 1.85 and 1.97.
-  The driver imposes strong Dirichlet rows; boundary conditions are not in the engine yet.
+  (`AnalyticalSolutions.LinearADR_Sine`, c = 0) on regular triangles (`AD_ND=2`) or tets (`AD_ND=3`): the engine, a plain
+  Proteus model (`ad_sine_p.py`, `ad_sine_n.py`) and the exact solution. Engine and Proteus nodal solutions agree to solver
+  round-off (< 5e-14) in every boundary mode and both dimensions; L2 error rates 1.85-2.00.
 
-Design, decisions and numbers: `ADRH_ENGINE_DESIGN.md` (kept beside the campaign notes, sections 11-13).
+## Boundary conditions (`boundary.h`)
+
+`exteriorFaceKernel` adds the exterior-face terms per face (residual and Jacobian, scattered by the driver), for
+segments of triangles and triangles of tetrahedra (`Coefficients::nSpace` = 2 or 3). A BC functor says, per component
+and boundary quadrature point, which of three kinds applies; each is the engine's version of a Proteus exterior
+numerical-flux treatment, and `ad_regression.py` compares against that treatment (`AD_BC=`, in 2D and 3D):
+
+| kind | what it does | Proteus equivalent | `AD_BC` |
+|---|---|---|---|
+| `Strong` | nothing on the face; the driver replaces the row | `StrongDirichletFactory`, `numericalFluxType=None` | `strong` |
+| `Flux` | prescribed diffusive flux, optionally affine in the trace (`gdiff + dgdiff u`); advective flux prescribed where given, else the interior `f(u).n` on outflow and zero on inflow | `StrongDirichletFactory({ci: 'outFlow'})` with `diffusiveFluxBoundaryConditions` / `advectiveFluxBoundaryConditions`; `RobinFactory` for `alpha u + beta q.n = g` | `mixed`, `inflow`, `robin` |
+| `WeakDirichlet` | upwinded advective flux against `f(g)`, `-a grad phi . n + (C/h^p)(u - g)`, adjoint `-sigma (u - g) (a^T n).grad w` | `Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior` (C = 2, 10, 100; sigma = 0, -1, +1) | `iipg`, `nipg`, `sipg` |
+
+The header's comment says where the engine deliberately differs from Proteus (the adjoint's tensor contraction, da/du in
+the Jacobian); neither shows with an isotropic constant diffusion tensor.
+
+Design, decisions and numbers: `ADRH_ENGINE_DESIGN.md` (kept beside the campaign notes, sections 11-15).

experimental/adrh/ad_regression.py

@@ -1,16 +1,26 @@
 """Advection-diffusion regression: ADRH engine vs the Proteus model (ad_sine_p/n.py) vs the exact solution.
 
-Steady div(b u - kappa grad u) + d(x) = 0 on the unit cube, exact solution sin(omega.x + omega0) (the
-verified AnalyticalSolutions.LinearADR_Sine with c = 0), P1 Galerkin on regular tetrahedra, strong
-Dirichlet data from the exact solution. The engine builds the element residual and Jacobian; the driver
-assembles, imposes the Dirichlet rows and solves with scipy. Checks:
-
-  1. engine and Proteus solve the same discrete problem: nodal solutions agree to solver round-off;
-     AD_MIXED=1 makes the face x = 1 a flux face (prescribed diffusive flux, advective outflow) instead of Dirichlet;
+Steady div(b u - kappa grad u) + d(x) = 0 on the unit square or cube, exact solution sin(omega.x + omega0) (the
+verified AnalyticalSolutions.LinearADR_Sine with c = 0), P1 Galerkin on regular triangles (AD_ND=2) or tetrahedra
+(AD_ND=3, the default). The engine builds the element and exterior-face residual and Jacobian; the driver
+assembles, imposes any strong Dirichlet rows and solves with scipy. Checks:
+
+  1. engine and Proteus solve the same discrete problem: nodal solutions agree to solver round-off.
+     AD_BC picks the boundary treatment on both sides (ad_sine_p.py has the Proteus side):
+       strong (default)  Dirichlet data everywhere, imposed strongly by the driver (Proteus: numericalFluxType=None)
+       mixed             the face x = 1 is a flux face: prescribed diffusive flux, advection out by the interior value
+                         (Proteus: StrongDirichletFactory({0: 'outFlow'}); AD_MIXED=1 means this)
+       inflow            as mixed, and the inflow face x = 0 has the advective flux (b . n) u and the diffusive flux
+                         prescribed (the engine's Flux kind with hasAdv)
+       robin             as mixed, but the face x = 1 carries alpha u + beta q.n = g (the engine's Flux kind with an
+                         affine diffusive flux, against Proteus's RobinFactory)
+       iipg, nipg, sipg  Dirichlet data everywhere, imposed weakly: upwinded advective flux against f(g), diffusive flux
+                         with penalty C/h (2, 10, 100) and the NIPG/SIPG adjoint term (the engine's WeakDirichlet kind
+                         against Proteus's Advection_DiagonalUpwind_Diffusion_*_exterior);
   2. both converge to the exact solution at the P1 rates (L2 ~ h^2, max-norm ~ h^2);
   3. the exact solution fed through the engine's residual is small and shrinks with h (consistency).
 
-usage: python ad_regression.py [nn ...]        default 5 9 17
+usage: [AD_ND=2|3] [AD_BC=mode] python ad_regression.py [nn ...]        default 5 9 17
 """
 import os
 import sys
@@ -27,16 +37,29 @@
 from proteus import default_s, default_so, NumericalSolution, Quadrature
 from proteus.defaults import load_numerics, load_physics
 from proteus.iproteus import opts as runopts
-from proteus import AnalyticalSolutions
 
-b = np.array([1.0, 0.5, 0.25]); kappa = 0.2; omega = np.array([3.0, 2.0, 4.0]); omega0 = 0.3
+ND = int(os.environ.get('AD_ND', '3'))
+os.environ['AD_ND'] = str(ND)
+nD = ND + 1                                                       # nodes per element
+b = np.array([1.0, 0.5, 0.25])[:ND]; kappa = 0.2; omega = np.array([3.0, 2.0, 4.0])[:ND]; omega0 = 0.3
 order = 4
 os.environ.setdefault('AD_QORDER', str(order))
-quad = Quadrature.SimplexGaussQuadrature(3, order)
-qp = np.array([list(p)[:3] for p in quad.points]); qw = np.array(quad.weights)
-# verified exact solution (LinearADR_Sine with c = 0): same u as the functor's source assumes
-exact_cls = AnalyticalSolutions.LinearADR_Sine(b=b, a=kappa * np.eye(3), c=0.0, omega=omega, omega0=omega0)
-exact = lambda X: np.sin(X @ omega + omega0)
+quad = Quadrature.SimplexGaussQuadrature(ND, order)
+qp = np.array([list(p)[:ND] for p in quad.points]); qw = np.array(quad.weights)
+quadf = Quadrature.SimplexGaussQuadrature(ND - 1, order)
+qpf = np.array([list(p)[:ND - 1] for p in quadf.points]); qwf = np.array(quadf.weights)
+# the verified exact solution (LinearADR_Sine with c = 0), on points with any leading shape and 2 or 3 coordinates
+exact = lambda X: np.sin(X[..., :ND] @ omega + omega0)
+
+BC = os.environ.get('AD_BC', 'mixed' if os.environ.get('AD_MIXED', '0') == '1' else 'strong')
+WEAK = {'iipg': (2.0, 1.0, 0.0), 'nipg': (10.0, 1.0, -1.0), 'sipg': (100.0, 1.0, 1.0)}   # C, p, sigma as in Proteus
+assert BC in ('strong', 'mixed', 'inflow', 'robin') or BC in WEAK, BC
+os.environ['AD_BC'] = BC
+ROBIN_ALPHA, ROBIN_BETA = 1.0, 2.0                                 # as ad_sine_p.py
+# the faces the engine's boundary kernel sees, by plane: (kind code of SineBC, axis, value); the rest is strong Dirichlet
+FACES = {'strong': [], 'mixed': [(1, 0, 1.0)], 'inflow': [(2, 0, 0.0), (1, 0, 1.0)], 'robin': [(4, 0, 1.0)]}
+for mode in WEAK:
+    FACES[mode] = [(3, ax, v) for ax in range(ND) for v in (0.0, 1.0)]
 
 
 def proteus_solution(nn):
@@ -53,33 +76,47 @@ def proteus_solution(nn):
     return np.array(lm.u[0].dof), np.array(lm.mesh.nodeArray), np.array(lm.mesh.elementNodesArray)
 
 
-MIXED = os.environ.get('AD_MIXED', '0') == '1'
-quadf = Quadrature.SimplexGaussQuadrature(2, order)
-qpf = np.array([list(p)[:2] for p in quadf.points]); qwf = np.array(quadf.weights)
+def plane_faces(nodes, conn, axis, value):
+    """Element faces (element, local face = the opposite node) whose nodes all lie on the plane x[axis] = value.
+    Per plane, because a face with its nodes on two different planes is an interior face."""
+    on = np.abs(nodes[conn][:, :, axis] - value) < 1e-12                       # (nE, nD)
+    return np.nonzero(on.sum(axis=1)[:, None] - on.astype(int) == ND)
 
 
 def engine_solution(nodes, conn):
     nE, nN = len(conn), len(nodes)
-    xn = np.ascontiguousarray(nodes[conn])
-    u0 = np.zeros((nE, 1, 4))
-    R, J = adrh_proto.advection_diffusion_sine(b, kappa, omega, omega0, xn, u0, qp, qw)
-    rows = np.repeat(conn, 4, axis=1).ravel(); cols = np.tile(conn, (1, 4)).ravel()
+    xn = np.ascontiguousarray(nodes[conn])                                      # (nE, nD, 3): Proteus keeps 3 coordinates
+    u0 = np.zeros((nE, 1, nD))
+    R, J = adrh_proto.advection_diffusion_sine(ND, b, kappa, omega, omega0, xn, u0, qp, qw)
+    rows = np.repeat(conn, nD, axis=1).ravel(); cols = np.tile(conn, (1, nD)).ravel()
     K = sp.coo_matrix((J[:, 0, :, 0, :].ravel(), (rows, cols)), shape=(nN, nN)).tocsr()
     rhs = -np.bincount(conn.ravel(), weights=R[:, 0, :].ravel(), minlength=nN)
-    if MIXED:
-        # strong Dirichlet on every boundary node except those strictly inside the face x = 1
-        onb = (nodes.min(axis=1) < 1e-12) | (nodes[:, 1:].max(axis=1) > 1 - 1e-12)
-        # flux faces: element faces whose three nodes all lie on x = 1; local face k is opposite node k
-        onx1 = nodes[conn][:, :, 0] > 1 - 1e-12                              # (nE, 4)
-        fe, fk = np.nonzero(onx1.sum(axis=1)[:, None] - onx1.astype(int) == 3)
+    if FACES[BC]:
+        fe, fk, kind = [], [], []
+        for code, ax, v in FACES[BC]:
+            e, k = plane_faces(nodes, conn, ax, v)
+            fe.append(e); fk.append(k); kind.append(np.full(len(e), code))
+        fe, fk, kind = (np.concatenate(a) for a in (fe, fk, kind))
+        pc, pp, sig = WEAK.get(BC, (2.0, 1.0, 0.0))
         face_conn = conn[fe]
-        Rf, Jf = adrh_proto.advection_diffusion_sine_flux_faces(b, kappa, omega, omega0, np.ascontiguousarray(fk, dtype='i'),
-                                                                np.ascontiguousarray(nodes[face_conn]), np.zeros((len(fe), 1, 4)), qpf, qwf)
-        rows_f = np.repeat(face_conn, 4, axis=1).ravel(); cols_f = np.tile(face_conn, (1, 4)).ravel()
+        Rf, Jf = adrh_proto.advection_diffusion_sine_faces(ND, b, kappa, omega, omega0, np.ascontiguousarray(fk, dtype='i'),
+                                                           np.ascontiguousarray(kind, dtype='i'),
+                                                           np.ascontiguousarray(nodes[face_conn]), np.zeros((len(fe), 1, nD)),
+                                                           qpf, qwf, pc, pp, sig, ROBIN_ALPHA, ROBIN_BETA)
+        rows_f = np.repeat(face_conn, nD, axis=1).ravel(); cols_f = np.tile(face_conn, (1, nD)).ravel()
         K = K + sp.coo_matrix((Jf[:, 0, :, 0, :].ravel(), (rows_f, cols_f)), shape=(nN, nN)).tocsr()
         rhs = rhs - np.bincount(face_conn.ravel(), weights=Rf[:, 0, :].ravel(), minlength=nN)
+    # the nodes whose rows the driver replaces by the Dirichlet value
+    X = nodes[:, :ND]
+    lo, hi = X.min(axis=1) < 1e-12, X.max(axis=1) > 1 - 1e-12
+    if BC in WEAK:
+        onb = np.zeros(nN, bool)
+    elif BC in ('mixed', 'robin'):
+        onb = lo | (X[:, 1:].max(axis=1) > 1 - 1e-12)
+    elif BC == 'inflow':
+        onb = (X[:, 1:].min(axis=1) < 1e-12) | (X[:, 1:].max(axis=1) > 1 - 1e-12)
     else:
-        onb = (nodes.min(axis=1) < 1e-12) | (nodes.max(axis=1) > 1 - 1e-12)
+        onb = lo | hi
     ub = exact(nodes)
     free = ~onb
     Kff = K[free][:, free]
@@ -93,24 +130,26 @@ def errors(u, nodes, conn, xn):
     lam = np.column_stack([1 - qp.sum(axis=1), qp])
     uh = u[conn] @ lam.T                                              # (nE, nQ)
     xq = np.einsum('qb,ebi->eqi', lam, xn)
-    d = np.einsum('ebi,bj->eij', xn[:, 1:, :] - xn[:, :1, :], np.eye(3))  # edge matrix (nE, 3, 3) rows = edges
-    vol = np.abs(np.linalg.det(xn[:, 1:, :] - xn[:, :1, :]))
+    vol = np.abs(np.linalg.det(xn[:, 1:, :ND] - xn[:, :1, :ND]))      # the reference weights carry the 1/d!
     l2 = np.sqrt(np.sum(vol[:, None] * qw[None, :] * (uh - exact(xq)) ** 2))
     return l2, np.abs(u - exact(nodes)).max()
 
 
 sizes = [int(a) for a in sys.argv[1:]] or [5, 9, 17]
 prev = None
+print("AD_ND=%d AD_BC=%s" % (ND, BC))
 print("%4s %9s %12s %12s %14s %14s %9s %9s %14s" % ("nn", "elements", "|eng-prot|", "eng-vs-ex L2", "prot-vs-ex L2", "eng max err", "L2 rate", "max rate", "res(exact) inf"))
 for nn in sizes:
     up, nodes, conn = proteus_solution(nn)
     ue, xn, R, J = engine_solution(nodes, conn)
     l2e, mxe = errors(ue, nodes, conn, xn)
     l2p, mxp = errors(up, nodes, conn, xn)
-    # the exact solution through the engine: element residual of u_exact (interior elements only, no boundary terms)
+    # the exact solution through the engine: element residual of u_exact at the interior nodes (no boundary terms)
     uex = exact(nodes)[conn][:, None, :]
-    Rx = adrh_proto.advection_diffusion_sine(b, kappa, omega, omega0, xn, np.ascontiguousarray(uex), qp, qw)[0]
-    resid = np.abs(np.bincount(conn.ravel(), weights=Rx[:, 0, :].ravel(), minlength=len(nodes))[(nodes.min(axis=1) > 1e-12) & (nodes.max(axis=1) < 1 - 1e-12)]).max()
+    Rx = adrh_proto.advection_diffusion_sine(ND, b, kappa, omega, omega0, xn, np.ascontiguousarray(uex), qp, qw)[0]
+    X = nodes[:, :ND]
+    interior = (X.min(axis=1) > 1e-12) & (X.max(axis=1) < 1 - 1e-12)
+    resid = np.abs(np.bincount(conn.ravel(), weights=Rx[:, 0, :].ravel(), minlength=len(nodes))[interior]).max()
     h = 1.0 / (nn - 1)
     rate = ("%9.2f %9.2f" % (np.log(prev[0] / l2e) / np.log(prev[1] / h), np.log(prev[2] / mxe) / np.log(prev[1] / h))) if prev else "%9s %9s" % ("-", "-")
     print("%4d %9d %12.2e %12.3e %14.3e %14.3e %s %14.3e" % (nn, len(conn), np.abs(ue - up).max(), l2e, l2p, mxe, rate, resid))

experimental/adrh/ad_sine_n.py

@@ -1,4 +1,4 @@
-"""Numerics for ad_sine_p.py: P1 Galerkin on regular tetrahedra, direct solve, no stabilization."""
+"""Numerics for ad_sine_p.py: P1 Galerkin on regular triangles or tetrahedra, direct solve, no stabilization."""
 import os
 from proteus import *
 from proteus.default_n import *
@@ -16,8 +16,19 @@
 nLevels = 1
 subgridError = None
 shockCapturing = None
-# strong Dirichlet (the Dirichlet unknowns are eliminated), advective outflow from the interior value, prescribed diffusive flux
-numericalFluxType = NumericalFlux.StrongDirichletFactory({0: 'outFlow'}) if mixed else None
+# strong Dirichlet (the Dirichlet unknowns are eliminated) with advective outflow by the interior value and prescribed
+# fluxes, or weak Dirichlet by the interior-penalty exterior fluxes (see ad_sine_p.py for AD_BC)
+_WEAK = {'iipg': NumericalFlux.Advection_DiagonalUpwind_Diffusion_IIPG_exterior,
+         'nipg': NumericalFlux.Advection_DiagonalUpwind_Diffusion_NIPG_exterior,
+         'sipg': NumericalFlux.Advection_DiagonalUpwind_Diffusion_SIPG_exterior}
+if BC == 'strong':
+    numericalFluxType = None
+elif BC in _WEAK:
+    numericalFluxType = _WEAK[BC]
+elif BC == 'robin':
+    numericalFluxType = NumericalFlux.RobinFactory({0: getRobin}, base=NumericalFlux.StrongDirichletFactory({0: 'outFlow'}))
+else:
+    numericalFluxType = NumericalFlux.StrongDirichletFactory({0: 'outFlow'})
 multilevelNonlinearSolver = Newton
 levelNonlinearSolver = Newton
 fullNewtonFlag = True