WIP: Uniform Cylindrical and Spherical Coordinates

[forrestglines]

Very WIP PR to support uniform cylindrical and spherical coordinates in Parthenon.

PR Summary

Many potential AthenaPK users need treatments of spherical and cylindrical coordinates that can preserve angular momentum. Most of the machinery for these coordinates should belong in Parthenon. This PR is an update of Shengtai's work to match changes to Parthenon's coordinate renaming.

PR Checklist

• Code passes cpplint
• New features are documented.
• Adds a test for any bugs fixed. Adds tests for new features.
• Code is formatted
• Changes are summarized in CHANGELOG.md
• CI has been triggered on Darwin for performance regression tests.
• Docs build
• (@lanl.gov employees) Update copyright on changed files

[forrestglines]

Outstanding Questions and To-do's:

1. Usage of x1s2, x1s3->Xs<dir,side>, etc. These are "area-averaged" terms but documentation in Athena++ is not specific what coordinate and over what area the coordinate is averaged. My best guess is that x1s2 is the "1" coordinate averaged over side "2". These terms are used for prolongation of field-centered variables and so might not be needed until a CT MHD method exists in AthenaPK.

2. Naming of h2v, h2f, dh2vd1, dh32vd2, etc. What are these terms? In Athena++, they seem to only be used for viscosity and conduction. Ideally we'd rename these with something more descriptive and replace the different functions with a handful of functions with template parameters like we've done with X<dir> and Xf<dir,face>. However, since their use is limited, we could spin these off into another friend object of the coordinate classes.

3. Do geometric source terms belong within Parthenon? Can these terms be abstracted or are they specific for AthenaPK's implementation?

4. How should coordinate terms be cached to save computation? Both the cylindrical and spherical coordinates in Athena++ make use of cached areas for coordinate systems. Some of these cached arrays may only be used in certain contexts (i.e. just for viscosity). Additionally, the current implementation of cached coordinate data in this PR is likely broken. For now, we could replace the cached data with functions so that for initial testing everything is computed on the fly.

5. Is there a simple test of coordinate systems that we can add to Parthenon? Maybe a disk advection problem?

[Yurlungur]

A few comments below. But a bigger picture comment also: would you be willing to chat with myself, Josh, and Luke about what's actually needed here? We were chatting today and we're pretty sure that (at least in a finite volume context) this interface can be simplified. Two things in particular:

1. The geometric source terms can be very likely expressed as differences in "areas" or "generalized volumes"
2. For prolongation, there are two ways to go. The way we copied from Athena++ uses differences of centroid positions, which requires storing centroids. However, you can also "densitize" the prolongated quantities before and after interpolation. For linear reconstructions, this is equivalent---there aren't that many ways do something physically correct with only 2 degrees of freedom. For higher-order reconstructions, the densitized approach is much cleaner.
3. Restriction doesn't need to be changed.

[Yurlungur]

if coordinate systems are going to compile-time only, maybe we drop coord from the input deck entirely?

[Yurlungur]

Hmm... is this actually something we need/want the Parthenon coords tooling to provide? My recommendation is one of the following:

1. All coords provide functionality to convert between themselves and Cartesian. For UniformCartesian, this is an identity transformation
2. Just make this not parthenon's problem.

Probably but not necessarily we want to go with option 1, as that's what's needed for vis tooling and initial conditions in some contexts.

[Yurlungur]

longer discussion below. Resolving this comment

[forrestglines]

From what I've done and what Shengtai included in his demo, there's a bigger use case to provide spherical and cylindrical conversions in Parthenon since we're likely to deal with physical systems where those coordinates make sense (i.e. galaxy clusters and stars for spherical setups, jets and disks for cylindrical setups). We could leave it up to each of these systems to implement these conversions on their own but I see them repeated so many times elsewhere it might be better to centralize that code in one place.

I'm not set on including these conversions in Parthenon at the moment but still thinking about it. We might also include conversions for vectors.

