Hi Matthew,

Thank you for your help on this. Modifying my above code with your suggestion I have:

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)
slice_3d = (slice(1, -1), slice(1, -1), slice(1, -1))
div_fac = 2.0
ds = div_fac * data['dx'].flat

slice_3dl = slice_3d[:axi] + (sl_left,) + slice_3d[axi+1:]
slice_3dr = slice_3d[:axi] + (sl_right,) + slice_3d[axi+1:]
ds = div_fac * data[ftype, "d%s" % ax]
new_field = np.zeros_like(data[grad_field], dtype=np.float64)
new_field = data.ds.arr(new_field, f.units)
new_field[slice_3d] = f
return new_field

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

So far this reproduces the original unwanted results. I believe the next step would be in regards to the ghost zones you mentioned. Unfortunately I don’t have any experience in this. I do know that the first argument for ValidateSpatial refers to the number of ghost zones and changing to accept 2 ghost zones results in an error. Any additional help would be greatly appreciated.

Best,
Michael

On Aug 6, 2020, at 11:15 AM, Matthew Turk <matthewturk@gmail.com> wrote:

Hi Michael,

Sure -- basically, I am suggesting wrapping the gradient field right into your calculation.  Basically, the gradient field is defined something like:

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]
new_field = np.zeros_like(data[grad_field], dtype=np.float64)
new_field = data.ds.arr(new_field, f.units)
new_field[slice_3d] = f
return new_field
return func

once for each axis (axi here).  I am suggesting, instead of requesting the field "temperature_gradient_x", compute the gradient yourself using code like this.  Please do let me know if this helps, or if it's unclear or not working.

On Thu, Aug 6, 2020 at 12:29 PM Michael Jennings <robertmjenningsjr@berkeley.edu> wrote:
Hi Matthew,

Thank you for your insight! I understand your explanation and suggestion for the most part but I am unsure about how to go about this. Is there anything in the documentation you could point me to that may help me? I’m somewhat new to this.

Thanks,
Michael

On Aug 5, 2020, at 11:56 AM, Matthew Turk <matthewturk@gmail.com> wrote:

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

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
f -= data["temperature_gradient_x"][sl_left ,1:-1,1:-1]/ds
if data.ds.dimensionality > 1:
ds = div_fac * data['dy'].flat
f -= data["temperature_gradient_y"][1:-1,sl_left ,1:-1]/ds
if data.ds.dimensionality > 2:
ds = div_fac * data['dz'].flat
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,

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/
_______________________________________________
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/