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) grad_field = 'temperature'
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[0]
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:] 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] 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
temperature_gradient_x = grad_func(0, 'x') temperature_gradient_y = grad_func(1, 'y') temperature_gradient_z = grad_func(2, 'z')
f = temperature_gradient_x[sl_right,1:-1,1:-1]/ds f -= temperature_gradient_x[sl_left ,1:-1,1:-1]/ds if data.ds.dimensionality > 1: ds = div_fac * data['dy'].flat[0] f += temperature_gradient_y[1:-1,sl_right,1:-1]/ds f -= temperature_gradient_y[1:-1,sl_left ,1:-1]/ds if data.ds.dimensionality > 2: ds = div_fac * data['dz'].flat[0] f += temperature_gradient_z[1:-1,1:-1,sl_right]/ds 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:
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] 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 funconce 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 mailto: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 mailto: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 mailto: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 http://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 mailto:yt-users@python.org To unsubscribe send an email to yt-users-leave@python.org mailto:yt-users-leave@python.org https://mail.python.org/mailman3/lists/yt-users.python.org/ https://mail.python.org/mailman3/lists/yt-users.python.org/ Member address: matthewturk@gmail.com mailto:matthewturk@gmail.com _______________________________________________ yt-users mailing list -- yt-users@python.org mailto:yt-users@python.org To unsubscribe send an email to yt-users-leave@python.org mailto:yt-users-leave@python.org https://mail.python.org/mailman3/lists/yt-users.python.org/ https://mail.python.org/mailman3/lists/yt-users.python.org/ Member address: robertmjenningsjr@berkeley.edu mailto:robertmjenningsjr@berkeley.edu
yt-users mailing list -- yt-users@python.org mailto:yt-users@python.org To unsubscribe send an email to yt-users-leave@python.org mailto:yt-users-leave@python.org https://mail.python.org/mailman3/lists/yt-users.python.org/ https://mail.python.org/mailman3/lists/yt-users.python.org/ Member address: matthewturk@gmail.com mailto:matthewturk@gmail.com _______________________________________________ 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: robertmjenningsjr@berkeley.edu