Clean up Coordinates() constructor

[johnomotani]

At the moment the Coordinates::Coordinates() constructor (plus some utility methods/functions) does a lot of interpolation, extrapolation and checking. This was useful to support new features, but has become more complicated and hard to maintain. It is also difficult for users to find and fix problems with grids if they arise. Extrapolation (especially) is still useful, because geqdsk files may not define the poloidal flux function far enough past the walls for a grid generator to create a grid with boundary points, or the grid may run into coils/singularities, etc.

A proposal for discussion:

1. Change Coordinates::Coordinates() so that it does no extrapolation or interpolation. Grid files are required to include boundary points, and if staggered Coordinates are created, the grid file must include metrics, etc. at the staggered grid locations. It should check for positivity, etc. as appropriate, and also require sane values in corner guard/boundary points.
2. Write a grid file validation and conversion tool. This would take a grid file that might be missing staggered quantities and/or boundary points and apply extrapolation/interpolation to generate them. It should also allow users to easily check grid files and fix problems.

Desired features of the conversion tool:

• A GUI similar to hypnotoad-gui would be ideal, that like hypnotoad can also create an input file for reproducible command-line conversions.
• Can take input grid files with/without staggered quantities and with/without boundary cells.
• To start with have at least the same extrapolation/interpolation algorithms currently implemented in coordinates.cxx, but would be easier to add new ones (e.g. to ensure positivity in problematic cases).
• Check for errors (same checks as Coordinates::Coordinates(), e.g. negative B, etc. Ideally provide some options to fix the problems.
• Would be nice to have functionality for further validation of grid files, e.g. a standard set of plots so users can check visually that there are no excessively large or small values, and that quantities are smooth. Probably a combination of 2d color-plots for all quantities and 1d line plots for selected quantities (e.g. dx, hthe/hy).
  • Edit would be nice to have a check of the 'no slip condition' from Dudson et al. (2009): Delta psi * Delta theta * d nu / d psi <= Delta zeta
• Possibly this could be implemented as an extra tab/mode in hypnotoad? The stuff above would certainly always be useful for validation of hypnotoad outputs.

See also #1472 - I think addressing that would be much nicer in a Python conversion tool than as more complication in coordinates.cxx.

[johnomotani]

Help wanted!

Unfortunately I don't see myself having any time to work on this any time soon. I'd be very happy to review some PRs though!

[ZedThree]

I'm definitely in favour of trying to reduce the amount of code in Coordinates -- it's quite complicated now!

I'm broadly in favour of your proposal, though I wonder if it might make sense to have an option to fall back to copying the last interior grid points into the boundaries? This would let you do something with a grid file without boundaries.

If the tool can take a grid file that works with BOUT++ v4.x and make it work with v5.x with minimal fuss, something like hypnotoad --to-v5 old_grid.nc, then that makes it much less painful for users.

Essentially just moving the code to hypnotoad seems like a good first step that does that conversion. I also don't have much time, but I will see if I can bash that out at some point.

[johnomotani]

it might make sense to have an option to fall back to copying the last interior grid points into the boundaries? This would let you do something with a grid file without boundaries.

Sounds like a good idea, especially for compatibility with existing tests/examples. I'd be in favour of the default option being to raise an error, so the user has to explicitly set the boundary cells this way.

[bendudson]

Agree the construction of Coordinates needs tidying up. One worry is that PR #2025 (3D metrics) also makes changes to Coordinates, so merging might be tricky...

To avoid breaking too many things, perhaps we can make all members of Coordinates const, so they can't be changed in inconsistent ways. This also means that the constructor can't really do anything except initialise the values.

To create a Coordinates, we could have classes for the separate pieces. It would be nice if these were grouped together in ways which corresponded to mathematical objects:

• MeshSpacing with dx, dy and dz, along with their derivatives (d1_dx, d1_dy)
• ContravariantMetric with g11, g22, ...
• CovariantMetric with g_11, g_22, ...
• DerivedMetrics or maybe ConnectionCoefficients with G1_11, ..., G1, G2
• Something with J and Bxy. Perhaps also in DerivedMetrics.

This breaks up the construction, which is currently one giant function, into separate pieces. Different combinations of these could be passed to the Coordinates constructor, e.g. if CovariantMetric isn't given then construct it from ContravariantMetric.

[dschwoerer]

Merging #2025 would solve that issue :-)

For the coordinates it might be easier to move everything to a pre-processing script, that writes out all the required quantities.

That way BOUT++ can concentrate on solving the differential equations, and the preprocessing tool can guess what the user might want to use - or maybe even ask in an interactive fashion. Also it wouldn't need to work in parallel, which would simplify some things (communicate before the mesh is initialised is as far as I remember a bit tricky)

[ZedThree]

I'm very much on-board with splitting it up into smaller sub-objects.

Rather than making the members const, I'd prefer to have getters/setters that can make sure everything is consistent. So calling coords->setContravariantMetric(contravariant_meric) would recalculate the covariant metric. const members make objects not behave terribly well in std containers for instance.

Getters/setters also allow us to put off calculation of some of the geometry until the first time we need it. For instance, G1_11 etc are not needed in many simulations, so we could just not calculate them until the first coords->getConnectionCoefficients()

A preprocessing tool would be nice, but not as a replacement for doing it all inside BOUT++.

[bendudson]

• Agree that const members would break things like assignment operators. It would avoid breaking code that uses coords->g_22 but I suppose if we're breaking things then getters would be ok. Calculating things only on demand is a nice benefit.
• I was trying to put all work in constructors, rather than modifying in place, but am not opposed to having nice setters which make sure everything is consistent
• Agree that doing this in preprocessing is nice: It allows debugging, plotting and sanity checking of the metrics before they go into the code. Indeed the new python hypnotoad does write most (all?) of these components to file.
• One common transformation of the metrics is normalisation: In almost every case the input metrics need to be normalised, usually by a quantity which is only known when BOUT++ runs (e.g. reference gyro-radius, B field). To normalise the metrics and derived quantities I think only needs:
  • The normalisation of length
  • The normalisation of dx, dy and dz
• A general question of style: Should we prefer a modify in place or create new instance? i.e.
  coords->setContravariantMetric(contravariant_meric) or
  coords = coords->withContravariantMetric(contravariant_meric) ?

[bendudson]

This Coordinates stuff is painful. As a plan of attack, how about:

1. Write tests for the Coordinates in next (PR Testing of Coordinates  #2167)
2. Merge next with the new tests into 3D coordinates branch (PR 3D metrics #2025), and make sure that we're not breaking anything. If all is ok, add more tests and merge 3D metrics #2025 into next.
3. Start pulling out parts of the construction of Coordinates into separate objects

[johnomotani]

To normalise the metrics and derived quantities I think only needs:

• The normalisation of length
• The normalisation of dx, dy and dz

Just a comment, I suspect it shouldn't make any difference in the end, but in STORM we actually normalise the metric not just with lengths - we normalise 'x' with something like B0/(rho_s**2) and then that feeds through the metrics... I'm not opposed to converting to just using lengths (that could be more natural anyway), but to keep doing what we're currently doing we'd need the option to do a more general normalisation.

[bendudson]

Yes I think that's quite common, to have dx a poloidal flux. I think that could be handled ok. We just need to be able to derive the normalisation of metrics. We know that

length^2 = g_ij dx^i dx^j

so by giving the normalisation of length, together with the normalisation of dx, dy, dz, all other things can (I think) be normalised consistently. Typically dy and dz are not normalised (since they're angles) so the normalisation might be:

coords->applyNormalisation(rho_s, B0 * rho_s**2, 1.0, 1.0)  // length, dx, dy, dz

[bendudson]

Running into issues around staggered coordinates, and calls to Mesh::recalculateStaggeredCoordinates() segfaulting (coordinates null in Mesh). Rather than having Mesh holding the different staggered Coordinates, and having to have this callback to recalculate, would it make sense to have the staggered versions of the metrics inside Coordinates, so they can be invalidated and recalculated?

• When creating a Coordinates from a Mesh (or Options), if data for the staggered metrics is included, then create them. If data is not included, then only calculate them on request.
• If the cell centred coordinate is modified, invalidate all the others (as now)
• The user can request staggered components e.g.

metric = coords->getCovariantMetric(CELL_YLOW);

[johnomotani]

Rather than having Mesh holding the different staggered Coordinates, and having to have this callback to recalculate, would it make sense to have the staggered versions of the metrics inside Coordinates, so they can be invalidated and recalculated?

I think the original motivation for the callback was that depending on what the user is doing we might or might not need to invalidate the staggered coordinates: if I have only CELL_CENTRE metrics in the grid file, then it makes sense to normalise them and then (re-)calculate everything else; on the other hand if I had the staggered metrics, then I probably want to normalise both the CELL_CENTRE and staggered coordinates, so don't want the staggered ones to be re-calculated.

If moving the staggered coordinates inside the 'standard' CELL_CENTRE ones simplifies this logic then it could be a good idea. Possible con is it adds a layer to the object hierarchy (mesh->coords_centre->staggered_coords instead of just mesh->coords). If we go in the end with the design I suggested at the top, where staggered and unstaggered are both always read from the grid then staggered Coordinates don't depend on the unstaggered ones - we get rid of all the code that interpolates/extrapolates - and then the current std::vector<Coordinates> in Mesh would seem like the natural thing to do. @bendudson is the refactoring you're doing urgent, or could we do something like: 1) merge 3d metrics; 2) simplify the Coordinates constructor by requiring more stuff to be calculated/saved by the grid generator; 3) refactor into more modular classes?