[Yurlungur]

I see. No I think that makes sense. As I comment below, we may also want it for vis tools. It would be better if, for example, Paraview could automatically visualize simulations done in spherical coordinates. So we might want to standardize some functions that always map from whatever coordinate system we're in to Cartesian.

[Yurlungur]

I have not had good success making the ternary operator vectorize. Maybe we either switch this to using masks or we split it into 3 functions?

[Yurlungur]

Also this appears to be for cell centers only?

[forrestglines]

This implementation was lifted from an inline function for swarms, I copied it nearly verbatim. I should change the function description for "cell centered index". I think this is sufficient for

My very outdated experience on Kepler GPUs was that ternary operators had less branch splitting, so I also default to using ternary operators.

[Yurlungur]

This implementation was lifted from an inline function for swarms, I copied it nearly verbatim. I should change the function description for "cell centered index". I think this is sufficient for

👍

My very outdated experience on Kepler GPUs was that ternary operators had less branch splitting, so I also default to using ternary operators.

It's not a huge deal either way, but I was thinking of CPUs when I made the comment, not GPUs. My experience is that GPUs are much better at handling ternary/branching than CPUs.

[Yurlungur]

not sure I like the name here.

[Yurlungur]

should these be if constexprs?

[Yurlungur]

These smell to me like coefficients in the christoffel symbols or something similar. I don't think we need them all for finite volumes-type discretizations... especially if we're restricting ourselves to cartesian, spherical, and cylindrical coordinates with uniform coordinate spacing.

[Yurlungur]

these are all 1D right? You could use ParArray1D here probably.

That said, as a first pass, these calculations are only being cached for performance. We could do them without and just compute inline as well. The only reason not to do this is that the trig functions are not transcendental.

[Yurlungur]

@forrestglines my apologies---After a longer discussion I think @lroberts36 @jdolence and I have all convinced ourselves that my above comments are not right:

The geometric source terms can be very likely expressed as differences in "areas" or "generalized volumes"

This is still mostly true, I think, but that doesn't mean all geometry output can be dropped. In particular, the face area argument holds for equations similar to the hydro equations, but not for all conservation laws. Equations with a diffusive, e.g., second derivative, term don't work out as nicely.

For prolongation, there are two ways to go. The way we copied from Athena++ uses differences of centroid positions, which requires storing centroids. However, you can also "densitize" the prolongated quantities before and after interpolation. For linear reconstructions, this is equivalent---there aren't that many ways do something physically correct with only 2 degrees of freedom. For higher-order reconstructions, the densitized approach is much cleaner.

I worked out the densitized interpolation vs centroids... I think there's reason to provide centroids---they are not exactly the same as the densitized thing, even in the linear case. My bad on this.

Regarding your questions/TODOs:

Usage of x1s2, x1s3->Xs<dir,side>, etc. These are "area-averaged" terms but documentation in Athena++ is not specific what coordinate and over what area the coordinate is averaged. My best guess is that x1s2 is the "1" coordinate averaged over side "2". These terms are used for prolongation of field-centered variables and so might not be needed until a CT MHD method exists in AthenaPK.

I think in the "topological element" framework that @lroberts36 developed, these become centroid values based on different topological elements, e.g., Face1 or Face2. We probably want to rethink the naming convention here... s for "centroid" and c for cell-center is in my opinion confusing. I would rather switch to X for coordinate location (templated on topological element) and Xc for centroid (templated on topological element).

Naming of h2v, h2f, dh2vd1, dh32vd2, etc. What are these terms? In Athena++, they seem to only be used for viscosity and conduction. Ideally we'd rename these with something more descriptive and replace the different functions with a handful of functions with template parameters like we've done with X

and Xf<dir,face>. However, since their use is limited, we could spin these off into another friend object of the coordinate classes.

