Mailman 3 August 2015 - yt-dev

Issue #1072: Add support to YTArray for remaining unsupported ufuncs (yt_analysis/yt)
by Nathan Goldbaum 24 Aug '15

24 Aug '15

New issue 1072: Add support to YTArray for remaining unsupported ufuncs https://bitbucket.org/yt_analysis/yt/issues/1072/add-support-to-ytarray-for… Nathan Goldbaum: Right now YTArray supports most numpy ufuncs, but it looks like there are a few that need to be implemented still. * bitwise_and * bitwise_or * invert * left_shift * right_shift * ldexp

1 0

PR Triage Hangout + Poll for Fall!
by Hilary Egan 24 Aug '15

24 Aug '15

Hey all, Here's this weeks hangout link, and a poll for what time we want these meetings to occur in the fall. Even if you haven't switched over to the fall semester quite yet, just treat this poll as though it were asking about your average week through the fall. Thanks! -Hilary *Poll:* http://doodle.com/mmyecmeweih7z382 *Hangout: * https://plus.google.com/b/107728486871834552760/events/c2vj4fn4oomve6lso1i4…

1 1

Issue #1071: ART does not auto-detect AGORA files (yt_analysis/yt)
by Kenza Arraki 23 Aug '15

23 Aug '15

New issue 1071: ART does not auto-detect AGORA files https://bitbucket.org/yt_analysis/yt/issues/1071/art-does-not-auto-detect-a… Kenza Arraki: ART cannot automatically load the AGORA files. This causes it to fail as it does not know the number of root grid cells. Data from: `AGORA-DISK-repository-for-use/Grackle+noSF/ART-I/IC` ``` #!python In [1]: import yt In [2]: ds = yt.load('AGORA_Galaxy_LOW.d') yt : [INFO ] 2015-08-23 11:21:45,490 Using root level of 09 yt : [INFO ] 2015-08-23 11:21:45,677 Max level is 08 --------------------------------------------------------------------------- KeyError Traceback (most recent call last) <ipython-input-2-661a3253e356> in <module>() ----> 1 ds = yt.load('AGORA_Galaxy_LOW.d') ... /global/u1/k/karraki/yt/yt/frontends/art/data_structures.pyc in _set_code_unit_attributes(self) 231 wmu = self.parameters["wmu"] 232 Om0 = self.parameters['Om0'] --> 233 ng = self.parameters['ng'] 234 boxh = self.parameters['boxh'] 235 aexpn = self.parameters["aexpn"] KeyError: 'ng' ``` Responsible: karraki

1 0

Backporting bugfixes and possible changes to development practices
by Nathan Goldbaum 21 Aug '15

21 Aug '15

Hi all, I spent some time this week coming up with a technical solution to our social problem of doing bugfix releases. This is embodied in a new script, "pr_backport.py", which I've created a pull request for here to include with the yt repository. https://bitbucket.org/yt_analysis/yt/pull-requests/1717/add-pr-backport-scr… Briefly, the script creates a temporary clone of yt_analysis/yt, figures out which commits on the "yt" branch need to be included on the stable branch, and then assigns those commits to a pull request, and then interactively suggests mercurial commands to run in the temporary repository to backport individual pull requests to the stable branch. This is all done "manually" by copy/pasting mercurial commands suggested by the script. I don't think this should be fully automated since merge conflicts can and likely will be triggered, particularly if the stable and yt branches begin to substantially diverge. My hope is that this will make the task of backporting bugfixes, which up until this point has been associated with such a high cognitive barrier that no one has wanted to sit down and acutally do it, much easier. I might be biased since I wrote the script, but using the final version that I pull requested it took me about 10 minutes or so to go through the list of pull requests and backport them to create this pull request: https://bitbucket.org/yt_analysis/yt/pull-requests/1716/backporting-bugfixe… Being able to backport bugfixes like this allows us to possibly make some changes to our development workflow and suggested best practices. 1. We can now suggest most users stay on the stable branch rather than switching to the dev branch to get bugfixes. We can also issue bugfix releases on a much quicker cadence (monthly? biweekly?). 2. Since the "stable" branch is more stable and less buggy and the dev branch is basically now just for new features, perhaps we can allow more experimental development on the main yt branch without an onerous requirement of documentation and review. This will help unblock the development of major features like support for unstructured mesh data and the new volume rendering interface. 3. Since there is less pressure to do a feature release just to get bugfixes out, we might be able to have longer intervals between feature releases, giving more time to shake out issues. I'm more or less just throwing these thoughts out there to get feedback from the development community. Does the easy availability of bugfix releases constitute a significant change for our development and release mindset? Thanks all for your time and input! -Nathan

