Inspecting gradient of a field at boundary
Hi people, I am trying to define curl of a vector field in yt, but first I needed to see if correct derivatives are taken at the boundary. I am starting with a mock data set and followed yt tutorial: import yt from yt import YTArray import numpy as np posx = np.linspace(-1.0,1.0,12) posy = np.linspace(-1.0,1.0,12) posz = np.linspace(-1.0,1.0,12) X, Y, Z = np.meshgrid(posx, posy, posz, indexing='ij') # From gradient output, yt seems to use ij indexing psi = X data_dict = {"psi":(psi,'kg/cm**3'), "number_of_particles":12, "particle_position_x":(posx, 'kg/cm**3'), "particle_position_y":(posy, 'kg/cm**3'), "particle_position_z":(posz, 'kg/cm**3')} bbox = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) y = yt.load_uniform_grid(data_dict,domain_dimensions=(12,12,12),bbox=bbox) grad_fields = y.add_gradient_fields(("stream","psi")) ad = y.all_data() print( ad["psi_gradient_x"] ) print( ad["psi_gradient_y"] ) print( ad["psi_gradient_z"] ) Output: yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: current_time = 0.0 yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: domain_dimensions = [12 12 12] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_left_edge = [-1. -1. -1.] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_right_edge = [ 1. 1. 1.] yt : [INFO ] 2017-05-17 19:40:16,412 Parameters: cosmological_simulation = 0.0 [-5.45454545 -5.45454545 -5.45454545 ..., -5.45454545 -5.45454545 -5.45454545] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4 Now my question is, since psi(x,y,z) = x, the gradient should be (1,0,0)...and it is everywhere except the boundaries of the box. Can you please explain to me what is going wrong here? I don't think YT is interpreting the bbox and grid points as I am. Best
Hi Tazkera,
I think what's happening is that the gradient is getting periodically
wrapped. It helps if you take a look at a SlicePlot of your 'psi' and
'psi_gradient_x' fields:
psi: http://i.imgur.com/EK73NNO.png
psi_gradient_x: http://i.imgur.com/wXPpqi9.png
The gradient at the left and right edges of the domain is being calculated
using the periodically wrapped values (i.e. the local difference between
neighboring values of the psi field shoots up from 1 in the middle of the
domain to 12 at the edges). This causes the gradient to have different
values at the edges of the domain.
I guess you are looking for a gradient field that uses forward differences
on the left side of the domain and backward differences on the right side
of the domain rather than being periodically wrapped? Unfortunately right
now yt doesn't have an option to specify that. That said, I think it would
be possible to do e.g. by adding a new keyword argument to
Datset.add_gradient_field and then adding the new difference logic to the
setup_gradient_fields function in yt/fields/fluid_fields.py.
Hope that helps,
Nathan
On Thu, May 18, 2017 at 10:15 AM, tazkera haque
Hi people,
I am trying to define curl of a vector field in yt, but first I needed to see if correct derivatives are taken at the boundary. I am starting with a mock data set and followed yt tutorial:
import yt from yt import YTArray import numpy as np
posx = np.linspace(-1.0,1.0,12) posy = np.linspace(-1.0,1.0,12) posz = np.linspace(-1.0,1.0,12)
X, Y, Z = np.meshgrid(posx, posy, posz, indexing='ij') # From gradient output, yt seems to use ij indexing
psi = X
data_dict = {"psi":(psi,'kg/cm**3'), "number_of_particles":12, "particle_position_x":(posx, 'kg/cm**3'), "particle_position_y":(posy, 'kg/cm**3'), "particle_position_z":(posz, 'kg/cm**3')}
bbox = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]])
y = yt.load_uniform_grid(data_dict,domain_dimensions=(12,12,12), bbox=bbox)
grad_fields = y.add_gradient_fields(("stream","psi"))
ad = y.all_data()
print( ad["psi_gradient_x"] ) print( ad["psi_gradient_y"] ) print( ad["psi_gradient_z"] )
Output: yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: current_time = 0.0 yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: domain_dimensions = [12 12 12] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_left_edge = [-1. -1. -1.] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_right_edge = [ 1. 1. 1.] yt : [INFO ] 2017-05-17 19:40:16,412 Parameters: cosmological_simulation = 0.0 [-5.45454545 -5.45454545 -5.45454545 ..., -5.45454545 -5.45454545 -5.45454545] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4
Now my question is, since psi(x,y,z) = x, the gradient should be (1,0,0)...and it is everywhere except the boundaries of the box. Can you please explain to me what is going wrong here? I don't think YT is interpreting the bbox and grid points as I am.
Best
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi Nathan,
Thanks for the explanation. I will follow up if I make any progress.
Best
Tazkera
On Thu, May 18, 2017 at 11:55 AM, Nathan Goldbaum
Hi Tazkera,
I think what's happening is that the gradient is getting periodically wrapped. It helps if you take a look at a SlicePlot of your 'psi' and 'psi_gradient_x' fields:
psi: http://i.imgur.com/EK73NNO.png psi_gradient_x: http://i.imgur.com/wXPpqi9.png
The gradient at the left and right edges of the domain is being calculated using the periodically wrapped values (i.e. the local difference between neighboring values of the psi field shoots up from 1 in the middle of the domain to 12 at the edges). This causes the gradient to have different values at the edges of the domain.
I guess you are looking for a gradient field that uses forward differences on the left side of the domain and backward differences on the right side of the domain rather than being periodically wrapped? Unfortunately right now yt doesn't have an option to specify that. That said, I think it would be possible to do e.g. by adding a new keyword argument to Datset.add_gradient_field and then adding the new difference logic to the setup_gradient_fields function in yt/fields/fluid_fields.py.
Hope that helps,
Nathan
On Thu, May 18, 2017 at 10:15 AM, tazkera haque
wrote: Hi people,
I am trying to define curl of a vector field in yt, but first I needed to see if correct derivatives are taken at the boundary. I am starting with a mock data set and followed yt tutorial:
import yt from yt import YTArray import numpy as np
posx = np.linspace(-1.0,1.0,12) posy = np.linspace(-1.0,1.0,12) posz = np.linspace(-1.0,1.0,12)
X, Y, Z = np.meshgrid(posx, posy, posz, indexing='ij') # From gradient output, yt seems to use ij indexing
psi = X
data_dict = {"psi":(psi,'kg/cm**3'), "number_of_particles":12, "particle_position_x":(posx, 'kg/cm**3'), "particle_position_y":(posy, 'kg/cm**3'), "particle_position_z":(posz, 'kg/cm**3')}
bbox = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]])
y = yt.load_uniform_grid(data_dict,domain_dimensions=(12,12,12), bbox=bbox)
grad_fields = y.add_gradient_fields(("stream","psi"))
ad = y.all_data()
print( ad["psi_gradient_x"] ) print( ad["psi_gradient_y"] ) print( ad["psi_gradient_z"] )
Output: yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: current_time = 0.0 yt : [INFO ] 2017-05-17 19:40:16,410 Parameters: domain_dimensions = [12 12 12] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_left_edge = [-1. -1. -1.] yt : [INFO ] 2017-05-17 19:40:16,411 Parameters: domain_right_edge = [ 1. 1. 1.] yt : [INFO ] 2017-05-17 19:40:16,412 Parameters: cosmological_simulation = 0.0 [-5.45454545 -5.45454545 -5.45454545 ..., -5.45454545 -5.45454545 -5.45454545] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4 [ 0. 0. 0. ..., 0. 0. 0.] kg/cm**4
Now my question is, since psi(x,y,z) = x, the gradient should be (1,0,0)...and it is everywhere except the boundaries of the box. Can you please explain to me what is going wrong here? I don't think YT is interpreting the bbox and grid points as I am.
Best
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (2)
-
Nathan Goldbaum
-
tazkera haque