-
Notifications
You must be signed in to change notification settings - Fork 107
Improvements for hermitespline interpolation #3298
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Changes from all commits
5502ce3
c74d330
afc68d9
84237c8
0662449
f77849e
820f0a5
1178243
3d6f587
265c8cd
16a162e
e5878c1
aee2713
b4f7be5
578ee8d
bf9c9eb
2155f06
0662247
28d4e28
e4507d9
d7ba297
f55f8e4
0e9f7d6
944adfb
6daec20
3444dd0
91bee68
d7aeae4
ac5cc70
885b465
a1c4033
1c25b0e
c7f8504
96aa2bd
65880e4
a91d3e7
f473a6b
f341037
999d4ca
614a727
d1cfb8a
5c0d41d
eaf35d9
a5392af
b8bfd3b
52373e7
c650881
3a62dc4
cb84b76
a81d3eb
4632df0
4b0ed1c
325dde3
8dc7b4f
db53852
29f9727
dcafb6c
424136b
500989b
cabf678
01e2263
6bd15d2
f2d860f
159ce2f
1de6329
47f93b1
58db33a
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -27,17 +27,18 @@ | |
| #include <bout/bout_types.hxx> | ||
| #include <bout/generic_factory.hxx> | ||
| #include <bout/mask.hxx> | ||
| #include <array> | ||
|
|
||
| #define USE_NEW_WEIGHTS 1 | ||
| #if BOUT_HAS_PETSC | ||
| #define HS_USE_PETSC 1 | ||
| #endif | ||
|
|
||
| #ifdef HS_USE_PETSC | ||
| #include "bout/petsclib.hxx" | ||
| #endif | ||
|
|
||
| namespace { | ||
| enum class implementation_type { new_weights, petsc, legacy }; | ||
| } | ||
|
|
||
| class Options; | ||
| class GlobalField3DAccess; | ||
|
|
||
| /// Interpolate a field onto a perturbed set of points | ||
| const Field3D interpolate(const Field3D& f, const Field3D& delta_x, | ||
|
|
@@ -135,14 +136,26 @@ public: | |
| } | ||
| }; | ||
|
|
||
| class XZHermiteSpline : public XZInterpolation { | ||
| /// Monotonic Hermite spline interpolator | ||
| /// | ||
| /// Similar to XZHermiteSpline, so uses most of the same code. | ||
| /// Forces the interpolated result to be in the range of the | ||
| /// neighbouring cell values. This prevents unphysical overshoots, | ||
| /// but also degrades accuracy near maxima and minima. | ||
| /// Perhaps should only impose near boundaries, since that is where | ||
| /// problems most obviously occur. | ||
| template <bool monotonic, implementation_type imp_type> | ||
| class XZHermiteSplineBase : public XZInterpolation { | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: class 'XZHermiteSplineBase' defines a non-default destructor but does not define a copy constructor, a copy assignment operator, a move constructor or a move assignment operator [cppcoreguidelines-special-member-functions] class XZHermiteSplineBase : public XZInterpolation {
^ |
||
| protected: | ||
| /// This is protected rather than private so that it can be | ||
| /// extended and used by HermiteSplineMonotonic | ||
|
|
||
| Tensor<SpecificInd<IND_TYPE::IND_3D>> i_corner; // index of bottom-left grid point | ||
| Tensor<int> k_corner; // z-index of bottom-left grid point | ||
|
|
||
| std::unique_ptr<GlobalField3DAccess> gf3daccess; | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: member variable 'gf3daccess' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes] std::unique_ptr<GlobalField3DAccess> gf3daccess;
^ |
||
| Tensor<std::array<int, 4>> g3dinds; | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: member variable 'g3dinds' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes] Tensor<std::array<int, 4>> g3dinds;
^
dschwoerer marked this conversation as resolved.
|
||
|
|
||
| // 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 | ||
|
|
@@ -160,23 +173,28 @@ protected: | |
|
|
||
| std::vector<Field3D> newWeights; | ||
|
|
||
| #if HS_USE_PETSC | ||
| #if BOUT_HAS_PETSC | ||
| PetscLib* petsclib; | ||
| bool isInit{false}; | ||
| Mat petscWeights; | ||
| Vec rhs, result; | ||
| Mat petscWeights{}; | ||
|
dschwoerer marked this conversation as resolved.
|
||
| Vec rhs{}, result{}; | ||
|
dschwoerer marked this conversation as resolved.
dschwoerer marked this conversation as resolved.
|
||
| #endif | ||
|
|
||
| /// Factors to allow for some wiggleroom | ||
| BoutReal abs_fac_monotonic{1e-2}; | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: member variable 'abs_fac_monotonic' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes] BoutReal abs_fac_monotonic{1e-2};
^ |
||
| BoutReal rel_fac_monotonic{1e-1}; | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: member variable 'rel_fac_monotonic' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes] BoutReal rel_fac_monotonic{1e-1};
^ |
||
|
|
||
| public: | ||
| XZHermiteSpline(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr) | ||
| : XZHermiteSpline(0, mesh) {} | ||
| XZHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr); | ||
| XZHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) | ||
| : XZHermiteSpline(y_offset, mesh) { | ||
| XZHermiteSplineBase(Mesh* mesh = nullptr, [[maybe_unused]] Options* options = nullptr) | ||
| : XZHermiteSplineBase(0, mesh, options) {} | ||
| XZHermiteSplineBase(int y_offset = 0, Mesh* mesh = nullptr, Options* options = nullptr); | ||
| XZHermiteSplineBase(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr, | ||
| Options* options = nullptr) | ||
| : XZHermiteSplineBase(y_offset, mesh, options) { | ||
| setRegion(regionFromMask(mask, localmesh)); | ||
| } | ||
| ~XZHermiteSpline() { | ||
| #if HS_USE_PETSC | ||
| ~XZHermiteSplineBase() override { | ||
| #if BOUT_HAS_PETSC | ||
| if (isInit) { | ||
| MatDestroy(&petscWeights); | ||
| VecDestroy(&rhs); | ||
|
|
@@ -205,61 +223,21 @@ public: | |
| getWeightsForYApproximation(int i, int j, int k, int yoffset) override; | ||
| }; | ||
|
|
||
| /// Monotonic Hermite spline interpolator | ||
| /// | ||
| /// Similar to XZHermiteSpline, so uses most of the same code. | ||
| /// Forces the interpolated result to be in the range of the | ||
| /// neighbouring cell values. This prevents unphysical overshoots, | ||
| /// but also degrades accuracy near maxima and minima. | ||
| /// Perhaps should only impose near boundaries, since that is where | ||
| /// problems most obviously occur. | ||
| /// | ||
| /// You can control how tight the clipping to the range of the neighbouring cell | ||
| /// values through ``rtol`` and ``atol``: | ||
| /// | ||
| /// diff = (max_of_neighours - min_of_neighours) * rtol + atol | ||
| /// | ||
| /// and the interpolated value is instead clipped to the range | ||
| /// ``[min_of_neighours - diff, max_of_neighours + diff]`` | ||
| class XZMonotonicHermiteSpline : public XZHermiteSpline { | ||
| /// Absolute tolerance for clipping | ||
| BoutReal atol = 0.0; | ||
| /// Relative tolerance for clipping | ||
| BoutReal rtol = 1.0; | ||
|
|
||
| public: | ||
| XZMonotonicHermiteSpline(Mesh* mesh = nullptr, Options* options = nullptr) | ||
| : XZHermiteSpline(0, mesh), | ||
| atol{(*options)["atol"] | ||
| .doc("Absolute tolerance for clipping overshoot") | ||
| .withDefault(0.0)}, | ||
| rtol{(*options)["rtol"] | ||
| .doc("Relative tolerance for clipping overshoot") | ||
| .withDefault(1.0)} { | ||
| if (localmesh->getNXPE() > 1) { | ||
| throw BoutException("Do not support MPI splitting in X"); | ||
| } | ||
| } | ||
| XZMonotonicHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr) | ||
| : XZHermiteSpline(y_offset, mesh) { | ||
| if (localmesh->getNXPE() > 1) { | ||
| throw BoutException("Do not support MPI splitting in X"); | ||
| } | ||
| } | ||
| XZMonotonicHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) | ||
| : XZHermiteSpline(mask, y_offset, mesh) { | ||
| if (localmesh->getNXPE() > 1) { | ||
| throw BoutException("Do not support MPI splitting in X"); | ||
| } | ||
| } | ||
|
|
||
| using XZHermiteSpline::interpolate; | ||
| /// Interpolate using precalculated weights. | ||
| /// This function is called by the other interpolate functions | ||
| /// in the base class XZHermiteSpline. | ||
| Field3D interpolate(const Field3D& f, | ||
| const std::string& region = "RGN_NOBNDRY") const override; | ||
| }; | ||
| using XZMonotonicHermiteSplineSerial = | ||
| XZHermiteSplineBase<true, implementation_type::new_weights>; | ||
| using XZHermiteSplineSerial = | ||
| XZHermiteSplineBase<false, implementation_type::new_weights>; | ||
| using XZMonotonicHermiteSplineLegacy = | ||
| XZHermiteSplineBase<true, implementation_type::legacy>; | ||
| using XZHermiteSplineLegacy = XZHermiteSplineBase<false, implementation_type::legacy>; | ||
| #if BOUT_HAS_PETSC | ||
| using XZMonotonicHermiteSpline = XZHermiteSplineBase<true, implementation_type::petsc>; | ||
| using XZHermiteSpline = XZHermiteSplineBase<false, implementation_type::petsc>; | ||
| #else | ||
| using XZMonotonicHermiteSpline = | ||
| XZHermiteSplineBase<true, implementation_type::new_weights>; | ||
| using XZHermiteSpline = XZHermiteSplineBase<false, implementation_type::new_weights>; | ||
| #endif | ||
|
|
||
| /// XZLagrange4pt interpolation class | ||
| /// | ||
|
|
||
| Original file line number | Diff line number | Diff line change | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -534,6 +534,8 @@ int BoutMesh::load() { | |||||||||
| PE_XIND = MYPE % NXPE; | ||||||||||
| PE_ZIND = 0; | ||||||||||
|
|
||||||||||
| ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND)); | ||||||||||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: no header providing "ASSERT2" is directly included [misc-include-cleaner] ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND));
^ |
||||||||||
|
|
||||||||||
| // Set the other grid sizes from nx, ny, nz | ||||||||||
| setDerivedGridSizes(); | ||||||||||
|
|
||||||||||
|
|
@@ -1015,6 +1017,8 @@ void BoutMesh::createXBoundaries() { | |||||||||
| } | ||||||||||
| } | ||||||||||
|
|
||||||||||
| int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; } | ||||||||||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
Suggested change
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: parameter 'Z' is unused [misc-unused-parameters]
Suggested change
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. warning: unused parameter 'Z' [clang-diagnostic-unused-parameter] int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }
^ |
||||||||||
|
|
||||||||||
| void BoutMesh::createYBoundaries() { | ||||||||||
| if (MYG <= 0) { | ||||||||||
| return; | ||||||||||
|
|
||||||||||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: enum 'implementation_type' uses a larger base type ('int', size: 4 bytes) than necessary for its value set, consider using 'std::uint8_t' (1 byte) as the base type to reduce its size [performance-enum-size]