Prevent NaN arrays in gradients

yt/fields/fluid_fields.py

@@ -208,11 +208,19 @@ def grad_func(axi, ax):
         slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:]
         slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:]
         def func(field, data):
-            ds = div_fac * data[ftype, "d%s" % ax]
-            f  = data[grad_field][slice_3dr]/ds[slice_3d]
-            f -= data[grad_field][slice_3dl]/ds[slice_3d]
+            ds = data[ftype, "d%s" % ax]
+            vl = data[grad_field][slice_3dl]
+            vc = data[grad_field][slice_3d]
+            vr = data[grad_field][slice_3dr]
+
+            okl = np.isfinite(vl)
+            okr = np.isfinite(vr)
+
+            f = np.where(okl & okr, (vr - vl) / 2,
+                np.where(okl, vc - vl, vr - vc)) / ds[slice_3d]
+
             new_field = np.zeros_like(data[grad_field], dtype=np.float64)
-            new_field = data.ds.arr(new_field, f.units)
+            new_field = data.ds.arr(new_field, vr.units / ds.units)
             new_field[slice_3d] = f
             return new_field
         return func

yt/fields/vector_operations.py

@@ -138,6 +138,7 @@ def create_vector_fields(registry, basename, field_units,
     if slice_info is None:
         sl_left = slice(None, -2, None)
         sl_right = slice(2, None, None)
+        sl_centre = slice(1, -1, None)
         div_fac = 2.0
     else:
         sl_left, sl_right, div_fac = slice_info