Redesigning the Coordinate classes

[tomidakn]

Currently, we use derived classes for Coordinate classes. This is a good C++ design as they can share the base functions and define necessary functions for each coordinate systems. However, there are two major issues with this design:

1. GPU compilers (NVIDIA in particular) are not compatible with virtual functions in GPU kernels.
2. Compiling all the coordinate files causes long compilation time.

Now I am talking to NVIDIA (Japan) and discussing how to port Athena++ to GPUs with minimum change (possibly using OpenACC), and they pointed out the first issue. The coordinate is selected at the configuration, and never changed at runtime. In that sense, it does not have to be a derived class.

So, I am thinking of redesigning the Coordinate classes. Like the Riemann solver files, we select one coordinate file in the configuration script. Each Coordinate class should have all the functions defined, but we can implement default ones using the "attribute((weak))" feature, which we already used for problem generators. How do you think? I would like to ask @felker and @c-white in particular, but everyone is welcome.

Here, I just want to make Athena++ compatible with GPUs (and CPUs) and I do not intend to get maximum performance, unlike AthenaK. Simply there are many features in Athena++ not ported to AthenaK, and if some of them can be used on GPUs with minimum change, I guess there will be a lot of demand.

[felker]

What about the virtual functions in bvals/? There are also minor uses of it in Multigrid, Chemistry, FFT, Outputs

[tomidakn]

I do not worry about bvals because typically they are not in GPU computing kernels. NVIDIA pointed out that the coordinate functions are called in the Riemann solvers and CT, which is a problem in GPU porting. (Note: nvcc is so primitive here, it does not return any compile nor runtime error but just call an empty function which does nothing.)

We need to think about other physics modules, but they can be done later.

[c-white]

If this only affects Coordinates, then this seems like a desirable thing if it allows running on a GPU.

If weak attributes work, then this seems straightforward.

Even if it doesn't, a little code duplication isn't the worst thing in this case. For example, currently someone reading the spherical polar file might wonder how Edge1Length is defined, or might not even know it exists, since the formula matches Cartesian coordinates and so is hidden in the base class.

[tomidakn]

Yes, it should be straightforward. I guess I will try implementing one or two coordinates and see how it works in a few weeks.

[tomidakn]

@c-white I noticed that "sinusoidal" and "tilted" options are still there in configure.py but there are no corresponding files in the code directory. Should I remove these options?

[c-white]

Yes, the options should have been removed some time ago.

[tomidakn]

OK, thanks.

[tomidakn]

...I will keep those options in athena_read.py just in case someone wants to load those old data with it.