From 6b8f2256ebc24eca0947b95563f5962d82cd0d0f Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Tue, 3 Mar 2026 14:54:15 +1100 Subject: [PATCH] Fix PETSc error in mesh.adapt() when extra MeshVariables exist When variables beyond the solver fields (e.g. strain rate, metric) existed on the mesh, adapt() would trigger "Invalid field number N; not in [0, M)" errors from PETSc. Two problems: 1. Old variable lvecs were destroyed one-at-a-time inside the reinit loop. When variable N's _setup_ds() ran its internal backup/restore loop over all mesh._vars, variables N+1, N+2 etc still held lvecs with stale field_ids from the pre-adaptation DM, causing createSubDM to fail. 2. The metric field was excluded from reinitialization and deleted after adaptation. As a user-created variable with external references, it should survive adaptation like any other variable. Fix: destroy all old lvecs upfront before reinitializing any variable, include the metric in the reinitialization set, and stop deleting it. Fixes #48 Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 48 ++++++++----------- 1 file changed, 20 insertions(+), 28 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 177c3a290..9c290392d 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2961,13 +2961,12 @@ def adapt(self, metric_field, verbose=False): print(f"[{uw.mpi.rank}] Notifying surface '{surface.name}' (marking distance stale)...", flush=True) surface._on_mesh_adapted(self) - # Capture current variable data, excluding only the metric field - # (which becomes invalid after adaptation) - # All other variables (including surface distance fields) are reinitialized + # Capture all user-supplied variables for reinitialization on the new mesh. + # The metric field is included — it's a user-created variable that may + # have external references and be reused in subsequent adaptation cycles. old_vars_data = {} - metric_name = metric_field.name if hasattr(metric_field, 'name') else None for var_name, var in self._vars.items(): - if var is not None and var_name != metric_name: + if var is not None: old_vars_data[var_name] = var # Stack boundary labels for adaptation @@ -3054,27 +3053,26 @@ def mesh_update_callback(array, change_context): # Rebuild coordinate navigation self.nuke_coords_and_rebuild(verbose=False) + # Destroy ALL old vectors upfront before reinitializing any variable. + # This is critical because _setup_ds() iterates mesh._vars to backup/restore + # data — if some variables still hold lvecs with stale field_ids from the + # pre-adaptation DM, createSubDM will fail on the new DM. (Fixes #48) + for old_var in old_vars_data.values(): + if old_var._lvec is not None: + old_var._lvec.destroy() + old_var._lvec = None + if old_var._gvec is not None: + old_var._gvec.destroy() + old_var._gvec = None + if hasattr(old_var, '_canonical_data'): + old_var._canonical_data = None + if hasattr(old_var, '_cached_data_array'): + old_var._cached_data_array = None + # Reinitialize MeshVariables on the new mesh # Note: Variables are reset to zero. Users should reinitialize with data. for var_name, old_var in old_vars_data.items(): try: - # Destroy old vectors - if old_var._lvec is not None: - old_var._lvec.destroy() - old_var._lvec = None - if old_var._gvec is not None: - old_var._gvec.destroy() - old_var._gvec = None - - # Eagerly invalidate cached data arrays. The .data property also - # self-validates via _lvec identity check, but clearing here avoids - # unnecessary recreation on next access. - if hasattr(old_var, '_canonical_data'): - old_var._canonical_data = None - if hasattr(old_var, '_cached_data_array'): - old_var._cached_data_array = None - - # Re-setup the variable on the new mesh old_var._setup_ds() old_var._set_vec(available=True) @@ -3094,12 +3092,6 @@ def mesh_update_callback(array, change_context): if verbose: print(f"[{uw.mpi.rank}] Solver marked for rebuild", flush=True) - # Remove only the metric field from mesh._vars - # (it was specific to the pre-adaptation mesh and is now invalid) - # Surface distance variables stay - they're just marked stale and will recompute - if metric_name and metric_name in self._vars: - del self._vars[metric_name] - # Clear caches self._evaluation_hash = None self._evaluation_interpolated_results = None