Hi Michael, sorry about the length reply time on this.

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

ds = yt.load('../pleiades/grf/grf_pr_2_3_nm5_3/TI.uov.00200.athdf',\
             units_override=units_override,unit_system="cgs")

grad_fields = ds.add_gradient_fields(("gas","temperature"))

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]
    f  = data["temperature_gradient_x"][sl_right,1:-1,1:-1]/ds
    f -= data["temperature_gradient_x"][sl_left ,1:-1,1:-1]/ds
    if data.ds.dimensionality > 1:
        ds = div_fac * data['dy'].flat[0]
        f += data["temperature_gradient_y"][1:-1,sl_right,1:-1]/ds
        f -= data["temperature_gradient_y"][1:-1,sl_left ,1:-1]/ds
    if data.ds.dimensionality > 2:
        ds = div_fac * data['dz'].flat[0]
        f += data["temperature_gradient_z"][1:-1,1:-1,sl_right]/ds
        f -= data["temperature_gradient_z"][1:-1,1:-1,sl_left ]/ds
    new_field = np.zeros(data["temperature_gradient_x"].shape, dtype='float64')
    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']

ds.add_field(('gas','therm_cond'), function=_therm_cond, units="erg/(g*s)", sampling_type="cell", take_log=False,
             display_name=r'$\nabla \cdot (\kappa\nabla T)/\rho$', force_override=True, validators = [ValidateSpatial(1,
                           ['temperature_gradient_x','temperature_gradient_y','temperature_gradient_z'])])

Has anyone come across this issue before? Is this a known bug or am I just implementing the divergence field incorrectly? Any help is greatly appreciated.

Thanks,
Michael
_______________________________________________
yt-users mailing list -- yt-users@python.org
To unsubscribe send an email to yt-users-leave@python.org
https://mail.python.org/mailman3/lists/yt-users.python.org/
Member address: matthewturk@gmail.com