My guess, which might be addressable, is that this is related to the way that the temperature gradient and the *additional* step of another gradient field interact.  I *think* what might be happening is that your gradient fields are being compute with 1 ghost zones, and then whe nyou get the additional ghost zone for the thermal condition, it does something strange (since it would need a net of two ghost zones, which I don't know that it will do.)

It is probable in my mind that you would be able to fix this by doing a new field that takes two ghost zones and then does both the temperature gradient and then applies your calculation.  I know that would be a bit more work, but it may fix the issue.

On Sun, Jul 19, 2020 at 9:31 PM Michael Jennings <robertmjenningsjr@berkeley.edu> wrote:
Hi all,

I am attempting to add a divergence field and am running into the issue apparent in the plot below:

The grid artifacts appear to be the meshblocks for mpi parallelization when running the simulation. This doesn’t occur for any other field.
I am working with an Athena++ athdf dataset and using yt version 3.7.dev0.
Here is the script I used to add the desired field:

import yt
from yt.fields.derived_field import ValidateSpatial
import yt.units as u

kappa0 = 9.68e4

units_override=units_override,unit_system="cgs")

def _therm_cond(field,data):
kappa = kappa0*u.erg/(u.cm*u.K*u.s)
sl_left = slice(None, -2, None)
sl_right = slice(2, None, None)
div_fac = 2.0
ds = div_fac * data['dx'].flat[0]
if data.ds.dimensionality > 1:
ds = div_fac * data['dy'].flat[0]
if data.ds.dimensionality > 2:
ds = div_fac * data['dz'].flat[0]
new_field = data.ds.arr(new_field,'K/cm**2')
new_field[1:-1,1:-1,1:-1] = f
return kappa*new_field/data['density']

display_name=r'$\nabla \cdot (\kappa\nabla T)/\rho$', force_override=True, validators = [ValidateSpatial(1,