Octree ghost zones — reloaded #2425

cphyc · 2020-01-23T18:11:41Z

PR Summary

This PR adds support for ghost zones (or a similar interface). This is another attempt than #2166 based on the OctVisitors machinery.

The good

import numpy as np
import yt
ds = yt.load('output_00080/info_00080.txt')

ad = ds.all_data()
center = ds.arr(ad.argmax('density'))
sp = ds.sphere(center, (100, 'kpc'))
ds.add_gradient_fields(('gas', 'density'))

center = np.array([.1, .5, .5])
w = np.array([.05, 1, 1])
reg = ds.box(center-w/2, center+w/2)
p = yt.ProjectionPlot(ds, 'x', ['density'] + [f'density_gradient_{k}' for k in 'xyz'], center=center, data_source=reg, weight_field='density')
p.set_axes_unit('Mpc')

for k in 'xyz':
    p.set_log(f'density_gradient_{k}', True, linthresh=1e-55)
    p.set_zlim(f'density_gradient_{k}', -1e-53, 1e-53)
    p.set_cmap(f'density_gradient_{k}', 'bwr')
p.save('frames/the_good/')

The bad

ad = ds.all_data()
center = ds.arr(ad.argmax('density'))
sp = ds.sphere(center, (100, 'kpc'))
ds.add_gradient_fields(('gas', 'density'))

p = yt.SlicePlot(ds, 'x', ['density'] + [f'density_gradient_{k}' for k in 'xyz'], width=sp.radius, center=center, data_source=sp)
p.set_axes_unit('kpc')

dens = p.frb['density']


p.annotate_clear()

for k in 'xyz':
    p.set_log(f'density_gradient_{k}', True, linthresh=1e-48)
    p.set_zlim(f'density_gradient_{k}', -1e-45, 1e-45)
    p.set_cmap(f'density_gradient_{k}', 'bwr')

p.annotate_domain_edges()
p.save('frames/the_bad/')

matthewturk · 2020-01-23T20:10:28Z

I dig the fstrings in your example.

cphyc · 2020-01-30T13:37:03Z

Current version:

The good

import numpy as np
import yt
ds = yt.load('output_00080/info_00080.txt')

ad = ds.all_data()
center = ds.arr(ad.argmax('density'))
sp = ds.sphere(center, (100, 'kpc'))
ds.add_gradient_fields(('gas', 'density'))

center = np.array([.1, .5, .5])
w = np.array([.05, 1, 1])
reg = ds.box(center-w/2, center+w/2)
p = yt.ProjectionPlot(ds, 'x', ['density'] + [f'density_gradient_{k}' for k in 'xyz'], center=center, data_source=reg, weight_field='density')
p.set_axes_unit('Mpc')

for k in 'xyz':
    p.set_log(f'density_gradient_{k}', True, linthresh=1e-55)
    p.set_zlim(f'density_gradient_{k}', -1e-53, 1e-53)
    p.set_cmap(f'density_gradient_{k}', 'bwr')
p.save('frames/the_good/')

The bad

ad = ds.all_data()
center = ds.arr(ad.argmax('density'))
sp = ds.sphere(center, (100, 'kpc'))
ds.add_gradient_fields(('gas', 'density'))

p = yt.SlicePlot(ds, 'x', ['density'] + [f'density_gradient_{k}' for k in 'xyz'], width=sp.radius, center=center, data_source=sp)
p.set_axes_unit('kpc')

dens = p.frb['density']


p.annotate_clear()

for k in 'xyz':
    p.set_log(f'density_gradient_{k}', True, linthresh=1e-48)
    p.set_zlim(f'density_gradient_{k}', -1e-45, 1e-45)
    p.set_cmap(f'density_gradient_{k}', 'bwr')

p.annotate_domain_edges()
p.save('frames/the_bad/')

yt/fields/fluid_fields.py

yt/data_objects/octree_subset.py

yt/fields/fluid_fields.py

matthewturk

This looks pretty good, but I'm not sure that I've hit all the corner cases. Nice!

setup.py

yt/fields/fluid_fields.py

yt/geometry/oct_container.pyx

cphyc · 2020-02-05T16:16:16Z

@yt-fido test this please

cphyc · 2020-02-05T17:06:32Z

This should be ready to go!

brittonsmith

This all looks good to me, but I think it would good to see the answer tests run on this before it's merged since the gradient fields are changed. We could either include this in the PR #2437 or wait until that's in. @matthewturk, what do you think?

