From 62b809d01d6d8e047f730442a61e330d136c3851 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 23 Oct 2026 15:27:40 +0200 Subject: [PATCH 1/2] ADRH prototype: boundary kinds (prescribed inflow flux, weak Dirichlet IIPG/NIPG/SIPG) against Proteus boundary.h has one exterior-face kernel, exteriorFaceKernel, and a BC functor that returns a kind per component and quadrature point: Strong (nothing on the face), Flux (prescribed diffusive flux; advective flux prescribed where given, else the interior f(u).n on outflow and zero on inflow, as StrongDirichlet's 'outFlow' does) and WeakDirichlet (upwinded advective flux against f(g), -a grad phi.n + (C/h^p)(u - g) with h the face's longest edge, and the adjoint -sigma (u - g)(a^T n).grad w), each the engine's version of a Proteus exterior numerical-flux treatment. WeakDirichletParams carries C, p, sigma (2/1/0 IIPG, 10/1/-1 NIPG, 100/1/+1 SIPG, Proteus's class defaults). The header's comment maps every term to the Proteus routine it follows and names the two deliberate differences (the adjoint's tensor contraction, da/du in the Jacobian), neither visible with an isotropic constant a. ad_regression.py gets AD_BC=strong|mixed|inflow|iipg|nipg|sipg, with the Proteus side in ad_sine_p.py/ad_sine_n.py (StrongDirichletFactory with flux boundary conditions, or Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior). Engine and Proteus nodal solutions agree to 7e-15 or better in every mode at nn = 5, 9, 17; L2 rates 1.85-1.97. Boundary faces are picked plane by plane: a face whose nodes lie on different planes is interior. Co-Authored-By: Claude Fable 5.1 --- experimental/adrh/README.md | 21 ++- experimental/adrh/ad_regression.py | 60 +++++-- experimental/adrh/ad_sine_n.py | 13 +- experimental/adrh/ad_sine_p.py | 31 +++- experimental/adrh/boundary.h | 278 ++++++++++++++++++++++++----- experimental/adrh/module.cpp | 19 +- 6 files changed, 347 insertions(+), 75 deletions(-) diff --git a/experimental/adrh/README.md b/experimental/adrh/README.md index b4a9993797..ca6a485ac9 100644 --- a/experimental/adrh/README.md +++ b/experimental/adrh/README.md @@ -30,7 +30,22 @@ against 0.39 s for the Python path's kernels alone with its quadrature arrays al 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. + `ad_sine_n.py`) and the exact solution. Engine and Proteus nodal solutions agree to solver round-off (< 1e-14) in every + boundary mode; L2 error rates 1.85-1.97. -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). 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=`): + +| 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; advective flux prescribed where given, else the interior `f(u).n` on outflow and zero on inflow | `StrongDirichletFactory({ci: 'outFlow'})` with `diffusiveFluxBoundaryConditions` / `advectiveFluxBoundaryConditions` | `mixed`, `inflow` | +| `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). diff --git a/experimental/adrh/ad_regression.py b/experimental/adrh/ad_regression.py index eb06c9ec7a..a21edfefed 100644 --- a/experimental/adrh/ad_regression.py +++ b/experimental/adrh/ad_regression.py @@ -5,12 +5,20 @@ 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; + 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) + 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_BC=mode] python ad_regression.py [nn ...] default 5 9 17 """ import os import sys @@ -53,9 +61,23 @@ 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' +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') or BC in WEAK, BC +os.environ['AD_BC'] = BC quadf = Quadrature.SimplexGaussQuadrature(2, order) qpf = np.array([list(p)[:2] for p in quadf.points]); qwf = np.array(quadf.weights) +# 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)]} +for mode in WEAK: + FACES[mode] = [(3, ax, v) for ax in range(3) for v in (0.0, 1.0)] + + +def plane_faces(nodes, conn, axis, value): + """Element faces (element, local face = the opposite node) whose three nodes 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, 4) + return np.nonzero(on.sum(axis=1)[:, None] - on.astype(int) == 3) def engine_solution(nodes, conn): @@ -66,20 +88,31 @@ def engine_solution(nodes, conn): rows = np.repeat(conn, 4, axis=1).ravel(); cols = np.tile(conn, (1, 4)).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) + Rf, Jf = adrh_proto.advection_diffusion_sine_faces(b, kappa, omega, omega0, np.ascontiguousarray(fk, dtype='i'), + np.ascontiguousarray(kind, dtype='i'), + np.ascontiguousarray(nodes[face_conn]), np.zeros((len(fe), 1, 4)), + qpf, qwf, pc, pp, sig) rows_f = np.repeat(face_conn, 4, axis=1).ravel(); cols_f = np.tile(face_conn, (1, 4)).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 + lo, hi = nodes.min(axis=1) < 1e-12, nodes.max(axis=1) > 1 - 1e-12 + if BC in WEAK: + onb = np.zeros(nN, bool) + elif BC == 'mixed': + onb = lo | (nodes[:, 1:].max(axis=1) > 1 - 1e-12) + elif BC == 'inflow': + onb = (nodes[:, 1:].min(axis=1) < 1e-12) | (nodes[:, 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] @@ -101,6 +134,7 @@ def errors(u, nodes, conn, xn): sizes = [int(a) for a in sys.argv[1:]] or [5, 9, 17] prev = None +print("AD_BC=%s" % 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) diff --git a/experimental/adrh/ad_sine_n.py b/experimental/adrh/ad_sine_n.py index 2640f7de39..78616e90db 100644 --- a/experimental/adrh/ad_sine_n.py +++ b/experimental/adrh/ad_sine_n.py @@ -16,8 +16,17 @@ 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] +else: + numericalFluxType = NumericalFlux.StrongDirichletFactory({0: 'outFlow'}) multilevelNonlinearSolver = Newton levelNonlinearSolver = Newton fullNewtonFlag = True diff --git a/experimental/adrh/ad_sine_p.py b/experimental/adrh/ad_sine_p.py index 66fb92b49c..871d4c2880 100644 --- a/experimental/adrh/ad_sine_p.py +++ b/experimental/adrh/ad_sine_p.py @@ -54,26 +54,43 @@ def evaluate(self, t, c): eps = 1.0e-8 -# AD_MIXED=1: Dirichlet data on every face except x = 1, where the diffusive flux -kappa grad u . n is prescribed -# (advection leaves by the interior value); the default is Dirichlet data on the whole boundary -mixed = os.environ.get('AD_MIXED', '0') == '1' +# AD_BC selects the boundary treatment, matching ad_regression.py: +# strong (default) Dirichlet data on the whole boundary, imposed strongly (numericalFluxType = None) +# mixed as strong except the face x = 1, where the diffusive flux -kappa grad u . n is prescribed and +# advection leaves by the interior value (StrongDirichletFactory({0: 'outFlow'}); AD_MIXED=1 is this) +# inflow as mixed, and on the face x = 0 (inflow, b . n < 0) the advective flux (b . n) u and the +# diffusive flux are both prescribed +# iipg, nipg, sipg Dirichlet data on the whole boundary, imposed weakly by +# Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior (the Dirichlet dofs stay unknowns) +BC = os.environ.get('AD_BC', 'mixed' if os.environ.get('AD_MIXED', '0') == '1' else 'strong') +assert BC in ('strong', 'mixed', 'inflow', 'iipg', 'nipg', 'sipg'), BC +mixed = BC in ('mixed', 'inflow') +not_dirichlet = {'mixed': [3], 'inflow': [0, 3]}.get(BC, []) # indices into distances below: 0 = x=0, 3 = x=1 def getDBC(x, flag): distances = [x[i] - x0[i] for i in range(3)] + [x0[i] + L[i] - x[i] for i in range(3)] - if mixed: - distances = distances[:3] + distances[4:] # not the face x = 1 (index 3) + distances = [d for i, d in enumerate(distances) if i not in not_dirichlet] if min(distances) < eps: return lambda x, t: exact(x) def getDiffusiveFlux(x, flag): + # -kappa grad u . n with n = +-e_x on the faces x = 1 and x = 0 if mixed and abs(x[0] - (x0[0] + L[0])) < eps: return lambda x, t: -kappa * omega[0] * np.cos(np.dot(omega, x[:3]) + omega0) + if BC == 'inflow' and abs(x[0] - x0[0]) < eps: + return lambda x, t: kappa * omega[0] * np.cos(np.dot(omega, x[:3]) + omega0) + + +def getAdvectiveFlux(x, flag): + # (b . n) u on the inflow face x = 0, n = -e_x + if BC == 'inflow' and abs(x[0] - x0[0]) < eps: + return lambda x, t: -b[0] * exact(x) dirichletConditions = {0: getDBC} -advectiveFluxBoundaryConditions = {} +advectiveFluxBoundaryConditions = {0: getAdvectiveFlux} if BC == 'inflow' else {} diffusiveFluxBoundaryConditions = {0: {0: getDiffusiveFlux}} if mixed else {0: {}} -fluxBoundaryConditions = {0: 'outFlow' if mixed else 'noFlow'} +fluxBoundaryConditions = {0: 'outFlow' if BC != 'strong' else 'noFlow'} analyticalSolution = {} diff --git a/experimental/adrh/boundary.h b/experimental/adrh/boundary.h index 97a0231f09..4d77d99d38 100644 --- a/experimental/adrh/boundary.h +++ b/experimental/adrh/boundary.h @@ -1,14 +1,40 @@ // Exterior-face terms of the weak form for P1 simplices (tetrahedra), as a separate kernel. // // The Galerkin form integrated by parts leaves, on every exterior face, -// + integral over the face of F . n w dS, F = f(u) - a grad(phi) (the total flux), -// which proteus adds to the element residual as (advectiveFlux + diffusiveFlux) w dS. A face carries -// one of the following (the BC functor says which, per face and quadrature point): -// strong Dirichlet: nothing here, the unknown is removed or its row replaced by the driver; -// flux face: the advective flux is the interior one, f(u) . n (outflow; the engine does not yet -// look at inflow), and the diffusive flux is prescribed, g = -a grad u . n, a function of x, t and n. -// Weak Dirichlet (SIPG/IIPG penalty and adjoint terms) and upwinding against a prescribed inflow value are -// the next steps; they are further BC functor kinds, not changes of this structure. +// + integral over the face of F . n w dS, F = f(u) - sum_k a^{ik} grad phi^k (the total flux), +// which proteus adds to the element residual as (advectiveFlux + diffusiveFlux) w dS +// (femIntegrals.c updateExteriorElementBoundaryFlux). What stands in for F . n at a boundary quadrature +// point is decided per component by the BC functor, one of the kinds below. The structure follows +// proteus's exterior numerical flux classes (NumericalFlux.py, numericalFlux.c) so the two can be +// compared term by term: +// +// Strong nothing here; the driver removes the unknown or replaces its row +// (StrongDirichletFactory, or numericalFluxType=None). +// Flux the total flux is given. Diffusive part: g = -a grad u . n, prescribed (a +// diffusiveFluxBoundaryCondition). Advective part: prescribed where a value is given (an +// advectiveFluxBoundaryCondition; no Jacobian, as proteus zeroes dadvectiveFlux_left there); +// where none is, the interior flux f(u) . n on outflow (df . n >= 0) and ZERO on inflow, +// which is what StrongDirichlet's 'outFlow' does (calculateExteriorNumericalAdvectiveFlux_NoBC +// drops the inflow flux). +// WeakDirichlet the value g_D is imposed weakly, as Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior +// do. Advective part: upwinded between f(u) . n and f(g_D) . n by the signs of df(u) . n and +// df(g_D) . n (calculateExteriorNumericalAdvectiveFlux: a shock speed decides the mixed case, a +// transonic rarefaction is averaged), with a Jacobian only where the interior flux is taken. +// Diffusive part: -a grad phi . n + (C / h^p) (u - g_D), h the face diameter (its longest +// edge); C, p = 2, 1 for IIPG, 10, 1 for NIPG, 100, 1 for SIPG. Adjoint part (sigma = 0 IIPG, +// -1 NIPG, +1 SIPG): - sigma (u - g_D) (a^T n) . grad w dS, in the residual and the Jacobian. +// +// nc > 1: the kind is per component i. The diffusive part of component i sums over the potentials k with +// a^{ik} nonzero; penalty and adjoint are on the jump of u_k against g_D[k], and zero where component k is +// not itself WeakDirichlet (proteus: bc_u is u where no Dirichlet value is set). The upwind decision uses +// df^{ii}, as proteus does; the Jacobian of the interior flux uses the whole row df^{ij}. +// +// Two places where the engine is not proteus, on purpose: (1) proteus's adjoint +// (femIntegrals.c updateExteriorElementBoundaryDiffusionAdjoint) contracts the tensor as +// sum_{I,J} n_I a_{JI} grad w_I, the column sums of a, which equals (a^T n) . grad w only for a diagonal a; +// the engine uses (a^T n) . grad w. (2) proteus's Jacobian of the adjoint ignores da/du; so does the engine's, +// but its diffusive-flux Jacobian carries the da term like the element kernel. The regressions have an +// isotropic, constant a, where none of this shows. // // Output is per face (not scattered): Rf (nF, nc, nD) and Jf (nF, nc, nD, nc, nD); the driver adds // them to the element rows of the face's element. @@ -19,21 +45,45 @@ namespace adrh { - // A flux boundary condition: g_i = prescribed diffusive flux of component i, -a grad u . n + enum BCKind : int { Strong = 0, Flux = 1, WeakDirichlet = 2 }; + + // What the BC functor fills in at one boundary quadrature point, per component. template - struct FluxValues + struct FaceBC + { + int kind[nc]; + double gD[nc]; // WeakDirichlet: the Dirichlet value + double gdiff[nc]; // Flux: the prescribed diffusive flux -a grad u . n (0 = no diffusive flux) + bool hasAdv[nc]; // Flux: whether the advective flux is prescribed too, + double fadv[nc]; // and its value f . n + void reset() + { + for (int i = 0; i < nc; ++i) + { + kind[i] = Strong; gD[i] = 0.0; gdiff[i] = 0.0; hasAdv[i] = false; fadv[i] = 0.0; + } + } + }; + + // Penalty and adjoint parameters of the weak Dirichlet kinds (proteus's defaults per class). + struct WeakDirichletParams { - double g[nc]; - void zero() { std::memset(g, 0, sizeof(g)); } + double penalty_constant, penalty_power, sigma; + static WeakDirichletParams IIPG() { return {2.0, 1.0, 0.0}; } + static WeakDirichletParams NIPG() { return {10.0, 1.0, -1.0}; } + static WeakDirichletParams SIPG() { return {100.0, 1.0, 1.0}; } }; + // The boundary kernel. bc.evaluate(t, x, n, F, FaceBC&) is called per face F and quadrature point; + // the face index lets a functor look a per-face kind up. template - void exteriorFluxFaceKernel(const Coefficients &C, const BC &bc, int nF, int nQf, - const int *faceLocal, // (nF) local face index k = the opposite node - const double *xn, // (nF, 4, 3) coordinates of the face's element - const double *u, // (nF, nc, 4) the element's dofs - const double *qpf, const double *qwf, // face quadrature on the reference triangle - double t, double *Rf, double *Jf) + void exteriorFaceKernel(const Coefficients &C, const BC &bc, const WeakDirichletParams &wd, + int nF, int nQf, + const int *faceLocal, // (nF) local face index k = the opposite node + const double *xn, // (nF, 4, 3) coordinates of the face's element + const double *u, // (nF, nc, 4) the element's dofs + const double *qpf, const double *qwf, // face quadrature on the reference triangle + double t, double *Rf, double *Jf) { using Co = Coefficients; constexpr int nc = Co::nc, nS = 3, nD = 4; @@ -59,8 +109,21 @@ namespace adrh double out = 0.0; // orient away from the opposite node for (int I = 0; I < 3; ++I) out += nrm[I] * (X[fn[0] * 3 + I] - X[k * 3 + I]); if (out < 0.0) for (int I = 0; I < 3; ++I) nrm[I] = -nrm[I]; - // gradients of u on the element are constant; the coefficient functor wants grad u (unused by the - // flux terms here but part of its interface), so compute them from the element's affine map + // face diameter = its longest edge (mesh.cpp computeGeometricInfo_tetrahedron), for the penalty + double h = 0.0; + for (int p = 0; p < 3; ++p) + for (int r = p + 1; r < 3; ++r) + { + double d2 = 0.0; + for (int I = 0; I < 3; ++I) + { + const double d = X[fn[p] * 3 + I] - X[fn[r] * 3 + I]; + d2 += d * d; + } + h = std::fmax(h, std::sqrt(d2)); + } + const double penalty = wd.penalty_constant / std::pow(h, wd.penalty_power); + // gradients of the basis and of u on the element (constant for P1) double gw[nD][nS], adet; affineP1(X, gw, adet); const double *ue = u + F * nc * nD; @@ -73,7 +136,7 @@ namespace adrh } double Rloc[nc][nD] = {}; double Jloc[nc][nD][nc][nD] = {}; - Terms T; + Terms T, Tb; for (int q = 0; q < nQf; ++q) { double w[nD] = {0.0, 0.0, 0.0, 0.0}; @@ -89,28 +152,144 @@ namespace adrh uq[j] = 0.0; for (int b = 0; b < nD; ++b) uq[j] += ue[j * nD + b] * w[b]; } + FaceBC B; + B.reset(); + bc.evaluate(t, x, nrm, F, B); + bool anyWeak = false; + for (int i = 0; i < nc; ++i) anyWeak = anyWeak || (B.kind[i] == WeakDirichlet); zeroMasked(T); C.evaluate(t, x, uq, gu, T); - FluxValues g; - g.zero(); - bc.diffusiveFlux(t, x, nrm, g.g); + // the boundary state: g_D where a component is WeakDirichlet, the interior value elsewhere + double ub[nc]; + for (int j = 0; j < nc; ++j) ub[j] = (B.kind[j] == WeakDirichlet) ? B.gD[j] : uq[j]; + if (anyWeak) + { + zeroMasked(Tb); + C.evaluate(t, x, ub, gu, Tb); + } + // grad phi^k = sum_j dphi[k][j] grad u_j + double gphi[nc][nS]; + for (int kk = 0; kk < nc; ++kk) + for (int I = 0; I < nS; ++I) + { + gphi[kk][I] = 0.0; + for (int j = 0; j < nc; ++j) gphi[kk][I] += T.dphi[kk][j] * gu[j][I]; + } const double dS = qwf[q] * len; constexpr_for<0, nc>([&](auto i_) { constexpr int i = decltype(i_)::value; - double Fn = g.g[i]; + if (B.kind[i] == Strong) return; + double Fn = 0.0; // the flux . n at this point + double dFn[nc] = {}; // d Fn / d u_j, the part that multiplies w_b (no gradients) + // advective part if constexpr (Co::has_advection && Co::f_nz(i)) - for (int I = 0; I < 3; ++I) Fn += T.f[i][I] * nrm[I]; - for (int a = 0; a < nD; ++a) Rloc[i][a] += dS * Fn * w[a]; - constexpr_for<0, nc>([&](auto j_) { - constexpr int j = decltype(j_)::value; - if constexpr (Co::has_advection && Co::df_nz(i, j)) + { + double fl = 0.0, dl = 0.0; + for (int I = 0; I < 3; ++I) fl += T.f[i][I] * nrm[I]; + if constexpr (Co::df_nz(i, i)) + for (int I = 0; I < 3; ++I) dl += T.df[i][i][I] * nrm[I]; + bool interior = false, prescribed = false; + if (B.kind[i] == Flux) { - double dfn = 0.0; - for (int I = 0; I < 3; ++I) dfn += T.df[i][j][I] * nrm[I]; - for (int a = 0; a < nD; ++a) - for (int b = 0; b < nD; ++b) Jloc[i][a][j][b] += dS * dfn * w[b] * w[a]; + if (B.hasAdv[i]) { prescribed = true; Fn += B.fadv[i]; } + else interior = (dl >= 0.0); // inflow without a value: dropped, as proteus does } - }); + else // WeakDirichlet: upwind against f(g_D) + { + double fr = 0.0, dr = 0.0; + for (int I = 0; I < 3; ++I) fr += Tb.f[i][I] * nrm[I]; + if constexpr (Co::df_nz(i, i)) + for (int I = 0; I < 3; ++I) dr += Tb.df[i][i][I] * nrm[I]; + if (dl >= 0.0 && dr >= 0.0) interior = true; + else if (dl <= 0.0 && dr <= 0.0) interior = false; + else if (dl >= 0.0 && dr <= 0.0) + { + const double uj = ub[i] - uq[i], fj = fr - fl; + const double s = (std::fabs(uj) > std::fabs(fj) * 1.0e-16) ? fj / uj : 0.0; + interior = (s >= 0.0); + } + else + { // transonic rarefaction: proteus averages the two + Fn += 0.5 * (fl + fr); + prescribed = true; // and takes half the interior speed as the Jacobian + constexpr_for<0, nc>([&](auto j_) { + constexpr int j = decltype(j_)::value; + if constexpr (Co::df_nz(i, j)) + for (int I = 0; I < 3; ++I) dFn[j] += 0.5 * T.df[i][j][I] * nrm[I]; + }); + } + if (!prescribed && !interior) Fn += fr; + } + if (interior) + { + Fn += fl; + constexpr_for<0, nc>([&](auto j_) { + constexpr int j = decltype(j_)::value; + if constexpr (Co::df_nz(i, j)) + for (int I = 0; I < 3; ++I) dFn[j] += T.df[i][j][I] * nrm[I]; + }); + } + } + // diffusive part + if constexpr (Co::has_diffusion) + { + if (B.kind[i] == Flux) + Fn += B.gdiff[i]; + else + { + constexpr_for<0, nc>([&](auto k_) { + constexpr int kk = decltype(k_)::value; + if constexpr (Co::a_nz(i, kk)) + { + double an[nS]; // (a^{ik})^T n + for (int J = 0; J < nS; ++J) + { + an[J] = 0.0; + for (int I = 0; I < nS; ++I) an[J] += nrm[I] * T.a[i][kk][I][J]; + } + const double jump = uq[kk] - ub[kk]; // zero unless component kk is WeakDirichlet + double qn = 0.0; + for (int J = 0; J < nS; ++J) qn -= an[J] * gphi[kk][J]; + Fn += qn + penalty * jump; + // Jacobian of -n . a grad phi^k + penalty (u_k - g) and of the adjoint, which carry grad w_b + constexpr_for<0, nc>([&](auto j_) { + constexpr int j = decltype(j_)::value; + for (int b = 0; b < nD; ++b) + { + double d = 0.0; + if constexpr (Co::dphi_nz(kk, j)) + for (int J = 0; J < nS; ++J) d -= an[J] * T.dphi[kk][j] * gw[b][J]; + if constexpr (Co::da_nz(i, kk, j)) + for (int I = 0; I < nS; ++I) + for (int J = 0; J < nS; ++J) d -= nrm[I] * T.da[i][kk][j][I][J] * gphi[kk][J] * w[b]; + if (j == kk) d += penalty * w[b]; + for (int a = 0; a < nD; ++a) Jloc[i][a][j][b] += dS * d * w[a]; + } + }); + if (wd.sigma != 0.0) + { + for (int a = 0; a < nD; ++a) + { + double angw = 0.0; + for (int J = 0; J < nS; ++J) angw += an[J] * gw[a][J]; + Rloc[i][a] -= dS * wd.sigma * jump * angw; + if (B.kind[kk] == WeakDirichlet) + for (int b = 0; b < nD; ++b) Jloc[i][a][kk][b] -= dS * wd.sigma * w[b] * angw; + } + } + } + }); + } + } + for (int a = 0; a < nD; ++a) + { + Rloc[i][a] += dS * Fn * w[a]; + constexpr_for<0, nc>([&](auto j_) { + constexpr int j = decltype(j_)::value; + if (dFn[j] != 0.0) + for (int b = 0; b < nD; ++b) Jloc[i][a][j][b] += dS * dFn[j] * w[b] * w[a]; + }); + } }); } std::memcpy(Rf + F * nc * nD, Rloc, sizeof(Rloc)); @@ -118,19 +297,30 @@ namespace adrh } } - // g = -kappa grad(u_exact) . n for u_exact = sin(omega.x + omega0): the flux of the sine solution's diffusion - struct SineDiffusiveFlux + // Boundary data of the sine solution u = sin(omega.x + omega0) for the regression. The kind of each face + // comes from an array: 1 = Flux (diffusive flux given; advection leaves by the interior value, inflow dropped), + // 2 = Flux with the advective flux (b.n) u given too, 3 = WeakDirichlet with u given; anything else = Strong. + struct SineBC { - double kappa, omega[3], omega0; - SineDiffusiveFlux(double kappa_, const double *omega_, double omega0_) : kappa(kappa_), omega0(omega0_) + double b[3], kappa, omega[3], omega0; + const int *faceKind; + SineBC(const double *b_, double kappa_, const double *omega_, double omega0_, const int *faceKind_) + : kappa(kappa_), omega0(omega0_), faceKind(faceKind_) { - for (int I = 0; I < 3; ++I) omega[I] = omega_[I]; + for (int I = 0; I < 3; ++I) { b[I] = b_[I]; omega[I] = omega_[I]; } } - void diffusiveFlux(double, const double x[3], const double n[3], double *g) const + void evaluate(double, const double x[3], const double n[3], int F, FaceBC<1> &B) const { - double arg = omega0, wn = 0.0; - for (int I = 0; I < 3; ++I) { arg += omega[I] * x[I]; wn += omega[I] * n[I]; } - g[0] = -kappa * wn * std::cos(arg); + double arg = omega0, wn = 0.0, bn = 0.0; + for (int I = 0; I < 3; ++I) { arg += omega[I] * x[I]; wn += omega[I] * n[I]; bn += b[I] * n[I]; } + const double ue = std::sin(arg), gradn = wn * std::cos(arg); + switch (faceKind[F]) + { + case 1: B.kind[0] = Flux; B.gdiff[0] = -kappa * gradn; break; + case 2: B.kind[0] = Flux; B.gdiff[0] = -kappa * gradn; B.hasAdv[0] = true; B.fadv[0] = bn * ue; break; + case 3: B.kind[0] = WeakDirichlet; B.gD[0] = ue; break; + default: B.kind[0] = Strong; + } } }; } // namespace adrh diff --git a/experimental/adrh/module.cpp b/experimental/adrh/module.cpp index 871c8c382b..569f14859a 100644 --- a/experimental/adrh/module.cpp +++ b/experimental/adrh/module.cpp @@ -50,17 +50,24 @@ py::tuple runMasked(const C &coef, Arr xn, Arr u, py::object uold, Arr qp, Arr q PYBIND11_MODULE(adrh_proto, m) { - // flux faces of the advection-diffusion sine problem: advective outflow f(u).n plus the prescribed diffusive flux - m.def("advection_diffusion_sine_flux_faces", [](Arr b, double kappa, Arr omega, double omega0, py::array_t faceLocal, - Arr xn, Arr u, Arr qpf, Arr qwf) { + // exterior faces of the advection-diffusion sine problem; faceKind (nF) per face: 1 = flux face (diffusive flux + // given, advection by the interior value on outflow), 2 = flux face with the advective flux given too, + // 3 = weak Dirichlet (penalty_constant / h^penalty_power, adjoint sigma: 0 IIPG, -1 NIPG, +1 SIPG) + m.def("advection_diffusion_sine_faces", [](Arr b, double kappa, Arr omega, double omega0, + py::array_t faceLocal, + py::array_t faceKind, + Arr xn, Arr u, Arr qpf, Arr qwf, + double penalty_constant, double penalty_power, double sigma) { adrh::AdvectionDiffusionSine c(b.data(), kappa, omega.data(), omega0); - adrh::SineDiffusiveFlux bc(kappa, omega.data(), omega0); + adrh::SineBC bc(b.data(), kappa, omega.data(), omega0, faceKind.data()); + adrh::WeakDirichletParams wd{penalty_constant, penalty_power, sigma}; const int nF = xn.shape(0), nQf = qwf.shape(0); py::array_t Rf({nF, 1, 4}), Jf({nF, 1, 4, 1, 4}); { py::gil_scoped_release release; - adrh::exteriorFluxFaceKernel( - c, bc, nF, nQf, faceLocal.data(), xn.data(), u.data(), qpf.data(), qwf.data(), 0.0, Rf.mutable_data(), Jf.mutable_data()); + adrh::exteriorFaceKernel( + c, bc, wd, nF, nQf, faceLocal.data(), xn.data(), u.data(), qpf.data(), qwf.data(), 0.0, + Rf.mutable_data(), Jf.mutable_data()); } return py::make_tuple(Rf, Jf); }); From 400144d216175cbc9b3079a8445eb3785d7ffff7 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Fri, 23 Oct 2026 23:01:48 +0200 Subject: [PATCH 2/2] ADRH prototype: boundary kernel for 2D segments, Robin fluxes; regression in 2D and 3D exteriorFaceKernel is now generic in Coefficients::nSpace (2 or 3): a face is the nS nodes other than the opposite one, the normal is the edge's perpendicular in 2D and the cross product in 3D, the penalty's h is the longest face edge (the edge itself in 2D), and the reference quadrature has nS - 1 coordinates. The Flux kind takes a diffusive flux affine in the trace, gdiff + dgdiff u, which is what RobinFactory's alpha u + beta q.n = g amounts to (q.n = g/beta - (alpha/beta) u, Jacobian -alpha/beta). AdvectionDiffusionSine and SineBC are templates on the dimension; the module's advection_diffusion_sine and advection_diffusion_sine_faces take nd first. ad_regression.py: AD_ND=2|3, and a seventh mode AD_BC=robin (face x = 1 carries alpha u + beta q.n = g with alpha, beta = 1, 2, against RobinFactory over StrongDirichletFactory({0: 'outFlow'})). Engine and Proteus nodal solutions agree in all seven modes: 3D (nn 5/9/17) to 7.1e-15, 2D (nn 9/17/33/65) to 4.6e-14; L2 rates 2.00 in 2D, 1.90-1.97 in 3D. Boundary faces are picked plane by plane with the plane count following nd. Co-Authored-By: Claude Fable 5.1 --- experimental/adrh/README.md | 15 ++-- experimental/adrh/ad_regression.py | 95 ++++++++++---------- experimental/adrh/ad_sine_n.py | 4 +- experimental/adrh/ad_sine_p.py | 49 +++++++---- experimental/adrh/boundary.h | 135 +++++++++++++++++------------ experimental/adrh/coefficients.h | 13 +-- experimental/adrh/module.cpp | 51 +++++++---- 7 files changed, 212 insertions(+), 150 deletions(-) diff --git a/experimental/adrh/README.md b/experimental/adrh/README.md index ca6a485ac9..973e6e4899 100644 --- a/experimental/adrh/README.md +++ b/experimental/adrh/README.md @@ -29,20 +29,21 @@ 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 solver round-off (< 1e-14) in every - boundary mode; L2 error rates 1.85-1.97. + (`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. ## Boundary conditions (`boundary.h`) -`exteriorFaceKernel` adds the exterior-face terms per face (residual and Jacobian, scattered by the driver). 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=`): +`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; advective flux prescribed where given, else the interior `f(u).n` on outflow and zero on inflow | `StrongDirichletFactory({ci: 'outFlow'})` with `diffusiveFluxBoundaryConditions` / `advectiveFluxBoundaryConditions` | `mixed`, `inflow` | +| `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 diff --git a/experimental/adrh/ad_regression.py b/experimental/adrh/ad_regression.py index a21edfefed..fb9fce27e5 100644 --- a/experimental/adrh/ad_regression.py +++ b/experimental/adrh/ad_regression.py @@ -1,9 +1,9 @@ """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: +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): @@ -12,13 +12,15 @@ (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: [AD_BC=mode] 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 @@ -35,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): @@ -61,31 +76,19 @@ def proteus_solution(nn): return np.array(lm.u[0].dof), np.array(lm.mesh.nodeArray), np.array(lm.mesh.elementNodesArray) -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') or BC in WEAK, BC -os.environ['AD_BC'] = BC -quadf = Quadrature.SimplexGaussQuadrature(2, order) -qpf = np.array([list(p)[:2] for p in quadf.points]); qwf = np.array(quadf.weights) -# 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)]} -for mode in WEAK: - FACES[mode] = [(3, ax, v) for ax in range(3) for v in (0.0, 1.0)] - - def plane_faces(nodes, conn, axis, value): - """Element faces (element, local face = the opposite node) whose three nodes lie on the plane x[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, 4) - return np.nonzero(on.sum(axis=1)[:, None] - on.astype(int) == 3) + 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 FACES[BC]: @@ -96,21 +99,22 @@ def engine_solution(nodes, conn): 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_faces(b, kappa, omega, omega0, np.ascontiguousarray(fk, dtype='i'), + 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, 4)), - qpf, qwf, pc, pp, sig) - rows_f = np.repeat(face_conn, 4, axis=1).ravel(); cols_f = np.tile(face_conn, (1, 4)).ravel() + 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 - lo, hi = nodes.min(axis=1) < 1e-12, nodes.max(axis=1) > 1 - 1e-12 + 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 == 'mixed': - onb = lo | (nodes[:, 1:].max(axis=1) > 1 - 1e-12) + elif BC in ('mixed', 'robin'): + onb = lo | (X[:, 1:].max(axis=1) > 1 - 1e-12) elif BC == 'inflow': - onb = (nodes[:, 1:].min(axis=1) < 1e-12) | (nodes[:, 1:].max(axis=1) > 1 - 1e-12) + onb = (X[:, 1:].min(axis=1) < 1e-12) | (X[:, 1:].max(axis=1) > 1 - 1e-12) else: onb = lo | hi ub = exact(nodes) @@ -126,25 +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_BC=%s" % BC) +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)) diff --git a/experimental/adrh/ad_sine_n.py b/experimental/adrh/ad_sine_n.py index 78616e90db..e618c7db63 100644 --- a/experimental/adrh/ad_sine_n.py +++ b/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 * @@ -25,6 +25,8 @@ 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 diff --git a/experimental/adrh/ad_sine_p.py b/experimental/adrh/ad_sine_p.py index 871d4c2880..8ccaa9376d 100644 --- a/experimental/adrh/ad_sine_p.py +++ b/experimental/adrh/ad_sine_p.py @@ -1,9 +1,10 @@ """Steady advection-diffusion with the manufactured sine solution, as a plain Proteus model: the Python-path reference for the ADRH engine regression (ad_regression.py). - div(b u - kappa grad u) + d(x) = 0 on the unit cube, u = sin(omega.x + omega0) on all of the boundary. + div(b u - kappa grad u) + d(x) = 0 on the unit square or cube (AD_ND = 2 or 3), + u = sin(omega.x + omega0) on the boundary, imposed as AD_BC says. -Galerkin: no stabilization, strong Dirichlet conditions, so it discretizes exactly what the engine does. +Galerkin: no stabilization, so it discretizes exactly what the engine does. Parameters come from the environment (AD_*) so the driver can vary them. """ import os @@ -12,19 +13,19 @@ from proteus.default_p import * from proteus.TransportCoefficients import TC_base -nd = 3 -x0 = (0.0, 0.0, 0.0) -L = (1.0, 1.0, 1.0) +nd = int(os.environ.get('AD_ND', '3')) +x0 = (0.0,) * nd +L = (1.0,) * nd domain = None initialConditions = None -b = np.array([float(v) for v in os.environ.get('AD_B', '1.0,0.5,0.25').split(',')]) +b = np.array([float(v) for v in os.environ.get('AD_B', '1.0,0.5,0.25').split(',')])[:nd] kappa = float(os.environ.get('AD_KAPPA', '0.2')) -omega = np.array([float(v) for v in os.environ.get('AD_OMEGA', '3.0,2.0,4.0').split(',')]) +omega = np.array([float(v) for v in os.environ.get('AD_OMEGA', '3.0,2.0,4.0').split(',')])[:nd] omega0 = float(os.environ.get('AD_OMEGA0', '0.3')) def exact(x): - return np.sin(np.dot(omega, x[:3]) + omega0) + return np.sin(np.dot(omega, x[:nd]) + omega0) class AdvectionDiffusionSine(TC_base): @@ -44,7 +45,7 @@ def evaluate(self, t, c): c[('da', 0, 0, 0)][...] = 0.0 c[('phi', 0)][...] = u c[('dphi', 0, 0)][...] = 1.0 - arg = x @ omega + omega0 + arg = x[..., :nd] @ omega + omega0 c[('r', 0)][...] = -np.dot(b, omega) * np.cos(arg) - kappa * np.dot(omega, omega) * np.sin(arg) c[('dr', 0, 0)][...] = 0.0 @@ -60,27 +61,33 @@ def evaluate(self, t, c): # advection leaves by the interior value (StrongDirichletFactory({0: 'outFlow'}); AD_MIXED=1 is this) # inflow as mixed, and on the face x = 0 (inflow, b . n < 0) the advective flux (b . n) u and the # diffusive flux are both prescribed +# robin as mixed, but the face x = 1 carries alpha u + beta q.n = g (RobinFactory over the same base) # iipg, nipg, sipg Dirichlet data on the whole boundary, imposed weakly by # Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior (the Dirichlet dofs stay unknowns) BC = os.environ.get('AD_BC', 'mixed' if os.environ.get('AD_MIXED', '0') == '1' else 'strong') -assert BC in ('strong', 'mixed', 'inflow', 'iipg', 'nipg', 'sipg'), BC -mixed = BC in ('mixed', 'inflow') -not_dirichlet = {'mixed': [3], 'inflow': [0, 3]}.get(BC, []) # indices into distances below: 0 = x=0, 3 = x=1 +assert BC in ('strong', 'mixed', 'inflow', 'robin', 'iipg', 'nipg', 'sipg'), BC +robin_alpha, robin_beta = 1.0, 2.0 +# indices into the distance list of getDBC: 0 = the plane x = 0, nd = the plane x = 1 +not_dirichlet = {'mixed': [nd], 'robin': [nd], 'inflow': [0, nd]}.get(BC, []) def getDBC(x, flag): - distances = [x[i] - x0[i] for i in range(3)] + [x0[i] + L[i] - x[i] for i in range(3)] + distances = [x[i] - x0[i] for i in range(nd)] + [x0[i] + L[i] - x[i] for i in range(nd)] distances = [d for i, d in enumerate(distances) if i not in not_dirichlet] if min(distances) < eps: return lambda x, t: exact(x) +def _qn(x, sign): + """-kappa grad u . n on a face x = const with n = sign * e_x.""" + return -kappa * sign * omega[0] * np.cos(np.dot(omega, x[:nd]) + omega0) + + def getDiffusiveFlux(x, flag): - # -kappa grad u . n with n = +-e_x on the faces x = 1 and x = 0 - if mixed and abs(x[0] - (x0[0] + L[0])) < eps: - return lambda x, t: -kappa * omega[0] * np.cos(np.dot(omega, x[:3]) + omega0) + if BC in ('mixed', 'inflow') and abs(x[0] - (x0[0] + L[0])) < eps: + return lambda x, t: _qn(x, 1.0) if BC == 'inflow' and abs(x[0] - x0[0]) < eps: - return lambda x, t: kappa * omega[0] * np.cos(np.dot(omega, x[:3]) + omega0) + return lambda x, t: _qn(x, -1.0) def getAdvectiveFlux(x, flag): @@ -89,8 +96,14 @@ def getAdvectiveFlux(x, flag): return lambda x, t: -b[0] * exact(x) +def getRobin(x, flag): + # alpha u + beta q.n = g on the face x = 1, with g from the exact solution + if BC == 'robin' and abs(x[0] - (x0[0] + L[0])) < eps: + return lambda x, t: (robin_alpha, robin_beta, robin_alpha * exact(x) + robin_beta * _qn(x, 1.0)) + + dirichletConditions = {0: getDBC} advectiveFluxBoundaryConditions = {0: getAdvectiveFlux} if BC == 'inflow' else {} -diffusiveFluxBoundaryConditions = {0: {0: getDiffusiveFlux}} if mixed else {0: {}} +diffusiveFluxBoundaryConditions = {0: {0: getDiffusiveFlux}} if BC in ('mixed', 'inflow') else {0: {}} fluxBoundaryConditions = {0: 'outFlow' if BC != 'strong' else 'noFlow'} analyticalSolution = {} diff --git a/experimental/adrh/boundary.h b/experimental/adrh/boundary.h index 4d77d99d38..81d8faff36 100644 --- a/experimental/adrh/boundary.h +++ b/experimental/adrh/boundary.h @@ -1,4 +1,5 @@ -// Exterior-face terms of the weak form for P1 simplices (tetrahedra), as a separate kernel. +// Exterior-face terms of the weak form for P1 simplices (segments of triangles, triangles of tetrahedra), +// as a separate kernel. // // The Galerkin form integrated by parts leaves, on every exterior face, // + integral over the face of F . n w dS, F = f(u) - sum_k a^{ik} grad phi^k (the total flux), @@ -11,18 +12,19 @@ // Strong nothing here; the driver removes the unknown or replaces its row // (StrongDirichletFactory, or numericalFluxType=None). // Flux the total flux is given. Diffusive part: g = -a grad u . n, prescribed (a -// diffusiveFluxBoundaryCondition). Advective part: prescribed where a value is given (an -// advectiveFluxBoundaryCondition; no Jacobian, as proteus zeroes dadvectiveFlux_left there); -// where none is, the interior flux f(u) . n on outflow (df . n >= 0) and ZERO on inflow, -// which is what StrongDirichlet's 'outFlow' does (calculateExteriorNumericalAdvectiveFlux_NoBC -// drops the inflow flux). +// diffusiveFluxBoundaryCondition), optionally affine in the trace, g0 + dg u (RobinFactory: +// alpha u + beta q.n = g gives q.n = g/beta - (alpha/beta) u, with the Jacobian -alpha/beta). +// Advective part: prescribed where a value is given (an advectiveFluxBoundaryCondition; no +// Jacobian, as proteus zeroes dadvectiveFlux_left there); where none is, the interior flux +// f(u) . n on outflow (df . n >= 0) and ZERO on inflow, which is what StrongDirichlet's +// 'outFlow' does (calculateExteriorNumericalAdvectiveFlux_NoBC drops the inflow flux). // WeakDirichlet the value g_D is imposed weakly, as Advection_DiagonalUpwind_Diffusion_{IIPG,NIPG,SIPG}_exterior // do. Advective part: upwinded between f(u) . n and f(g_D) . n by the signs of df(u) . n and // df(g_D) . n (calculateExteriorNumericalAdvectiveFlux: a shock speed decides the mixed case, a // transonic rarefaction is averaged), with a Jacobian only where the interior flux is taken. // Diffusive part: -a grad phi . n + (C / h^p) (u - g_D), h the face diameter (its longest -// edge); C, p = 2, 1 for IIPG, 10, 1 for NIPG, 100, 1 for SIPG. Adjoint part (sigma = 0 IIPG, -// -1 NIPG, +1 SIPG): - sigma (u - g_D) (a^T n) . grad w dS, in the residual and the Jacobian. +// edge; the edge itself in 2D); C, p = 2, 1 for IIPG, 10, 1 for NIPG, 100, 1 for SIPG. Adjoint part +// (sigma = 0 IIPG, -1 NIPG, +1 SIPG): - sigma (u - g_D) (a^T n) . grad w dS, residual and Jacobian. // // nc > 1: the kind is per component i. The diffusive part of component i sums over the potentials k with // a^{ik} nonzero; penalty and adjoint are on the jump of u_k against g_D[k], and zero where component k is @@ -36,8 +38,8 @@ // but its diffusive-flux Jacobian carries the da term like the element kernel. The regressions have an // isotropic, constant a, where none of this shows. // -// Output is per face (not scattered): Rf (nF, nc, nD) and Jf (nF, nc, nD, nc, nD); the driver adds -// them to the element rows of the face's element. +// Output is per face (not scattered): Rf (nF, nc, nD) and Jf (nF, nc, nD, nc, nD), nD = nS + 1; the driver +// adds them to the element rows of the face's element. #pragma once #include #include @@ -53,14 +55,15 @@ namespace adrh { int kind[nc]; double gD[nc]; // WeakDirichlet: the Dirichlet value - double gdiff[nc]; // Flux: the prescribed diffusive flux -a grad u . n (0 = no diffusive flux) + double gdiff[nc]; // Flux: the prescribed diffusive flux -a grad u . n at u = 0 (0 = no diffusive flux), + double dgdiff[nc]; // plus dgdiff * u_i (Robin); 0 for a flux that does not depend on the trace bool hasAdv[nc]; // Flux: whether the advective flux is prescribed too, double fadv[nc]; // and its value f . n void reset() { for (int i = 0; i < nc; ++i) { - kind[i] = Strong; gD[i] = 0.0; gdiff[i] = 0.0; hasAdv[i] = false; fadv[i] = 0.0; + kind[i] = Strong; gD[i] = 0.0; gdiff[i] = 0.0; dgdiff[i] = 0.0; hasAdv[i] = false; fadv[i] = 0.0; } } }; @@ -75,44 +78,55 @@ namespace adrh }; // The boundary kernel. bc.evaluate(t, x, n, F, FaceBC&) is called per face F and quadrature point; - // the face index lets a functor look a per-face kind up. + // the face index lets a functor look a per-face kind up. nS = Coefficients::nSpace = 2 or 3. template void exteriorFaceKernel(const Coefficients &C, const BC &bc, const WeakDirichletParams &wd, int nF, int nQf, const int *faceLocal, // (nF) local face index k = the opposite node - const double *xn, // (nF, 4, 3) coordinates of the face's element - const double *u, // (nF, nc, 4) the element's dofs - const double *qpf, const double *qwf, // face quadrature on the reference triangle + const double *xn, // (nF, nS+1, 3) coordinates of the face's element + const double *u, // (nF, nc, nS+1) the element's dofs + const double *qpf, const double *qwf, // (nQf, nS-1), (nQf): quadrature on the reference face double t, double *Rf, double *Jf) { using Co = Coefficients; - constexpr int nc = Co::nc, nS = 3, nD = 4; - static_assert(Co::nSpace == 3, "tetrahedra"); + constexpr int nc = Co::nc, nS = Co::nSpace, nD = nS + 1, nFN = nS; // nFN = nodes of a face + static_assert(nS == 2 || nS == 3, "segments of triangles or triangles of tetrahedra"); for (int F = 0; F < nF; ++F) { const double *X = xn + F * nD * 3; const int k = faceLocal[F]; - int fn[3], m = 0; + int fn[nFN], m = 0; for (int b = 0; b < nD; ++b) if (b != k) fn[m++] = b; - double e1[3], e2[3], cr[3]; - for (int I = 0; I < 3; ++I) + // outward unit normal and the face measure factor len (|edge| in 2D, |e1 x e2| = 2 area in 3D) + double nrm[3] = {0.0, 0.0, 0.0}, len; + if constexpr (nS == 3) { - e1[I] = X[fn[1] * 3 + I] - X[fn[0] * 3 + I]; - e2[I] = X[fn[2] * 3 + I] - X[fn[0] * 3 + I]; + double e1[3], e2[3], cr[3]; + for (int I = 0; I < 3; ++I) + { + e1[I] = X[fn[1] * 3 + I] - X[fn[0] * 3 + I]; + e2[I] = X[fn[2] * 3 + I] - X[fn[0] * 3 + I]; + } + cr[0] = e1[1] * e2[2] - e1[2] * e2[1]; + cr[1] = e1[2] * e2[0] - e1[0] * e2[2]; + cr[2] = e1[0] * e2[1] - e1[1] * e2[0]; + len = std::sqrt(cr[0] * cr[0] + cr[1] * cr[1] + cr[2] * cr[2]); + for (int I = 0; I < 3; ++I) nrm[I] = cr[I] / len; + } + else + { + const double ex = X[fn[1] * 3 + 0] - X[fn[0] * 3 + 0], ey = X[fn[1] * 3 + 1] - X[fn[0] * 3 + 1]; + len = std::sqrt(ex * ex + ey * ey); + nrm[0] = ey / len; nrm[1] = -ex / len; } - cr[0] = e1[1] * e2[2] - e1[2] * e2[1]; - cr[1] = e1[2] * e2[0] - e1[0] * e2[2]; - cr[2] = e1[0] * e2[1] - e1[1] * e2[0]; - const double len = std::sqrt(cr[0] * cr[0] + cr[1] * cr[1] + cr[2] * cr[2]); - double nrm[3] = {cr[0] / len, cr[1] / len, cr[2] / len}; double out = 0.0; // orient away from the opposite node - for (int I = 0; I < 3; ++I) out += nrm[I] * (X[fn[0] * 3 + I] - X[k * 3 + I]); - if (out < 0.0) for (int I = 0; I < 3; ++I) nrm[I] = -nrm[I]; - // face diameter = its longest edge (mesh.cpp computeGeometricInfo_tetrahedron), for the penalty + for (int I = 0; I < nS; ++I) out += nrm[I] * (X[fn[0] * 3 + I] - X[k * 3 + I]); + if (out < 0.0) for (int I = 0; I < nS; ++I) nrm[I] = -nrm[I]; + // face diameter = its longest edge (mesh.cpp computeGeometricInfo_*), for the penalty double h = 0.0; - for (int p = 0; p < 3; ++p) - for (int r = p + 1; r < 3; ++r) + for (int p = 0; p < nFN; ++p) + for (int r = p + 1; r < nFN; ++r) { double d2 = 0.0; for (int I = 0; I < 3; ++I) @@ -139,10 +153,13 @@ namespace adrh Terms T, Tb; for (int q = 0; q < nQf; ++q) { - double w[nD] = {0.0, 0.0, 0.0, 0.0}; - w[fn[0]] = 1.0 - qpf[q * 2] - qpf[q * 2 + 1]; - w[fn[1]] = qpf[q * 2]; - w[fn[2]] = qpf[q * 2 + 1]; + double w[nD] = {}; + w[fn[0]] = 1.0; + for (int mm = 0; mm < nS - 1; ++mm) + { + w[fn[mm + 1]] = qpf[q * (nS - 1) + mm]; + w[fn[0]] -= qpf[q * (nS - 1) + mm]; + } double x[3] = {0.0, 0.0, 0.0}; for (int b = 0; b < nD; ++b) for (int I = 0; I < 3; ++I) x[I] += w[b] * X[b * 3 + I]; @@ -185,9 +202,9 @@ namespace adrh if constexpr (Co::has_advection && Co::f_nz(i)) { double fl = 0.0, dl = 0.0; - for (int I = 0; I < 3; ++I) fl += T.f[i][I] * nrm[I]; + for (int I = 0; I < nS; ++I) fl += T.f[i][I] * nrm[I]; if constexpr (Co::df_nz(i, i)) - for (int I = 0; I < 3; ++I) dl += T.df[i][i][I] * nrm[I]; + for (int I = 0; I < nS; ++I) dl += T.df[i][i][I] * nrm[I]; bool interior = false, prescribed = false; if (B.kind[i] == Flux) { @@ -197,9 +214,9 @@ namespace adrh else // WeakDirichlet: upwind against f(g_D) { double fr = 0.0, dr = 0.0; - for (int I = 0; I < 3; ++I) fr += Tb.f[i][I] * nrm[I]; + for (int I = 0; I < nS; ++I) fr += Tb.f[i][I] * nrm[I]; if constexpr (Co::df_nz(i, i)) - for (int I = 0; I < 3; ++I) dr += Tb.df[i][i][I] * nrm[I]; + for (int I = 0; I < nS; ++I) dr += Tb.df[i][i][I] * nrm[I]; if (dl >= 0.0 && dr >= 0.0) interior = true; else if (dl <= 0.0 && dr <= 0.0) interior = false; else if (dl >= 0.0 && dr <= 0.0) @@ -215,7 +232,7 @@ namespace adrh constexpr_for<0, nc>([&](auto j_) { constexpr int j = decltype(j_)::value; if constexpr (Co::df_nz(i, j)) - for (int I = 0; I < 3; ++I) dFn[j] += 0.5 * T.df[i][j][I] * nrm[I]; + for (int I = 0; I < nS; ++I) dFn[j] += 0.5 * T.df[i][j][I] * nrm[I]; }); } if (!prescribed && !interior) Fn += fr; @@ -226,7 +243,7 @@ namespace adrh constexpr_for<0, nc>([&](auto j_) { constexpr int j = decltype(j_)::value; if constexpr (Co::df_nz(i, j)) - for (int I = 0; I < 3; ++I) dFn[j] += T.df[i][j][I] * nrm[I]; + for (int I = 0; I < nS; ++I) dFn[j] += T.df[i][j][I] * nrm[I]; }); } } @@ -234,7 +251,10 @@ namespace adrh if constexpr (Co::has_diffusion) { if (B.kind[i] == Flux) - Fn += B.gdiff[i]; + { + Fn += B.gdiff[i] + B.dgdiff[i] * uq[i]; + dFn[i] += B.dgdiff[i]; + } else { constexpr_for<0, nc>([&](auto k_) { @@ -297,28 +317,33 @@ namespace adrh } } - // Boundary data of the sine solution u = sin(omega.x + omega0) for the regression. The kind of each face - // comes from an array: 1 = Flux (diffusive flux given; advection leaves by the interior value, inflow dropped), - // 2 = Flux with the advective flux (b.n) u given too, 3 = WeakDirichlet with u given; anything else = Strong. + // Boundary data of the sine solution u = sin(omega.x + omega0) for the regression, in nS dimensions. The kind + // of each face comes from an array: 1 = Flux (diffusive flux given; advection leaves by the interior value, + // inflow dropped), 2 = Flux with the advective flux (b.n) u given too, 3 = WeakDirichlet with u given, + // 4 = Robin alpha u + beta q.n = g with g from the exact solution; anything else = Strong. + template struct SineBC { - double b[3], kappa, omega[3], omega0; + double b[3], kappa, omega[3], omega0, alpha, beta; const int *faceKind; - SineBC(const double *b_, double kappa_, const double *omega_, double omega0_, const int *faceKind_) - : kappa(kappa_), omega0(omega0_), faceKind(faceKind_) + SineBC(const double *b_, double kappa_, const double *omega_, double omega0_, const int *faceKind_, + double alpha_ = 1.0, double beta_ = 1.0) + : kappa(kappa_), omega0(omega0_), alpha(alpha_), beta(beta_), faceKind(faceKind_) { - for (int I = 0; I < 3; ++I) { b[I] = b_[I]; omega[I] = omega_[I]; } + for (int I = 0; I < 3; ++I) { b[I] = I < nS ? b_[I] : 0.0; omega[I] = I < nS ? omega_[I] : 0.0; } } void evaluate(double, const double x[3], const double n[3], int F, FaceBC<1> &B) const { double arg = omega0, wn = 0.0, bn = 0.0; - for (int I = 0; I < 3; ++I) { arg += omega[I] * x[I]; wn += omega[I] * n[I]; bn += b[I] * n[I]; } - const double ue = std::sin(arg), gradn = wn * std::cos(arg); + for (int I = 0; I < nS; ++I) { arg += omega[I] * x[I]; wn += omega[I] * n[I]; bn += b[I] * n[I]; } + const double ue = std::sin(arg), qn = -kappa * wn * std::cos(arg); // exact u and -kappa grad u . n switch (faceKind[F]) { - case 1: B.kind[0] = Flux; B.gdiff[0] = -kappa * gradn; break; - case 2: B.kind[0] = Flux; B.gdiff[0] = -kappa * gradn; B.hasAdv[0] = true; B.fadv[0] = bn * ue; break; + case 1: B.kind[0] = Flux; B.gdiff[0] = qn; break; + case 2: B.kind[0] = Flux; B.gdiff[0] = qn; B.hasAdv[0] = true; B.fadv[0] = bn * ue; break; case 3: B.kind[0] = WeakDirichlet; B.gD[0] = ue; break; + case 4: // g = alpha u_exact + beta q_exact + B.kind[0] = Flux; B.gdiff[0] = (alpha * ue + beta * qn) / beta; B.dgdiff[0] = -alpha / beta; break; default: B.kind[0] = Strong; } } diff --git a/experimental/adrh/coefficients.h b/experimental/adrh/coefficients.h index 3fed15bb2e..606887bd59 100644 --- a/experimental/adrh/coefficients.h +++ b/experimental/adrh/coefficients.h @@ -111,9 +111,12 @@ namespace adrh // div(b u - kappa grad u) + d(x) = 0, u = sin(omega.x + omega0), // d(x) = -(b.omega) cos(omega.x + omega0) - kappa |omega|^2 sin(omega.x + omega0). // In proteus's form f = b u, a = kappa I, phi = u, r = d(x). The parameters are data of the functor. + // Steady advection-diffusion with the manufactured solution u = sin(omega.x + omega0) in nS = 2 or 3 dimensions: + // div(b u - kappa grad u) + r(x) = 0 with r chosen so u is exact (AnalyticalSolutions.LinearADR_Sine, c = 0). + template struct AdvectionDiffusionSine : DenseMasks { - static constexpr int nc = 1, nSpace = 3; + static constexpr int nc = 1, nSpace = nS_; static constexpr bool has_mass = false, has_hamiltonian = false; static constexpr bool da_nz(int, int, int) { return false; } // a does not depend on u static constexpr bool dr_nz(int, int) { return false; } // nor does the source @@ -121,15 +124,15 @@ namespace adrh AdvectionDiffusionSine(const double *b_, double kappa_, const double *omega_, double omega0_) : kappa(kappa_), omega0(omega0_) { - for (int I = 0; I < 3; ++I) { b[I] = b_[I]; omega[I] = omega_[I]; } + for (int I = 0; I < 3; ++I) { b[I] = I < nS_ ? b_[I] : 0.0; omega[I] = I < nS_ ? omega_[I] : 0.0; } } - void evaluate(double, const double x[3], const double u[1], const double gu[1][3], Terms<1, 3> &T) const + void evaluate(double, const double x[3], const double u[1], const double gu[1][nS_], Terms<1, nS_> &T) const { T.phi[0] = u[0]; T.dphi[0][0] = 1.0; double bw = 0.0, ww = 0.0, arg = omega0; - for (int I = 0; I < 3; ++I) { bw += b[I] * omega[I]; ww += omega[I] * omega[I]; arg += omega[I] * x[I]; } - for (int I = 0; I < 3; ++I) + for (int I = 0; I < nS_; ++I) { bw += b[I] * omega[I]; ww += omega[I] * omega[I]; arg += omega[I] * x[I]; } + for (int I = 0; I < nS_; ++I) { T.f[0][I] = b[I] * u[0]; T.df[0][0][I] = b[I]; diff --git a/experimental/adrh/module.cpp b/experimental/adrh/module.cpp index 569f14859a..bee9b23e98 100644 --- a/experimental/adrh/module.cpp +++ b/experimental/adrh/module.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "coefficients.h" #include "boundary.h" @@ -50,30 +51,42 @@ py::tuple runMasked(const C &coef, Arr xn, Arr u, py::object uold, Arr qp, Arr q PYBIND11_MODULE(adrh_proto, m) { - // exterior faces of the advection-diffusion sine problem; faceKind (nF) per face: 1 = flux face (diffusive flux - // given, advection by the interior value on outflow), 2 = flux face with the advective flux given too, - // 3 = weak Dirichlet (penalty_constant / h^penalty_power, adjoint sigma: 0 IIPG, -1 NIPG, +1 SIPG) - m.def("advection_diffusion_sine_faces", [](Arr b, double kappa, Arr omega, double omega0, + // exterior faces of the advection-diffusion sine problem in nd = 2 (segments) or 3 (triangles) dimensions; + // faceKind (nF) per face: 1 = flux face (diffusive flux given, advection by the interior value on outflow), + // 2 = flux face with the advective flux given too, 3 = weak Dirichlet (penalty_constant / h^penalty_power, + // adjoint sigma: 0 IIPG, -1 NIPG, +1 SIPG), 4 = Robin alpha u + beta q.n = g with (alpha, beta) given + m.def("advection_diffusion_sine_faces", [](int nd, Arr b, double kappa, Arr omega, double omega0, py::array_t faceLocal, py::array_t faceKind, Arr xn, Arr u, Arr qpf, Arr qwf, - double penalty_constant, double penalty_power, double sigma) { - adrh::AdvectionDiffusionSine c(b.data(), kappa, omega.data(), omega0); - adrh::SineBC bc(b.data(), kappa, omega.data(), omega0, faceKind.data()); - adrh::WeakDirichletParams wd{penalty_constant, penalty_power, sigma}; - const int nF = xn.shape(0), nQf = qwf.shape(0); - py::array_t Rf({nF, 1, 4}), Jf({nF, 1, 4, 1, 4}); + double penalty_constant, double penalty_power, double sigma, + double robin_alpha, double robin_beta) { + auto run = [&](auto dim) { + constexpr int nS = decltype(dim)::value; + adrh::AdvectionDiffusionSine c(b.data(), kappa, omega.data(), omega0); + adrh::SineBC bc(b.data(), kappa, omega.data(), omega0, faceKind.data(), robin_alpha, robin_beta); + adrh::WeakDirichletParams wd{penalty_constant, penalty_power, sigma}; + const int nF = xn.shape(0), nQf = qwf.shape(0); + py::array_t Rf({nF, 1, nS + 1}), Jf({nF, 1, nS + 1, 1, nS + 1}); + { + py::gil_scoped_release release; + adrh::exteriorFaceKernel, adrh::SineBC>( + c, bc, wd, nF, nQf, faceLocal.data(), xn.data(), u.data(), qpf.data(), qwf.data(), 0.0, + Rf.mutable_data(), Jf.mutable_data()); + } + return py::make_tuple(Rf, Jf); + }; + if (nd == 2) return run(std::integral_constant{}); + return run(std::integral_constant{}); + }); + m.def("advection_diffusion_sine", [](int nd, Arr b, double kappa, Arr omega, double omega0, Arr xn, Arr u, Arr qp, Arr qw) { + if (nd == 2) { - py::gil_scoped_release release; - adrh::exteriorFaceKernel( - c, bc, wd, nF, nQf, faceLocal.data(), xn.data(), u.data(), qpf.data(), qwf.data(), 0.0, - Rf.mutable_data(), Jf.mutable_data()); + adrh::AdvectionDiffusionSine<2> c(b.data(), kappa, omega.data(), omega0); + return runMasked, 2>(c, xn, u, py::none(), qp, qw, 0.0, 0.0); } - return py::make_tuple(Rf, Jf); - }); - m.def("advection_diffusion_sine", [](Arr b, double kappa, Arr omega, double omega0, Arr xn, Arr u, Arr qp, Arr qw) { - adrh::AdvectionDiffusionSine c(b.data(), kappa, omega.data(), omega0); - return runMasked(c, xn, u, py::none(), qp, qw, 0.0, 0.0); + adrh::AdvectionDiffusionSine<3> c(b.data(), kappa, omega.data(), omega0); + return runMasked, 3>(c, xn, u, py::none(), qp, qw, 0.0, 0.0); }); m.def("poisson3d_masked", [](Arr xn, Arr u, Arr qp, Arr qw) { return runMasked(adrh::Poisson3D(), xn, u, py::none(), qp, qw, 0.0, 0.0);