Re: [yt-users] yt-users Digest, Vol 57, Issue 7
Thanks, Andrew! It works now, I altered my yt-x86_64/src/yt-hg/yt/frontends/orion/fields.py to match yours. Changing add_orion_field to add_field corrects the problem. I'm not quite sure why though. On Nov 13, 2012, at 5:04 AM, yt-users-request@lists.spacepope.org wrote:
Send yt-users mailing list submissions to yt-users@lists.spacepope.org
To subscribe or unsubscribe via the World Wide Web, visit http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org or, via email, send a message with subject or body 'help' to yt-users-request@lists.spacepope.org
You can reach the person managing the list at yt-users-owner@lists.spacepope.org
When replying, please edit your Subject line so it is more specific than "Re: Contents of yt-users digest..."
Today's Topics:
1. Orion Derived Fields (Anna Rosen) 2. Re: Orion Derived Fields (Andrew Myers) 3. Re: Strange Behaviour in ProjectionPlot (Patrick Rieser)
----------------------------------------------------------------------
Message: 1 Date: Mon, 12 Nov 2012 16:04:00 -0800 From: Anna Rosen <rosen@ucolick.org> To: yt-users@lists.spacepope.org Subject: [yt-users] Orion Derived Fields Message-ID: <64BC5460-361B-441A-A676-CE433376AF3B@ucolick.org> Content-Type: text/plain; charset=us-ascii
Hi all,
I am unable to access the derived fields written for Orion such as temperature, pressure, etc. The definitions are these fields are located here: yt-x86_64/src/yt-hg/yt/frontends/orion/fields.py
This is the error that I get if I try to access temperature:
In [14]: mi,ma= dd.quantities["Extrema"]("Temperature")[0] --------------------------------------------------------------------------- KeyError Traceback (most recent call last) <ipython-input-14-c8df716340ab> in <module>() ----> 1 mi,ma= dd.quantities["Extrema"]("Temperature")[0]
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in __call__(self, *args, **kwargs) 90 self._data_source.pf.h.io) 91 if lazy_reader and not self.force_unlazy: ---> 92 return self._call_func_lazy(args, kwargs) 93 else: 94 return self._call_func_unlazy(args, kwargs)
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in _call_func_lazy(self, args, kwargs) 97 self.retvals = [ [] for i in range(self.n_ret)] 98 for gi,g in enumerate(self._get_grids()): ---> 99 rv = self.func(GridChildMaskWrapper(g, self._data_source), *args, **kwargs) 100 if not iterable(rv): rv = (rv,) 101 for i in range(self.n_ret): self.retvals[i].append(rv[i])
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in _Extrema(data, fields, non_zero, filter) 586 mins, maxs = [], [] 587 for field in fields: --> 588 if data[field].size < 1: 589 mins.append(1e90) 590 maxs.append(-1e90)
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in __getitem__(self, item) 58 def __getitem__(self, item): 59 if item not in self.local_cache: ---> 60 data = self.data_source._get_data_from_grid(self.grid, item) 61 self.local_cache[item] = data 62 return self.local_cache[item]
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in save_state(self, grid, field, *args, **kwargs) 89 old_keys = grid.field_data.keys() 90 grid.field_parameters = self.field_parameters ---> 91 tr = func(self, grid, field, *args, **kwargs) 92 grid.field_parameters = old_params 93 grid.field_data = YTFieldData( [(k, grid.field_data[k]) for k in old_keys] )
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _get_data_from_grid(self, grid, field) 2555 return na.array([f[i,:][pointI] for i in range(3)]) 2556 else: -> 2557 tr = grid[field] 2558 if tr.size == 1: # dx, dy, dz, cellvolume 2559 tr = tr * na.ones(grid.ActiveDimensions, dtype='float64')
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in __getitem__(self, key) 155 """ 156 if not self.field_data.has_key(key): --> 157 self.get_data(key) 158 return self.field_data[key] 159
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in get_data(self, field, convert) 198 else: raise 199 else: --> 200 self._generate_field(field) 201 return self.field_data[field] 202
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in _generate_field(self, field) 145 self[field] = self.pf.field_info[field](self) 146 else: # Can't find the field, try as it might --> 147 raise exceptions.KeyError(field) 148 149 def has_key(self, key):
KeyError: 'Temperature'
Thanks!
Anna Rosen Graduate Student UCSC, Astronomy & Astrophysics rosen@ucolick.org
------------------------------
Message: 2 Date: Mon, 12 Nov 2012 18:30:14 -0800 From: Andrew Myers <atmyers@berkeley.edu> To: Discussion of the yt analysis package <yt-users@lists.spacepope.org> Subject: Re: [yt-users] Orion Derived Fields Message-ID: <CAB7F9=PN4DL9NfHnj0gBaUBdBkewVGhx7+_g1o5ewwn7Mm+TXg@mail.gmail.com> Content-Type: text/plain; charset="iso-8859-1"
Hi Anna,
Check out this<https://bitbucket.org/atmyers/yt/changeset/53fa42c3dbc9c0228203900a1b1a169e732e9a8c>changeset, which I think fixes the problem. I'm not 100% sure that I understand the distinction between "fallback" and "known" fields, but I can confirm that with those changes all the basic orion fields work.
-Andrew
On Mon, Nov 12, 2012 at 4:04 PM, Anna Rosen <rosen@ucolick.org> wrote:
Hi all,
I am unable to access the derived fields written for Orion such as temperature, pressure, etc. The definitions are these fields are located here: yt-x86_64/src/yt-hg/yt/frontends/orion/fields.py
This is the error that I get if I try to access temperature:
In [14]: mi,ma= dd.quantities["Extrema"]("Temperature")[0] --------------------------------------------------------------------------- KeyError Traceback (most recent call last) <ipython-input-14-c8df716340ab> in <module>() ----> 1 mi,ma= dd.quantities["Extrema"]("Temperature")[0]
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in __call__(self, *args, **kwargs) 90 self._data_source.pf.h.io) 91 if lazy_reader and not self.force_unlazy: ---> 92 return self._call_func_lazy(args, kwargs) 93 else: 94 return self._call_func_unlazy(args, kwargs)
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in _call_func_lazy(self, args, kwargs) 97 self.retvals = [ [] for i in range(self.n_ret)] 98 for gi,g in enumerate(self._get_grids()): ---> 99 rv = self.func(GridChildMaskWrapper(g, self._data_source), *args, **kwargs) 100 if not iterable(rv): rv = (rv,) 101 for i in range(self.n_ret): self.retvals[i].append(rv[i])
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in _Extrema(data, fields, non_zero, filter) 586 mins, maxs = [], [] 587 for field in fields: --> 588 if data[field].size < 1: 589 mins.append(1e90) 590 maxs.append(-1e90)
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.pyc in __getitem__(self, item) 58 def __getitem__(self, item): 59 if item not in self.local_cache: ---> 60 data = self.data_source._get_data_from_grid(self.grid, item) 61 self.local_cache[item] = data 62 return self.local_cache[item]
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in save_state(self, grid, field, *args, **kwargs) 89 old_keys = grid.field_data.keys() 90 grid.field_parameters = self.field_parameters ---> 91 tr = func(self, grid, field, *args, **kwargs) 92 grid.field_parameters = old_params 93 grid.field_data = YTFieldData( [(k, grid.field_data[k]) for k in old_keys] )
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.pyc in _get_data_from_grid(self, grid, field) 2555 return na.array([f[i,:][pointI] for i in range(3)]) 2556 else: -> 2557 tr = grid[field] 2558 if tr.size == 1: # dx, dy, dz, cellvolume 2559 tr = tr * na.ones(grid.ActiveDimensions, dtype='float64')
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in __getitem__(self, key) 155 """ 156 if not self.field_data.has_key(key): --> 157 self.get_data(key) 158 return self.field_data[key] 159
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in get_data(self, field, convert) 198 else: raise 199 else: --> 200 self._generate_field(field) 201 return self.field_data[field] 202
/home1/arosen2/yt-x86_64/src/yt-hg/yt/data_objects/grid_patch.pyc in _generate_field(self, field) 145 self[field] = self.pf.field_info[field](self) 146 else: # Can't find the field, try as it might --> 147 raise exceptions.KeyError(field) 148 149 def has_key(self, key):
KeyError: 'Temperature'
Thanks!
Anna Rosen Graduate Student UCSC, Astronomy & Astrophysics rosen@ucolick.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
Hi, everybody and Matthew: I happened to encounter the blog of Matthew in the internet, and left some notes in the blog. I was led to this mailing list with the following request: I like Python so much that I would like to do my experimentation with radiative transfer with Python environment. So, yt is a perfect environment. But it does not support a full radiative transfer equation yet. By a full radiative transfer equation I mean: (1) dot (w, del L (x, w) ) = LeV_0 (x,w) + sigma_s (x) [ int_{ Sp} f_p( w, x, w' ) dw' ] - simga_t(x) L(x, w) for all x in a given volume V where V_0 is the "free " part of the volume V in which objects are not occupied. Here del is the spatial gradient operator, L(x,w) is the radiance at point x in the direction of w. LeV_0 (x,w) is the emitting radiance at point x in the direction of w. sigma_s(x) is the scattering coefficient at x. Sp is the set of all solid angles. f_p(w, x, w') is the phase function at x, from direction w' to direction w. sigma_t (x) is the extinction or attenuation coefficient, and sigma_t = sigma_s + sigma_a, where sigma_a is the absorption coefficient. (2) boundary conditions: for all x on surfaces of objects within the volume V: L(x,w) = LeS (x, w) + int_{Sp} f_s (w, x, w') L(x, w') dot ( N(x), w' ) dw' Here LeS(x,w) is the emitting radiace at x on a surface to direction w. f_s(w, x, w') is the reflection function at x from direction w' to direction w. By combining (1) and (2), we get a single integral equation: (3) L(x, w) = rau( x_S, w) L (x_S, w) + int_{x_S, x} tau(x', x) [ LeV_0 (x', w) + sigma_s(x') int_{Sp } f_p( w, x', w') dw' ] dx' Here tau(x', x) is the optical depth from x' to x: tau(x', x) = exp [ -int_{x, x'} sigma_t (x'') dx'' ] It is not easy to find the radiance function L(x,w) that satisfies the integral equation of (3), where L occurs in both sides of the equation. In computer graphics [within computer science], this problem has been solved mainly by means of Monte Carlo Integration technique, in particular "path tracing" techniques. In these techniques, we randomly generate a set of paths from the camera to the light sources, and compute the solution L by averaging the radiances carried by each path. A path is a sequence of positions in space. The essence of the techniques comes down to how to sample good paths which contributes much to the final radiance at the camera. So may paths would not connect the camera and the light sources and thus would not contribute to the final radiance at the camera. I use the term "radiance" to mean the radiant energy per unit area per solid angle, and is the same as the "intensity" often used in the literature. NOW: As far as I gathered information from various sources, the yt volume renderer implements only a special case of a radiative transfer equation. It seems that it does not consider the second term ( the integration) in equation (3). Or it simplifies the integration int_{x_S, x} tau(x', x) [ LeV_0 (x', w) + sigma_s(x') int_{Sp } f_p( w, x', w') dw' ] dx' to int_{x_S, x} tau(x', x) [ LeV_0 (x', w) ] dx'. That is, it does not support the phase function f_p(w, x', w'). But I would like to use the phase function. In that case, we may try a "single scattering" version or a "multiple scattering" version. In the first version, each sampled position in the volume scatters only the radiance coming directly from the light source. In the second version, each sampled position in the volume scatters the radiances scattered by other positions. In this case, two or three multiple scattering would be good enough. But for me the single scattering version would be just fine. ----- The project I am interested in is rainbow simulation. The phase functions of a collection of water drops will be computed by means of Mie theory or some other variants. The implementation in C++ of path tracers are available in open source in the following site: http://www.mitsuba-renderer.org/devblog/?p=473 This implementation is reliable and tested. I am willing to add the relevant parts of the code into the yt volume renderer. But as is well known, it is not easy to find where to hack the existing code.... I think I need to cooperate with some one who knows the internal structure of the volume rendering code. Sincerely Moon R. Jung Professor Dept of Media Technology School of Media Sogang University. Science is the refinement of everyday thinking. Zen is your everyday life Life is at root playing
Dear Moon-Ryul, Sorry for the delay in replying. Your email took a bit of time to think about and absorb, and I wanted to give the best answer I could. As background, in this paper we go into detail on the main method by which yt volume renders: http://adsabs.harvard.edu/abs/2011ApJS..192....9T The underlying code has changed since that paper has been published, but I believe the methods in there are still accurate. For starters, as it stands, the yt method of integration currently starts from fixed positions and fixed vectors and then integrates outward. This means that regardless of how the values change as a ray moves over a cell, the ray always carries with it a source point and a vector. The next item that should be considered is that there *is* an experimental routine (likely non-functional) that dynamically adds and removes rays during integration; this also carries with it a modification of the *direction* of a ray, but not the source. This can be found in the "adaptive healpix" code, but again I suspect it is non-functional. (Non-adaptive healpix *is* functional, however.) So with that out of the way, I see only one reason that as a ray moves across the domain can't change both direction and source. The reason is the integration order of the volumetric patches (we call them "bricks") that it traverses. As it moves, if a ray changes source or direction, it could potentially re-enter a brick it has previously traversed. Sam Skillman's streamlines module nicely addresses this, but for a different purpose. The scalability is also different, because the streamlines seek to retain all traversed values, whereas volume rendering is a strictly accumulating process. So I think what you are describing may be *possible* with yt, but it may be difficult. That being said, there has been considerable interest in adding in a monte carlo module to yt; I've explored this a bit in the past and I think that it would be quite feasible and potentially quite scalable. On Thu, Nov 15, 2012 at 8:56 PM, Moon-Ryul Jung <moon@sogang.ac.kr> wrote:
Hi, everybody and Matthew:
I happened to encounter the blog of Matthew in the internet, and left some notes in the blog. I was led to this mailing list with the following request:
I like Python so much that I would like to do my experimentation with radiative transfer with Python environment. So, yt is a perfect environment. But it does not support a full radiative transfer equation yet.
By a full radiative transfer equation I mean:
(1) dot (w, del L (x, w) ) = LeV_0 (x,w) + sigma_s (x) [ int_{ Sp} f_p( w, x, w' ) dw' ] - simga_t(x) L(x, w)
for all x in a given volume V where V_0 is the "free " part of the volume V in which objects are not occupied.
Here del is the spatial gradient operator, L(x,w) is the radiance at point x in the direction of w. LeV_0 (x,w) is the emitting radiance at point x in the direction of w. sigma_s(x) is the scattering coefficient at x. Sp is the set of all solid angles. f_p(w, x, w') is the phase function at x, from direction w' to direction w. sigma_t (x) is the extinction or attenuation coefficient, and sigma_t = sigma_s + sigma_a, where sigma_a is the absorption coefficient.
(2) boundary conditions: for all x on surfaces of objects within the volume V:
L(x,w) = LeS (x, w) + int_{Sp} f_s (w, x, w') L(x, w') dot ( N(x), w' ) dw'
Here LeS(x,w) is the emitting radiace at x on a surface to direction w. f_s(w, x, w') is the reflection function at x from direction w' to direction w.
By combining (1) and (2), we get a single integral equation:
(3) L(x, w) = rau( x_S, w) L (x_S, w) + int_{x_S, x} tau(x', x) [ LeV_0 (x', w) + sigma_s(x') int_{Sp } f_p( w, x', w') dw' ] dx'
Here tau(x', x) is the optical depth from x' to x:
tau(x', x) = exp [ -int_{x, x'} sigma_t (x'') dx'' ]
It is not easy to find the radiance function L(x,w) that satisfies the integral equation of (3), where L occurs in both sides of the equation.
In computer graphics [within computer science], this problem has been solved mainly by means of Monte Carlo Integration technique, in particular "path tracing" techniques. In these techniques, we randomly generate a set of paths from the camera to the light sources, and compute the solution L by averaging the radiances carried by each path. A path is a sequence of positions in space. The essence of the techniques comes down to how to sample good paths which contributes much to the final radiance at the camera. So may paths would not connect the camera and the light sources and thus would not contribute to the final radiance at the camera.
I use the term "radiance" to mean the radiant energy per unit area per solid angle, and is the same as the "intensity" often used in the literature.
NOW: As far as I gathered information from various sources, the yt volume renderer implements only a special case of a radiative transfer equation. It seems that it does not consider the second term ( the integration) in equation (3). Or it simplifies the integration int_{x_S, x} tau(x', x) [ LeV_0 (x', w) + sigma_s(x') int_{Sp } f_p( w, x', w') dw' ] dx'
to
int_{x_S, x} tau(x', x) [ LeV_0 (x', w) ] dx'.
That is, it does not support the phase function f_p(w, x', w').
Yes, that is correct -- we do not phase shift.
But I would like to use the phase function. In that case, we may try a "single scattering" version or a "multiple scattering" version. In the first version, each sampled position in the volume scatters only the radiance coming directly from the light source. In the second version, each sampled position in the volume scatters the radiances scattered by other positions. In this case, two or three multiple scattering would be good enough. But for me the single scattering version would be just fine.
We may be able to support this, but I am not entirely sure. If I am to understand correctly, this would imply that at a certain critical value or threshold, a given ray emanating from a source would be scattered and be sent in a new direction. I'm thinking about this and I think that if the ray vectors are constant in number (i.e., we do not add new rays) then this could be done quite simply by swapping out the integrator in the Streamlines, which would be simpler than the method of using the volume integrator. What this says to me is that we could actually gain quite a bit of interesting functionality if the streamline integrator were refactored to fit better with the volume rendering -- separation of accumulation steps from the integration steps, and allowing both to be swapped out independently. Sam, what do you think of this?
----- The project I am interested in is rainbow simulation. The phase functions of a collection of water drops will be computed by means of Mie theory or some other variants.
The implementation in C++ of path tracers are available in open source in the following site:
http://www.mitsuba-renderer.org/devblog/?p=473
This implementation is reliable and tested. I am willing to add the relevant parts of the code into the yt volume renderer.
WOW! That is quite cool. This owuld be far, far simpler than what I mentioned above, with the ray tracing. In fact, reading the Mitsuba page, it looks like it already accepts hierarchical grids, which is what yt provides. This would be amazing ,and I think mapping the data structures could be straightforward. I think the "gridvolume" data source (p116 of the documentation) would be a good match. There are already python bindings, so I think what would work best would be: 1) Using python (no yt) supply a single grid to Mitsuba and render this 2) Using yt, use the amr kD-Tree to decompose the data 3) Using yt, supply *many* bricks from the kD-tree to Mitsuba and render these This would be an amazing addition to the code. I'm happy to help you with identifying where the right places are to hook up Mitsuba and yt. I'd never heard of Mitsuba before, but it's really impressive.
But as is well known, it is not easy to find where to hack the existing code....
I think I need to cooperate with some one who knows the internal structure of the volume rendering code.
Sure! I'm happy to serve as that person. I'd be eager to help out with this. Let me know if there's anything else I can do to help you get started, and please let us know how it goes if you decide to take this on! -Matt
Sincerely
Moon R. Jung Professor Dept of Media Technology School of Media Sogang University.
Science is the refinement of everyday thinking. Zen is your everyday life Life is at root playing
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (3)
-
Anna Rosen -
Matthew Turk -
Moon-Ryul Jung