Skip to content

projection: unit-aware smoothing_length API on the four Projection classes#200

Open
lmoresi wants to merge 1 commit into
developmentfrom
feature/projection-smoothing-length
Open

projection: unit-aware smoothing_length API on the four Projection classes#200
lmoresi wants to merge 1 commit into
developmentfrom
feature/projection-smoothing-length

Conversation

@lmoresi
Copy link
Copy Markdown
Member

@lmoresi lmoresi commented May 21, 2026

@benknight — this is the bit you asked me about: a length-scale-aware way to set the smoother on Projection (and friends) instead of the unit-less smoothing knob. Splitting it out of a larger meshing branch so it can land independently.

Summary

Adds an L-valued, unit-aware smoothing_length accessor to the four projection classes:

  • SNES_Projection
  • SNES_Vector_Projection
  • SNES_Tensor_Projection
  • SNES_MultiComponent_Projection

The existing smoothing knob is preserved unchanged (units of length², historically used as a tiny numerical regulariser); smoothing_length is a parallel view with the physical convention smoothing = smoothing_length². The two stay consistent: setting one updates both.

Why — the mathematics

The L2-projection solvers all minimise a Tikhonov-regularised misfit

$$\min_u; \tfrac12!\int_\Omega (u-\tilde f)^2 ;+; \tfrac{\alpha}{2}!\int_\Omega |\nabla u|^2 ,$$

whose Euler–Lagrange equation is the screened-Poisson equation

$$u ;-; \nabla!\cdot!\bigl(\alpha,\nabla u\bigr) ;=; \tilde f .$$

The free-space Green's function of $\bigl(1 - \alpha\nabla^2\bigr)$ decays as $\exp(-r/L)$ where

$$L ;=; \sqrt{\alpha},$$

so the projection acts as a Gaussian-like convolution of width $L$, obtained implicitly by one elliptic solve — no kernel is ever assembled. In Fourier space the transfer function is

$$\widehat{u}(k) ;=; \frac{1}{1 + k^{2} L^{2}},\widehat{\tilde f}(k),$$

so features at scale $\ll L$ are damped, features at $\gg L$ pass through, and features at $\sim L$ are attenuated by ≈ 1/2.

This makes $L$ — not $\alpha$ — the natural physical knob: you ask for "smooth this field at scale $L$" and the solver delivers it.

API

proj = uw.systems.Projection(mesh, V)

# L-valued, unit-aware
proj.smoothing_length = 0.05                        # plain ND scalar
proj.smoothing_length = 5 * uw.scaling.units.km     # Pint Quantity → ND'd

print(proj.smoothing)         # α = L²  (sympy)
print(proj.smoothing_length)  # Pint Quantity in m when a scaling
                              # context is set; otherwise √α

Identical surface on Vector_Projection, Tensor_Projection, and MultiComponent_Projection. The class docstrings now carry the screened-Poisson derivation (with guidance: L ≳ h is needed for any real filtering; L ≈ 1–2·h is a useful light de-noising default).

Test plan

tests/test_0505_projection_smoothing_length.py is new and locks:

  • L → α = L² round-trip on all four classes
  • α → smoothing_length = √α round-trip
  • Unit-aware getter returns a Pint Quantity in length units
  • End-to-end behaviour: smoothing a step function at successively larger L monotonically shrinks the projected field's peak-to-peak range (the screened-Poisson low-pass behaviour the docstring promises)
  • Existing tests/test_0504_projections.py and tests/test_multicomponent_projection.py pass unchanged
tests/test_0505_projection_smoothing_length.py ......       [100%]
6 passed in 6.4s
13 passed in 40.8s   (existing projection tests)

Underworld development team with AI support from Claude Code

…asses

The L2-projection solvers solve

    u  -  ∇·(α ∇u)  =  ũ

i.e. a screened-Poisson smoother whose Green's function decays as
exp(-r/L) with L = √α — equivalent in effect to a Gaussian convolution
of width L, obtained by one elliptic solve without ever forming the
kernel.

