Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Optimized radius calculation. #4079

Merged
merged 4 commits into from
Aug 16, 2022
Merged

ENH: Optimized radius calculation. #4079

Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e6d15b3
Optimized get_radius function.
yipihey Aug 14, 2022
699c2ad
Clean up comments
yipihey Aug 14, 2022
3aaeac6
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Aug 14, 2022
81b0ac8
ENH: save an unused array copy
neutrinoceros Aug 14, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 14 additions & 17 deletions yt/fields/field_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,37 +7,34 @@


def get_radius(data, field_prefix, ftype):
unit_system = data.ds.unit_system
center = data.get_field_parameter("center").in_base(unit_system.name)
DW = (data.ds.domain_right_edge - data.ds.domain_left_edge).in_base(
unit_system.name
)
# This is in cm**2 so it can be the destination for our r later.
center = data.get_field_parameter("center").to("code_length")
DW = (data.ds.domain_right_edge - data.ds.domain_left_edge).to("code_length")
# This is in code_length so it can be the destination for our r later.
radius2 = data.ds.arr(
np.zeros(data[ftype, field_prefix + "x"].shape, dtype="float64"), "cm**2"
np.zeros(data[ftype, field_prefix + "x"].shape, dtype="float64"), "code_length"
)
r = radius2.copy()

r = radius2.v
neutrinoceros marked this conversation as resolved.
Show resolved Hide resolved
if any(data.ds.periodicity):
rdw = radius2.copy()
rdw = radius2.v
for i, ax in enumerate("xyz"):
# This will coerce the units, so we don't need to worry that we copied
# it from a cm**2 array.
np.subtract(
data[ftype, f"{field_prefix}{ax}"].in_base(unit_system.name),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The removal of the unit conversion (i.e., in_base(unit_system.name)) has changed radius values for frontends where the default unit is not "code_length". I'm not sure how much performance is lost by adding it back, but I think it is necessary unless we mandate position units always come back in "code_length", which I don't think we do.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes perhaps we are mixing two concepts.
One is "code_units" and one is the "original_units", i.e. which units are used in the file we are using. I think it is advantageous to require frontends to not modify the numbers it reads from disk.
I think of two different use cases of yt.

  1. Debugging a code that generates the data we are reading in. It is desirable that yt attempts to not modify the data it reads so we can have some certainty we are debugging the code that generates the data.
  2. If I use yt to convert from one data format to another. E.g. I read all particle positions in enzo and then store them as a simpler hdf5 file where the positions are a simple array. In this case I want to make sure that I'm just changing the shape of the data but do not change any of the numbers.

The positive side effect that it also avoid unnecessary computation which saves a little time.

Is there already the concept of "original_units" ? I.e. give us a simple way to check in what units the data was provided?

center[i],
data[ftype, f"{field_prefix}{ax}"].d,
center[i].d,
r,
)
if data.ds.periodicity[i]:
np.abs(r, r)
np.subtract(r, DW[i], rdw)
np.subtract(r, DW.d[i], rdw)
np.abs(rdw, rdw)
np.minimum(r, rdw, r)
np.multiply(r, r, r)
np.add(radius2, r, radius2)
np.add(radius2.d, r, radius2.d)
if data.ds.dimensionality < i + 1:
break
# Now it's cm.
np.sqrt(radius2, radius2)
# Using the views into the array is not changing units and as such keeps
# from having to do symbolic manipulations
np.sqrt(radius2.d, radius2.d)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like this.

# Alias it, just for clarity.
radius = radius2
return radius
Expand Down