diff --git a/docs/advanced/cylindrical-geometry-correction.md b/docs/advanced/cylindrical-geometry-correction.md new file mode 100644 index 000000000..78a2c148a --- /dev/null +++ b/docs/advanced/cylindrical-geometry-correction.md @@ -0,0 +1,336 @@ +# 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. 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. + +## 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 +# 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 +``` + +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. + +## 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. + +### 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 + +$$ +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. + +## 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 +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, 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. + +## 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. + 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. + +- 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/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..930734892 100644 --- a/docs/developer/guides/boundary-condition-rulings.md +++ b/docs/developer/guides/boundary-condition-rulings.md @@ -87,3 +87,31 @@ 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. + +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/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..68f1489be 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8652,6 +8652,427 @@ 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 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 + ----- + 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. + + 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 + ----- + 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. + + 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 + ----- + 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 + 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 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_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. + + 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. + + 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 + 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 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. + + 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 + ----- + 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 + 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_dirichlet.py b/src/underworld3/utilities/_cylindrical_dirichlet.py new file mode 100644 index 000000000..3c14e0d38 --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_dirichlet.py @@ -0,0 +1,356 @@ +"""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) + # 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) + 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/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/src/underworld3/utilities/_cylindrical_pressure.py b/src/underworld3/utilities/_cylindrical_pressure.py new file mode 100644 index 000000000..9c5b31435 --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_pressure.py @@ -0,0 +1,244 @@ +"""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. + +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 + +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..e18f3b8eb --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_reactions.py @@ -0,0 +1,183 @@ +"""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 + +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 +import underworld3 as uw +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 + + +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 _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. + + 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() + 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 + 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_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/src/underworld3/utilities/_cylindrical_wall.py b/src/underworld3/utilities/_cylindrical_wall.py new file mode 100644 index 000000000..734ca8a37 --- /dev/null +++ b/src/underworld3/utilities/_cylindrical_wall.py @@ -0,0 +1,429 @@ +"""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. + +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 +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..49aaf9a2e --- /dev/null +++ b/tests/parallel/test_1084_rotated_datum_collectives.py @@ -0,0 +1,45 @@ +"""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 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"] + # 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) + + +@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) + 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_dirichlet.py b/tests/test_1089_cylindrical_dirichlet.py new file mode 100644 index 000000000..6de82bb04 --- /dev/null +++ b/tests/test_1089_cylindrical_dirichlet.py @@ -0,0 +1,373 @@ +"""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) + + +@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. 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 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 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