I looked a bit more carefully at this... In conduction and viscosity, h*v shows up in calculations of derivatives in combinations with grid spacings that indicate they're used for physical lengths. I.e., the line element in an integral so you can compute the length in a derivative.

I believe these are also expressed more descriptively in other parts of the coordinate machinery. For example, the term

(pco_->h2v(i)*pco_->dx2v(j-1))

that shows up in the x2 flux of conduction in Athena++ can be re-expressed as VolCenter2Length in the Athena++ coordinates class. I wonder if there's a lack of code cleanup here where essentially identical code paths were not unified. It might be worth asking the Athena++ folks.

The dh*vd* versions are the derivatives of these lengths with respect to another direction. And v vs f is for volume (i.e., cell center) vs., cell face. This also shows up in equations with second derivatives, but in particular, it shows up when differentiating vector quantities in those second derivative terms. These are coming from the parallel transport of vectors in curved space. I.e., they represent that a radial vector does not point in the same direction at $\phi = 0$ and $\phi=\pi/2$.

There is a mathematical formalism for describing these relevant quantities. In particular the h* terms are formally described as part of the metric tensor, which describes the physical magnitudes of vectors in curved space and, as a corollary, describes lengths of paths, which are the integrals of the magnitudes of tangent vectors. The dh*vd* type terms are formally described by the Christoffel symbols these are essentially rotation coefficients describing how a vector rotates as you move in a coordinate direction (for example $\theta$).

In full generality the metric is a symmetric 3x3 matrix, so it has six possible nonzero terms. In Cartesian, spherical, and cylindrical coordinates, it is a diagonal matrix, so it's 3 nonzero numbers. In full generality, the Christoffel symbols are three symmetric matrices, so 18 nonzero numbers. However for the Cartesian, spherical, and cylindrical coordinates, there are at most two nonzero terms.

Do geometric source terms belong within Parthenon? Can these terms be abstracted or are they specific for AthenaPK's implementation?

I could go either way on this. I think one option would be to use the above geometric description, using the metric tensor and christoffel symbols. Then to get to the AthenaPK version of those equations, when you get there, it would be straightforward to write a translation layer to represent them as these h and dh type terms. That would be the most general solution from the perspective of a library. Alternatively we could provide them as-is with new names, or not provide them. I lean towards the geometric description. Curious what other people think. (And if folks are interested I could prepare a very brief primer on the fomalism and how it relates specifically to the coordinates stuff in a finite volumes context.)

One thing I wanted to ask about---What is the team's perspective on what parthenon should support broadly? Should parthenon only support Cartesian, Cylindrical, and Spherical, but provide enough hooks for downstream codes to add more generic objects? Or should we design with more generality in mind? Or should we just restrict ourselves to things that are no more complex than Spherical? One reason I ask is that the geometric language I describe above is dramatically simpler for the special cases we are considering here than in full generality. My lenaing is to support only Cartesian, cylindrical, and spherical in parthenon proper, but try to design with letting downstream codes do more general things if they want.

How should coordinate terms be cached to save computation?

This may depend on CPU vs GPU. IMO the reason to cache is that sin and cos are transcendental functions, and thus expensive to compute. But for these simple coordinate systems, we may discover that memory bandwidth vs flops tradeoffs mean that the cache is actually slower.

Additionally, the current implementation of cached coordinate data in this PR is likely broken. For now, we could replace the cached data with functions so that for initial testing everything is computed on the fly.

Agreed. Let's just do on-the-fly and see if we need caching at all.

Is there a simple test of coordinate systems that we can add to Parthenon? Maybe a disk advection problem?

I'd be in favor of that. A more strenuous test might be can we advect a Gaussian blob (with a constant linear velocity) through a spherical/cylindrical coordinate system. I.e., it should still go in a straight line, despite the coordinate system being curved.

[Yurlungur]

For completeness, the metric tensor for cylindrical coordinates is:
$$\gamma = diag(1, r^2, 1)$$
and for spherical coordinates:
$$\gamma = diag(1, r^2, r^2 \sin^2\theta)$$

