Fix thermal convection units tutorial

[ss2098]

This updates docs/examples/Tutorial_Thermal_Convection_Units.py so that it runs with the current unit-aware array and strict-units behavior.

The previous tutorial failed in several places under the current API:

1. The buoyancy expression mixed a MeshVariable symbolic field with a unit-bearing reference temperature, producing a dimensional compatibility error.
2. The initial temperature assignment attempted to assign unit-aware evaluate() output directly into MeshVariable storage, which required explicit conversion to a plain numeric array.
3. The semi-Lagrangian advection step encountered incompatible coordinate and timestep units unless the velocity field carried physical velocity units.
4. The final reporting attempted to convert model-unit temperature quantities directly to offset Celsius units, which raised a Pint OffsetUnitCalculusError.
5. The initial temperature profile could become incorrectly scaled because unit-aware coordinates needed to be normalized before constructing the model-unit temperature field.

The updated tutorial keeps the physical inputs in units, but uses explicit model-unit numeric values in solver expressions. It also converts unit-aware arrays before assignment, constructs the initial temperature field from normalized coordinates, and reports original physical parameters for interpretation.

No solver behavior is changed.

Tested with:

pixi run python docs/examples/Tutorial_Thermal_Convection_Units.py

Result: script completes successfully.

[lmoresi]

Holding this one — it exposes a real units + semi-Lagrangian bug

