Fix geographic mesh units: ellipsoid quantities, checkpoint, Darcy sign

[lmoresi]

Summary

Fixes for geographic meshes with units active.

1. Ellipsoid as uw.quantity — geographic.py stores uw.quantity(a, "km") so symbolic depth/basis vector expressions carry units through JIT nondimensionalisation
2. Checkpoint round-trip — ellipsoid quantities serialized to HDF5 metadata and restored on load
3. Darcy velocity sign — projection was +flux, corrected to -flux
4. surfaces.py — extract km magnitude from quantity for from_trace()

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes several unit-handling and restart issues for GEOGRAPHIC meshes by carrying ellipsoid parameters as unit-aware quantities through coordinate expressions and checkpoint metadata, and corrects the Darcy velocity projection sign.

Changes:

• Store ellipsoid semi-axes as uw.quantity(..., "km") when units are active, and nondimensionalise them only where numeric gmsh inputs are required.
• Persist/restore ellipsoid metadata in mesh HDF5 checkpoints and plumb it into CoordinateSystem initialisation on reload.
• Correct Darcy velocity projection to use v = -darcy_flux, and update surface tracing to safely extract km magnitudes from quantity ellipsoids.

Reviewed changes

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

Show a summary per file

File	Description
src/underworld3/systems/solvers.py	Fixes Darcy velocity projection sign (v = -flux).
src/underworld3/meshing/surfaces.py	Ensures from_trace() converts ellipsoid quantities to km floats for geographic conversions.
src/underworld3/meshing/geographic.py	Stores ellipsoid a/b as quantities under units; nondimensionalises for gmsh mesh generation.
src/underworld3/discretisation/discretisation_mesh.py	Adds checkpoint serialization/deserialization for ellipsoid metadata and passes it into coordinate system init.
src/underworld3/coordinates.py	Updates geographic coordinate computations and symbolic expressions to work with quantity ellipsoid parameters and checkpoint restoration.

---

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

[Copilot]

The regions metadata is loaded from the checkpoint and assigned to a local regions variable, but it is never used afterwards (unlike boundaries and coordinate_system_type). Either wire this through (e.g., store it on self and/or use it to rebuild labels) or remove the dead code to avoid confusion about whether region restoration is supported.

Suggested change

	regions = None
	try:
	json_str = f["metadata"].attrs["regions"]
	rgn_dict = json.loads(json_str)
	regions = Enum("Regions", rgn_dict)
	except KeyError:
	pass

[Copilot]

This sign change affects the meaning/direction of the projected Darcy velocity field (self.v) and could easily regress again because the existing Darcy solver tests primarily validate pressure, not velocity direction. Add/extend a test to assert that the projected velocity matches -darcy_flux (e.g., check sign/direction at a few sample points for a known pressure gradient).

[Copilot]

Ellipsoid serialization/deserialization is new checkpoint behavior for geographic meshes (including quantity handling). Please add a checkpoint round-trip test that: creates a GEOGRAPHIC mesh with units active, writes it, reloads it, and asserts the restored ellipsoid values/units and coordinate conversions (geo↔cartesian) are consistent. This will protect the HDF5 metadata schema and the quantity reconstruction logic from future breakage.

Suggested change

	# Save ellipsoid metadata for geographic meshes
	if hasattr(self.CoordinateSystem, "ellipsoid"):
	ellipsoid_ser = {}
	for k, v in self.CoordinateSystem.ellipsoid.items():
	if hasattr(v, "to"): # uw.quantity
	ellipsoid_ser[k] = {
	"value": float(v.magnitude),
	"unit": str(v.units),
	}
	else:
	ellipsoid_ser[k] = v
	def _serialise_checkpoint_metadata_value(value):
	if (
	hasattr(value, "magnitude")
	and hasattr(value, "units")
	and hasattr(value, "to")
	):
	magnitude = value.magnitude
	if isinstance(magnitude, numpy.ndarray):
	magnitude = magnitude.tolist()
	elif isinstance(magnitude, numpy.generic):
	magnitude = magnitude.item()
	return {
	"__uw_quantity__": True,
	"value": magnitude,
	"unit": str(value.units),
	}
	if isinstance(value, dict):
	return {
	key: _serialise_checkpoint_metadata_value(val)
	for key, val in value.items()
	}
	if isinstance(value, (list, tuple)):
	return [
	_serialise_checkpoint_metadata_value(val) for val in value
	]
	if isinstance(value, numpy.ndarray):
	return value.tolist()
	if isinstance(value, numpy.generic):
	return value.item()
	return value
	# Save ellipsoid metadata for geographic meshes
	if hasattr(self.CoordinateSystem, "ellipsoid"):
	ellipsoid_ser = {
	k: _serialise_checkpoint_metadata_value(v)
	for k, v in self.CoordinateSystem.ellipsoid.items()
	}

[Copilot]

GeographicCoordinateAccessor.depth still dimensionalises nondimensional depth using self.cs.ellipsoid.get('L_ref_km', 1000), but RegionalGeographicBox no longer stores L_ref_km (or a_nd/b_nd). This will silently default to 1000 km and return incorrect dimensional depths when units are active. Consider deriving the length scale from the mesh/model (e.g., mesh.length_scale converted to km, or model.get_scale_for_dimensionality('[length]')) instead of relying on an ellipsoid dict field, or reintroduce a reliably-populated reference length entry during mesh creation/checkpoint restore.

Suggested change

	self._nondimensional = True
	self._nondimensional = True
	# Preserve the reference length used for nondimensionalisation so
	# downstream dimensional depth conversion does not silently fall
	# back to an incorrect default.
	if "L_ref_km" not in ellipsoid:
	ref_scales_km = []
	if a != 0.0:
	ref_scales_km.append(float(a_raw.to("kilometer").magnitude) / a)
	if b != 0.0:
	ref_scales_km.append(float(b_raw.to("kilometer").magnitude) / b)
	if ref_scales_km:
	ellipsoid["L_ref_km"] = sum(ref_scales_km) / len(ref_scales_km)