diff --git a/biosteam/_unit.py b/biosteam/_unit.py
index 6e3363ab..0277375c 100644
--- a/biosteam/_unit.py
+++ b/biosteam/_unit.py
@@ -52,11 +52,17 @@ def phenomena_oriented_run(self):
     Ts = [i.T for i in outs]
     Q = self.Hnet
     Unit.run(self)
-    self._dmol = sum(
+    new = sum(
         [i.mol for i in outs],
         -sum([i.mol for i in ins], 0)
     )
-    self._dmol[np.abs(self._dmol) < 1e-9] = 0.
+    if hasattr(self, '_dmol'):
+        old = self._dmol
+        f = bst.PhasePartition.dmol_relaxation_factor
+        self._dmol = dmol = f * old + (1 - f) * new
+    else:
+        self._dmol = dmol = new
+    self._dmol[np.abs(dmol) < 1e-9] = 0.
     self._duty = self.Hnet
     Ts_new = [i.T for i in outs]
     if all([i == j for i, j in zip(Ts_new, Ts)]): # T constant
diff --git a/biosteam/plots/plots.py b/biosteam/plots/plots.py
index 93d07ad6..1f657f06 100644
--- a/biosteam/plots/plots.py
+++ b/biosteam/plots/plots.py
@@ -1043,7 +1043,7 @@ def generate_contour_data(
         if filterwarnings:
             from warnings import filterwarnings
             filterwarnings('ignore')
-        data0 = z_at_xy(x0, y0, *args)
+        data0 = np.asarray(z_at_xy(x0, y0, *args))
         shape = data0.shape
         if len(shape) == 1:
             shape = f"({shape[0]})"