matthewturk · 2020-02-06T13:47:16Z

I suspect it will be fine to include it now. I will work on getting answer tests in #2437 working for the RAMSES datasets before merging this into that PR.

brittonsmith · 2020-02-06T14:17:12Z

I think as long as answer tests are run at some point, it's ok with me. If that happens as part of the yt-4.0 PR, that's fine, too.

yt/fields/fluid_fields.py

yt/frontends/ramses/data_structures.py

cphyc · 2020-02-12T16:21:46Z

Since @brittonsmith approved the PR, I refactored parts of it to reduce code repetition and make it more readable. Sorry for updating a PR post-review!

brittonsmith

New changes look ok to me!

cphyc · 2020-05-05T12:39:56Z

For information, I just rebased the branch on yt-4.0. @matthewturk, you mentioned #2437 (in #2425 (comment)), is it still a blocker to merge this PR?

munkm · 2020-05-28T19:23:31Z

We're temporarily blocking all merges to the yt-4.0 branch until #2437 gets merged. Is that what you're asking?

cphyc · 2020-05-29T08:47:26Z

Yes, thanks for the reply @munkm ;) In the meantime is there anything I can do to prepare this PR for a future merge into yt-4?

munkm · 2020-05-29T13:40:34Z

I'll do my best to review it today! I can't think of anything at the moment. 🙂

cphyc · 2020-06-11T15:01:04Z

yt/geometry/oct_visitors.pyx

+        for i in range(-1, 3):
+            ishift[0] = i
+            for j in range(-1, 3):
+                ishift[1] = j
+                for k in range(-1, 3):
+                    ishift[2] = k
+                    self.set_neighbour_info(o, ishift)
+


This piece of code is very much optimizable. Here are possible routes to achieve much better performance:

Work oct by oct, which would reduce the number of neighbor lookup from 4³=64 to 3³=27,

