uw.function.evaluate returns negative values for a squared DERIVATIVE expression

[lmoresi]

Summary

uw.function.evaluate returns negative values for a squared derivative, and negative values for sqrt(...) of one. A square cannot be negative, so any invariant built the obvious way — strain-rate second invariant, stress invariant, |grad T| — is silently wrong.

Expressions that do not contain a derivative are exact, so this is specific to the derivative path.

Reproducer

Plain box mesh, no cuts, no adaptation, ~20 lines:

import numpy as np, sympy, underworld3 as uw

mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(0., 0.), maxCoords=(1., 1.), cellSize=1/8,
    regular=False, qdegree=2)
x, y = mesh.X
f = uw.discretisation.MeshVariable("f", mesh, 1, degree=2)
c = np.asarray(f.coords)
f.data[:, 0] = np.sin(6.0 * c[:, 0]) * np.cos(6.0 * c[:, 1])

pts = np.random.default_rng(0).uniform(0.05, 0.95, size=(4000, 2))

d  = np.asarray(uw.function.evaluate(f.sym[0].diff(x), pts)).ravel()
d2 = np.asarray(uw.function.evaluate(f.sym[0].diff(x)**2, pts)).ravel()
s  = np.asarray(uw.function.evaluate(sympy.sqrt(f.sym[0].diff(x)**2), pts)).ravel()

print("evaluate(df/dx)          min", d.min(),  "max", d.max())
print("evaluate(df/dx)**2       min", (d**2).min())          # correct
print("evaluate((df/dx)**2)     min", d2.min(), "negatives", (d2 < 0).sum())
print("evaluate(sqrt((df/dx)**2)) min", s.min(), "negatives", (s < 0).sum())
print("max |evaluate(d**2) - evaluate(d)**2| =", np.abs(d2 - d**2).max())

Observed

evaluate(df/dx)            min -5.6261   max +6.0884
evaluate(df/dx)**2         min +0.0000                 <- correct
evaluate((df/dx)**2)       min -2.6386   negatives 614 of 4000
evaluate(sqrt((df/dx)**2)) min -0.0622   negatives 25
max |evaluate(d**2) - evaluate(d)**2| = 3.3521

sqrt((df/dx)**2) is |df/dx| and cannot be negative.

Expected

evaluate(expr(D)) should agree with composing the evaluated derivative pointwise, or at least never violate the sign of a square.

Contrast: no derivative, no problem

The same test on the field itself is exact:

evaluate(f**2) vs evaluate(f)**2 : max difference 0.00000, 0 negatives

Likely cause

evaluate appears to project internally. For a non-derivative expression the projection is exact enough to be invisible; for a derivative of a P2 field (a P1-per-cell, discontinuous quantity) the square is formed on a projected/continuous representation and undershoots.

Impact

Found while computing a strain-rate invariant next to a sharp viscosity contrast. There the sqrt argument reached -143.7, and the derived field was unusable near the interface — while looking plausible in the bulk (medians were within 0.4% of correct, so it is easy to miss).

Workaround

Evaluate the derivative components separately and compose the invariant in numpy:

exx = np.asarray(uw.function.evaluate(E[0, 0], pts)).ravel()
eyy = np.asarray(uw.function.evaluate(E[1, 1], pts)).ravel()
exy = np.asarray(uw.function.evaluate(E[0, 1], pts)).ravel()
inv = np.sqrt(0.5 * (exx**2 + eyy**2 + 2 * exy**2))

For fields next to a material discontinuity, better still to L2-project the stress components to a continuous field first and form invariants from the recovered components.

Underworld development team with AI support from Claude Code