DDt semi-Lagrangian: unit mismatch on geographic meshes with units active

[lmoresi]

Summary

The SemiLagrangian_DDt.update_pre_solve() in systems/ddt.py fails with unit incompatibility errors when used on geographic meshes with the units system active. The AdvDiffusion solver is unusable on geographic meshes.

Error

ValueError: Cannot subtract arrays with incompatible units: 'meter' and 'meter / second'

at ddt.py:1450:

mid_pt_coords = coords - v_at_node_pts * (0.5 * dt_for_calc)

Root Cause

Three quantities with inconsistent units:

• coords = psi_star.coords → UnitAwareArray in meters (dimensional)
• v_at_node_pts → dimensionalised to m/s (correct)
• dt_for_calc → bare float or incorrectly handled (should be seconds)

When dt is passed as a bare float (nondimensional), line 1348-1349 treats it as dimensionless rather than converting to seconds:

# If dt is already a dimensionless number, treat it as seconds
dt_for_calc = dt  # BUG: doesn't actually convert to seconds

When dt is passed as a uw.quantity in seconds, the JIT compilation appears to hang (>30 min with no output) — likely due to the complex symbolic expressions from geographic basis vectors interacting with the dimensional dt.

Reproduction

import underworld3 as uw

uw.reset_default_model()
model = uw.Model()
model.set_reference_quantities(
    length=uw.quantity(1000, "km"),
    viscosity=uw.quantity(1e21, "Pa.s"),
    diffusivity=uw.quantity(1e-6, "m**2/s"),
)

mesh = uw.meshing.RegionalGeographicBox(
    lon_range=(135.5, 137.5), lat_range=(-34.5, -33.0),
    depth_range=(uw.quantity(0, "km"), uw.quantity(50, "km")),
    ellipsoid="WGS84", numElements=(8, 8, 4), simplex=True,
)

v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=1,
                                    continuous=True, units="m/s")
C = uw.discretisation.MeshVariable("C", mesh, 1, degree=2)

adv_diff = uw.systems.AdvDiffusion(mesh, u_Field=C, V_fn=v)
adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel
adv_diff.constitutive_model.Parameters.diffusivity = 1e-4
adv_diff.add_essential_bc(1.0, "Bottom")

adv_diff.solve(timestep=1e-4)  # Fails with unit mismatch

Suggested Fix

The DDt code should either:

1. Work in nondimensional DM space internally — strip units from coords, nondimensionalise velocity and dt, do all arithmetic in nondimensional, then re-wrap
2. Consistently dimensionalise everything — ensure coords (m), velocity (m/s), and dt (s) are all in compatible SI units before arithmetic

Option 1 is simpler and matches how the PETSc solver backend works.

Additional Notes

• Velocity variables must have units= set (e.g. units="m/s") for DDt to detect the dimensional mode. Without units, DDt gets V_fn units: dimensionless and the mismatch is different (dimensionless vs meters).
• This also surfaces a broader issue: MeshVariable units should be required when the units system is active, not optional. Currently a missing units= causes silent failures in downstream solvers.

Context

Discovered while building the H2Ex fault-controlled flow workflow on the feature/fault-system-workflow branch. The Darcy velocity field needs to be advected through an AdvDiffusion solver for concentration field accumulation mapping.

[lmoresi]

Additional Finding: DMInterpolation performance on adapted meshes

The unit mismatch issue was partially a misdiagnosis — the AdvDiffusion solver does work on geographic meshes when:

• Concentration field has no units (dimensionless)
• Velocity has units="m/s"
• Timestep is a uw.quantity in seconds

However, on adapted meshes (~6000 P1 / ~44k P2 nodes), the DDt semi-Lagrangian evaluation hangs in DMInterpolation point location:

ddt.py:1452 → uw.function.global_evaluate()
  → functions_unit_system.py:530 → _global_evaluate_nd()
    → dminterpolation_cache.py:88 → store_structure()

This takes >3 minutes (alarm-killed) for a single advection step on a 6058-node adapted tet mesh. The same solve completes instantly on a small (67 node) mesh.

The issue is DMInterpolation point location for ~44k P2 query points on an adapted (non-uniform) tet mesh. The RBF evaluation path (mode="fast") handles this efficiently (~0.03s per 484-point batch), but global_evaluate in DDt uses the heavier DMInterpolation path.

Possible fixes:

1. Use RBF path in DDt for the semi-Lagrangian upstream evaluation
2. Use P1 for the transported variable (fewer query points)
3. Cache the DMInterpolation point location between timesteps (points don't move much)

Reproduction on adapted geographic mesh:

# After loading adapted mesh from checkpoint + darcy velocity
C = uw.discretisation.MeshVariable("C", mesh, 1, degree=2)  # P2 = ~44k DOFs
adv = uw.systems.AdvDiffusion(mesh, u_Field=C, V_fn=v_darcy)
adv.constitutive_model = uw.constitutive_models.DiffusionModel
adv.constitutive_model.Parameters.diffusivity = 1e-5
adv.add_essential_bc(1.0, "Bottom")
adv.add_essential_bc(0.0, "Surface")
adv.solve(timestep=uw.quantity(1e14, "s"))  # Hangs in DMInterpolation

Works instantly on a 4×4×2 mesh (67 nodes / ~500 P2 DOFs).

[lmoresi]

Correction: units work correctly, DMInterpolation is the real bottleneck

The unit mismatch was caused by missing units= on the velocity variable, not a DDt bug. With proper setup (velocity units="m/s", dimensionless C, timestep from estimate_dt()), the AdvDiffusion solve works on small geographic meshes.

The actual problem is DMInterpolation performance on adapted tet meshes:

• Small mesh (67 P1 nodes): instant
• Adapted mesh (6058 P1 nodes): hangs >90s per step in dminterpolation_cache.store_structure

This occurs even with tiny timesteps (dt_nd = 8e-8, well within CFL). The global_evaluate → _global_evaluate_nd → DMInterpolation point location is O(slow) on irregular adapted meshes.

The RBF evaluation path (mode="fast") handles the same mesh in ~0.03s per batch. The DDt semi-Lagrangian should use this path for upstream velocity evaluation.