Use faster neighbor lookup method(s). For now, all searches are started from the root mesh down to leaf nodes, but we could instead go up the tree from the central oct then down to find all neighbors (see e.g. https://geidav.wordpress.com/2017/12/02/advanced-octrees-4-finding-neighbor-nodes/).

Pre-compute the face-neighbors of all octs.

Note that for the last point, algorithms exist that generate the neighbors of all octs in O(1) time (https://link.springer.com/article/10.1007/s13319-015-0060-9) during the octree construction.
Another possible solution would be to keep a unordered hash map of all the octs indexed by their (3-integers) position. With such structure, finding a neighbor takes O(1) time. This could even come as a replacement of the current pointer-based octree structure.

I think this should be in a comment. It obviously doesn't need to be implemented right now.

cphyc · 2020-06-30T19:58:46Z

@yt-fido test this please

cphyc · 2020-07-01T09:16:17Z

So the new tests answer tests are hidden here https://tests.yt-project.org/job/yt_py38_git/423/testReport/yt.frontends.ramses.tests (they pass).

Xarthisius · 2020-07-01T14:36:13Z

So the new tests answer tests are hidden here https://tests.yt-project.org/job/yt_py38_git/423/testReport/yt.frontends.ramses.tests (they pass).

The test that you've added won't work as intended. There's a few things that need to be fulfilled for a method to work as an answer test:

It has to be decorated with requires_ds (not requires_file)
It has to yield test classes based on AnswerTestingTest e.g. FieldValuesTest
It has to be explicitly enumerated in tests/tests.yaml as yt/frontends/ramses/tests/test_outputs.py:test_field_accession
In case of name collision (which would happen here with test_output_00800 should you yield some additional FieldValuesTest) a unique suffix needs to be added to the test name.

Ultimately it would look roughly like this:

@requires_ds(output_00080)
def test_field_accession():
    ds = data_dir_load(output_00080)
    fields = [ 
        ('gas', 'density'),  # basic ones
        ('gas' ,'pressure'),
        ('gas', 'pressure_gradient_magnitude'), # requires ghost zones
    )   
    dso = ( 
        None,
        ("sphere", ([.1]*3, (0.01, 'unitary'))),
        ("sphere", ([.5]*3, (0.05, 'unitary'))),
        ("box", ([.1]*3, [.2]*3)),
    )   
    for dobj_name in dso:
        for field in fields:
            test = FieldValuesTest(ds, field, dobj_name)
            test.suffix = "accession"
            test_field_accession.__name__ = test.description + test.suffix
            yield test

Having said all that, I think most of the things from this tests are already done in the test_output_00800. It tests density and pressure gradient magnitude, it tests entire domain and one sphere. If you would add ('gas' ,'pressure') to L20, and additional objects in L28 it would do the same (and more) for you. For the sake of efficiency though I'd evaluate if test_output_00800 is sufficient to test your work even without any additions.

neutrinoceros

I'm not competent enough to approve this, but I left a few comments and questions, I hope it helps making this even better. In any case there's a lot going on here and you're doing great ! Keep up the awesome PRs

yt/data_objects/data_containers.py

yt/fields/fluid_fields.py

yt/fields/geometric_fields.py

yt/geometry/oct_container.pyx

yt/frontends/ramses/tests/test_outputs.py

cphyc · 2020-07-01T18:16:42Z

So the new tests answer tests are hidden here https://tests.yt-project.org/job/yt_py38_git/423/testReport/yt.frontends.ramses.tests (they pass).

The test that you've added won't work as intended. There's a few things that need to be fulfilled for a method to work as an answer test:
1. It has to be decorated with `requires_ds` (not `requires_file`)

2. It has to yield test classes based on `AnswerTestingTest` e.g. `FieldValuesTest`

3. It has to be explicitly enumerated in `tests/tests.yaml` as `yt/frontends/ramses/tests/test_outputs.py:test_field_accession`

4. In case of name collision (which would happen here with `test_output_00800` should you yield some additional FieldValuesTest) a unique suffix needs to be added to the test name.
Ultimately it would look roughly like this:
@requires_ds(output_00080)
def test_field_accession():
    ds = data_dir_load(output_00080)
    fields = [ 
        ('gas', 'density'),  # basic ones
        ('gas' ,'pressure'),
        ('gas', 'pressure_gradient_magnitude'), # requires ghost zones
    )   
    dso = ( 
        None,
        ("sphere", ([.1]*3, (0.01, 'unitary'))),
        ("sphere", ([.5]*3, (0.05, 'unitary'))),
        ("box", ([.1]*3, [.2]*3)),
    )   
    for dobj_name in dso:
        for field in fields:
            test = FieldValuesTest(ds, field, dobj_name)
            test.suffix = "accession"
            test_field_accession.__name__ = test.description + test.suffix
            yield test
Having said all that, I think most of the things from this tests are already done in the test_output_00800. It tests density and pressure gradient magnitude, it tests entire domain and one sphere. If you would add ('gas' ,'pressure') to L20, and additional objects in L28 it would do the same (and more) for you. For the sake of efficiency though I'd evaluate if test_output_00800 is sufficient to test your work even without any additions.

Thanks for the detailed answer. The answer tests should now be included in the already-existing test_output_00080 answer tests. I also added some other unit tests, tell me if you're ok with them.

neutrinoceros

I changed my mind, here's my approval 🥇

yt/data_objects/octree_subset.py

cphyc · 2020-07-13T21:33:30Z

@yt-fido test this please.

matthewturk

I think this is really good. Only a few comments.

matthewturk · 2020-07-16T20:02:40Z

yt/data_objects/data_containers.py

-                            self.selector,
-                            gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz],
-                            rv, ind)
+                    data = gz[field][ngz:-ngz, ngz:-ngz, ngz:-ngz]


I like getting rid of the conditional! Nice.

yt/frontends/ramses/data_structures.py

yt/geometry/oct_container.pyx

matthewturk · 2020-07-16T20:10:51Z

yt/geometry/oct_visitors.pxd

+cdef class BaseNeighbourVisitor(OctVisitor):
+    cdef int idim      # 0,1,2 for x,y,z
+    cdef int direction # +1 for +x, -1 for -x
+    cdef np.int8_t[:] neigh_ind


two quick things -- should neigh_ind be a 3d view? and, it's a bit faster if we declare ::1 in them so that it knows they have stride 1.

neigh_ind has size 3. I actually just changed it to use a C array instead of a memory view.

matthewturk · 2020-07-16T20:10:57Z

yt/geometry/oct_visitors.pxd

+    cdef np.int64_t[:] domain_inds
+
+cdef class NeighbourCellVisitor(BaseNeighbourVisitor):
+    cdef np.uint8_t[:] levels


same here about stride

matthewturk · 2020-07-16T20:11:54Z

yt/geometry/oct_visitors.pyx

