diff --git a/docs/developer/index.md b/docs/developer/index.md index 82f67deb..a8c1f512 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -190,6 +190,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 00000000..542ec9f7 --- /dev/null +++ b/docs/developer/subsystems/integration-point-variables.md @@ -0,0 +1,262 @@ +# 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` 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 +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). + +## 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.AdvDiffusionSLCN(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. + +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 + +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, 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`, `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 +(`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. 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 + +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/setup.py b/setup.py index 8d46fcac..924d753b 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_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f56fbe0f..aa22e699 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 new file mode 100644 index 00000000..efbbd553 --- /dev/null +++ b/src/underworld3/cython/petsc_quadrature_fe.pyx @@ -0,0 +1,233 @@ +# 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*) + 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() + 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 + + +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/cython/uw_delta_space.h b/src/underworld3/cython/uw_delta_space.h new file mode 100644 index 00000000..65214c16 --- /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/src/underworld3/discretisation/__init__.py b/src/underworld3/discretisation/__init__.py index 687b3451..a60150b4 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 b0844701..7f045997 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -5484,14 +5484,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 5d8eb7d6..4a277851 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 4ff8861c..385f60cc 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 21fcf641..6a144f25 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/systems/__init__.py b/src/underworld3/systems/__init__.py index 17097bbf..c47a3383 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -93,6 +93,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 from .ddt import EulerianSUPG as EulerianSUPG_DDt diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index f82f58c4..add33178 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -738,6 +738,56 @@ 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_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()) + degs = [fn.meshvar().degree for fn in varfns] + return max(degs) if degs else 2 + + def _make_velocity_level(self, tag): + """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}] }} }}", + 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 (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[...] + def bdf(self, order: Optional[int] = None): r"""Backward differentiation approximation of the time-derivative of :math:`\psi`. @@ -1951,6 +2001,7 @@ def __init__( monotone_mode: Optional[str] = None, theta: float = 0.5, old_frame_traceback: bool = False, + midtime_velocity: bool = True, ): super().__init__() @@ -1962,6 +2013,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" @@ -2662,12 +2716,46 @@ 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`. :math:`v^{n-1}` is + ``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) + if level is None or not getattr(self, "_v_prev_valid", False): + return None + return self._V_matrix() * sympy.Rational(3, 2) - level["expr"] * sympy.Rational(1, 2) + + def _record_velocity_history(self): + """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 + self._copy_velocity_level(self._v_prev_level) + 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. @@ -2699,15 +2787,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 @@ -2968,12 +3057,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 @@ -3185,6 +3279,11 @@ 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. + 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 # plain trace-back. If multiple adapts happen before the next @@ -3879,3 +3978,321 @@ 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_fn`` may be any expression (``-v``, ``v/2``, ``c(t) v``); the + velocity history caches it by evaluation at each time level. + """ + + 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=None, + order: int = 1, + theta: float = 0.5, + monotone_mode: Optional[str] = None, + **_unsupported, + ): + super().__init__() + if vtype != VarType.SCALAR: + raise NotImplementedError( + "IntegrationPointSemiLagrangian: scalar histories only for now" + ) + self.monotone_mode = monotone_mode + self.mesh = mesh + 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 + 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._init_history_tracking(order) + self._check_rule_oversampling(degree) + + if varsymbol is None: + 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) + ] + # 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}) }} }}", + units=psi_units, + ) + 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 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) + + 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. + + The solve fits the sampled departure-point values to the continuous + 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( + 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: + 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. 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( + f"IntegrationPointSemiLagrangian: {Nq} rule points per cell for " + 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, + ) + + 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): + 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(_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(_to_nondim_ndarray(v, units=self._velocity_units())) + if v.ndim == 3: + v = v[:, 0, :] + return v.reshape(coords.shape[0], self.mesh.dim) + + 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_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_start_sym, X, evalf) + Xm = X - 0.5 * dt * v0 + if clamp is not None: + Xm = clamp(Xm) + vm = self._velocity_at(v_mid_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() + 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. + # 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. + V = [lvl["expr"] for lvl in self.v_levels] + if k == 0: + v_start = V[0] + v_mid = V[0] * sympy.Rational(3, 2) - V[1] * half + else: + 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 + ) + 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 + ``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[...] + 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 = 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 + + 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[...] + for k in range(self._n_v - 1, 0, -1): + self._copy_velocity_level(self.v_levels[k], self.v_levels[k - 1]) + 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/src/underworld3/utilities/_jitextension.py b/src/underworld3/utilities/_jitextension.py index 395cc723..5301b43f 100644 --- a/src/underworld3/utilities/_jitextension.py +++ b/src/underworld3/utilities/_jitextension.py @@ -903,7 +903,21 @@ def ccode_patch_fns(varlist, prefix_str, component_offsets=None): 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 component_offsets is not None: u_i = component_offsets[var.field_id] u_x_i = u_i * mesh.cdim @@ -922,7 +936,7 @@ def ccode_patch_fns(varlist, prefix_str, component_offsets=None): 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 @@ -941,7 +955,7 @@ def ccode_patch_fns(varlist, prefix_str, component_offsets=None): 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_0064_quadrature_point_fe.py b/tests/test_0064_quadrature_point_fe.py new file mode 100644 index 00000000..e989750a --- /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 diff --git a/tests/test_0065_integration_point_variable.py b/tests/test_0065_integration_point_variable.py new file mode 100644 index 00000000..4f334436 --- /dev/null +++ b/tests/test_0065_integration_point_variable.py @@ -0,0 +1,164 @@ +"""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) + + +@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") + 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 (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.allclose(own, h.data[:, 0], rtol=0, atol=1e-14) + + # 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 = 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.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. + 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) diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py new file mode 100644 index 00000000..18520e44 --- /dev/null +++ b/tests/test_0066_integration_point_slcn.py @@ -0,0 +1,216 @@ +"""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=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)) + 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.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V, DuDt=DuDt, order=1) + else: + adv = uw.systems.AdvDiffusionSLCN(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()) + 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(): + """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=3 + ) + 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 + + +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 + 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 + # V_fn is symbolic by design: any expression of the variable must work. + # "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) + else: + 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) + ddt.update_post_solve(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 * 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", "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) + # The trace is exact to round-off. Integration-point: the tolerance + # 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-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 + + +@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 diff --git a/tests/test_1056_units_slcn_traceback.py b/tests/test_1056_units_slcn_traceback.py index de2e29a2..c8ae585d 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}"