3 2

Issue #1070: Cannot access field that needs ghost zones from a covering grid (yt_analysis/yt)
by Nathan Goldbaum 21 Aug '15

21 Aug '15

New issue 1070: Cannot access field that needs ghost zones from a covering grid https://bitbucket.org/yt_analysis/yt/issues/1070/cannot-access-field-that-n… Nathan Goldbaum: The following test script: ``` #!python import yt ds = yt.load('GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0100') ad = ds.all_data() ml = ad.quantities.min_location('gpot') center = (ml[2], ml[3], ml[4]) dx = ds.index.get_smallest_dx() width = 256*dx LE = center-width.in_units('code_length')/2. CN = int(width/dx) Roi = ds.covering_grid(level=ds.index.max_level, left_edge=LE, dims=[CN, CN, CN]) Roi['gas', 'vorticity_x'] ``` Fails with the following traceback ``` Traceback (most recent call last): File "test2.py", line 20, in <module> Roi['gas', 'vorticity_x'] File "/Users/goldbaum/Documents/yt-hg/yt/data_objects/data_containers.py", line 253, in __getitem__ self.get_data(f) File "/Users/goldbaum/Documents/yt-hg/yt/data_objects/construction_data_containers.py", line 598, in get_data fill, gen, part, alias = self._split_fields(fields_to_get) File "/Users/goldbaum/Documents/yt-hg/yt/data_objects/construction_data_containers.py", line 616, in _split_fields finfo.check_available(self) File "/Users/goldbaum/Documents/yt-hg/yt/fields/derived_field.py", line 150, in check_available validator(data) File "/Users/goldbaum/Documents/yt-hg/yt/fields/derived_field.py", line 293, in __call__ raise NeedsGridType(self.ghost_zones,self.fields) yt.fields.field_exceptions.NeedsGridType: (1, [('gas', 'velocity_x'), ('gas', 'velocity_y'), ('gas', 'velocity_z')]) ``` Thanks to Markus Haider for reporting this on the mailing list.

1 0

Issue #1069: Data objects used to create a profile in a for loop get garbage collected (yt_analysis/yt)
by Nathan Goldbaum 20 Aug '15

20 Aug '15

New issue 1069: Data objects used to create a profile in a for loop get garbage collected https://bitbucket.org/yt_analysis/yt/issues/1069/data-objects-used-to-creat… Nathan Goldbaum: Illustrated by a test script Mateusz Ruszkowski reported on the mailing list: ``` #!python import yt profiles = [] labels = [] for i in range(0,2,1): print i ds = yt.load("Sedov_3D/sedov_hdf5_chk_%04i" % i) my_data = ds.sphere(ds.domain_center, 1.0) sim_time = ds.current_time profiles.append(yt.create_profile(my_data, "radius", fields=["density"], weight_field="cell_volume", n_bins=32)) labels.append("t= %.0f" % sim_time) plot = yt.ProfilePlot.from_profiles(profiles, labels=labels) plot.set_xlim([0.1, 1.0]) plot.save() ``` This uses the Sedov3D test dataset on yt-project.org/data.

1 0

Issue #1068: Cannot use ipython notebook with anaconda once yt is installed (yt_analysis/yt)
by ggoodwin52 20 Aug '15

20 Aug '15