The h* terms in Athena++ are the so-called scale factors which are the square root of the diagonal components of the metric---they're only defined for diagonal metrics, which is why they're only in the Cartesian, cylindrical, and spherical coordinate objects. I would be supportive of including them over the metric if they are called scale factor rather than h as I think the latter is ambiguous.

The derivatives of the scale factors contain similar information to the christoffel symbols. The Christoffel symbols in cylindrical coordinates have essentially one nonzero term:
$$\Gamma^r_{\theta\theta} = - \Gamma^\theta_{r r} = -\frac{1}{r}$$
In Spherical coordinates, they have the following nonzero terms:
$$\Gamma^r_{\theta\theta} = \Gamma^r_{\phi\phi} = - \Gamma^\theta_{r r} = -\Gamma^\phi_{r\phi} = -\frac{1}{r}$$
and
$$\Gamma^\theta_{\theta\phi} = -\Gamma^\phi_{\theta\theta} = \frac{\cot\theta}{r}$$
I would be open to including the Christoffel symbols themselves or the gradients of the scale factors if they're named descriptively.

[Yurlungur]

Oh one final thought: We might want to add functions that transform coordinate positions to Cartesian positions for the purpose of vis tools... And we may also want Jacobians to transform to Cartesian coordinates for visualizing vectors. Those are just for outputs though... not needed in an actual simulation. And they could be encoded instead in, e.g., yt.

[forrestglines]

A few comments below. But a bigger picture comment also: would you be willing to chat with myself, Josh, and Luke about what's actually needed here? We were chatting today and we're pretty sure that (at least in a finite volume context) this interface can be simplified.

Thank you so much for looking over this PR. Still early stages and we're open to changes/suggestions.

I'd be open to chatting if we think there's a need. Talking others in T-2 though, our initial guiding principle is "Make it work like Athena++ first". Since we've generally followed Athena++'s AMR strategy and AthenaPK's hydro implementation follows Athena++, following Athena++'s curvilinear strategy hopefully should be the minimal set of changes.

Ideally in the future we'd like swap linear momentum for angular momentum in the conserved variables. I'm not sure how many changes that would require to Parthenon and AthenaPK and might end up in a code separate to AthenaPK.

[forrestglines]

There is a mathematical formalism for describing these relevant quantities. In particular the h* terms are formally described as part of the metric tensor, which describes the physical magnitudes of vectors in curved space and, as a corollary, describes lengths of paths, which are the integrals of the magnitudes of tangent vectors. The dhvd type terms are formally described by the Christoffel symbols these are essentially rotation coefficients describing how a vector rotates as you move in a coordinate direction (for example ).

The h* terms in Athena++ are the so-called scale factors which are the square root of the diagonal components of the metric

Ah, thank you! This was the insight I was hoping for, seeing as I'm newer to using different coordinate systems.

Since they're only used in specific use cases, it might be better to spin these off into different objects templated on coordinate system i.e. ScaleFactors<class Coords> and Christoffels<class Coords> as we need them. I'm not sure how this would be implemented with cached coordinate data if we decide to do that later (make them member objects of Coordinate classes?) but we can figure that out later.

[forrestglines]

I think in the "topological element" framework that @lroberts36 developed, these become centroid values based on different topological elements, e.g., Face1 or Face2. We probably want to rethink the naming convention here... s for "centroid" and c for cell-center is in my opinion confusing. I would rather switch to X for coordinate location (templated on topological element) and Xc for centroid (templated on topological element).

Bit of clarification, in this PR Xc becomes "X-centriod" (center of volume) instead of just "X-center" (midpoint of coordinates). For cylindrical the centroid r coordinate differs from the midpoint r and for spherical the centroid r and centroid theta (azimuthal?) differ from their midpoints. As far as I can tell, midpoints don't appear in Athena++, centoid positions are always used. For the topological element framework X function that @lroberts36 made, I think cell-centers should return centoids since we apparently don't use midpoints.