+            # Assuming periodicity
+            if fcoords[i] < 0:
+                fcoords[i] += 1
+            elif fcoords[i] > 1:


can you make this explicitly a float? and, shouldn't this be checking if it's greater than (and above, less than) the domain width?

actually i may be wrong on that, but want confirmation.

Thing is, we are computing here the center of each cell, which cannot be at the edge. For example, the cell positions at level 2 in 1D would be [0.125, 0.375, 0.625, 0.875] to which we add ±0.25, so that should cause no problems.

However if you prefer, I can do

Suggested change

elif fcoords[i] > 1:

elif fcoords[i] >= 1:

matthewturk · 2020-07-16T20:13:13Z

yt/geometry/oct_visitors.pyx

+        for i in range(-1, 3):
+            ishift[0] = i
+            for j in range(-1, 3):
+                ishift[1] = j
+                for k in range(-1, 3):
+                    ishift[2] = k
+                    self.set_neighbour_info(o, ishift)
+


I think this should be in a comment. It obviously doesn't need to be implemented right now.

cphyc

Comments addressed! Thanks for them

yt/frontends/ramses/data_structures.py

cphyc · 2020-07-17T08:30:38Z

yt/frontends/ramses/data_structures.py

+
+        if num_ghost_zones > 0:
+            if not all(ds.periodicity):
+                warnings.warn('Ghost zones will wrongly assume the domain to be periodic.')


No, I changed it to raise a RuntimeError unless there's a better one to raise here?

cphyc · 2020-07-17T08:31:00Z

