From bf81ed97122d6dea7626c12dfe223f8e168612d8 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 23 Oct 2026 20:32:14 +1100 Subject: [PATCH 1/6] Add opt-in cylindrical pressure and circular-wall corrections Expose solver-owned pressure extraction, numerical Taylor wall updates, circular CBF topography recovery and modal trace integration for affine isoviscous P2/continuous-P1 annuli. Keep five public methods backed by private lazy helpers, explicit configuration and numerical failure histories. Document technical sources and the required projected pressure gauge. Include the registered-boundary collective datum scheduling prerequisite from open PR #855; no analytical, geoid, Poisson, spherical or child-submesh changes are included. Add independent mathematical regressions and one explicit CI entry. Validation: executable kernel/method/test ASTs match d2c05279 tested on the integration branch (14 serial UW3, 189 annulus components and 14 tests on 16 MPI ranks). Four standalone h32 migrations reproduce calibrated fields and signed coefficients within 1e-14. Coverage and style checks pass. A matching upstream-base native installation and focused serial/MPI gates are next before opening the review PR. Underworld development team with AI support from Claude Code --- .../cylindrical-geometry-correction.md | 146 ++++++ docs/api/solvers.md | 13 + .../guides/boundary-condition-rulings.md | 18 + scripts/test.sh | 1 + .../cython/petsc_generic_snes_solvers.pyx | 203 +++++++++ .../utilities/_cylindrical_pressure.py | 241 ++++++++++ .../utilities/_cylindrical_reactions.py | 143 ++++++ .../utilities/_cylindrical_wall.py | 426 ++++++++++++++++++ src/underworld3/utilities/rotated_bc.py | 38 +- .../test_1084_rotated_datum_collectives.py | 40 ++ tests/test_1089_cylindrical_geometry.py | 334 ++++++++++++++ 11 files changed, 1593 insertions(+), 10 deletions(-) create mode 100644 docs/advanced/cylindrical-geometry-correction.md create mode 100644 src/underworld3/utilities/_cylindrical_pressure.py create mode 100644 src/underworld3/utilities/_cylindrical_reactions.py create mode 100644 src/underworld3/utilities/_cylindrical_wall.py create mode 100644 tests/parallel/test_1084_rotated_datum_collectives.py create mode 100644 tests/test_1089_cylindrical_geometry.py diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md new file mode 100644 index 000000000..966838f55 --- /dev/null +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -0,0 +1,146 @@ +# Circular geometry correction on affine annulus meshes + +These explicitly selected Stokes methods support two-dimensional, affine +triangular annuli with continuous P2 velocity, continuous P1 regular pressure +and constant viscosity. They correct source integration, wall data and reaction +recovery; they do not replace straight elements with curved elements. Numerical +validation currently targets the isoviscous annulus benchmarks. They do not +establish support for spherical shells, variable viscosity or other FE pairs. + +## Pressure extraction for a prescribed circular sheet + +For a radial sheet load at $r=R$, write + +$$ +p=q+p_J,\qquad p_J=-A H(r-R)b(r)\cos(n\theta),\qquad +b(r)=\frac{r_o-r}{r_o-R}. +$$ + +In $-\nabla\cdot\tau+\nabla p=f$, the distributional gradient of $p_J$ +cancels $f=-A\cos(n\theta)e_r\delta(r-R)$. The regular equation uses +$A H(r-R)\nabla[b\cos(n\theta)]$. Circular-cut integration constructs six +P2 load moments on each outer triangle. A DG P2 field represents these moments +under the native quadrature; it is a force representation, not a new velocity +or pressure element pair. + +```python +lift = stokes.cylindrical_pressure_lifting( + interface_radius=2.0, outer_radius=2.22, wavenumber=2, + amplitude=1.0, quadrature_points=12, +) +stokes.bodyforce = lift.regular_force.sym +stokes.add_natural_bc(lift.natural_datum, "Upper") +physical_pressure = stokes.p.sym[0] + lift.pressure_jump +``` + +The natural datum is $-p_J\Gamma$ on the polygonal outer wall. Although +$b(r_o)=0$, $b$ is nonzero inside each outer chord. Dropping this term changes +the regular problem. `pressure_jump_integral` permits physical pressure-gauge +calibration; it is the owned-cell integral of the prescribed extraction. +Smooth density problems use their original volume force and physical P1 pressure; +they do not require a pressure lift. + +## Numerical transfer of circular free-slip conditions + +At a polygon wall radius $r_h$, let $d=R-r_h$. Lagged numerical traces impose +$u_r(r_h)=-d\partial_r u_r$ and replace shear stress by +$\sigma_{rt}(r_h)=-d\partial_r\sigma_{rt}$. Diagonal stresses remain numerical. +The force moments are projected tangentially before inversion of the consistent +P2 line mass. Discrete normal flux and rigid-rotation torque are projected to zero. +No analytical Stokes solution supplies the wall data. + +```python +stokes.petsc_use_pressure_nullspace = True +walls = {} +for boundary, radius, sign in (("Lower", 1.22, -1), ("Upper", 2.22, 1)): + walls[boundary] = stokes.add_cylindrical_wall_transfer( + boundary=boundary, radius=radius, sign=sign, viscosity=1.0, + interface_radius=2.0, outer_radius=2.22, wavenumber=2, + amplitude=1.0, quadrature=12, + ) +history = stokes.solve_cylindrical_wall_updates(walls=walls, updates=4) +``` + +Create wall fields before the first solve and before another free-slip condition +is registered on these boundaries. Use `amplitude=0` for smooth physical pressure. +Four updates produce five accepted solves, with final fields retained from the +last solve. A converged, zero-increment warm start is accepted. Failed histories +remain on `stokes.cylindrical_wall_history`; file output belongs to the caller. +Enable `stokes.petsc_use_pressure_nullspace = True` before the first solve. +The wall history reads the projected pressure-gauge residual produced by that +configuration; other pressure-gauge configurations are outside this helper's +supported contract. The helper does not install the pressure projection or +velocity nullspace itself. Configure the rigid-rotation velocity modes separately +when required by the enclosing free-slip system. + +## Spatial circular reaction and topography + +The circular reaction moments are + +$$ +R_{\rm circle}=R_{\rm CBF}+L_{\rm circle}(\sigma_h)-L_{\rm polygon}(\sigma_h). +$$ + +The circular stress uses the first-order radial Taylor transfer. Each geometric +load has its own mean removed; the circular consistent line mass determines the +spatial trace. The original assembled CBF contribution is preserved, including +its numerical residual, and its solver cache is not overwritten. + +```python +checks = stokes.recover_cylindrical_topography( + wall=walls["Upper"], field=outer_topography, buoyancy_scale=-1.0, +) +``` + +The target must be a continuous scalar P2 field created before solving. The +method writes only this wall and returns mass-moment, mean and length checks. +Direct-field Poisson consumes the spatial field. Spectral gravity subsequently +extracts a mode using + +```python +coefficient = stokes.cylindrical_boundary_coefficient( + wall=walls["Upper"], field=outer_topography, wavenumber=2, +) +``` + +Coefficient extraction removes the circle-weighted mean and uses the exact +radial-projection line Jacobian. It is not part of spatial reaction recovery. + +## Parallel scope and reproducibility + +Call the methods collectively. Only owned boundary-cell traces and boundary +nodes are replicated; full volume fields are never gathered. Small SciPy sparse +solves invert the replicated one-dimensional boundary mass systems. The Stokes +operator and volume solves remain distributed PETSc. Wall handles cache geometry +and are valid only for their original mesh, fields and fixed parameters. + +Record pressure representation, quadrature, source amplitude, wall-update history, +reaction checks and the target actually supplied to gravity. Keep complete spatial +errors distinct from modal errors. A correction changes the numerical problem +on an unchanged affine mesh; both baseline and corrected cases must be solved +independently for a scientific before/after comparison. + +## References and attribution boundaries + +- Chang, K. S., and Kwak, D. Y. (2009), *Extraction method for Stokes Flow with + jumps in the pressure*, [arXiv:0911.4769](https://arxiv.org/abs/0911.4769). + Section 2, equations (4)-(6), motivates regular/singular pressure splitting. + UW3 uses its own prescribed radial profile, exact-circle cap integration, + DG P2 load moments and outer-chord natural datum, rather than the paper's + harmonic extension and CR discretisation. +- Atallah, N. M., Canuto, C., and Scovazzi, G. (2020), *The second-generation + Shifted Boundary Method and its numerical analysis*, Computer Methods in + Applied Mechanics and Engineering 372, 113341, + [doi:10.1016/j.cma.2020.113341](https://doi.org/10.1016/j.cma.2020.113341), + [author preprint](https://arxiv.org/abs/2004.10584). Section 2.2, equations + (8)-(13), supplies the Taylor-shift motivation. Its Section 4.2 uses a + Nitsche formulation and assumes body-fitted Neumann boundaries. UW3's lagged + strong free-slip, tangential load projection and compatibility corrections + are separate constructions; that analysis provides no stability proof here. +- Zhong, S., Gurnis, M., and Hulbert, G. (1993), *Accurate determination of + surface normal stress in viscous flow from a consistent boundary flux + method*, Physics of the Earth and Planetary Interiors 78(1-2), 1-8, + [doi:10.1016/0031-9201(93)90078-N](https://doi.org/10.1016/0031-9201(93)90078-N). + This is the CBF stress-recovery reference. UW3's circular moment transfer, + Taylor reaction correction and weighted coefficient integration extend the + calculation and are not results established by that paper. diff --git a/docs/api/solvers.md b/docs/api/solvers.md index ac7f51a53..3eca569d4 100644 --- a/docs/api/solvers.md +++ b/docs/api/solvers.md @@ -10,6 +10,19 @@ :show-inheritance: ``` +### Affine-annulus circular geometry correction + +The Stokes methods `cylindrical_pressure_lifting`, +`add_cylindrical_wall_transfer`, `solve_cylindrical_wall_updates`, +`recover_cylindrical_topography` and `cylindrical_boundary_coefficient` +provide explicitly selected source, wall and reaction corrections for +constant-viscosity affine 2D P2/continuous-P1 annuli. They lazily load private +implementations; no benchmark package is required. The +[circular geometry guide](../advanced/cylindrical-geometry-correction.md) +contains the mathematical contracts and examples. Four wall updates mean five +accepted solves. Spatial topography recovery is independent of the optional +spectral coefficient calculation. + ### SNES_VE_Stokes Viscoelastic extension of the Stokes solver. diff --git a/docs/developer/guides/boundary-condition-rulings.md b/docs/developer/guides/boundary-condition-rulings.md index b5b7f9605..df385384b 100644 --- a/docs/developer/guides/boundary-condition-rulings.md +++ b/docs/developer/guides/boundary-condition-rulings.md @@ -87,3 +87,21 @@ given as `None` (or `sympy.oo`) are left free. Three things to know: - A boundary condition on a field that is not the solver's unknown (a pressure datum on a Stokes solver) is accepted by the solver but not by the adjoint, which refuses rather than drop it. + +## Explicit circular-wall correction on affine annuli + +For the restricted constant-viscosity, continuous P2/P1 annulus workflow, +`stokes.add_cylindrical_wall_transfer(...)` registers numerical normal and +natural data together; `stokes.solve_cylindrical_wall_updates(...)` alternates +accepted solves and lagged boundary updates. Do not also register homogeneous +rotated conditions on the same walls. The correction follows the physical +circle while retaining affine elements, so it is explicitly selected rather +than a change to the default rotated condition. + +`stokes.cylindrical_pressure_lifting(...)` constructs a prescribed circular-sheet +pressure extraction, including the nonzero outer-chord natural datum. It leaves +force and datum registration visible to the caller. +`stokes.recover_cylindrical_topography(...)` transfers CBF moments into a spatial +P2 target without changing the raw reaction cache or fitting harmonics. See the +[circular geometry guide](../../advanced/cylindrical-geometry-correction.md) +for equations, limitations and source/pressure conventions. diff --git a/scripts/test.sh b/scripts/test.sh index b6caaedda..9f81360d3 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -193,6 +193,7 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # hardened, fast, and exactly what CI exists to run. test_1064 is the slow one # and it is excluded by WHAT IT IS (tier_c) below, not by its number. "${PYTEST[@]}" tests/test_106*py tests/test_107*py || status=1 + "${PYTEST[@]}" tests/test_1089_cylindrical_geometry.py || status=1 # # test_1072 used to be pulled forward out of test_107* by name, because that # band was not batched and it is the ONLY guard on the 3-D free-surface sign diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 5e3c23dda..23065e675 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8652,6 +8652,209 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): J1.destroy(); J2.destroy() return rel > tol + def cylindrical_pressure_lifting(self, *, interface_radius, outer_radius, + wavenumber, amplitude, quadrature_points): + r"""Construct an exact-circle pressure lift and moment-equivalent force. + + Parameters + ---------- + interface_radius, outer_radius : float + Prescribed circular sheet and outer wall radii, with outer > sheet. + wavenumber : int + Azimuthal cosine mode of the prescribed radial sheet load. + amplitude : float + Sheet amplitude in :math:`f=-A\cos(n\theta)e_r\delta(r-R)`. + quadrature_points : int + Points per direction in circular-cut quadrature. + + Returns + ------- + named tuple + ``pressure_jump``, ``regular_force`` (DG P2 MeshVariable), + ``natural_datum``, ``pressure_jump_integral`` and ``checks``. + Reconstruct physical pressure as :math:`p=q+p_J`. Assign the + force and register ``natural_datum`` on the outer boundary; the + latter is :math:`-p_J\Gamma` and does not vanish on outer chords. + + Notes + ----- + Explicitly selected correction for affine 2D P2/continuous-P1 annuli, + constant viscosity and a fitted circular sheet. No analytical Stokes + solution enters construction. Call collectively on every rank. + + References + ---------- + K. S. Chang and D. Y. Kwak (2009), "Extraction method for Stokes + Flow with jumps in the pressure", arXiv:0911.4769, + https://arxiv.org/abs/0911.4769, Section 2, equations (4)-(6). + The pressure split motivates this method. The radial lift, circular-cut + load integration and outer-chord natural datum are UW3 adaptations, + not the paper's harmonic-extension/CR scheme. + """ + from underworld3.utilities._cylindrical_pressure import _create_circular_pressure_lift + if self.u.degree != 2 or not self.u.continuous or self.p.degree != 1 or not self.p.continuous: + raise ValueError("Circular pressure lifting requires continuous P2/P1 fields.") + return _create_circular_pressure_lift(self.mesh, interface_radius, outer_radius, + wavenumber, amplitude, quadrature_points) + + def add_cylindrical_wall_transfer(self, *, boundary, radius, sign, viscosity, + interface_radius, outer_radius, wavenumber, + amplitude, quadrature): + r"""Register numerical Taylor data for a circular free-slip wall. + + Parameters + ---------- + boundary : str + Exterior boundary label, before registering another condition on it. + radius : float + Physical circular wall radius. + sign : {-1, 1} + Outward radial orientation: -1 for inner, +1 for outer. + viscosity : float + Constant viscosity used by this Stokes system. + interface_radius, outer_radius : float + Radii of the prescribed pressure lift and outer wall. + wavenumber : int + Prescribed cosine sheet mode. + amplitude : float + Pressure-lift amplitude; use zero for smooth physical pressure. + quadrature : int + Number of line quadrature points. + + Returns + ------- + wall + Opaque numerical wall handle, with ``normal_field`` and + ``traction_field``. Pass it to wall updates and reaction recovery. + + Notes + ----- + Restricted to affine 2D P2/continuous-P1 annuli and constant viscosity. + Registers rotated normal and tangential natural data. Only owned + boundary-cell traces are replicated; no volume field is gathered. + + References + ---------- + N. M. Atallah, C. Canuto and G. Scovazzi (2020), "The + second-generation Shifted Boundary Method and its numerical analysis", + Computer Methods in Applied Mechanics and Engineering 372, 113341, + https://doi.org/10.1016/j.cma.2020.113341; arXiv:2004.10584, + Section 2.2, equations (8)-(13). + These equations motivate Taylor transfer. Strong free-slip constraints, + lagged numerical traction and flux/torque projection are UW3 additions. + This is not that paper's Nitsche scheme; its Section 4.2 assumes + body-fitted Neumann boundaries and does not prove this method stable. + """ + from underworld3.utilities._cylindrical_wall import _add_cylindrical_wall_transfer + return _add_cylindrical_wall_transfer( + self, boundary=boundary, radius=radius, sign=sign, viscosity=viscosity, + interface_radius=interface_radius, outer_radius=outer_radius, + wavenumber=wavenumber, amplitude=amplitude, quadrature=quadrature) + + def solve_cylindrical_wall_updates(self, *, walls, updates): + """Alternate accepted Stokes solves and numerical wall-data updates. + + Parameters + ---------- + walls : dict + Boundary labels mapped to handles returned by + :meth:`add_cylindrical_wall_transfer`. + updates : int + Number of updates; four produces five solves. + + Returns + ------- + list of dict + Convergence and boundary-update histories. Partial history remains + available as ``cylindrical_wall_history`` if an iterate fails. This + method performs no file I/O. A residual-accepted zero-increment + warm start counts as an accepted solve. + + Notes + ----- + Before calling, set ``stokes.petsc_use_pressure_nullspace = True``. + The accepted-solve history reads the rotated solver's projected + pressure-gauge residual, which is recorded only with this projection + enabled. This method does not configure the gauge automatically; + other pressure-gauge configurations are outside its supported contract. + + This lagged solve/update schedule is a UW3 implementation choice. + The Taylor motivation is Atallah, Canuto and Scovazzi (2020), + "The second-generation Shifted Boundary Method and its numerical + analysis", CMAME 372, 113341, + https://doi.org/10.1016/j.cma.2020.113341, Section 2.2. + Four updates and their convergence are verified by benchmark evidence; + they are not guarantees provided by that paper. + """ + from underworld3.utilities._cylindrical_wall import _solve_with_wall_updates + return _solve_with_wall_updates(self, walls, updates) + + def recover_cylindrical_topography(self, *, wall, field, buoyancy_scale): + r"""Transfer CBF reaction moments to a circular spatial P2 topography. + + Parameters + ---------- + wall : wall + Handle returned by :meth:`add_cylindrical_wall_transfer`. + field : scalar continuous P2 MeshVariable + Spatial target; only this boundary's values are written. + buoyancy_scale : float + Nonzero signed :math:`\Delta\rho g`, with :math:`h=-R/(\Delta\rho g)`. + + Returns + ------- + dict + Circular length, moment, mean and correction checks. The original + assembled CBF reaction cache remains unchanged. No coefficients or + reference fields enter this transfer. + + References + ---------- + S. Zhong, M. Gurnis and G. Hulbert (1993), "Accurate determination + of surface normal stress in viscous flow from a consistent boundary + flux method", Physics of the Earth and Planetary Interiors 78(1-2), + 1-8, https://doi.org/10.1016/0031-9201(93)90078-N. + Cited for momentum-residual CBF stress recovery. The circular-minus- + polygonal moment transfer and Taylor correction here are UW3 + adaptations; this reference does not analyse those corrections. + """ + from underworld3.utilities._cylindrical_reactions import _recover_circular_topography + return _recover_circular_topography(self, wall, field, buoyancy_scale) + + def cylindrical_boundary_coefficient(self, *, wall, field, wavenumber): + r"""Integrate a spatial P2 wall field's cosine coefficient on the circle. + + Parameters + ---------- + wall : wall + Handle returned by :meth:`add_cylindrical_wall_transfer`. + field : scalar continuous P2 MeshVariable + Spatial field on the same mesh and P2 nodes as the wall. + wavenumber : int + Nonzero cosine harmonic. The spatial mean is removed first. + + Returns + ------- + float + Circle-weighted least-squares cosine coefficient. This operation + is for spectral postprocessing; Direct-field Poisson uses the + spatial field directly. Call collectively on every rank. + + Notes + ----- + This circle-weighted integration of an existing spatial field is a UW3 + postprocessing construction. The underlying CBF provenance is Zhong, + Gurnis and Hulbert (1993), "Accurate determination of surface normal + stress in viscous flow from a consistent boundary flux method", + Physics of the Earth and Planetary Interiors 78(1-2), 1-8, + https://doi.org/10.1016/0031-9201(93)90078-N. That paper does not + define the circle-transfer Jacobian or this coefficient integration. + """ + from underworld3.utilities._cylindrical_reactions import _circular_field_coefficient + if wavenumber <= 0: + raise ValueError("A mean-removed boundary coefficient requires wavenumber > 0.") + return _circular_field_coefficient(wall, field, wavenumber) + def boundary_normal_traction(self, boundary, mass="auto"): r"""Return the boundary normal traction :math:`\sigma_{nn}` on a rotated-free-slip ``boundary`` as the constraint reaction from the last diff --git a/src/underworld3/utilities/_cylindrical_pressure.py b/src/underworld3/utilities/_cylindrical_pressure.py new file mode 100644 index 000000000..212c97e24 --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_pressure.py @@ -0,0 +1,241 @@ +"""Private circular source extraction for affine two-dimensional Stokes meshes. + +The caller selects this mathematical correction explicitly. No reference Stokes +solution enters the source. Owned-cell checks are reduced once; force values +also populate local ghost cells. Public entry: Stokes.cylindrical_pressure_lifting. + +Pressure decomposition is motivated by Chang and Kwak (2009), "Extraction method +for Stokes Flow with jumps in the pressure", arXiv:0911.4769, +https://arxiv.org/abs/0911.4769, Section 2, equations (4)-(6). The prescribed +radial lift, exact-circle cap integration, DG P2 moment representation and +outer-chord natural datum are this UW3 implementation. It does not implement +the paper's Section 3 harmonic extension or Section 4 CR discretisation. +""" + +from typing import NamedTuple + +import numpy as np +import sympy as sp +import underworld3 as uw +from mpi4py import MPI +from petsc4py import PETSc + + +class _PressureLift(NamedTuple): + pressure_jump: object + regular_force: object + natural_datum: object + pressure_jump_integral: float + checks: dict + + +def _p2_shapes(barycentric: np.ndarray) -> np.ndarray: + """P2 values in vertex order 0,1,2 followed by edges 01,12,20.""" + barycentric = np.asarray(barycentric) + return np.stack( + [barycentric[..., i] * (2 * barycentric[..., i] - 1) for i in range(3)] + + [4 * barycentric[..., i] * barycentric[..., j] for i, j in ((0, 1), (1, 2), (2, 0))], + axis=-1, + ) + +def _profile_values(radius: np.ndarray, interface_radius: float, width: float): + """Linear radial lifting and derivative, extended below R for cap integration.""" + s = (np.asarray(radius) - interface_radius) / width + return 1 - s, np.full_like(s, -1 / width) + +def _extended_regular_force(points, interface_radius, width, wavenumber, amplitude): + """Smooth outer branch of A*grad(b*cos(n*theta)), including its cap extension.""" + radius = np.linalg.norm(points, axis=-1) + theta = np.arctan2(points[..., 1], points[..., 0]) + radial = points / radius[..., None] + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + value, derivative = _profile_values(radius, interface_radius, width) + return amplitude * ( + (derivative * np.cos(wavenumber * theta))[..., None] * radial + - (wavenumber * value * np.sin(wavenumber * theta) / radius)[..., None] * tangent + ) + +def _load_moments(vertices, interface_radius, width, wavenumber, amplitude, quadrature_points): + """P2 force moments on outer fitted triangles, excluding each circular cap. + + Interior quadrature integrates the smooth outer branch on each triangle. + Polar cap quadrature removes the region rnj", bary, vertices) + force = _extended_regular_force(point, interface_radius, width, wavenumber, amplitude) + moments += ( + (area * 2 * x * wx * wy)[:, None, None] + * _p2_shapes(bary)[None, :, None] + * force[:, None, :] + ) + radius = np.linalg.norm(point, axis=1) + lift_integral -= ( + area + * 2 + * x + * wx + * wy + * amplitude + * (1 - (radius - interface_radius) / width) + * np.cos(wavenumber * np.arctan2(point[:, 1], point[:, 0])) + ) + + cap_area = np.zeros(len(vertices)) + chord_count = np.zeros(len(vertices), dtype=int) + radii = np.linalg.norm(vertices, axis=-1) + for first, second in ((0, 1), (1, 2), (2, 0)): + selected = np.flatnonzero( + np.isclose(radii[:, first], interface_radius, rtol=0, atol=1e-10) + & np.isclose(radii[:, second], interface_radius, rtol=0, atol=1e-10) + ) + if not len(selected): + continue + left, right = vertices[selected, first], vertices[selected, second] + normal = left + right + normal /= np.linalg.norm(normal, axis=1)[:, None] + tangent = np.column_stack((-normal[:, 1], normal[:, 0])) + distance = np.sum(left * normal, axis=1) + half_angle = np.arctan2(np.linalg.norm(right - left, axis=1) / 2, distance) + measured_area = np.zeros(len(selected)) + for t, wt in zip(nodes, weights): + angle = (2 * t - 1) * half_angle + direction = np.cos(angle)[:, None] * normal + np.sin(angle)[:, None] * tangent + lower = distance / np.cos(angle) + thickness = interface_radius - lower + for s, ws in zip(nodes, weights): + radius = lower + s * thickness + point = radius[:, None] * direction + xi = np.einsum("nij,nj->ni", inverse[selected], point - vertices[selected, 0]) + bary = np.column_stack((1 - xi.sum(axis=1), xi)) + if bary.min() < -1e-9: + raise ValueError("Circular cap leaves its adjacent affine triangle.") + weight = radius * thickness * ws * 2 * half_angle * wt + force = _extended_regular_force( + point, interface_radius, width, wavenumber, amplitude + ) + moments[selected] -= ( + weight[:, None, None] * _p2_shapes(bary)[..., None] * force[:, None, :] + ) + lift_integral[selected] += ( + weight + * amplitude + * (1 - (radius - interface_radius) / width) + * np.cos(wavenumber * np.arctan2(point[:, 1], point[:, 0])) + ) + measured_area += weight + expected_area = interface_radius**2 * half_angle - distance**2 * np.tan(half_angle) + if not np.allclose(measured_area, expected_area, rtol=1e-8, atol=1e-13): + raise ArithmeticError("Polar quadrature did not reproduce circular-cap areas.") + cap_area[selected] += measured_area + chord_count[selected] += 1 + return moments, cap_area, chord_count, lift_integral + +def _create_circular_pressure_lift( + mesh, interface_radius, outer_radius, wavenumber, amplitude, quadrature_points +): + """Return exact-circle p_jump, a moment-equivalent DG P2 force, and global checks. + + A linear radial profile avoids introducing a narrow transition layer. + Local cell connectivity includes ghosts; each DG cell is populated locally, + then owned-cell diagnostics are summed once across MPI ranks. + + References + ---------- + K. S. Chang and D. Y. Kwak (2009), "Extraction method for Stokes + Flow with jumps in the pressure", arXiv:0911.4769, + https://arxiv.org/abs/0911.4769, Section 2, equations (4)-(6). + The regular-plus-singular pressure split motivates this construction; + its radial profile, cut loads and boundary datum are UW3 adaptations. + """ + if mesh.dim != 2 or mesh.cdim != 2: + raise ValueError("Circular pressure extraction requires a two-dimensional affine mesh.") + if not 0 < interface_radius < outer_radius: + raise ValueError("Require 0 < interface_radius < outer_radius.") + if quadrature_points < 2: + raise ValueError("Circular-cut quadrature requires at least two points.") + width = outer_radius - interface_radius + force_field = uw.discretisation.MeshVariable( + "CircleRegularForce", mesh, 2, degree=2, continuous=False + ) + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, 3) + vertices = np.asarray(mesh._get_coords_for_basis(1, True))[vertex_rows] + outer = np.linalg.norm(vertices.mean(axis=1), axis=1) > interface_radius + moments, cap_area, chord_count, lift_integral = _load_moments( + vertices[outer], interface_radius, width, wavenumber, amplitude, quadrature_points + ) + edges = np.stack((vertices[:, 1] - vertices[:, 0], vertices[:, 2] - vertices[:, 0]), axis=-1) + area = np.abs(np.linalg.det(edges)) / 2 + + # Match native P2 weak loads, rather than interpolate pointwise source values. + fe = PETSc.FE().createDefault(2, 1, True, mesh.qdegree, "circle_force_mass_") + points, weights = fe.getQuadrature().getData() + xy = (np.asarray(points).reshape(-1, 2) + 1) / 2 + bary = np.column_stack((1 - xy.sum(axis=1), xy)) + shapes = _p2_shapes(bary) + mass = np.einsum("qi,qj,q->ij", shapes, shapes, np.asarray(weights).reshape(-1) / 2) + fe.destroy() + if np.linalg.matrix_rank(mass) != 6: + raise ValueError("Quadrature must integrate the full P2 mass matrix (qdegree >= 2).") + nodal = np.zeros((len(vertices), 6, 2)) + nodal[outer] = np.einsum("ij,njk->nik", np.linalg.inv(mass), moments / area[outer, None, None]) + reproduced = area[outer, None, None] * np.einsum("ij,njk->nik", mass, nodal[outer]) + moment_error = float(np.max(np.abs(reproduced - moments))) if len(moments) else 0.0 + if moment_error > 1e-11: + raise ArithmeticError(f"DG P2 load-moment representation failed: {moment_error}") + + rows = mesh._cell_node_indices(2, False).reshape(-1, 6) + coordinates = np.asarray(force_field.coords)[rows] + xi = np.einsum("nij,nkj->nki", np.linalg.inv(edges), coordinates - vertices[:, None, 0]) + bary = np.concatenate((1 - xi.sum(axis=-1, keepdims=True), xi), axis=-1) + force_field.array[rows, 0, :] = np.einsum("nki,nij->nkj", _p2_shapes(bary), nodal) + owned = mesh._get_owned_cells_mask()[outer] + expected_force = uw.mpi.comm.allreduce(moments[owned].sum(axis=(0, 1)), op=MPI.SUM) + actual_force = np.array( + [uw.maths.Integral(mesh, force_field.sym[i]).evaluate() for i in range(2)] + ) + force_error = float(np.max(np.abs(actual_force - expected_force))) + if force_error > 1e-10: + raise ArithmeticError(f"Native DG force integral differs from the cut load: {force_error}") + checks = dict( + width=width, + quadrature_points=quadrature_points, + circular_segment_area=float(uw.mpi.comm.allreduce(cap_area[owned].sum(), op=MPI.SUM)), + internal_chord_count=int(uw.mpi.comm.allreduce(chord_count[owned].sum(), op=MPI.SUM)), + maximum_load_moment_error=float(uw.mpi.comm.allreduce(moment_error, op=MPI.MAX)), + pressure_lift_integral=float(uw.mpi.comm.allreduce(lift_integral[owned].sum(), op=MPI.SUM)), + native_force_integral_error=force_error, + ) + radius, theta = mesh.CoordinateSystem.xR + interface_radius, width, amplitude = [sp.nsimplify(value, rational=True) for value in (interface_radius, width, amplitude)] + s = (radius - interface_radius) / width + lift = sp.Piecewise( + (-amplitude * (1 - s) * sp.cos(wavenumber * theta), radius > interface_radius), + (0, True), + ) + return _PressureLift(lift, force_field, -lift * mesh.Gamma, checks["pressure_lift_integral"], checks) + + +def _pressure_jump_values(points, interface_radius, outer_radius, wavenumber, amplitude): + """Prescribed pressure extraction and its smooth one-sided Cartesian gradient.""" + radius = np.linalg.norm(points, axis=-1) + above = radius > interface_radius + width = outer_radius - interface_radius + profile, _ = _profile_values(radius, interface_radius, width) + pressure = -amplitude * above * profile * np.cos( + wavenumber * np.arctan2(points[..., 1], points[..., 0])) + gradient = -(above[..., None] * _extended_regular_force( + points, interface_radius, width, wavenumber, amplitude)) + return pressure, gradient diff --git a/src/underworld3/utilities/_cylindrical_reactions.py b/src/underworld3/utilities/_cylindrical_reactions.py new file mode 100644 index 000000000..0ec63323d --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_reactions.py @@ -0,0 +1,143 @@ +"""Private spatial reaction transfer from polygonal to circular annulus walls. + +Only boundary traces are replicated. Sparse solves are one-dimensional line +mass systems; no volume field is gathered. The spatial target is independent +of Fourier fitting and the original solver reaction cache remains unchanged. + +The consistent boundary flux (CBF) principle recovers surface stress from the +discrete momentum balance [1]_. The circular-minus-polygonal stress moments and +Taylor shift below are UW3 geometry adaptations, not a scheme or stability +result established in that paper. + +References +---------- +.. [1] Zhong, S., Gurnis, M., and Hulbert, G. (1993). Accurate determination + of surface normal stress in viscous flow from a consistent boundary flux + method. Physics of the Earth and Planetary Interiors, 78(1-2), 1-8. + https://doi.org/10.1016/0031-9201(93)90078-N +""" + +import numpy as np +import underworld3 as uw +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.spatial import cKDTree + +from ._cylindrical_wall import _trace_rule, _numerical_trace + + +def _circular_reaction_field(edges, trace, raw_reaction, radius, sign, quadrature): + """Transfer CBF nodal moments to a circular P2 trace without Fourier fitting. + + R_circle = R_CBF + L_circle(sigma_h) - L_polygon(sigma_h). Remove each + traction mean on its own geometry, then invert the circular consistent + line mass. The circular stress uses the same first-order Taylor shift as + the coefficient diagnostic. Keep the assembled CBF residual contribution. + + The residual-based stress recovery follows the CBF principle of Zhong, + Gurnis and Hulbert (1993), cited in the module. The additive moment transfer + is our extension to a physical circular trace; its accuracy is verified + independently and no published stability proof is claimed for it. + """ + rule = _trace_rule(edges, sign, 2, quadrature) + t = (np.polynomial.legendre.leggauss(quadrature)[0] + 1) / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + edge_nodes = np.stack((edges[:, 0], edges[:, 1], edges.mean(axis=1)), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + normal_nodes = sign * edge_nodes / np.linalg.norm(edge_nodes, axis=-1)[..., None] + r = np.linalg.norm(rule["points"], axis=-1) + radial = rule["points"] / r[..., None] + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + nr = np.sum(rule["normal"][:, None, :] * radial, axis=-1) + nt = np.sum(rule["normal"][:, None, :] * tangent, axis=-1) + traction = (trace["rr"] * nr + trace["rt"] * nt)[..., None] * radial + ( + trace["rt"] * nr + trace["tt"] * nt + )[..., None] * tangent + polygon_load = -np.einsum("eq,qi,eqd,eid->ei", rule["weight"], shape, traction, normal_nodes) + circle_weight = rule["weight"] * radius / r * np.abs(nr) + shifted_stress = trace["rr"] + (radius - r) * trace["rr_derivative"] + circle_load = -np.einsum("eq,qi,eq->ei", circle_weight, shape, shifted_stress) + masses, loads = [], [] + for weight, load in ((rule["weight"], polygon_load), (circle_weight, circle_load)): + local = np.einsum("eq,qi,qj->eij", weight, shape, shape) + rows = np.broadcast_to(index[:, :, None], local.shape).ravel() + columns = np.broadcast_to(index[:, None, :], local.shape).ravel() + mass = coo_matrix((local.ravel(), (rows, columns)), shape=(len(nodes), len(nodes))).tocsr() + total = np.zeros(len(nodes)) + np.add.at(total, index.ravel(), load.ravel()) + uniform = mass @ np.ones(len(nodes)) + total -= total.sum() / uniform.sum() * uniform + masses.append(mass) + loads.append(total) + uniform = masses[0] @ np.ones(len(nodes)) + raw = np.asarray(raw_reaction).reshape(-1) + raw_mean = float(uniform @ raw / uniform.sum()) + target = masses[0] @ (raw - raw_mean) + loads[1] - loads[0] + corrected = spsolve(masses[1], target) + checks = dict( + raw_reaction_mean_removed=raw_mean, + maximum_mass_moment_error=float(np.max(np.abs(masses[1] @ corrected - target))), + circular_reaction_integral=float(np.sum(masses[1] @ corrected)), + mapped_circle_length=float(circle_weight.sum()), + correction_dof_l2=float(np.linalg.norm(corrected - raw)), + ) + np.testing.assert_allclose(checks["mapped_circle_length"], 2 * np.pi * radius, rtol=1e-12) + if max(checks["maximum_mass_moment_error"], abs(checks["circular_reaction_integral"])) > 1e-11: + raise ArithmeticError(f"Circular reaction field moment/gauge failure: {checks}") + return nodes, corrected, checks + +def _validate_spatial_target(wall, field): + """Boundary masks are valid only on the wall's own scalar P2 coordinates.""" + if (field.mesh is not wall.mesh or field.degree != 2 + or not field.continuous or field.num_components != 1): + raise ValueError("Circular wall fields require a continuous scalar P2 target on the same mesh.") + + +def _recover_circular_topography(stokes, wall, field, buoyancy_scale): + r"""Write circular-wall topography from CBF reactions and numerical moments. + + Recover polygonal normal traction from the discrete Stokes momentum + residual using a consistent boundary mass, following the CBF principle + of Zhong, Gurnis and Hulbert (1993; module reference). Apply our circular + moment transfer, then use the signed buoyancy balance + :math:`h=-\sigma_{nn}/(\Delta\rho g)`. No projected volume stress or + analytical Stokes field replaces the assembled reaction. + """ + if buoyancy_scale == 0: + raise ValueError("Topography requires nonzero buoyancy_scale.") + _validate_spatial_target(wall, field) + coords, values = stokes.boundary_normal_traction(wall.name, mass="consistent") + gathered_coords = np.concatenate(uw.mpi.comm.allgather(np.asarray(coords).reshape(-1, 2))) + gathered_values = np.concatenate(uw.mpi.comm.allgather(np.asarray(values).reshape(-1))) + distance, index = cKDTree(gathered_coords).query(wall.nodes) + if distance.max() > 1e-10: + raise ValueError(f"CBF nodes do not cover circular wall {wall.name}.") + velocity, pressure = wall.cell_values() + trace = _numerical_trace(wall.vertices, velocity, pressure, wall.rule["points"], *wall.physics) + nodes, reaction, checks = _circular_reaction_field( + wall.edges, trace, gathered_values[index], wall.radius, wall.sign, wall.quadrature + ) + np.testing.assert_array_equal(nodes, wall.nodes) + distance, index = cKDTree(nodes).query(np.asarray(field.coords)) + boundary = distance < 1e-10 + field.array[boundary, 0, 0] = -reaction[index[boundary]] / buoyancy_scale + return checks + +def _circular_field_coefficient(wall, field, mode): + """Integrate a saved spatial P2 wall field on the physical circle (spectral only).""" + _validate_spatial_target(wall, field) + boundary = wall.local_boundary + coords = np.concatenate(uw.mpi.comm.allgather(np.asarray(field.coords)[boundary])) + values = np.concatenate(uw.mpi.comm.allgather(np.asarray(field.array)[boundary, 0, 0])) + distance, index = cKDTree(coords).query(wall.nodes) + np.testing.assert_allclose(distance, 0, atol=1e-10) + t = (np.polynomial.legendre.leggauss(wall.quadrature)[0] + 1) / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + samples = np.einsum("qi,ei->eq", shape, values[index][wall.node_index]) + r = np.linalg.norm(wall.rule["points"], axis=-1) + nr = np.sum(wall.rule["normal"][:, None, :] * wall.rule["points"] / r[..., None], axis=-1) + weights = wall.rule["weight"] * wall.radius / r * np.abs(nr) + harmonic = np.cos(mode * np.arctan2(wall.rule["points"][..., 1], wall.rule["points"][..., 0])) + samples -= np.sum(weights * samples) / weights.sum() + return float(np.sum(weights * samples * harmonic) / np.sum(weights * harmonic**2)) diff --git a/src/underworld3/utilities/_cylindrical_wall.py b/src/underworld3/utilities/_cylindrical_wall.py new file mode 100644 index 000000000..305899c9b --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_wall.py @@ -0,0 +1,426 @@ +"""Private lagged wall transfer on affine P2/continuous-P1 annuli. + +Only owned boundary-cell traces are replicated across ranks. The sparse SciPy +solves below are bounded one-dimensional boundary mass systems, not volume +operators; the Stokes system remains distributed PETSc. Parallel collectives +are internal implementation details, never part of the notebook contract. +No analytical Stokes data or coefficient fitting enters a numerical wall update. + +Taylor boundary transfer is motivated by N. M. Atallah, C. Canuto and G. Scovazzi +(2020), "The second-generation Shifted Boundary Method and its numerical +analysis", Computer Methods in Applied Mechanics and Engineering 372, 113341, +https://doi.org/10.1016/j.cma.2020.113341, Section 2.2, equations (8)-(13). +This lagged strong free-slip/traction projection is a UW3 adaptation, not their +Nitsche formulation. Their Section 4.2 assumes body-fitted Neumann boundaries; +its stability/convergence analysis does not cover this implementation. +""" + +import numpy as np +import underworld3 as uw +from mpi4py import MPI +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.spatial import cKDTree + +from ._cylindrical_pressure import _pressure_jump_values + + +def _trace_rule(edges: np.ndarray, sign: int, mode: int, quadrature: int) -> dict: + """Return P2-interpolated nodal radial tests and polygonal line quadrature. + + The weighted test is I_h(n_node*cos(mode*theta_node)); its normalization + is the integral of the continuous cos² over the polygonal wall, as in UW3. + """ + t, weight = np.polynomial.legendre.leggauss(quadrature) + t, weight = (t + 1) / 2, weight / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + midpoint = edges.mean(axis=1) + nodes = np.stack([edges[:, 0], edges[:, 1], midpoint], axis=1) + node_normal = sign * nodes / np.linalg.norm(nodes, axis=2)[..., None] + node_angle = np.arctan2(nodes[..., 1], nodes[..., 0]) + node_cosine = np.cos(mode * node_angle) + node_tangent = sign * np.stack([-node_normal[..., 1], node_normal[..., 0]], axis=-1) + points = (1 - t)[None, :, None] * edges[:, None, 0] + t[None, :, None] * edges[:, None, 1] + return dict( + points=points, + normal=sign * midpoint / np.linalg.norm(midpoint, axis=1)[:, None], + weight=np.linalg.norm(edges[:, 1] - edges[:, 0], axis=1)[:, None] * weight, + harmonic=np.cos(mode * np.arctan2(points[..., 1], points[..., 0])), + test=np.einsum("qi,eid->eqd", shape, node_normal * node_cosine[..., None]), + test_one=np.einsum("qi,eid->eqd", shape, node_normal), + tangent_test=np.einsum( + "qi,eid->eqd", shape, node_tangent * np.sin(mode * node_angle)[..., None] + ), + ) + +def _tangent_load_density(edges: np.ndarray, traction: np.ndarray, sign: int, quadrature: int): + """Project integrated loads tangentially at P2 nodes; preserve their moments. + + Project loads before inverting the mass matrix. Projecting pointwise density + instead would generally inject normal nodal loads into the reaction. + """ + rule = _trace_rule(edges, sign, 2, quadrature) + t = (np.polynomial.legendre.leggauss(quadrature)[0] + 1) / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + points = np.stack((edges[:, 0], edges[:, 1], edges.mean(axis=1)), axis=1) + nodes, index = np.unique(points.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + local_mass = np.einsum("eq,qi,qj->eij", rule["weight"], shape, shape) + local_load = np.einsum("eq,qi,eqd->eid", rule["weight"], shape, traction) + row = np.broadcast_to(index[:, :, None], local_mass.shape).ravel() + column = np.broadcast_to(index[:, None, :], local_mass.shape).ravel() + mass = coo_matrix((local_mass.ravel(), (row, column)), shape=(len(nodes), len(nodes))).tocsr() + load = np.zeros_like(nodes) + np.add.at(load, index.ravel(), local_load.reshape(-1, 2)) + normal = sign * nodes / np.linalg.norm(nodes, axis=1)[:, None] + tangent_load = load - np.sum(load * normal, axis=1)[:, None] * normal + density = spsolve(mass, tangent_load) + reconstructed = mass @ density + checks = dict( + nodes=len(nodes), + edges=len(edges), + maximum_mass_moment_error=float(np.max(np.abs(reconstructed - tangent_load))), + maximum_normal_nodal_load=float(np.max(np.abs(np.sum(reconstructed * normal, axis=1)))), + tangential_load_norm=float(np.linalg.norm(tangent_load)), + ) + if checks["maximum_mass_moment_error"] > 1e-12 or checks["maximum_normal_nodal_load"] > 1e-12: + raise ArithmeticError(f"Tangential boundary-load projection failed: {checks}") + return nodes, density, checks + +def _wall_solve_succeeded(info): + """Accept a converged residual, including a valid zero-increment warm start.""" + return bool( + info["converged"] + and ( + info["ksp_reason"] > 0 + or (info["ksp_reason"] == 0 and info["nonlinear_iterations"] == 0) + ) + ) + +def _solve_with_wall_updates(stokes, walls, updates): + """Alternate converged Stokes solves and numerical wall updates. + + The final solve uses the last recorded boundary fields. A residual-accepted + zero-increment warm start is valid. Partial history is stored on the solver + before raising so the caller can record a failed step without accepting it. + Requires the caller's ``petsc_use_pressure_nullspace=True`` configuration: + accepted history reads the projected pressure-gauge residual. This helper + does not install a pressure gauge or support another gauge configuration. + The fixed-point schedule and choice of four updates are UW3 benchmark + choices; they are not an iteration theorem from the shifted-boundary paper + cited in this module's documentation. + """ + if updates < 0 or int(updates) != updates: + raise ValueError("Wall updates must be a non-negative integer.") + velocity = stokes.u + history = [] + stokes.cylindrical_wall_history = history + for iteration in range(updates + 1): + previous = np.asarray(velocity.array).copy() + stokes.solve() + info = stokes._rotated_freeslip_info + step = dict( + iteration=iteration, + velocity_change_maximum=uw.mpi.comm.allreduce( + float(np.max(np.abs(velocity.array - previous), initial=0.0)), op=MPI.MAX + ), + ksp_iterations=info["ksp_its"], + ksp_reason=int(info["ksp_reason"]), + nonlinear_iterations=int(info["nonlinear_iterations"]), + converged=bool(info["converged"]), + residual_norm=float(info["rnorm"]), + pressure_gauge_residual=float(stokes._rotated_pressure_gauge_residual), + boundary_updates={}, + ) + history.append(step) + if not _wall_solve_succeeded(info): + raise RuntimeError( + f"Numerical wall iteration {iteration} failed: " + f"KSP={info['ksp_reason']}, residual={info['rnorm']}, reference={info['rnorm0']}" + ) + if iteration < updates: + step["boundary_updates"] = {name: wall.update() for name, wall in walls.items()} + uw.pprint( + f"Numerical wall iteration {iteration}: " + f"velocity change={step['velocity_change_maximum']:.3e}" + ) + return history + +def _zero_torque_density(edges, nodes, density): + """Remove the M-weighted rigid-rotation load, preserving radial nodal loads. + + For affine P2 line traces M*w=diag(int(N_i))*w for w=(-y,x). + Subtracting this linear density therefore changes only tangential loads. + The radial sheet body force and true free-slip walls exert zero torque. + """ + edge_nodes = np.stack((edges[:, 0], edges[:, 1], edges.mean(axis=1)), axis=1) + distance, index = cKDTree(nodes).query(edge_nodes) + if distance.max() > 1e-10: + raise ValueError("Boundary density nodes do not match P2 edge nodes.") + length = np.linalg.norm(edges[:, 1] - edges[:, 0], axis=1) + weights = np.zeros(len(nodes)) + np.add.at(weights, index.ravel(), (length[:, None] * np.array([1 / 6, 1 / 6, 2 / 3])).ravel()) + rotation = np.column_stack((-nodes[:, 1], nodes[:, 0])) + torque = float(np.sum(weights[:, None] * rotation * density)) + coefficient = torque / np.sum(weights * np.sum(rotation**2, axis=1)) + balanced = density - coefficient * rotation + return balanced, dict( + raw_boundary_torque=torque, + removed_rotation_density_coefficient=float(coefficient), + corrected_boundary_torque=float(np.sum(weights[:, None] * rotation * balanced)), + ) + +def _numerical_trace( + vertices, velocity, pressure, points, viscosity, sheet_radius, outer_radius, mode, amplitude +): + """Evaluate physical stress and radial derivatives from affine-cell FE data. + + Velocity nodes are vertices 012 then edge midpoints 01,12,20; pressure nodes + are the three vertices. Spatial derivatives are Cartesian and one-sided in + each boundary-adjacent triangle. Radial differentiation holds angle fixed. + """ + inverse = np.linalg.inv( + np.stack((vertices[:, 1] - vertices[:, 0], vertices[:, 2] - vertices[:, 0]), axis=-1) + ) + xi = np.einsum("eij,eqj->eqi", inverse, points - vertices[:, None, 0]) + bary = np.concatenate((1 - xi.sum(axis=-1, keepdims=True), xi), axis=-1) + grad = np.stack((-inverse.sum(axis=1), inverse[:, 0], inverse[:, 1]), axis=1) + shapes_grad = np.stack( + [(4 * bary[..., i, None] - 1) * grad[:, None, i] for i in range(3)] + + [ + 4 * (bary[..., i, None] * grad[:, None, j] + bary[..., j, None] * grad[:, None, i]) + for i, j in ((0, 1), (1, 2), (2, 0)) + ], + axis=2, + ) + hessian = np.stack( + [4 * np.einsum("ei,ej->eij", grad[:, i], grad[:, i]) for i in range(3)] + + [ + 4 + * ( + np.einsum("ei,ej->eij", grad[:, i], grad[:, j]) + + np.einsum("ei,ej->eij", grad[:, j], grad[:, i]) + ) + for i, j in ((0, 1), (1, 2), (2, 0)) + ], + axis=1, + ) + velocity_grad = np.einsum("eni,eqnj->eqij", velocity, shapes_grad) + velocity_hessian = np.einsum("eni,enjk->eijk", velocity, hessian) + radius = np.linalg.norm(points, axis=-1) + radial = points / radius[..., None] + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + lifting = np.zeros_like(radius) + lifting_grad = np.zeros_like(points) + if amplitude != 0: + lifting, lifting_grad = _pressure_jump_values( + points, sheet_radius, outer_radius, mode, amplitude) + physical_pressure = np.einsum("en,eqn->eq", pressure, bary) + lifting + pressure_grad = np.einsum("en,enj->ej", pressure, grad)[:, None, :] + lifting_grad + derivative = np.einsum("eqi,eqij,eqj->eq", radial, velocity_grad, radial) + rr = 2 * viscosity * derivative - physical_pressure + tt = ( + 2 * viscosity * np.einsum("eqi,eqij,eqj->eq", tangent, velocity_grad, tangent) + - physical_pressure + ) + rt = viscosity * ( + np.einsum("eqi,eqij,eqj->eq", radial, velocity_grad, tangent) + + np.einsum("eqi,eqij,eqj->eq", tangent, velocity_grad, radial) + ) + rt_derivative = viscosity * ( + np.einsum("eqi,eijk,eqj,eqk->eq", radial, velocity_hessian, tangent, radial) + + np.einsum("eqi,eijk,eqj,eqk->eq", tangent, velocity_hessian, radial, radial) + ) + rr_derivative = 2 * viscosity * np.einsum( + "eqi,eijk,eqj,eqk->eq", radial, velocity_hessian, radial, radial + ) - np.sum(radial * pressure_grad, axis=-1) + return dict( + radial_velocity_derivative=derivative, + rr=rr, + tt=tt, + rt=rt, + rt_derivative=rt_derivative, + rr_derivative=rr_derivative, + ) + +class _NumericalWall: + """Cache one wall's numerical boundary traces and lagged Taylor data. + + First-order transfer follows the Taylor identity in Atallah, Canuto and + Scovazzi (2020), Section 2.2; see the full reference in this module. + Strong normal constraints and compatible tangential loads are UW3's + formulation, with no claim to inherit their Nitsche stability analysis. + """ + + def __init__( + self, + mesh, + velocity, + pressure, + radius, + sign, + viscosity, + sheet_radius, + outer_radius, + mode, + amplitude, + name, + quadrature, + ): + if mesh.dim != 2 or mesh.cdim != 2 or velocity.degree != 2 or not velocity.continuous: + raise ValueError("Wall transfer requires a two-dimensional continuous affine P2 velocity.") + if pressure.degree != 1 or not pressure.continuous: + raise ValueError("Wall transfer requires continuous affine P1 regular pressure.") + if sign not in (-1, 1) or radius <= 0 or viscosity <= 0 or quadrature < 2: + raise ValueError("Require positive radius/viscosity, sign +/-1 and quadrature >= 2.") + self.mesh, self.velocity, self.pressure = mesh, velocity, pressure + self.radius, self.sign, self.quadrature = radius, sign, quadrature + self.name = name + self.physics = (viscosity, sheet_radius, outer_radius, mode, amplitude) + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, 3) + vertices = np.asarray(mesh._get_coords_for_basis(1, True))[vertex_rows] + owned = mesh._get_owned_cells_mask() + cells, edges = [], [] + for first, second in ((0, 1), (1, 2), (2, 0)): + chosen = np.flatnonzero( + owned + & np.all( + np.isclose( + np.linalg.norm(vertices[:, [first, second]], axis=2), + radius, + rtol=0, + atol=1e-10, + ), + axis=1, + ) + ) + cells.extend(chosen) + edges.extend(vertices[chosen][:, [first, second]]) + self.cells = np.array(cells, dtype=int) + self.p_rows = vertex_rows[self.cells] + self.v_rows = mesh._cell_node_indices(2, True).reshape(-1, 6)[self.cells] + local_vertices = vertices[self.cells] + standard = np.concatenate( + (local_vertices, (local_vertices[:, [0, 1, 2]] + local_vertices[:, [1, 2, 0]]) / 2), + axis=1, + ) + distance = np.linalg.norm( + np.asarray(velocity.coords)[self.v_rows][:, :, None, :] - standard[:, None, :, :], + axis=-1, + ) + self.v_order = distance.argmin(axis=1) + if distance.size and np.max(distance.min(axis=1)) > 1e-10: + raise ValueError("Wall transfer requires affine P2 velocity nodes.") + np.testing.assert_allclose( + np.asarray(pressure.coords)[self.p_rows], local_vertices, rtol=0, atol=1e-10 + ) + self.vertices = np.concatenate(uw.mpi.comm.allgather(local_vertices)) + self.edges = np.concatenate(uw.mpi.comm.allgather(np.asarray(edges).reshape(-1, 2, 2))) + self.rule = _trace_rule(self.edges, sign, mode, quadrature) + self.edge_nodes = np.stack( + (self.edges[:, 0], self.edges[:, 1], self.edges.mean(axis=1)), axis=1 + ) + self.nodes, self.node_index = np.unique( + self.edge_nodes.reshape(-1, 2), axis=0, return_inverse=True + ) + self.node_index = self.node_index.reshape(-1, 3) + self.normal_field = uw.discretisation.MeshVariable( + "NormalShift" + name, mesh, 1, degree=2, continuous=True + ) + self.traction_field = uw.discretisation.MeshVariable( + "TangentShift" + name, mesh, 2, degree=2, continuous=True + ) + np.testing.assert_array_equal(self.normal_field.coords, self.traction_field.coords) + distance, index = cKDTree(self.nodes).query(np.asarray(self.normal_field.coords)) + self.local_boundary, self.local_index = distance < 1e-10, index + self.normal_values = np.zeros(len(self.nodes)) + self.density = np.zeros_like(self.nodes) + + def cell_values(self): + """Gather only owned wall-cell traces in the cached geometry order.""" + local_velocity = np.take_along_axis( + np.asarray(self.velocity.array)[:, 0, :][self.v_rows], self.v_order[..., None], axis=1 + ) + velocity = np.concatenate(uw.mpi.comm.allgather(local_velocity)) + pressure = np.concatenate( + uw.mpi.comm.allgather(np.asarray(self.pressure.array)[self.p_rows, 0, 0]) + ) + return velocity, pressure + + def update(self): + """Rebuild Taylor data and project normal flux and tangential torque. + + The first-order boundary shift is motivated by Atallah, Canuto and + Scovazzi (2020), Section 2.2, equations (8)-(13). Nodal load projection, + zero-flux/zero-torque corrections and lagged solves are UW3 additions. + """ + velocity, pressure = self.cell_values() + points = self.rule["points"] + trace = _numerical_trace(self.vertices, velocity, pressure, points, *self.physics) + radius = np.linalg.norm(points, axis=-1) + radial = points / radius[..., None] + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + shear = -(self.radius - radius) * trace["rt_derivative"] + nr = np.sum(self.rule["normal"][:, None, :] * radial, axis=-1) + nt = np.sum(self.rule["normal"][:, None, :] * tangent, axis=-1) + traction = (trace["rr"] * nr + shear * nt)[..., None] * radial + ( + shear * nr + trace["tt"] * nt + )[..., None] * tangent + nodes, density, checks = _tangent_load_density( + self.edges, traction, self.sign, self.quadrature + ) + np.testing.assert_array_equal(nodes, self.nodes) + density, torque_checks = _zero_torque_density(self.edges, nodes, density) + checks.update(torque_checks) + node_trace = _numerical_trace( + self.vertices, velocity, pressure, self.edge_nodes, *self.physics + ) + datum = ( + -self.sign + * (self.radius - np.linalg.norm(self.edge_nodes, axis=-1)) + * node_trace["radial_velocity_derivative"] + ) + length = np.linalg.norm(self.edges[:, 1] - self.edges[:, 0], axis=1) + weights = np.zeros(len(nodes)) + normal = np.zeros(len(nodes)) + np.add.at(normal, self.node_index.ravel(), (length[:, None] * datum).ravel()) + np.add.at( + weights, self.node_index.ravel(), np.broadcast_to(length[:, None], datum.shape).ravel() + ) + normal /= weights + # Integrate I_h(g_i*n_i).n_h exactly, rather than averaging nodal values. + t = (np.polynomial.legendre.leggauss(self.quadrature)[0] + 1) / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + normals = self.sign * self.edge_nodes / np.linalg.norm(self.edge_nodes, axis=-1)[..., None] + flux_weight = np.einsum( + "eq,qi,eid,ed->ei", self.rule["weight"], shape, normals, self.rule["normal"] + ) + raw_flux = np.sum(flux_weight * normal[self.node_index]) + removed_mean = raw_flux / np.sum(flux_weight) + normal -= removed_mean + checks.update( + normal_datum_change=float(np.max(np.abs(normal - self.normal_values))), + traction_density_change=float(np.max(np.abs(density - self.density))), + normal_datum_maximum=float(np.max(np.abs(normal))), + raw_normal_flux=float(raw_flux), + removed_normal_datum_mean=float(removed_mean), + corrected_normal_flux=float(np.sum(flux_weight * normal[self.node_index])), + ) + self.normal_values, self.density = normal, density + selected, index = self.local_boundary, self.local_index[self.local_boundary] + self.normal_field.array[selected, 0, 0] = normal[index] + self.traction_field.array[selected, 0, :] = density[index] + return checks + + +def _add_cylindrical_wall_transfer(stokes, *, boundary, radius, sign, viscosity, + interface_radius, outer_radius, wavenumber, + amplitude, quadrature): + wall = _NumericalWall(stokes.mesh, stokes.u, stokes.p, radius, sign, + viscosity, interface_radius, outer_radius, wavenumber, + amplitude, boundary, quadrature) + radial = stokes.mesh.CoordinateSystem.unit_e_0 + stokes.add_rotated_freeslip_bc(wall.normal_field.sym[0], boundary, + normal=sign * radial) + stokes.add_natural_bc(-wall.traction_field.sym, boundary) + return wall diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 33b5bdad6..4fd412a90 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -368,6 +368,26 @@ def _eval_datum(dspec, coords, solver): return np.asarray(uw.function.evaluate(dspec, np.ascontiguousarray(coords))).flatten() +def _datum_batches(boundaries, datum_specs, datum_nodes, coordinates, dim): + """Schedule datum evaluation by registered boundary, including empty ranks. + + FE datum evaluation updates a mesh local vector collectively. Scheduling + by locally owned nodes (or local expression identities) can give ranks + different call counts/orders and deadlock before the first Stokes solve. + Registration order is common to all ranks; coordinate arrays may be empty. + """ + if datum_specs is None: + return + for spec in boundaries: + name, _ = _boundary_spec(spec) + dspec = datum_specs.get(name) + if dspec is None: + continue + items = [(row, q, sign) for row, q, boundary, sign in datum_nodes if boundary == name] + coords = np.asarray([coordinates[q] for _, q, _ in items], dtype=float).reshape(-1, dim) + yield dspec, coords, items + + def build_rotation(solver, boundaries, datum_specs=None): """Global sparse rotation Q on the composite saddle vector: identity except a per-node (normal,tangential) block at each velocity node of `boundaries`. @@ -394,7 +414,7 @@ def build_rotation(solver, boundaries, datum_specs=None): for q, nrm in _boundary_velocity_nodes(solver, name, normal=normal): node_normals.setdefault(q, []).append(nrm) if dspec is not None: - node_dspec.setdefault(q, dspec) + node_dspec.setdefault(q, name) # Split-fault contact pairs (fault_contact.add_frictionless_fault_bc): # their 2·dim mean/jump blocks widen the widest row of Q. @@ -559,19 +579,17 @@ def build_rotation(solver, boundaries, datum_specs=None): Q.assemble() Qt = Q.transpose(PETSc.Mat()) - # Evaluate the prescribed normal datum at each simple constrained node, batched per - # distinct datum spec (usually one — the surface). Rotated row grows[0] == v_n. + # Every rank evaluates every registered datum, even with no owned nodes. + # FE expressions contain collectives; local-only batching can deadlock. datum_map = {} - if datum_nodes: + if datum_specs is not None: csec = dm.getCoordinateSection() cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) v0, v1 = dm.getDepthStratum(0) - by_spec = {} - for grow0, q, dspec, sgn in datum_nodes: - by_spec.setdefault(id(dspec), (dspec, []))[1].append((grow0, q, sgn)) - for dspec, items in by_spec.values(): - coords = np.array([_point_coord(dm, dim, cvec, csec, v0, v1, q) - for _, q, _ in items]) + coordinates = {q: _point_coord(dm, dim, cvec, csec, v0, v1, q) + for _, q, _, _ in datum_nodes} + for dspec, coords, items in _datum_batches( + boundaries, datum_specs, datum_nodes, coordinates, dim): vals = _eval_datum(dspec, coords, solver) for (grow0, _, sgn), val in zip(items, vals): # A zero datum IS free-slip; leave those rows out so an all-zero datum diff --git a/tests/parallel/test_1084_rotated_datum_collectives.py b/tests/parallel/test_1084_rotated_datum_collectives.py new file mode 100644 index 000000000..a0191790f --- /dev/null +++ b/tests/parallel/test_1084_rotated_datum_collectives.py @@ -0,0 +1,40 @@ +"""Prescribed FE wall data must be evaluated in the same order on empty ranks. + +Run serial or with mpiexec. This exercises the datum schedule/evaluator only; +no Stokes model or spherical mesh is run. +""" + +import numpy as np +import underworld3 as uw +from mpi4py import MPI + +from underworld3.utilities.rotated_bc import _datum_batches, _eval_datum + + +def test_registered_boundaries_determine_collective_schedule_not_local_nodes(): + rank = uw.mpi.rank + boundaries = ["Upper", "Lower"] + # All ranks have the same registration, but zero, one or two local groups. + owned = {0: [(10, 1, "Upper", 1.0), (11, 2, "Lower", -1.0)], + 1: [(20, 2, "Lower", 1.0)]}.get(rank, []) + coordinates = {1: [2.0, 0.0], 2: [1.0, 0.0]} + batches = list(_datum_batches(boundaries, {"Upper": 2.0, "Lower": 3.0}, owned, coordinates, 2)) + assert len(batches) == 2 + for index, (datum, coords, items) in enumerate(batches): + assert uw.mpi.comm.allgather(datum) == [2.0+index]*uw.mpi.size + assert coords.shape == (len(items), 2) + if rank > 1: + assert coords.shape == (0, 2) + + +def test_fe_datum_evaluation_accepts_empty_coordinate_batches_on_other_ranks(): + mesh = uw.meshing.Annulus(radiusInner=1.22, radiusOuter=2.22, cellSize=1/32, qdegree=2) + field = uw.discretisation.MeshVariable("DatumCollectiveP2", mesh, 1, degree=2, continuous=True) + field.data[:, 0] = 0.13*field.coords[:, 0]-0.2*field.coords[:, 1] + # Only one rank queries nodes; every rank must still enter FE evaluation. + coords = np.asarray(field.coords[:3]) if uw.mpi.rank == 0 else np.empty((0, 2)) + for _ in range(2): + actual = _eval_datum(field.sym[0], coords, None) + expected = 0.13*coords[:, 0]-0.2*coords[:, 1] + np.testing.assert_allclose(actual, expected, atol=1e-12) + assert uw.mpi.comm.allreduce(1, op=MPI.SUM) == uw.mpi.size diff --git a/tests/test_1089_cylindrical_geometry.py b/tests/test_1089_cylindrical_geometry.py new file mode 100644 index 000000000..d65cb2a4a --- /dev/null +++ b/tests/test_1089_cylindrical_geometry.py @@ -0,0 +1,334 @@ +"""Hard mathematical baselines for affine annulus geometry corrections. + +No Stokes model is run. Independent polar integration verifies circular cuts; +symbolic polynomial stresses verify derivatives; exact line masses verify +reaction and compatibility moments. Run serial or on multiple MPI ranks. +""" + +from types import SimpleNamespace + +import numpy as np +import pytest +import sympy as sp + +from underworld3.utilities._cylindrical_pressure import ( + _extended_regular_force as extended_regular_force, + _load_moments as load_moments, + _p2_shapes as p2_shapes, + _profile_values as profile_values, +) +from underworld3.utilities._cylindrical_wall import ( + _numerical_trace as numerical_trace, + _solve_with_wall_updates, + _wall_solve_succeeded as wall_solve_succeeded, + _zero_torque_density as zero_torque_density, +) +from underworld3.utilities._cylindrical_reactions import ( + _circular_reaction_field as circular_reaction_field, +) + +pytestmark = pytest.mark.level_1 + + +@pytest.mark.tier_b +def test_profile_has_prescribed_jump_and_vanishes_at_outer_radius(): + value, derivative = profile_values(np.array([2.0, 2.3, 2.4]), 2.0, 0.3) + np.testing.assert_allclose(value, [1, 0, -1 / 3], atol=1e-14) + np.testing.assert_allclose(derivative, -1 / 0.3, atol=1e-14) + + +@pytest.mark.tier_b +def test_regular_force_is_positive_gradient_of_extracted_profile(): + points = np.array([[2.04, 0.03], [1.98, 0.5], [-2.05, 0.1]]) + actual = extended_regular_force(points, 2.0, 0.3, 2, 1.7) + numerical = np.empty_like(actual) + for component in range(2): + displacement = np.eye(2)[component] * 1e-6 + values = [] + for shifted in (points + displacement, points - displacement): + radius = np.linalg.norm(shifted, axis=1) + profile, _ = profile_values(radius, 2.0, 0.3) + values.append(1.7 * profile * np.cos(2 * np.arctan2(shifted[:, 1], shifted[:, 0]))) + numerical[:, component] = (values[0] - values[1]) / 2e-6 + np.testing.assert_allclose(actual, numerical, rtol=1e-8, atol=1e-8) + + +@pytest.mark.tier_b +def test_triangle_load_against_independent_polar_domain_integration(): + radius, alpha = 2.0, 0.05 + vertices = np.array( + [ + [radius * np.cos(alpha), -radius * np.sin(alpha)], + [radius * np.cos(alpha), radius * np.sin(alpha)], + [2.15, 0], + ] + ) + moments, cap_area, chords, _ = load_moments(vertices[None], radius, 0.4, 2, 1.0, 12) + edges = np.column_stack((vertices[1] - vertices[0], vertices[2] - vertices[0])) + inverse = np.linalg.inv(edges) + nodes, weights = np.polynomial.legendre.leggauss(24) + nodes, weights = (nodes + 1) / 2, weights / 2 + independent = np.zeros((6, 2)) + pressure_divergence = 0.0 + for sign in (-1, 1): + endpoint = vertices[0 if sign < 0 else 1] + edge = vertices[2] - endpoint + determinant = endpoint[0] * vertices[2, 1] - endpoint[1] * vertices[2, 0] + for t, wt in zip(nodes, weights): + theta = sign * alpha * t + direction = np.array([np.cos(theta), np.sin(theta)]) + upper = determinant / (direction[0] * edge[1] - direction[1] * edge[0]) + for s, ws in zip(nodes, weights): + r = radius + s * (upper - radius) + point = r * direction + xi = inverse @ (point - vertices[0]) + bary = np.r_[1 - xi.sum(), xi] + force = extended_regular_force(point, radius, 0.4, 2, 1.0) + independent += ( + r * (upper - radius) * alpha * wt * ws * p2_shapes(bary)[:, None] * force + ) + profile, _ = profile_values(r, radius, 0.4) + pressure_divergence += ( + r * (upper - radius) * alpha * wt * ws * profile * np.cos(2 * theta) + ) + np.testing.assert_allclose(moments[0], independent, rtol=1e-9, atol=1e-12) + np.testing.assert_allclose(cap_area, radius**2 * (alpha - np.sin(2 * alpha) / 2), rtol=1e-10) + np.testing.assert_array_equal(chords, [1]) + # Independently verify the weak sign and boundary term using v=(x,0). + boundary_flux = 0.0 + for first, second in ((1, 2), (2, 0)): + edge = vertices[second] - vertices[first] + normal_ds = np.array([-edge[1], edge[0]]) # clockwise triangle + for t, wt in zip(nodes, weights): + point = vertices[first] + t * edge + profile, _ = profile_values(np.linalg.norm(point), radius, 0.4) + boundary_flux += ( + wt * profile * np.cos(2 * np.arctan2(point[1], point[0])) * point[0] * normal_ds[0] + ) + for t, wt in zip(nodes, weights): + theta = alpha * (2 * t - 1) + boundary_flux -= 2 * alpha * wt * radius**2 * np.cos(2 * theta) * np.cos(theta) ** 2 + x_dofs = np.r_[vertices[:, 0], (vertices[[0, 1, 2], 0] + vertices[[1, 2, 0], 0]) / 2] + np.testing.assert_allclose( + moments[0, :, 0] @ x_dofs, boundary_flux - pressure_divergence, rtol=1e-9, atol=1e-12 + ) + + +@pytest.mark.tier_b +def test_converged_zero_increment_warm_start_is_not_a_ksp_failure(): + assert wall_solve_succeeded(dict(converged=True, ksp_reason=0, nonlinear_iterations=0)) + assert wall_solve_succeeded(dict(converged=True, ksp_reason=2, nonlinear_iterations=1)) + assert not wall_solve_succeeded(dict(converged=False, ksp_reason=0, nonlinear_iterations=0)) + assert not wall_solve_succeeded(dict(converged=True, ksp_reason=0, nonlinear_iterations=1)) + assert not wall_solve_succeeded(dict(converged=True, ksp_reason=-3, nonlinear_iterations=1)) + + +@pytest.mark.tier_b +def test_p2_p1_wall_derivatives_against_symbolic_physical_stress(): + x, y = sp.symbols("x y", real=True) + radius = sp.sqrt(x * x + y * y) + radial = sp.Matrix([x, y]) / radius + tangent = sp.Matrix([-y, x]) / radius + velocity = sp.Matrix([x * x + 2 * x * y - y * y + 3 * x, -x * x - 2 * x * y - y * y - 3 * y]) + regular_pressure = 1 + 2 * x - sp.Rational(7, 10) * y + viscosity = 1.3 + for offset in (1.1, 2.02): + vertices = np.array([[[offset, 0.02], [offset + 0.12, 0.02], [offset + 0.04, 0.16]]]) + nodes = np.concatenate( + (vertices, (vertices[:, [0, 1, 2]] + vertices[:, [1, 2, 0]]) / 2), axis=1 + ) + v = np.stack( + [ + sp.lambdify((x, y), value, "numpy")(nodes[..., 0], nodes[..., 1]) + for value in velocity + ], + axis=-1, + ) + p = sp.lambdify((x, y), regular_pressure, "numpy")(vertices[..., 0], vertices[..., 1]) + bary = np.array([[0.2, 0.3, 0.5], [0.4, 0.2, 0.4]]) + points = np.einsum("qi,eid->eqd", bary, vertices) + result = numerical_trace(vertices, v, p, points, viscosity, 2.0, 2.22, 2, 0.8) + jump = ( + -0.8 * (2.22 - radius) / 0.22 * (x * x - y * y) / (x * x + y * y) if offset > 2 else 0 + ) + gradient = velocity.jacobian([x, y]) + stress = viscosity * (gradient + gradient.T) - (regular_pressure + jump) * sp.eye(2) + rr = (radial.T * stress * radial)[0] + rt = (radial.T * stress * tangent)[0] + expressions = dict( + radial_velocity_derivative=(radial.T * gradient * radial)[0], + rr=rr, + tt=(tangent.T * stress * tangent)[0], + rt=rt, + rt_derivative=radial[0] * sp.diff(rt, x) + radial[1] * sp.diff(rt, y), + rr_derivative=radial[0] * sp.diff(rr, x) + radial[1] * sp.diff(rr, y), + ) + for key, expression in expressions.items(): + expected = sp.lambdify((x, y), expression, "numpy")(points[..., 0], points[..., 1]) + np.testing.assert_allclose(result[key], expected, rtol=1e-11, atol=1e-11, err_msg=key) + + +@pytest.mark.tier_b +def test_rigid_rotation_has_no_stress_or_radial_velocity_transfer(): + vertices = np.array([[[1.1, 0.1], [1.2, 0.1], [1.1, 0.2]]]) + nodes = np.concatenate( + (vertices, (vertices[:, [0, 1, 2]] + vertices[:, [1, 2, 0]]) / 2), axis=1 + ) + velocity = 2.3 * np.stack((-nodes[..., 1], nodes[..., 0]), axis=-1) + points = np.einsum("qi,eid->eqd", np.array([[0.2, 0.3, 0.5], [0.1, 0.8, 0.1]]), vertices) + result = numerical_trace(vertices, velocity, np.zeros((1, 3)), points, 1.0, 2.0, 2.22, 2, 0.0) + for value in result.values(): + np.testing.assert_allclose(value, 0, atol=1e-11) + + +@pytest.mark.tier_b +def test_smooth_pressure_needs_no_sheet_radius_or_jump(): + vertices = np.array([[[1.1, 0.1], [1.2, 0.1], [1.1, 0.2]]]) + points = vertices.mean(axis=1)[:, None, :] + trace = numerical_trace( + vertices, + np.zeros((1, 6, 2)), + np.ones((1, 3)) * 2.3, + points, + 1.0, + 2.22, + 2.22, + 2, + 0.0, + ) + np.testing.assert_allclose(trace["rr"], -2.3) + np.testing.assert_allclose(trace["tt"], -2.3) + np.testing.assert_allclose(trace["rr_derivative"], 0, atol=1e-13) + + +@pytest.mark.tier_b +def test_torque_projection_preserves_normal_nodal_loads_on_irregular_trace(): + angles = np.linspace(0, 2 * np.pi, 41)[:-1] + angles += 0.025 * np.sin(3 * angles) + vertices = 1.7 * np.column_stack((np.cos(angles), np.sin(angles))) + edges = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + edge_nodes = np.stack((edges[:, 0], edges[:, 1], edges.mean(axis=1)), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + rotation = np.column_stack((-nodes[:, 1], nodes[:, 0])) + density = 0.4 * rotation + nodes * np.sin(2 * np.arctan2(nodes[:, 1], nodes[:, 0]))[:, None] + balanced, checks = zero_torque_density(edges, nodes, density) + loads = [] + for field in (density, balanced): + value = np.zeros_like(nodes) + for edge, rows in zip(edges, index): + mass = ( + np.linalg.norm(edge[1] - edge[0]) + / 30 + * np.array([[4, -1, 2], [-1, 4, 2], [2, 2, 16]]) + ) + np.add.at(value, rows, mass @ field[rows]) + loads.append(value) + np.testing.assert_allclose( + np.sum(loads[0] * rotation), checks["raw_boundary_torque"], atol=1e-13 + ) + np.testing.assert_allclose(np.sum(loads[1] * rotation), 0, atol=1e-13) + normal = nodes / np.linalg.norm(nodes, axis=1)[:, None] + np.testing.assert_allclose(np.sum((loads[1] - loads[0]) * normal, axis=1), 0, atol=1e-13) + + +@pytest.mark.tier_b +def test_spatial_reaction_transfer_preserves_cbf_residual_and_removes_pressure_gauge(): + radius, count = 1.7, 48 + theta = np.arange(count) * 2 * np.pi / count + theta += 0.015 * np.sin(3 * theta) + vertices = radius * np.column_stack((np.cos(theta), np.sin(theta))) + edges = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + edge_nodes = np.stack((edges[:, 0], edges[:, 1], edges.mean(axis=1)), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + z, w = np.polynomial.legendre.leggauss(16) + t, w = (z + 1) / 2, w / 2 + shape = np.array([(1 - t) * (1 - 2 * t), t * (2 * t - 1), 4 * t * (1 - t)]).T + points = np.einsum("qi,eid->eqd", shape, edge_nodes) + r = np.linalg.norm(points, axis=-1) + harmonic = (points[..., 0] ** 2 - points[..., 1] ** 2) / r**2 + normal = edges.mean(axis=1) / np.linalg.norm(edges.mean(axis=1), axis=1)[:, None] + length = np.linalg.norm(edges[:, 1] - edges[:, 0], axis=1) + polygon_weights = length[:, None] * w + circle_weights = ( + polygon_weights + * radius + / r + * np.abs(np.sum(normal[:, None, :] * points / r[..., None], axis=-1)) + ) + masses = [] + for weights in (polygon_weights, circle_weights): + mass = np.zeros((len(nodes), len(nodes))) + for e, rows in enumerate(index): + mass[np.ix_(rows, rows)] += np.einsum("q,qi,qj->ij", weights[e], shape, shape) + masses.append(mass) + for sign in (-1, 1): + for gauge in (0, 2.3): + pressure = r * harmonic + gauge + normal_nodes = sign * edge_nodes / np.linalg.norm(edge_nodes, axis=-1)[..., None] + traction = -pressure[..., None] * sign * normal[:, None, :] + loads = np.zeros(len(nodes)) + for e, rows in enumerate(index): + load = -np.einsum( + "q,qi,qd,id->i", polygon_weights[e], shape, traction[e], normal_nodes[e] + ) + np.add.at(loads, rows, load) + raw = np.linalg.solve(masses[0], loads) + trace = dict(rr=-pressure, tt=-pressure, rt=np.zeros_like(r), rr_derivative=-harmonic) + _, corrected, checks = circular_reaction_field(edges, trace, raw, radius, sign, 16) + expected = np.zeros(len(nodes)) + for e, rows in enumerate(index): + np.add.at( + expected, + rows, + np.einsum("q,qi,q->i", circle_weights[e], shape, radius * harmonic[e]), + ) + expected -= expected.sum() / masses[1].sum() * (masses[1] @ np.ones(len(nodes))) + np.testing.assert_allclose(masses[1] @ corrected, expected, atol=1e-13) + # A numerical CBF residual must survive the geometric correction. + residual = np.cos(5 * np.arctan2(nodes[:, 1], nodes[:, 0])) + residual -= (masses[0] @ residual).sum() / masses[0].sum() + _, changed, _ = circular_reaction_field(edges, trace, raw + residual, radius, sign, 16) + np.testing.assert_allclose( + masses[1] @ (changed - corrected), masses[0] @ residual, atol=1e-13 + ) + assert checks["maximum_mass_moment_error"] < 1e-13 + + +@pytest.mark.tier_b +def test_wall_update_sequence_accepts_zero_increment_and_retains_failure(): + class Wall: + count = 0 + + def update(self): + self.count += 1 + return {"update": self.count} + + class Stokes: + count = 0 + u = SimpleNamespace(array=np.zeros((0, 1, 2))) + _rotated_pressure_gauge_residual = 0.0 + + def solve(self): + self.count += 1 + self._rotated_freeslip_info = dict( + converged=True, ksp_reason=2 if self.count == 1 else 0, + nonlinear_iterations=1 if self.count == 1 else 0, + ksp_its=[3] if self.count == 1 else [], rnorm=1e-13, rnorm0=1.0) + + wall, stokes = Wall(), Stokes() + history = _solve_with_wall_updates(stokes, {"Upper": wall}, 4) + assert stokes.count == 5 and wall.count == 4 + assert len(history) == 5 and history[-1]["boundary_updates"] == {} + assert history[-1]["ksp_reason"] == 0 + assert history[-1]["velocity_change_maximum"] == 0 + failed_info = dict(converged=False, ksp_reason=-3, nonlinear_iterations=1, + ksp_its=[9], rnorm=0.02, rnorm0=1.0) + failed = SimpleNamespace(solve=lambda: None, u=stokes.u, + _rotated_freeslip_info=failed_info, + _rotated_pressure_gauge_residual=0.0) + with pytest.raises(RuntimeError, match="iteration 0 failed"): + _solve_with_wall_updates(failed, {}, 4) + assert len(failed.cylindrical_wall_history) == 1 + assert failed.cylindrical_wall_history[0]["ksp_reason"] == -3 From 80243e81065b31abf531c6ceb25fa02dee7eca55 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Fri, 23 Oct 2026 20:55:28 +1100 Subject: [PATCH 2/6] docs: state nondimensional cylindrical correction contract Constrain the public examples and helper docstrings to nondimensional inputs; mark prerequisite datum tests level_1/tier_b. No numerical implementation or test-body changes. Validation: Gadi182128805 built bf81ed97 against the existing int32 image, verified native/Python origins and passed all 12 tests serially and on every one of 16 MPI ranks. Executable ASTs remain identical after these documentation/selection changes; coverage355, style and diff gates pass. No native/model reruns are affected. Next: narrow upstream review, retaining separate integration-campaign provenance. Underworld development team with AI support from Claude Code --- .../advanced/cylindrical-geometry-correction.md | 4 +++- .../cython/petsc_generic_snes_solvers.pyx | 17 +++++++++++++++++ .../utilities/_cylindrical_pressure.py | 3 +++ .../utilities/_cylindrical_reactions.py | 3 +++ src/underworld3/utilities/_cylindrical_wall.py | 3 +++ .../test_1084_rotated_datum_collectives.py | 5 +++++ 6 files changed, 34 insertions(+), 1 deletion(-) diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md index 966838f55..5929484bb 100644 --- a/docs/advanced/cylindrical-geometry-correction.md +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -2,7 +2,9 @@ These explicitly selected Stokes methods support two-dimensional, affine triangular annuli with continuous P2 velocity, continuous P1 regular pressure -and constant viscosity. They correct source integration, wall data and reaction +and constant viscosity. Use nondimensional mesh coordinates, fields and input +parameters throughout; active reference-unit scaling is unsupported. +They correct source integration, wall data and reaction recovery; they do not replace straight elements with curved elements. Numerical validation currently targets the isoviscous annulus benchmarks. They do not establish support for spherical shells, variable viscosity or other FE pairs. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 23065e675..fbf1a352e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8678,6 +8678,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Notes ----- + Only nondimensional mesh coordinates, fields and input scalars are + supported; active reference-unit scaling is unsupported. + Explicitly selected correction for affine 2D P2/continuous-P1 annuli, constant viscosity and a fitted circular sheet. No analytical Stokes solution enters construction. Call collectively on every rank. @@ -8729,6 +8732,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Notes ----- + Only nondimensional mesh coordinates, fields and input scalars are + supported; active reference-unit scaling is unsupported. + Restricted to affine 2D P2/continuous-P1 annuli and constant viscosity. Registers rotated normal and tangential natural data. Only owned boundary-cell traces are replicated; no volume field is gathered. @@ -8772,6 +8778,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Notes ----- + Only nondimensional mesh coordinates, fields and input scalars are + supported; active reference-unit scaling is unsupported. + Before calling, set ``stokes.petsc_use_pressure_nullspace = True``. The accepted-solve history reads the rotated solver's projected pressure-gauge residual, which is recorded only with this projection @@ -8808,6 +8817,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): assembled CBF reaction cache remains unchanged. No coefficients or reference fields enter this transfer. + Notes + ----- + Only nondimensional mesh coordinates, fields and input scalars are + supported; active reference-unit scaling is unsupported. + References ---------- S. Zhong, M. Gurnis and G. Hulbert (1993), "Accurate determination @@ -8842,6 +8856,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Notes ----- + Only nondimensional mesh coordinates, fields and input scalars are + supported; active reference-unit scaling is unsupported. + This circle-weighted integration of an existing spatial field is a UW3 postprocessing construction. The underlying CBF provenance is Zhong, Gurnis and Hulbert (1993), "Accurate determination of surface normal diff --git a/src/underworld3/utilities/_cylindrical_pressure.py b/src/underworld3/utilities/_cylindrical_pressure.py index 212c97e24..9c5b31435 100644 --- a/src/underworld3/utilities/_cylindrical_pressure.py +++ b/src/underworld3/utilities/_cylindrical_pressure.py @@ -10,6 +10,9 @@ radial lift, exact-circle cap integration, DG P2 moment representation and outer-chord natural datum are this UW3 implementation. It does not implement the paper's Section 3 harmonic extension or Section 4 CR discretisation. + +Only nondimensional mesh coordinates, fields and input scalars are supported. +Active reference-unit scaling is outside this helper's supported contract. """ from typing import NamedTuple diff --git a/src/underworld3/utilities/_cylindrical_reactions.py b/src/underworld3/utilities/_cylindrical_reactions.py index 0ec63323d..8b1bb54e9 100644 --- a/src/underworld3/utilities/_cylindrical_reactions.py +++ b/src/underworld3/utilities/_cylindrical_reactions.py @@ -15,6 +15,9 @@ of surface normal stress in viscous flow from a consistent boundary flux method. Physics of the Earth and Planetary Interiors, 78(1-2), 1-8. https://doi.org/10.1016/0031-9201(93)90078-N + +Only nondimensional mesh coordinates, fields and input scalars are supported. +Active reference-unit scaling is outside this helper's supported contract. """ import numpy as np diff --git a/src/underworld3/utilities/_cylindrical_wall.py b/src/underworld3/utilities/_cylindrical_wall.py index 305899c9b..734ca8a37 100644 --- a/src/underworld3/utilities/_cylindrical_wall.py +++ b/src/underworld3/utilities/_cylindrical_wall.py @@ -13,6 +13,9 @@ This lagged strong free-slip/traction projection is a UW3 adaptation, not their Nitsche formulation. Their Section 4.2 assumes body-fitted Neumann boundaries; its stability/convergence analysis does not cover this implementation. + +Only nondimensional mesh coordinates, fields and input scalars are supported. +Active reference-unit scaling is outside this helper's supported contract. """ import numpy as np diff --git a/tests/parallel/test_1084_rotated_datum_collectives.py b/tests/parallel/test_1084_rotated_datum_collectives.py index a0191790f..49aaf9a2e 100644 --- a/tests/parallel/test_1084_rotated_datum_collectives.py +++ b/tests/parallel/test_1084_rotated_datum_collectives.py @@ -5,12 +5,16 @@ """ import numpy as np +import pytest import underworld3 as uw from mpi4py import MPI from underworld3.utilities.rotated_bc import _datum_batches, _eval_datum +pytestmark = pytest.mark.level_1 + +@pytest.mark.tier_b def test_registered_boundaries_determine_collective_schedule_not_local_nodes(): rank = uw.mpi.rank boundaries = ["Upper", "Lower"] @@ -27,6 +31,7 @@ def test_registered_boundaries_determine_collective_schedule_not_local_nodes(): assert coords.shape == (0, 2) +@pytest.mark.tier_b def test_fe_datum_evaluation_accepts_empty_coordinate_batches_on_other_ranks(): mesh = uw.meshing.Annulus(radiusInner=1.22, radiusOuter=2.22, cellSize=1/32, qdegree=2) field = uw.discretisation.MeshVariable("DatumCollectiveP2", mesh, 1, degree=2, continuous=True) From 6649a60942bc3964227cef5fa26df4e6b2d18e28 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 25 Oct 2026 14:14:55 +1100 Subject: [PATCH 3/6] Add conservative spatial annulus trace reconstruction and gauge consistency Recover a separate circular P2 trace using endpoint-weighted local polynomial convolution, exact breakpoint-resolved circular projection, and constraints on the physical mean and Cartesian radial resultant. Keep the original CBF trace and cache unchanged; replicate only boundary data. Explicit support and degree control accuracy; do not claim spline/DG SIAC superconvergence. Align transfer pressure with the cached reaction gauge after constant pressure calibration. Distributed pressure-subvector differences must be uniform; nonuniform stale pressure changes are rejected. Constant pressure regression proves the original defect before the fix; its actual benchmark effect is small compared with the combined FE pressure/velocity moment error. Independent pure-helper/adapter contracts pass serially and on two ranks: polynomial moments, closed P6 kernel, direct nested integration, conservation, rotation/scale/permutation, graded meshes, empty boundary ranks and duplicated ghost values. Native wrapper and fresh-install gates are pending the batched container build; no model results from this commit are yet accepted. Next: build immutable versioned Mac/Gadi installs, run full serial/MPI gates and matched annulus refinement before production acceptance. Underworld development team with AI support from Claude Code (cherry picked from commit a882c34030de443923406f4fa9516f96be9785ac) --- .../cylindrical-geometry-correction.md | 73 +++++- .../guides/boundary-condition-rulings.md | 10 + .../cython/petsc_generic_snes_solvers.pyx | 42 +++ .../utilities/_cylindrical_reactions.py | 37 +++ .../utilities/_cylindrical_reconstruction.py | 194 ++++++++++++++ tests/test_1089_cylindrical_pressure_gauge.py | 104 ++++++++ tests/test_1089_cylindrical_reconstruction.py | 240 ++++++++++++++++++ 7 files changed, 699 insertions(+), 1 deletion(-) create mode 100644 src/underworld3/utilities/_cylindrical_reconstruction.py create mode 100644 tests/test_1089_cylindrical_pressure_gauge.py create mode 100644 tests/test_1089_cylindrical_reconstruction.py diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md index 5929484bb..7a8cd4308 100644 --- a/docs/advanced/cylindrical-geometry-correction.md +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -108,6 +108,69 @@ coefficient = stokes.cylindrical_boundary_coefficient( Coefficient extraction removes the circle-weighted mean and uses the exact radial-projection line Jacobian. It is not part of spatial reaction recovery. +## Spatial trace reconstruction + +Retain the circular CBF trace and reconstruct into a **distinct** continuous +scalar P2 field created alongside the other mesh variables before solving: + +```python +checks = stokes.reconstruct_cylindrical_trace( + wall=walls["Upper"], source=cbf_topography, target=reconstructed_topography, + degree=6, halfwidth_edges=5, +) +``` + +This is an explicitly selected spatial postprocessor. It does not change the +Stokes solution or the cached assembled reaction. The compact local-L2 kernel +uses the endpoint weight $w(z)=(1-z^2)^2$ and reproduces evaluation at zero for +polynomials through the requested even degree (4, 6 or 8). Degree 6 has the +normalized kernel + +$$ +K_6(z)=\frac{(1-z^2)^2}{8192} +(24255-315315z^2+945945z^4-765765z^6),\qquad |z|\leq1, +$$ + +and is zero outside its support. The physical support halfwidth is +`halfwidth_edges` times the mean circular edge arc length, not a prescribed +source wavelength. The currently measured $n=2$ annulus choice is degree 6 and +halfwidth 5. The degree describes the kernel moments; the output remains P2. + +Convolution integrates the existing P2 trace using its radial projection onto +the true circle. Integration is split at original and shifted source vertices, +including on nonuniform walls, before projection onto the circular consistent +P2 mass. A minimum circular-L2 correction preserves +$\int f\,ds$, $\int f\cos\theta\,ds$ and $\int f\sin\theta\,ds$. +The latter two are Cartesian components of the radial resultant, defined by +geometry rather than a fitted benchmark wavenumber. + +The reconstructed field **changes fine-scale nodal CBF moments**. Keeping every +such moment in the same P2 space would uniquely retain the original CBF trace. +Keep both fields, record which is supplied to later calculations, and assess +complete spatial errors as well as coefficients. Polynomial reproduction of the +kernel is not a proof of superconvergence of the complete annular Stokes chain. +The spline/DG SIAC results discussed by Docampo-Sanchez et al. (2020) do not +establish this weighted polynomial operator's accuracy on circular or graded +CBF traces. + +The finite support imposes a bandwidth limit: genuine rapid angular variations +within that support can be attenuated or amplified. For example, the degree-6 kernel with +halfwidth 5 has continuous response about 0.8213 for an eighth harmonic on a +32-edge circle. No automatic tuning to a source wavenumber is performed. The +caller must assess supported physical scales as well as unwanted spatial errors. +Returned checks record actual fine-scale nodal-moment changes, circular-L2 change +and physical-moment conservation, alongside the kernel degree and support. + +### Pressure calibration and the reaction cache + +Pressure calibration may change the solved pressure by a constant while the +cached Cartesian reaction retains the solve's earlier gauge. Circular recovery +checks that the difference is constant and evaluates its numerical Taylor trace +in the cached-reaction gauge. It leaves the calibrated pressure unchanged. +A nonconstant pressure change makes the cache stale and is rejected; solve +Stokes again before recovery. No analytical pressure or wall +reference is used for this synchronization. + ## Parallel scope and reproducibility Call the methods collectively. Only owned boundary-cell traces and boundary @@ -117,7 +180,7 @@ operator and volume solves remain distributed PETSc. Wall handles cache geometry and are valid only for their original mesh, fields and fixed parameters. Record pressure representation, quadrature, source amplitude, wall-update history, -reaction checks and the target actually supplied to gravity. Keep complete spatial +reaction checks, reconstruction degree/support and the target actually supplied to gravity. Keep complete spatial errors distinct from modal errors. A correction changes the numerical problem on an unchanged affine mesh; both baseline and corrected cases must be solved independently for a scientific before/after comparison. @@ -146,3 +209,11 @@ independently for a scientific before/after comparison. This is the CBF stress-recovery reference. UW3's circular moment transfer, Taylor reaction correction and weighted coefficient integration extend the calculation and are not results established by that paper. + +- Docampo-Sanchez, J., Jacobs, G. B., Li, X., and Ryan, J. K. (2020), + *Enhancing accuracy with a convolution filter: What works and why!*, + Computers & Fluids 213, 104727, + [doi:10.1016/j.compfluid.2020.104727](https://doi.org/10.1016/j.compfluid.2020.104727). + This supplies the accuracy-enhancing convolution context. UW3's endpoint-weighted + polynomial kernel, circular P2 projection and physical-moment restoration are + distinct from its spline/DG operator; its convergence results are not inherited. diff --git a/docs/developer/guides/boundary-condition-rulings.md b/docs/developer/guides/boundary-condition-rulings.md index df385384b..930734892 100644 --- a/docs/developer/guides/boundary-condition-rulings.md +++ b/docs/developer/guides/boundary-condition-rulings.md @@ -105,3 +105,13 @@ force and datum registration visible to the caller. P2 target without changing the raw reaction cache or fitting harmonics. See the [circular geometry guide](../../advanced/cylindrical-geometry-correction.md) for equations, limitations and source/pressure conventions. + +For an explicitly selected reconstructed spatial trace, use +`stokes.reconstruct_cylindrical_trace(wall=..., source=..., target=..., degree=6, +halfwidth_edges=5)` with distinct P2 fields. This leaves Stokes and the original +CBF trace intact, preserves the circular mean and Cartesian radial resultants, +and changes fine-scale nodal moments. It uses no prescribed source harmonic. +Circular recovery also synchronizes constant post-solve pressure calibration +with the cached reaction gauge; nonconstant changes invalidate that cache. +The [circular geometry guide](../../advanced/cylindrical-geometry-correction.md) +states the reconstruction's mathematical contract and limitations. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index fbf1a352e..9f4afa40b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8835,6 +8835,48 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): from underworld3.utilities._cylindrical_reactions import _recover_circular_topography return _recover_circular_topography(self, wall, field, buoyancy_scale) + def reconstruct_cylindrical_trace(self, *, wall, source, target, degree, halfwidth_edges): + r"""Reconstruct a circular spatial trace into a distinct continuous P2 field. + + Parameters + ---------- + wall : wall + Handle returned by :meth:`add_cylindrical_wall_transfer`. + source, target : scalar continuous P2 MeshVariable + Existing circular trace and distinct reconstructed target on this mesh. + Only this wall's target values are written; source and CBF cache remain intact. + degree : int + Polynomial reproduction degree of the compact weighted local-L2 + kernel, one of 4, 6 or 8. This does not change the P2 trace space. + halfwidth_edges : float + Support halfwidth in mean boundary arc lengths. The measured n=2 + annulus choice is degree 6 and halfwidth 5; no source wavenumber is used. + + Returns + ------- + dict + Support, actual nodal-moment changes and physical-moment checks. + + Notes + ----- + Call collectively. This postprocessing changes fine-scale nodal CBF + moments while preserving the circular mean and Cartesian radial-force + moments. It is not an unchanged CBF trace or a Stokes solve. Compact + support can attenuate or amplify genuine rapid variations; retain + the source and assess the physical bandwidth. Inputs must be + nondimensional; only boundary data are replicated. + + Accuracy-enhancing convolution is discussed by Docampo-Sanchez et al. + (2020), Computers & Fluids 213, 104727, + https://doi.org/10.1016/j.compfluid.2020.104727. This weighted polynomial, + circular projection and physical-moment correction are UW3 adaptations; + the spline/DG SIAC convergence theorems do not establish their accuracy. + """ + from underworld3.utilities._cylindrical_reconstruction import _reconstruct_cylindrical_trace + return _reconstruct_cylindrical_trace( + self, wall=wall, source=source, target=target, + degree=degree, halfwidth_edges=halfwidth_edges) + def cylindrical_boundary_coefficient(self, *, wall, field, wavenumber): r"""Integrate a spatial P2 wall field's cosine coefficient on the circle. diff --git a/src/underworld3/utilities/_cylindrical_reactions.py b/src/underworld3/utilities/_cylindrical_reactions.py index 8b1bb54e9..e18f3b8eb 100644 --- a/src/underworld3/utilities/_cylindrical_reactions.py +++ b/src/underworld3/utilities/_cylindrical_reactions.py @@ -25,6 +25,7 @@ from scipy.sparse import coo_matrix from scipy.sparse.linalg import spsolve from scipy.spatial import cKDTree +from petsc4py import PETSc from ._cylindrical_wall import _trace_rule, _numerical_trace @@ -97,6 +98,39 @@ def _validate_spatial_target(wall, field): raise ValueError("Circular wall fields require a continuous scalar P2 target on the same mesh.") +def _reaction_pressure_offset(stokes): + """Express current pressure in the gauge of the cached momentum reaction. + + Subtracting a constant after a solve is a valid pressure calibration. The + cached reaction and the polygon/circle stress moments must nevertheless use + the same constant: nodal radial normals vary over each affine edge, so an + inconsistent constant becomes a nonconstant scalar reaction after mass + inversion. Distributed PETSc vectors detect the single global offset; no + pressure field is gathered or modified. Nonconstant changes require a new + solve and cannot be interpreted as a gauge. + """ + pressure_is, pressure_dm = stokes._subdict["pressure"] + state = stokes._rotated_freeslip_info["U"] + cached = state.getSubVector(pressure_is) + current, difference = cached.duplicate(), cached.duplicate() + try: + current.set(0.) + pressure_dm.localToGlobal(stokes.p.vec, current) + cached.copy(difference) + difference.axpy(-1., current) + offset = float(difference.sum()/difference.getSize()) + difference.shift(-offset) + infinity = PETSc.NormType.NORM_INFINITY + tolerance = 64*np.finfo(float).eps*(1+cached.norm(infinity)+current.norm(infinity)) + if difference.norm(infinity) > tolerance: + raise ValueError("Pressure changed nonuniformly after the cached Stokes reaction; solve again before recovery.") + finally: + current.destroy() + difference.destroy() + state.restoreSubVector(pressure_is, cached) + return offset + + def _recover_circular_topography(stokes, wall, field, buoyancy_scale): r"""Write circular-wall topography from CBF reactions and numerical moments. @@ -117,10 +151,13 @@ def _recover_circular_topography(stokes, wall, field, buoyancy_scale): if distance.max() > 1e-10: raise ValueError(f"CBF nodes do not cover circular wall {wall.name}.") velocity, pressure = wall.cell_values() + pressure_offset = _reaction_pressure_offset(stokes) + pressure = pressure+pressure_offset trace = _numerical_trace(wall.vertices, velocity, pressure, wall.rule["points"], *wall.physics) nodes, reaction, checks = _circular_reaction_field( wall.edges, trace, gathered_values[index], wall.radius, wall.sign, wall.quadrature ) + checks["reaction_pressure_gauge_offset"] = pressure_offset np.testing.assert_array_equal(nodes, wall.nodes) distance, index = cKDTree(nodes).query(np.asarray(field.coords)) boundary = distance < 1e-10 diff --git a/src/underworld3/utilities/_cylindrical_reconstruction.py b/src/underworld3/utilities/_cylindrical_reconstruction.py new file mode 100644 index 000000000..3f4fe12cf --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_reconstruction.py @@ -0,0 +1,194 @@ +"""Private spatial reconstruction of circular P2 boundary traces. + +The compact local-L2 evaluation kernel reproduces polynomials with weight +``(1-z**2)**2`` and vanishes smoothly at its support endpoints. Convolve the +existing spatial trace, project onto the same circular P2 space, and preserve +its mean and Cartesian radial-resultant moments. No prescribed harmonic or +analytical solution enters the reconstruction. + +Accuracy-enhancing convolution is discussed by Docampo-Sanchez, J., Jacobs, +G. B., Li, X., and Ryan, J. K. (2020), "Enhancing accuracy with a convolution +filter: What works and why!", Computers & Fluids 213, 104727, +https://doi.org/10.1016/j.compfluid.2020.104727. The weighted polynomial kernel, +circular projection and physical-moment constraints here are UW3 adaptations. +This is not that paper's spline/DG operator and inherits no SIAC convergence +theorem for Stokes reactions or nonuniform circular meshes. + +Only boundary nodes are replicated. Sparse one-dimensional SciPy mass solves +are used, as in the existing circular reaction transfer; no volume field or +Stokes operator is gathered. Nondimensional inputs are required. +""" + +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.special import eval_jacobi +from scipy.spatial import cKDTree +import underworld3 as uw + + +def _reconstruction_kernel(points, degree): + r"""Evaluate the normalized weighted polynomial-reproducing kernel. + + The Jacobi polynomials :math:`P_j^{(2,2)}` are orthogonal with weight + :math:`w=(1-z^2)^2`. Thus :math:`w(z)\sum_j P_j(0)P_j(z)/\|P_j\|_w^2` + reproduces evaluation at zero through ``degree``. Its support is [-1, 1]. + """ + if degree not in (4, 6, 8): + raise ValueError("Circular trace reconstruction supports degrees 4, 6 and 8.") + points = np.asarray(points, dtype=float) + inside = np.abs(points) <= 1 + result = np.zeros_like(points) + z = points[inside] + polynomial = np.zeros_like(z) + for order in range(0, degree+1, 2): + norm = 32*(order+1)*(order+2)/((2*order+5)*(order+3)*(order+4)) + polynomial += eval_jacobi(order, 2, 2, 0)*eval_jacobi(order, 2, 2, z)/norm + result[inside] = (1-z*z)**2*polynomial + return result + + +def _recover_circular_trace(edges, nodal_values, radius, degree, halfwidth_edges, quadrature): + """Reconstruct a complete circular P2 trace from boundary arrays. + + Input values follow lexicographically unique endpoints and midpoints of + ``edges``, matching cylindrical reaction recovery. The support halfwidth + is ``halfwidth_edges`` times the mean physical boundary arc length. Split + integration at original and shifted source vertices, also on graded meshes. + This relaxes fine-scale CBF moments; the original trace is not overwritten. + """ + edges = np.asarray(edges, dtype=float) + if edges.ndim != 3 or edges.shape[1:] != (2, 2) or not np.all(np.isfinite(edges)): + raise ValueError("Circular trace edges require finite (edges, 2, 2) coordinates.") + if radius <= 0 or not np.isfinite(radius): + raise ValueError("Circular trace radius must be positive and finite.") + if degree not in (4, 6, 8): + raise ValueError("Circular trace reconstruction supports degrees 4, 6 and 8.") + if not np.isfinite(halfwidth_edges) or halfwidth_edges <= 0 or quadrature < 2: + raise ValueError("Require positive finite halfwidth and quadrature >= 2.") + vertices = np.unique(edges.reshape(-1, 2), axis=0) + if not np.allclose(np.linalg.norm(vertices, axis=1), radius, rtol=1e-10, atol=0): + raise ValueError("Boundary vertices must lie on the specified circle.") + angles = np.mod(np.arctan2(vertices[:, 1], vertices[:, 0]), 2*np.pi) + vertices = vertices[np.argsort(angles)] + ordered = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + original_pairs = {tuple(sorted(map(tuple, pair))) for pair in edges} + ordered_pairs = {tuple(sorted(map(tuple, pair))) for pair in ordered} + if len(edges) != len(vertices) or original_pairs != ordered_pairs: + raise ValueError("Input edges must form one complete circular wall.") + edge_nodes = np.concatenate((ordered, ordered.mean(axis=1)[:, None]), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + values = np.asarray(nodal_values, dtype=float).reshape(-1) + if len(values) != len(nodes) or not np.all(np.isfinite(values)): + raise ValueError("Require one finite value per unique P2 boundary node.") + starts = np.unwrap(np.mod(np.arctan2(ordered[:, 0, 1], ordered[:, 0, 0]), 2*np.pi)) + origin, period, count = starts[0], 2*np.pi, len(ordered) + ends = np.r_[starts[1:], origin+period] + if np.any(ends-starts >= np.pi): + raise ValueError("Circular trace arcs must each span less than half a circle.") + width = halfwidth_edges*period/count + if not 0 < width < np.pi: + raise ValueError("The recovery halfwidth must span less than half a circle.") + gauss, weight = np.polynomial.legendre.leggauss(quadrature) + + def basis(theta): + query = (theta-origin) % period+origin + selected = np.searchsorted(starts, query, side="right")-1 + unit = np.stack((np.cos(query), np.sin(query)), axis=-1) + first = ordered[selected, 0] + direction = ordered[selected, 1]-first + t = -(first[..., 0]*unit[..., 1]-first[..., 1]*unit[..., 0])/( + direction[..., 0]*unit[..., 1]-direction[..., 1]*unit[..., 0]) + shape = np.stack(((1-t)*(1-2*t), t*(2*t-1), 4*t*(1-t)), axis=-1) + return selected, shape + + def unwrapped_cell(theta): + return (np.floor((theta-origin)/period).astype(int)*count + + np.searchsorted(starts, (theta-origin) % period+origin, side="right")-1) + + def convolution(theta): + first, last = unwrapped_cell(theta-width), unwrapped_cell(theta+width) + result = np.zeros_like(theta) + for offset in range(int(np.max(last-first))+1): + cell = first+offset + wrap, turn = cell % count, cell//count + lo = np.maximum(starts[wrap]+turn*period, theta-width) + hi = np.minimum(ends[wrap]+turn*period, theta+width) + keep = (cell <= last) & (hi > lo) + points = ((lo[keep, None]+hi[keep, None])/2 + + (hi[keep, None]-lo[keep, None])*gauss/2) + selected, shape = basis(points) + samples = np.sum(shape*values[index[selected]], axis=-1) + kernel = _reconstruction_kernel((points-theta[keep, None])/width, degree) + result[keep] += np.sum( + (hi[keep, None]-lo[keep, None])*weight/2*kernel*samples/width, axis=1) + return result + + # Exact circular arclength defines the projection and physical moments. + points = ((starts[:, None]+ends[:, None])/2 + + (ends[:, None]-starts[:, None])*gauss/2).ravel() + weights = (radius*(ends[:, None]-starts[:, None])*weight/2).ravel() + selected, shape = basis(points) + local = np.einsum("q,qi,qj->qij", weights, shape, shape) + rows = np.broadcast_to(index[selected, :, None], local.shape).ravel() + columns = np.broadcast_to(index[selected, None, :], local.shape).ravel() + mass = coo_matrix((local.ravel(), (rows, columns)), shape=(len(nodes), len(nodes))).tocsr() + constraints = np.zeros((3, len(nodes))) + for row, radial_component in enumerate((np.ones_like(points), np.cos(points), np.sin(points))): + for j in range(3): + np.add.at(constraints[row], index[selected, j], weights*shape[:, j]*radial_component) + + cuts = (np.r_[starts, starts-width, starts+width]-origin) % period+origin + cuts = np.sort(np.r_[cuts, origin, origin+period]) + cuts = cuts[np.r_[True, np.diff(cuts) > 1e-13]] + cuts[-1] = origin+period + points = ((cuts[:-1, None]+cuts[1:, None])/2 + + (cuts[1:, None]-cuts[:-1, None])*gauss/2).ravel() + weights = (radius*(cuts[1:, None]-cuts[:-1, None])*weight/2).ravel() + selected, shape = basis(points) + samples = convolution(points) + load = np.zeros(len(nodes)) + for j in range(3): + np.add.at(load, index[selected, j], weights*shape[:, j]*samples) + recovered = spsolve(mass, load) + correction_basis = spsolve(mass, constraints.T) + recovered += correction_basis@np.linalg.solve( + constraints@correction_basis, constraints@(values-recovered)) + change = recovered-values + moment_change = constraints@change + checks = dict(degree=degree, halfwidth_edges=halfwidth_edges, + physical_arc_halfwidth=radius*width, quadrature=quadrature, + maximum_conserved_moment_error=float(np.max(np.abs(moment_change))), + maximum_nodal_moment_change=float(np.max(np.abs(mass@change))), + circular_l2_change=float(np.sqrt(change@(mass@change))), + circular_mean_change=float(moment_change[0]/(2*np.pi*radius)), + cartesian_radial_resultant_change=moment_change[1:].tolist(), + fine_scale_cbf_moments_preserved=False) + return nodes, recovered, checks + + +def _reconstruct_cylindrical_trace(stokes, *, wall, source, target, degree, halfwidth_edges): + """Write a separate reconstructed trace; replicate boundary nodes only.""" + from ._cylindrical_reactions import _validate_spatial_target + + if wall.mesh is not stokes.mesh: + raise ValueError("The wall handle must belong to this Stokes mesh.") + if source is target: + raise ValueError("Source and reconstructed target must be distinct fields.") + _validate_spatial_target(wall, source) + _validate_spatial_target(wall, target) + boundary = wall.local_boundary + coordinates = np.concatenate(uw.mpi.comm.allgather(np.asarray(source.coords)[boundary])) + values = np.concatenate(uw.mpi.comm.allgather(np.asarray(source.array)[boundary, 0, 0])) + unique, first, inverse = np.unique(coordinates, axis=0, return_index=True, return_inverse=True) + if not np.allclose(values, values[first][inverse], rtol=1e-12, atol=0): + raise ValueError("Shared boundary nodes have inconsistent source values.") + distance, index = cKDTree(unique).query(wall.nodes) + if np.max(distance) > 1e-10: + raise ValueError("The source field does not cover the circular wall nodes.") + nodes, reconstructed, checks = _recover_circular_trace( + wall.edges, values[first][index], wall.radius, degree, halfwidth_edges, 12) + np.testing.assert_array_equal(nodes, wall.nodes) + target.array[boundary, 0, 0] = reconstructed[wall.local_index[boundary]] + return checks diff --git a/tests/test_1089_cylindrical_pressure_gauge.py b/tests/test_1089_cylindrical_pressure_gauge.py new file mode 100644 index 000000000..1e9c7aef2 --- /dev/null +++ b/tests/test_1089_cylindrical_pressure_gauge.py @@ -0,0 +1,104 @@ +"""Constant-pressure invariance of circular CBF transfer after calibration. + +No Stokes solve is required: independently integrated constant-pressure loads +are the cached reaction, while the current pressure field is shifted by a gauge. +The physical mean-free circular stress must remain identically zero. +""" + +from types import SimpleNamespace + +import numpy as np +import pytest +from petsc4py import PETSc +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve + +from underworld3.utilities._cylindrical_reactions import _recover_circular_topography + + +pytestmark = pytest.mark.level_1 + + +class _PressureDM: + """Preserve the native distributed pressure-vector contract in this fixture.""" + + def localToGlobal(self, source, target): + source.copy(target) + + +def _constant_pressure_fixture(shift: float): + """Build exact P2 line loads on a polygon and a separately shifted pressure.""" + angle = np.linspace(0, 2*np.pi, 33)[:-1] + first = np.column_stack((np.cos(angle), np.sin(angle))) + edges = np.stack((first, np.roll(first, -1, axis=0)), axis=1) + vertices = np.concatenate((edges, .9*edges.mean(axis=1)[:, None]), axis=1) + edge_nodes = np.concatenate((edges, edges.mean(axis=1)[:, None]), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + index = index.reshape(-1, 3) + t, weights = np.polynomial.legendre.leggauss(12) + t, weights = (t+1)/2, weights/2 + shape = np.array([(1-t)*(1-2*t), t*(2*t-1), 4*t*(1-t)]).T + delta = edges[:, 1]-edges[:, 0] + length = np.linalg.norm(delta, axis=1) + normal = np.column_stack((delta[:, 1], -delta[:, 0]))/length[:, None] + nodal_normal = edge_nodes/np.linalg.norm(edge_nodes, axis=-1)[..., None] + local_mass = np.einsum("e,q,qi,qj->eij", length, weights, shape, shape) + rows = np.broadcast_to(index[:, :, None], local_mass.shape).ravel() + columns = np.broadcast_to(index[:, None, :], local_mass.shape).ravel() + mass = coo_matrix((local_mass.ravel(), (rows, columns)), shape=(len(nodes), len(nodes))).tocsr() + # sigma=-I gives the opposite solver-reaction convention +n_i.n_facet. + loads = np.zeros(len(nodes)) + local_load = np.einsum("e,q,qi,eid,ed->ei", length, weights, shape, nodal_normal, normal) + np.add.at(loads, index.ravel(), local_load.ravel()) + reaction = spsolve(mass, loads) + uniform = mass@np.ones(len(nodes)) + reaction -= uniform@reaction/uniform.sum() + cached = PETSc.Vec().createMPI((len(vertices)*3, None), comm=PETSc.COMM_WORLD) + cached.set(1) + current = cached.duplicate() + current.set(1+shift) + start, end = cached.getOwnershipRange() + pressure_is = PETSc.IS().createStride(end-start, first=start, step=1, comm=PETSc.COMM_WORLD) + mesh = object() + wall = SimpleNamespace(mesh=mesh, name="Upper", nodes=nodes, vertices=vertices, + edges=edges, radius=1., sign=1, quadrature=12, physics=(1., .5, 1., 2, 0.), + cell_values=lambda: (np.zeros((len(vertices), 6, 2)), np.full((len(vertices), 3), 1+shift))) + wall.rule = dict(points=edges[:, :1]+t[None, :, None]*delta[:, None]) + field = SimpleNamespace(mesh=mesh, degree=2, continuous=True, num_components=1, + coords=nodes, array=np.zeros((len(nodes), 1, 1))) + solver = SimpleNamespace(p=SimpleNamespace(vec=current), + _subdict={"pressure": (pressure_is, _PressureDM())}, + _rotated_freeslip_info={"U": cached}, + boundary_normal_traction=lambda boundary, mass: (nodes, reaction)) + return solver, wall, field, (cached, current, pressure_is) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("shift", [0., -1., .37]) +def test_circular_recovery_keeps_constant_pressure_gauge_out_of_shape(shift): + solver, wall, field, owned = _constant_pressure_fixture(shift) + cached_values, current_values = owned[0].array.copy(), owned[1].array.copy() + try: + for repeat in range(2): + checks = _recover_circular_topography(solver, wall, field, buoyancy_scale=1.) + np.testing.assert_allclose(field.array, 0, rtol=0, atol=4e-14, + err_msg="A constant pressure gauge must not become a circular facet-scale stress.") + np.testing.assert_allclose(checks["reaction_pressure_gauge_offset"], -shift, rtol=0, atol=2e-15) + np.testing.assert_array_equal(owned[0].array, cached_values) + np.testing.assert_array_equal(owned[1].array, current_values) + finally: + for item in owned: + item.destroy() + + +@pytest.mark.tier_b +def test_nonconstant_pressure_change_is_not_accepted_as_a_gauge(): + solver, wall, field, owned = _constant_pressure_fixture(0.) + try: + # Change the first owned pressure entry on each rank, rather than a constant. + solver.p.vec.array[0] += .01 + with pytest.raises(ValueError, match="Pressure changed nonuniformly"): + _recover_circular_topography(solver, wall, field, buoyancy_scale=1.) + finally: + for item in owned: + item.destroy() diff --git a/tests/test_1089_cylindrical_reconstruction.py b/tests/test_1089_cylindrical_reconstruction.py new file mode 100644 index 000000000..c85576ed8 --- /dev/null +++ b/tests/test_1089_cylindrical_reconstruction.py @@ -0,0 +1,240 @@ +"""Mathematical contracts for spatial circular trace reconstruction. + +No Stokes solve is required. Independent kernel moments, circular integration, +frame transformations and conservation define the baselines. +""" + +import numpy as np +import pytest +import underworld3 as uw + +pytestmark = pytest.mark.level_1 + + +@pytest.mark.tier_b +def test_stokes_exposes_explicit_spatial_trace_reconstruction(): + """The solver capability must be reachable without a private utility import.""" + assert callable(getattr(uw.systems.Stokes, "reconstruct_cylindrical_trace", None)), ( + "Stokes lacks reconstruct_cylindrical_trace(*, wall, source, target, degree, halfwidth_edges)." + ) + + +def _ring(count, graded, radius): + steps = 1 + (0.7*np.sin(0.7*np.arange(count)) if graded else np.zeros(count)) + angles = 2*np.pi*np.r_[0, np.cumsum(steps[:-1])]/steps.sum() + vertices = radius*np.stack((np.cos(angles), np.sin(angles)), axis=-1) + edges = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + nodes = np.unique(np.concatenate((edges, edges.mean(axis=1)[:, None]), axis=1).reshape(-1, 2), axis=0) + return edges, nodes + + +def _reference_reconstruction(edges, values, radius, degree, halfwidth_edges): + """Direct scalar quadrature with an independent angular/chord coordinate map.""" + from underworld3.utilities._cylindrical_reconstruction import _reconstruction_kernel + + vertices = np.unique(edges.reshape(-1, 2), axis=0) + starts = np.sort(np.mod(np.arctan2(vertices[:, 1], vertices[:, 0]), 2*np.pi)) + ends = np.r_[starts[1:], starts[0]+2*np.pi] + ordered_vertices = radius*np.stack((np.cos(starts), np.sin(starts)), axis=-1) + ordered_edges = np.stack((ordered_vertices, np.roll(ordered_vertices, -1, axis=0)), axis=1) + nodes = np.unique(np.concatenate((edges, edges.mean(axis=1)[:, None]), axis=1).reshape(-1, 2), axis=0) + expected_nodes = np.concatenate((ordered_edges, ordered_edges.mean(axis=1)[:, None]), axis=1) + index = np.argmin(np.linalg.norm(expected_nodes[:, :, None]-nodes[None, None], axis=-1), axis=-1) + width = halfwidth_edges*2*np.pi/len(edges) + gauss, weight = np.polynomial.legendre.leggauss(24) + + def shapes(theta): + wrapped = (theta-starts[0]) % (2*np.pi)+starts[0] + cell = np.searchsorted(starts, wrapped, side="right")-1 + half_arc = (ends[cell]-starts[cell])/2 + midpoint = (ends[cell]+starts[cell])/2 + chord = (np.tan(wrapped-midpoint)/np.tan(half_arc)+1)/2 + return cell, np.stack(((1-chord)*(1-2*chord), chord*(2*chord-1), 4*chord*(1-chord)), axis=-1) + + cuts = np.sort(np.r_[0, 2*np.pi, starts, (starts-width) % (2*np.pi), (starts+width) % (2*np.pi)]) + cuts = cuts[np.r_[True, np.diff(cuts)>1e-13]] + mass = np.zeros((len(nodes), len(nodes))) + load = np.zeros(len(nodes)) + conserved = np.zeros((3, len(nodes))) + for lo, hi in zip(cuts[:-1], cuts[1:]): + for theta, outer_weight in zip((lo+hi)/2+(hi-lo)*gauss/2, radius*(hi-lo)*weight/2): + cell, basis = shapes(np.asarray(theta)) + source_cuts = np.sort(np.r_[theta-width, theta+width, + np.tile(starts, 3)+np.repeat([-2*np.pi, 0, 2*np.pi], len(starts))]) + source_cuts = source_cuts[(source_cuts>=theta-width) & (source_cuts<=theta+width)] + integral = 0.0 + for first, last in zip(source_cuts[:-1], source_cuts[1:]): + points = (first+last)/2+(last-first)*gauss/2 + source_cell, source_basis = shapes(points) + source_value = np.sum(source_basis*values[index[source_cell]], axis=-1) + integral += np.sum((last-first)*weight/2*source_value + *_reconstruction_kernel((points-theta)/width, degree)/width) + rows = index[cell] + mass[np.ix_(rows, rows)] += outer_weight*np.outer(basis, basis) + load[rows] += outer_weight*basis*integral + conserved[:, rows] += outer_weight*np.outer([1, np.cos(theta), np.sin(theta)], basis) + result = np.linalg.solve(mass, load) + correction = np.linalg.solve(mass, conserved.T) + result += correction@np.linalg.solve(conserved@correction, conserved@(values-result)) + return result, conserved + + +@pytest.mark.tier_b +@pytest.mark.parametrize("degree", [4, 6, 8]) +def test_kernel_reproduces_polynomial_moments(degree): + from underworld3.utilities._cylindrical_reconstruction import _reconstruction_kernel + + points, weights = np.polynomial.legendre.leggauss(20) + kernel = _reconstruction_kernel(points, degree) + for power in range(degree+1): + expected = 1 if power == 0 else 0 + np.testing.assert_allclose(np.sum(weights*kernel*points**power), expected, rtol=0, atol=2e-13) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("graded", [False, True]) +def test_constant_and_multimode_trace_against_direct_variational_integral(graded): + from underworld3.utilities._cylindrical_reconstruction import _recover_circular_trace + + edges, nodes = _ring(24, graded, 1.7) + _, constant, _ = _recover_circular_trace(edges, np.full(len(nodes), 1.234), 1.7, 6, 5, 12) + np.testing.assert_allclose(constant, 1.234, rtol=0, atol=3e-12) + theta = np.arctan2(nodes[:, 1], nodes[:, 0]) + values = 1+.2*np.cos(theta)-.3*np.sin(theta)+.4*np.cos(3*theta)-.1*np.sin(7*theta) + original = values.copy() + returned_nodes, actual, _ = _recover_circular_trace(edges, values, 1.7, 6, 5, 12) + reference, conserved = _reference_reconstruction(edges, values, 1.7, 6, 5) + np.testing.assert_array_equal(returned_nodes, nodes) + np.testing.assert_allclose(actual, reference, rtol=2e-10, atol=3e-11) + np.testing.assert_allclose(conserved@actual, conserved@values, rtol=2e-12, atol=2e-12) + np.testing.assert_array_equal(values, original) + + +@pytest.mark.tier_b +def test_edge_permutation_rotation_and_physical_scaling(): + from underworld3.utilities._cylindrical_reconstruction import _recover_circular_trace + + edges, nodes = _ring(40, True, 1.22) + theta = np.arctan2(nodes[:, 1], nodes[:, 0]) + values = 1+.3*np.cos(3*theta)-.1*np.sin(5*theta)+.2*np.cos(8*theta) + _, expected, _ = _recover_circular_trace(edges, values, 1.22, 6, 5, 12) + reordered = edges[np.random.default_rng(223).permutation(len(edges)), ::-1] + _, permuted, _ = _recover_circular_trace(reordered, values, 1.22, 6, 5, 12) + np.testing.assert_allclose(permuted, expected, rtol=0, atol=2e-12) + _, scaled, _ = _recover_circular_trace(edges*1.7, values, 1.22*1.7, 6, 5, 12) + np.testing.assert_allclose(scaled, expected, rtol=0, atol=2e-12) + angle = .324 + rotation = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) + rotated_edges = edges@rotation.T + rotated_nodes = np.unique(np.concatenate((rotated_edges, rotated_edges.mean(axis=1)[:, None]), axis=1).reshape(-1, 2), axis=0) + index = np.argmin(np.linalg.norm(rotated_nodes[:, None]-(nodes@rotation.T)[None], axis=-1), axis=1) + _, rotated, _ = _recover_circular_trace(rotated_edges, values[index], 1.22, 6, 5, 12) + np.testing.assert_allclose(rotated, expected[index], rtol=0, atol=2e-12) + + +@pytest.mark.tier_b +def test_boundary_only_field_adapter_preserves_source_and_other_nodes(): + from types import SimpleNamespace + from underworld3.utilities._cylindrical_reconstruction import _reconstruct_cylindrical_trace + + edges, nodes = _ring(24, False, 1.7) + # Empty-boundary ranks participate in the same collective adapter path. + selected = np.arange(len(nodes)) if uw.mpi.rank == 0 else np.empty(0, dtype=int) + coordinates = np.vstack((nodes[selected], [[0, 0]])) + mesh = SimpleNamespace(dim=2, cdim=2) + wall = SimpleNamespace(mesh=mesh, edges=edges, nodes=nodes, radius=1.7, quadrature=12, + local_boundary=np.r_[np.ones(len(selected), dtype=bool), False], + local_index=np.r_[selected, 0]) + source = SimpleNamespace(mesh=mesh, degree=2, continuous=True, num_components=1, + coords=coordinates, array=np.full((len(coordinates), 1, 1), 2.3)) + source.array[-1, 0, 0] = 777 + target = SimpleNamespace(mesh=mesh, degree=2, continuous=True, num_components=1, + coords=coordinates, array=np.full((len(coordinates), 1, 1), -999.0)) + original = source.array.copy() + _reconstruct_cylindrical_trace(SimpleNamespace(mesh=mesh), wall=wall, + source=source, target=target, degree=6, halfwidth_edges=5) + np.testing.assert_allclose(target.array[:-1, 0, 0], 2.3, rtol=0, atol=3e-12) + assert target.array[-1, 0, 0] == -999 + np.testing.assert_array_equal(source.array, original) + with pytest.raises(ValueError, match="distinct"): + _reconstruct_cylindrical_trace(SimpleNamespace(mesh=mesh), wall=wall, + source=source, target=source, degree=6, halfwidth_edges=5) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("invalid", ["missing_edge", "wrong_radius", "wrong_values", "degree", "width"]) +def test_incomplete_or_unsupported_circular_trace_fails_explicitly(invalid): + from underworld3.utilities._cylindrical_reconstruction import _recover_circular_trace + + edges, nodes = _ring(24, False, 1.7) + values, radius, degree, width = np.ones(len(nodes)), 1.7, 6, 5 + if invalid == "missing_edge": + edges = edges[:-1] + elif invalid == "wrong_radius": + radius = 1.9 + elif invalid == "wrong_values": + values = values[:-1] + elif invalid == "degree": + degree = 3 + else: + width = 12 + with pytest.raises(ValueError): + _recover_circular_trace(edges, values, radius, degree, width, 12) + + +@pytest.mark.tier_b +def test_degree_six_kernel_matches_exact_weighted_reproducing_polynomial(): + from underworld3.utilities._cylindrical_reconstruction import _reconstruction_kernel + + z = np.linspace(-1, 1, 31) + expected = (1-z*z)**2*(24255-315315*z*z+945945*z**4-765765*z**6)/8192 + np.testing.assert_allclose(_reconstruction_kernel(z, 6), expected, rtol=0, atol=2e-14) + np.testing.assert_array_equal(_reconstruction_kernel(np.array([-2, 2]), 6), [0, 0]) + + +@pytest.mark.tier_b +def test_degree_six_kernel_has_expected_finite_bandwidth_response(): + """At E=32, n=8, halfwidth=5 the analytic response is about 0.8213.""" + from underworld3.utilities._cylindrical_reconstruction import _reconstruction_kernel + + z, weights = np.polynomial.legendre.leggauss(64) + frequency = 5*np.pi/2 + actual = np.sum(weights*_reconstruction_kernel(z, 6)*np.cos(frequency*z)) + # Exact integral of the degree-ten polynomial times cos(a*z), sin(a)=1. + expected = (3465/frequency**3-2848230/frequency**5 + +265540275/frequency**7-5037157125/frequency**9 + +10854718875/frequency**11)/16 + np.testing.assert_allclose(actual, expected, rtol=0, atol=3e-13) + + +@pytest.mark.tier_b +def test_partitioned_nonconstant_trace_with_shared_boundary_ghosts(): + from types import SimpleNamespace + from underworld3.utilities._cylindrical_reconstruction import ( + _reconstruct_cylindrical_trace, _recover_circular_trace, + ) + + edges, nodes = _ring(24, True, 1.7) + theta = np.arctan2(nodes[:, 1], nodes[:, 0]) + values = 1+.2*np.cos(theta)-.3*np.sin(theta)+.4*np.cos(3*theta)-.1*np.sin(7*theta) + _, reference, _ = _recover_circular_trace(edges, values, 1.7, 6, 5, 12) + # Distributed owned entries plus the same five boundary ghosts on every rank. + owned = np.arange(len(nodes))[uw.mpi.rank::uw.mpi.size] + selected = np.unique(np.r_[owned, np.arange(5)]) + coordinates = np.vstack((nodes[selected], [[0, 0]])) + mesh = SimpleNamespace(dim=2, cdim=2) + wall = SimpleNamespace(mesh=mesh, edges=edges, nodes=nodes, radius=1.7, quadrature=12, + local_boundary=np.r_[np.ones(len(selected), dtype=bool), False], + local_index=np.r_[selected, 0]) + source = SimpleNamespace(mesh=mesh, degree=2, continuous=True, num_components=1, + coords=coordinates, array=np.r_[values[selected], 777].reshape(-1, 1, 1)) + target = SimpleNamespace(mesh=mesh, degree=2, continuous=True, num_components=1, + coords=coordinates, array=np.full((len(coordinates), 1, 1), -999.0)) + original = source.array.copy() + checks = _reconstruct_cylindrical_trace(SimpleNamespace(mesh=mesh), wall=wall, + source=source, target=target, degree=6, halfwidth_edges=5) + np.testing.assert_allclose(target.array[:-1, 0, 0], reference[selected], rtol=0, atol=2e-12) + assert target.array[-1, 0, 0] == -999 + np.testing.assert_array_equal(source.array, original) + assert checks["fine_scale_cbf_moments_preserved"] is False + assert checks["maximum_conserved_moment_error"] < 2e-12 From a00d9111243bb1badefa7017ff997fd61af291c0 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 25 Oct 2026 15:10:16 +1100 Subject: [PATCH 4/6] Add numerical full-vector transfer for zero-slip circular walls Extend one-sided affine P2 velocity polynomials to true circles and register distinct essential vector fields once. Enforce one joint polygon flux constraint with the minimum-L2 normal midpoint-bubble correction, retaining vertex and tangential data. Replicate only owned boundary-cell traces, including empty-rank/ghost-safe label extraction; report physical-circle residuals and actual standard-solver verdicts at every fixed-point step. Cite the exact Taylor-Hood trial/test precedent in Ruas 2019 alongside shifted-boundary foundations, while documenting symmetric-stress and flux-projection differences. Clarify that pressure-lift natural data are omitted where full velocity test traces vanish. Seventeen independent polynomial, flux, geometry and adapter tests pass serially and on each MPI2 rank as development helper checks. Fresh compiled/native gates and model acceptance remain pending; no all-case accuracy claim yet. Next: versioned Mac/Gadi builds and separate smooth-wall/delta-pressure pilots. Underworld development team with AI support from Claude Code --- .../cylindrical-geometry-correction.md | 53 +++ .../cython/petsc_generic_snes_solvers.pyx | 89 ++++- .../utilities/_cylindrical_dirichlet.py | 354 ++++++++++++++++++ tests/test_1089_cylindrical_dirichlet.py | 317 ++++++++++++++++ 4 files changed, 811 insertions(+), 2 deletions(-) create mode 100644 src/underworld3/utilities/_cylindrical_dirichlet.py create mode 100644 tests/test_1089_cylindrical_dirichlet.py diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md index 7a8cd4308..41de27f98 100644 --- a/docs/advanced/cylindrical-geometry-correction.md +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -31,6 +31,7 @@ lift = stokes.cylindrical_pressure_lifting( amplitude=1.0, quadrature_points=12, ) stokes.bodyforce = lift.regular_force.sym +# Register only on a natural/free-slip wall; omit for full velocity Dirichlet. stokes.add_natural_bc(lift.natural_datum, "Upper") physical_pressure = stokes.p.sym[0] + lift.pressure_jump ``` @@ -75,6 +76,58 @@ supported contract. The helper does not install the pressure projection or velocity nullspace itself. Configure the rigid-rotation velocity modes separately when required by the enclosing free-slip system. +## Full-vector zero-slip wall transfer + +The separate zero-slip API registers one vector P2 essential datum per wall. +For each boundary node $x_h$, project to $P(x_h)=R x_h/|x_h|$ and update + +$$g_h(x_h)=u_h(x_h)-u_h(P(x_h)).$$ + +This evaluates the complete one-sided P2 polynomial, including the quadratic +Taylor term. The radial shift has the same coordinate sign at both walls; +the inner/outer normal sign enters only polygon-flux integration. No exact +Stokes field, pressure trace or harmonic coefficient supplies the update. + +```python +walls = { + "Lower": stokes.add_cylindrical_zero_slip_transfer( + boundary="Lower", radius=1.22, sign=-1, quadrature=12), + "Upper": stokes.add_cylindrical_zero_slip_transfer( + boundary="Upper", radius=2.22, sign=1, quadrature=12), +} +stokes.petsc_use_pressure_nullspace = True +history = stokes.solve_cylindrical_zero_slip_updates(walls=walls, updates=4) +``` + +Use standard Stokes with both velocity components prescribed; there is no +rigid-rotation nullspace. Creating the handles registers the essential fields +once. Later updates change their values before each standard solve. Pressure +extraction is independent: zero-slip case3 uses the same regular force and +physical $q+p_J$, with no outer natural datum because its test trace is zero. + +Before each solve, remove the combined polygon flux using the normal P2 +midpoint bubbles on both walls. For total chord length $L$ and raw net flux $F$, +subtract $c n_h$ at every midpoint, where $c=F/(2L/3)$. Vertex and midpoint +tangential values remain unchanged. This is the minimum integrated-L2 change +within that bubble span. It enforces one global compatibility condition; +individual wall fluxes are recorded rather than fitted independently. + +`cylindrical_zero_slip_history` retains accepted standard-solver reports, +flux before/after projection, its L2 change, and true-circle nodal/full-Gauss +velocity residuals. At a fixed point a nonzero $c$ perturbs the true-circle +constraint; verify that defect and its refinement behaviour. Four updates are +an explicit experiment schedule, not a convergence guarantee. + +Ruas (2019), DOI10.1016/j.camwa.2018.05.029, sections2/4 and equations29–30, +provides the close P2 true-boundary trial/canonical polygon-test precedent. +UW3 uses symmetric stress and a global flux projection instead of that exact +formulation. Atallah et al. (2020), section4, and Burman et al. (2018), +https://arxiv.org/abs/1801.07463, motivate shifted Dirichlet data through different +weak formulations. None of these results is claimed as a proof of this strong +numerical fixed point, especially for the delta case's piecewise regularity. +Independent polynomial/flux/parallel tests pass; model acceptance remains pending +in the benchmark investigation. + ## Spatial circular reaction and topography The circular reaction moments are diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 9f4afa40b..c547f2710 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8673,8 +8673,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ``pressure_jump``, ``regular_force`` (DG P2 MeshVariable), ``natural_datum``, ``pressure_jump_integral`` and ``checks``. Reconstruct physical pressure as :math:`p=q+p_J`. Assign the - force and register ``natural_datum`` on the outer boundary; the - latter is :math:`-p_J\Gamma` and does not vanish on outer chords. + force and register ``natural_datum`` on the natural outer-wall + components. This term is :math:`-p_J\Gamma` and does not vanish on + outer chords. Fully prescribed velocity has zero test trace, so + no outer natural datum is registered for zero-slip walls. Notes ----- @@ -8798,6 +8800,89 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): from underworld3.utilities._cylindrical_wall import _solve_with_wall_updates return _solve_with_wall_updates(self, walls, updates) + def add_cylindrical_zero_slip_transfer(self, *, boundary, radius, sign, quadrature): + r"""Register numerical full-velocity transfer to a zero-slip circle. + + Parameters + ---------- + boundary : str + Complete circular wall label on the affine annulus mesh. + radius : float + Positive physical-circle radius through the wall vertices. + sign : int + Outward-domain normal sign: -1 at the inner wall, +1 at the outer. + It enters flux integration, not the radial geometric displacement. + quadrature : int + Gauss points per edge for numerical trace/flux audits. + + Returns + ------- + wall + Boundary-only handle with a distinct vector P2 auxiliary ``field``. + Its essential condition is registered once. Use both wall handles + with :meth:`solve_cylindrical_zero_slip_updates`. + + Notes + ----- + Supported scope is nondimensional affine 2D P2/continuous-P1 Stokes + with zero physical velocity at both circular walls. Numerical P2 + polynomial continuation supplies both velocity components; no reference + Stokes field or prescribed harmonic enters the update. Global polygon + flux compatibility changes only normal edge-midpoint data and keeps + physical-circle vertex data zero. No rotational nullspace is removed. + + References + ---------- + Atallah, Canuto and Scovazzi (2020), "The second-generation Shifted + Boundary Method and its numerical analysis", CMAME 372, 113341, + https://doi.org/10.1016/j.cma.2020.113341, Section 4; Burman, Hansbo and + Larson (2018), "A cut finite element method with boundary value + correction for the incompressible Stokes equations", + https://arxiv.org/abs/1801.07463, equation (16). + These motivate geometric Taylor correction. Strong nodal data, a + lagged fixed point and midpoint flux projection are UW3 adaptations; + neither cited Nitsche/stabilized formulation proves this method stable. + Ruas (2019), "Accuracy enhancement for non-isoparametric finite-element + simulations in curved domains; application to fluid flow", + DOI10.1016/j.camwa.2018.05.029, Sections 2/4, equations (29)-(30), uses + closely related true-boundary P2 trials and canonical polygon tests. + Symmetric stress and the global flux projection remain adaptations. + """ + from underworld3.utilities._cylindrical_dirichlet import _add_cylindrical_zero_slip_transfer + return _add_cylindrical_zero_slip_transfer( + self, boundary=boundary, radius=radius, sign=sign, quadrature=quadrature) + + def solve_cylindrical_zero_slip_updates(self, *, walls, updates): + """Alternate standard Stokes solves and full-vector zero-slip updates. + + Parameters + ---------- + walls : dict + Both boundary labels mapped to handles returned by + :meth:`add_cylindrical_zero_slip_transfer`. + updates : int + Nonnegative number of data updates; the solve count is updates+1. + + Returns + ------- + list of dict + Standard-solver convergence, prescribed polygon flux, update + changes and true-circle nodal/full-trace velocity residuals. + Partial history remains in ``cylindrical_zero_slip_history`` on failure. + + Notes + ----- + Set ``petsc_use_pressure_nullspace=True`` before calling. The supported + homogeneous physical data require zero combined polygon flux; a joint + midpoint-bubble projection enforces this compatibility. Fixed iteration + counts do not guarantee wall convergence: inspect the returned residuals. + This is a strong numerical Taylor iteration, not the weak formulation + analysed by Atallah et al. (2020), DOI10.1016/j.cma.2020.113341. + No field files or analytical solutions are read by this method. + """ + from underworld3.utilities._cylindrical_dirichlet import _solve_cylindrical_zero_slip_updates + return _solve_cylindrical_zero_slip_updates(self, walls=walls, updates=updates) + def recover_cylindrical_topography(self, *, wall, field, buoyancy_scale): r"""Transfer CBF reaction moments to a circular spatial P2 topography. diff --git a/src/underworld3/utilities/_cylindrical_dirichlet.py b/src/underworld3/utilities/_cylindrical_dirichlet.py new file mode 100644 index 000000000..930ee99a1 --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_dirichlet.py @@ -0,0 +1,354 @@ +"""Lagged zero-slip transfer from affine P2 chords to true circular walls. + +Only owned boundary-cell velocity coefficients are replicated. The essential +datum is a separate P2 vector field, registered once and refreshed before each +standard Stokes solve. No pressure, analytical Stokes solution, harmonic fit or +traction recovery enters this correction. Coordinates and fields must be +nondimensional; the mesh must remain fixed after registration. + +The shifted Dirichlet relation is motivated by Atallah, Canuto & Scovazzi +(2020), CMAME 372, 113341, https://doi.org/10.1016/j.cma.2020.113341, +Section 2.2, equations (8)-(13), and its Section 4 Stokes formulation; +see also Burman, Hansbo & Larson (2018), +"A cut finite element method with boundary value correction for the +incompressible Stokes' equations", https://arxiv.org/abs/1801.07463. +Here the complete affine-cell P2 polynomial is continued radially and its +boundary defect is used in a lagged strong essential condition. This differs +from those papers' variational Nitsche formulations; their stability and +convergence results do not establish convergence of this fixed-point scheme. +The joint zero-flux projection below is an additional compatibility condition. + +A closer trial/test-space precedent is Ruas (2019), "Accuracy enhancement for +non-isoparametric finite-element simulations in curved domains; application to +fluid flow", Computers & Mathematics with Applications 77(6), 1756-1769, +https://doi.org/10.1016/j.camwa.2018.05.029, Sections 2/4, equations (29)-(30). +Its P2 velocity trials obey true-boundary point constraints; tests retain zero +polygon traces, with continuous P1 pressure. UW3's symmetric-stress form and +joint flux projection differ. A nonzero projection coefficient perturbs the +true-circle nodal constraint, so the physical residual must be reported. The +paper's regularity and stability results are not claimed for this adaptation. +""" + +import numpy as np +import underworld3 as uw +from mpi4py import MPI +from scipy.spatial import cKDTree + +from .boundary_flux import _boundary_stratum_is + + +def _edge_key(edge, tolerance): + """Coordinate key independent of facet orientation and partition numbering.""" + return tuple(sorted(tuple(row) for row in np.rint(edge / tolerance).astype(np.int64))) + + +def _labelled_edge_keys(mesh, boundary, radius, tolerance): + """Collect only labelled boundary facets, including ghosts on empty ranks.""" + dm = mesh.dm + facets = _boundary_stratum_is(dm, mesh, boundary) + indices = np.empty(0, dtype=int) + if facets is not None: + indices = np.asarray(facets.indices).copy() + facets.destroy() + f0, f1 = dm.getDepthStratum(1) + v0, v1 = dm.getDepthStratum(0) + section = dm.getCoordinateSection() + coordinates = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + local, unsupported = [], False + for facet in indices: + if f0 <= facet < f1: + closure = dm.getTransitiveClosure(int(facet))[0] + vertices = [int(point) for point in closure if v0 <= point < v1] + if len(vertices) != 2: + unsupported = True + else: + local.append([coordinates[section.getOffset(point) // 2] for point in vertices]) + if uw.mpi.comm.allreduce(unsupported, op=MPI.LOR): + raise ValueError("Zero-slip transfer supports two-vertex affine boundary facets only.") + edges = np.concatenate(uw.mpi.comm.allgather(np.asarray(local).reshape(-1, 2, 2))) + if not len(edges) or not np.allclose(np.linalg.norm(edges, axis=-1), radius, + rtol=0, atol=tolerance): + raise ValueError("The labelled boundary must lie on the requested centered circle.") + return {_edge_key(edge, tolerance) for edge in edges} + + +def _p2_velocity(vertices, velocity, points): + """Evaluate Cartesian P2 cell polynomials, including one-sided continuation. + + Array shapes are ``(cells,3,2)``, ``(cells,6,2)`` and ``(cells,points,2)``. + Velocity coefficients are ordered vertices 012, then midpoints 01,12,20. + No radial-frame transformation or angular averaging is applied. + """ + inverse = np.linalg.inv(np.stack( + (vertices[:, 1] - vertices[:, 0], vertices[:, 2] - vertices[:, 0]), axis=-1)) + xi = np.einsum("eij,eqj->eqi", inverse, points - vertices[:, None, 0]) + bary = np.concatenate((1 - xi.sum(axis=-1, keepdims=True), xi), axis=-1) + shape = np.stack( + [bary[..., i] * (2 * bary[..., i] - 1) for i in range(3)] + + [4 * bary[..., i] * bary[..., j] for i, j in ((0, 1), (1, 2), (2, 0))], + axis=-1) + return np.einsum("eqn,end->eqd", shape, velocity) + + +def _zero_slip_datum(vertices, velocity, edge_nodes, radius): + r"""Return ``u_h(x)-u_h(R*x/|x|)`` at endpoint, endpoint, midpoint nodes. + + The prescribed physical velocity is zero. For an affine P2 polynomial this + is the complete quadratic Taylor shift, including its second derivative. + The displacement is radially outward at both walls; outward-domain normal + signs belong only to the flux calculation, never to this vector identity. + """ + radii = np.linalg.norm(edge_nodes, axis=-1) + true_points = radius * edge_nodes / radii[..., None] + datum = (_p2_velocity(vertices, velocity, edge_nodes) + - _p2_velocity(vertices, velocity, true_points)) + # Circle vertices have exactly zero shift; retain that essential value + # despite floating-point error in the radial map. + datum[:, :2] = 0.0 + return datum + + +def _facet_normals(edges, sign): + """Use actual chord normals, oriented radially independently of edge order.""" + direction = edges[:, 1] - edges[:, 0] + normal = np.column_stack((direction[:, 1], -direction[:, 0])) + normal /= np.linalg.norm(normal, axis=1)[:, None] + orientation = np.sign(np.sum(normal * edges.mean(axis=1), axis=1)) + return sign * orientation[:, None] * normal + + +def _p2_facet_flux(edges, values, sign): + """Integrate the P2 vector datum against the constant polygon outward normal. + + ``values`` has shape ``(edges,3,2)`` in endpoint,endpoint,midpoint order. + Simpson weights are exact here. ``sign`` is -1 inner and +1 outer. + """ + normal = _facet_normals(edges, sign) + length = np.linalg.norm(edges[:, 1] - edges[:, 0], axis=1) + average = np.einsum("i,eid->ed", np.array([1, 1, 4]) / 6, values) + return float(np.sum(length * np.sum(average * normal, axis=1))) + + +def _project_joint_flux(walls): + r"""Remove total polygon flux with the smallest L2 normal-midpoint change. + + Each correction is ``-c*N_mid*n_facet``. Its flux is ``-c*2L/3`` and its + squared L2 norm is ``c²*8L/15`` on an edge of length L. Consequently one + coefficient shared by every edge of both walls minimizes the integrated + L2 change in this bubble span. Vertex values and tangential midpoint data + are unchanged. Individual wall fluxes are reported, not independently fit. + + Handles provide ``edges,nodes,node_index,values,sign``. Only ``values`` is + changed, so the mathematical projection is testable without a solver. + """ + raw, lengths = {}, {} + for name, wall in sorted(walls.items()): + raw[name] = _p2_facet_flux(wall.edges, wall.values[wall.node_index], wall.sign) + lengths[name] = float(np.linalg.norm(wall.edges[:, 1] - wall.edges[:, 0], axis=1).sum()) + total_length = sum(lengths.values()) + if not np.isfinite(total_length) or total_length <= 0: + raise ValueError("Joint flux projection requires nonempty circular walls.") + raw_total = sum(raw.values()) + coefficient = raw_total / (2 * total_length / 3) + corrected = {} + for name, wall in sorted(walls.items()): + mid = wall.node_index[:, 2] + if len(np.unique(mid)) != len(mid): + raise ValueError("Each circular boundary midpoint must belong to one edge.") + normal = _facet_normals(wall.edges, wall.sign) + wall.values[mid] -= coefficient * normal + corrected[name] = _p2_facet_flux(wall.edges, wall.values[wall.node_index], wall.sign) + corrected_total = sum(corrected.values()) + tolerance = 128 * np.finfo(float).eps * max(1.0, abs(raw_total)) + if abs(corrected_total) > tolerance: + raise ArithmeticError("Joint P2 zero-slip flux projection failed.") + return dict( + raw_total_flux=float(raw_total), corrected_total_flux=float(corrected_total), + removed_midpoint_normal_value=float(coefficient), + l2_correction_norm=float(abs(coefficient) * np.sqrt(8 * total_length / 15)), + boundaries={name: dict(raw_flux=raw[name], corrected_flux=corrected[name]) + for name in sorted(walls)}, + ) + + +class _ZeroSlipWall: + """Fixed-mesh boundary geometry and one lagged vector essential datum.""" + + def __init__(self, stokes, boundary, radius, sign, quadrature): + mesh, velocity = stokes.mesh, stokes.u + if (mesh.dim != 2 or mesh.cdim != 2 or not mesh.isSimplex + or velocity.degree != 2 or not velocity.continuous): + raise ValueError("Zero-slip transfer requires a two-dimensional continuous affine P2 velocity.") + if stokes.p.degree != 1 or not stokes.p.continuous: + raise ValueError("Zero-slip transfer requires continuous affine P1 regular pressure.") + if uw.is_nondimensional_scaling_active() or uw.get_default_model().has_units_active(): + raise ValueError("Zero-slip transfer supports nondimensional coordinates and fields only.") + if (not np.isfinite(radius) or radius <= 0 or sign not in (-1, 1) + or int(quadrature) != quadrature or quadrature < 3): + raise ValueError("Require finite positive radius, sign +/-1 and integer quadrature >= 3.") + self.stokes, self.mesh, self.velocity = stokes, mesh, velocity + self.name, self.radius, self.sign = boundary, float(radius), sign + self.quadrature = int(quadrature) + tolerance = 1e-10 * radius + boundary_edges = _labelled_edge_keys(mesh, boundary, radius, tolerance) + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, 3) + vertices = np.asarray(mesh._get_coords_for_basis(1, True))[vertex_rows] + owned = mesh._get_owned_cells_mask() + cells, edges = [], [] + for first, second in ((0, 1), (1, 2), (2, 0)): + chosen = np.flatnonzero(owned & np.all(np.isclose( + np.linalg.norm(vertices[:, [first, second]], axis=2), radius, + rtol=0, atol=tolerance), axis=1)) + for cell in chosen: + edge = vertices[cell, [first, second]] + if _edge_key(edge, tolerance) in boundary_edges: + cells.append(cell) + edges.append(edge) + self.cells = np.asarray(cells, dtype=int) + self.v_rows = mesh._cell_node_indices(2, True).reshape(-1, 6)[self.cells] + local_vertices = vertices[self.cells] + standard = np.concatenate((local_vertices, + (local_vertices[:, [0, 1, 2]] + local_vertices[:, [1, 2, 0]]) / 2), axis=1) + distance = np.linalg.norm(np.asarray(velocity.coords)[self.v_rows][:, :, None] + - standard[:, None], axis=-1) + self.v_order = distance.argmin(axis=1) + error = float(np.max(distance.min(axis=1), initial=0.0)) + if uw.mpi.comm.allreduce(error, op=MPI.MAX) > tolerance: + raise ValueError("Zero-slip transfer requires straight-sided affine P2 cells.") + self.vertices = np.concatenate(uw.mpi.comm.allgather(local_vertices)) + self.edges = np.concatenate(uw.mpi.comm.allgather(np.asarray(edges).reshape(-1, 2, 2))) + if not len(self.edges): + raise ValueError("No owned facets lie on the requested circular wall.") + if {_edge_key(edge, tolerance) for edge in self.edges} != boundary_edges: + raise ValueError("Owned boundary-cell traces do not cover the labelled circular wall.") + if np.any(np.linalg.norm(self.edges.mean(axis=1), axis=1) <= tolerance): + raise ValueError("Circular boundary edges must span less than a semicircle.") + self.edge_nodes = np.stack((self.edges[:, 0], self.edges[:, 1], self.edges.mean(axis=1)), axis=1) + self.nodes, index = np.unique(self.edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + self.node_index = index.reshape(-1, 3) + # An exterior closed wall has one owned cell per edge and two edges at + # each vertex. Reject internal circles and incomplete circular arcs. + _, counts = np.unique(np.sort(self.node_index[:, :2], axis=1), axis=0, return_counts=True) + _, incidence = np.unique(self.node_index[:, :2], return_counts=True) + if np.any(counts != 1) or np.any(incidence != 2): + raise ValueError("Transfer requires a complete exterior circular wall, not an interface or arc.") + self.field = uw.discretisation.MeshVariable( + "ZeroSlipShift" + boundary, mesh, 2, degree=2, continuous=True) + distance, index = cKDTree(self.nodes).query(np.asarray(self.field.coords)) + self.local_boundary, self.local_index = distance < tolerance, index + self.values = np.zeros_like(self.nodes) + t, weight = np.polynomial.legendre.leggauss(self.quadrature) + t, weight = (t + 1) / 2, weight / 2 + points = ((1 - t)[None, :, None] * self.edges[:, None, 0] + + t[None, :, None] * self.edges[:, None, 1]) + direction = self.edges[:, 1] - self.edges[:, 0] + cross = np.abs(points[..., 0] * direction[:, None, 1] + - points[..., 1] * direction[:, None, 0]) + self.rule = dict(points=points, + circle_points=radius * points / np.linalg.norm(points, axis=-1)[..., None], + circle_weight=radius * cross / np.sum(points**2, axis=-1) * weight) + + def cell_values(self): + """Collect only owned wall-cell P2 coefficients, including empty ranks.""" + local = np.take_along_axis( + np.asarray(self.velocity.array)[:, 0, :][self.v_rows], + self.v_order[..., None], axis=1) + return np.concatenate(uw.mpi.comm.allgather(local)) + + def update(self): + """Replace raw lagged datum; joint flux projection must follow both walls.""" + datum = _zero_slip_datum(self.vertices, self.cell_values(), self.edge_nodes, self.radius) + length = np.linalg.norm(self.edges[:, 1] - self.edges[:, 0], axis=1) + value, weight = np.zeros_like(self.values), np.zeros(len(self.nodes)) + np.add.at(value, self.node_index.ravel(), (length[:, None, None] * datum).reshape(-1, 2)) + np.add.at(weight, self.node_index.ravel(), np.broadcast_to(length[:, None], datum.shape[:2]).ravel()) + value /= weight[:, None] + change = float(np.max(np.linalg.norm(value - self.values, axis=1))) + self.values = value + return dict(datum_change_maximum=change, + datum_maximum=float(np.max(np.linalg.norm(value, axis=1)))) + + def write(self): + """Write the current compatible datum; the registered expression is fixed.""" + selected = self.local_boundary + self.field.array[selected, 0, :] = self.values[self.local_index[selected]] + + def audit(self): + """Measure the solved polynomial's zero-slip defect on the true circle. + + Full Gauss RMS uses the radial-map circular arclength Jacobian. Nodal + maximum and the polygon essential residual are independent diagnostics; + no pressure-gauge adjustment or trace filtering is applied. + """ + velocity = self.cell_values() + circle = self.radius * self.edge_nodes / np.linalg.norm(self.edge_nodes, axis=-1)[..., None] + nodal = _p2_velocity(self.vertices, velocity, circle) + gauss = _p2_velocity(self.vertices, velocity, self.rule["circle_points"]) + square = float(np.sum(self.rule["circle_weight"] * np.sum(gauss**2, axis=-1))) + measure = float(self.rule["circle_weight"].sum()) + polygon = _p2_velocity(self.vertices, velocity, self.edge_nodes) + return dict( + true_circle_nodal_maximum=float(np.max(np.linalg.norm(nodal, axis=-1))), + true_circle_gauss_maximum=float(np.max(np.linalg.norm(gauss, axis=-1))), + true_circle_velocity_l2=float(np.sqrt(square)), + true_circle_velocity_rms=float(np.sqrt(square / measure)), + true_circle_measure=measure, + polygon_essential_residual_maximum=float(np.max(np.linalg.norm( + polygon - self.values[self.node_index], axis=-1))), + ) + + +def _add_cylindrical_zero_slip_transfer(stokes, *, boundary, radius, sign, quadrature): + """Register one auxiliary P2 vector essential datum before the first solve.""" + if getattr(stokes, "_rotated_freeslip_bcs", None): + raise ValueError("Zero-slip transfer uses standard Stokes, not rotated free slip.") + wall = _ZeroSlipWall(stokes, boundary, radius, sign, quadrature) + stokes.add_dirichlet_bc(wall.field.sym, boundary) + return wall + + +def _solve_cylindrical_zero_slip_updates(stokes, *, walls, updates): + """Solve with jointly compatible, lagged full-vector circular-wall data. + + Both inner and outer handles must be supplied. ``updates`` is an explicit + experimental iteration count, not a convergence guarantee. Each history + row retains the standard solver verdict and the physical-circle velocity + defects so callers can assess fixed-point convergence. No rigid-rotation + removal or torque projection is added here. The caller must enable the + standard pressure nullspace; this helper does not install another gauge. + """ + if int(updates) != updates or updates < 0: + raise ValueError("Zero-slip updates must be a non-negative integer.") + if not stokes.petsc_use_pressure_nullspace: + raise ValueError("Zero-slip updates require petsc_use_pressure_nullspace=True.") + if (len(walls) != 2 or {wall.sign for wall in walls.values()} != {-1, 1} + or any(wall.stokes is not stokes for wall in walls.values())): + raise ValueError("Supply this Stokes system's inner and outer zero-slip handles together.") + radii = {wall.sign: wall.radius for wall in walls.values()} + if not radii[-1] < radii[1]: + raise ValueError("The inner zero-slip radius must be smaller than the outer radius.") + history = [] + stokes.cylindrical_zero_slip_history = history + for iteration in range(int(updates) + 1): + flux = _project_joint_flux(walls) + for _, wall in sorted(walls.items()): + wall.write() + previous = np.asarray(stokes.u.array).copy() + stokes.solve() + report = stokes.solve_report + step = dict(iteration=iteration, flux_projection=flux, + converged=bool(report.converged), reason=int(report.reason), + reason_string=report.reason_str, nonlinear_iterations=int(report.nl_its), + ksp_iterations=int(report.ksp_its), residual_norm=float(report.fnorm), + velocity_change_maximum=uw.mpi.comm.allreduce( + float(np.max(np.abs(stokes.u.array - previous), initial=0.0)), op=MPI.MAX), + boundary_velocity={name: wall.audit() for name, wall in sorted(walls.items())}, + boundary_updates={}) + history.append(step) + if not report.converged or report.bounded or report.deadline_expired or not np.isfinite(report.fnorm): + raise RuntimeError(f"Zero-slip wall iteration {iteration} failed: {report}") + if iteration < updates: + step["boundary_updates"] = {name: wall.update() for name, wall in sorted(walls.items())} + uw.pprint(f"Zero-slip wall iteration {iteration}: velocity change=" + f"{step['velocity_change_maximum']:.3e}") + return history diff --git a/tests/test_1089_cylindrical_dirichlet.py b/tests/test_1089_cylindrical_dirichlet.py new file mode 100644 index 000000000..2591a4064 --- /dev/null +++ b/tests/test_1089_cylindrical_dirichlet.py @@ -0,0 +1,317 @@ +"""Independent polynomial and flux contracts for zero-slip circle correction. + +These controls require no Stokes solve or analytical benchmark input. They test +quadratic continuation and the one global compatibility constraint directly. +""" + +from types import SimpleNamespace + +import numpy as np +import pytest +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] + + +def _ring(count, radius, graded=False): + spacing = 1 + (0.6*np.sin(np.arange(count)*0.73) if graded else np.zeros(count)) + theta = 2*np.pi*np.r_[0, np.cumsum(spacing[:-1])]/spacing.sum() + vertices = radius*np.stack((np.cos(theta), np.sin(theta)), axis=-1) + edges = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + nodes = np.concatenate((vertices, edges.mean(axis=1))) + index = np.stack((np.arange(count), (np.arange(count)+1)%count, + np.arange(count)+count), axis=-1) + return SimpleNamespace(edges=edges, nodes=nodes, node_index=index, + values=np.zeros_like(nodes)) + + +def _cells(wall, sign): + midpoint = wall.edges.mean(axis=1) + radius = np.linalg.norm(wall.edges[:, 0], axis=-1) + interior = midpoint/np.linalg.norm(midpoint, axis=-1)[:, None]*(radius*(1-sign*.18))[:, None] + vertices = np.concatenate((wall.edges, interior[:, None]), axis=1) + nodes = np.concatenate((vertices, (vertices[:, [0, 1, 2]]+vertices[:, [1, 2, 0]])/2), axis=1) + return vertices, nodes + + +def _quadratic(points): + x, y = points[..., 0], points[..., 1] + return np.stack((.7+x*x+2*x*y+3*y, -.2-x*x-2*x*y-y*y+2*x), axis=-1) + + +def _flux_reference(edges, edge_values, sign): + """Exact Simpson integral of a vector quadratic against constant facet normal.""" + tangent = edges[:, 1]-edges[:, 0] + length = np.linalg.norm(tangent, axis=-1) + normal = np.stack((tangent[:, 1], -tangent[:, 0]), axis=-1)/length[:, None] + normal *= (sign*np.sign(np.einsum("ei,ei->e", normal, edges.mean(axis=1))))[:, None] + integral = np.einsum("q,eqi->ei", np.array([1., 1., 4.])/6, edge_values) + return float(np.sum(length*np.einsum("ei,ei->e", integral, normal))) + + +def test_quadratic_continuation_matches_global_cartesian_polynomial_outside_cells(): + from underworld3.utilities._cylindrical_dirichlet import _p2_velocity + + vertices = np.array([[[1.1, .1], [1.4, .1], [1.2, .4]], + [[1.8, -.2], [2.2, -.2], [1.9, .3]]]) + nodes = np.concatenate((vertices, (vertices[:, [0, 1, 2]]+vertices[:, [1, 2, 0]])/2), axis=1) + values = _quadratic(nodes) + original = values.copy() + points = np.array([[[1.22, -.1], [1.8, .7]], [[2.22, .2], [1.6, -.4]]]) + np.testing.assert_allclose(_p2_velocity(vertices, values, points), _quadratic(points), + rtol=2e-13, atol=2e-13) + np.testing.assert_array_equal(values, original) + + +@pytest.mark.parametrize("radius,sign", [(1.22, -1), (2.22, 1)]) +@pytest.mark.parametrize("graded", [False, True]) +def test_zero_slip_datum_is_exact_quadratic_circle_continuation(radius, sign, graded): + from underworld3.utilities._cylindrical_dirichlet import _zero_slip_datum + + wall = _ring(17, radius, graded) + vertices, nodes = _cells(wall, sign) + edge_nodes = wall.nodes[wall.node_index] + values = _quadratic(nodes) + original = values.copy() + actual = _zero_slip_datum(vertices, values, edge_nodes, radius) + physical = radius*edge_nodes/np.linalg.norm(edge_nodes, axis=-1)[..., None] + expected = _quadratic(edge_nodes)-_quadratic(physical) + np.testing.assert_allclose(actual, expected, rtol=0, atol=3e-13) + np.testing.assert_allclose(actual[:, :2], 0, rtol=0, atol=3e-13) + np.testing.assert_array_equal(values, original) + # An independent global quadratic vanishing on this circle has precisely + # its chord value as the shifted datum; no analytical Stokes field is used. + vanishing = (radius**2-np.sum(nodes**2, axis=-1))[..., None]*np.array([.4, -.7]) + exact = (radius**2-np.sum(edge_nodes**2, axis=-1))[..., None]*np.array([.4, -.7]) + np.testing.assert_allclose(_zero_slip_datum(vertices, vanishing, edge_nodes, radius), + exact, rtol=0, atol=3e-13) + + +def test_translation_and_rigid_rotation_have_known_shifted_data(): + from underworld3.utilities._cylindrical_dirichlet import _zero_slip_datum + + wall = _ring(19, 1.7, True) + vertices, nodes = _cells(wall, 1) + edge_nodes = wall.nodes[wall.node_index] + translation = np.broadcast_to([2.3, -.7], nodes.shape).copy() + np.testing.assert_allclose(_zero_slip_datum(vertices, translation, edge_nodes, 1.7), + 0, rtol=0, atol=2e-13) + rotation = np.array([[0., -1.3], [1.3, 0.]]) + physical = 1.7*edge_nodes/np.linalg.norm(edge_nodes, axis=-1)[..., None] + expected = (edge_nodes-physical)@rotation.T + np.testing.assert_allclose(_zero_slip_datum(vertices, nodes@rotation.T+translation, + edge_nodes, 1.7), expected, rtol=0, atol=3e-13) + + +def test_datum_respects_edge_reversal_rotation_scaling_and_cell_permutation(): + from underworld3.utilities._cylindrical_dirichlet import _zero_slip_datum + + wall = _ring(23, 2.22, True) + vertices, nodes = _cells(wall, 1) + edge_nodes = wall.nodes[wall.node_index] + values = _quadratic(nodes) + expected = _zero_slip_datum(vertices, values, edge_nodes, 2.22) + order = np.random.default_rng(31).permutation(len(vertices)) + reversed_cells = [1, 0, 2] + reversed_p2 = [1, 0, 2, 3, 5, 4] + actual = _zero_slip_datum(vertices[order][:, reversed_cells], values[order][:, reversed_p2], + edge_nodes[order][:, [1, 0, 2]], 2.22) + np.testing.assert_allclose(actual, expected[order][:, [1, 0, 2]], rtol=0, atol=4e-13) + angle = .391 + rotation = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) + actual = _zero_slip_datum(vertices@rotation.T, values@rotation.T, + edge_nodes@rotation.T, 2.22) + np.testing.assert_allclose(actual, expected@rotation.T, rtol=0, atol=4e-13) + actual = _zero_slip_datum(vertices*1.7, values, edge_nodes*1.7, 2.22*1.7) + np.testing.assert_allclose(actual, expected, rtol=0, atol=4e-13) + + +def test_exact_facet_flux_obeys_polygon_divergence_theorem(): + from underworld3.utilities._cylindrical_dirichlet import _p2_facet_flux + + actual = 0. + expected = 0. + for count, radius, sign in [(21, 1.22, -1), (29, 2.22, 1)]: + wall = _ring(count, radius, True) + x, y = wall.nodes[..., 0], wall.nodes[..., 1] + values = np.stack((x*x+x, y*y+y), axis=-1)[wall.node_index] + actual += _p2_facet_flux(wall.edges, values, sign) + np.testing.assert_allclose(_p2_facet_flux(wall.edges, values, sign), + _flux_reference(wall.edges, values, sign), rtol=2e-13, atol=2e-13) + a, b = wall.edges[:, 0], wall.edges[:, 1] + cross = a[:, 0]*b[:, 1]-a[:, 1]*b[:, 0] + area = cross.sum()/2 + moment = np.sum((a+b)*cross[:, None], axis=0)/6 + expected += sign*(2*area+2*moment.sum()) # div(x²+x,y²+y)=2x+2y+2. + order = np.random.default_rng(7).permutation(count) + np.testing.assert_allclose(_p2_facet_flux(wall.edges[order, ::-1], values[order][:, [1, 0, 2]], sign), + _flux_reference(wall.edges, values, sign), rtol=2e-13, atol=2e-13) + np.testing.assert_allclose(actual, expected, rtol=2e-13, atol=2e-13) + + +def test_joint_flux_projection_is_one_global_midpoint_constraint_not_wallwise(): + from underworld3.utilities._cylindrical_dirichlet import _project_joint_flux + + walls = dict(inner=_ring(19, 1.22, True), outer=_ring(27, 2.22, True)) + for name, wall in walls.items(): + wall.sign = -1 if name == "inner" else 1 + wall.values[len(wall.edges):] = (.7 if name == "inner" else -.2)*wall.nodes[len(wall.edges):] + originals = {name: wall.values.copy() for name, wall in walls.items()} + raw = sum(_flux_reference(w.edges, w.values[w.node_index], w.sign) for w in walls.values()) + length = sum(np.linalg.norm(w.edges[:, 1]-w.edges[:, 0], axis=-1).sum() for w in walls.values()) + coefficient = raw/(2*length/3) + check = _project_joint_flux(walls) + np.testing.assert_allclose(check["raw_total_flux"], raw, rtol=2e-13, atol=2e-13) + np.testing.assert_allclose(check["removed_midpoint_normal_value"], coefficient, rtol=2e-13, atol=2e-13) + np.testing.assert_allclose(check["l2_correction_norm"], abs(coefficient)*np.sqrt(8*length/15), + rtol=2e-13, atol=2e-13) + for name, wall in walls.items(): + count = len(wall.edges) + expected = originals[name].copy() + midpoint = wall.nodes[count:] + expected[count:] -= coefficient*wall.sign*midpoint/np.linalg.norm(midpoint, axis=-1)[:, None] + np.testing.assert_allclose(wall.values, expected, rtol=0, atol=2e-13) + np.testing.assert_array_equal(wall.values[:count], 0) + assert abs(_flux_reference(wall.edges, wall.values[wall.node_index], wall.sign)) > 1e-4 + corrected = sum(_flux_reference(w.edges, w.values[w.node_index], w.sign) for w in walls.values()) + np.testing.assert_allclose(corrected, 0, rtol=0, atol=2e-13) + np.testing.assert_allclose(check["corrected_total_flux"], 0, rtol=0, atol=2e-13) + + +def test_zero_shift_and_zero_joint_flux_need_no_field_correction(): + from underworld3.utilities._cylindrical_dirichlet import _project_joint_flux, _zero_slip_datum + + wall = _ring(16, 1.22) + vertices, nodes = _cells(wall, -1) + # Nodes placed on the true circle have exactly zero displacement. + points = wall.nodes[wall.node_index] + points = 1.22*points/np.linalg.norm(points, axis=-1)[..., None] + np.testing.assert_allclose(_zero_slip_datum(vertices, _quadratic(nodes), points, 1.22), + 0, rtol=0, atol=3e-13) + wall.sign = -1 + original = wall.values.copy() + check = _project_joint_flux(dict(inner=wall)) + np.testing.assert_array_equal(wall.values, original) + assert check["raw_total_flux"] == 0 and check["corrected_total_flux"] == 0 + + +@pytest.mark.parametrize("layout", ["empty_other_ranks", "partitioned_with_ghosts"]) +def test_boundary_adapter_collects_owned_coefficients_and_preserves_velocity(layout, monkeypatch): + from underworld3.utilities import _cylindrical_dirichlet as module + _ZeroSlipWall = module._ZeroSlipWall + + radius = 2.22 + ring = _ring(23, radius, True) + # A valid polygon triangulation with three circle vertices in each cell + # includes internal diagonals joining circle vertices. Only labelled exterior + # edges belong to this wall; a radius-only selector would incorrectly fail. + polygon = ring.edges[:, 0] + all_vertices = np.stack([np.stack((polygon[0], polygon[j], polygon[j+1])) + for j in range(1, len(polygon)-1)]) + all_nodes = np.concatenate((all_vertices, + (all_vertices[:, [0, 1, 2]]+all_vertices[:, [1, 2, 0]])/2), axis=1) + if layout == "empty_other_ranks": + selected = np.arange(len(all_vertices)) if uw.mpi.rank == 0 else np.array([], dtype=int) + owned = np.ones(len(selected), dtype=bool) + else: + chosen = np.arange(len(all_vertices))[uw.mpi.rank::uw.mpi.size] + ghosts = np.setdiff1d(np.arange(3), chosen) + selected = np.r_[chosen, ghosts] + owned = np.r_[np.ones(len(chosen), dtype=bool), np.zeros(len(ghosts), dtype=bool)] + local_vertices, local_nodes = all_vertices[selected], all_nodes[selected] + coordinates, inverse = np.unique(local_nodes.reshape(-1, 2), axis=0, return_inverse=True) + rows = inverse.reshape(-1, 6) + # Deliberately use a noncanonical local P2 ordering. The adapter must discover + # it geometrically while excluding every ghost cell from the gather. + rows = rows[:, [3, 1, 5, 2, 0, 4]] + p1_coordinates, p1_inverse = np.unique(local_vertices.reshape(-1, 2), axis=0, return_inverse=True) + mesh = SimpleNamespace(dim=2, cdim=2, isSimplex=True, boundaries=[SimpleNamespace(name="Outer")], + _cell_node_indices=lambda degree, continuous: p1_inverse.reshape(-1, 3) if degree == 1 else rows, + _get_coords_for_basis=lambda degree, continuous: p1_coordinates, + _get_owned_cells_mask=lambda: owned) + velocity = SimpleNamespace(mesh=mesh, degree=2, continuous=True, + coords=coordinates, array=_quadratic(coordinates)[:, None, :]) + original = np.asarray(velocity.array).copy() + + def auxiliary(*args, **kwargs): + return SimpleNamespace(mesh=mesh, coords=coordinates, + array=np.full_like(velocity.array, 777.), sym=None) + + # PETSc label extraction has a separate real-mesh test; here supply its + # known fixture boundary so the owned/ghost coefficient adapter is isolated. + monkeypatch.setattr(module, "_labelled_edge_keys", lambda mesh, name, radius, tolerance: + {module._edge_key(edge, tolerance) for edge in ring.edges}) + monkeypatch.setattr(uw.discretisation, "MeshVariable", auxiliary) + monkeypatch.setattr(uw, "is_nondimensional_scaling_active", lambda: False) + monkeypatch.setattr(uw, "get_default_model", lambda: SimpleNamespace(has_units_active=lambda: False)) + wall = _ZeroSlipWall(SimpleNamespace(mesh=mesh, u=velocity, p=SimpleNamespace(degree=1, continuous=True)), "Outer", radius, 1, 12) + assert len(wall.edges) == len(ring.edges) + np.testing.assert_allclose(wall.cell_values(), _quadratic(np.concatenate((wall.vertices, + (wall.vertices[:, [0, 1, 2]]+wall.vertices[:, [1, 2, 0]])/2), axis=1)), rtol=0, atol=1e-13) + wall.update() + circle = radius*wall.nodes/np.linalg.norm(wall.nodes, axis=-1)[:, None] + expected = _quadratic(wall.nodes)-_quadratic(circle) + is_vertex = np.isclose(np.linalg.norm(wall.nodes, axis=-1), radius, rtol=0, atol=1e-12) + expected[is_vertex] = 0 + np.testing.assert_allclose(wall.values, expected, rtol=0, atol=5e-13) + wall.write() + np.testing.assert_allclose(wall.field.array[wall.local_boundary, 0], + expected[wall.local_index[wall.local_boundary]], rtol=0, atol=5e-13) + np.testing.assert_array_equal(wall.field.array[~wall.local_boundary, 0], 777.) + np.testing.assert_array_equal(velocity.array, original) + check = wall.audit() + theta = 2*np.pi*np.arange(32)/32 + exact = _quadratic(radius*np.stack((np.cos(theta), np.sin(theta)), axis=-1)) + rms = np.sqrt(np.mean(np.sum(exact**2, axis=-1))) + np.testing.assert_allclose(check["true_circle_velocity_rms"], rms, rtol=3e-13, atol=3e-13) + np.testing.assert_allclose(check["true_circle_measure"], 2*np.pi*radius, rtol=3e-13, atol=3e-13) + np.testing.assert_array_equal(velocity.array, original) + + +def test_real_mesh_labelled_walls_continue_exact_polynomial_without_a_solve(): + from underworld3.utilities._cylindrical_dirichlet import _ZeroSlipWall, _project_joint_flux + + mesh = uw.meshing.Annulus(radiusInner=.5, radiusOuter=1., cellSize=1/8, qdegree=2) + velocity = uw.discretisation.MeshVariable("DirichletPolynomialWitness", mesh, 2, degree=2, continuous=True) + velocity.array[:, 0, :] = _quadratic(np.asarray(velocity.coords)) + original = np.asarray(velocity.array).copy() + stokes = SimpleNamespace(mesh=mesh, u=velocity, p=SimpleNamespace(degree=1, continuous=True)) + walls = {name: _ZeroSlipWall(stokes, name, radius, sign, 12) + for name, radius, sign in [("Lower", .5, -1), ("Upper", 1., 1)]} + for wall in walls.values(): + wall.update() + np.testing.assert_allclose(wall.cell_values(), _quadratic(np.concatenate((wall.vertices, + (wall.vertices[:, [0, 1, 2]]+wall.vertices[:, [1, 2, 0]])/2), axis=1)), rtol=0, atol=3e-13) + check = wall.audit() + theta = 2*np.pi*np.arange(32)/32 + exact = _quadratic(wall.radius*np.stack((np.cos(theta), np.sin(theta)), axis=-1)) + np.testing.assert_allclose(check["true_circle_velocity_rms"], + np.sqrt(np.mean(np.sum(exact**2, axis=-1))), rtol=4e-13, atol=4e-13) + projection = _project_joint_flux(walls) + assert abs(projection["corrected_total_flux"]) < 2e-13 + for wall in walls.values(): + wall.write() + np.testing.assert_allclose(wall.field.array[wall.local_boundary, 0], + wall.values[wall.local_index[wall.local_boundary]], rtol=0, atol=3e-13) + np.testing.assert_array_equal(velocity.array, original) + + +@pytest.mark.parametrize("degree,continuous", [(0, False), (1, False), (2, True)]) +def test_adapter_rejects_unsupported_regular_pressure_before_collectives(degree, continuous): + from underworld3.utilities._cylindrical_dirichlet import _ZeroSlipWall + + mesh = SimpleNamespace(dim=2, cdim=2, isSimplex=True) + stokes = SimpleNamespace(mesh=mesh, u=SimpleNamespace(degree=2, continuous=True), + p=SimpleNamespace(degree=degree, continuous=continuous)) + with pytest.raises(ValueError, match="continuous affine P1"): + _ZeroSlipWall(stokes, "Outer", 1., 1, 12) + + +def test_lagged_solve_requires_the_standard_pressure_nullspace(): + from underworld3.utilities._cylindrical_dirichlet import _solve_cylindrical_zero_slip_updates + + # Fail before any field modification, collective or solve; full Dirichlet + # Stokes has the standard constant-pressure kernel. + stokes = SimpleNamespace(petsc_use_pressure_nullspace=False) + with pytest.raises(ValueError, match="petsc_use_pressure_nullspace=True"): + _solve_cylindrical_zero_slip_updates(stokes, walls={}, updates=1) From b198af164f8311528c599e794973f616c16d773f Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 25 Oct 2026 15:31:10 +1100 Subject: [PATCH 5/6] Guard null PETSc boundary index sets on empty wall ranks The first16-rank zero-slip installation gate exposed a real null-handle dereference in labelled edge extraction: PETSc may return a false IS wrapper rather than PythonNone for an absent local stratum. Use the established PETSc truthiness guard before indices/destroy. Existing stratum consumers already use this contract; no generic boundary-flux behavior or numerical wall data changes. Add real-null, safe false-sentinel and valid zero-length index-set regressions. The exact prior commit fails the safe sentinel at the same indices access; corrected20-test helper suite passes serially and on every MPI2/MPI16 rank through explicit development overlays. Preserve failed Gadi build182210837 with no validated-install marker and no model submissions. Next: fresh versioned single-revision builds and complete native gates before MPI model use. Underworld development team with AI support from Claude Code --- .../utilities/_cylindrical_dirichlet.py | 4 +- tests/test_1089_cylindrical_dirichlet.py | 56 +++++++++++++++++++ 2 files changed, 59 insertions(+), 1 deletion(-) diff --git a/src/underworld3/utilities/_cylindrical_dirichlet.py b/src/underworld3/utilities/_cylindrical_dirichlet.py index 930ee99a1..3c14e0d38 100644 --- a/src/underworld3/utilities/_cylindrical_dirichlet.py +++ b/src/underworld3/utilities/_cylindrical_dirichlet.py @@ -47,7 +47,9 @@ def _labelled_edge_keys(mesh, boundary, radius, tolerance): dm = mesh.dm facets = _boundary_stratum_is(dm, mesh, boundary) indices = np.empty(0, dtype=int) - if facets is not None: + # PETSc may return a false IS wrapper with a null C handle on empty ranks. + # Checking only Python None would dereference that null handle at .indices. + if facets: indices = np.asarray(facets.indices).copy() facets.destroy() f0, f1 = dm.getDepthStratum(1) diff --git a/tests/test_1089_cylindrical_dirichlet.py b/tests/test_1089_cylindrical_dirichlet.py index 2591a4064..6de82bb04 100644 --- a/tests/test_1089_cylindrical_dirichlet.py +++ b/tests/test_1089_cylindrical_dirichlet.py @@ -315,3 +315,59 @@ def test_lagged_solve_requires_the_standard_pressure_nullspace(): stokes = SimpleNamespace(petsc_use_pressure_nullspace=False) with pytest.raises(ValueError, match="petsc_use_pressure_nullspace=True"): _solve_cylindrical_zero_slip_updates(stokes, walls={}, updates=1) + + +@pytest.mark.parametrize("local_stratum", ["null_petsc_is", "false_sentinel", "valid_empty_is"]) +def test_absent_local_boundary_is_not_dereferenced_before_global_edge_gather(local_stratum, monkeypatch): + """PETSc.IS() has a false null handle, distinct from a valid zero-length IS. + + The Gadi MPI16 gate first exposed the real null-handle crash. Never inspect + its .indices here; a separate false sentinel safely makes any such access + an assertion failure instead of triggering the native segmentation fault. + """ + from petsc4py import PETSc + from underworld3.utilities import _cylindrical_dirichlet as module + + class FalseHandle: + def __bool__(self): + return False + + @property + def indices(self): + raise AssertionError("A false PETSc handle must not be dereferenced.") + + def destroy(self): + raise AssertionError("A false PETSc handle must not be destroyed.") + + if local_stratum == "null_petsc_is": + facets = PETSc.IS() + assert not facets + elif local_stratum == "valid_empty_is": + facets = PETSc.IS().createGeneral(np.empty(0, dtype=PETSc.IntType), comm=PETSc.COMM_SELF) + assert bool(facets) # Handle validity, not the length of its index array. + else: + facets = FalseHandle() + radius = 1.7 + remote_edges = _ring(11, radius, True).edges + + class OtherRankOwnsWall: + def allreduce(self, value, op): + assert value is False + return False + + def allgather(self, local): + np.testing.assert_array_equal(local.shape, (0, 2, 2)) + return [local, remote_edges.copy()] + + empty_dm = SimpleNamespace( + getDepthStratum=lambda depth: (0, 0), + getCoordinateSection=lambda: SimpleNamespace(), + getCoordinatesLocal=lambda: SimpleNamespace(array=np.empty(0)), + ) + monkeypatch.setattr(module, "_boundary_stratum_is", lambda dm, mesh, boundary: facets) + monkeypatch.setattr(uw.mpi, "comm", OtherRankOwnsWall()) + tolerance = 1e-10*radius + actual = module._labelled_edge_keys(SimpleNamespace(dm=empty_dm), "Outer", radius, tolerance) + assert actual == {module._edge_key(edge, tolerance) for edge in remote_edges} + if local_stratum == "valid_empty_is": + assert not facets # The valid empty IS is released; null handles are untouched. From 1e2d0a5e96a72c8628d912f4a06d817838836308 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sun, 25 Oct 2026 16:04:46 +1100 Subject: [PATCH 6/6] Recover full zero-slip annulus traction on physical circles Transfer freshly assembled Cartesian regular-pressure reactions using physical-circle and polygon stress moments, then solve a rotating normal/tangent P2 trace system. Retain circular and weighted-reconstructed fields separately, without wall pressure shifts. Add exact constant pressure, affine stress/lift, conservation, gauge, and real empty-owner contracts. Python helper gates pass serial/MPI2/MPI4 and saved-state prototype parity; fresh native wrapper builds and full refinement acceptance remain pending. Underworld development team with AI support from Claude Code --- .../cylindrical-geometry-correction.md | 64 +++++ .../cython/petsc_generic_snes_solvers.pyx | 74 +++++ .../_cylindrical_dirichlet_traction.py | 193 +++++++++++++ ...est_1089_cylindrical_dirichlet_traction.py | 269 ++++++++++++++++++ 4 files changed, 600 insertions(+) create mode 100644 src/underworld3/utilities/_cylindrical_dirichlet_traction.py create mode 100644 tests/test_1089_cylindrical_dirichlet_traction.py diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md index 41de27f98..78a2c148a 100644 --- a/docs/advanced/cylindrical-geometry-correction.md +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -128,6 +128,63 @@ numerical fixed point, especially for the delta case's piecewise regularity. Independent polynomial/flux/parallel tests pass; model acceptance remains pending in the benchmark investigation. +### Full Dirichlet traction recovery + +Standard Stokes can recover the full Cartesian traction from essential-boundary +reactions with `boundary_flux(..., mass="consistent", remove_mean=False)`. +Its fresh volume-residual assembly uses the current pressure gauge and registered +essential data. For saved-state replay, restore both the solved velocity/pressure +and the last `ZeroSlipShiftLower/Upper` auxiliary fields; replacing those data by +zero or by another wall update changes the state being measured. + +`recover_cylindrical_zero_slip_traction` transfers the regular-pressure reaction +to the physical circle using one-sided numerical stress: + +$$R_c=R_q+L_c(\sigma_{q+p_J})-L_h(\sigma_q).$$ + +The polygon subtraction must use **regular** stress when the solved unknown is +$q$. A prescribed lift can vanish on the true wall while remaining nonzero on +the chords. Using physical stress in that subtraction without first converting +the regular reaction leaves a spurious load. Circle continuation uses the complete +affine-cell P2/P1 polynomial, rather than an analytical Stokes solution. + +The line system has Cartesian P2 test components and rotating trial vectors +$[n,e_\theta]$ multiplying scalar P2 normal/tangential coefficients. Thus it +preserves every transferred Cartesian weak moment and represents a constant +normal pressure exactly. This avoids contracting reactions with nodal normals. +The public shear target is $\sigma_{r\theta}$; the outward tangential traction +has the additional inner/outer `sign`. No wall mean is removed. + +```python +checks = stokes.recover_cylindrical_zero_slip_traction( + wall=walls["Upper"], normal_field=normal_circle, shear_field=shear_circle, + recovered_normal_field=normal_recovered, recovered_shear_field=shear_recovered, + viscosity=1.0, pressure_jump=pressure_jump, degree=6, + halfwidth_edges=5, quadrature=24, +) +``` + +All four targets are distinct continuous scalar P2 fields on the Stokes mesh. +Call after the global pressure calibration. The optional spatial reconstruction +is applied separately to the two rotating scalar traces, preserving each mean +and first angular moments and therefore their combined Cartesian resultant. +It changes finer transferred moments; the circular CBF fields remain available. +Returned checks distinguish the original polygon resultant, geometry-transfer +change, transferred resultant, and subsequent reconstruction changes. The +kernel's physical bandwidth limits described below still apply to shear. + +For an exact stationary incompressible zero-slip solution, $u=0$ along the +circle implies $\partial_r u_r=0$, so $p=-\sigma_{nn}$. This gives a derived +**boundary pressure estimate**, not an improved FE pressure field. Keep the +original pressure trace and the defect $2\mu\partial_r u_r$ in verification; +weak incompressibility does not make that identity exact for the discrete field. + +Brunátová et al. (2025), Section 2.3.2, equations (9)-(10), gives a direct +variational precedent for full Dirichlet traction and its tangential decomposition. +Its Stokes example uses gradient pseudostress, while UW3 uses symmetric stress. +Neither it nor the CBF reference proves this circular-transfer/rotating-trial/ +weighted-reconstruction adaptation stable or superconvergent. + ## Spatial circular reaction and topography The circular reaction moments are @@ -240,6 +297,13 @@ independently for a scientific before/after comparison. ## References and attribution boundaries +- Brunátová, J., Dokken, J. S., Valen-Sendstad, K., and Hron, J. (2025), + *On the Numerical Evaluation of Wall Shear Stress Using the Finite Element + Method*, International Journal for Numerical Methods in Biomedical Engineering + 41(9), e70086, [doi:10.1002/cnm.70086](https://doi.org/10.1002/cnm.70086). + Section 2.3.2 supplies the full Dirichlet variational recovery class and + normal/tangential decomposition. UW3's circular moments, symmetric Stokes + form, rotating trial basis and spatial reconstruction are separate adaptations. - Chang, K. S., and Kwak, D. Y. (2009), *Extraction method for Stokes Flow with jumps in the pressure*, [arXiv:0911.4769](https://arxiv.org/abs/0911.4769). Section 2, equations (4)-(6), motivates regular/singular pressure splitting. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c547f2710..68f1489be 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8883,6 +8883,80 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): from underworld3.utilities._cylindrical_dirichlet import _solve_cylindrical_zero_slip_updates return _solve_cylindrical_zero_slip_updates(self, walls=walls, updates=updates) + def recover_cylindrical_zero_slip_traction(self, *, wall, normal_field, shear_field, + recovered_normal_field, recovered_shear_field, viscosity, pressure_jump, + degree, halfwidth_edges, quadrature): + r"""Recover physical normal/shear stress at a fully prescribed circle. + + Parameters + ---------- + wall + Handle returned by :meth:`add_cylindrical_zero_slip_transfer`. + Its registered auxiliary datum must still describe the solved state. + normal_field, shear_field : MeshVariable + Distinct continuous scalar P2 targets for the circular CBF + :math:`\sigma_{nn}` and polar shear :math:`\sigma_{r\theta}`. + recovered_normal_field, recovered_shear_field : MeshVariable + Two further scalar P2 targets for spatially reconstructed traces. + The preceding circular CBF fields are retained unchanged. + viscosity : float + Positive constant viscosity matching the Stokes constitutive law. + Augmentation and DEVSS must be disabled. + pressure_jump : sympy expression + Prescribed source-pressure lift :math:`p_J`, with physical pressure + :math:`p=q+p_J`. Pass exact zero when the solved pressure is physical. + This is a source profile, never an analytical Stokes pressure. + degree : {4, 6, 8} + Weighted reconstruction-kernel degree; targets remain P2. + halfwidth_edges : float + Kernel support halfwidth in mean circular edge arc lengths. + quadrature : int + Gauss points per edge, at least three. Assess quadrature convergence. + + Returns + ------- + dict + Cartesian weak-moment, absolute pressure mean, resultant and reconstruction + checks. No wall mean is removed. Fine CBF moments change only in + the two reconstructed targets; mean/resultants are preserved. + + Notes + ----- + Call collectively after pressure calibration. Standard ``boundary_flux`` + reassembles the current regular-pressure vector reaction. Numerical + one-sided P2/P1 stress moments transfer it from chords to true circles: + :math:`R_c=R_q+L_c(\sigma_{q+p_J})-L_h(\sigma_q)`. + Cartesian tests and rotating normal/tangent P2 trials preserve constant + normal pressure without nodal-normal contraction. Only boundary traces + are replicated; small sparse line systems are solved on each rank. + + Supported scope is fixed nondimensional affine 2D continuous P2/P1 + constant-viscosity, fully Dirichlet Stokes. No Stokes solve, analytical + reference, pressure gauge change or volume-stress projection is performed. + At an exact stationary incompressible zero-slip wall, :math:`p=-\sigma_{nn}`; + this derived boundary pressure estimate does not change the FE pressure. + + References + ---------- + Zhong, Gurnis and Hulbert (1993), PEPI 78, 1-8, + DOI10.1016/0031-9201(93)90078-N, provides the CBF principle. + Brunatova, Dokken, Valen-Sendstad and Hron (2025), + DOI10.1002/cnm.70086, Section 2.3.2, equations (9)-(10), describes + variational recovery for fully Dirichlet boundaries and tangential + decomposition. Its Stokes example uses pseudostress; our Stokes form + uses full symmetric stress. Circular transfer, rotating trials and + weighted reconstruction are UW3 adaptations, without an inherited + convergence theorem. Reconstruction context: Docampo-Sanchez et al. + (2020), DOI10.1016/j.compfluid.2020.104727. + """ + from underworld3.utilities._cylindrical_dirichlet_traction import _recover_cylindrical_zero_slip_traction + return _recover_cylindrical_zero_slip_traction(self, wall=wall, + normal_field=normal_field, shear_field=shear_field, + recovered_normal_field=recovered_normal_field, + recovered_shear_field=recovered_shear_field, viscosity=viscosity, + pressure_jump=pressure_jump, degree=degree, + halfwidth_edges=halfwidth_edges, quadrature=quadrature) + def recover_cylindrical_topography(self, *, wall, field, buoyancy_scale): r"""Transfer CBF reaction moments to a circular spatial P2 topography. diff --git a/src/underworld3/utilities/_cylindrical_dirichlet_traction.py b/src/underworld3/utilities/_cylindrical_dirichlet_traction.py new file mode 100644 index 000000000..f0f9f1ace --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_dirichlet_traction.py @@ -0,0 +1,193 @@ +"""Physical circular traction from standard fully Dirichlet Stokes reactions. + +The vector consistent boundary flux (CBF) is retained as a weak load, corrected +by circle-minus-polygon numerical stress moments, then represented in a rotating +normal/tangent P2 trial basis. No pressure mean is removed and no reference +Stokes solution is used. Only boundary traces are replicated; small SciPy sparse +line systems follow the existing circular reaction/reconstruction implementation. + +Zhong, Gurnis & Hulbert (1993), doi:10.1016/0031-9201(93)90078-N, motivates CBF. +Brunatova, Dokken, Valen-Sendstad & Hron (2025), doi:10.1002/cnm.70086, +Section 2.3.2, equations (9)-(10), describes full Dirichlet variational traction +recovery and normal/tangential decomposition. Those works do not establish +the circular moment transfer, rotating trial basis or reconstruction used here. +""" + +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.spatial import cKDTree +import sympy as sp +import underworld3 as uw +from underworld3.systems.solvers import _penalty_value + +from ._cylindrical_dirichlet import _ZeroSlipWall, _facet_normals +from ._cylindrical_reactions import _validate_spatial_target +from ._cylindrical_reconstruction import _recover_circular_trace +from ._cylindrical_wall import _numerical_trace, _trace_rule + + +def _traction_rule(edges, radius, sign, quadrature): + """Radially map polygon P2 tests onto the physical circle, with exact ds.""" + polygon = _trace_rule(edges, sign, 0, quadrature) + polygon_normal = _facet_normals(edges, sign) + points = polygon["points"] + radial = points/np.linalg.norm(points, axis=-1)[..., None] + circle = radius*radial + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + t = (np.polynomial.legendre.leggauss(quadrature)[0]+1)/2 + shape = np.stack(((1-t)*(1-2*t), t*(2*t-1), 4*t*(1-t)), axis=-1) + edge_nodes = np.concatenate((edges, edges.mean(axis=1)[:, None]), axis=1) + nodes, index = np.unique(edge_nodes.reshape(-1, 2), axis=0, return_inverse=True) + normal = sign*radial + circle_weight = polygon["weight"]*radius/np.linalg.norm(points, axis=-1)*np.abs( + np.einsum("ed,eqd->eq", polygon_normal, radial)) + return dict(points=points, circle_points=circle, nodes=nodes, index=index.reshape(-1, 3), + shape=shape, polygon_weight=polygon["weight"], circle_weight=circle_weight, + polygon_normal=polygon_normal, normal=normal, + frame=np.stack((normal, tangent), axis=-1), radius=radius, sign=sign) + + +def _regular_cartesian_stress(vertices, velocity, pressure, points, viscosity): + """Continue the affine-cell P2/P1 regular-pressure stress one-sidedly.""" + trace = _numerical_trace(vertices, velocity, pressure, points, viscosity, + sheet_radius=0, outer_radius=1, mode=0, amplitude=0) + radial = points/np.linalg.norm(points, axis=-1)[..., None] + tangent = np.stack((-radial[..., 1], radial[..., 0]), axis=-1) + return (trace["rr"][..., None, None]*np.einsum("eqi,eqj->eqij", radial, radial) + + trace["tt"][..., None, None]*np.einsum("eqi,eqj->eqij", tangent, tangent) + + trace["rt"][..., None, None]*(np.einsum("eqi,eqj->eqij", radial, tangent) + + np.einsum("eqi,eqj->eqij", tangent, radial))) + + +def _transfer_circle_traction(rule, raw_traction, polygon_stress, circle_stress): + r"""Preserve every transferred Cartesian moment in a rotating P2 trial space. + + ``raw_traction`` is the physical-sign vector CBF for the regular pressure q. + Subtract polygon regular-stress moments and add circle physical-stress + moments: :math:`R_c=R_q+L_c(\sigma_{q+p_J})-L_h(\sigma_q)`. + The two unknowns are :math:`\sigma_{nn}` and :math:`t\cdot e_\theta`; + the latter equals ``sign*sigma_rtheta``. No mean or modal fit is applied. + """ + index, shape, size = rule["index"], rule["shape"], len(rule["nodes"]) + local_mass = np.einsum("eq,qi,qj->eij", rule["polygon_weight"], shape, shape) + rows = np.broadcast_to(index[:, :, None], local_mass.shape).ravel() + columns = np.broadcast_to(index[:, None, :], local_mass.shape).ravel() + mass = coo_matrix((local_mass.ravel(), (rows, columns)), shape=(size, size)).tocsr() + polygon_load = np.einsum("eq,qi,eqcd,ed->eic", rule["polygon_weight"], shape, + polygon_stress, rule["polygon_normal"]) + circle_load = np.einsum("eq,qi,eqcd,eqd->eic", rule["circle_weight"], shape, + circle_stress, rule["normal"]) + raw_load = mass@raw_traction + target = raw_load.copy() + np.add.at(target, index.ravel(), (circle_load-polygon_load).reshape(-1, 2)) + # Rows are Cartesian test components; columns are normal/tangent trials. + local = np.einsum("eq,qi,qj,eqck->eicjk", rule["circle_weight"], shape, shape, rule["frame"]) + rows = np.broadcast_to(2*index[:, :, None, None, None]+np.arange(2)[None, None, :, None, None], local.shape) + columns = np.broadcast_to(2*index[:, None, None, :, None]+np.arange(2)[None, None, None, None, :], local.shape) + circular_mass = coo_matrix((local.ravel(), (rows.ravel(), columns.ravel())), shape=(2*size, 2*size)).tocsr() + values = spsolve(circular_mass, target.ravel()).reshape(-1, 2) + if not np.all(np.isfinite(values)): + raise ArithmeticError("Circular rotating-frame traction solve returned nonfinite values.") + sampled = np.einsum("qi,eik->eqk", shape, values[index]) + traction = np.einsum("eqk,eqck->eqc", sampled, rule["frame"]) + resultant = np.einsum("eq,eqc->c", rule["circle_weight"], traction) + checks = dict(maximum_cartesian_moment_residual=float(np.max(np.abs( + circular_mass@values.ravel()-target.ravel()))), + cartesian_resultant=resultant.tolist(), cartesian_resultant_target=target.sum(axis=0).tolist(), + raw_polygon_cartesian_resultant=raw_load.sum(axis=0).tolist(), + circle_transfer_resultant_change=(target-raw_load).sum(axis=0).tolist(), + polygon_numerical_stress_resultant=polygon_load.sum(axis=(0, 1)).tolist(), + circle_numerical_stress_resultant=circle_load.sum(axis=(0, 1)).tolist(), + cartesian_resultant_residual=float(np.max(np.abs(resultant-target.sum(axis=0)))), + circular_length=float(rule["circle_weight"].sum()), pressure_mean_removed=False, + normal_mean=float(np.sum(rule["circle_weight"]*sampled[..., 0])/(2*np.pi*rule["radius"])), + tangential_mean=float(np.sum(rule["circle_weight"]*sampled[..., 1])/(2*np.pi*rule["radius"]))) + return values, checks + + +def _trace_change_checks(rule, source, target): + """Measure changed nodal loads separately from conserved physical moments.""" + delta = np.einsum("qi,eik->eqk", rule["shape"], (target-source)[rule["index"]]) + traction = np.einsum("eqk,eqck->eqc", delta, rule["frame"]) + loads = np.zeros((len(source), 2)) + local = np.einsum("eq,qi,eqc->eic", rule["circle_weight"], rule["shape"], traction) + np.add.at(loads, rule["index"].ravel(), local.reshape(-1, 2)) + return dict(maximum_changed_cartesian_nodal_moment=float(np.max(np.abs(loads))), + maximum_cartesian_resultant_change=float(np.max(np.abs(np.einsum( + "eq,eqc->c", rule["circle_weight"], traction)))), + scalar_integral_changes=np.einsum("eq,eqk->k", rule["circle_weight"], delta).tolist(), + fine_scale_cbf_moments_preserved=False) + + +def _boundary_pressure_values(stokes, wall): + """Gather owned boundary-cell P1 values, including ranks with no wall cells.""" + rows = stokes.mesh._cell_node_indices(1, True).reshape(-1, 3)[wall.cells] + local = np.asarray(stokes.p.array)[rows, 0, 0] + return np.concatenate(uw.mpi.comm.allgather(local)) + + +def _recover_cylindrical_zero_slip_traction(stokes, *, wall, normal_field, shear_field, + recovered_normal_field, recovered_shear_field, viscosity, pressure_jump, + degree, halfwidth_edges, quadrature): + """Write distinct raw-circle/reconstructed physical normal and shear fields.""" + if not isinstance(wall, _ZeroSlipWall) or wall.stokes is not stokes: + raise ValueError("Supply this Stokes system's zero-slip wall handle.") + if getattr(stokes, "_rotated_freeslip_bcs", None): + raise ValueError("Full Dirichlet traction recovery requires standard Stokes.") + if not np.isfinite(viscosity) or viscosity <= 0: + raise ValueError("Traction recovery requires finite positive constant viscosity.") + constitutive_viscosity = uw.function.unwrap(stokes.constitutive_model.K, + keep_constants=False, return_self=False) + if (_penalty_value(stokes.penalty) != 0 or stokes.devss_viscosity is not None + or constitutive_viscosity.is_number is not True + or float(constitutive_viscosity) != viscosity): + raise ValueError("Require matching constant viscosity and zero augmentation/DEVSS.") + if int(quadrature) != quadrature or quadrature < 3: + raise ValueError("Traction quadrature requires an integer number of points >= 3.") + targets = (normal_field, shear_field, recovered_normal_field, recovered_shear_field) + if len({id(field) for field in targets}) != 4: + raise ValueError("Circular and reconstructed normal/shear targets must be four distinct fields.") + for field in targets: + _validate_spatial_target(wall, field) + rule = _traction_rule(wall.edges, wall.radius, wall.sign, quadrature) + np.testing.assert_array_equal(rule["nodes"], wall.nodes) + velocity, pressure = wall.cell_values(), _boundary_pressure_values(stokes, wall) + polygon_stress = _regular_cartesian_stress(wall.vertices, velocity, pressure, rule["points"], viscosity) + circle_stress = _regular_cartesian_stress(wall.vertices, velocity, pressure, rule["circle_points"], viscosity) + jump = np.asarray(uw.function.evaluate(pressure_jump, rule["circle_points"].reshape(-1, 2))).reshape( + rule["circle_weight"].shape) + circle_stress -= jump[..., None, None]*np.eye(2) + # Standard CBF freshly reassembles in the current pressure/essential-data gauge. + coordinates, traction = stokes.boundary_flux(wall.name, mass="consistent", remove_mean=False) + # Native assembly first wires constitutive placeholder gradients to this solver. + stress_difference = stokes.stress_deviator - 2*sp.nsimplify(viscosity, rational=True)*stokes.strainrate + if any(sp.simplify(uw.function.unwrap(value, keep_constants=False, return_self=False)) != 0 + for value in stress_difference): + raise ValueError("Traction recovery requires the symmetric Newtonian stress 2*viscosity*strainrate.") + coordinates = np.concatenate(uw.mpi.comm.allgather(np.asarray(coordinates).reshape(-1, 2))) + traction = np.concatenate(uw.mpi.comm.allgather(np.asarray(traction).reshape(-1, 2))) + distance, index = cKDTree(coordinates).query(wall.nodes) + if np.max(distance) > 1e-10*wall.radius: + raise ValueError("The standard CBF trace does not cover this circular wall.") + values, checks = _transfer_circle_traction(rule, traction[index], polygon_stress, circle_stress) + recovered, reconstruction = [], [] + for component in range(2): + nodes, value, component_checks = _recover_circular_trace(wall.edges, values[:, component], + wall.radius, degree, halfwidth_edges, quadrature) + np.testing.assert_array_equal(nodes, wall.nodes) + recovered.append(value) + reconstruction.append(component_checks) + recovered = np.stack(recovered, axis=-1) + # The public shear field is sigma_rtheta, not outward-normal tangential traction. + for field, value in zip(targets, (values[:, 0], wall.sign*values[:, 1], + recovered[:, 0], wall.sign*recovered[:, 1])): + selected = wall.local_boundary + field.array[selected, 0, 0] = value[wall.local_index[selected]] + checks.update(quadrature=quadrature, viscosity=viscosity, + pressure_transfer="R_regular+L_circle(sigma_physical)-L_polygon(sigma_regular)", + physical_circle_pressure_jump_maximum=float(np.max(np.abs(jump))), + reconstruction=_trace_change_checks(rule, values, recovered), + component_reconstruction_checks=reconstruction, + pressure_proxy="minus normal stress is a zero-slip boundary estimate, not the FE pressure field") + return checks diff --git a/tests/test_1089_cylindrical_dirichlet_traction.py b/tests/test_1089_cylindrical_dirichlet_traction.py new file mode 100644 index 000000000..cab966522 --- /dev/null +++ b/tests/test_1089_cylindrical_dirichlet_traction.py @@ -0,0 +1,269 @@ +"""Hard manufactured-stress contracts for full Dirichlet circular CBF recovery. + +No Stokes solve or benchmark reference is required. Constant pressure fixes +the sign/gauge; affine Cartesian stress fixes every transferred vector moment. +""" + +from enum import Enum +from types import SimpleNamespace + +import numpy as np +import pytest +from petsc4py import PETSc +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import spsolve +from scipy.spatial import cKDTree +import sympy as sp +import underworld3 as uw +from underworld3.coordinates import CoordinateSystemType + +pytestmark = pytest.mark.level_1 + + +def _ring(count, radius, graded): + spacing = 1+(.55*np.sin(np.arange(count)*.61) if graded else np.zeros(count)) + theta = 2*np.pi*np.r_[0, np.cumsum(spacing[:-1])]/spacing.sum() + nodes = radius*np.stack((np.cos(theta), np.sin(theta)), axis=-1) + return np.stack((nodes, np.roll(nodes, -1, axis=0)), axis=1) + + +def _stress(points): + x, y = points[..., 0], points[..., 1] + pressure = .3+.2*x-.4*y + result = np.empty(points.shape[:-1]+(2, 2)) + result[..., 0, 0], result[..., 1, 1] = 2*x-pressure, -2*x-pressure + result[..., 0, 1] = result[..., 1, 0] = -y + return result + + +def _polygon_cbf(rule, tensor): + """Independent boundary Galerkin solve for a known Cartesian stress.""" + index, phi = rule["index"], rule["shape"] + size = len(rule["nodes"]) + local_mass = np.einsum("eq,qi,qj->eij", rule["polygon_weight"], phi, phi) + row = np.broadcast_to(index[:, :, None], local_mass.shape).ravel() + col = np.broadcast_to(index[:, None, :], local_mass.shape).ravel() + mass = coo_matrix((local_mass.ravel(), (row, col)), shape=(size, size)).tocsr() + local_load = np.einsum("eq,qi,eqcd,ed->eic", rule["polygon_weight"], phi, tensor, rule["polygon_normal"]) + load = np.zeros((size, 2)) + np.add.at(load, index.ravel(), local_load.reshape(-1, 2)) + return spsolve(mass, load), load + + +@pytest.mark.tier_b +@pytest.mark.parametrize("radius,sign", [(1.22, -1), (2.22, 1)]) +@pytest.mark.parametrize("graded", [False, True]) +def test_constant_pressure_keeps_absolute_gauge_and_zero_shear(radius, sign, graded): + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule, _transfer_circle_traction + + rule = _traction_rule(_ring(29, radius, graded), radius, sign, 32) + stress = np.broadcast_to(-.7*np.eye(2), rule["points"].shape[:-1]+(2, 2)) + cbf, _ = _polygon_cbf(rule, stress) + actual, checks = _transfer_circle_traction(rule, cbf, stress, stress) + np.testing.assert_allclose(actual[:, 0], -.7, rtol=0, atol=3e-13) + np.testing.assert_allclose(actual[:, 1], 0, rtol=0, atol=3e-13) + assert checks["maximum_cartesian_moment_residual"] < 3e-14 + assert checks["cartesian_resultant_residual"] < 3e-13 + + +@pytest.mark.tier_b +@pytest.mark.parametrize("sign", [-1, 1]) +def test_affine_stress_preserves_every_circle_vector_moment_and_quadrature(sign): + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule, _transfer_circle_traction + + edges, radius = _ring(23, 1.7, True), 1.7 + traces = [] + for quadrature in (24, 48): + rule = _traction_rule(edges, radius, sign, quadrature) + cbf, _ = _polygon_cbf(rule, _stress(rule["points"])) + actual, checks = _transfer_circle_traction(rule, cbf, _stress(rule["points"]), _stress(rule["circle_points"])) + sampled = np.einsum("qi,eik->eqk", rule["shape"], actual[rule["index"]]) + traction = np.einsum("eqk,eqck->eqc", sampled, rule["frame"]) + exact = np.einsum("eqcd,eqd->eqc", _stress(rule["circle_points"]), rule["normal"]) + defect = np.zeros((len(actual), 2)) + local = np.einsum("eq,qi,eqc->eic", rule["circle_weight"], rule["shape"], traction-exact) + np.add.at(defect, rule["index"].ravel(), local.reshape(-1, 2)) + np.testing.assert_allclose(defect, 0, rtol=0, atol=2e-13) + expected_resultant = sign*np.pi*radius**2*np.array([.8, .4]) + np.testing.assert_allclose(checks["cartesian_resultant"], expected_resultant, rtol=3e-13, atol=3e-13) + traces.append(actual) + np.testing.assert_allclose(traces[0], traces[1], rtol=3e-13, atol=3e-13) + + +@pytest.mark.tier_b +def test_regular_pressure_reaction_needs_regular_polygon_stress(): + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule, _transfer_circle_traction + + radius = 2.22 + rule = _traction_rule(_ring(27, radius, True), radius, 1, 32) + regular = _stress(rule["points"]) + cbf, _ = _polygon_cbf(rule, regular) + x = rule["circle_points"] + lift_circle = (radius-np.linalg.norm(x, axis=-1))*.4 + physical_circle = _stress(x)-lift_circle[..., None, None]*np.eye(2) + actual, _ = _transfer_circle_traction(rule, cbf, regular, physical_circle) + reference, _ = _transfer_circle_traction(rule, cbf, regular, _stress(x)) + np.testing.assert_allclose(actual, reference, rtol=0, atol=2e-13) + # The prescribed lift vanishes at the true wall, but not at chord points. + assert np.max(radius-np.linalg.norm(rule["points"], axis=-1)) > 1e-3 + + +@pytest.mark.tier_b +def test_constant_pressure_uses_actual_facet_normals_with_tolerated_vertex_roundoff(): + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule, _transfer_circle_traction + + radius = 1.7 + edges = _ring(31, radius, True) + vertices = edges[:, 0]*(1+4e-11*np.sin(np.arange(len(edges))))[:, None] + edges = np.stack((vertices, np.roll(vertices, -1, axis=0)), axis=1) + rule = _traction_rule(edges, radius, 1, 32) + stress = np.broadcast_to(-.7*np.eye(2), rule["points"].shape[:-1]+(2, 2)) + cbf, _ = _polygon_cbf(rule, stress) + actual, _ = _transfer_circle_traction(rule, cbf, stress, stress) + np.testing.assert_allclose(actual[:, 0], -.7, rtol=0, atol=3e-13) + np.testing.assert_allclose(actual[:, 1], 0, rtol=0, atol=3e-13) + + +@pytest.mark.tier_b +def test_reconstruction_preserves_mean_and_combined_cartesian_resultants(): + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule, _trace_change_checks + from underworld3.utilities._cylindrical_reconstruction import _recover_circular_trace + + edges, radius = _ring(39, 1.7, True), 1.7 + rule = _traction_rule(edges, radius, -1, 48) + theta = np.arctan2(rule["nodes"][:, 1], rule["nodes"][:, 0]) + source = np.stack((.4+np.cos(3*theta)+.03*np.cos(15*theta), + -.2+np.sin(2*theta)+.02*np.sin(17*theta)), axis=-1) + target = np.stack([_recover_circular_trace(edges, source[:, k], radius, 6, 5, 16)[1] + for k in range(2)], axis=-1) + checks = _trace_change_checks(rule, source, target) + assert checks["maximum_cartesian_resultant_change"] < 2e-12 + np.testing.assert_allclose(checks["scalar_integral_changes"], 0, rtol=0, atol=2e-12) + assert checks["maximum_changed_cartesian_nodal_moment"] > 1e-5 + + +@pytest.mark.tier_b +@pytest.mark.parametrize("pressure_profile", ["zero", "affine"]) +def test_real_mesh_native_reassembly_preserves_gauge_and_separate_targets(pressure_profile): + """Exercise labelled facets, essential DOFs and shared MPI nodes without a solve.""" + mesh = uw.meshing.Annulus(radiusInner=.5, radiusOuter=1., cellSize=.18, qdegree=2) + velocity = uw.discretisation.MeshVariable("TractionWitnessU", mesh, 2, degree=2) + pressure = uw.discretisation.MeshVariable("TractionWitnessP", mesh, 1, degree=1, continuous=True) + velocity.array[:, 0, :] = 0 + pressure.array[:, 0, 0] = .7 + stokes = uw.systems.Stokes(mesh, velocityField=velocity, pressureField=pressure) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.viscosity = 1 + stokes.bodyforce = sp.Matrix([0, 0]) + wall = stokes.add_cylindrical_zero_slip_transfer(boundary="Upper", radius=1., sign=1, quadrature=16) + stokes.add_dirichlet_bc(sp.Matrix([[0, 0]]), "Lower") + targets = [uw.discretisation.MeshVariable("TractionWitness"+name, mesh, 1, degree=2) + for name in ("N", "T", "Nr", "Tr")] + for target in targets: + target.array[:, 0, 0] = 777 + jump = (sp.Rational(3, 10)+mesh.N.x/5-2*mesh.N.y/5 + if pressure_profile == "affine" else sp.S.Zero) + checks = stokes.recover_cylindrical_zero_slip_traction(wall=wall, + normal_field=targets[0], shear_field=targets[1], + recovered_normal_field=targets[2], recovered_shear_field=targets[3], + viscosity=1., pressure_jump=jump, degree=6, halfwidth_edges=3, quadrature=24) + for target, exact in zip(targets, [-.7, 0, -.7, 0]): + if pressure_profile == "zero": + np.testing.assert_allclose(target.array[wall.local_boundary, 0, 0], exact, rtol=0, atol=2e-12) + np.testing.assert_array_equal(target.array[~wall.local_boundary, 0, 0], 777) + np.testing.assert_array_equal(pressure.array[:, 0, 0], .7) + np.testing.assert_array_equal(velocity.array[:, 0, :], 0) + assert checks["reconstruction"]["maximum_cartesian_resultant_change"] < 2e-12 + if pressure_profile == "affine": + from underworld3.utilities._cylindrical_dirichlet_traction import _traction_rule + + rule = _traction_rule(wall.edges, 1., 1, 48) + # Gather only trace values; no interior field participates in this audit. + coords = np.concatenate(uw.mpi.comm.allgather(targets[0].coords[wall.local_boundary])) + values = np.concatenate(uw.mpi.comm.allgather(np.stack([ + field.array[wall.local_boundary, 0, 0] for field in targets[:2]], axis=-1))) + distance, index = cKDTree(coords).query(rule["nodes"]) + np.testing.assert_allclose(distance, 0, rtol=0, atol=1e-12) + sample = np.einsum("qi,eik->eqk", rule["shape"], values[index][rule["index"]]) + traction = np.einsum("eqk,eqck->eqc", sample, rule["frame"]) + points = rule["circle_points"] + exact_pressure = 1+points[..., 0]/5-2*points[..., 1]/5 + error = traction+exact_pressure[..., None]*rule["normal"] + defect = np.zeros((len(rule["nodes"]), 2)) + np.add.at(defect, rule["index"].ravel(), np.einsum( + "eq,qi,eqc->eic", rule["circle_weight"], rule["shape"], error).reshape(-1, 2)) + np.testing.assert_allclose(defect, 0, rtol=0, atol=2e-12) + np.testing.assert_allclose(checks["cartesian_resultant"], -np.pi*np.array([.2, -.4]), + rtol=2e-12, atol=2e-12) + + +@pytest.mark.tier_b +def test_pressure_adapter_gathers_boundary_cells_only_with_empty_local_owner(monkeypatch): + from underworld3.utilities import _cylindrical_dirichlet_traction as module + + class BoundaryOnlyComm: + def allgather(self, local): + np.testing.assert_array_equal(local.shape, (0, 3)) + return [local, np.array([[.2, .3, .4]])] + + mesh = SimpleNamespace(_cell_node_indices=lambda degree, continuous: np.empty((0, 3), dtype=int)) + stokes = SimpleNamespace(mesh=mesh, p=SimpleNamespace(array=np.zeros((500, 1, 1)))) + wall = SimpleNamespace(cells=np.empty(0, dtype=int)) + monkeypatch.setattr(uw.mpi, "comm", BoundaryOnlyComm()) + np.testing.assert_array_equal(module._boundary_pressure_values(stokes, wall), [[.2, .3, .4]]) + + +@pytest.mark.tier_b +@pytest.mark.mpi(min_size=4) +@pytest.mark.skipif(uw.mpi.size < 4, reason="Three labelled facets require >=4 ranks to force an empty owner.") +def test_real_three_facet_wall_recovery_with_ranks_owning_no_boundary_cells(): + """A complete real DMPlex wall has fewer facets than MPI ranks by construction.""" + angle = 2*np.pi*np.arange(3)/3 + coordinates = np.concatenate([radius*np.stack((np.cos(angle), np.sin(angle)), axis=-1) + for radius in (.5, .75, 1.)]) + cells = [] + for layer in range(2): + for i in range(3): + a, b = 3*layer+i, 3*layer+(i+1)%3 + c, d = 3*(layer+1)+i, 3*(layer+1)+(i+1)%3 + cells.extend(((a, c, d), (a, d, b))) + dm = PETSc.DMPlex().createFromCellList(2, np.asarray(cells, dtype=PETSc.IntType), coordinates) + + class Boundaries(Enum): + Lower = 1 + Upper = 2 + + for name in ("Lower", "Upper"): + dm.createLabel(name) + f0, f1 = dm.getDepthStratum(1) + v0, v1 = dm.getDepthStratum(0) + section = dm.getCoordinateSection() + coords = dm.getCoordinatesLocal().array.reshape(-1, 2) + for facet in range(f0, f1): + vertices = [point for point in dm.getTransitiveClosure(facet)[0] if v0 <= point < v1] + points = coords[[section.getOffset(point)//2 for point in vertices]] + for name, radius, tag in (("Lower", .5, 1), ("Upper", 1., 2)): + if np.allclose(np.linalg.norm(points, axis=-1), radius, rtol=0, atol=1e-12): + dm.setLabelValue(name, facet, tag) + mesh = uw.discretisation.Mesh(dm, degree=1, qdegree=2, simplex=True, + boundaries=Boundaries, markVertices=True, coordinate_system_type=CoordinateSystemType.CYLINDRICAL2D) + velocity = uw.discretisation.MeshVariable("EmptyTractionU", mesh, 2, degree=2) + pressure = uw.discretisation.MeshVariable("EmptyTractionP", mesh, 1, degree=1) + velocity.array[:, 0, :] = 0 + pressure.array[:, 0, 0] = .7 + stokes = uw.systems.Stokes(mesh, velocityField=velocity, pressureField=pressure) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.viscosity = 1 + stokes.bodyforce = sp.Matrix([0, 0]) + wall = stokes.add_cylindrical_zero_slip_transfer(boundary="Lower", radius=.5, sign=-1, quadrature=48) + stokes.add_dirichlet_bc(sp.Matrix([[0, 0]]), "Upper") + counts = uw.mpi.comm.allgather(len(wall.cells)) + assert sum(counts) == 3 and min(counts) == 0 + fields = [uw.discretisation.MeshVariable("EmptyTraction"+str(k), mesh, 1, degree=2) for k in range(4)] + checks = stokes.recover_cylindrical_zero_slip_traction(wall=wall, + normal_field=fields[0], shear_field=fields[1], recovered_normal_field=fields[2], + recovered_shear_field=fields[3], viscosity=1., pressure_jump=sp.S.Zero, + degree=6, halfwidth_edges=.5, quadrature=48) + for field, value in zip(fields, [-.7, 0, -.7, 0]): + np.testing.assert_allclose(field.array[wall.local_boundary, 0, 0], value, rtol=0, atol=2e-12) + assert checks["cartesian_resultant_residual"] < 2e-12