Xs is a "face-averaged" position which is apparently different from the centoid position for only the specific case of the r coordinate in the spherical coordinate system and is only used for prolongation of face-centered fields. I chose Xs because it matches Athena++'s convention of x1s2 and Xfa would confuse it with the "function argument" convention we're already using. I'd be open to suggestions of better names for "face-averaged X"

[Yurlungur]

Bit of clarification, in this PR Xc becomes "X-centriod" (center of volume) instead of just "X-center" (midpoint of coordinates). For cylindrical the centroid r coordinate differs from the midpoint r and for spherical the centroid r and centroid theta (azimuthal?) differ from their midpoints. As far as I can tell, midpoints don't appear in Athena++, centoid positions are always used.

Interesting... I'm not sure, but I think we may want to retain the ability to use coordinate midpoint values in addition to the centroid values. I can see that being valuable in some contexts---for example the volume-weighted interpolation (densitized) if anyone wants to use that. That said maybe that philosophy and the centroid philosophy shouldn't meet.

Curious what other folks think about that.

Xs is a "face-averaged" position which is apparently different than the centoid position for only the specific case of the r coordinate in the spherical coordinate system and is only used for prolongation of face-centered fields.

I'm not sure I understand this---my thinking was that a "face-averaged" position would be just the centroid position based at the face topological element. I was thinking there was sort of a generic generalization of the cell centroids for each topological element which matches that definition. What am I missing here?

[forrestglines]

@Yurlungur here's the two lines where the centroids and so called "face-averaged quantities" differ in spherical coordinates. I'm headed to a lunch meeting so I'll have think about why these two are different later

[Yurlungur]

Ah I see---the former is volume integrated, the latter surface-integrated.

I guess then the question is philosophical/API: If I ask for X1 "centroid" value with topological element "face2" do I ever want the volume-averaged version? Or do I always want the face-averaged version? I think the answer is we always want the face-averaged version?

[forrestglines]

I like that logic, the face position should probably be the face-averaged position.

The other place I can think of where face positions are relevant in AthenaPK is when we define a vector potential on faces then take a centered finite difference curl to get the magnetic field. We'd like the face position to correspond to a divergenceless curl and match the magnetic field we'd like to set as close as possible. I'm not sure if that would mean face-area averaged or cell-volume-averaged for the proper position on the face.

[Yurlungur]

The other place I can think of where face positions are relevant in AthenaPK is when we define a vector potential on faces then take a centered finite difference curl to get the magnetic field.

You don't define it on edges? I would think to get a divergence free B-field, you would take the vector potential on edges, curl it, and get B-field on faces?

We'd like the face position to correspond to a divergenceless curl and match the magnetic field we'd like to set as close as possible. I'm not sure if that would mean face-area averaged or cell-volume-averaged for the proper position on the face.

I don't know off the top of my head. Another thing to consider is reconstructions... There you want to use the volume-averaged centroid positions to interpolate to the coordinate centers.

[lroberts36]

I think for conservative prolongation on volumes, faces, and edges, you always want the centroid coordinates of the topological element itself (not the volume average for faces, for instance). At least at linear order choosing to model your function in the cell as $f(\vec{x}) = f_\textrm{centroid} + a_i (x^i - \bar x^i)$ with

$$\bar x^i = \frac{\int_\text{element} d^nx \sqrt{\gamma} x^i}{\int_\text{element} d^nx \sqrt{\gamma}}$$

where $\gamma_{ij}$ is the induced metric on the element, results in the volume averaged value of the prolongation function being equal to $f_\textrm{centroid}$.

@Yurlungur: Why do you necessarily want to reconstruct to the face coordinate centers instead of the face centroids?

[Yurlungur]

Maybe I just wasn't thinking it through. Your argument makes sense to me @lroberts36