New issue 1068: Cannot use ipython notebook with anaconda once yt is installed https://bitbucket.org/yt_analysis/yt/issues/1068/cannot-use-ipython-noteboo… ggoodwin52: I installed the most updated version of anaconda for linux-64bit and was able to open ipython, ipython notebook, spyder, etc. without any issue. Then, I installed yt using the "conda install yt" command. The install was successful, but now when I try and start ipython notebook or spyder to use yt I get an error message. The following script from my terminal shows successful starting of the ipython notebook via anaconda, closing the ipython notebook, installing yt, and then trying to re-open the ipython notebook and getting an error. ``` #!bash [ggoodwin@Omoikane ~]$ ipython notebook [I 13:36:38.296 NotebookApp] Using MathJax from CDN: https://cdn.mathjax.org/mathjax/latest/MathJax.js [I 13:36:38.376 NotebookApp] The port 8888 is already in use, trying another random port. [I 13:36:38.376 NotebookApp] The port 8889 is already in use, trying another random port. [I 13:36:38.377 NotebookApp] The port 8890 is already in use, trying another random port. [I 13:36:38.385 NotebookApp] Serving notebooks from local directory: /home/ggoodwin [I 13:36:38.385 NotebookApp] 0 active kernels [I 13:36:38.386 NotebookApp] The IPython Notebook is running at: http://localhost:8891/ [I 13:36:38.386 NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation). ^C[I 13:36:40.027 NotebookApp] interrupted Serving notebooks from local directory: /home/ggoodwin 0 active kernels The IPython Notebook is running at: http://localhost:8891/ Shutdown this notebook server (y/[n])? y [C 13:36:41.525 NotebookApp] Shutdown confirmed [I 13:36:41.525 NotebookApp] Shutting down kernels [ggoodwin@Omoikane ~]$ conda install yt Fetching package metadata: .... Solving package specifications: . Package plan for installation in environment /home/ggoodwin/anaconda: The following packages will be downloaded: package | build ---------------------------|----------------- openblas-0.2.14 | 0 6.6 MB conda-env-2.4.2 | py27_0 24 KB decorator-4.0.2 | py27_0 11 KB ipython_genutils-0.1.0 | py27_0 32 KB numpy-1.9.2 | py27_1 5.7 MB path.py-7.6 | py27_0 41 KB pexpect-3.3 | py27_0 60 KB pickleshare-0.5 | py27_0 8 KB setuptools-18.1 | py27_0 341 KB simplegeneric-0.8.1 | py27_0 6 KB traitlets-4.0.0 | py27_0 88 KB wheel-0.24.0 | py27_0 75 KB conda-3.16.0 | py27_0 172 KB ipython-4.0.0 | py27_0 915 KB pip-7.1.0 | py27_1 1.4 MB yt-3.2 | np19py27_0 8.0 MB ------------------------------------------------------------ Total: 23.3 MB The following NEW packages will be INSTALLED: ipython_genutils: 0.1.0-py27_0 openblas: 0.2.14-0 path.py: 7.6-py27_0 pexpect: 3.3-py27_0 pickleshare: 0.5-py27_0 simplegeneric: 0.8.1-py27_0 traitlets: 4.0.0-py27_0 wheel: 0.24.0-py27_0 yt: 3.2-np19py27_0 The following packages will be UPDATED: conda: 3.14.1-py27_0 --> 3.16.0-py27_0 conda-env: 2.2.3-py27_0 --> 2.4.2-py27_0 decorator: 3.4.2-py27_0 --> 4.0.2-py27_0 ipython: 3.2.0-py27_0 --> 4.0.0-py27_0 numpy: 1.9.2-py27_0 --> 1.9.2-py27_1 pip: 7.0.3-py27_0 --> 7.1.0-py27_1 setuptools: 17.1.1-py27_0 --> 18.1-py27_0 Proceed ([y]/n)? y Fetching packages ... openblas-0.2.1 100% |################################| Time: 0:00:00 18.01 MB/s conda-env-2.4. 100% |################################| Time: 0:00:00 212.81 kB/s decorator-4.0. 100% |################################| Time: 0:00:00 99.67 kB/s ipython_genuti 100% |################################| Time: 0:00:00 274.77 kB/s numpy-1.9.2-py 100% |################################| Time: 0:00:00 7.59 MB/s path.py-7.6-py 100% |################################| Time: 0:00:00 2.65 MB/s pexpect-3.3-py 100% |################################| Time: 0:00:00 164.63 kB/s pickleshare-0. 100% |################################| Time: 0:00:00 68.05 kB/s setuptools-18. 100% |################################| Time: 0:00:00 9.75 MB/s simplegeneric- 100% |################################| Time: 0:00:00 4.05 MB/s traitlets-4.0. 100% |################################| Time: 0:00:00 508.57 kB/s wheel-0.24.0-p 100% |################################| Time: 0:00:00 3.54 MB/s conda-3.16.0-p 100% |################################| Time: 0:00:00 1.66 MB/s ipython-4.0.0- 100% |################################| Time: 0:00:00 5.93 MB/s pip-7.1.0-py27 100% |################################| Time: 0:00:00 5.45 MB/s yt-3.2-np19py2 100% |################################| Time: 0:00:00 8.83 MB/s Extracting packages ... [ COMPLETE ]|###################################################| 100% Unlinking packages ... [ COMPLETE ]|###################################################| 100% Linking packages ... [ COMPLETE ]|###################################################| 100% [ggoodwin@Omoikane ~]$ ipython notebook Traceback (most recent call last): File "/home/ggoodwin/anaconda/bin/ipython", line 6, in <module> sys.exit(start_ipython()) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/IPython/__init__.py", line 118, in start_ipython return launch_new_instance(argv=argv, **kwargs) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 591, in launch_instance app.initialize(argv) File "<string>", line 2, in initialize File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 75, in catch_config_error return method(app, *args, **kwargs) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/IPython/terminal/ipapp.py", line 305, in initialize super(TerminalIPythonApp, self).initialize(argv) File "<string>", line 2, in initialize File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 75, in catch_config_error return method(app, *args, **kwargs) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/IPython/core/application.py", line 386, in initialize self.parse_command_line(argv) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/IPython/terminal/ipapp.py", line 300, in parse_command_line return super(TerminalIPythonApp, self).parse_command_line(argv) File "<string>", line 2, in parse_command_line File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 75, in catch_config_error return method(app, *args, **kwargs) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 487, in parse_command_line return self.initialize_subcommand(subc, subargv) File "<string>", line 2, in initialize_subcommand File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 75, in catch_config_error return method(app, *args, **kwargs) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/traitlets/config/application.py", line 418, in initialize_subcommand subapp = import_item(subapp) File "/home/ggoodwin/anaconda/lib/python2.7/site-packages/ipython_genutils/importstring.py", line 31, in import_item module = __import__(package, fromlist=[obj]) ImportError: No module named notebook.notebookapp [ggoodwin@Omoikane ~]$ ```