yt/frontends/ramses/data_structures.py

        fill_hydro(fd, file_handler.offset,
-                   file_handler.level_count, levels, cell_inds,
+                   file_handler.level_count,
+                   [self.domain_id-1],


I made the name of the variable explicit here and below (the next fill_hydro occurrence).

yt/frontends/ramses/data_structures.py

cphyc · 2020-07-17T08:35:05Z

yt/frontends/ramses/data_structures.py

+        oinds = oinds.reshape(-1, 64)
+        cinds = cinds.reshape(-1, 64)


Correct, I changed the names.

cphyc · 2020-07-17T08:47:57Z

yt/geometry/oct_container.pyx

            for i in range(levels.shape[0]):
                if levels[i] != level: continue


cphyc · 2020-07-17T08:48:47Z

yt/geometry/oct_container.pyx

+        -------
+        oct_inds : int64 ndarray (nocts*8, )
+            The on-domain index of the octs containing each cell
+        cell_inds : uint8 array (nocts*8, )


cphyc · 2020-07-17T08:51:24Z

yt/geometry/oct_visitors.pxd

+cdef class BaseNeighbourVisitor(OctVisitor):
+    cdef int idim      # 0,1,2 for x,y,z
+    cdef int direction # +1 for +x, -1 for -x
+    cdef np.int8_t[:] neigh_ind


neigh_ind has size 3. I actually just changed it to use a C array instead of a memory view.

cphyc · 2020-07-17T08:57:38Z

yt/geometry/oct_visitors.pyx

+            # Assuming periodicity
+            if fcoords[i] < 0:
+                fcoords[i] += 1
+            elif fcoords[i] > 1:


Thing is, we are computing here the center of each cell, which cannot be at the edge. For example, the cell positions at level 2 in 1D would be [0.125, 0.375, 0.625, 0.875] to which we add ±0.25, so that should cause no problems.

However if you prefer, I can do

Suggested change

elif fcoords[i] > 1:

elif fcoords[i] >= 1:

cphyc · 2020-07-17T09:00:00Z

yt/geometry/oct_visitors.pyx

+        for i in range(-1, 3):
+            ishift[0] = i
+            for j in range(-1, 3):
+                ishift[1] = j
+                for k in range(-1, 3):
+                    ishift[2] = k
+                    self.set_neighbour_info(o, ishift)
+


cphyc commented Jan 30, 2020

View reviewed changes

yt/fields/fluid_fields.py Outdated Show resolved Hide resolved

matthewturk reviewed Jan 30, 2020

View reviewed changes

yt/data_objects/octree_subset.py Outdated Show resolved Hide resolved

yt/fields/fluid_fields.py Outdated Show resolved Hide resolved

matthewturk reviewed Jan 30, 2020

View reviewed changes

setup.py Outdated Show resolved Hide resolved

yt/fields/fluid_fields.py Outdated Show resolved Hide resolved

yt/geometry/oct_container.pyx Show resolved Hide resolved

cphyc changed the title ~~[WIP] Octree ghost zones — reloaded~~ Octree ghost zones — reloaded Feb 5, 2020

brittonsmith approved these changes Feb 6, 2020

View reviewed changes

Xarthisius reviewed Feb 6, 2020

View reviewed changes

yt/fields/fluid_fields.py Show resolved Hide resolved

yt/frontends/ramses/data_structures.py Outdated Show resolved Hide resolved

cphyc mentioned this pull request Feb 10, 2020

[WIP] Octree ray tracing #2442

Closed

5 tasks

brittonsmith approved these changes Feb 13, 2020

View reviewed changes

munkm mentioned this pull request Apr 3, 2020

Ghost zones should only fill boundary #363

Open

munkm added index: octree new feature Something fun and new! labels Apr 24, 2020

cphyc force-pushed the neighbour-search branch from f5c2464 to fd3f5ed Compare May 5, 2020 10:45

cphyc mentioned this pull request May 27, 2020

Volume rendering for octrees #2610

Merged

8 tasks

cphyc mentioned this pull request Jun 4, 2020

Support vertex centred data #2615

Merged

cphyc commented Jun 11, 2020

View reviewed changes

cphyc added the yt-4.0 feature targeted for the yt-4.0 release label Jun 11, 2020

cphyc mentioned this pull request Jun 12, 2020

Octree dataset don't support generating ghost zones #1051

Open

cphyc added 2 commits June 16, 2020 19:55

Add support for fast neighbour search

333436c

More generic code

6ce0416

neutrinoceros reviewed Jul 1, 2020

View reviewed changes

cphyc added 10 commits July 2, 2020 10:38

Add kwa to select calls

6c318dc

ds → dx

236134b

Return "arr" not "tmp" variable

550e887

Each line for its argument

0fc5131

Minor nits in docstrings

f7d9e4a

Using python-friendly syntax

9ca0d0e

Remove duplicates in testing file

a92d477

Merge

58f4384

_block_reorder -> _block_order

0bf76b2

Field detector does not have a block_order argument

19058d8

neutrinoceros approved these changes Jul 2, 2020

View reviewed changes

cphyc commented Jul 2, 2020

View reviewed changes

yt/data_objects/octree_subset.py Outdated Show resolved Hide resolved

Reorder if necessary

6d2b439

Merge branch 'master' into neighbour-search

f21e058

matthewturk approved these changes Jul 16, 2020

View reviewed changes

Take comments into account

3c8b75c

cphyc commented Jul 17, 2020

View reviewed changes

cphyc added 5 commits July 17, 2020 10:06

Inform Cython about strides

13e5dfa

Merge

d97eb12

tr -> data

b199540

Make isort happy

8ecce95

Fix flaking

c033f18

munkm merged commit 6daf9cd into yt-project:master Jul 19, 2020

cphyc deleted the neighbour-search branch July 20, 2020 09:51

		for i in range(levels.shape[0]):
		if levels[i] != level: continue

Octree ghost zones — reloaded #2425

Octree ghost zones — reloaded #2425

Conversation

cphyc commented Jan 23, 2020

PR Summary

The good

The bad

matthewturk commented Jan 23, 2020

cphyc commented Jan 30, 2020

The good

The bad

matthewturk left a comment

Choose a reason for hiding this comment

cphyc commented Feb 5, 2020

cphyc commented Feb 5, 2020

brittonsmith left a comment

Choose a reason for hiding this comment

matthewturk commented Feb 6, 2020

brittonsmith commented Feb 6, 2020

cphyc commented Feb 12, 2020

brittonsmith left a comment

Choose a reason for hiding this comment

cphyc commented May 5, 2020

munkm commented May 28, 2020

cphyc commented May 29, 2020 • edited Loading

munkm commented May 29, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cphyc commented Jun 30, 2020

cphyc commented Jul 1, 2020

Xarthisius commented Jul 1, 2020 • edited Loading

neutrinoceros left a comment

Choose a reason for hiding this comment

cphyc commented Jul 1, 2020

neutrinoceros left a comment

Choose a reason for hiding this comment

cphyc commented Jul 13, 2020

matthewturk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cphyc left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cphyc commented May 29, 2020 •

edited

Loading

Xarthisius commented Jul 1, 2020 •

edited

Loading