From d1cc6581a27245c3ca77458631f9d71b88a83389 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 13:19:51 -0700 Subject: [PATCH 01/15] Quadrature-point (delta) finite element as a PETSc plugin space Adds create_delta_fe(quad, polytope): a PetscFE whose basis is the identity on the mesh quadrature rule, built from a UW3-registered PetscSpace type "uwdelta" (uw_delta_space.h, PetscSpaceRegister) and a PETSCDUALSPACESIMPLE dual space with one point-evaluation functional per rule point. A field of this type in the auxiliary DM is read by the pointwise functions as a[] with no interpolation: the intended carrier for the semi-Lagrangian history (values injected at the integration points, never sampled elsewhere). All dofs sit on the cell, so the local vector is (ncells, Nq). PETSc's own PETSCSPACEPOINT is the same idea but errors unless tabulated at exactly its own points in its own order, which breaks PetscFESetUp (one point per functional), face tabulation in PetscDSSetUp and boundary integrals over auxiliary fields; the plugin type returns zeros off its rule instead and needs no PETSc patch, so stock conda PETSc works. tests/test_0064: identity tabulation on all four cell types, cell-only section layout, zeros off-rule and permuted identity as negative controls. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- setup.py | 8 + .../cython/petsc_quadrature_fe.pyx | 192 ++++++++++++++++++ src/underworld3/cython/uw_delta_space.h | 157 ++++++++++++++ tests/test_0064_quadrature_point_fe.py | 102 ++++++++++ 4 files changed, 459 insertions(+) create mode 100644 src/underworld3/cython/petsc_quadrature_fe.pyx create mode 100644 src/underworld3/cython/uw_delta_space.h create mode 100644 tests/test_0064_quadrature_point_fe.py diff --git a/setup.py b/setup.py index 8d46fcac3..924d753b3 100644 --- a/setup.py +++ b/setup.py @@ -194,6 +194,14 @@ def configure(): extra_compile_args=extra_compile_args, **conf, ), + Extension( + "underworld3.cython.petsc_quadrature_fe", + sources=[ + "src/underworld3/cython/petsc_quadrature_fe.pyx", + ], + extra_compile_args=extra_compile_args, + **conf, + ), Extension( "underworld3.cython.petsc_maths", sources=[ diff --git a/src/underworld3/cython/petsc_quadrature_fe.pyx b/src/underworld3/cython/petsc_quadrature_fe.pyx new file mode 100644 index 000000000..1497c3f56 --- /dev/null +++ b/src/underworld3/cython/petsc_quadrature_fe.pyx @@ -0,0 +1,192 @@ +# cython: language_level=3 +r""" +Quadrature-point finite element (the "delta space"). + +A ``PetscFE`` whose basis functions are Kronecker deltas at the points of a +quadrature rule and whose dual space is point evaluation at those same +points. Tabulated on its own rule the basis is the identity matrix, so a +field of this type that is read by the assembler as an auxiliary field +(``a[]`` in the pointwise functions) delivers the stored value at each +quadrature point with no interpolation at all. The dofs all sit on the cell +interior, so the local vector is laid out cell-major, point-minor. + +Use it for values that are *injected* at the integration points (a +semi-Lagrangian history, a per-point material property reconstructed from a +swarm). It cannot be *sampled* anywhere else: the derivative tabulation is +zero and evaluation at points off the rule returns zeros. + +Built from a UW3-registered prime space ``uwdelta`` (``uw_delta_space.h``, +a PETSc plugin type: PETSc's own ``PETSCSPACEPOINT`` cannot be tabulated +anywhere but its own points, which breaks ``PetscFESetUp``, face tabulation +and boundary integrals) and a ``PETSCDUALSPACESIMPLE`` dual space, through +``PetscFECreateFromSpaces``. Works on any PETSc build. petsc4py cannot +construct the one-point delta functionals itself (``Quad`` has no +``setData``), which is why this helper is Cython. + +Scalar (one-component) elements only. +""" + +from petsc4py import PETSc +from petsc4py.PETSc cimport FE, PetscFE, Quad, PetscQuadrature, DM, PetscDM +from petsc4py.PETSc cimport PetscSpace, PetscDualSpace, PetscObject, MPI_Comm +from petsc4py.PETSc cimport CHKERR as CHKERRQ +from underworld3.cython.petsc_types cimport PetscInt, PetscReal, PetscErrorCode + +import numpy as np + + +cdef extern from "petsc.h" nogil: + MPI_Comm PETSC_COMM_SELF + ctypedef int DMPolytopeType + + PetscErrorCode PetscSpaceCreate(MPI_Comm, PetscSpace*) + PetscErrorCode PetscSpaceSetType(PetscSpace, const char*) + PetscErrorCode PetscSpaceSetNumVariables(PetscSpace, PetscInt) + PetscErrorCode PetscSpaceSetNumComponents(PetscSpace, PetscInt) + PetscErrorCode PetscSpaceSetUp(PetscSpace) + + PetscErrorCode PetscDualSpaceCreate(MPI_Comm, PetscDualSpace*) + PetscErrorCode PetscDualSpaceSetType(PetscDualSpace, const char*) + PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace, PetscDM) + PetscErrorCode PetscDualSpaceSetNumComponents(PetscDualSpace, PetscInt) + PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace, PetscInt) + PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace, PetscInt, PetscQuadrature) + PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace) + + PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm, DMPolytopeType, PetscDM*) + PetscErrorCode DMDestroy(PetscDM*) + + PetscErrorCode PetscQuadratureCreate(MPI_Comm, PetscQuadrature*) + PetscErrorCode PetscQuadratureSetData(PetscQuadrature, PetscInt, PetscInt, PetscInt, const PetscReal*, const PetscReal*) + PetscErrorCode PetscQuadratureGetData(PetscQuadrature, PetscInt*, PetscInt*, PetscInt*, const PetscReal**, const PetscReal**) + PetscErrorCode PetscQuadratureDestroy(PetscQuadrature*) + + PetscErrorCode PetscFECreateFromSpaces(PetscSpace, PetscDualSpace, PetscQuadrature, PetscQuadrature, PetscFE*) + PetscErrorCode PetscObjectReference(PetscObject) + PetscErrorCode PetscObjectSetName(PetscObject, const char*) + PetscErrorCode PetscMalloc(size_t, void**) + + ctypedef struct _n_PetscTabulation: + PetscInt K + PetscInt Nr + PetscInt Np + PetscInt Nb + PetscInt Nc + PetscInt cdim + PetscReal **T + ctypedef _n_PetscTabulation* PetscTabulation + PetscErrorCode PetscFECreateTabulation(PetscFE, PetscInt, PetscInt, const PetscReal*, PetscInt, PetscTabulation*) + PetscErrorCode PetscTabulationDestroy(PetscTabulation*) + +cdef extern from "uw_delta_space.h" nogil: + PetscErrorCode UWDeltaSpaceRegister() + PetscErrorCode UWDeltaSpaceSetPoints(PetscSpace, PetscQuadrature) + + +# Register the plugin space type once, at import. +CHKERRQ(UWDeltaSpaceRegister()) + + +def create_delta_fe(Quad quad, int polytope, name="quadrature_point_fe"): + r"""Build the scalar quadrature-point element on ``quad``. + + Parameters + ---------- + quad : petsc4py.PETSc.Quad + The cell rule the element's points coincide with. Take it from an + existing field, ``fe.getQuadrature()``, so it is the mesh's rule. + polytope : int + The reference cell type (``dm.getCellType(cStart)``) for the dual + space's reference cell. + name : str + PETSc object name. + + Returns + ------- + petsc4py.PETSc.FE + Element of dimension ``Nq`` (points in the rule), one component, + with ``quad`` as its cell quadrature and no face quadrature. + """ + cdef PetscInt qdim = 0, qNc = 0, Nq = 0, i, d + cdef const PetscReal *points = NULL + cdef const PetscReal *weights = NULL + cdef PetscReal *fpts = NULL + cdef PetscReal *fwts = NULL + cdef PetscQuadrature functional = NULL + cdef PetscSpace P = NULL + cdef PetscDualSpace Q = NULL + cdef PetscDM refcell = NULL + cdef PetscFE cfe = NULL + cdef FE pyfe + + CHKERRQ(PetscQuadratureGetData(quad.quad, &qdim, &qNc, &Nq, &points, &weights)) + if qNc != 1: + raise ValueError("create_delta_fe: the rule must have one component") + + # Prime space: deltas at the rule's points. + CHKERRQ(PetscSpaceCreate(PETSC_COMM_SELF, &P)) + CHKERRQ(PetscSpaceSetType(P, b"uwdelta")) + CHKERRQ(PetscSpaceSetNumVariables(P, qdim)) + CHKERRQ(PetscSpaceSetNumComponents(P, 1)) + CHKERRQ(UWDeltaSpaceSetPoints(P, quad.quad)) + CHKERRQ(PetscSpaceSetUp(P)) + + # Dual space: one point-evaluation functional per rule point, all on the + # cell interior of the reference cell. + CHKERRQ(DMPlexCreateReferenceCell(PETSC_COMM_SELF, polytope, &refcell)) + CHKERRQ(PetscDualSpaceCreate(PETSC_COMM_SELF, &Q)) + CHKERRQ(PetscDualSpaceSetType(Q, b"simple")) + CHKERRQ(PetscDualSpaceSetDM(Q, refcell)) + CHKERRQ(PetscDualSpaceSetNumComponents(Q, 1)) + CHKERRQ(PetscDualSpaceSimpleSetDimension(Q, Nq)) + for i in range(Nq): + # PetscQuadratureSetData takes ownership: arrays must be PetscMalloc'd. + CHKERRQ(PetscMalloc(sizeof(PetscReal) * qdim, &fpts)) + CHKERRQ(PetscMalloc(sizeof(PetscReal), &fwts)) + for d in range(qdim): + fpts[d] = points[i * qdim + d] + fwts[0] = 1.0 + CHKERRQ(PetscQuadratureCreate(PETSC_COMM_SELF, &functional)) + CHKERRQ(PetscQuadratureSetData(functional, qdim, 1, 1, fpts, fwts)) + # SimpleSetFunctional duplicates; release ours. + CHKERRQ(PetscDualSpaceSimpleSetFunctional(Q, i, functional)) + CHKERRQ(PetscQuadratureDestroy(&functional)) + CHKERRQ(PetscDualSpaceSetUp(Q)) + CHKERRQ(DMDestroy(&refcell)) + + # PetscFECreateFromSpaces consumes P, Q and the quadrature: keep the + # caller's Quad alive by taking a reference first. No face quadrature. + CHKERRQ(PetscObjectReference(quad.quad)) + CHKERRQ(PetscFECreateFromSpaces(P, Q, quad.quad, NULL, &cfe)) + CHKERRQ(PetscObjectSetName(cfe, name.encode())) + + pyfe = FE() + pyfe.fe = cfe + return pyfe + + +def tabulate(FE fe, points, int K=0): + r"""Tabulate ``fe``'s basis at reference-cell ``points``. + + Returns the value tabulation as an array shaped ``(Np, Nb, Nc)``. + Exposed for tests: on its own rule the delta element returns the + identity. + """ + cdef PetscTabulation T = NULL + cdef PetscInt Np, Nb, Nc, p, b, c + pts = np.ascontiguousarray(points, dtype=np.float64) + if pts.ndim != 2: + raise ValueError("points must be (Np, dim)") + cdef double[:, ::1] pv = pts + Np = pts.shape[0] + CHKERRQ(PetscFECreateTabulation(fe.fe, 1, Np, &pv[0, 0], K, &T)) + Nb = T.Nb + Nc = T.Nc + out = np.empty((Np, Nb, Nc), dtype=np.float64) + cdef double[:, :, ::1] ov = out + for p in range(Np): + for b in range(Nb): + for c in range(Nc): + ov[p, b, c] = T.T[0][(p * Nb + b) * Nc + c] + CHKERRQ(PetscTabulationDestroy(&T)) + return out diff --git a/src/underworld3/cython/uw_delta_space.h b/src/underworld3/cython/uw_delta_space.h new file mode 100644 index 000000000..65214c164 --- /dev/null +++ b/src/underworld3/cython/uw_delta_space.h @@ -0,0 +1,157 @@ +/* + * UWDELTA: a PetscSpace of Kronecker deltas at the points of a quadrature rule. + * + * Registered as a PetscSpace type ("uwdelta") through PETSc's plugin API, so + * it works on any PETSc build, stock conda packages included. It is the + * prime space of the quadrature-point finite element: tabulated on its own + * rule the basis is the identity, tabulated anywhere else (a face rule, a + * different cell rule) it is zero, and derivatives are always zero. + * + * PETSc's own PETSCSPACEPOINT is the same idea but (as of 3.25) it errors + * unless asked for exactly its own points in its own order, which breaks + * PetscFESetUp (one point per functional), face tabulation in PetscDSSetUp + * and boundary integrals over auxiliary fields. That is why this type exists. + * + * Header-only: include from exactly one extension module. + */ +#ifndef UW_DELTA_SPACE_H +#define UW_DELTA_SPACE_H + +#include +#include + +#define UWDELTA_TOL 1.0e-10 + +typedef struct { + PetscQuadrature quad; /* the rule whose points carry the deltas */ +} UWDeltaSpace; + +static PetscErrorCode UWDeltaSpace_Destroy(PetscSpace sp) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + + PetscFunctionBegin; + PetscCall(PetscQuadratureDestroy(&dl->quad)); + PetscCall(PetscFree(dl)); + PetscFunctionReturn(PETSC_SUCCESS); +} + +static PetscErrorCode UWDeltaSpace_SetUp(PetscSpace sp) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + + PetscFunctionBegin; + PetscCheck(dl->quad, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "UWDELTA space has no points: call UWDeltaSpaceSetPoints() first"); + PetscFunctionReturn(PETSC_SUCCESS); +} + +static PetscErrorCode UWDeltaSpace_View(PetscSpace sp, PetscViewer viewer) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + PetscBool isascii; + PetscInt Nq = 0; + + PetscFunctionBegin; + PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); + if (isascii) { + if (dl->quad) PetscCall(PetscQuadratureGetData(dl->quad, NULL, NULL, &Nq, NULL, NULL)); + PetscCall(PetscViewerASCIIPrintf(viewer, "UWDELTA space in dimension %" PetscInt_FMT " on %" PetscInt_FMT " points\n", sp->Nv, Nq)); + } + PetscFunctionReturn(PETSC_SUCCESS); +} + +static PetscErrorCode UWDeltaSpace_GetDimension(PetscSpace sp, PetscInt *dim) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + PetscInt Nq = 0; + + PetscFunctionBegin; + if (dl->quad) PetscCall(PetscQuadratureGetData(dl->quad, NULL, NULL, &Nq, NULL, NULL)); + *dim = Nq; + PetscFunctionReturn(PETSC_SUCCESS); +} + +/* B is laid out [point][basis][component] as for every PetscSpace. Basis i is + the delta at the space's point i: a requested point coincident with point i + gives e_i, any other point gives zero. All components share the basis. */ +static PetscErrorCode UWDeltaSpace_Evaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[]) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + const PetscInt dim = sp->Nv, Nc = sp->Nc; + const PetscReal *qp; + PetscInt pdim = 0, p, i, d, c; + + PetscFunctionBegin; + PetscCheck(dl->quad, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "UWDELTA space has no points"); + PetscCall(PetscQuadratureGetData(dl->quad, NULL, NULL, &pdim, &qp, NULL)); + if (B) { + PetscCall(PetscArrayzero(B, npoints * pdim * Nc)); + for (p = 0; p < npoints; ++p) { + for (i = 0; i < pdim; ++i) { + for (d = 0; d < dim; ++d) { + if (PetscAbsReal(points[p * dim + d] - qp[i * dim + d]) > UWDELTA_TOL) break; + } + if (d >= dim) { + for (c = 0; c < Nc; ++c) B[(p * pdim + i) * Nc + c] = 1.0; + break; + } + } + } + } + if (D) PetscCall(PetscArrayzero(D, npoints * pdim * Nc * dim)); + if (H) PetscCall(PetscArrayzero(H, npoints * pdim * Nc * dim * dim)); + PetscFunctionReturn(PETSC_SUCCESS); +} + +static PetscErrorCode PetscSpaceCreate_UWDelta(PetscSpace sp) +{ + UWDeltaSpace *dl; + + PetscFunctionBegin; + PetscCall(PetscNew(&dl)); + dl->quad = NULL; + sp->data = dl; + sp->maxDegree = PETSC_INT_MAX; + + sp->ops->setfromoptions = NULL; + sp->ops->setup = UWDeltaSpace_SetUp; + sp->ops->view = UWDeltaSpace_View; + sp->ops->destroy = UWDeltaSpace_Destroy; + sp->ops->getdimension = UWDeltaSpace_GetDimension; + sp->ops->evaluate = UWDeltaSpace_Evaluate; + sp->ops->getheightsubspace = NULL; + PetscFunctionReturn(PETSC_SUCCESS); +} + +/* Idempotent: PETSc's function list rejects nothing on re-registration, but + registering once per process keeps the list clean. */ +static PetscErrorCode UWDeltaSpaceRegister(void) +{ + static PetscBool registered = PETSC_FALSE; + + PetscFunctionBegin; + if (!registered) { + PetscCall(PetscSpaceRegister("uwdelta", PetscSpaceCreate_UWDelta)); + registered = PETSC_TRUE; + } + PetscFunctionReturn(PETSC_SUCCESS); +} + +/* Set the rule whose points carry the deltas (duplicated; caller keeps its own). */ +static PetscErrorCode UWDeltaSpaceSetPoints(PetscSpace sp, PetscQuadrature q) +{ + UWDeltaSpace *dl = (UWDeltaSpace *)sp->data; + PetscBool isdelta; + PetscInt qdim; + + PetscFunctionBegin; + PetscCall(PetscObjectTypeCompare((PetscObject)sp, "uwdelta", &isdelta)); + PetscCheck(isdelta, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Space is not of type uwdelta"); + PetscCall(PetscQuadratureGetData(q, &qdim, NULL, NULL, NULL, NULL)); + PetscCheck(qdim == sp->Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_INCOMP, "Rule dimension %" PetscInt_FMT " != space variables %" PetscInt_FMT, qdim, sp->Nv); + PetscCall(PetscQuadratureDestroy(&dl->quad)); + PetscCall(PetscQuadratureDuplicate(q, &dl->quad)); + PetscFunctionReturn(PETSC_SUCCESS); +} + +#endif /* UW_DELTA_SPACE_H */ diff --git a/tests/test_0064_quadrature_point_fe.py b/tests/test_0064_quadrature_point_fe.py new file mode 100644 index 000000000..e989750a1 --- /dev/null +++ b/tests/test_0064_quadrature_point_fe.py @@ -0,0 +1,102 @@ +"""Quadrature-point ("delta") finite element. + +The element's basis is the identity on the mesh quadrature rule, its dofs +all live on the cell, and it tabulates to zero at any point off its rule. Those three properties are what +let a field of this type carry pre-evaluated values straight into the +pointwise functions as ``a[]``. +""" + +import numpy as np +import pytest +from petsc4py import PETSc + +import underworld3 as uw +from underworld3.cython.petsc_quadrature_fe import create_delta_fe, tabulate + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +CELLS = [(2, True), (3, True), (2, False), (3, False)] +IDS = ["triangle", "tetrahedron", "quadrilateral", "hexahedron"] + + +def _box(dim, simplex): + """A bare clone of a UW3 mesh DM (no fields) and its cell polytope.""" + if simplex: + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, cellSize=0.5, qdegree=2, + ) + else: + mesh = uw.meshing.StructuredQuadBox(elementRes=(2,) * dim, qdegree=2) + dm = mesh.dm.clone() + cStart, _ = dm.getHeightStratum(0) + return dm, dm.getCellType(cStart) + + +def _rule(dim, simplex, qdegree): + ref = PETSc.FE().createDefault(dim, 1, simplex, qdegree, "ref_", PETSc.COMM_SELF) + quad = ref.getQuadrature() + pts = np.array(quad.getData()[0]).reshape(-1, dim) + return ref, quad, pts + + +@pytest.mark.parametrize("dim,simplex", CELLS, ids=IDS) +@pytest.mark.parametrize("qdegree", [1, 2]) +def test_identity_on_own_rule(dim, simplex, qdegree): + _, quad, pts = _rule(dim, simplex, qdegree) + _, polytope = _box(dim, simplex) + fe = create_delta_fe(quad, polytope) + + assert fe.getDimension() == len(pts) + assert fe.getNumComponents() == 1 + B = tabulate(fe, pts)[:, :, 0] + assert np.array_equal(B, np.eye(len(pts))) + # Derivatives are zero by construction. + D = tabulate(fe, pts, K=1) + assert D.shape[0] == len(pts) + + +@pytest.mark.parametrize("dim,simplex", CELLS, ids=IDS) +def test_dofs_live_on_the_cell(dim, simplex): + """Local layout is (ncells, Nq): every dof on the cell, none elsewhere.""" + dm, polytope = _box(dim, simplex) + _, quad, pts = _rule(dim, simplex, 2) + fe = create_delta_fe(quad, polytope) + dm.setNumFields(1) + dm.setField(0, fe) + # PetscDSSetUp asks for the face tabulation; the delta space answers zeros. + dm.createDS() + section = dm.getLocalSection() + cStart, cEnd = dm.getHeightStratum(0) + for c in range(cStart, cEnd): + assert section.getDof(c) == len(pts) + pStart, pEnd = dm.getChart() + for p in range(pStart, pEnd): + if not (cStart <= p < cEnd): + assert section.getDof(p) == 0 + assert section.getStorageSize() == (cEnd - cStart) * len(pts) + + +def test_off_rule_points_are_zero(): + """Negative control: a point that is not on the rule contributes + nothing, and the same points in a different order give the permuted + identity (PETSc's own point space compares point p only with its own p).""" + dim, simplex = 2, True + _, quad, pts = _rule(dim, simplex, 2) + _, polytope = _box(dim, simplex) + fe = create_delta_fe(quad, polytope) + Nq = len(pts) + perm = np.arange(Nq)[::-1] + off = pts[:2] + 0.05 + + assert np.all(tabulate(fe, off) == 0.0) + Bperm = tabulate(fe, pts[perm])[:, :, 0] + assert np.array_equal(Bperm, np.eye(Nq)[perm]) + + +def test_rule_is_the_mesh_rule(): + """The element must be built on the rule every other field uses, and + that rule is fixed by the quadrature degree, not the field degree.""" + for degree in (1, 2, 3): + fe = PETSc.FE().createDefault(2, 1, True, 2, f"p{degree}_", PETSc.COMM_SELF) + pts = np.array(fe.getQuadrature().getData()[0]).reshape(-1, 2) + assert len(pts) == 6 From 178edbb5c47dce259e1da576b37d39e999f6fded Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 14:07:28 -0700 Subject: [PATCH 02/15] IntegrationPointVariable: a public field stored at the mesh integration points A peer of MeshVariable and swarm variables, built on the quadrature-point element (create_delta_fe) on mesh.integration_rule, so the pointwise functions read its values at the integration points with no interpolation. One value per rule point per cell; cell_data views (ncells, Nq, 1); coords are the assembler's own integration points (DMPlexComputeCellGeometryFEM). evaluate() is defined as the nearest integration point of the owning cell - the only extension under which a query agrees with what the assembler used at that point; exterior points take the nearest point on the rank. The JIT refuses derivatives of the symbol (the tabulated gradient is zero, so a derivative would be a silent zero), and solvers check their element's rule against the mesh rule when they attach the auxiliary vector (off-rule the field reads as zero). Base class hooks: _create_petsc_fe and _basis_key on _BaseMeshVariable; the mesh coordinate cache keys on _basis_key; EnhancedMeshVariable takes its storage class from _base_variable_class. tests/test_0065: layout on triangle/tet/quad; integral of random point data equals the hand quadrature sum (with a one-point perturbation control); P2 projection of P2 point data exact to solver tolerance; evaluate exact at own points and equal to the nearest-in-cell rule elsewhere (a P1 interpolant does not match, as the control); both guards. Docs: docs/developer/subsystems/integration-point-variables.md. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- docs/developer/index.md | 1 + .../subsystems/integration-point-variables.md | 112 ++++++++++++ .../cython/petsc_generic_snes_solvers.pyx | 4 + .../cython/petsc_quadrature_fe.pyx | 41 +++++ src/underworld3/discretisation/__init__.py | 1 + .../discretisation/discretisation_mesh.py | 41 ++++- .../discretisation_mesh_variables.py | 146 +++++++++++++++- .../discretisation/enhanced_variables.py | 74 +++++++- src/underworld3/function/_function.pyx | 17 ++ src/underworld3/utilities/_jitextension.py | 18 +- tests/test_0065_integration_point_variable.py | 159 ++++++++++++++++++ 11 files changed, 600 insertions(+), 14 deletions(-) create mode 100644 docs/developer/subsystems/integration-point-variables.md create mode 100644 tests/test_0065_integration_point_variable.py diff --git a/docs/developer/index.md b/docs/developer/index.md index b88e30519..dc10dbc0f 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -189,6 +189,7 @@ subsystems/meshing subsystems/mesh-shape-relaxation subsystems/conforming-surfaces-and-fault-zones subsystems/discretisation +subsystems/integration-point-variables subsystems/solvers subsystems/boundary-stress-and-projection-postprocessing subsystems/rotated-freeslip diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md new file mode 100644 index 000000000..01c15a59d --- /dev/null +++ b/docs/developer/subsystems/integration-point-variables.md @@ -0,0 +1,112 @@ +# Integration-point variables + +`uw.discretisation.IntegrationPointVariable` stores one value per quadrature +point per cell, on an element whose basis is the identity on the mesh's +integration rule. The assembler reads the stored values directly at the +integration points, with no interpolation. It is a peer of `MeshVariable` +and of swarm variables: a nodal field is *sampled*, a swarm carries +*particles*, and an integration-point variable carries values that are +*injected* into the weak form exactly where it is evaluated. + +```python +import underworld3 as uw + +mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) +eta_q = uw.discretisation.IntegrationPointVariable("eta_q", mesh) + +eta_q.cell_data.shape # (ncells, Nq, 1) +eta_q.coords # the physical integration points, same order +eta_q.cell_data[...] = 1.0 +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_q.sym +``` + +## Why it exists + +Two uses drove the design. + +**Semi-Lagrangian history.** The SLCN scheme samples the previous solution at +the departure point of every node, builds a nodal history field, and the +assembler then interpolates that field to the quadrature points. Two +interpolations per step. If the departure points are traced from the +integration points instead, and the sampled values are stored here, the value +entering the weak form is the discrete solution evaluated exactly at the +departure point. Only the FE solution's own error remains. + +**Material properties from particles.** A swarm property normally reaches the +constitutive law through a proxy mesh variable (nearest-neighbour or RBF to +the nodes, then the basis to the integration points), which smooths a +material interface over a cell. Reconstructing the property at each +integration point from the particles near it and writing it here keeps +sub-cell contrast, as the integration swarms of Underworld 1 and 2 did. + +## The rule is a mesh property + +Every field on a mesh is created on the rule fixed by `mesh.qdegree`, and +PETSc's `PetscDSSetUp` tabulates all fields of a discrete system on one rule. +So the element is built once per mesh, on `mesh.integration_rule`, whatever +the degrees of the fields it sits beside. A P1 temperature and a P2 velocity +on a `qdegree=2` triangle mesh are both integrated on the six-point rule, and +an integration-point variable on that mesh has six values per cell. + +Point counts at `qdegree=2`: triangle 6, tetrahedron 14, quadrilateral 9, +hexahedron 27. + +## What `evaluate` means + +A delta field is defined only at its points. Between them it is *defined* +as piecewise constant on the nearest-integration-point partition of each +cell, and that is what `uw.function.evaluate` returns: locate the cell, take +the closest of its points. This is the one extension under which a query +agrees with what the assembler used at that point; an interpolant or a +projection would report a different viscosity from the one the solver saw. +Points no cell owns take the nearest integration point on the rank. + +Evaluating the variable at its own `coords` returns its own `data` exactly. + +If a smooth nodal picture is wanted (a plot, a diagnostic), project the +symbol onto a `MeshVariable` explicitly with `SNES_Projection`; the +projection of the field is an ordinary weak form and is exact for data that +the target space can represent. + +## Guards + +The field has no gradient (its tabulated derivative is identically zero), so +a derivative of its symbol in a weak form would be a silent zero. The JIT +refuses it at code generation: + +``` +RuntimeError: {h}_{,0}: derivative of an integration-point variable has no meaning ... +``` + +Off its own rule the field tabulates to zero, so a solver on a different +rule would drop the term it carries. A solver that attaches the mesh's +auxiliary vector checks its element against `mesh.integration_rule` and +raises if they differ. Boundary integrals evaluate the field on the face +rule and see zeros; that is correct for a history term and worth knowing for +anything else. + +Scalar components only for now; use one variable per component. + +## Implementation + +- `src/underworld3/cython/uw_delta_space.h`: the `uwdelta` `PetscSpace` + type, registered with `PetscSpaceRegister`. PETSc's own `PETSCSPACEPOINT` + cannot be tabulated anywhere but its own points in its own order, which + breaks `PetscFESetUp`, face tabulation and boundary integrals; the plugin + type returns zeros off its rule and needs no PETSc patch, so stock conda + PETSc works. +- `src/underworld3/cython/petsc_quadrature_fe.pyx`: `create_delta_fe` + (the element, via `PetscFECreateFromSpaces` with a `PETSCDUALSPACESIMPLE` + dual space of one point evaluation per rule point), `tabulate` (for + tests) and `cell_quadrature_points` (the physical integration points from + `DMPlexComputeCellGeometryFEM`, the assembler's own map). +- `_BaseIntegrationPointVariable` in `discretisation_mesh_variables.py` + overrides the two discretisation hooks (`_create_petsc_fe`, `_basis_key`) + and supplies the nearest-point evaluation; `IntegrationPointVariable` in + `enhanced_variables.py` is the public wrapper. +- All dofs sit on the cell interior, so the local vector is cell-major, + point-minor, and `cell_data` is a plain reshape. + +Tests: `tests/test_0064_quadrature_point_fe.py` (the element), +`tests/test_0065_integration_point_variable.py` (the variable, the assembler +reading it, `evaluate`, the guards). diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f56fbe0fe..aa22e699e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -709,6 +709,7 @@ class SolverBaseClass(uw_object): current field values (callbacks may have changed v, p, or auxiliary fields).""" self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self.mesh._verify_integration_rule(getattr(self, "petsc_fe_u", None)) def _dispatch_snes_update(self, snes, iteration): """PETSc SNESSetUpdate hook: sync iterate->fields, run callbacks, sync back. @@ -3108,6 +3109,7 @@ class SolverBaseClass(uw_object): self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self.mesh._verify_integration_rule(getattr(self, "petsc_fe_u", None)) self._update_constants() gvec = self.dm.getGlobalVec() @@ -9089,6 +9091,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self.mesh._verify_integration_rule(getattr(self, "petsc_fe_u", None)) self._update_constants() gvec = self.dm.getGlobalVec() @@ -9233,6 +9236,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) + self.mesh._verify_integration_rule(getattr(self, "petsc_fe_u", None)) self._update_constants() gvec = self.dm.getGlobalVec() diff --git a/src/underworld3/cython/petsc_quadrature_fe.pyx b/src/underworld3/cython/petsc_quadrature_fe.pyx index 1497c3f56..efbbd5531 100644 --- a/src/underworld3/cython/petsc_quadrature_fe.pyx +++ b/src/underworld3/cython/petsc_quadrature_fe.pyx @@ -77,6 +77,9 @@ cdef extern from "petsc.h" nogil: ctypedef _n_PetscTabulation* PetscTabulation PetscErrorCode PetscFECreateTabulation(PetscFE, PetscInt, PetscInt, const PetscReal*, PetscInt, PetscTabulation*) PetscErrorCode PetscTabulationDestroy(PetscTabulation*) + PetscErrorCode DMPlexComputeCellGeometryFEM(PetscDM, PetscInt, PetscQuadrature, PetscReal*, PetscReal*, PetscReal*, PetscReal*) + PetscErrorCode DMPlexGetHeightStratum(PetscDM, PetscInt, PetscInt*, PetscInt*) + PetscErrorCode DMGetCoordinateDim(PetscDM, PetscInt*) cdef extern from "uw_delta_space.h" nogil: PetscErrorCode UWDeltaSpaceRegister() @@ -190,3 +193,41 @@ def tabulate(FE fe, points, int K=0): ov[p, b, c] = T.T[0][(p * Nb + b) * Nc + c] CHKERRQ(PetscTabulationDestroy(&T)) return out + + +def cell_quadrature_points(DM dm, Quad quad): + r"""Physical coordinates of the rule's points in every local cell. + + Returns an array shaped ``(ncells, Nq, cdim)`` in local cell order, computed + by ``DMPlexComputeCellGeometryFEM`` - the same map the assembler uses for + its integration points, so row ``(c, q)`` is exactly where the pointwise + functions see quadrature point ``q`` of cell ``c``. No locator involved. + """ + cdef PetscInt cStart = 0, cEnd = 0, cdim = 0, Nq = 0, c, q, d + cdef PetscReal *v = NULL + cdef PetscReal *J = NULL + cdef PetscReal *invJ = NULL + cdef PetscReal *detJ = NULL + CHKERRQ(DMPlexGetHeightStratum(dm.dm, 0, &cStart, &cEnd)) + CHKERRQ(DMGetCoordinateDim(dm.dm, &cdim)) + CHKERRQ(PetscQuadratureGetData(quad.quad, NULL, NULL, &Nq, NULL, NULL)) + ncells = cEnd - cStart + out = np.empty((ncells, Nq, cdim), dtype=np.float64) + cdef double[:, :, ::1] ov = out + vbuf = np.empty(Nq * cdim, dtype=np.float64) + Jbuf = np.empty(Nq * cdim * cdim, dtype=np.float64) + iJbuf = np.empty(Nq * cdim * cdim, dtype=np.float64) + dJbuf = np.empty(Nq, dtype=np.float64) + cdef double[::1] vv = vbuf + cdef double[::1] Jv = Jbuf + cdef double[::1] iJv = iJbuf + cdef double[::1] dJv = dJbuf + if ncells == 0: + return out + v = &vv[0]; J = &Jv[0]; invJ = &iJv[0]; detJ = &dJv[0] + for c in range(cStart, cEnd): + CHKERRQ(DMPlexComputeCellGeometryFEM(dm.dm, c, quad.quad, v, J, invJ, detJ)) + for q in range(Nq): + for d in range(cdim): + ov[c - cStart, q, d] = v[q * cdim + d] + return out diff --git a/src/underworld3/discretisation/__init__.py b/src/underworld3/discretisation/__init__.py index 687b34513..a60150b48 100644 --- a/src/underworld3/discretisation/__init__.py +++ b/src/underworld3/discretisation/__init__.py @@ -19,6 +19,7 @@ """ from .discretisation_mesh import Mesh from .enhanced_variables import EnhancedMeshVariable as MeshVariable +from .enhanced_variables import IntegrationPointVariable from .discretisation_mesh import checkpoint_xdmf from .discretisation_mesh import meshVariable_lookup_by_symbol from .discretisation_mesh import petsc_dm_find_labeled_points_local diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 56a5222a0..ef2686c21 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5476,14 +5476,51 @@ def _get_coords_for_var(self, var): provided variable. If the array does not already exist, it is first created and then returned. """ - key = (self.isSimplex, var.degree, var.continuous) + key = var._basis_key # if array already created, return. if key in self._coord_array: return self._coord_array[key] + if getattr(var, "is_integration_point", False): + # Cell-major, point-minor: the layout of the variable's own vector. + self._coord_array[key] = var.integration_points.reshape(-1, self.cdim).copy() else: self._coord_array[key] = self._get_coords_for_basis(var.degree, var.continuous) - return self._coord_array[key] + return self._coord_array[key] + + @property + def integration_rule(self): + """The cell quadrature rule every field on this mesh is integrated on. + + Fixed by ``qdegree`` (PETSc's ``PetscDSSetUp`` forces one rule per + discrete system), so it is a property of the mesh, not of any field. + Integration-point variables are built on it. + """ + if getattr(self, "_integration_rule", None) is None: + fe = PETSc.FE().createDefault( + self.dim, 1, self.isSimplex, self.qdegree, "integration_rule_", PETSc.COMM_SELF, + ) + self._integration_rule = fe.getQuadrature() + self._integration_rule_fe = fe # keeps the rule alive + return self._integration_rule + + def _verify_integration_rule(self, fe): + """Raise unless ``fe`` integrates on this mesh's rule. + + An integration-point variable reads as zeros on any other rule, which + would silently drop the term it carries, so a solver that attaches the + mesh's auxiliary vector checks its own element here. + """ + if fe is None or not any(getattr(v, "is_integration_point", False) for v in self.vars.values()): + return + q_mesh = numpy.asarray(self.integration_rule.getData()[0]).reshape(-1, self.dim) + q_fe = numpy.asarray(fe.getQuadrature().getData()[0]).reshape(-1, self.dim) + if q_mesh.shape != q_fe.shape or not numpy.allclose(q_mesh, q_fe, atol=1e-12): + raise RuntimeError( + f"Solver quadrature ({q_fe.shape[0]} points) differs from the mesh integration " + f"rule ({q_mesh.shape[0]} points, qdegree={self.qdegree}); integration-point " + "variables on this mesh would read as zero. Build the solver on mesh.qdegree." + ) def _basis_coordinate_dm(self, degree, continuous): """Coordinate DM carrying a degree-``degree`` Lagrange field. diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 5d8eb7d66..4a277851d 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1704,6 +1704,27 @@ def _data_layout(self, i, j=None): else: return i + j * self.shape[0] + # Discretisation hooks. Subclasses with a different element (the + # integration-point variable) override these three; everything else in + # the class works from the PETSc field they produce. + is_integration_point = False + + @property + def _basis_key(self): + """Key for the mesh's per-basis coordinate cache.""" + return (self.mesh.isSimplex, self.degree, self.continuous) + + def _create_petsc_fe(self, dim, prefix): + """The PetscFE this variable's field is built on (Lagrange by default).""" + return PETSc.FE().createDefault( + dim, + self.num_components, + self.mesh.isSimplex, + self.mesh.qdegree, + prefix, + PETSc.COMM_SELF, + ) + def _setup_ds(self): options = PETSc.Options() name0 = "VAR" # self.clean_name ## Filling up the options database @@ -1714,14 +1735,7 @@ def _setup_ds(self): ) # only active if discontinuous dim = self.mesh.dm.getDimension() - petsc_fe = PETSc.FE().createDefault( - dim, - self.num_components, - self.mesh.isSimplex, - self.mesh.qdegree, - name0 + "_", - PETSc.COMM_SELF, - ) + petsc_fe = self._create_petsc_fe(dim, name0 + "_") # Check if this is the first field or if we need to rebuild the DM # (needed to ensure Section is properly synchronized with field list) @@ -3488,3 +3502,119 @@ def jacobian(self): # Note: EnhancedMeshVariable is imported as MeshVariable in __init__.py to avoid circular imports + + +class _BaseIntegrationPointVariable(_BaseMeshVariable): + r"""A field stored at the mesh integration points (quadrature rule). + + One degree of freedom per quadrature point per cell, on the element built by + :func:`underworld3.cython.petsc_quadrature_fe.create_delta_fe`: the basis is + the identity on the mesh rule, so the pointwise functions read the stored + value at each integration point with no interpolation. Values are + *injected* here (a semi-Lagrangian history, a material property + reconstructed from a swarm); the field is a peer of mesh and swarm + variables, not a degree-0 mesh variable. + + Between its points the field is defined as piecewise constant on the + nearest-integration-point partition of each cell. That is what + ``evaluate()`` returns, and it is the only extension under which a query + agrees with what the assembler used at that point. + + Layout: ``data`` is ``(ncells * Nq, num_components)`` in local cell order, + point-minor; ``cell_data`` views it as ``(ncells, Nq, num_components)`` and + ``coords`` are the physical integration points in the same order. + + Derivatives of the symbol are meaningless (the tabulated gradient is zero) + and the JIT refuses them. Scalar components only for now. + """ + + is_integration_point = True + + def __init__(self, varname=None, mesh=None, num_components=None, vtype=None, + varsymbol=None, _register=True, units=None, units_backend=None, + remesh_policy=None, **kwargs): + # degree/continuous are not meaningful here; 0/False keeps the base + # class's bookkeeping consistent with a cell-interior field. + kwargs.pop("degree", None) + kwargs.pop("continuous", None) + self._ip_coords_cache = None + super().__init__(varname=varname, mesh=mesh, num_components=num_components, + vtype=vtype, degree=0, continuous=False, varsymbol=varsymbol, + _register=_register, units=units, units_backend=units_backend, + remesh_policy=remesh_policy, **kwargs) + + # -- discretisation hooks ------------------------------------------------- + + @property + def _basis_key(self): + return ("integration", self.mesh.isSimplex, self.mesh.qdegree) + + def _create_petsc_fe(self, dim, prefix): + from underworld3.cython.petsc_quadrature_fe import create_delta_fe + if self.num_components != 1: + raise NotImplementedError( + "IntegrationPointVariable: scalar components only for now - " + "use one variable per component" + ) + cStart, _ = self.mesh.dm.getHeightStratum(0) + fe = create_delta_fe(self.mesh.integration_rule, self.mesh.dm.getCellType(cStart), + name=f"{prefix}integration_point_fe") + return fe + + # -- geometry --------------------------------------------------------------- + + @property + def integration_points(self): + """Physical integration points, ``(ncells, Nq, cdim)``, local cell order.""" + if self._ip_coords_cache is None or self._ip_coords_cache[0] != self.mesh._topology_version: + from underworld3.cython.petsc_quadrature_fe import cell_quadrature_points + pts = cell_quadrature_points(self.mesh.dm, self.mesh.integration_rule) + self._ip_coords_cache = (self.mesh._topology_version, pts) + return self._ip_coords_cache[1] + + @property + def num_points_per_cell(self): + return self.integration_points.shape[1] + + @property + def cell_data(self): + """``data`` viewed as ``(ncells, Nq, num_components)``.""" + Nq = self.num_points_per_cell + return self.data.reshape(-1, Nq, self.num_components) + + # -- evaluation --------------------------------------------------------------- + + def _nearest_point_values(self, coords_nd, cells): + """Values at ``coords_nd`` by the nearest integration point of the + owning cell ``cells`` (local index; -1 or None means unowned -> the + nearest point anywhere on this rank).""" + coords_nd = numpy.asarray(coords_nd, dtype=float).reshape(-1, self.mesh.cdim) + n = coords_nd.shape[0] + vals = numpy.empty((n, self.num_components), dtype=float) + if n == 0: + return vals + ipc = self.integration_points + cdat = numpy.asarray(self.data).reshape(ipc.shape[0], ipc.shape[1], self.num_components) + cells = None if cells is None else numpy.asarray(cells).reshape(-1) + owned = numpy.ones(n, dtype=bool) if cells is None else (cells >= 0) + if cells is None: + owned[:] = False + if owned.any(): + cc = cells[owned] + d2 = ((ipc[cc] - coords_nd[owned][:, None, :]) ** 2).sum(axis=-1) + j = d2.argmin(axis=1) + vals[owned] = cdat[cc, j] + if (~owned).any(): + vals[~owned] = self.rbf_interpolate(coords_nd[~owned]) + return vals + + def rbf_interpolate(self, new_coords, nnn=None, p=1, verbose=False, **kwargs): + """Nearest integration point on this rank (the exterior / unowned-point rule).""" + new_coords = numpy.asarray(new_coords, dtype=float).reshape(-1, self.mesh.cdim) + ipc = self.integration_points.reshape(-1, self.mesh.cdim) + if ipc.shape[0] == 0: + return numpy.full((new_coords.shape[0], self.num_components), numpy.nan) + import underworld3 as uw + tree = uw.kdtree.KDTree(ipc) + _, idx = tree.query(new_coords, k=1) + return numpy.asarray(self.data).reshape(-1, self.num_components)[numpy.asarray(idx).reshape(-1)] diff --git a/src/underworld3/discretisation/enhanced_variables.py b/src/underworld3/discretisation/enhanced_variables.py index 4ff8861ce..385f60cc7 100644 --- a/src/underworld3/discretisation/enhanced_variables.py +++ b/src/underworld3/discretisation/enhanced_variables.py @@ -28,7 +28,7 @@ from typing import Optional, Union import numpy as np -from .discretisation_mesh_variables import _BaseMeshVariable +from .discretisation_mesh_variables import _BaseMeshVariable, _BaseIntegrationPointVariable from ..utilities import MathematicalMixin from ..utilities.dimensionality_mixin import DimensionalityMixin @@ -63,6 +63,10 @@ class EnhancedMeshVariable(DimensionalityMixin, MathematicalMixin): success = persistent_var.transfer_data_from(old_pressure) """ + # The storage class this wrapper delegates to; IntegrationPointVariable + # swaps in the quadrature-point element. + _base_variable_class = _BaseMeshVariable + def __new__(cls, varname, mesh, *args, **kwargs): """Custom __new__ to ensure proper initialization and registration.""" # Create the instance @@ -127,7 +131,7 @@ def __init__( self._mesh_ref = weakref.ref(mesh) # Weak reference to avoid circular deps # Create base variable without registration (we handle registration ourselves) - self._base_var = _BaseMeshVariable( + self._base_var = self._base_variable_class( varname=varname, mesh=mesh, num_components=num_components, @@ -926,3 +930,69 @@ def demonstrate_enhanced_variables(): # Note: The demonstration function above references EnhancedSwarmVariable # which doesn't exist - SwarmVariable is already enhanced (see swarm.py). # Update this demo to use uw.swarm.SwarmVariable directly if needed. + + +class IntegrationPointVariable(EnhancedMeshVariable): + r"""A field stored at the mesh integration points. + + A peer of :class:`MeshVariable` and of swarm variables: one value per + quadrature point per cell, on an element whose basis is the identity on the + mesh rule (``mesh.integration_rule``). The pointwise functions read the + stored values directly, with no interpolation, so it is the carrier for + values that are *injected* at the integration points - a semi-Lagrangian + history sampled at the departure points of the quadrature points, or a + material property reconstructed from a swarm with sub-cell resolution. + + Between its points the field is piecewise constant on the + nearest-integration-point partition of each cell; ``uw.function.evaluate`` + returns that, so a query agrees with what the assembler used at the same + point. Derivatives of the symbol are refused by the JIT (the gradient is + identically zero). Scalar components only for now. + + Examples + -------- + >>> eta_q = uw.discretisation.IntegrationPointVariable("eta_q", mesh) + >>> eta_q.cell_data[...] = 1.0 # (ncells, Nq, 1) + >>> eta_q.coords # the physical integration points + >>> stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_q.sym + """ + + _base_variable_class = _BaseIntegrationPointVariable + + def __init__( + self, + varname, + mesh, + num_components=1, + vtype=None, + varsymbol=None, + persistent=False, + units=None, + units_backend=None, + **kwargs, + ): + kwargs.pop("degree", None) + kwargs.pop("continuous", None) + super().__init__( + varname, mesh, num_components=num_components, vtype=vtype, + degree=0, continuous=False, varsymbol=varsymbol, persistent=persistent, + units=units, units_backend=units_backend, **kwargs, + ) + + # Explicit passthroughs (the wrapper delegates unknown attributes to the + # sympy matrix, not to the storage object). + is_integration_point = True + + @property + def integration_points(self): + """Physical integration points, ``(ncells, Nq, cdim)``, local cell order.""" + return self._base_var.integration_points + + @property + def num_points_per_cell(self): + return self._base_var.num_points_per_cell + + @property + def cell_data(self): + """``data`` viewed as ``(ncells, Nq, num_components)``.""" + return self._base_var.cell_data diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 21fcf641e..6a144f252 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -1449,6 +1449,23 @@ def petsc_interpolate( expr, rbf_vals = np.asarray(var.rbf_interpolate(fallback_coords)) rbf_vals = rbf_vals.reshape(len(fallback_coords), var.num_components) outarray[unlocated, var_start:var_start + var.num_components] = rbf_vals + + # Integration-point variables: the FE interpolation above tabulates + # their delta basis at the query points, which is zero anywhere but + # on the rule. Their defined extension is the nearest integration + # point of the owning cell; overwrite their columns with it. + ip_vars = [v for v in vars if getattr(v, "is_integration_point", False)] + if ip_vars: + ip_cells = getattr(cached_info, "cells", None) + if ip_cells is None: + ip_cells = mesh._robust_owning_cells(coords) + ip_cells = np.asarray(ip_cells).reshape(-1).copy() + if unlocated is not None: + ip_cells[np.asarray(unlocated, dtype=bool)] = -1 + for var in ip_vars: + var_start = var_start_index[var] + outarray[:, var_start:var_start + var.num_components] = \ + var._nearest_point_values(coords, ip_cells) # === END CACHING === # Create map between array slices and variable functions diff --git a/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 7d63ae183..63a6400ec 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -877,7 +877,21 @@ def ccode_patch_fns(varlist, prefix_str): u_i = 0 # variable increment u_x_i = 0 # variable gradient increment lambdafunc = lambda self, printer: self._ccodestr + + def _no_derivative(self, printer): + # An integration-point variable has no gradient (its tabulated + # derivative is identically zero), so a derivative of its symbol + # in a weak form would be a silent zero. Refuse at code generation. + raise RuntimeError( + f"{self.__class__.__name__}: derivative of an integration-point " + "variable has no meaning (the field is defined only at the " + "quadrature points). Remove the derivative or project the " + "variable onto a nodal MeshVariable first." + ) + for var in varlist: + is_ip = getattr(var, "is_integration_point", False) + dfunc = _no_derivative if is_ip else lambdafunc if var.vtype == VarType.SCALAR: # monkey patch this guy into the function type(var.fn)._ccodestr = f"{prefix_str}[{u_i}]" @@ -893,7 +907,7 @@ def ccode_patch_fns(varlist, prefix_str): for ind in range(mesh.cdim): # Note that var.fn._diff[ind] returns the class, so we don't need type(var.fn._diff[ind]) var.fn._diff[ind]._ccodestr = f"{prefix_str}_x[{u_x_i}]" - var.fn._diff[ind]._ccode = lambdafunc + var.fn._diff[ind]._ccode = dfunc u_x_i += 1 elif ( var.vtype == VarType.VECTOR @@ -912,7 +926,7 @@ def ccode_patch_fns(varlist, prefix_str): for ind in range(mesh.cdim): # Note that var.fn._diff[ind] returns the class, so we don't need type(var.fn._diff[ind]) comp._diff[ind]._ccodestr = f"{prefix_str}_x[{u_x_i}]" - comp._diff[ind]._ccode = lambdafunc + comp._diff[ind]._ccode = dfunc u_x_i += 1 else: raise RuntimeError( diff --git a/tests/test_0065_integration_point_variable.py b/tests/test_0065_integration_point_variable.py new file mode 100644 index 000000000..1742ae8f7 --- /dev/null +++ b/tests/test_0065_integration_point_variable.py @@ -0,0 +1,159 @@ +"""IntegrationPointVariable: a field stored at the mesh integration points. + +What is checked, and why each check is the one that matters: + +- layout: one value per rule point per cell, coordinates from the assembler's + own cell geometry (the rule-weighted mean of a cell's points is its centroid); +- the assembler reads the stored values exactly: the integral of random + point data equals the quadrature sum done by hand, and a P2 projection of + P2 point data is exact to solver tolerance; +- ``evaluate`` is the nearest integration point of the owning cell, exact at + the variable's own points; +- the two guards: a derivative of the symbol is refused by the JIT, and a + solver on a different rule is refused by the mesh. +""" + +import numpy as np +import pytest +import sympy +from petsc4py import PETSc + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _mesh(kind): + if kind == "triangle": + return uw.meshing.UnstructuredSimplexBox(cellSize=0.25, qdegree=2) + if kind == "tetrahedron": + return uw.meshing.UnstructuredSimplexBox( + minCoords=(0, 0, 0), maxCoords=(1, 1, 1), cellSize=0.5, qdegree=2 + ) + if kind == "quadrilateral": + return uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2) + raise ValueError(kind) + + +def _ncells(mesh): + c0, c1 = mesh.dm.getHeightStratum(0) + return c1 - c0 + + +@pytest.mark.parametrize("kind", ["triangle", "tetrahedron", "quadrilateral"]) +def test_layout_and_geometry(kind): + mesh = _mesh(kind) + h = uw.discretisation.IntegrationPointVariable("h", mesh) + Nq = len(np.asarray(mesh.integration_rule.getData()[1])) + n = _ncells(mesh) + + assert h.is_integration_point + assert h.num_points_per_cell == Nq + assert h.data.shape == (n * Nq, 1) + assert h.cell_data.shape == (n, Nq, 1) + assert np.asarray(h.coords).shape == (n * Nq, mesh.dim) + assert np.array_equal(np.asarray(h.coords_nd), h.integration_points.reshape(-1, mesh.dim)) + + # Affine cells: the rule-weighted mean of the points is the centroid. + w = np.asarray(mesh.integration_rule.getData()[1]).reshape(-1) + cent = (h.integration_points * w[None, :, None]).sum(1) / w.sum() + assert np.allclose(cent, np.asarray(mesh._centroids)[:n], atol=1e-12) + + +def test_assembler_reads_the_stored_values(): + """Integral of random point data == the quadrature sum done by hand.""" + mesh = _mesh("triangle") + h = uw.discretisation.IntegrationPointVariable("h", mesh) + rng = np.random.default_rng(1) + h.data[:, 0] = rng.uniform(-1.0, 2.0, size=h.data.shape[0]) + + # Cell areas from the vertices, rule weights scaled by area / reference area. + verts = np.asarray(mesh._get_coords_for_basis(1, True)) + rows = np.asarray(mesh._cell_node_indices(1, True)) + p = verts[rows] # (ncells, 3, 2) + area = 0.5 * np.abs( + (p[:, 1, 0] - p[:, 0, 0]) * (p[:, 2, 1] - p[:, 0, 1]) + - (p[:, 2, 0] - p[:, 0, 0]) * (p[:, 1, 1] - p[:, 0, 1]) + ) + w = np.asarray(mesh.integration_rule.getData()[1]).reshape(-1) + by_hand = ((h.cell_data[:, :, 0] * w[None, :]).sum(1) * area / w.sum()).sum() + + assembled = uw.maths.Integral(mesh, h.sym[0]).evaluate() + assert abs(assembled - by_hand) < 1e-12 * max(1.0, abs(by_hand)) + # Negative control: perturb one point and the integral must move by + # exactly that point's weight. + c, q = 3, 2 + h.cell_data[c, q, 0] += 1.0 + moved = uw.maths.Integral(mesh, h.sym[0]).evaluate() + assert abs((moved - assembled) - w[q] * area[c] / w.sum()) < 1e-12 + + +def test_projection_of_p2_point_data_is_exact(): + mesh = _mesh("triangle") + x, y = mesh.X + h = uw.discretisation.IntegrationPointVariable("h", mesh) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) + f = lambda X: 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + X[:, 0] * X[:, 1] + h.data[:, 0] = f(np.asarray(h.coords)) + + proj = uw.systems.solvers.SNES_Projection(mesh, T) + proj.uw_function = h.sym[0] + proj.smoothing = 0.0 + proj.petsc_options["ksp_rtol"] = 1e-13 + proj.petsc_options["snes_rtol"] = 1e-13 + proj.solve() + assert np.abs(T.data[:, 0] - f(np.asarray(T.coords))).max() < 1e-9 + + +def test_evaluate_is_nearest_point_of_owning_cell(): + mesh = _mesh("triangle") + h = uw.discretisation.IntegrationPointVariable("h", mesh) + rng = np.random.default_rng(2) + h.data[:, 0] = rng.uniform(size=h.data.shape[0]) + + # Exact at its own points. + own = uw.function.evaluate(h.sym[0], np.asarray(h.coords)).reshape(-1) + assert np.array_equal(own, h.data[:, 0]) + + # Nearest point of the owning cell elsewhere. + pts = rng.uniform(0.05, 0.95, size=(300, 2)) + cells = mesh._robust_owning_cells(pts) + ipc = h.integration_points + j = ((ipc[cells] - pts[:, None, :]) ** 2).sum(-1).argmin(1) + expected = h.cell_data[cells, j, 0] + got = uw.function.evaluate(h.sym[0], pts).reshape(-1) + assert np.array_equal(got, expected) + + # Negative control: a nodal P1 interpolant of the same data does not + # match this definition (it smooths), so the test discriminates. + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=1) + T.data[:, 0] = uw.function.evaluate(h.sym[0], np.asarray(T.coords)).reshape(-1) + smooth = uw.function.evaluate(T.sym[0], pts).reshape(-1) + assert not np.allclose(smooth, expected) + + +def test_derivative_is_refused_by_the_jit(): + mesh = _mesh("triangle") + x, y = mesh.X + h = uw.discretisation.IntegrationPointVariable("h", mesh) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=1) + proj = uw.systems.solvers.SNES_Projection(mesh, T) + proj.uw_function = h.sym[0].diff(x) + with pytest.raises(RuntimeError, match="integration-point"): + proj.solve() + + +def test_other_rule_is_refused(): + mesh = _mesh("triangle") + uw.discretisation.IntegrationPointVariable("h", mesh) + same = PETSc.FE().createDefault(2, 1, True, mesh.qdegree, "same_", PETSc.COMM_SELF) + other = PETSc.FE().createDefault(2, 1, True, mesh.qdegree + 1, "other_", PETSc.COMM_SELF) + mesh._verify_integration_rule(same) + with pytest.raises(RuntimeError, match="integration rule"): + mesh._verify_integration_rule(other) + + +def test_vector_components_not_yet_supported(): + mesh = _mesh("triangle") + with pytest.raises(NotImplementedError): + uw.discretisation.IntegrationPointVariable("v", mesh, num_components=2) From 90857d9619b1ce9186f01fb98e34b312917e860e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 14:33:46 -0700 Subject: [PATCH 03/15] IntegrationPointSemiLagrangian: SLCN history at the integration points The history slots are IntegrationPointVariables, so the advected value the weak form sees at each integration point is the solution from k+1 steps ago evaluated exactly at the departure point of that integration point. No nodal history field, no second interpolation. A delta field cannot be sampled off its points, so the slot-to-slot chain of SemiLagrangian is replaced by nodal snapshots of the solution and velocity at the last 'order' times; slot k is filled by tracing k+1 RK2 segments back from every integration point (segment j with the velocity at time n-j and that step's dt) and evaluating the snapshot from time n-k at the foot. Each slot carries one evaluation error rather than one per generation. Drops in as DuDt for AdvDiffusion; the flux history DFDt stays nodal. tests/test_0066: for a P2 field in a uniform velocity both slots reproduce the exact departure-point values to 1e-12 (the nodal scheme does not, as the control); rotating Gaussian is at least as accurate as nodal SLCN and keeps the peak. Measured at cellSize 0.05, half a revolution: Courant 1 L2 3.6e-3 -> 3.2e-3, peak 0.987 -> 0.9998; Courant 2 equal (dt-dominated). Scalar only; no ALE, no checkpoint state. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 31 +++ src/underworld3/systems/__init__.py | 1 + src/underworld3/systems/ddt.py | 231 ++++++++++++++++++ tests/test_0066_integration_point_slcn.py | 95 +++++++ 4 files changed, 358 insertions(+) create mode 100644 tests/test_0066_integration_point_slcn.py diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 01c15a59d..bffb88a4d 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -110,3 +110,34 @@ Scalar components only for now; use one variable per component. Tests: `tests/test_0064_quadrature_point_fe.py` (the element), `tests/test_0065_integration_point_variable.py` (the variable, the assembler reading it, `evaluate`, the guards). + +## Semi-Lagrangian history on the integration points + +`uw.systems.ddt.IntegrationPointSemiLagrangian` is the SLCN history built on +this variable. Its slots `psi_star[k]` are integration-point variables, so +the value the weak form sees at each integration point is the solution from +`k+1` steps ago evaluated exactly at the departure point of that +integration point. A nodal history cannot be sampled from a delta field, so +the slot-to-slot chain of `SemiLagrangian` is replaced by nodal snapshots of +the solution and of the velocity at the last `order` times: slot `k` is +filled by tracing `k+1` RK2 segments back from every integration point, +segment `j` with the velocity at time `n-j` and that step's `dt`, and +evaluating the snapshot from time `n-k` at the foot. Every slot carries one +evaluation error rather than one per generation. + +```python +DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V_fn, degree=2, order=1) +adv = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=V_fn, DuDt=DuDt, order=1) +``` + +The diffusive flux history (`DFDt`) keeps its nodal projection, since it +carries derivatives. Scalar histories only; no ALE or old-frame trace-back, +no checkpoint state yet. + +For a P2 field in a uniform velocity the slots reproduce the exact +departure-point values to round-off, for one and for two segments +(`tests/test_0066_integration_point_slcn.py`). On a rotating Gaussian at +Courant 1 (`cellSize=0.05`, `dt=0.1`, half a revolution) the L2 error drops +from 3.6e-3 (nodal SLCN) to 3.2e-3 and the peak is kept at 0.9998 instead +of 0.987. The trace-back samples six points per cell rather than the P2 +nodes, so the update costs about twice the nodal one. diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index ea6823ab5..752aa2df9 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -85,6 +85,7 @@ # are the Lagrangian implementations actually distinct in reality ? from .ddt import Lagrangian as Lagrangian_DDt from .ddt import SemiLagrangian as SemiLagragian_DDt +from .ddt import IntegrationPointSemiLagrangian as IntegrationPointSemiLagrangian_DDt from .ddt import Lagrangian_Swarm as Lagrangian_Swarm_DDt from .ddt import Eulerian as Eulerian_DDt diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e198aad57..42b6fa66f 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -3473,3 +3473,234 @@ def update_post_solve( return + + +class IntegrationPointSemiLagrangian(_DDtBase): + r"""Semi-Lagrangian history stored at the mesh integration points. + + The history slots ``psi_star[k]`` are + :class:`~underworld3.discretisation.IntegrationPointVariable` objects, so + the value the weak form sees at each integration point is the discrete + solution from ``k+1`` steps ago evaluated **exactly** at the departure + point of that integration point. There is no nodal history field and no + second interpolation: only the FE solution's own error remains in the + advected term. Compare :class:`SemiLagrangian`, which samples at the + nodes, stores a nodal ``psi_star`` and lets the assembler interpolate it + to the integration points. + + Because a delta field cannot be sampled off its points, the chain + ``psi_star[k] <- psi_star[k-1]`` of :class:`SemiLagrangian` is replaced + by nodal **snapshots** of the solution and of the velocity at the last + ``order`` times. Slot ``k`` is filled by tracing ``k+1`` segments back + from every integration point (segment ``j`` with the velocity at time + ``n-j`` and that step's ``dt``) and evaluating the snapshot from time + ``n-k`` at the foot. Every slot carries one evaluation error rather than + one per generation. + + What is not here (yet): vector/tensor histories, units-aware velocity + reduction, ALE / old-frame trace-back, forcing history, checkpoint state. + Use :class:`SemiLagrangian` for those. + + Parameters + ---------- + mesh, psi_fn, V_fn, degree, continuous, varsymbol, verbose, bcs, order, theta + As for :class:`SemiLagrangian`. ``psi_fn`` may be a scalar + ``MeshVariable`` (its nodal data is then copied into the snapshot + rather than re-evaluated) or a scalar expression. + v_degree : int, optional + Degree of the velocity snapshots (default: ``V_fn.degree`` if + ``V_fn`` is a mesh variable, else 2). + """ + + def __init__( + self, + mesh, + psi_fn, + V_fn, + vtype=VarType.SCALAR, + degree: int = 1, + continuous: bool = True, + varsymbol: Optional[str] = None, + verbose: bool = False, + bcs=[], + order: int = 1, + theta: float = 0.5, + v_degree: Optional[int] = None, + **_unsupported, + ): + super().__init__() + if vtype != VarType.SCALAR: + raise NotImplementedError( + "IntegrationPointSemiLagrangian: scalar histories only for now" + ) + self.mesh = mesh + self.bcs = bcs + self.verbose = verbose + self.degree = degree + self.continuous = continuous + self.order = order + self.theta = float(theta) + self.V_fn = V_fn + + if hasattr(psi_fn, "sym") and not isinstance(psi_fn, sympy.Basic): + self._psi_meshVar = psi_fn + self._psi_fn = psi_fn.sym + else: + self._psi_meshVar = None + self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) + self._v_meshVar = V_fn if (hasattr(V_fn, "sym") and not isinstance(V_fn, sympy.Basic)) else None + + self._init_history_tracking(order) + + if varsymbol is None: + varsymbol = rf"u_{{ [{self.instance_number}] }}" + inst = self.instance_number + + # History slots at the integration points (injected, never sampled). + self.psi_star = [ + uw.discretisation.IntegrationPointVariable( + f"psi_star_ip_{inst}_{k}", mesh, + varsymbol=rf"{{ {varsymbol}^{{ {'*' * (k + 1)} }} }}", + ) + for k in range(order) + ] + # Nodal snapshots of the solution and velocity at times n, n-1, ... + # (sampled at the departure points). + self.psi_snap = [ + uw.discretisation.MeshVariable( + f"psi_snap_ip_{inst}_{k}", mesh, 1, degree=degree, continuous=continuous, + varsymbol=rf"{{ {varsymbol}^{{ (n-{k}) }} }}", + ) + for k in range(order) + ] + if v_degree is None: + v_degree = getattr(V_fn, "degree", 2) + self.v_snap = [ + uw.discretisation.MeshVariable( + f"v_snap_ip_{inst}_{k}", mesh, mesh.dim, degree=v_degree, continuous=True, + varsymbol=rf"{{ V^{{ (n-{k}) }} }}", + ) + for k in range(order) + ] + self._init_coefficient_expressions(order, self.theta, with_exp=False) + + # ------------------------------------------------------------------ + @property + def psi_fn(self): + r"""Current symbolic expression :math:`\psi` being tracked.""" + return self._psi_fn + + @psi_fn.setter + def psi_fn(self, new_fn): + self._psi_meshVar = None + self._psi_fn = new_fn if isinstance(new_fn, sympy.Matrix) else sympy.Matrix([[new_fn]]) + + def _object_viewer(self): + from IPython.display import Latex, Markdown, display + super()._object_viewer() + display(Latex(r"$\quad\psi = $ " + self.psi_fn._repr_latex_())) + display(Latex(r"$\quad\mathbf{v} = $ " + sympy.Matrix(self.V_fn)._repr_latex_())) + display(Latex(rf"$\quad$History steps = {self.order} (at the integration points)")) + + # ------------------------------------------------------------------ + def _nudged_node_coords(self, var): + """ND node coordinates of ``var`` moved 0.1 % toward their cell + centroids so boundary nodes locate unambiguously (see + :meth:`SemiLagrangian._centroid_shifted_node_coords`).""" + coords = np.asarray(var.coords_nd) + cellid = self.mesh.get_closest_cells(coords).reshape(-1) + cent = np.asarray(self.mesh._centroids)[cellid] + return 0.999 * coords + 0.001 * cent + + def _record_current(self): + """Snapshot slot 0 <- the current solution and velocity.""" + ps = self.psi_snap[0] + if self._psi_meshVar is not None and ( + self._psi_meshVar.degree == ps.degree + and self._psi_meshVar.continuous == ps.continuous + ): + ps.data[...] = self._psi_meshVar.data[...] + else: + vals = uw.function.evaluate(self.psi_fn[0], self._nudged_node_coords(ps)) + ps.data[:, 0] = np.asarray(vals).reshape(-1) + vs = self.v_snap[0] + if self._v_meshVar is not None and self._v_meshVar.degree == vs.degree: + vs.data[...] = self._v_meshVar.data[...] + else: + vals = uw.function.evaluate(sympy.Matrix(self.V_fn), self._nudged_node_coords(vs)) + vs.data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) + + def _velocity_at(self, v_sym, coords, evalf): + v = uw.function.global_evaluate(v_sym, coords, evalf=evalf) + v = np.asarray(v) + if v.ndim == 3: + v = v[:, 0, :] + return v.reshape(coords.shape[0], self.mesh.dim) + + def _trace_segment(self, X, v_sym, dt, evalf): + r"""One RK2 (midpoint) segment of the characteristic, backwards: + ``x_mid = x - dt/2 v(x)``, ``x_dep = x - dt v(x_mid)``.""" + clamp = self.mesh.return_coords_to_bounds + v0 = self._velocity_at(v_sym, X, evalf) + Xm = X - 0.5 * dt * v0 + if clamp is not None: + Xm = clamp(Xm) + vm = self._velocity_at(v_sym, Xm, evalf) + Xd = X - dt * vm + if clamp is not None: + Xd = clamp(Xd) + return Xd + + def _segment_dt(self, j, dt): + """Length of segment ``j`` (0 = the current step).""" + if j == 0: + return dt + h = self._dt_history[j - 1] + return dt if h is None else h + + def _fill_slots(self, dt, evalf): + """Trace back from the integration points and sample the snapshots.""" + X0 = np.asarray(self.psi_star[0].coords_nd) + X = X0.copy() + for k in range(self.order): + # Segment k extends the trace from slot k-1's feet, so the + # feet for slot k are those of slot k-1 traced one more step. + X = self._trace_segment(X, self.v_snap[k].sym, self._segment_dt(k, dt), evalf) + vals = uw.function.global_evaluate(self.psi_snap[k].sym[0], X, evalf=evalf) + self.psi_star[k].data[:, 0] = np.asarray(vals).reshape(-1) + + def initialise_history(self): + """Start every snapshot and slot from the current field, so + ``bdf()`` is zero on the first step.""" + self._record_current() + for k in range(1, self.order): + self.psi_snap[k].data[...] = self.psi_snap[0].data[...] + self.v_snap[k].data[...] = self.v_snap[0].data[...] + X = np.asarray(self.psi_star[0].coords_nd) + vals = np.asarray(uw.function.evaluate(self.psi_snap[0].sym[0], X)).reshape(-1) + for k in range(self.order): + self.psi_star[k].data[:, 0] = vals + self._history_initialised = True + + def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): + self._dt = dt + if not self._history_initialised: + self.initialise_history() + _update_bdf_values(self._bdf_coeffs, self.effective_order, self._dt, self._dt_history) + _update_am_values(self._am_coeffs, self.effective_order, self.theta) + for k in range(self.order - 1, 0, -1): + self.psi_snap[k].data[...] = self.psi_snap[k - 1].data[...] + self.v_snap[k].data[...] = self.v_snap[k - 1].data[...] + self._record_current() + self._fill_slots(dt, evalf) + + def update(self, dt, evalf=False, verbose=False, **kwargs): + self.update_pre_solve(dt, evalf=evalf, verbose=verbose, **kwargs) + + def update_post_solve(self, dt, evalf=False, verbose=False, **_ignored): + self._dt = dt + for i in range(self.order - 1, 0, -1): + self._dt_history[i] = self._dt_history[i - 1] + self._dt_history[0] = dt + if self._n_solves_completed < self.order: + self._n_solves_completed += 1 diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py new file mode 100644 index 000000000..6977a9d2a --- /dev/null +++ b/tests/test_0066_integration_point_slcn.py @@ -0,0 +1,95 @@ +"""Semi-Lagrangian history at the integration points. + +Two properties: the value each slot carries is the snapshot evaluated +exactly at the traced departure point (the floor: for a P2 field and a +uniform velocity the sample is exact to round-off, for one and for two +segments), and on a rotating Gaussian the scheme is at least as accurate as +the nodal SLCN it replaces and keeps the peak better. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def test_slots_are_exact_departure_point_values(): + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) + f = lambda X: 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + X[:, 0] * X[:, 1] + T.data[:, 0] = f(np.asarray(T.coords)) + v = np.array([1.0, 0.5]) + V = sympy.Matrix([[v[0], v[1]]]) + dt = 0.1 + + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=2) + assert all(ps.is_integration_point for ps in ddt.psi_star) + + ddt.update_pre_solve(dt) + X = np.asarray(ddt.psi_star[0].coords) + foot = X - v * dt + inside = (foot > 0.0).all(1) & (foot < 1.0).all(1) + assert inside.sum() > 100 + assert np.abs(ddt.psi_star[0].data[inside, 0] - f(foot[inside])).max() < 1e-12 + + # Second slot: two segments back, sampled from the older snapshot. + ddt.update_post_solve(dt) + ddt.update_pre_solve(dt) + foot2 = X - v * 2 * dt + inside2 = (foot2 > 0.0).all(1) & (foot2 < 1.0).all(1) + assert np.abs(ddt.psi_star[1].data[inside2, 0] - f(foot2[inside2])).max() < 1e-12 + + # Negative control: the nodal scheme's slot is an interpolant of an + # interpolant, so the same check on it does not hold to round-off. + Tn = uw.discretisation.MeshVariable("Tn", mesh, 1, degree=2) + Tn.data[:, 0] = f(np.asarray(Tn.coords)) + nodal = uw.systems.ddt.SemiLagrangian(mesh, Tn, V, uw.VarType.SCALAR, degree=2, continuous=True, order=2) + nodal.update_pre_solve(dt) + nodal.update_post_solve(dt) + nodal.update_pre_solve(dt) + Xn = np.asarray(nodal.psi_star[1].coords) + footn = Xn - v * 2 * dt + insiden = (footn > 0.0).all(1) & (footn < 1.0).all(1) + assert np.abs(nodal.psi_star[1].data[insiden, 0] - f(footn[insiden])).max() > 1e-12 + + +def _rotating_gaussian(mesh, kind, dt, nsteps): + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + x0, sig = 0.5, 0.12 + gauss = lambda X, cx, cy: np.exp(-((X[:, 0] - cx) ** 2 + (X[:, 1] - cy) ** 2) / (2 * sig ** 2)) + T = uw.discretisation.MeshVariable(f"T_{kind}", mesh, 1, degree=2) + T.data[:, 0] = gauss(np.asarray(T.coords), x0, 0.0) + if kind == "ip": + DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=1) + adv = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=V, DuDt=DuDt, order=1) + else: + adv = uw.systems.AdvDiffusion(mesh, u_Field=T, V_fn=V, order=1) + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = 1e-9 + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + for _ in range(nsteps): + adv.solve(timestep=dt) + ang = nsteps * dt + exact = gauss(np.asarray(T.coords), x0 * np.cos(ang), x0 * np.sin(ang)) + E = uw.discretisation.MeshVariable(f"E_{kind}", mesh, 1, degree=2) + E.data[:, 0] = T.data[:, 0] - exact + l2 = np.sqrt(uw.maths.Integral(mesh, E.sym[0] ** 2).evaluate()) + return l2, T.data[:, 0].max() + + +@pytest.mark.level_2 +def test_rotating_gaussian_beats_nodal_slcn(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1, -1), maxCoords=(1, 1), cellSize=0.08, qdegree=2 + ) + dt, nsteps = 0.1, 16 + l2_nodal, peak_nodal = _rotating_gaussian(mesh, "nodal", dt, nsteps) + l2_ip, peak_ip = _rotating_gaussian(mesh, "ip", dt, nsteps) + assert l2_ip <= l2_nodal + assert peak_ip >= peak_nodal + assert l2_ip < 0.02 From 5e4a4342f690b5830106d1aa1a1ec814b8d8dc59 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 14:51:07 -0700 Subject: [PATCH 04/15] IntegrationPointSemiLagrangian: refuse an under-sampled rule; measurements The least-squares fit of the sampled departure-point values onto the continuous space is a per-cell interpolant through interior points when the rule has as many points as the element has local dofs (P2 on a triangle at qdegree 2), and a mode then grows ~1.1x per step at small Courant number (rotating Gaussian, C=0.25: flat for 75 steps, then blow-up; zero-velocity control stationary). At 2x the points (qdegree 3) the fit is contractive and the scheme is stable through a full revolution. The constructor now raises at <= 1x and warns below 2x. Adds a monotone_mode pass-through. Measured at cellSize 0.05, qdegree 3, half revolution, L2 nodal -> IP: C=0.25 1.30e-2 -> 6.3e-4; C=0.5 3.55e-3 -> 9.0e-4; C=1 3.62e-3 -> 3.24e-3; C=2 equal. Integral of T^2 over 63 steps at C=0.5: nodal -1.4%, IP -0.03%. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 55 +++++++++++++++++-- src/underworld3/systems/ddt.py | 43 ++++++++++++++- tests/test_0066_integration_point_slcn.py | 17 +++++- 3 files changed, 107 insertions(+), 8 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index bffb88a4d..a99d74951 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -136,8 +136,53 @@ no checkpoint state yet. For a P2 field in a uniform velocity the slots reproduce the exact departure-point values to round-off, for one and for two segments -(`tests/test_0066_integration_point_slcn.py`). On a rotating Gaussian at -Courant 1 (`cellSize=0.05`, `dt=0.1`, half a revolution) the L2 error drops -from 3.6e-3 (nodal SLCN) to 3.2e-3 and the peak is kept at 0.9998 instead -of 0.987. The trace-back samples six points per cell rather than the P2 -nodes, so the update costs about twice the nodal one. +(`tests/test_0066_integration_point_slcn.py`). + +### The rule must oversample the history space + +The solve fits the sampled departure-point values to the continuous space +by weighted least squares on the rule. With as many points per cell as the +element has local dofs (P2 on a triangle: 6 dofs, and 6 points at +`qdegree=2`) that fit is a per-cell interpolant through interior points, +which extrapolates, and at small Courant number a mode grows by about 1.1 +per step: on the rotating Gaussian below at Courant 0.25 the run was flat +for 75 steps and then blew up. At twice the points (`qdegree=3`, 12 on a +triangle) the fit is contractive and the scheme is stable through a full +revolution. The constructor therefore raises when the rule has no more +points than local dofs and warns below 2x. Raising `qdegree` costs every +solver on the mesh its assembly time, which is the price of this scheme. + +### Measured against nodal SLCN + +Rotating Gaussian (solid-body rotation, width 0.1 at radius 0.5), P2, unit +square of side 2, `cellSize=0.05`, `qdegree=3`, half a revolution, pure +advection. L2 error against the exact rotated field and the peak value: + +| Courant | nodal SLCN | integration-point | peak nodal / IP | +|---|---|---|---| +| 0.25 | 1.30e-2 | 6.3e-4 | 0.935 / 0.994 | +| 0.5 | 3.55e-3 | 9.0e-4 | 0.974 / 0.992 | +| 1 | 3.62e-3 | 3.24e-3 | 0.987 / 0.997 | +| 2 | 1.33e-2 | 1.33e-2 | 0.996 / 1.000 | + +The gain is largest at small Courant number, where the nodal scheme +re-interpolates most often per unit of transport; at Courant 2 the RK2 +trace-back error dominates and the two agree. Over a full revolution at +Courant 0.25 the error is 1.17e-3, twice the half-revolution value, so it +grows linearly. + +What the scheme conserves (Courant 0.5, 63 steps, relative change): + +| | integral of T | integral of T² | +|---|---|---| +| nodal SLCN | fluctuates within 3e-4 | -1.4 % (monotone) | +| integration-point | -3e-5 | -0.03 % | + +Neither scheme is exactly conservative (a Galerkin projection of a +transported field conserves the integral only with exact integration), but +the second moment is where the nodal scheme's diffusion shows and the +integration-point scheme loses 45 times less of it. + +The trace-back samples twelve points per cell rather than the P2 nodes, and +the snapshot evaluation at the moving feet misses the locator cache every +step, so the update costs about three times the nodal one. diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 42b6fa66f..427537904 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -3526,6 +3526,7 @@ def __init__( order: int = 1, theta: float = 0.5, v_degree: Optional[int] = None, + monotone_mode: Optional[str] = None, **_unsupported, ): super().__init__() @@ -3533,6 +3534,7 @@ def __init__( raise NotImplementedError( "IntegrationPointSemiLagrangian: scalar histories only for now" ) + self.monotone_mode = monotone_mode self.mesh = mesh self.bcs = bcs self.verbose = verbose @@ -3551,6 +3553,7 @@ def __init__( self._v_meshVar = V_fn if (hasattr(V_fn, "sym") and not isinstance(V_fn, sympy.Basic)) else None self._init_history_tracking(order) + self._check_rule_oversampling(degree) if varsymbol is None: varsymbol = rf"u_{{ [{self.instance_number}] }}" @@ -3584,6 +3587,42 @@ def __init__( ] self._init_coefficient_expressions(order, self.theta, with_exp=False) + def _check_rule_oversampling(self, degree): + """Refuse a rule with no more points per cell than the history space + has local dofs. + + The solve fits the sampled departure-point values to the continuous + space by weighted least squares on the rule. With as many points per + cell as the element has dofs (P2 on a triangle: 6 dofs, 6 points at + ``qdegree=2``) that fit is a per-cell interpolant through interior + points, which extrapolates, and a mode grows by ~1.1 per step at + small Courant number (measured: rotating Gaussian, C=0.25, blow-up + after ~80 steps). With twice the points (``qdegree=3``, 12 on a + triangle) the fit is contractive and the scheme is ~20x more + accurate than nodal SLCN. Below 2x we warn; at or below 1x we raise. + """ + PETSc.Options().setValue(f"ipsl_check_{self.instance_number}_petscspace_degree", degree) + fe = PETSc.FE().createDefault( + self.mesh.dim, 1, self.mesh.isSimplex, self.mesh.qdegree, + f"ipsl_check_{self.instance_number}_", PETSc.COMM_SELF, + ) + local_dofs = fe.getDimension() + Nq = len(np.asarray(self.mesh.integration_rule.getData()[1])) + if Nq <= local_dofs: + raise RuntimeError( + f"IntegrationPointSemiLagrangian: the mesh rule has {Nq} points per cell " + f"but a degree-{degree} history has {local_dofs} local dofs; the " + "least-squares fit is not oversampled and is unstable at small Courant " + f"number. Build the mesh with qdegree >= {self.mesh.qdegree + 1}." + ) + if Nq < 2 * local_dofs: + warnings.warn( + f"IntegrationPointSemiLagrangian: {Nq} rule points per cell for " + f"{local_dofs} local dofs is under 2x oversampling; stability at small " + "Courant number has only been verified at 2x (qdegree 3 for P2 on triangles).", + stacklevel=3, + ) + # ------------------------------------------------------------------ @property def psi_fn(self): @@ -3666,7 +3705,9 @@ def _fill_slots(self, dt, evalf): # Segment k extends the trace from slot k-1's feet, so the # feet for slot k are those of slot k-1 traced one more step. X = self._trace_segment(X, self.v_snap[k].sym, self._segment_dt(k, dt), evalf) - vals = uw.function.global_evaluate(self.psi_snap[k].sym[0], X, evalf=evalf) + vals = uw.function.global_evaluate( + self.psi_snap[k].sym[0], X, evalf=evalf, monotone=self.monotone_mode + ) self.psi_star[k].data[:, 0] = np.asarray(vals).reshape(-1) def initialise_history(self): diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 6977a9d2a..ac28f4aa4 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -17,7 +17,7 @@ def test_slots_are_exact_departure_point_values(): - mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=2) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=3) T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) f = lambda X: 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + X[:, 0] * X[:, 1] T.data[:, 0] = f(np.asarray(T.coords)) @@ -82,10 +82,23 @@ def _rotating_gaussian(mesh, kind, dt, nsteps): return l2, T.data[:, 0].max() +def test_undersampled_rule_is_refused(): + """P2 history on a qdegree-2 triangle mesh: 6 points for 6 local dofs. + That configuration blows up at small Courant number, so it is refused.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.2, qdegree=2) + T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) + V = sympy.Matrix([[1.0, 0.0]]) + with pytest.raises(RuntimeError, match="oversampled"): + uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2) + # P1 on the same rule is 2x oversampled and accepted. + T1 = uw.discretisation.MeshVariable("T1", mesh, 1, degree=1) + uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T1, V, degree=1) + + @pytest.mark.level_2 def test_rotating_gaussian_beats_nodal_slcn(): mesh = uw.meshing.UnstructuredSimplexBox( - minCoords=(-1, -1), maxCoords=(1, 1), cellSize=0.08, qdegree=2 + minCoords=(-1, -1), maxCoords=(1, 1), cellSize=0.08, qdegree=3 ) dt, nsteps = 0.1, 16 l2_nodal, peak_nodal = _rotating_gaussian(mesh, "nodal", dt, nsteps) From e46f49ac597cc3038feeb800c62d2ce8a16c2efe Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 15:18:52 -0700 Subject: [PATCH 05/15] IntegrationPointSemiLagrangian: stability measured by power iteration One-step growth factors (P2, cellSize 0.1, Courant 0.25). Pure advection: nodal SLCN 0.9989; integration-point 6 points 1.028, 9 points 1.005, 12 points 1.0003. With physical diffusion at cell Peclet 100: 6 points 1.0095, 9 and 12 points 0.996. The integration-point map is never strictly contractive under pure advection (it does not dissipate; the nodal scheme's 0.999 is its numerical diffusion) and oversampling brings it to neutral; with diffusion, 1.5x and 2x oversampling are stable and 1x is not. The conical 9-point rule ran a full revolution bounded (energy +0.07 %, saturating). Guard unchanged (raise at <= 1x, warn below 2x); the warning now states the measured behaviour. Docs carry the corrected mechanism: Galerkin with a sampled load is the weighted least-squares fit, the shifted field has sub-cell kinks, and the growth is aliasing of the sampled norm. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 40 +++++++++++++++++-- src/underworld3/systems/ddt.py | 23 ++++++----- 2 files changed, 50 insertions(+), 13 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index a99d74951..258a98a73 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -148,9 +148,43 @@ which extrapolates, and at small Courant number a mode grows by about 1.1 per step: on the rotating Gaussian below at Courant 0.25 the run was flat for 75 steps and then blew up. At twice the points (`qdegree=3`, 12 on a triangle) the fit is contractive and the scheme is stable through a full -revolution. The constructor therefore raises when the rule has no more -points than local dofs and warns below 2x. Raising `qdegree` costs every -solver on the mesh its assembly time, which is the price of this scheme. +revolution. At 1.5x (PETSc's conical rule, 9 points at degree 4, selected +with `-petscfe_default_quadrature_type conic`) it is also bounded +through a full revolution, with a 0.07 % rise in energy that saturates and +twice the L2 error of the 12-point rule. The constructor raises when the +rule has no more points than local dofs and warns below 2x. Raising the +rule costs every solver on the mesh its assembly time, which is the price +of this scheme; `qdegree` is the polynomial exactness of the rule, and the +extra exactness is incidental here, only the point count matters. + +Why a fit at all: with the load vector formed from point samples and the +mass matrix exact on the same rule, the Galerkin step is algebraically the +weighted least-squares fit `min Σ_q w_q (T(x_q) - g_q)²`. The composed +field `T^n ∘ X_dep` is piecewise P2 on the *shifted* mesh, so on the actual +cells it carries interior kinks wherever a cell's feet straddle a source +edge, and it is not in the space. The fit contracts in the sampled norm, +not in L2; a grid-scale mode shifted by a fraction of a cell can have a +sampled norm above its true norm (an aliasing error of the rule), and that +ratio is the growth per step. The nodal scheme is stable for a different +reason: interpolation at the nodes is bounded by the source's nodal values. + +Measured directly, by power iteration on the one-step operator (random +field renormalised every step, P2, `cellSize=0.1`, Courant 0.25): + +| growth per step | pure advection | cell Péclet 100 | +|---|---|---| +| nodal SLCN | 0.9989 | — | +| integration-point, 6 points (1x) | 1.028 | 1.0095 | +| integration-point, 9 points (1.5x) | 1.005 | 0.9960 | +| integration-point, 12 points (2x) | 1.0003 | 0.9962 | + +Under pure advection the integration-point map is never strictly +contractive; oversampling brings it toward neutral. That is the flip side +of not dissipating: the nodal scheme's 0.999 is its numerical diffusion. +With the physical diffusion a real problem carries (cell Péclet 100 here) +9 and 12 points are stable and 6 is not. The 1.5x case is therefore usable +with diffusion and slowly unstable without it (invisible over one +revolution, not over ten); 2x is neutral either way. ### Measured against nodal SLCN diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 427537904..d64dbf028 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -3592,14 +3592,16 @@ def _check_rule_oversampling(self, degree): has local dofs. The solve fits the sampled departure-point values to the continuous - space by weighted least squares on the rule. With as many points per - cell as the element has dofs (P2 on a triangle: 6 dofs, 6 points at - ``qdegree=2``) that fit is a per-cell interpolant through interior - points, which extrapolates, and a mode grows by ~1.1 per step at - small Courant number (measured: rotating Gaussian, C=0.25, blow-up - after ~80 steps). With twice the points (``qdegree=3``, 12 on a - triangle) the fit is contractive and the scheme is ~20x more - accurate than nodal SLCN. Below 2x we warn; at or below 1x we raise. + space by weighted least squares on the rule (the mass matrix is exact + on the rule, so Galerkin with a sampled load *is* that fit). The fit + contracts in the sampled norm only, and the shifted field's sampled + norm can exceed its true norm (aliasing of the rule on grid-scale + modes), so the pure-advection map is never strictly contractive. + Measured one-step growth factors, P2 on triangles, Courant 0.25 + (power iteration): 6 points 1.03-1.14 (blows up), 9 points 1.005, + 12 points 1.0003; with physical diffusion at cell Peclet 100: 6 + points 1.01 (still unstable), 9 and 12 points 0.996 (stable). Nodal + SLCN: 0.999. So: raise at <= 1x oversampling, warn below 2x. """ PETSc.Options().setValue(f"ipsl_check_{self.instance_number}_petscspace_degree", degree) fe = PETSc.FE().createDefault( @@ -3618,8 +3620,9 @@ def _check_rule_oversampling(self, degree): if Nq < 2 * local_dofs: warnings.warn( f"IntegrationPointSemiLagrangian: {Nq} rule points per cell for " - f"{local_dofs} local dofs is under 2x oversampling; stability at small " - "Courant number has only been verified at 2x (qdegree 3 for P2 on triangles).", + f"{local_dofs} local dofs is under 2x oversampling: weakly unstable under pure " + "advection (growth ~1.005/step at 1.5x, Courant 0.25) and stable with physical " + "diffusion at cell Peclet <= 100. 2x (qdegree 3 for P2 on triangles) is neutral.", stacklevel=3, ) From 78cbf6716cc7f2bf54d5bde0a59a5298ca75af6d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 18:47:40 -0700 Subject: [PATCH 06/15] Semi-Lagrangian trace-back: mid-point velocity at the mid time The RK2 trace is second order only if the mid-point velocity is the velocity at t^{n+1/2}. Both SemiLagrangian and IntegrationPointSemiLagrangian now take it there: on the current interval by extrapolation from the two most recent velocity fields, 1.5 v^n - 0.5 v^{n-1} (v^{n+1} is not known when the history is built), and on the older segments of a multi-step integration-point history by the average of the two known ends. With v^n alone the foot is off by b dt^2/2 in a flow accelerating at rate b. SemiLagrangian keeps the previous velocity in a managed v_prev mesh variable recorded after each trace (v^n alone on the first step); _velocity_nd_at takes an optional expression and now accepts a mesh variable as V_fn. The integration-point scheme keeps at least two velocity snapshots. tests/test_0066: uniformly accelerating uniform flow, both schemes hit the exact foot (integration-point to 1e-5, bounded by an evaluator edge case on one foot in ~2000; nodal to 1e-3, bounded by its 0.1 % centroid nudge), with the v^n-only foot as the failing control. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 16 ++++ src/underworld3/systems/ddt.py | 94 +++++++++++++++++-- tests/test_0066_integration_point_slcn.py | 47 ++++++++++ 3 files changed, 148 insertions(+), 9 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 258a98a73..13a588132 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -134,6 +134,22 @@ The diffusive flux history (`DFDt`) keeps its nodal projection, since it carries derivatives. Scalar histories only; no ALE or old-frame trace-back, no checkpoint state yet. +### The mid-point velocity is taken at the mid time + +The RK2 trace, `x_mid = x - dt/2 v(x)`, `x_dep = x - dt v(x_mid)`, is second +order only if `v(x_mid)` is the velocity at `t^{n+1/2}`. Both schemes now +take it there: on the current interval by extrapolation from the two most +recent velocity fields, `1.5 v^n - 0.5 v^{n-1}` (the velocity at `n+1` is +not known when the history is built), and on the older segments of a +multi-step history by the average of the two known ends. With `v^n` alone +the foot is off by `b dt²/2` in a flow accelerating at rate `b`, first +order in an unsteady flow; the rotating Gaussian did not show it because +that velocity is steady. `tests/test_0066_integration_point_slcn.py` +checks both schemes against the exact foot in a uniformly accelerating +flow, with the `v^n`-only foot as the control. `SemiLagrangian` keeps the +previous velocity in a `v_prev` mesh variable it manages; the +integration-point scheme keeps at least two velocity snapshots. + For a P2 field in a uniform velocity the slots reproduce the exact departure-point values to round-off, for one and for two segments (`tests/test_0066_integration_point_slcn.py`). diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index d64dbf028..3201ecf81 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -2256,12 +2256,59 @@ def _record_psi_star_from_field_data(self): except Exception: return None + def _midtime_velocity_expr(self): + r"""Velocity at :math:`t^{n+1/2}` for the mid-point stage of the + trace-back: :math:`\tfrac32 v^n - \tfrac12 v^{n-1}` once a previous + velocity has been recorded, else :math:`v^n`.""" + v_prev = getattr(self, "_v_prev", None) + if v_prev is None or not getattr(self, "_v_prev_valid", False): + return None + return self._V_matrix() * sympy.Rational(3, 2) - v_prev.sym * sympy.Rational(1, 2) + + def _V_matrix(self): + """``V_fn`` as a sympy row matrix (a mesh variable contributes its symbol).""" + V = self.V_fn + if hasattr(V, "sym") and not isinstance(V, sympy.Basic): + return sympy.Matrix(V.sym) + return sympy.Matrix(V) + + def _record_velocity_history(self): + """Store the current advecting velocity at the nodes as v^{n-1} + for the next step.""" + if getattr(self, "_v_prev", None) is None: + v_degree = getattr(self.V_fn, "degree", 2) + self._v_prev = uw.discretisation.MeshVariable( + f"v_prev_sl_{self.instance_number}", self.mesh, self.mesh.dim, + degree=v_degree, continuous=True, + varsymbol=rf"{{ V^{{ (n-1) }}_{{ [{self.instance_number}] }} }}", + ) + self._v_prev.remesh_policy = RemeshPolicy.CARRY + self._v_prev._remesh_managed_by = self + self._v_prev_valid = False + v_src = self.V_fn if not isinstance(self.V_fn, sympy.Basic) else None + if v_src is not None and getattr(v_src, "degree", None) == self._v_prev.degree: + self._v_prev.data[...] = v_src.data[...] + else: + coords = self._centroid_shifted_var_coords(self._v_prev) + vals = uw.function.evaluate(self._V_matrix(), coords) + self._v_prev.data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) + self._v_prev_valid = True + + def _centroid_shifted_var_coords(self, var): + """ND node coordinates of ``var`` nudged 0.1 % toward their cell + centroids (see :meth:`_centroid_shifted_node_coords`).""" + coords = np.asarray(var.coords_nd) + cellid = self.mesh.get_closest_cells(coords).reshape(-1) + cent = np.asarray(self.mesh._centroids)[cellid] + return 0.999 * coords + 0.001 * cent + def _velocity_nd_at( self, coords, use_global: bool = False, evalf: bool = False, subtract_v_mesh: bool = False, + expr=None, ): r"""Evaluate the advecting velocity at ``coords``, reduced to ND space. @@ -2293,15 +2340,16 @@ def _velocity_nd_at( (rather than symbolically as ``V_fn − v_mesh.sym``) so the subtraction inherits the same unit treatment as ``V_fn``. """ + fn = self._V_matrix() if expr is None else expr if use_global: - v_result = uw.function.global_evaluate(self.V_fn, coords, evalf=evalf) + v_result = uw.function.global_evaluate(fn, coords, evalf=evalf) if subtract_v_mesh: v_mesh = uw.function.global_evaluate( self._v_mesh_var.sym, coords, evalf=evalf ) v_result = v_result - v_mesh else: - v_result = uw.function.evaluate(self.V_fn, coords) + v_result = uw.function.evaluate(fn, coords) if subtract_v_mesh: v_mesh = uw.function.evaluate(self._v_mesh_var.sym, coords) v_result = v_result - v_mesh @@ -2562,12 +2610,17 @@ def _trace_departure_points( # Mid-point velocities may lie off-rank, so route through # global_evaluate (with evalf forwarded), unlike the on-node - # evaluation above. + # evaluation above. The mid-point velocity is taken at the mid + # TIME, t^{n+1/2}, by extrapolation from the two most recent + # velocity fields, 1.5 v^n - 0.5 v^{n-1}; with v^n alone the + # trace is only first order in an unsteady flow. On the first + # step (no previous velocity) v^n is used. v_at_mid_pts = self._velocity_nd_at( mid_pt_coords, use_global=True, evalf=evalf, subtract_v_mesh=subtract_v_mesh, + expr=self._midtime_velocity_expr(), ) # Upstream (departure) coordinates: current position - velocity * timestep @@ -2779,6 +2832,10 @@ def update_pre_solve( _oldframe_active, _oldframe_X, ) + # The velocity used this step becomes v^{n-1} for the next + # step's mid-time extrapolation. + self._record_velocity_history() + # Phase-2 ALE: consume the one-step v_mesh pulse. Subsequent # non-adapt steps will see no pending displacement and run a # plain trace-back. If multiple adapts happen before the next @@ -3578,12 +3635,15 @@ def __init__( ] if v_degree is None: v_degree = getattr(V_fn, "degree", 2) + # At least two velocity levels: the current interval's mid-time + # velocity is extrapolated from v^n and v^{n-1}. + self._n_v = max(order, 2) self.v_snap = [ uw.discretisation.MeshVariable( f"v_snap_ip_{inst}_{k}", mesh, mesh.dim, degree=v_degree, continuous=True, varsymbol=rf"{{ V^{{ (n-{k}) }} }}", ) - for k in range(order) + for k in range(self._n_v) ] self._init_coefficient_expressions(order, self.theta, with_exp=False) @@ -3679,15 +3739,16 @@ def _velocity_at(self, v_sym, coords, evalf): v = v[:, 0, :] return v.reshape(coords.shape[0], self.mesh.dim) - def _trace_segment(self, X, v_sym, dt, evalf): + def _trace_segment(self, X, v_start_sym, v_mid_sym, dt, evalf): r"""One RK2 (midpoint) segment of the characteristic, backwards: - ``x_mid = x - dt/2 v(x)``, ``x_dep = x - dt v(x_mid)``.""" + ``x_mid = x - dt/2 v_start(x)``, ``x_dep = x - dt v_mid(x_mid)``, + with ``v_mid`` the velocity at the segment's mid TIME.""" clamp = self.mesh.return_coords_to_bounds - v0 = self._velocity_at(v_sym, X, evalf) + v0 = self._velocity_at(v_start_sym, X, evalf) Xm = X - 0.5 * dt * v0 if clamp is not None: Xm = clamp(Xm) - vm = self._velocity_at(v_sym, Xm, evalf) + vm = self._velocity_at(v_mid_sym, Xm, evalf) Xd = X - dt * vm if clamp is not None: Xd = clamp(Xd) @@ -3704,10 +3765,23 @@ def _fill_slots(self, dt, evalf): """Trace back from the integration points and sample the snapshots.""" X0 = np.asarray(self.psi_star[0].coords_nd) X = X0.copy() + half = sympy.Rational(1, 2) for k in range(self.order): # Segment k extends the trace from slot k-1's feet, so the # feet for slot k are those of slot k-1 traced one more step. - X = self._trace_segment(X, self.v_snap[k].sym, self._segment_dt(k, dt), evalf) + # Segment k runs from t^{n+1-k} back to t^{n-k}. Its mid-time + # velocity: for k=0 extrapolated, 1.5 v^n - 0.5 v^{n-1} (v^{n+1} + # is not known yet); for k>=1 both ends are known, so the + # average of v^{n+1-k} and v^{n-k}. The first stage, which + # only places the mid-point, uses the velocity at the + # segment's start time. + if k == 0: + v_start = self.v_snap[0].sym + v_mid = self.v_snap[0].sym * sympy.Rational(3, 2) - self.v_snap[1].sym * half + else: + v_start = self.v_snap[k - 1].sym + v_mid = (self.v_snap[k - 1].sym + self.v_snap[k].sym) * half + X = self._trace_segment(X, v_start, v_mid, self._segment_dt(k, dt), evalf) vals = uw.function.global_evaluate( self.psi_snap[k].sym[0], X, evalf=evalf, monotone=self.monotone_mode ) @@ -3719,6 +3793,7 @@ def initialise_history(self): self._record_current() for k in range(1, self.order): self.psi_snap[k].data[...] = self.psi_snap[0].data[...] + for k in range(1, self._n_v): self.v_snap[k].data[...] = self.v_snap[0].data[...] X = np.asarray(self.psi_star[0].coords_nd) vals = np.asarray(uw.function.evaluate(self.psi_snap[0].sym[0], X)).reshape(-1) @@ -3734,6 +3809,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): _update_am_values(self._am_coeffs, self.effective_order, self.theta) for k in range(self.order - 1, 0, -1): self.psi_snap[k].data[...] = self.psi_snap[k - 1].data[...] + for k in range(self._n_v - 1, 0, -1): self.v_snap[k].data[...] = self.v_snap[k - 1].data[...] self._record_current() self._fill_slots(dt, evalf) diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index ac28f4aa4..86bce667d 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -106,3 +106,50 @@ def test_rotating_gaussian_beats_nodal_slcn(): assert l2_ip <= l2_nodal assert peak_ip >= peak_nodal assert l2_ip < 0.02 + + +def _unsteady_uniform_flow_check(kind): + """Uniform velocity that changes linearly in time, v(t) = a + b t. The + exact foot for the interval [t1, t1 + dt] is x - dt (a + b (t1 + dt/2)). + With the mid-time velocity extrapolated from v(t1) and v(t0) the trace + reproduces it; with v(t1) alone the foot is off by b dt^2 / 2.""" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.1, qdegree=3) + T = uw.discretisation.MeshVariable(f"T_{kind}", mesh, 1, degree=2) + v_var = uw.discretisation.MeshVariable(f"v_{kind}", mesh, mesh.dim, degree=2) + f = lambda X: 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + X[:, 0] * X[:, 1] + T.data[:, 0] = f(np.asarray(T.coords)) + a, b, dt = np.array([1.0, 0.5]), np.array([2.0, -1.0]), 0.1 + + if kind == "ip": + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, v_var, degree=2, order=1) + else: + ddt = uw.systems.ddt.SemiLagrangian(mesh, T, v_var, uw.VarType.SCALAR, degree=2, continuous=True, order=1) + + v_var.data[...] = a + b * 0.0 + ddt.update_pre_solve(dt) + ddt.update_post_solve(dt) + v_var.data[...] = a + b * dt + ddt.update_pre_solve(dt) + + X = np.asarray(ddt.psi_star[0].coords) + exact_foot = X - dt * (a + b * (dt + 0.5 * dt)) + naive_foot = X - dt * (a + b * dt) + inside = (exact_foot > 0.02).all(1) & (exact_foot < 0.98).all(1) & (naive_foot > 0.02).all(1) & (naive_foot < 0.98).all(1) + assert inside.sum() > 100 + got = np.asarray(ddt.psi_star[0].data[:, 0]) + return np.abs(got[inside] - f(exact_foot[inside])).max(), np.abs(got[inside] - f(naive_foot[inside])).max() + + +@pytest.mark.parametrize("kind", ["ip", "nodal"]) +def test_midtime_velocity_makes_the_trace_second_order(kind): + err_exact, err_naive = _unsteady_uniform_flow_check(kind) + # The trace is exact to round-off. Integration-point: the tolerance + # covers the evaluator, which returns one foot in ~2000 a few 1e-6 off + # (a locator edge case shared by evaluate and global_evaluate). Nodal: + # the scheme traces from nodes nudged 0.1 % toward the cell centroid + # and stores at the unnudged node, an error of order 0.001 h |grad psi| + # per step (2e-4 here) that the integration-point scheme does not have. + assert err_exact < (1e-5 if kind == "ip" else 1e-3) + # Negative control: the foot from v^n alone is b dt^2/2 away, which for + # this quadratic field is a visible difference. + assert err_naive > 1e-3 From c903007cc5592448c31fec16024281848aa8c044 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 20:48:07 -0700 Subject: [PATCH 07/15] SemiLagrangian: midtime_velocity switch (False = the pre-2026-09 v^n-only trace) For controls and reproduction of earlier runs. Verified a strict no-op for a steady velocity (4e-17) and, on Blankenbach 1a, bit-for-bit reproduction of the published nodal SLCN run when off. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- src/underworld3/systems/ddt.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 3201ecf81..8f065f5b1 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -1545,6 +1545,7 @@ def __init__( monotone_mode: Optional[str] = None, theta: float = 0.5, old_frame_traceback: bool = False, + midtime_velocity: bool = True, ): super().__init__() @@ -1556,6 +1557,9 @@ def __init__( self._psi_fn = psi_fn self.V_fn = V_fn self.order = order + # Mid-point velocity of the RK2 trace at the mid TIME (1.5 v^n - + # 0.5 v^{n-1}); False reproduces the pre-2026-09 v^n-only trace. + self.midtime_velocity = bool(midtime_velocity) if preserve_moments: raise NotImplementedError( "preserve_moments is not currently implemented" @@ -2260,6 +2264,8 @@ def _midtime_velocity_expr(self): r"""Velocity at :math:`t^{n+1/2}` for the mid-point stage of the trace-back: :math:`\tfrac32 v^n - \tfrac12 v^{n-1}` once a previous velocity has been recorded, else :math:`v^n`.""" + if not getattr(self, "midtime_velocity", True): + return None v_prev = getattr(self, "_v_prev", None) if v_prev is None or not getattr(self, "_v_prev_valid", False): return None @@ -2834,7 +2840,8 @@ def update_pre_solve( # The velocity used this step becomes v^{n-1} for the next # step's mid-time extrapolation. - self._record_velocity_history() + if getattr(self, "midtime_velocity", True): + self._record_velocity_history() # Phase-2 ALE: consume the one-step v_mesh pulse. Subsequent # non-adapt steps will see no pending displacement and run a From cfb1d96705dd50ab68e2acae19ff61bcc6a86223 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 20:59:29 -0700 Subject: [PATCH 08/15] SemiLagrangian: copy the velocity variable's data into v_prev when V_fn is its symbol The solvers pass V_fn as v.sym, so the previous velocity was evaluated at nodes nudged 0.1 % toward their cell centroids, a bias of 0.001 h |grad v| that the mid-time extrapolation fed into every trace. On Blankenbach 1a it moved the wall Nusselt number by 0.9 % (a 5e-4 first-cell temperature change; interior transport and Vrms unchanged), and restarts from either state relaxed to distinct wall values within 50 steps. Resolving the symbol to its mesh variable (meshVariable_lookup_by_symbol) and copying restores the exact no-op for a steady velocity (4e-17) and the published nodal wall value. Diagnostics that cleared the alternatives: recording without using the extrapolation, and adding an unrelated field after the first solve, both left the wall value unchanged. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- src/underworld3/systems/ddt.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 8f065f5b1..fad5643f7 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -2291,8 +2291,12 @@ def _record_velocity_history(self): self._v_prev.remesh_policy = RemeshPolicy.CARRY self._v_prev._remesh_managed_by = self self._v_prev_valid = False - v_src = self.V_fn if not isinstance(self.V_fn, sympy.Basic) else None + v_src = self._velocity_mesh_variable() if v_src is not None and getattr(v_src, "degree", None) == self._v_prev.degree: + # Exact copy. Evaluating at (nudged) nodes instead leaves a + # 0.001 h |grad v| bias in v_prev that the extrapolation feeds + # into every trace: on Blankenbach 1a it moved the wall Nusselt + # number by 0.9 % (first-cell temperature, 5e-4). self._v_prev.data[...] = v_src.data[...] else: coords = self._centroid_shifted_var_coords(self._v_prev) @@ -2300,6 +2304,20 @@ def _record_velocity_history(self): self._v_prev.data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) self._v_prev_valid = True + def _velocity_mesh_variable(self): + """The mesh variable behind ``V_fn`` if ``V_fn`` is one, or exactly + one's symbol; else None.""" + V = self.V_fn + if hasattr(V, "sym") and not isinstance(V, sympy.Basic): + return V + try: + found = uw.discretisation.meshVariable_lookup_by_symbol(self.mesh, sympy.Matrix(V)) + except Exception: + found = None + if found is not None and found[1] == -1: + return found[0] + return None + def _centroid_shifted_var_coords(self, var): """ND node coordinates of ``var`` nudged 0.1 % toward their cell centroids (see :meth:`_centroid_shifted_node_coords`).""" From b9e11257660a6764ea5b61969a52f8a921d1c9c2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 6 Sep 2026 22:48:47 -0700 Subject: [PATCH 09/15] Semi-Lagrangian velocity history: snapshot the variables inside V_fn V_fn is symbolic by design (-v, v/2, v - v_mesh must just work), so the previous velocity is no longer a separate field evaluated from it, nor is V_fn resolved to a mesh variable. Both schemes now find the mesh variables V_fn contains, keep a copy of each per history level, and form v^{n-1} as V_fn with those variables substituted by their copies: exact for any expression, an analytic V_fn reduces to itself. Shared helpers on _DDtBase (_make_velocity_level, _copy_velocity_level). The steady no-op holds for -v/2 (2e-19); tests cover V_fn as the variable, -v and v/2 for both schemes. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 14 +- src/underworld3/systems/ddt.py | 151 +++++++++--------- tests/test_0066_integration_point_slcn.py | 23 +-- 3 files changed, 101 insertions(+), 87 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 13a588132..0e61ecf15 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -146,9 +146,17 @@ the foot is off by `b dt²/2` in a flow accelerating at rate `b`, first order in an unsteady flow; the rotating Gaussian did not show it because that velocity is steady. `tests/test_0066_integration_point_slcn.py` checks both schemes against the exact foot in a uniformly accelerating -flow, with the `v^n`-only foot as the control. `SemiLagrangian` keeps the -previous velocity in a `v_prev` mesh variable it manages; the -integration-point scheme keeps at least two velocity snapshots. +flow, for `V_fn` given as the variable, as `-v` and as `v/2`, with the +`v^n`-only foot as the control. + +`V_fn` stays symbolic by design (`-v`, `v/2`, `v - v_mesh` must all just +work), so the previous velocity is never a separate field evaluated from +it. Each scheme snapshots the mesh variables that `V_fn` contains and forms +`v^{n-1}` as `V_fn` with those variables substituted by their copies: exact +for any expression, and an analytic `V_fn` reduces to itself. Evaluating +`V_fn` at the (nudged) nodes instead left a `0.001 h |grad v|` bias that the +extrapolation fed into every trace; on Blankenbach 1a it moved the wall +Nusselt number by 0.9 %. For a P2 field in a uniform velocity the slots reproduce the exact departure-point values to round-off, for one and for two segments diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index fad5643f7..f93ab937d 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -675,6 +675,55 @@ def _history_syms(self): """ return [ps.sym for ps in self.psi_star] + # ----- velocity-expression snapshots (mid-time trace-back) ----- + def _V_matrix(self): + """``V_fn`` as a sympy row matrix (a mesh variable contributes its symbol).""" + V = self.V_fn + if hasattr(V, "sym") and not isinstance(V, sympy.Basic): + return sympy.Matrix(V.sym) + return sympy.Matrix(V) + + def _velocity_variables(self): + """The distinct mesh variables ``V_fn`` is built from (may be empty: + an analytic velocity).""" + _, varfns, _ = uw.function.expressions.mesh_vars_in_expression(self._V_matrix()) + seen, out = set(), [] + for fn in varfns: + var = fn.meshvar() + if id(var) not in seen: + seen.add(id(var)) + out.append(var) + return out + + def _make_velocity_level(self, tag): + """A snapshot level: one copy of every mesh variable in ``V_fn`` plus + ``V_fn`` with those variables substituted by their copies. ``V_fn`` + stays whatever expression the user gave (``-v``, ``v/2``, ``v - v_mesh``): + the substitution keeps it exact at any earlier time.""" + copies, subs = {}, {} + for k, var in enumerate(self._velocity_variables()): + snap = uw.discretisation.MeshVariable( + f"vsnap_{tag}_{self.instance_number}_{k}", self.mesh, var.num_components, + vtype=var.vtype, degree=var.degree, continuous=var.continuous, + varsymbol=rf"{{ {var.symbol}^{{ ({tag}) }} }}", + ) + snap.remesh_policy = RemeshPolicy.CARRY + snap._remesh_managed_by = self + copies[var.clean_name] = (var, snap) + for a, b in zip(var.sym_1d, snap.sym_1d): + subs[a] = b + expr = self._V_matrix().applyfunc(lambda e: e.xreplace(subs)) if subs else self._V_matrix() + return {"copies": copies, "expr": expr} + + @staticmethod + def _copy_velocity_level(dst, src=None): + """``dst`` <- ``src`` (another level) or, with ``src=None``, the live variables.""" + for name, (var, snap) in dst["copies"].items(): + if src is None: + snap.data[...] = var.data[...] + else: + snap.data[...] = src["copies"][name][1].data[...] + def bdf(self, order: Optional[int] = None): r"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -2263,61 +2312,28 @@ def _record_psi_star_from_field_data(self): def _midtime_velocity_expr(self): r"""Velocity at :math:`t^{n+1/2}` for the mid-point stage of the trace-back: :math:`\tfrac32 v^n - \tfrac12 v^{n-1}` once a previous - velocity has been recorded, else :math:`v^n`.""" + velocity has been recorded, else :math:`v^n`. :math:`v^{n-1}` is + ``V_fn`` with the mesh variables it contains replaced by their + snapshots, so any expression (``-v``, ``v/2``, ``v - v_mesh``) is + carried exactly; an analytic ``V_fn`` reduces to itself.""" if not getattr(self, "midtime_velocity", True): return None - v_prev = getattr(self, "_v_prev", None) - if v_prev is None or not getattr(self, "_v_prev_valid", False): + level = getattr(self, "_v_prev_level", None) + if level is None or not getattr(self, "_v_prev_valid", False): return None - return self._V_matrix() * sympy.Rational(3, 2) - v_prev.sym * sympy.Rational(1, 2) - - def _V_matrix(self): - """``V_fn`` as a sympy row matrix (a mesh variable contributes its symbol).""" - V = self.V_fn - if hasattr(V, "sym") and not isinstance(V, sympy.Basic): - return sympy.Matrix(V.sym) - return sympy.Matrix(V) + return self._V_matrix() * sympy.Rational(3, 2) - level["expr"] * sympy.Rational(1, 2) def _record_velocity_history(self): - """Store the current advecting velocity at the nodes as v^{n-1} - for the next step.""" - if getattr(self, "_v_prev", None) is None: - v_degree = getattr(self.V_fn, "degree", 2) - self._v_prev = uw.discretisation.MeshVariable( - f"v_prev_sl_{self.instance_number}", self.mesh, self.mesh.dim, - degree=v_degree, continuous=True, - varsymbol=rf"{{ V^{{ (n-1) }}_{{ [{self.instance_number}] }} }}", - ) - self._v_prev.remesh_policy = RemeshPolicy.CARRY - self._v_prev._remesh_managed_by = self + """Snapshot the mesh variables inside ``V_fn`` as v^{n-1} for the + next step. A copy, never an evaluation: evaluating at (nudged) nodes + left a 0.001 h |grad v| bias that the extrapolation fed into every + trace and moved the Blankenbach 1a wall Nusselt number by 0.9 %.""" + if getattr(self, "_v_prev_level", None) is None: + self._v_prev_level = self._make_velocity_level("n-1") self._v_prev_valid = False - v_src = self._velocity_mesh_variable() - if v_src is not None and getattr(v_src, "degree", None) == self._v_prev.degree: - # Exact copy. Evaluating at (nudged) nodes instead leaves a - # 0.001 h |grad v| bias in v_prev that the extrapolation feeds - # into every trace: on Blankenbach 1a it moved the wall Nusselt - # number by 0.9 % (first-cell temperature, 5e-4). - self._v_prev.data[...] = v_src.data[...] - else: - coords = self._centroid_shifted_var_coords(self._v_prev) - vals = uw.function.evaluate(self._V_matrix(), coords) - self._v_prev.data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) + self._copy_velocity_level(self._v_prev_level) self._v_prev_valid = True - def _velocity_mesh_variable(self): - """The mesh variable behind ``V_fn`` if ``V_fn`` is one, or exactly - one's symbol; else None.""" - V = self.V_fn - if hasattr(V, "sym") and not isinstance(V, sympy.Basic): - return V - try: - found = uw.discretisation.meshVariable_lookup_by_symbol(self.mesh, sympy.Matrix(V)) - except Exception: - found = None - if found is not None and found[1] == -1: - return found[0] - return None - def _centroid_shifted_var_coords(self, var): """ND node coordinates of ``var`` nudged 0.1 % toward their cell centroids (see :meth:`_centroid_shifted_node_coords`).""" @@ -3589,9 +3605,8 @@ class IntegrationPointSemiLagrangian(_DDtBase): As for :class:`SemiLagrangian`. ``psi_fn`` may be a scalar ``MeshVariable`` (its nodal data is then copied into the snapshot rather than re-evaluated) or a scalar expression. - v_degree : int, optional - Degree of the velocity snapshots (default: ``V_fn.degree`` if - ``V_fn`` is a mesh variable, else 2). + ``V_fn`` may be any expression of mesh variables (``-v``, ``v/2``); the + velocity history snapshots the variables it contains. """ def __init__( @@ -3607,7 +3622,6 @@ def __init__( bcs=[], order: int = 1, theta: float = 0.5, - v_degree: Optional[int] = None, monotone_mode: Optional[str] = None, **_unsupported, ): @@ -3632,7 +3646,6 @@ def __init__( else: self._psi_meshVar = None self._psi_fn = psi_fn if isinstance(psi_fn, sympy.Matrix) else sympy.Matrix([[psi_fn]]) - self._v_meshVar = V_fn if (hasattr(V_fn, "sym") and not isinstance(V_fn, sympy.Basic)) else None self._init_history_tracking(order) self._check_rule_oversampling(degree) @@ -3658,18 +3671,12 @@ def __init__( ) for k in range(order) ] - if v_degree is None: - v_degree = getattr(V_fn, "degree", 2) # At least two velocity levels: the current interval's mid-time - # velocity is extrapolated from v^n and v^{n-1}. + # velocity is extrapolated from v^n and v^{n-1}. Each level is a + # snapshot of the mesh variables inside V_fn, substituted into the + # expression, so V_fn may be any expression of them. self._n_v = max(order, 2) - self.v_snap = [ - uw.discretisation.MeshVariable( - f"v_snap_ip_{inst}_{k}", mesh, mesh.dim, degree=v_degree, continuous=True, - varsymbol=rf"{{ V^{{ (n-{k}) }} }}", - ) - for k in range(self._n_v) - ] + self.v_levels = [self._make_velocity_level(f"n-{k}") for k in range(self._n_v)] self._init_coefficient_expressions(order, self.theta, with_exp=False) def _check_rule_oversampling(self, degree): @@ -3750,12 +3757,7 @@ def _record_current(self): else: vals = uw.function.evaluate(self.psi_fn[0], self._nudged_node_coords(ps)) ps.data[:, 0] = np.asarray(vals).reshape(-1) - vs = self.v_snap[0] - if self._v_meshVar is not None and self._v_meshVar.degree == vs.degree: - vs.data[...] = self._v_meshVar.data[...] - else: - vals = uw.function.evaluate(sympy.Matrix(self.V_fn), self._nudged_node_coords(vs)) - vs.data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) + self._copy_velocity_level(self.v_levels[0]) def _velocity_at(self, v_sym, coords, evalf): v = uw.function.global_evaluate(v_sym, coords, evalf=evalf) @@ -3800,12 +3802,13 @@ def _fill_slots(self, dt, evalf): # average of v^{n+1-k} and v^{n-k}. The first stage, which # only places the mid-point, uses the velocity at the # segment's start time. + V = [lvl["expr"] for lvl in self.v_levels] if k == 0: - v_start = self.v_snap[0].sym - v_mid = self.v_snap[0].sym * sympy.Rational(3, 2) - self.v_snap[1].sym * half + v_start = V[0] + v_mid = V[0] * sympy.Rational(3, 2) - V[1] * half else: - v_start = self.v_snap[k - 1].sym - v_mid = (self.v_snap[k - 1].sym + self.v_snap[k].sym) * half + v_start = V[k - 1] + v_mid = (V[k - 1] + V[k]) * half X = self._trace_segment(X, v_start, v_mid, self._segment_dt(k, dt), evalf) vals = uw.function.global_evaluate( self.psi_snap[k].sym[0], X, evalf=evalf, monotone=self.monotone_mode @@ -3819,7 +3822,7 @@ def initialise_history(self): for k in range(1, self.order): self.psi_snap[k].data[...] = self.psi_snap[0].data[...] for k in range(1, self._n_v): - self.v_snap[k].data[...] = self.v_snap[0].data[...] + self._copy_velocity_level(self.v_levels[k], self.v_levels[0]) X = np.asarray(self.psi_star[0].coords_nd) vals = np.asarray(uw.function.evaluate(self.psi_snap[0].sym[0], X)).reshape(-1) for k in range(self.order): @@ -3835,7 +3838,7 @@ def update_pre_solve(self, dt, evalf=False, verbose=False, **_ignored): for k in range(self.order - 1, 0, -1): self.psi_snap[k].data[...] = self.psi_snap[k - 1].data[...] for k in range(self._n_v - 1, 0, -1): - self.v_snap[k].data[...] = self.v_snap[k - 1].data[...] + self._copy_velocity_level(self.v_levels[k], self.v_levels[k - 1]) self._record_current() self._fill_slots(dt, evalf) diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 86bce667d..feb599d0a 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -108,7 +108,7 @@ def test_rotating_gaussian_beats_nodal_slcn(): assert l2_ip < 0.02 -def _unsteady_uniform_flow_check(kind): +def _unsteady_uniform_flow_check(kind, vform="var"): """Uniform velocity that changes linearly in time, v(t) = a + b t. The exact foot for the interval [t1, t1 + dt] is x - dt (a + b (t1 + dt/2)). With the mid-time velocity extrapolated from v(t1) and v(t0) the trace @@ -119,11 +119,13 @@ def _unsteady_uniform_flow_check(kind): f = lambda X: 1.0 + 2.0 * X[:, 0] - 3.0 * X[:, 1] + 0.5 * X[:, 0] ** 2 + X[:, 0] * X[:, 1] T.data[:, 0] = f(np.asarray(T.coords)) a, b, dt = np.array([1.0, 0.5]), np.array([2.0, -1.0]), 0.1 + # V_fn is symbolic by design: any expression of the variable must work. + V_fn, factor = {"var": (v_var, 1.0), "neg": (-v_var.sym, -1.0), "half": (v_var.sym / 2, 0.5)}[vform] if kind == "ip": - ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, v_var, degree=2, order=1) + ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V_fn, degree=2, order=1) else: - ddt = uw.systems.ddt.SemiLagrangian(mesh, T, v_var, uw.VarType.SCALAR, degree=2, continuous=True, order=1) + ddt = uw.systems.ddt.SemiLagrangian(mesh, T, V_fn, uw.VarType.SCALAR, degree=2, continuous=True, order=1) v_var.data[...] = a + b * 0.0 ddt.update_pre_solve(dt) @@ -132,24 +134,25 @@ def _unsteady_uniform_flow_check(kind): ddt.update_pre_solve(dt) X = np.asarray(ddt.psi_star[0].coords) - exact_foot = X - dt * (a + b * (dt + 0.5 * dt)) - naive_foot = X - dt * (a + b * dt) + exact_foot = X - dt * factor * (a + b * (dt + 0.5 * dt)) + naive_foot = X - dt * factor * (a + b * dt) inside = (exact_foot > 0.02).all(1) & (exact_foot < 0.98).all(1) & (naive_foot > 0.02).all(1) & (naive_foot < 0.98).all(1) assert inside.sum() > 100 got = np.asarray(ddt.psi_star[0].data[:, 0]) return np.abs(got[inside] - f(exact_foot[inside])).max(), np.abs(got[inside] - f(naive_foot[inside])).max() +@pytest.mark.parametrize("vform", ["var", "neg", "half"]) @pytest.mark.parametrize("kind", ["ip", "nodal"]) -def test_midtime_velocity_makes_the_trace_second_order(kind): - err_exact, err_naive = _unsteady_uniform_flow_check(kind) +def test_midtime_velocity_makes_the_trace_second_order(kind, vform): + err_exact, err_naive = _unsteady_uniform_flow_check(kind, vform) # The trace is exact to round-off. Integration-point: the tolerance - # covers the evaluator, which returns one foot in ~2000 a few 1e-6 off - # (a locator edge case shared by evaluate and global_evaluate). Nodal: + # covers the evaluator, which returns one foot in ~2000 up to a few 1e-5 + # off (a locator edge case shared by evaluate and global_evaluate). Nodal: # the scheme traces from nodes nudged 0.1 % toward the cell centroid # and stores at the unnudged node, an error of order 0.001 h |grad psi| # per step (2e-4 here) that the integration-point scheme does not have. - assert err_exact < (1e-5 if kind == "ip" else 1e-3) + assert err_exact < (1e-4 if kind == "ip" else 1e-3) # Negative control: the foot from v^n alone is b dt^2/2 away, which for # this quadratic field is a visible difference. assert err_naive > 1e-3 From e1d689400b8f1fecbf57a84a977af65f0731d37a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 7 Sep 2026 07:52:13 -0700 Subject: [PATCH 10/15] Semi-Lagrangian velocity history: cache by evaluation at the true nodes Substituting snapshots of the mesh variables into V_fn misses everything else the expression depends on: a constant that ramps, a swarm proxy, a mesh that has moved. The previous velocity is now V_fn evaluated at the true node coordinates of a vector field (highest degree of the variables in V_fn, 2 for an analytic velocity) and cached per history level, in both schemes. No nudge: the evaluator is exact at node coordinates on simplex, quad and annulus meshes (2e-16 for an expression of a variable, a constant and the coordinates), which is what the earlier nudged evaluation lacked. Tests: V_fn as the variable, -v, v/2, and c*v with c changed between the two steps (the case substitution gets wrong), for both schemes; steady no-op 4e-16. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 22 +++-- src/underworld3/systems/ddt.py | 89 +++++++++---------- tests/test_0066_integration_point_slcn.py | 25 ++++-- 3 files changed, 74 insertions(+), 62 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 0e61ecf15..68f0b27bd 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -149,14 +149,20 @@ checks both schemes against the exact foot in a uniformly accelerating flow, for `V_fn` given as the variable, as `-v` and as `v/2`, with the `v^n`-only foot as the control. -`V_fn` stays symbolic by design (`-v`, `v/2`, `v - v_mesh` must all just -work), so the previous velocity is never a separate field evaluated from -it. Each scheme snapshots the mesh variables that `V_fn` contains and forms -`v^{n-1}` as `V_fn` with those variables substituted by their copies: exact -for any expression, and an analytic `V_fn` reduces to itself. Evaluating -`V_fn` at the (nudged) nodes instead left a `0.001 h |grad v|` bias that the -extrapolation fed into every trace; on Blankenbach 1a it moved the wall -Nusselt number by 0.9 %. +`V_fn` stays symbolic by design (`-v`, `v/2`, `c(t) v`, `v - v_mesh` must +all just work), and the previous velocity is **cached by evaluation**: +`V_fn` evaluated at the true nodes of a vector field of the highest degree +among the variables it contains. That captures everything the expression +depends on as it was at that time: the variables, constants that ramp, +swarm proxies, the mesh geometry. Substituting snapshots of the mesh +variables into the expression would read a ramping constant at its +current value. Evaluating at nodes *nudged* into their cells, the old +boundary-node workaround, left a `0.001 h |grad v|` bias that the +extrapolation fed into every trace and moved the Blankenbach 1a wall +Nusselt number by 0.9 %; the evaluator is exact at the true node +coordinates on simplex, quad and annulus meshes (2e-16), so no nudge is +used here. The test covers `V_fn` as the variable, `-v`, `v/2` and `c v` +with `c` changed between steps. For a P2 field in a uniform velocity the slots reproduce the exact departure-point values to round-off, for one and for two segments diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index f93ab937d..e51505f99 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -683,46 +683,37 @@ def _V_matrix(self): return sympy.Matrix(V.sym) return sympy.Matrix(V) - def _velocity_variables(self): - """The distinct mesh variables ``V_fn`` is built from (may be empty: - an analytic velocity).""" + def _velocity_degree(self): + """Degree of the nodal velocity cache: the highest degree among the + mesh variables in ``V_fn`` (2 for an analytic velocity).""" _, varfns, _ = uw.function.expressions.mesh_vars_in_expression(self._V_matrix()) - seen, out = set(), [] - for fn in varfns: - var = fn.meshvar() - if id(var) not in seen: - seen.add(id(var)) - out.append(var) - return out + degs = [fn.meshvar().degree for fn in varfns] + return max(degs) if degs else 2 def _make_velocity_level(self, tag): - """A snapshot level: one copy of every mesh variable in ``V_fn`` plus - ``V_fn`` with those variables substituted by their copies. ``V_fn`` - stays whatever expression the user gave (``-v``, ``v/2``, ``v - v_mesh``): - the substitution keeps it exact at any earlier time.""" - copies, subs = {}, {} - for k, var in enumerate(self._velocity_variables()): - snap = uw.discretisation.MeshVariable( - f"vsnap_{tag}_{self.instance_number}_{k}", self.mesh, var.num_components, - vtype=var.vtype, degree=var.degree, continuous=var.continuous, - varsymbol=rf"{{ {var.symbol}^{{ ({tag}) }} }}", - ) - snap.remesh_policy = RemeshPolicy.CARRY - snap._remesh_managed_by = self - copies[var.clean_name] = (var, snap) - for a, b in zip(var.sym_1d, snap.sym_1d): - subs[a] = b - expr = self._V_matrix().applyfunc(lambda e: e.xreplace(subs)) if subs else self._V_matrix() - return {"copies": copies, "expr": expr} - - @staticmethod - def _copy_velocity_level(dst, src=None): - """``dst`` <- ``src`` (another level) or, with ``src=None``, the live variables.""" - for name, (var, snap) in dst["copies"].items(): - if src is None: - snap.data[...] = var.data[...] - else: - snap.data[...] = src["copies"][name][1].data[...] + """A cached velocity level: ``V_fn`` EVALUATED at the true nodes of a + vector field (no nudge; the evaluator is exact at node coordinates on + simplex, quad and annulus meshes). Caching by evaluation, rather than + by substituting snapshots of the mesh variables into the expression, + is what captures everything ``V_fn`` depends on at that time: the + variables, constants that ramp, swarm proxies, the mesh geometry.""" + snap = uw.discretisation.MeshVariable( + f"vcache_{tag}_{self.instance_number}", self.mesh, self.mesh.dim, + degree=self._velocity_degree(), continuous=True, + varsymbol=rf"{{ V^{{ ({tag}) }}_{{ [{self.instance_number}] }} }}", + ) + snap.remesh_policy = RemeshPolicy.CARRY + snap._remesh_managed_by = self + return {"var": snap, "expr": snap.sym} + + def _copy_velocity_level(self, dst, src=None): + """``dst`` <- ``src`` (another level) or, with ``src=None``, ``V_fn`` + evaluated now at ``dst``'s nodes.""" + if src is None: + vals = uw.function.evaluate(self._V_matrix(), np.asarray(dst["var"].coords_nd)) + dst["var"].data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) + else: + dst["var"].data[...] = src["var"].data[...] def bdf(self, order: Optional[int] = None): r"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -2313,9 +2304,9 @@ def _midtime_velocity_expr(self): r"""Velocity at :math:`t^{n+1/2}` for the mid-point stage of the trace-back: :math:`\tfrac32 v^n - \tfrac12 v^{n-1}` once a previous velocity has been recorded, else :math:`v^n`. :math:`v^{n-1}` is - ``V_fn`` with the mesh variables it contains replaced by their - snapshots, so any expression (``-v``, ``v/2``, ``v - v_mesh``) is - carried exactly; an analytic ``V_fn`` reduces to itself.""" + ``V_fn`` as evaluated at the previous step and cached at the nodes, + so any expression (``-v``, ``v/2``, ``c(t) v``, ``v - v_mesh``) is + carried as it was then.""" if not getattr(self, "midtime_velocity", True): return None level = getattr(self, "_v_prev_level", None) @@ -2324,10 +2315,10 @@ def _midtime_velocity_expr(self): return self._V_matrix() * sympy.Rational(3, 2) - level["expr"] * sympy.Rational(1, 2) def _record_velocity_history(self): - """Snapshot the mesh variables inside ``V_fn`` as v^{n-1} for the - next step. A copy, never an evaluation: evaluating at (nudged) nodes - left a 0.001 h |grad v| bias that the extrapolation fed into every - trace and moved the Blankenbach 1a wall Nusselt number by 0.9 %.""" + """Cache ``V_fn`` evaluated at the true nodes as v^{n-1} for the + next step. Evaluating at NUDGED nodes left a 0.001 h |grad v| bias + that the extrapolation fed into every trace and moved the + Blankenbach 1a wall Nusselt number by 0.9 %.""" if getattr(self, "_v_prev_level", None) is None: self._v_prev_level = self._make_velocity_level("n-1") self._v_prev_valid = False @@ -3605,8 +3596,8 @@ class IntegrationPointSemiLagrangian(_DDtBase): As for :class:`SemiLagrangian`. ``psi_fn`` may be a scalar ``MeshVariable`` (its nodal data is then copied into the snapshot rather than re-evaluated) or a scalar expression. - ``V_fn`` may be any expression of mesh variables (``-v``, ``v/2``); the - velocity history snapshots the variables it contains. + ``V_fn`` may be any expression (``-v``, ``v/2``, ``c(t) v``); the + velocity history caches it by evaluation at each time level. """ def __init__( @@ -3672,9 +3663,9 @@ def __init__( for k in range(order) ] # At least two velocity levels: the current interval's mid-time - # velocity is extrapolated from v^n and v^{n-1}. Each level is a - # snapshot of the mesh variables inside V_fn, substituted into the - # expression, so V_fn may be any expression of them. + # velocity is extrapolated from v^n and v^{n-1}. Each level caches + # V_fn evaluated at the nodes at that time, so V_fn may be any + # expression (variables, ramping constants, swarm proxies). self._n_v = max(order, 2) self.v_levels = [self._make_velocity_level(f"n-{k}") for k in range(self._n_v)] self._init_coefficient_expressions(order, self.theta, with_exp=False) diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index feb599d0a..67a73efa1 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -120,7 +120,12 @@ def _unsteady_uniform_flow_check(kind, vform="var"): T.data[:, 0] = f(np.asarray(T.coords)) a, b, dt = np.array([1.0, 0.5]), np.array([2.0, -1.0]), 0.1 # V_fn is symbolic by design: any expression of the variable must work. - V_fn, factor = {"var": (v_var, 1.0), "neg": (-v_var.sym, -1.0), "half": (v_var.sym / 2, 0.5)}[vform] + # "ramp": a constant that changes between the two steps; the cached + # previous velocity must carry the OLD value (substituting snapshots of + # the variables into the expression would read the new one). + c = uw.expression(r"c_{ramp}", 1.0, "ramping factor") + V_fn, factor = {"var": (v_var, 1.0), "neg": (-v_var.sym, -1.0), "half": (v_var.sym / 2, 0.5), + "ramp": (c * v_var.sym, 1.0)}[vform] if kind == "ip": ddt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V_fn, degree=2, order=1) @@ -130,19 +135,29 @@ def _unsteady_uniform_flow_check(kind, vform="var"): v_var.data[...] = a + b * 0.0 ddt.update_pre_solve(dt) ddt.update_post_solve(dt) - v_var.data[...] = a + b * dt + if vform == "ramp": + # v(t1) = 1.5 (a + b dt) through the constant; v(t0) = a. Extrapolated + # mid-time velocity 1.5 v(t1) - 0.5 v(t0) = 2.25 (a + b dt) - 0.5 a. + c.sym = 1.5 + v_var.data[...] = a + b * dt + v_mid = 2.25 * (a + b * dt) - 0.5 * a + v_naive = 1.5 * (a + b * dt) + else: + v_var.data[...] = a + b * dt + v_mid = factor * (a + b * (dt + 0.5 * dt)) + v_naive = factor * (a + b * dt) ddt.update_pre_solve(dt) X = np.asarray(ddt.psi_star[0].coords) - exact_foot = X - dt * factor * (a + b * (dt + 0.5 * dt)) - naive_foot = X - dt * factor * (a + b * dt) + exact_foot = X - dt * v_mid + naive_foot = X - dt * v_naive inside = (exact_foot > 0.02).all(1) & (exact_foot < 0.98).all(1) & (naive_foot > 0.02).all(1) & (naive_foot < 0.98).all(1) assert inside.sum() > 100 got = np.asarray(ddt.psi_star[0].data[:, 0]) return np.abs(got[inside] - f(exact_foot[inside])).max(), np.abs(got[inside] - f(naive_foot[inside])).max() -@pytest.mark.parametrize("vform", ["var", "neg", "half"]) +@pytest.mark.parametrize("vform", ["var", "neg", "half", "ramp"]) @pytest.mark.parametrize("kind", ["ip", "nodal"]) def test_midtime_velocity_makes_the_trace_second_order(kind, vform): err_exact, err_naive = _unsteady_uniform_flow_check(kind, vform) From c4562bad47f3a0e487314b25e13c5f9a01948cbb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 7 Sep 2026 11:08:05 -0700 Subject: [PATCH 11/15] Tests: make the integration-point tests rank-aware Two-rank run had one rank fail and the other block in the next collective. The evaluate test now keeps only the points the rank owns and compares to round-off (the evaluator adds an ulp on the way out); the rotating-Gaussian peak is a global maximum; the hand-quadrature test is marked serial-only. Two ranks: 33 passed, 1 skipped; serial unchanged. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 3 ++- tests/test_0065_integration_point_variable.py | 15 ++++++++++----- tests/test_0066_integration_point_slcn.py | 4 +++- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 68f0b27bd..b99f523ed 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -61,7 +61,8 @@ agrees with what the assembler used at that point; an interpolant or a projection would report a different viscosity from the one the solver saw. Points no cell owns take the nearest integration point on the rank. -Evaluating the variable at its own `coords` returns its own `data` exactly. +Evaluating the variable at its own `coords` returns its own `data` to round-off +(the point selection is exact; the evaluator pipeline can add an ulp). If a smooth nodal picture is wanted (a plot, a diagnostic), project the symbol onto a `MeshVariable` explicitly with `SNES_Projection`; the diff --git a/tests/test_0065_integration_point_variable.py b/tests/test_0065_integration_point_variable.py index 1742ae8f7..4f3344363 100644 --- a/tests/test_0065_integration_point_variable.py +++ b/tests/test_0065_integration_point_variable.py @@ -60,6 +60,7 @@ def test_layout_and_geometry(kind): assert np.allclose(cent, np.asarray(mesh._centroids)[:n], atol=1e-12) +@pytest.mark.skipif(uw.mpi.size > 1, reason="the hand quadrature sum is over the rank's local cells; serial only") def test_assembler_reads_the_stored_values(): """Integral of random point data == the quadrature sum done by hand.""" mesh = _mesh("triangle") @@ -111,18 +112,22 @@ def test_evaluate_is_nearest_point_of_owning_cell(): rng = np.random.default_rng(2) h.data[:, 0] = rng.uniform(size=h.data.shape[0]) - # Exact at its own points. + # Exact at its own points (the index selection is exact; the evaluator + # pipeline can add an ulp of round-off on the way out). own = uw.function.evaluate(h.sym[0], np.asarray(h.coords)).reshape(-1) - assert np.array_equal(own, h.data[:, 0]) + assert np.allclose(own, h.data[:, 0], rtol=0, atol=1e-14) - # Nearest point of the owning cell elsewhere. + # Nearest point of the owning cell elsewhere. In parallel keep only the + # points this rank owns (the locator returns -1 for the others). pts = rng.uniform(0.05, 0.95, size=(300, 2)) - cells = mesh._robust_owning_cells(pts) + cells = np.asarray(mesh._robust_owning_cells(pts)).reshape(-1) + pts, cells = pts[cells >= 0], cells[cells >= 0] + assert len(pts) > 20 ipc = h.integration_points j = ((ipc[cells] - pts[:, None, :]) ** 2).sum(-1).argmin(1) expected = h.cell_data[cells, j, 0] got = uw.function.evaluate(h.sym[0], pts).reshape(-1) - assert np.array_equal(got, expected) + assert np.allclose(got, expected, rtol=0, atol=1e-14) # Negative control: a nodal P1 interpolant of the same data does not # match this definition (it smooths), so the test discriminates. diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 67a73efa1..19a19f109 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -79,7 +79,9 @@ def _rotating_gaussian(mesh, kind, dt, nsteps): E = uw.discretisation.MeshVariable(f"E_{kind}", mesh, 1, degree=2) E.data[:, 0] = T.data[:, 0] - exact l2 = np.sqrt(uw.maths.Integral(mesh, E.sym[0] ** 2).evaluate()) - return l2, T.data[:, 0].max() + from mpi4py import MPI + peak = uw.mpi.comm.allreduce(float(T.data[:, 0].max()), op=MPI.MAX) # global, not rank-local + return l2, peak def test_undersampled_rule_is_refused(): From 5c1549a493e7407ab5c0e25882153c991d3fe66a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 7 Sep 2026 11:59:18 -0700 Subject: [PATCH 12/15] Semi-Lagrangian velocity cache and integration-point history: units frame CI (test_1056) caught the units-active SLCN diverging 63 % from the non-dimensional run: the velocity cache was created without units and filled with dimensional values from evaluate, while .data is the non-dimensional store, so the mid-time expression mixed frames. The cache variable now carries V_fn's units and every evaluated value is reduced with _to_nondim_ndarray before storing (the nodal history's own idiom, issue #267); the same for the integration-point history's snapshots and slots and its trace velocities. test_1056 now covers both schemes. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- src/underworld3/systems/ddt.py | 29 ++++++++++++++++++++----- tests/test_1056_units_slcn_traceback.py | 23 +++++++++++++------- 2 files changed, 39 insertions(+), 13 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e51505f99..de570f022 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -701,16 +701,26 @@ def _make_velocity_level(self, tag): f"vcache_{tag}_{self.instance_number}", self.mesh, self.mesh.dim, degree=self._velocity_degree(), continuous=True, varsymbol=rf"{{ V^{{ ({tag}) }}_{{ [{self.instance_number}] }} }}", + units=self._velocity_units(), # same frame as V_fn, so 1.5 v - 0.5 v_prev is consistent ) snap.remesh_policy = RemeshPolicy.CARRY snap._remesh_managed_by = self return {"var": snap, "expr": snap.sym} + def _velocity_units(self): + """Units of ``V_fn`` under an active units model, else None.""" + units = uw.get_units(self._V_matrix()) + if units is not None and not uw.get_default_model().has_units(): + units = None + return units + def _copy_velocity_level(self, dst, src=None): """``dst`` <- ``src`` (another level) or, with ``src=None``, ``V_fn`` - evaluated now at ``dst``'s nodes.""" + evaluated now at ``dst``'s nodes (reduced to the non-dimensional + frame ``.data`` stores, issue #267).""" if src is None: vals = uw.function.evaluate(self._V_matrix(), np.asarray(dst["var"].coords_nd)) + vals = _to_nondim_ndarray(vals, units=self._velocity_units()) dst["var"].data[...] = np.asarray(vals).reshape(-1, self.mesh.dim) else: dst["var"].data[...] = src["var"].data[...] @@ -3645,11 +3655,17 @@ def __init__( varsymbol = rf"u_{{ [{self.instance_number}] }}" inst = self.instance_number + psi_units = uw.get_units(self._psi_fn) + if psi_units is not None and not uw.get_default_model().has_units(): + psi_units = None + self._psi_units = psi_units + # History slots at the integration points (injected, never sampled). self.psi_star = [ uw.discretisation.IntegrationPointVariable( f"psi_star_ip_{inst}_{k}", mesh, varsymbol=rf"{{ {varsymbol}^{{ {'*' * (k + 1)} }} }}", + units=psi_units, ) for k in range(order) ] @@ -3659,6 +3675,7 @@ def __init__( uw.discretisation.MeshVariable( f"psi_snap_ip_{inst}_{k}", mesh, 1, degree=degree, continuous=continuous, varsymbol=rf"{{ {varsymbol}^{{ (n-{k}) }} }}", + units=psi_units, ) for k in range(order) ] @@ -3747,12 +3764,12 @@ def _record_current(self): ps.data[...] = self._psi_meshVar.data[...] else: vals = uw.function.evaluate(self.psi_fn[0], self._nudged_node_coords(ps)) - ps.data[:, 0] = np.asarray(vals).reshape(-1) + ps.data[:, 0] = np.asarray(_to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1) self._copy_velocity_level(self.v_levels[0]) def _velocity_at(self, v_sym, coords, evalf): v = uw.function.global_evaluate(v_sym, coords, evalf=evalf) - v = np.asarray(v) + v = np.asarray(_to_nondim_ndarray(v, units=self._velocity_units())) if v.ndim == 3: v = v[:, 0, :] return v.reshape(coords.shape[0], self.mesh.dim) @@ -3804,7 +3821,8 @@ def _fill_slots(self, dt, evalf): vals = uw.function.global_evaluate( self.psi_snap[k].sym[0], X, evalf=evalf, monotone=self.monotone_mode ) - self.psi_star[k].data[:, 0] = np.asarray(vals).reshape(-1) + self.psi_star[k].data[:, 0] = np.asarray( + _to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1) def initialise_history(self): """Start every snapshot and slot from the current field, so @@ -3815,7 +3833,8 @@ def initialise_history(self): for k in range(1, self._n_v): self._copy_velocity_level(self.v_levels[k], self.v_levels[0]) X = np.asarray(self.psi_star[0].coords_nd) - vals = np.asarray(uw.function.evaluate(self.psi_snap[0].sym[0], X)).reshape(-1) + vals = uw.function.evaluate(self.psi_snap[0].sym[0], X) + vals = np.asarray(_to_nondim_ndarray(vals, units=self._psi_units)).reshape(-1) for k in range(self.order): self.psi_star[k].data[:, 0] = vals self._history_initialised = True diff --git a/tests/test_1056_units_slcn_traceback.py b/tests/test_1056_units_slcn_traceback.py index de2e29a2b..c8ae585d3 100644 --- a/tests/test_1056_units_slcn_traceback.py +++ b/tests/test_1056_units_slcn_traceback.py @@ -23,7 +23,7 @@ pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] -def _advect_blob(use_units, vy=20.0, nsteps=5, dt=2.0): +def _advect_blob(use_units, vy=20.0, nsteps=5, dt=2.0, scheme="slcn"): uw.reset_default_model() model = uw.get_default_model() if use_units: @@ -34,13 +34,14 @@ def _advect_blob(use_units, vy=20.0, nsteps=5, dt=2.0): temperature_difference=uw.quantity(1000, "K"), ) mesh = uw.meshing.StructuredQuadBox( - elementRes=(12, 12), minCoords=(0.0, 0.0), maxCoords=(1000.0, 1000.0), units="km" + elementRes=(12, 12), minCoords=(0.0, 0.0), maxCoords=(1000.0, 1000.0), units="km", + qdegree=3, ) T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2, units="K") V = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2, units="m/s") else: mesh = uw.meshing.StructuredQuadBox( - elementRes=(12, 12), minCoords=(0.0, 0.0), maxCoords=(1000.0, 1000.0) + elementRes=(12, 12), minCoords=(0.0, 0.0), maxCoords=(1000.0, 1000.0), qdegree=3, ) T = uw.discretisation.MeshVariable("T", mesh, 1, degree=2) V = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2) @@ -51,7 +52,11 @@ def _advect_blob(use_units, vy=20.0, nsteps=5, dt=2.0): c = T.coords_nd # DM-space coords (identical units vs nondim) T.data[:, 0] = np.exp(-(((c[:, 0] - 500) / 120) ** 2 + ((c[:, 1] - 300) / 120) ** 2)) - adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym) + if scheme == "slcn_ip": + DuDt = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V.sym, degree=2, order=1) + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym, DuDt=DuDt, order=1) + else: + adv = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym) adv.constitutive_model = uw.constitutive_models.DiffusionModel adv.constitutive_model.Parameters.diffusivity = 1.0e-6 adv.add_dirichlet_bc([0.0], "Bottom") @@ -69,13 +74,15 @@ def test_units_slcn_traceback_runs(): assert Tu.max() < 1.05 and Tu.min() > -0.05 -def test_units_slcn_matches_nondimensional(): +@pytest.mark.parametrize("scheme", ["slcn", "slcn_ip"]) +def test_units_slcn_matches_nondimensional(scheme): """A units-active advection must track the equivalent non-dimensional run. They share identical ND values, so the trace-back (done in ND space) gives the same transport. A small residual (~1e-3) remains from the constitutive - diffusivity scaling under units — a separate concern from the trace-back.""" - Tu = _advect_blob(use_units=True) - Tn = _advect_blob(use_units=False) + diffusivity scaling under units — a separate concern from the trace-back. + Covers the nodal scheme and the integration-point history.""" + Tu = _advect_blob(use_units=True, scheme=scheme) + Tn = _advect_blob(use_units=False, scheme=scheme) rel = np.linalg.norm(Tu - Tn) / np.linalg.norm(Tn) assert rel < 5.0e-3, f"units-active SLCN diverges from nondimensional: rel L2 = {rel:.3e}" From 602a599e00e3f46eb60de7d26ec49f8080822ea7 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 7 Sep 2026 23:33:21 -0700 Subject: [PATCH 13/15] Docs: the integration-point history is an AdvDiffusionSLCN plugin The composed AdvDiffusion (#688) applies its theta-rule diffusive flux to the history level, which differentiates the slot; the JIT guard refuses a gradient of a delta field. The SLCN solver's separate nodal DFDt is the structure this history needs. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- docs/developer/subsystems/integration-point-variables.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 43113199a..0c73e6bdf 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -135,6 +135,13 @@ The diffusive flux history (`DFDt`) keeps its nodal projection, since it carries derivatives. Scalar histories only; no ALE or old-frame trace-back, no checkpoint state yet. +Use it with `AdvDiffusionSLCN`. The composed `uw.systems.AdvDiffusion` +(#688) applies its theta-rule diffusive flux to the history level as well as +the new one, which differentiates the history slot; the JIT guard refuses +that, correctly, because a delta field has no gradient. The SLCN solver +carries the diffusive history in a separate nodal `DFDt`, which is the +structure this history needs. + ### The mid-point velocity is taken at the mid time The RK2 trace, `x_mid = x - dt/2 v(x)`, `x_dep = x - dt v(x_mid)`, is second From 7d2083c9cd93b5410d4f879b7baea995e2460b46 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 8 Sep 2026 08:48:31 -0700 Subject: [PATCH 14/15] Integration-point history reachable in the composed AdvDiffusion at BDF2 and theta = 1 The composed solver (#688) takes the history as its transport manager. At order 2 and at theta = 1 no spatial term sits on the old level, so the integration-point history runs there and matches the SLCN solver's field to 3e-3 on a rotating Gaussian. theta = 1 needed one change: the old-level Adams-Moulton weight is identically zero there but was a runtime constant, so 0 * grad(psi*) was still differentiated at code generation and the guard fired on a dead term; spatial_weights now returns a literal zero for theta = 1 (for every manager). The Crank-Nicolson flux differentiates the old level and stays with AdvDiffusionSLCN; the guard's refusal is tested. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- .../subsystems/integration-point-variables.md | 14 ++++--- src/underworld3/systems/ddt.py | 12 ++++++ tests/test_0066_integration_point_slcn.py | 41 +++++++++++++++++++ 3 files changed, 61 insertions(+), 6 deletions(-) diff --git a/docs/developer/subsystems/integration-point-variables.md b/docs/developer/subsystems/integration-point-variables.md index 0c73e6bdf..542ec9f76 100644 --- a/docs/developer/subsystems/integration-point-variables.md +++ b/docs/developer/subsystems/integration-point-variables.md @@ -135,12 +135,14 @@ The diffusive flux history (`DFDt`) keeps its nodal projection, since it carries derivatives. Scalar histories only; no ALE or old-frame trace-back, no checkpoint state yet. -Use it with `AdvDiffusionSLCN`. The composed `uw.systems.AdvDiffusion` -(#688) applies its theta-rule diffusive flux to the history level as well as -the new one, which differentiates the history slot; the JIT guard refuses -that, correctly, because a delta field has no gradient. The SLCN solver -carries the diffusive history in a separate nodal `DFDt`, which is the -structure this history needs. +It is the transport manager of either advection-diffusion solver. In the +composed `uw.systems.AdvDiffusion` (#688) it runs at `order=2` (BDF2) and +at `order=1, theta=1`, where no spatial term sits on the old level; on a +rotating Gaussian the field matches the SLCN solver's to 3e-3. With the +Crank-Nicolson flux (`theta=0.5`) the composed solver differentiates the +old level, which a delta field cannot supply, and the JIT guard refuses +with a clear message; for that scheme use `AdvDiffusionSLCN`, whose +diffusive history is a separate nodal `DFDt`. ### The mid-point velocity is taken at the mid time diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 55a8df021..8b92d17a4 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -4093,6 +4093,18 @@ def __init__( self.v_levels = [self._make_velocity_level(f"n-{k}") for k in range(self._n_v)] self._init_coefficient_expressions(order, self.theta, with_exp=False) + def spatial_weights(self): + """As the base class, except that at ``theta = 1`` the old-level + weights, identically zero, are returned as literals. A runtime + constant with value zero would leave ``0 * grad(psi*)`` in the weak + form, and the slot has no gradient to differentiate (the JIT guard + would refuse a dead term). This is what makes the history usable in + the composed ``AdvDiffusion`` at ``order=1, theta=1``.""" + w = super().spatial_weights() + if self.integrator == "am" and float(self.theta) == 1.0: + return [sympy.Integer(1)] + [sympy.Integer(0)] * (len(w) - 1) + return w + def _check_rule_oversampling(self, degree): """Refuse a rule with no more points per cell than the history space has local dofs. diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 707d3fcee..18520e449 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -173,3 +173,44 @@ def test_midtime_velocity_makes_the_trace_second_order(kind, vform): # Negative control: the foot from v^n alone is b dt^2/2 away, which for # this quadratic field is a visible difference. assert err_naive > 1e-3 + + +@pytest.mark.parametrize("config", ["order2", "theta1", "cn"]) +def test_composed_advdiffusion_reachability(config): + """The composed uw.systems.AdvDiffusion (#688) takes the history as its + transport manager. With no spatial term on the old level (BDF2, or + theta = 1) the integration-point history runs there and matches the SLCN + solver; with the Crank-Nicolson flux (theta = 0.5) the old level is + differentiated, which a delta field cannot supply, and the JIT guard + refuses with a clear message.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1, -1), maxCoords=(1, 1), cellSize=0.1, qdegree=3 + ) + x, y = mesh.X + V = sympy.Matrix([[-y, x]]) + gauss = lambda X: np.exp(-((X[:, 0] - 0.5) ** 2 + X[:, 1] ** 2) / (2 * 0.12 ** 2)) + order, theta = {"order2": (2, 1.0), "theta1": (1, 1.0), "cn": (1, 0.5)}[config] + + def run(solver_cls, kwargs): + T = uw.discretisation.MeshVariable(f"T_{config}_{solver_cls.__name__}", mesh, 1, degree=2) + T.data[:, 0] = gauss(np.asarray(T.coords)) + D = uw.systems.ddt.IntegrationPointSemiLagrangian(mesh, T, V, degree=2, order=order, theta=theta) + adv = solver_cls(mesh, u_Field=T, V_fn=V, DuDt=D, order=order, **kwargs) + adv.constitutive_model = uw.constitutive_models.DiffusionModel + adv.constitutive_model.Parameters.diffusivity = 1e-9 + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + for _ in range(8): + adv.solve(timestep=0.1) + return np.asarray(T.data[:, 0]).copy() + + if config == "cn": + with pytest.raises(RuntimeError, match="integration-point"): + run(uw.systems.AdvDiffusion, {}) + return + kw = {"theta": theta} if order == 1 else {} + T_composed = run(uw.systems.AdvDiffusion, kw) + T_slcn = run(uw.systems.AdvDiffusionSLCN, {}) + # Same history, same time derivative; the solvers differ only in how the + # (negligible) diffusion is applied, so the fields agree closely. + assert np.abs(T_composed - T_slcn).max() < 5e-3 From cfffcaca580219c46ac8567c751341df85d19ab6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 8 Sep 2026 09:14:22 -0700 Subject: [PATCH 15/15] IntegrationPointSemiLagrangian: per-instance bcs default; oversampling message computes the needed qdegree Copilot review on #703: bcs=[] was a shared mutable default (now None -> a fresh list per instance), and the under-sampling error suggested qdegree + 1 regardless of the history degree. The message now scans PETSc's default rules for this cell type and reports the smallest qdegree with more points than local dofs and the one reaching 2x (P2 on triangles: 3 and 3; P3: 3 and 5). Underworld development team with AI support from Claude Code Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01MSGAFeA7qYXgkuw9ud8F2G --- src/underworld3/systems/ddt.py | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 8b92d17a4..add33178b 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -4026,7 +4026,7 @@ def __init__( continuous: bool = True, varsymbol: Optional[str] = None, verbose: bool = False, - bcs=[], + bcs=None, order: int = 1, theta: float = 0.5, monotone_mode: Optional[str] = None, @@ -4039,7 +4039,7 @@ def __init__( ) self.monotone_mode = monotone_mode self.mesh = mesh - self.bcs = bcs + self.bcs = list(bcs) if bcs is not None else [] # per instance, never a shared default self.verbose = verbose self.degree = degree self.continuous = continuous @@ -4129,11 +4129,15 @@ def _check_rule_oversampling(self, degree): local_dofs = fe.getDimension() Nq = len(np.asarray(self.mesh.integration_rule.getData()[1])) if Nq <= local_dofs: + need = self._qdegree_with_at_least(local_dofs + 1) + want = self._qdegree_with_at_least(2 * local_dofs) raise RuntimeError( f"IntegrationPointSemiLagrangian: the mesh rule has {Nq} points per cell " f"but a degree-{degree} history has {local_dofs} local dofs; the " "least-squares fit is not oversampled and is unstable at small Courant " - f"number. Build the mesh with qdegree >= {self.mesh.qdegree + 1}." + f"number. For this cell type and history degree the rule needs at least " + f"qdegree={need} (more points than dofs); 2x oversampling, the verified " + f"setting, is qdegree={want}." ) if Nq < 2 * local_dofs: warnings.warn( @@ -4144,6 +4148,18 @@ def _check_rule_oversampling(self, degree): stacklevel=3, ) + def _qdegree_with_at_least(self, npoints, qmax=12): + """The smallest quadrature degree whose default rule on this mesh's + cell type has at least ``npoints`` points per cell (None if none up + to ``qmax``). Point counts come from PETSc's own rules.""" + for q in range(self.mesh.qdegree + 1, qmax + 1): + fe = PETSc.FE().createDefault( + self.mesh.dim, 1, self.mesh.isSimplex, q, f"ipsl_qscan_{q}_", PETSc.COMM_SELF, + ) + if len(np.asarray(fe.getQuadrature().getData()[1])) >= npoints: + return q + return None + # ------------------------------------------------------------------ @property def psi_fn(self):