1 0

PR to remove Theia
by Nathan Goldbaum 20 Aug '15

20 Aug '15

Hi all, As part of my ongoing efforts to clean up the codebase, I just opened a pull request to remove the Theia hardware volume rendering interface from yt: https://bitbucket.org/yt_analysis/yt/pull-requests/1703/removing-theia-hard… I'd only like to see this pull request merged if there's consensus that this is the way to go. I'd especially like to hear from anyone that is actually using this code, or if anyone is interested in maintaining it and adding tests. -Nathan

4 3

Issue #1067: ds.h.find_max fail with cylindrical data (yt_analysis/yt)
by Suoqing Ji 19 Aug '15

19 Aug '15

New issue 1067: ds.h.find_max fail with cylindrical data https://bitbucket.org/yt_analysis/yt/issues/1067/dshfind_max-fail-with-cyli… Suoqing Ji: Test [data](http://yt-project.org/data/WDMerger_hdf5_chk_1000.tar.gz) and script: ``` #!python import yt ds = yt.load("WDmerger_hdf5_chk_1000/WDmerger_hdf5_chk_1000.hdf5") v, c = ds.h.find_max('density') ``` Error message: ``` #!bash v, c = ds.h.find_max('density') File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/static_output.py", line 583, in find_max source.quantities.max_location(field) File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 540, in __call__ rv = super(MaxLocation, self).__call__(field) File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 58, in __call__ sto.result = self.process_chunk(ds, *args, **kwargs) File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/derived_quantities.py", line 547, in process_chunk mx = array_like_field(data, -1, "x") File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/units/yt_array.py", line 1271, in array_like_field field = data._determine_fields(field)[0] File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/data_containers.py", line 555, in _determine_fields finfo = self.ds._get_field_info("unknown", fname) File "/Users/jisuoqing/Workspace/code/yt-x86_64/src/yt-hg/yt/data_objects/static_output.py", line 543, in _get_field_info raise YTFieldNotFound((ftype, fname), self) yt.utilities.exceptions.YTFieldNotFound: Could not find field '('all', 'x')' in WDmerger_hdf5_chk_1000/WDmerger_hdf5_chk_1000.hdf5. ``` Although I'm kind of know what causes this error, someone else might have a better and more global idea about this. Also I'm in no hurry to use this function -- it just appears and I report it here.