`smoothing` (= α, units length²) was the only existing knob and is
historically used as a tiny *numerical* regulariser (e.g. 1e-6), which
gives sub-grid L and no physical smoothing. This change adds a
parallel L-valued, unit-aware `smoothing_length` accessor on:

  - SNES_Projection
  - SNES_Vector_Projection
  - SNES_Tensor_Projection (inherits)
  - SNES_MultiComponent_Projection

Setter accepts a plain float, a Pint Quantity with length units, or
any unit-aware object understood by uw.non_dimensionalise; the value
is non-dimensionalised, squared and stored in `self._smoothing`, so
the two views stay consistent. Getter returns a Pint Quantity in
length units when a scaling context is set, else the plain ND float.

Class docstrings have been rewritten to put `smoothing_length` in
its mathematical context (screened-Poisson form, Green's function,
attenuation as 1/(1+k²L²), guidance on choosing L relative to the
cell size). Vector / Tensor / MultiComponent classes defer the full
derivation to SNES_Projection.

A new tests/test_0505_projection_smoothing_length.py locks:

  - L → α = L² round-trip on all four classes
  - α → smoothing_length = √α round-trip
  - End-to-end: smoothing a step function at increasing L shrinks
    the peak-to-peak range monotonically (screened-Poisson low-pass
    behaviour)

Existing projection tests pass unchanged.

Underworld development team with AI support from Claude Code
Copilot AI review requested due to automatic review settings May 21, 2026 05:21
Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR adds a new, physically meaningful unit-aware smoothing_length accessor across the Projection solver family, intended as a length-scale view of the existing smoothing (α) coefficient with the convention α = smoothing_length². It also expands class/property docstrings to document the screened-Poisson interpretation and introduces a new test module to lock the API.

Changes:

  • Adds smoothing_length property (getter/setter) to SNES_Projection, SNES_Vector_Projection, and SNES_MultiComponent_Projection (and therefore to SNES_Tensor_Projection via inheritance), keeping smoothing and smoothing_length consistent.
  • Updates projection class/property docstrings with a screened-Poisson / Helmholtz derivation and guidance on choosing the filter scale.
  • Adds a new test module to check round-trips and end-to-end monotone smoothing behavior.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 5 comments.

File Description
src/underworld3/systems/solvers.py Adds smoothing_length API and substantial docstring updates for projection solvers.
tests/test_0505_projection_smoothing_length.py New tests for round-trip behavior and a basic end-to-end smoothing monotonicity check.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +2216 to +2224
except Exception:
# Fall back to magnitude-or-float coercion if the
# value doesn't carry/expect units.
if hasattr(L, "magnitude"):
L_nd = L.magnitude
else:
L_nd = L
self._smoothing = sympify(L_nd) ** 2

except (TypeError, ValueError):
return sympy.sqrt(s)
if sval < 0:
return None
Comment on lines +27 to +38
def test_scalar_smoothing_length_roundtrip(mesh):
"""L=0.05 ⇒ α=L²=0.0025; reading back gives L (unit-aware)."""
V = uw.discretisation.MeshVariable(
"V_sl", mesh, vtype=uw.VarType.SCALAR, degree=2, continuous=True)
proj = sys.Projection(mesh, V)
proj.smoothing_length = 0.05
assert float(proj.smoothing) == pytest.approx(0.0025)
L = proj.smoothing_length
# In an active scaling context this is a Pint Quantity in metres;
# without one it's a plain float — both should round-trip.
L_num = getattr(L, "magnitude", L)
assert float(L_num) == pytest.approx(0.05)
Comment on lines +2432 to +2433
self._smoothing = sympify(L_nd) ** 2

Comment on lines +2768 to +2776
try:
L_nd = uw.non_dimensionalise(L)
except Exception:
if hasattr(L, "magnitude"):
L_nd = L.magnitude
else:
L_nd = L
self._smoothing = sympify(L_nd) ** 2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants