Proof of concept for using variant metric elements

[ZedThree]

I've made a proof of concept for using variant<BoutReal, Field2D, Field3D> for the Coordinate metric elements, grid spacing and so on. There are two main reasons why this approach is attractive:

1. runtime choice of metric dimension
2. constant metric elements could be easily detected

The second reason is interesting for turning off identically-zero terms, for example the g13 * ::D2DXDZ term in Delp2, or for enabling FFT in any direction. Set dx = 1.0, and we would know that it's constant and FFTs could be used.

The timings are very promising! These are the results of the average of 50 loops of auto result = DDX(f), with 128x128x128 points:

Method	Avg. time (s)
current DDX	0.351670
variant DDX with BoutReal	0.371770
variant DDX with Field2D	0.422765
variant DDX with Field3D	0.410582
hardcoded DDX with BoutReal	0.371798
hardcoded DDX with Field2D	0.422805
hardcoded DDX with Field3D	0.410404

With both the variant and hardcoded DDX, I'm passing in dx.

There's basically no difference between using the variant and using a hardcoded type. The current DDX is a bit faster because it actually uses result /= coords->dx, and we have some optimisations for Field3D::operator/=(Field2D), whereas all the other results above use result = DDX(f) / dx. This isn't a barrier to a real implementation, just what was the fastest to get implemented.

The implementation looks something like this:

using metric_type = bout::utils::variant<BoutReal, Field2D, Field3D>;

template <class T>
T operator/(const T &field, const metric_type &metric) {
  metric_type lhs = field;
  auto result = bout::utils::visit(
      [](const auto &lhs, const auto &rhs) -> metric_type { return lhs / rhs; },
      lhs, metric);
  return bout::utils::get<T>(result);
}

// Just for demonstration, real implementation follows current DDX more closely
Field3D variant_DDX(const Field3D &f, const metric_type &dx) {
  return bout::derivatives::index::DDX(f, outloc, method, region) / dx;
}

Thanks to C++14's generic lambdas, this was surprisingly easy to implement!

One downside of this approach is if things go wrong, the error messages are generally not very helpful to the average user.

A second issue, and I think this is actually common to all the 3D coords approaches we've looked at, is that what should the return type of DDX(Field2D) be? The current approach from @bshanahan and @d7919 sets this to be Coordinates::metric_type, which is set at configure-time to be either Field2D or Field3D.

With the variant approach, this could either remain as Field2D and throw an exception if the metric is 3D, which is what I've done in the implementation above. I think this makes a certain amount of sense to do as you probably don't want such operations with a 3D metric. There may be some circumstances where you do want them however, and this could be a problem.

The other option is to return a variant, but then this has a downside of "infecting" everything downstream of it, as it were, and I'd worry about lots of seemingly nonsense errors for the user.

Lastly, there's an issue with indexing the variant. This could be solved by wrapping up the variant in a class, and providing operator(), operator[], which would Do The Right Thing.

[bshanahan]

This looks great! The only issue I can see at the moment is with something like DDX(FieldPerp), which would likely need 3D metrics if available.

[bendudson]

Nice! Thanks @ZedThree this looks good. I think I agree that a Field2D input implies some kind of optimisation is going on, and it should throw an exception rather than returning a variant. Otherwise this goes more and more towards dynamic typing with the overhead that usually comes with it (and with worse error messages).

A wrapper with indexing could be nice, but I think it might be better not to do this: Usually the indexing is used when there is a need to optimise, and having a switch inside operator[] in the inner loop would probably be very slow. Perhaps wrapper functions to make the process of testing the type of the metric would be nice, though std::holds_alternative isn't too bad. This would allow code to switch between versions, depending on the metric type.

It would be nice (for some value of nice) to have some combination of template and/or macro magic, which would compile multiple versions of a block of code, and put in a switch. Inside each instance of the code, the type would be BoutReal, Field2D, or Field3D. Probably with BoutReal we could define a thin wrapper with an operator[] which could be optimised away by the compiler.

e.g. argument list is list of variant types. Code is compiled for each combination (warning: combinatorial explosion ahead...)

VARIANT_CODE(dx, dy, ...) {
   // Here use dx, dx[i] etc. 
}

[ZedThree]

@bshanahan I don't see a DDX(FieldPerp)? At any rate, it's possible to check exactly which type the variant is holding at runtime, so it's possible to throw if the metric isn't the correct type.

@bendudson That could probably be done with visitors? At least the implementation of such a macro maybe.

This might be a sufficient wrapper for BoutReal:

struct BoutScalar {
  BoutReal value;
  constexpr BoutReal operator[](Ind3D) const { return value; }
  constexpr BoutReal operator()(int, int) const { return value; }
  constexpr BoutReal operator()(int, int, int) const { return value; }
}

I think the current coords 3d approach is #ifdef out bits that won't work with 3D metrics, and to index everything with either Ind3D or three indices. This would probably work with the above, except instead of #ifdef, if (holds_alternative<Field3D>(dx)) throw.

[bendudson]

Yes, possibly a visitor with some template parameters. Nested visitors if multiple metric components are being used.