PETSc XZHermiteSpline #2917
PETSc XZHermiteSpline #2917
Conversation
Reduce duplication by moving some overloads to the base class.
Optimised build was omitting XZInterpolation factory registrations.
* next: (194 commits) invert_laplace: Fix outer boundary metrics Apply clang-format changes Apply clang-format changes Apply suggestions from code review Fix missed adios Prefer ADIOS2 over ADIOS in more cases Update tests/unit/test_extras.hxx Apply clang-format changes Silence warning about assignment in `if` statement Fix some short identifiers Make a couple of implicit casts explicit CI: Suppress some warnings about short identifiers Move some static functions to anonymous namespace CI: Remove deprecated clang-tidy option Include header for `BoutReal`; remove unused header Apply clang-format Add options to set Butcher tables by name Simplify SUNDIALS cmake Switch to [[maybe_unused]] Apply suggestions from code review ...
Based on @dschwoerer's implementation in #2651
Currently gives pretty bad results, not sure why
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
|
Ok, I have managed to convert this to use the So I've implemented a new |
|
clang-tidy review says "All clean, LGTM! 👍" |
Although `getGlobal` is only valid for local indices, unfortunately `isLocal` calls `getGlobal`, so we can't actually check!
|
clang-tidy review says "All clean, LGTM! 👍" |
There was a problem hiding this comment.
A few things are indeed much nicer then with the abstraction.
I am not sure whether z-splitting will be tricky. Right now it works because only x is split, and x is the slowest index.
One advantage with the alternative implementation is, that the monotonic version is implementation for all hermite spline interpolation, although we probably want to keep at most the non-petsc version, although I am not sure of that, we might just say: Get PETSc if you want to do FCI ...
Have you by chance did a performance measurement? This implementation should copy slightly less data ...
| // Actually, I'm not 100% sure what this does, but this appears to the minimum | ||
| // necessary to get things to work |
There was a problem hiding this comment.
I am not to sure either.
I have to say I am worried, if no-one understands how this works ...
Why do we not need .xp(2).zp(2) and similar? Does that include the corners, or just the boundaries?
| using RegisterUnavailableXZInterpolation = | ||
| XZInterpolationFactory::RegisterUnavailableInFactory; | ||
|
|
||
| namespace { | ||
| const inline RegisterXZInterpolation<XZHermiteSpline> registerinterphermitespline{ | ||
| "hermitespline"}; | ||
| const inline RegisterXZInterpolation<XZMonotonicHermiteSpline> | ||
| registerinterpmonotonichermitespline{"monotonichermitespline"}; | ||
| const inline RegisterXZInterpolation<XZLagrange4pt> registerinterplagrange4pt{ | ||
| "lagrange4pt"}; | ||
| const inline RegisterXZInterpolation<XZBilinear> registerinterpbilinear{"bilinear"}; | ||
|
|
||
| #if BOUT_HAS_PETSC | ||
| const inline RegisterXZInterpolation<PetscXZHermiteSpline> registerpetschermitespline{ | ||
| "petschermitespline"}; | ||
| #else | ||
| const inline RegisterUnavailableXZInterpolation register_unavailable_petsc_hermite_spline{ | ||
| "petschermitespline", "BOUT++ was not configured with PETSc"}; | ||
| #endif | ||
| } // namespace | ||
|
|
There was a problem hiding this comment.
Should they really be in the header?
| const int xend = (localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE() | ||
| + localmesh->xstart - 1; | ||
|
|
||
| // TODO: work out why using `getRegion(region)` directly causes |
There was a problem hiding this comment.
is this still an issue?
| // Basis functions for cubic Hermite spline interpolation | ||
| // see http://en.wikipedia.org/wiki/Cubic_Hermite_spline | ||
| // The h00 and h01 basis functions are applied to the function itself | ||
| // and the h10 and h11 basis functions are applied to its derivative | ||
| // along the interpolation direction. | ||
| Field3D h00_x; | ||
| Field3D h01_x; | ||
| Field3D h10_x; | ||
| Field3D h11_x; | ||
| Field3D h00_z; | ||
| Field3D h01_z; | ||
| Field3D h10_z; | ||
| Field3D h11_z; |
There was a problem hiding this comment.
They should only be used in the constructor, so no need to keep them around
| // Need to convert from global indices to local | ||
| const auto i_c = Field3D::ind_type( | ||
| (((localmesh->getLocalXIndex(i_corn) * ny) + (y + y_offset)) * nz | ||
| + localmesh->getLocalZIndex(k_corner(x, y, z))), |
There was a problem hiding this comment.
Note that this seems like it would work with Z-splitting, but it does not!
This is inherently incompatible with z-splitting. It only works if ny and nz are the same as the local ...
| const PetscVector<Field3D> f_result = weights * f_petsc; | ||
| // TODO: we should only consider the passed-in region | ||
| // const auto region2 = y_offset == 0 ? region : fmt::format("RGN_YPAR_{:+d}", y_offset); | ||
| return f_result.toField(); |
There was a problem hiding this comment.
This does not set the region - is this set in fci.cxx/fci.hxx?
Otherwise this will not work with high check level and 3d metrics, afaik.
Continues #2858 pulling out a separate class for
PetscXZHermiteSplinethat was started in #2651This does have a working implementation, but I've tried to switch it the existing
PetscMatrixwrapper, and unfortunately I've gone wrong somewhere -- probably at the boundaries.Double unfortunately, I really don't have much time to spend digging into this. @dschwoerer Maybe you want to take a look?