From 35a63e84cfc7739ae34efb3b1a0f395b78690d10 Mon Sep 17 00:00:00 2001
From: Saurabh Shukla <136088049+ss2098@users.noreply.github.com>
Date: Wed, 29 Jul 2026 09:55:55 -0800
Subject: [PATCH] Update anisotropic convection demo for current APIs

---
 scripts/aniso_convection_demo.py | 30 +++++++++++++++++++++++-------
 1 file changed, 23 insertions(+), 7 deletions(-)

diff --git a/scripts/aniso_convection_demo.py b/scripts/aniso_convection_demo.py
index 474b07684..8cd1f1b7f 100644
--- a/scripts/aniso_convection_demo.py
+++ b/scripts/aniso_convection_demo.py
@@ -31,6 +31,23 @@
 import matplotlib.pyplot as plt
 import matplotlib.tri as mtri
 import underworld3 as uw
+
+
+DT_SAFETY = 0.1
+
+
+def scalar_dt(value):
+    """Convert estimate_dt() output to a plain scalar timestep."""
+    try:
+        dt = float(value)
+    except TypeError:
+        arr = np.asarray(value, dtype=float)
+        dt = float(np.nanmin(arr))
+
+    if not np.isfinite(dt) or dt <= 0.0:
+        raise ValueError(f"Bad timestep from estimate_dt(): {value!r}")
+
+    return DT_SAFETY * dt
 from underworld3.meshing import smooth_mesh_interior
 from underworld3.meshing.smoothing import _tri_cells, _signed_areas
 
@@ -68,8 +85,7 @@ def run_convection(mesh, r, th, v, P, t_soln, cellsize):
     stokes.penalty = 0.0
     unit_r = mesh.CoordinateSystem.unit_e_0
     stokes.add_essential_bc((0.0, 0.0), mesh.boundaries.Lower.name)
-    stokes.add_natural_bc(1.0e6 * v.sym.dot(unit_r) * unit_r,
-                          mesh.boundaries.Upper.name)
+    stokes.add_essential_bc((0.0, 0.0), mesh.boundaries.Upper.name)
     T_cond = (r_o - r) / (r_o - r_inner)
     stokes.bodyforce = RA * (t_soln.sym[0] - T_cond) * unit_r
 
@@ -91,9 +107,9 @@ def run_convection(mesh, r, th, v, P, t_soln, cellsize):
 
     t_sim = 0.0
     for s in range(N_STEPS):
-        dt = adv_diff.estimate_dt()
+        dt = scalar_dt(adv_diff.estimate_dt())
         adv_diff.solve(timestep=dt, zero_init_guess=False)
-        stokes.solve(zero_init_guess=False)
+        stokes.solve(zero_init_guess=True)
         t_sim += dt
         tt = t_soln.data[:, 0]
         print(f"  step {s+1:2d}: t={t_sim:.4f} Δt={dt:.2e} "
@@ -112,6 +128,7 @@ def run_convection(mesh, r, th, v, P, t_soln, cellsize):
     print(f"=== Ra={RA:.0e} annulus convection, {N_STEPS} steps, "
           f"res-{RES} (FIXED mesh) ===")
     run_convection(mesh, r, th, v, P, t_soln, cellsize)
+    os.makedirs(os.path.dirname(CACHE), exist_ok=True)
     np.savez(CACHE, T=np.asarray(t_soln.data),
              V=np.asarray(v.data),
              Xc=np.asarray(mesh.X.coords))
@@ -150,9 +167,8 @@ def run_convection(mesh, r, th, v, P, t_soln, cellsize):
 
 # --- 3. refine the mesh on the T-gradient metric -------------------
 A0 = np.abs(_signed_areas(X_orig, tris))
-print(f"=== refine: method='anisotropic' on ρ∝|∇T| ===")
-smooth_mesh_interior(mesh, metric=metric, method="anisotropic",
-                     verbose=True)
+print("=== refine: node_redistribution on ρ∝|∇T| ===")
+uw.meshing.node_redistribution(mesh, metric)
 X_ref = np.asarray(mesh.X.coords).copy()
 A1 = np.abs(_signed_areas(X_ref, tris))
 print(f"minA/meanA  before={A0.min()/A0.mean():.4f}  "