Subdomain extraction

[Adnan-Ali-Ahmad]

This PR is linked to #52 and is in anticipation of an upcoming PR that will handle cartopy projections of surfaces (#53).

A spatial.py file containing common coordinate transforms alongside matrix rotations has been added. I am still open for discussion on how these transformations should be handled/called, but for the moment the current routines are generic enough for us to easily modify in the future.

Git should handle the diff if #81 and #85 are merged.

[nvaytet]

I appreciate you wanting to make a start on this earlier rather than later, but as I've said before I think it's worth taking our time to think properly about how we want to implement this.

We have to ask ourselves what is the goal here.
If we want to be able to e.g. easily compute a radius, and make a radial profile, then these functions would definitely be useful.

If we wish to make maps that depend on theta and phi (which was your original motivation), we also need to consider how this would then be used to make the maps. Maybe this approach does work, I am not sure, I have not thought enough about it.

Another approach could have been to give the Array additional properties, such as .r, .theta, and .phi, which would compute the radius and angles on-the-fly. I am thinking that the coordinates are inherently cartesian, because they are as such inside Ramses. So they should be loaded as such, and am therefore not sure that converting from spherical back to cartesian would be something that is ever needed?

In addition, as you pointed out when you first suggested this idea, velocity or B field vectors could also have radial, theta and phi components (and those calculations would be a little different from that of the positions), which would certainly be useful.

This last point, combined with the slicing inconsistency of vectors in the other PR, is making me wonder whether we want a slightly different class for vector arrays (that would inherit from Array).

[Adnan-Ali-Ahmad]

Another approach could have been to give the Array additional properties, such as .r, .theta, and .phi, which would compute the radius and angles on-the-fly.

In addition, as you pointed out when you first suggested this idea, velocity or B field vectors could also have radial, theta and phi components (and those calculations would be a little different from that of the positions), which would certainly be useful.

Maybe we can group these into .get_spherical_components() & .get_cylindrical_components() methods? These shouldn't be automatically computed by osyris upon loading the dataset because otherwise it comes with an extra RAM burden, but they can provide a quick and easy way to gather all components.

This last point, combined with the slicing inconsistency of vectors in the other PR, is making me wonder whether we want a slightly different class for vector arrays (that would inherit from Array).

Instead of having a different class for vectors, I would rather have the type be an attribute of the Array class (eg. a .vector boolean attribute).

[nvaytet]

In your original suggestion, you also had vector bases as part of the arguments. I think these would still be needed here?

[Adnan-Ali-Ahmad]

Indeed, I just wanted to keep it simple and generic for now, but expressing your coordinates in a new basis requires you to simply rotate them (which is the purpose of the rotation_matrix function).

[nvaytet]

Maybe thinking about how one would use this in a notebook would be a good way to figure out the API.

You have a disk somewhere in your simulation. You want to basically rotate everything so that the z axis is aligned with the angular momentum (is that all you need to do?)

You could have the following steps:

• find angular momentum vector
• rotate all coords and vectors into a new basis (still in cartesian coords)
• then access the radial or angular components of your vectors (this could be done via the .r and .theta properties?

The problem with just having .r property on a vector is that it works for position vectors, but not for other vectors such as velocity. If you want the radial velocity, you need to take the scalar product of the velocity vector with the position vector, so the conversion needs both position vectors and velocity vectors as inputs.

Maybe that is an argument for having functions that turn a whole Datagroup into spherical or cylindrical components?
Something like

spherical_hydro = osyris.spatial.get_spherical_components(data['hydro'])

Which would return a new Datagroup (i think it might be safer to return a new group rather than doing operations in-place).

But then the get_spherical_components also needs the data['amr']['position'] array...
Not sure if sending the entire Dataset to be converted is good, because you may be doubling your RAM usage.

Maybe the function should just accept a single Array, and an optional positions argument, which could just be None for computing the spherical components of the cell positions?

[Adnan-Ali-Ahmad]

You have a disk somewhere in your simulation. You want to basically rotate everything so that the z axis is aligned with the angular momentum (is that all you need to do?)

Yes, that's all there is to it.

The problem with just having .r property on a vector is that it works for position vectors, but not for other vectors such as velocity. If you want the radial velocity, you need to take the scalar product of the velocity vector with the position vector, so the conversion needs both position vectors and velocity vectors as inputs.

Yeah that's what's been causing me a bit of a headache. It's difficult to come up with a clean API for something like this.

Not sure if sending the entire Dataset to be converted is good, because you may be doubling your RAM usage.

That is why in #52 I suggested that we make the coordinate system an attribute of the Dataset. That way with data.transform(system="spherical", basis="top"), you do not return a new dataset, but simply convert all components in it to spherical and keep the old basis stored somewhere so that the user can revert his changes. However, this solution is going to be quite costly in development time.

What I suggest is to have a get_spherical_components function that takes a Dataset as the argument but also a dictionary containing all the groups and arrays that the user would like to convert. Something like:

osyris.spatial.get_spherical_components(data, {"amr":"position", "hydro": ["velocity", "B_field"]})

It should be easier in this case to make a few checks on the user's request (eg. if they request spherical velocity components, we can access the position array since the whole dataset has been passed as an argument).

[Adnan-Ali-Ahmad]

I forgot to say that this suggestion would return a new Dataset containing all the converted coordinates and .r, .theta and .phi components for each array.

Eg. :

spherical_coords = osyris.spatial.get_spherical_components(data, {"amr":"position", "hydro": ["velocity", "B_field"]})

spherical_coords["amr"]["position"].r # radius
spherical_coords["amr"]["position"].theta # colatitude
spherical_coords["amr"]["position"].phi # azimuth

spherical_coords[“hydro"]["velocity"].r # radial velocity
spherical_coords[“hydro"]["velocity"].phi # azimuthal velocity
etc...

And spherical_coords["amr"]["position"].x returns an attribute error.

[nvaytet]

Trying to get this going again:

After doing some more thinking about this, about how it would be used, I think what you would want is to be able to

• rotate your entire domain into a new basis (could be to align to angular momentum)
• some mechanism to obtain spherical or cylindrical components of vectors (position or other)

That's probably what it boils down to (I think?).

I suggest the following:

1. Make a function that would rotate your domain, supplying a normal vector (we could use the top and side orientations like for the maps, or make the function that computes the angular momentum vector more accessible, or both)
2. Make .r, .theta, .phi properties on the vectors that would each compute and return a new Array (you can always store those in your Datagroup if you don't want to compute them over and over again)
3. If the .theta, .phi is requested on a non-position vector, get the corresponding position vector from the .parent reference to compute the components (for hydro velocities, we would have to go a look into the .parent.parent['amr'], but we are already doing this for plotting maps)
4. If needed, make wrapper functions that can convert a whole vector to spherical or cylindrical (which calls the .r, .theta, .phi internally)

We could try this as a first step, see how usable it is, and make wrapper functions that convert an entire DataGroup in a later step, if we feel that is also needed?

[Adnan-Ali-Ahmad]

PR is indeed trying to do two things, but this keeps things simpler for us to review the code together.

Checklist:

• Subdomain extraction using spatial criterion (cubic & spherical)
• Numpy indexing for subdomain extraction
• Spherical coordinate transforms
• Wrappers in Vector for spherical coordinate transforms
• Cylindrical coordinate transforms
• Change of origins
• Change of basis using rotations
• Add matmul operator overload for Vector to make code cleaner and more user friendly (Matmul operator overload for osyris Arrays #85)
• Spatial toolkit to be accessed by users with osyris.spatial
• Add examples to doc
• Add coordinate transform unit tests to test suite

[Adnan-Ali-Ahmad]

Current workflow looks like this:

# data load
data = osyris.Dataset(8,path=path).load()

ind = np.argmax(data["hydro"]["density"])
center = data["amr"]["position"][ind.values] # (or data["sink"]["position"])

#proceed to subdomain extraction
subdomain = osyris.Subdomain(data, selection="spherical", dr=1e3*osyris.units("au"), origin=center, basis="top", dr_L=5e2*osyris.units("au"))

#do a histogram
fig,ax = plt.subplots()
osyris.histogram2d(subdomain["amr"]["position"].r.to("au"), subdomain["hydro"]["velocity"].r.to("km/s"),
                   norm="log", logx=True, logy=False, cmap="viridis", ax = ax)
fig.show()

Basis can be changed with spatial.change_basis, and the origin with spatial.change_origin.
Accessing spherical components is fast and easy in subdomains, with the .r, .theta and .phi components using a rotation matrix, allowing for any arbitrary cartesian vector to be converted into its spherical counterpart.

[Adnan-Ali-Ahmad]

This is however messing up in slices. I never had this problem in my own simulations, presumably because they're much higher resolution, but I'm getting some weird effects as can be seen here:

osyris.map({"data": subdomain["hydro"]["density"], "norm":"log"},
                     {"data": subdomain["sink"]["position"], "mode": "scatter", "c": "white"},
             dx=2e3 * osyris.units("au"),
             direction="z",plot=True, resolution=256)

This doesn't seem to occur when the chosen basis for the subdomain is `[0,1,0]` or `[1,0,0]`:

I'm still trying to track the issue down, but I'm certain from the spherical components that things are working as intended. Maybe there's something I've misunderstood with slices.

[nvaytet]

This is however messing up in slices. I never had this problem in my own simulations, presumably because they're much higher resolution, but I'm getting some weird effects as can be seen here:

The effect on maps is indeed weird. You can actually see the little squares from the cells that are rotated (around the edges).
I suspect it comes from the fact that when making a map, we try to find for each image pixel the cell that contains it. And there must be somewhere in there an assumption on the orientation of the grid.

[nvaytet]

Have a look at the changes in #102 , maybe if we get that merged it might make some things easier/more compact here?

[Adnan-Ali-Ahmad]

Have a look at the changes in #102 , maybe if we get that merged it might make some things easier/more compact here?

It should reduce the number of .values calls.

[Adnan-Ali-Ahmad]

Thanks for all your comments. I'll work on this again and push the changes sometime this week.

[Adnan-Ali-Ahmad]

PR now simply extracts subdomains. Ordinary coordinate transforms will be added after this PR is merged.

[Adnan-Ali-Ahmad]

@nvaytet I think this PR is ready to be merged, and we can move on to discuss the changes in #106 ?