[forrestglines]

@lroberts36 So if I understand correctly, we'd likely always want the position on the face to be the centroid of the face-area, not a centroid based on cell-volume. That's consistent with how Athena++ usually treats face positions although I'm suspicious that they've mixed them up in some place.

@Yurlungur Since we're cell centered the magnetic fields are currently cell-centered. That won't be the case once a CT-MHD scheme makes it into AthenaPK

[lroberts36]

@forrestglines: For prolongation, I think it is always true that we want the centroids of the topological element that we are prolongating on (at least for the default prolongation operator).

For reconstruction, we will have cell values and want to reconstruct them to faces (or for some upwind MHD algorithms have face values and want to reconstruct them on edges). The cell averaged value is probably most naturally interpreted as a point-wise value at the cell centroid. Then the reconstructed values on the face will be at two of the cell centroid coordinates and one of the face centroid coordinates, e.g. if I am reconstructing on the face normal to the $\theta$ direction the reconstructed position on the face will be $(r_\textrm{cell-centroid}, \theta_\textrm{face}, \phi_\textrm{cell-centroid})$. This maybe poses some issues for the order of the method, since for second order accuracy we would probably like to solve point-wise Riemann problems at the face centroids[^1], but maybe in practice this isn't a large source of error.

All that being said, at least for Cartesian, cylindrical, and spherical coordinates where the radial position only depends on the $i$ index, etc., all of this information can easily be recovered from a single centroid coordinate template function

Coordinates::Xc<CoordinateDirection dir, TopologicalElement el>(int k, int j, int i);

so that e.g. the reconstructed position of the lower face of a zone normal to the $\theta$ direction is given by {coords.Xc<X1DIR, TE::CC>(k, j, i), coords.Xc<X2DIR, TE::F2>(k, j, i), coords.Xc<X3DIR, TE::CC>(k, j, i)}.

[1:] I would be interested to know how/if Athena++ deals with this.

[Yurlungur]

I'd be open to chatting if we think there's a need. Talking others in T-2 though, our initial guiding principle is "Make it work like Athena++ first". Since we've generally followed Athena++'s AMR strategy and AthenaPK's hydro implementation follows Athena++, following Athena++'s curvilinear strategy hopefully should be the minimal set of changes.

I think it might be worth touching base at some point. I think your strategy makes sense for AthenaPK. From the parthenon perspective, I think we want the minimal feature set that spans AthenaPK's needs and the needs of other downstream codes that might want to use this machinery.

Ideally in the future we'd like swap linear momentum for angular momentum in the conserved variables. I'm not sure how many changes that would require to Parthenon and AthenaPK and might end up in a code separate to AthenaPK.

Interesting... That probably requires coordinate transformation matrices (i.e., the Jacobian). There may be some need for that anyway in pgens and vis tools though.

Since they're only used in specific use cases, it might be better to spin these off into different objects templated on coordinate system i.e. ScaleFactors and Christoffels as we need them. I'm not sure how this would be implemented with cached coordinate data if we decide to do that later (make them member objects of Coordinate classes?) but we can figure that out later.

That's an interesting idea. Broadly I think this makes sense. We should think about how that fits into the meshblock and pack machinery though to make sure they're accessible where they're needed.

Alternatively, one could imagine caching them on the mesh as variables. Maybe some automated machinery for computing scale factors/christoffels/etc centered at relevant topological elements on the grid? We could just treat them as actual mesh variables? That might cause them to have pretty nonstandard shapes per meshblock though.

[forrestglines]

I'm currently working on a private branch of AthenaPK using these constructs in Parthenon. Once I have something partially working I can share that private branch to interested LANL employees and make an advection example in Parthenon.

I think we might have enough for cell-centered MHD in this PR already. I'm finding some of the functions in the coordinate class (like the Coord_src terms) might better belong downstream. If we decide to cache these values we'll have to show an example in Parthenon on how to do that.