1 0

PR prematurely merged
by Nathan Goldbaum 19 Aug '15

19 Aug '15

Hey Ben, What happened here? Did you mean to merge this PR? I think it was still marked as work in progress, e.g. [WIP], and Billi had more commits to push to it. -Nathan ---------- Forwarded message ---------- From: <commits-noreply(a)bitbucket.org> Date: Wed, Aug 19, 2015 at 12:09 PM Subject: [yt-svn] commit/yt: bwkeller: Merged in qobilidop/yt (pull request #1670) To: yt-svn(a)lists.spacepope.org 1 new commit in yt: https://bitbucket.org/yt_analysis/yt/commits/cfc98fddd9d7/ Changeset: cfc98fddd9d7 Branch: yt User: bwkeller Date: 2015-08-19 17:09:19+00:00 Summary: Merged in qobilidop/yt (pull request #1670) Alternative Smoothing Kernels Affected #: 11 files diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/construction_data_containers.py --- a/yt/data_objects/construction_data_containers.py +++ b/yt/data_objects/construction_data_containers.py @@ -682,11 +682,13 @@ def RightEdge(self): return self.right_edge - def deposit(self, positions, fields = None, method = None): + def deposit(self, positions, fields = None, method = None, + kernel_name = 'cubic'): cls = getattr(particle_deposit, "deposit_%s" % method, None) if cls is None: raise YTParticleDepositionNotImplemented(method) - op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs + # We allocate number of zones, not number of octs + op = cls(self.ActiveDimensions.prod(), kernel_name) op.initialize() op.process_grid(self, positions, fields) vals = op.finalize() @@ -1787,5 +1789,3 @@ else: mylog.error("Problem uploading.") return upload_id - - diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/grid_patch.py --- a/yt/data_objects/grid_patch.py +++ b/yt/data_objects/grid_patch.py @@ -330,12 +330,14 @@ def particle_operation(self, *args, **kwargs): raise NotImplementedError - def deposit(self, positions, fields = None, method = None): + def deposit(self, positions, fields = None, method = None, + kernel_name = 'cubic'): # Here we perform our particle deposition. cls = getattr(particle_deposit, "deposit_%s" % method, None) if cls is None: raise YTParticleDepositionNotImplemented(method) - op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs + # We allocate number of zones, not number of octs + op = cls(self.ActiveDimensions.prod(), kernel_name) op.initialize() op.process_grid(self, positions, fields) vals = op.finalize() diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/octree_subset.py --- a/yt/data_objects/octree_subset.py +++ b/yt/data_objects/octree_subset.py @@ -141,7 +141,8 @@ self._domain_ind = di return self._domain_ind - def deposit(self, positions, fields = None, method = None): + def deposit(self, positions, fields = None, method = None, + kernel_name='cubic'): r"""Operate on the mesh, in a particle-against-mesh fashion, with exclusively local input. @@ -176,7 +177,8 @@ raise YTParticleDepositionNotImplemented(method) nz = self.nz nvals = (nz, nz, nz, (self.domain_ind >= 0).sum()) - op = cls(nvals) # We allocate number of zones, not number of octs + # We allocate number of zones, not number of octs + op = cls(nvals, kernel_name) op.initialize() mylog.debug("Depositing %s (%s^3) particles into %s Octs", positions.shape[0], positions.shape[0]**0.3333333, nvals[-1]) @@ -192,7 +194,8 @@ return np.asfortranarray(vals) def smooth(self, positions, fields = None, index_fields = None, - method = None, create_octree = False, nneighbors = 64): + method = None, create_octree = False, nneighbors = 64, + kernel_name = 'cubic'): r"""Operate on the mesh, in a particle-against-mesh fashion, with non-local input. @@ -258,7 +261,7 @@ nz = self.nz mdom_ind = self.domain_ind nvals = (nz, nz, nz, (mdom_ind >= 0).sum()) - op = cls(nvals, len(fields), nneighbors) + op = cls(nvals, len(fields), nneighbors, kernel_name) op.initialize() mylog.debug("Smoothing %s particles into %s Octs", positions.shape[0], nvals[-1]) @@ -280,7 +283,7 @@ return vals def particle_operation(self, positions, fields = None, - method = None, nneighbors = 64): + method = None, nneighbors = 64, kernel_name = 'cubic'): r"""Operate on particles, in a particle-against-particle fashion. This uses the octree indexing system to call a "smoothing" operation @@ -335,7 +338,7 @@ nz = self.nz mdom_ind = self.domain_ind nvals = (nz, nz, nz, (mdom_ind >= 0).sum()) - op = cls(nvals, len(fields), nneighbors) + op = cls(nvals, len(fields), nneighbors, kernel_name) op.initialize() mylog.debug("Smoothing %s particles into %s Octs", positions.shape[0], nvals[-1]) diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/data_objects/unstructured_mesh.py --- a/yt/data_objects/unstructured_mesh.py +++ b/yt/data_objects/unstructured_mesh.py @@ -105,13 +105,15 @@ def select_tcoords(self, dobj): raise NotImplementedError - def deposit(self, positions, fields = None, method = None): + def deposit(self, positions, fields = None, method = None, + kernel_name = 'cubic'): raise NotImplementedError # Here we perform our particle deposition. cls = getattr(particle_deposit, "deposit_%s" % method, None) if cls is None: raise YTParticleDepositionNotImplemented(method) - op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs + # We allocate number of zones, not number of octs + op = cls(self.ActiveDimensions.prod(), kernel_name) op.initialize() op.process_grid(self, positions, fields) vals = op.finalize() @@ -200,4 +202,3 @@ else: self._last_count = mask.sum() return mask - diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/fields/particle_fields.py --- a/yt/fields/particle_fields.py +++ b/yt/fields/particle_fields.py @@ -766,8 +766,12 @@ def add_volume_weighted_smoothed_field(ptype, coord_name, mass_name, smoothing_length_name, density_name, smoothed_field, registry, - nneighbors = None): - field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field)) + nneighbors = None, kernel_name = 'cubic'): + if kernel_name == 'cubic': + field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field)) + else: + field_name = ("deposit", "%s_%s_smoothed_%s" % (ptype, kernel_name, + smoothed_field)) field_units = registry[ptype, smoothed_field].units def _vol_weight(field, data): pos = data[ptype, coord_name] @@ -789,7 +793,8 @@ # volume_weighted smooth operations return lists of length 1. rv = data.smooth(pos, [mass, hsml, dens, quan], method="volume_weighted", - create_octree=True)[0] + create_octree=True, + kernel_name=kernel_name)[0] rv[np.isnan(rv)] = 0.0 # Now some quick unit conversions. rv = data.apply_units(rv, field_units) @@ -816,8 +821,12 @@ units = "code_length") return [field_name] -def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64): - field_name = (ptype, "smoothed_density") +def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64, + kernel_name = 'cubic'): + if kernel_name == 'cubic': + field_name = (ptype, "smoothed_density") + else: + field_name = (ptype, "%s_smoothed_density" % (kernel_name)) field_units = registry[ptype, mass_name].units def _nth_neighbor(field, data): pos = data[ptype, coord_name] @@ -827,7 +836,8 @@ densities = mass * 0.0 data.particle_operation(pos, [mass, densities], method="density", - nneighbors = nneighbors) + nneighbors = nneighbors, + kernel_name = kernel_name) ones = pos.prod(axis=1) # Get us in code_length**3 ones[:] = 1.0 densities /= ones diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/frontends/artio/data_structures.py --- a/yt/frontends/artio/data_structures.py +++ b/yt/frontends/artio/data_structures.py @@ -133,7 +133,8 @@ tr = dict((field, v) for field, v in zip(fields, tr)) return tr - def deposit(self, positions, fields = None, method = None): + def deposit(self, positions, fields = None, method = None, + kernel_name = 'cubic'): # Here we perform our particle deposition. if fields is None: fields = [] cls = getattr(particle_deposit, "deposit_%s" % method, None) @@ -141,7 +142,8 @@ raise YTParticleDepositionNotImplemented(method) nz = self.nz nvals = (nz, nz, nz, self.ires.size) - op = cls(nvals) # We allocate number of zones, not number of octs + # We allocate number of zones, not number of octs + op = cls(nvals, kernel_name) op.initialize() mylog.debug("Depositing %s (%s^3) particles into %s Root Mesh", positions.shape[0], positions.shape[0]**0.3333333, nvals[-1]) diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/frontends/artio/setup.py --- a/yt/frontends/artio/setup.py +++ b/yt/frontends/artio/setup.py @@ -19,7 +19,8 @@ depends=artio_sources + ["yt/utilities/lib/fp_utils.pxd", "yt/geometry/oct_container.pxd", - "yt/geometry/selection_routines.pxd"]) + "yt/geometry/selection_routines.pxd", + "yt/geometry/particle_deposit.pxd"]) config.make_config_py() # installs __config__.py #config.make_svn_version_py() return config diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_deposit.pxd --- a/yt/geometry/particle_deposit.pxd +++ b/yt/geometry/particle_deposit.pxd @@ -38,7 +38,7 @@ # Standard SPH kernel for use with the Grid method # #################################################### -cdef inline np.float64_t sph_kernel(np.float64_t x) nogil: +cdef inline np.float64_t sph_kernel_cubic(np.float64_t x) nogil: cdef np.float64_t kernel if x <= 0.5: kernel = 1.-6.*x*x*(1.-x) @@ -48,8 +48,55 @@ kernel = 0. return kernel +######################################################## +# Alternative SPH kernels for use with the Grid method # +######################################################## + +# quartic spline +cdef inline np.float64_t sph_kernel_quartic(np.float64_t x): + cdef np.float64_t kernel + cdef np.float64_t C = 5.**6/512/np.pi + if x < 1: + kernel = (1.-x)**4 + if x < 3./5: + kernel -= 5*(3./5-x)**4 + if x < 1./5: + kernel += 10*(1./5-x)**4 + else: + kernel = 0. + return kernel * C + +# quintic spline +cdef inline np.float64_t sph_kernel_quintic(np.float64_t x): + cdef np.float64_t kernel + cdef np.float64_t C = 3.**7/40/np.pi + if x < 1: + kernel = (1.-x)**5 + if x < 2./3: + kernel -= 6*(2./3-x)**5 + if x < 1./3: + kernel += 15*(1./3-x)**5 + else: + kernel = 0. + return kernel * C + +# I don't know the way to use a dict in a cdef class. +# So in order to mimic a registry functionality, +# I manually created a function to lookup the kernel functions. +ctypedef np.float64_t (*kernel_func) (np.float64_t) +cdef inline kernel_func get_kernel_func(str kernel_name): + if kernel_name == 'cubic': + return sph_kernel_cubic + elif kernel_name == 'quartic': + return sph_kernel_quartic + elif kernel_name == 'quintic': + return sph_kernel_quintic + else: + raise NotImplementedError + cdef class ParticleDepositOperation: # We assume each will allocate and define their own temporary storage + cdef kernel_func sph_kernel cdef public object nvals cdef public int update_values cdef void process(self, int dim[3], np.float64_t left_edge[3], diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_deposit.pyx --- a/yt/geometry/particle_deposit.pyx +++ b/yt/geometry/particle_deposit.pyx @@ -25,9 +25,10 @@ OctreeContainer, OctInfo cdef class ParticleDepositOperation: - def __init__(self, nvals): + def __init__(self, nvals, kernel_name): self.nvals = nvals self.update_values = 0 # This is the default + self.sph_kernel = get_kernel_func(kernel_name) def initialize(self, *args): raise NotImplementedError @@ -227,7 +228,7 @@ dist = idist[0] + idist[1] + idist[2] # Calculate distance in multiples of the smoothing length dist = sqrt(dist) / fields[0] - self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist) + self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist) kernel_sum += self.temp[gind(i,j,k,dim) + offset] # Having found the kernel, deposit accordingly into gdata for i from ib0[0] <= i <= ib1[0]: @@ -493,4 +494,3 @@ return self.onnfield deposit_nearest = NNParticleField - diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_smooth.pxd --- a/yt/geometry/particle_smooth.pxd +++ b/yt/geometry/particle_smooth.pxd @@ -22,7 +22,7 @@ from fp_utils cimport * from oct_container cimport Oct, OctAllocationContainer, OctreeContainer -from .particle_deposit cimport sph_kernel, gind +from .particle_deposit cimport kernel_func, get_kernel_func, gind cdef extern from "platform_dep.h": void *alloca(int) @@ -34,6 +34,7 @@ cdef class ParticleSmoothOperation: # We assume each will allocate and define their own temporary storage + cdef kernel_func sph_kernel cdef public object nvals cdef np.float64_t DW[3] cdef int nfields diff -r d5816147380b81418aec9ccc7ac0f551e8ddce21 -r cfc98fddd9d791464330adf166d4229749dab22c yt/geometry/particle_smooth.pyx --- a/yt/geometry/particle_smooth.pyx +++ b/yt/geometry/particle_smooth.pyx @@ -74,7 +74,7 @@ opos[2] = ipos[2] cdef class ParticleSmoothOperation: - def __init__(self, nvals, nfields, max_neighbors): + def __init__(self, nvals, nfields, max_neighbors, kernel_name): # This is the set of cells, in grids, blocks or octs, we are handling. cdef int i self.nvals = nvals @@ -83,6 +83,7 @@ self.neighbors = <NeighborList *> malloc( sizeof(NeighborList) * self.maxn) self.neighbor_reset() + self.sph_kernel = get_kernel_func(kernel_name) def initialize(self, *args): raise NotImplementedError @@ -630,7 +631,7 @@ # Usually this density has been computed dens = fields[2][pn] if dens == 0.0: continue - weight = mass * sph_kernel(sqrt(r2) / hsml) / dens + weight = mass * self.sph_kernel(sqrt(r2) / hsml) / dens # Mass of the particle times the value for fi in range(self.nfields - 3): val = fields[fi + 3][pn] @@ -756,7 +757,7 @@ for pn in range(self.curn): mass = fields[0][self.neighbors[pn].pn] r2 = self.neighbors[pn].r2 - lw = sph_kernel(sqrt(r2) / hsml) + lw = self.sph_kernel(sqrt(r2) / hsml) dens += mass * lw weight = (4.0/3.0) * 3.1415926 * hsml**3 fields[1][offset] = dens/weight Repository URL: https://bitbucket.org/yt_analysis/yt/ -- This is a commit notification from bitbucket.org. You are receiving this because you have the service enabled, addressing the recipient of this email. _______________________________________________ yt-svn mailing list yt-svn(a)lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-svn-spacepope.org

1 0