Thanks @ss2098 — the setup half of the tutorial is correctly updated (it gets all the way through the Rayleigh-number / scaling section). But I ran it end-to-end on current development (note: CI does not execute docs/examples/, so green checks don't cover this), and it fails at the time-stepping loop:

Running 10 coupled Stokes + thermal steps...
  thermal.solve(timestep=dt, ...)
  -> DuDt.update_pre_solve(...)            # ddt.py:2651
  ValueError: Cannot subtract arrays with incompatible units: 'meter' and 'meter / second'

Root cause (core, not tutorial)

The SL trace-back midpoint is coords - v * (0.5*dt) (ddt.py:2651). In the units-active path, coords is UnitAwareArray[meter] and v is [meter/second], but dt arrives as a plain float and is not dimensionalised (ddt.py:2519):

if has_units:
    if hasattr(dt, "to"):
        dt_for_calc = dt.to("second")
    else:
        dt_for_calc = dt          # comment says "treat as seconds" but leaves it unitless

so v[m/s] * dt[dimensionless] = [m/s], and [m] - [m/s] raises. Simply attaching seconds here would fix the units but be numerically wrong — the tutorial's CFL dt is in nondimensional model-time, not SI seconds, so it needs the model time-scale (the not has_units branch just below does the consistent thing: nondimensionalise coords, velocity and dt together).

So the tutorial can't run until the units-active SL trace-back consistently reconciles dt with the dimensional coords/velocity. This is a genuine UW3 bug the tutorial usefully surfaced, not something to patch in the example. I'll file it separately and link it here; let's hold this PR until that lands (then this tutorial should be the natural regression test for it).

Everything else in the PR (the setup/units fixes) looks good and can ride along once the core fix is in.

[lmoresi]

Filed the underlying bug: #267 — holding this PR until that lands.

[ss2098]

Thanks Professor @lmoresi , this makes sense. I will leave PR #263 open and on hold until #267 is fixed.

I agree that this should be handled in the semi-Lagrangian units path rather than patched only in the tutorial. Once the core fix lands, I can rebase the tutorial branch on current development, rerun Tutorial_Thermal_Convection_Units.py end-to-end, and update the PR if needed.

Thanks for identifying and filing the underlying issue.

[lmoresi]

Unblocked: the underlying SL trace-back bug is fixed on development (#277 / #267). The time-loop crash this tutorial hit is resolved. Please rebase this PR on current development and re-run the tutorial end-to-end — the time-stepping section should now complete. Once it runs clean, this is good to merge (and it'll serve as the natural integration test for #277).

[ss2098]

Thanks Professor @lmoresi. I rebased this branch on current development after #277 landed and re-ran the tutorial end-to-end.

Tested with:

pixi run python docs/examples/Tutorial_Thermal_Convection_Units.py

The time-stepping section now completes successfully.

Result:

✅ Time-stepping completed!
🎉 SUCCESS: Physical thermal convection tutorial completed!

I have pushed the rebased branch.

[lmoresi]

Re-ran end-to-end on current development (which now has the #267 trace-back fix): good progress — the time-loop crash is gone, the tutorial gets all the way to the buoyancy force. It now hits a different, distinct units bug at buoyancy = ... * (temperature - reference_temp) (line ~202): EnhancedMeshVariable - UWQuantity fails even for compatible units. Filed as #282.

So this PR is unblocked on #267 but blocked on #282. Two options:

1. wait for EnhancedMeshVariable ± UWQuantity fails for compatible units (blocks units-active buoyancy T − T_ref) #282 (the proper fix — meshvar ± quantity should just work), or
2. work around it in the tutorial now by dropping to .sym + a model-unit reference value (e.g. temperature.sym[0,0] - reference_temp_model), which I verified works.

Happy to push option (2) so the tutorial merges now, or hold for #282 — your call / @ss2098's.

[ss2098]

Thanks Professor @lmoresi. I think option (2) makes sense here so the tutorial can merge now.

Since the tutorial is already using physical quantities for setup and then explicit model-unit numeric values in solver expressions, using temperature.sym[0,0] - reference_temp_model for the buoyancy term seems consistent with the current tutorial design. The deeper EnhancedMeshVariable - UWQuantity issue can then remain tracked separately in #282 as the proper core fix.

Please feel free to push that workaround, or I can update the branch on my side if you prefer.

[lmoresi]

Units-active convection now runs — the system blockers are fixed

Re-ran on current development. The simulation reaches ✅ Time-stepping completed! (10 coupled Stokes + thermal steps). The three underlying system bugs this tutorial exposed are all fixed and merged:

1. Semi-Lagrangian trace-back is unit-inconsistent under active units (dt not dimensionalised in has_units path) #267 — semi-Lagrangian trace-back under units (fix(ddt): units-active semi-Lagrangian trace-back (#267) #277)
2. EnhancedMeshVariable ± UWQuantity fails for compatible units (blocks units-active buoyancy T − T_ref) #282 — meshvar ± UWQuantity arithmetic (fix(units): meshvar arithmetic with a UWQuantity operand (#282) #283) — the T − T_ref buoyancy
3. bodyforce setter rejecting a UWQuantity component (fix(stokes): bodyforce setter accepts a UWQuantity component #284) — stokes.bodyforce = [0, buoyancy]

So the units-active convection path works now. What remains is tutorial-side:

• Line 388 (and 78/79): T_model.to('degC') raises OffsetUnitCalculusError — Pint can't convert a model-scaled temperature to the offset unit °C. Use .to('K') for display, or compute °C manually as K − 273.15. (These are diagnostic prints, after the sim.)
• Stability: mean_T grows over the 10 steps (−6 → −8873), i.e. the convection is numerically unstable — looks like the timestep (dt ≈ 4.5e13 model-time) / Rayleigh scaling needs tuning so the run stays bounded. Worth checking the CFL dt and the buoyancy magnitude against the intended Ra.

@ss2098 — once the °C prints use K (or manual −273.15) and the timestep/scaling is tuned so the run stays stable, this should be ready. The hard part (the three system bugs) is done. Happy to help with the stability tuning if useful.

[ss2098]

Professor @lmoresi, I also tested docs/examples/porous_flow/advanced/Ex_viscousFingering.py on current development and found that the migrated script currently fails in the swarm path. First swarm.populate(fill_param=4) hits a NameError for all_local_cells, and after bypassing populate the later swarm.advection path hits additional missing internal bookkeeping such as particle_X_orig and cellid.

I made a local runnable version by removing the swarm-material route and using the existing AdvDiffusionSLCN mesh-field material advection instead. That version runs end-to-end and finishes with:

Viscous fingering example complete.

However, this changes the example from swarm-tracked material advection to mesh-field advection-diffusion, so I do not want to open it as a simple bug-fix PR without checking. Would you prefer that I open this as a script-safe migration PR, or should the swarm/populate/advection issues be treated as a separate core bug first?

[ss2098]

Thanks Professor @lmoresi. I updated the tutorial-side issues after rebasing on current development.

Changes made:

• Replaced direct degC conversions with Kelvin display plus manual Celsius conversion, avoiding Pint OffsetUnitCalculusError for offset units.
• Kept the thermal timestep unit-aware so the semi-Lagrangian trace-back receives a proper time quantity.
• Added a conservative timestep cap for the tutorial run so the 10-step demonstration remains bounded.

Tested with:

pixi run python docs/examples/Tutorial_Thermal_Convection_Units.py

Result:

✅ Time-stepping completed!
Temperature field range: 0.200 to 1.242 model units
Mean temperature: 1.047 model units
🎉 SUCCESS: Physical thermal convection tutorial completed!

I pushed the updated branch.