diff --git a/CMakeLists.txt b/CMakeLists.txt index f57a78a14a..83dde0f1b6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -270,6 +270,7 @@ set(BOUT_SOURCES ./src/mesh/interpolation/interpolation_z.cxx ./src/mesh/interpolation/lagrange_4pt_xz.cxx ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx + ./src/mesh/interpolation/petsc_hermite_spline_xz.cxx ./src/mesh/mesh.cxx ./src/mesh/parallel/fci.cxx ./src/mesh/parallel/fci.hxx diff --git a/include/bout/globalindexer.hxx b/include/bout/globalindexer.hxx index e756ead3b2..76c60e1cb6 100644 --- a/include/bout/globalindexer.hxx +++ b/include/bout/globalindexer.hxx @@ -17,21 +17,25 @@ using IndexerPtr = std::shared_ptr>; using InterpolationWeights = std::vector; -/*! - * An object which accepts index objects produced by iterating over - * fields and returns a global index. This index can be used when - * constructing PETSc arrays. Guard regions used for communication - * between processes will have the indices of the part of the interior - * region they are mirroring. Boundaries required by the stencil will - * have unique indices. If no stencil is provided then only the - * interior region will be assigned global indices. By default, the - * indexer is fully initialised so that guard cells are communicated - * (ensuring they hold the appropriate global indices). However, by - * passing `autoInitialise = false` this behaviour can be turned off - * and the user can then manually call the `initialise()` method - * later. This can be useful for mocking/faking the class when - * testing. - */ +/// Convert local indices on current process to consistent, unique +/// global flat indices across all processes. +/// +/// Note that this can only convert indices both local to *and owned by* the +/// current process! +/// +/// An object which accepts index objects produced by iterating over +/// fields and returns a global index. This index can be used when +/// constructing PETSc arrays. Guard regions used for communication +/// between processes will have the indices of the part of the interior +/// region they are mirroring. Boundaries required by the stencil will +/// have unique indices. If no stencil is provided then only the +/// interior region will be assigned global indices. By default, the +/// indexer is fully initialised so that guard cells are communicated +/// (ensuring they hold the appropriate global indices). However, by +/// passing `autoInitialise = false` this behaviour can be turned off +/// and the user can then manually call the `initialise()` method +/// later. This can be useful for mocking/faking the class when +/// testing. template class GlobalIndexer { public: @@ -98,7 +102,7 @@ public: } } - virtual ~GlobalIndexer() {} + virtual ~GlobalIndexer() = default; /// Call this immediately after construction when running unit tests. void initialiseTest() {} @@ -123,7 +127,7 @@ public: /// Convert the local index object to a global index which can be /// used in PETSc vectors and matrices. - int getGlobal(const ind_type& ind) const { + virtual int getGlobal(const ind_type& ind) const { return static_cast(std::round(indices[ind])); } @@ -224,7 +228,8 @@ private: Array int_indices; /// The first and last global index on this processor (inclusive in /// both cases) - int globalStart, globalEnd; + int globalStart = 0; + int globalEnd = 0; /// Stencil for which this indexer has been configured OperatorStencil stencils; @@ -237,4 +242,49 @@ private: mutable std::vector numDiagonal, numOffDiagonal; }; +namespace bout { +/// A simplified local-to-global index converter: includes all boundaries, but +/// not guard cells, and assumes a simple rectangular division of processes. +/// +/// The conversion indices are not precalculated or cached, so this may be less +/// efficient than `GlobalIndexer`. However, it can convert "local" indices +/// outside of the range owned by the current process, for example a local index +/// with `x = -1` which would be owned by the process to the left. +class SimpleGlobalIndexer : public GlobalIndexer { + // Stencil to ensure x-boundaries (to a depth of 2) are included in the base + // indexer + // + // Actually, I'm not 100% sure what this does, but this appears to the minimum + // necessary to get things to work + static auto XZstencil(const Mesh& mesh) { + using ind_type = Field3D::ind_type; + const IndexOffset zero; + OperatorStencil stencil; + + stencil.add( + [&mesh](ind_type ind) -> bool { + return (mesh.xstart <= ind.x() and ind.x() <= mesh.xend) + and (mesh.ystart <= ind.y() and ind.y() <= mesh.yend) + and (mesh.zstart <= ind.z() and ind.z() <= mesh.zend); + }, + {zero.xm(2), zero.xp(2), zero.zm(2), zero.zp(2)}); + + return stencil; + } + +public: + SimpleGlobalIndexer() = default; + explicit SimpleGlobalIndexer(Mesh* localmesh, bool autoInitialise = true) + : GlobalIndexer(localmesh, XZstencil(*localmesh), autoInitialise) {} + + int getGlobal(const ind_type& ind) const override { + const auto& mesh = *getMesh(); + return ((mesh.getGlobalXIndex(ind.x()) * mesh.GlobalNy) + + mesh.getGlobalYIndex(ind.y())) + * mesh.GlobalNz + + mesh.getGlobalZIndex(ind.z()); + } +}; +} // namespace bout + #endif // BOUT_GLOBALINDEXER_H diff --git a/include/bout/interpolation_xz.hxx b/include/bout/interpolation_xz.hxx index 52dc38f174..0ea5efa7ac 100644 --- a/include/bout/interpolation_xz.hxx +++ b/include/bout/interpolation_xz.hxx @@ -1,8 +1,7 @@ /************************************************************************** - * Copyright 2010-2020 B.D.Dudson, S.Farley, P. Hill, J.T. Omotani, J.T. Parker, - * M.V.Umansky, X.Q.Xu + * Copyright 2010-2024 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -24,7 +23,12 @@ #ifndef BOUT_INTERP_XZ_H #define BOUT_INTERP_XZ_H +#include "bout/build_defines.hxx" +#include "bout/field3d.hxx" +#include "bout/globalindexer.hxx" #include "bout/mask.hxx" +#include "bout/petsc_interface.hxx" +#include "bout/petsclib.hxx" class Options; @@ -95,20 +99,30 @@ public: ASSERT1(has_region); return getRegion(); } + /// Calculate weights using given offsets in X and Z virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") = 0; - virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") = 0; + void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, + const std::string& region = "RGN_NOBNDRY") { + setMask(mask); + calcWeights(delta_x, delta_z, region); + } + /// Use pre-calculated weights virtual Field3D interpolate(const Field3D& f, const std::string& region = "RGN_NOBNDRY") const = 0; - virtual Field3D interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, - const std::string& region = "RGN_NOBNDRY") = 0; - virtual Field3D interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") = 0; + + /// Calculate weights then interpolate + Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, + const std::string& region = "RGN_NOBNDRY") { + calcWeights(delta_x, delta_z, region); + return interpolate(f, region); + } + Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, + const BoutMask& mask, const std::string& region = "RGN_NOBNDRY") { + calcWeights(delta_x, delta_z, mask, region); + return interpolate(f, region); + } // Interpolate using the field at (x,y+y_offset,z), rather than (x,y,z) void setYOffset(int offset) { y_offset = offset; } @@ -162,18 +176,10 @@ public: void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") override; - void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") override; - // Use precalculated weights Field3D interpolate(const Field3D& f, const std::string& region = "RGN_NOBNDRY") const override; - // Calculate weights and interpolate - Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const std::string& region = "RGN_NOBNDRY") override; - Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") override; + std::vector getWeightsForYApproximation(int i, int j, int k, int yoffset) override; }; @@ -218,21 +224,12 @@ public: void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") override; - void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") override; + + using XZInterpolation::interpolate; // Use precalculated weights Field3D interpolate(const Field3D& f, const std::string& region = "RGN_NOBNDRY") const override; - // Calculate weights and interpolate - Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const std::string& region = "RGN_NOBNDRY") override; - Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") override; - BoutReal lagrange_4pt(BoutReal v2m, BoutReal vm, BoutReal vp, BoutReal v2p, - BoutReal offset) const; - BoutReal lagrange_4pt(const BoutReal v[], BoutReal offset) const; }; class XZBilinear : public XZInterpolation { @@ -251,19 +248,68 @@ public: void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") override; - void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, + + using XZInterpolation::interpolate; + + // Use precalculated weights + Field3D interpolate(const Field3D& f, + const std::string& region = "RGN_NOBNDRY") const override; +}; + +#if BOUT_HAS_PETSC +class PetscXZHermiteSpline : public XZInterpolation { + PetscLib petsclib; + IndexerPtr indexer; + PetscMatrix weights; + + Tensor i_corner; // x-index of bottom-left grid point + Tensor k_corner; // z-index of bottom-left grid point + + // 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; + +public: + PetscXZHermiteSpline(int y_offset = 0, Mesh* mesh = nullptr); + PetscXZHermiteSpline(Mesh* mesh = nullptr) : PetscXZHermiteSpline(0, mesh) {} + PetscXZHermiteSpline(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) + : PetscXZHermiteSpline(y_offset, mesh) { + region = regionFromMask(mask, localmesh); + } + PetscXZHermiteSpline(const PetscXZHermiteSpline& other) = delete; + PetscXZHermiteSpline(PetscXZHermiteSpline&& other) = default; + PetscXZHermiteSpline& operator=(const PetscXZHermiteSpline& other) = delete; + PetscXZHermiteSpline& operator=(PetscXZHermiteSpline&& other) = delete; + ~PetscXZHermiteSpline() override = default; + + void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") override; + void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, + const std::string& region = "RGN_NOBNDRY"); // Use precalculated weights Field3D interpolate(const Field3D& f, const std::string& region = "RGN_NOBNDRY") const override; // Calculate weights and interpolate Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const std::string& region = "RGN_NOBNDRY") override; + const std::string& region = "RGN_NOBNDRY"); Field3D interpolate(const Field3D& f, const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, - const std::string& region = "RGN_NOBNDRY") override; + const BoutMask& mask, const std::string& region = "RGN_NOBNDRY"); + + std::vector + getWeightsForYApproximation(int i, int j, int k, int yoffset) override; }; +#endif // BOUT_HAS_PETSC class XZInterpolationFactory : public Factory { @@ -279,11 +325,30 @@ public: ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const { return Factory::create(type, nullptr); } - - static void ensureRegistered(); }; template using RegisterXZInterpolation = XZInterpolationFactory::RegisterInFactory; +using RegisterUnavailableXZInterpolation = + XZInterpolationFactory::RegisterUnavailableInFactory; + +namespace { +const inline RegisterXZInterpolation registerinterphermitespline{ + "hermitespline"}; +const inline RegisterXZInterpolation + registerinterpmonotonichermitespline{"monotonichermitespline"}; +const inline RegisterXZInterpolation registerinterplagrange4pt{ + "lagrange4pt"}; +const inline RegisterXZInterpolation registerinterpbilinear{"bilinear"}; + +#if BOUT_HAS_PETSC +const inline RegisterXZInterpolation registerpetschermitespline{ + "petschermitespline"}; +#else +const inline RegisterUnavailableXZInterpolation register_unavailable_petsc_hermite_spline{ + "petschermitespline", "BOUT++ was not configured with PETSc"}; +#endif +} // namespace + #endif // BOUT_INTERP_XZ_H diff --git a/src/mesh/interpolation/bilinear_xz.cxx b/src/mesh/interpolation/bilinear_xz.cxx index 8445764a8f..4c0603f983 100644 --- a/src/mesh/interpolation/bilinear_xz.cxx +++ b/src/mesh/interpolation/bilinear_xz.cxx @@ -83,12 +83,6 @@ void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, } } -void XZBilinear::calcWeights(const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, const std::string& region) { - setMask(mask); - calcWeights(delta_x, delta_z, region); -} - Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) const { ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; @@ -113,16 +107,3 @@ Field3D XZBilinear::interpolate(const Field3D& f, const std::string& region) con } return f_interp; } - -Field3D XZBilinear::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const std::string& region) { - calcWeights(delta_x, delta_z, region); - return interpolate(f, region); -} - -Field3D XZBilinear::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const BoutMask& mask, - const std::string& region) { - calcWeights(delta_x, delta_z, mask, region); - return interpolate(f, region); -} diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index a8e5df7cdf..e2ea540a1f 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -111,12 +111,6 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z } } -void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, const std::string& region) { - setMask(mask); - calcWeights(delta_x, delta_z, region); -} - /*! * Return position and weight of points needed to approximate the function value at the * point that the field line through (i,j,k) meets the (j+1)-plane. For the case where @@ -223,16 +217,3 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region } return f_interp; } - -Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const std::string& region) { - calcWeights(delta_x, delta_z, region); - return interpolate(f, region); -} - -Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const BoutMask& mask, - const std::string& region) { - calcWeights(delta_x, delta_z, mask, region); - return interpolate(f, region); -} diff --git a/src/mesh/interpolation/lagrange_4pt_xz.cxx b/src/mesh/interpolation/lagrange_4pt_xz.cxx index 92c14ecfd5..9d37bf501b 100644 --- a/src/mesh/interpolation/lagrange_4pt_xz.cxx +++ b/src/mesh/interpolation/lagrange_4pt_xz.cxx @@ -20,11 +20,29 @@ * **************************************************************************/ +#include "bout/bout_types.hxx" #include "bout/globals.hxx" #include "bout/interpolation_xz.hxx" #include "bout/mesh.hxx" -#include +#include + +namespace { +// 4-point Lagrangian interpolation +// offset must be between 0 and 1 +BoutReal lagrange_4pt(const BoutReal v2m, const BoutReal v1m, const BoutReal v1p, + const BoutReal v2p, const BoutReal offset) { + return -offset * (offset - 1.0) * (offset - 2.0) * v2m / 6.0 + + 0.5 * (offset * offset - 1.0) * (offset - 2.0) * v1m + - 0.5 * offset * (offset + 1.0) * (offset - 2.0) * v1p + + offset * (offset * offset - 1.0) * v2p / 6.0; +} +// Convenience helper +BoutReal lagrange_4pt(const std::array& x, const BoutReal offset) { + return lagrange_4pt(x[0], x[1], x[2], x[3], offset); +} + +} // namespace XZLagrange4pt::XZLagrange4pt(int y_offset, Mesh* mesh) : XZInterpolation(y_offset, mesh), t_x(localmesh), t_z(localmesh) { @@ -75,12 +93,6 @@ void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, } } -void XZLagrange4pt::calcWeights(const Field3D& delta_x, const Field3D& delta_z, - const BoutMask& mask, const std::string& region) { - setMask(mask); - calcWeights(delta_x, delta_z, region); -} - Field3D XZLagrange4pt::interpolate(const Field3D& f, const std::string& region) const { ASSERT1(f.getMesh() == localmesh); @@ -107,55 +119,25 @@ Field3D XZLagrange4pt::interpolate(const Field3D& f, const std::string& region) const int jz2mnew = (jz - 1 + ncz) % ncz; // Interpolate in Z first - BoutReal xvals[4]; - const int y_next = y + y_offset; - xvals[0] = lagrange_4pt(f(jx2mnew, y_next, jz2mnew), f(jx2mnew, y_next, jz), - f(jx2mnew, y_next, jzpnew), f(jx2mnew, y_next, jz2pnew), - t_z(x, y, z)); + const std::array xvals{ + lagrange_4pt(f(jx2mnew, y_next, jz2mnew), f(jx2mnew, y_next, jz), + f(jx2mnew, y_next, jzpnew), f(jx2mnew, y_next, jz2pnew), + t_z(x, y, z)), - xvals[1] = lagrange_4pt(f(jx, y_next, jz2mnew), f(jx, y_next, jz), - f(jx, y_next, jzpnew), f(jx, y_next, jz2pnew), t_z(x, y, z)); + lagrange_4pt(f(jx, y_next, jz2mnew), f(jx, y_next, jz), f(jx, y_next, jzpnew), + f(jx, y_next, jz2pnew), t_z(x, y, z)), - xvals[2] = lagrange_4pt(f(jxpnew, y_next, jz2mnew), f(jxpnew, y_next, jz), - f(jxpnew, y_next, jzpnew), f(jxpnew, y_next, jz2pnew), t_z(x, y, z)); - - xvals[3] = lagrange_4pt(f(jx2pnew, y_next, jz2mnew), f(jx2pnew, y_next, jz), - f(jx2pnew, y_next, jzpnew), f(jx2pnew, y_next, jz2pnew), - t_z(x, y, z)); + f(jxpnew, y_next, jzpnew), f(jxpnew, y_next, jz2pnew), t_z(x, y, z)), + lagrange_4pt(f(jx2pnew, y_next, jz2mnew), f(jx2pnew, y_next, jz), + f(jx2pnew, y_next, jzpnew), f(jx2pnew, y_next, jz2pnew), + t_z(x, y, z)), + }; // Then in X f_interp(x, y_next, z) = lagrange_4pt(xvals, t_x(x, y, z)); } return f_interp; } - -Field3D XZLagrange4pt::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const std::string& region) { - calcWeights(delta_x, delta_z, region); - return interpolate(f, region); -} - -Field3D XZLagrange4pt::interpolate(const Field3D& f, const Field3D& delta_x, - const Field3D& delta_z, const BoutMask& mask, - const std::string& region) { - calcWeights(delta_x, delta_z, mask, region); - return interpolate(f, region); -} - -// 4-point Lagrangian interpolation -// offset must be between 0 and 1 -BoutReal XZLagrange4pt::lagrange_4pt(const BoutReal v2m, const BoutReal vm, - const BoutReal vp, const BoutReal v2p, - const BoutReal offset) const { - return -offset * (offset - 1.0) * (offset - 2.0) * v2m / 6.0 - + 0.5 * (offset * offset - 1.0) * (offset - 2.0) * vm - - 0.5 * offset * (offset + 1.0) * (offset - 2.0) * vp - + offset * (offset * offset - 1.0) * v2p / 6.0; -} - -BoutReal XZLagrange4pt::lagrange_4pt(const BoutReal v[], const BoutReal offset) const { - return lagrange_4pt(v[0], v[1], v[2], v[3], offset); -} diff --git a/src/mesh/interpolation/petsc_hermite_spline_xz.cxx b/src/mesh/interpolation/petsc_hermite_spline_xz.cxx new file mode 100644 index 0000000000..7f5323a39d --- /dev/null +++ b/src/mesh/interpolation/petsc_hermite_spline_xz.cxx @@ -0,0 +1,240 @@ +/************************************************************************** + * Copyright 2024 The BOUT++ Team + * + * Contact: Ben Dudson, bd512@york.ac.uk + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/field3d.hxx" +#include "bout/globalindexer.hxx" +#include "bout/globals.hxx" +#include "bout/index_derivs_interface.hxx" +#include "bout/interpolation_xz.hxx" +#include "bout/operatorstencil.hxx" +#include "bout/paralleltransform.hxx" +#include "bout/petsc_interface.hxx" +#include "bout/petsclib.hxx" +#include "bout/region.hxx" + +#include +#include +#include +#include + +PetscXZHermiteSpline::PetscXZHermiteSpline(int y_offset, Mesh* mesh) + : XZInterpolation(y_offset, mesh), + petsclib(&Options::root()["mesh:paralleltransform:xzinterpolation:hermitespline"]), + indexer(std::make_shared(localmesh)), + weights(indexer, false), h00_x(localmesh), h01_x(localmesh), h10_x(localmesh), + h11_x(localmesh), h00_z(localmesh), h01_z(localmesh), h10_z(localmesh), + h11_z(localmesh) { + + // Index arrays contain guard cells in order to get subscripts right + i_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz); + k_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz); + + // Initialise in order to avoid 'uninitialized value' errors from Valgrind when using + // guard-cell values + k_corner = -1; + + // Allocate Field3D members + h00_x.allocate(); + h01_x.allocate(); + h10_x.allocate(); + h11_x.allocate(); + h00_z.allocate(); + h01_z.allocate(); + h10_z.allocate(); + h11_z.allocate(); + + // The stencil has 16 elements, so maximum number of columns in any + // given row (whether on local or remote process) is 16 + // NOLINTNEXTLINE(misc-include-cleaner) + MatMPIAIJSetPreallocation(*weights.get(), 16, nullptr, 16, nullptr); +} + +void PetscXZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, + const std::string& region) { + + const int ny = localmesh->LocalNy; + const int nz = localmesh->LocalNz; + const int xend = (localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE() + + localmesh->xstart - 1; + + // TODO: work out why using `getRegion(region)` directly causes + // issues on more than one core + const auto actual_region = getRegion(region); + BOUT_FOR(i, actual_region) { + const int x = i.x(); + const int y = i.y(); + const int z = i.z(); + + // The integer part of xt_prime, zt_prime are the indices of the cell + // containing the field line end-point + int i_corn = static_cast(std::floor(delta_x(x, y, z))); + k_corner(x, y, z) = static_cast(std::floor(delta_z(x, y, z))); + + // t_x, t_z are the normalised coordinates \in [0,1) within the cell + // calculated by taking the remainder of the floating point index + BoutReal t_x = delta_x(x, y, z) - static_cast(i_corn); + const BoutReal t_z = delta_z(x, y, z) - static_cast(k_corner(x, y, z)); + + // NOTE: A (small) hack to avoid one-sided differences + if (i_corn >= xend) { + i_corn = xend - 1; + t_x = 1.0; + } + if (i_corn < localmesh->xstart) { + i_corn = localmesh->xstart; + t_x = 0.0; + } + + k_corner(x, y, z) = ((k_corner(x, y, z) % nz) + nz) % nz; + + // Check that t_x and t_z are in range + if ((t_x < 0.0) || (t_x > 1.0)) { + throw BoutException( + "t_x={:e} out of range at ({:d},{:d},{:d}) (delta_x={:e}, i_corn={:d})", t_x, x, + y, z, delta_x(x, y, z), i_corn); + } + + if ((t_z < 0.0) || (t_z > 1.0)) { + throw BoutException( + "t_z={:e} out of range at ({:d},{:d},{:d}) (delta_z={:e}, k_corner={:d})", t_z, + x, y, z, delta_z(x, y, z), k_corner(x, y, z)); + } + + h00_x[i] = (2. * t_x * t_x * t_x) - (3. * t_x * t_x) + 1.; + h00_z[i] = (2. * t_z * t_z * t_z) - (3. * t_z * t_z) + 1.; + + h01_x[i] = (-2. * t_x * t_x * t_x) + (3. * t_x * t_x); + h01_z[i] = (-2. * t_z * t_z * t_z) + (3. * t_z * t_z); + + h10_x[i] = t_x * (1. - t_x) * (1. - t_x); + h10_z[i] = t_z * (1. - t_z) * (1. - t_z); + + h11_x[i] = (t_x * t_x * t_x) - (t_x * t_x); + h11_z[i] = (t_z * t_z * t_z) - (t_z * t_z); + + // 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))), + ny, nz); + + // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; + // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; + // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; + // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; + + // Quite a few duplicated terms, could pre-calculate + weights(i, i_c.xm()) = (h10_x[i] * h11_z[i] / 4) - (h10_x[i] * h00_z[i] / 2); + weights(i, i_c.xm().zm()) = (h10_x[i] * h10_z[i] / 4); + weights(i, i_c.xm().zp()) = -(h10_x[i] * h01_z[i] / 2) - (h10_x[i] * h10_z[i] / 4); + weights(i, i_c.xm().zp(2)) = -(h10_x[i] * h11_z[i] / 4); + weights(i, i_c.xp()) = (h01_x[i] * h00_z[i]) + (h10_x[i] * h00_z[i] / 2) + - (h01_x[i] * h11_z[i] / 2) - (h10_x[i] * h11_z[i] / 4); + weights(i, i_c.xp().zm()) = -(h01_x[i] * h10_z[i] / 2) - (h10_x[i] * h10_z[i] / 4); + weights(i, i_c.xp().zp()) = (h01_x[i] * h01_z[i]) + (h01_x[i] * h10_z[i] / 2) + + (h10_x[i] * h01_z[i] / 2) + (h10_x[i] * h10_z[i] / 4); + weights(i, i_c.xp().zp(2)) = (h01_x[i] * h11_z[i] / 2) + (h10_x[i] * h11_z[i] / 4); + weights(i, i_c.xp(2)) = (h11_x[i] * h00_z[i] / 2) - (h11_x[i] * h11_z[i] / 4); + weights(i, i_c.xp(2).zm()) = -(h11_x[i] * h10_z[i] / 4); + weights(i, i_c.xp(2).zp()) = (h11_x[i] * h01_z[i] / 2) + (h11_x[i] * h10_z[i] / 4); + weights(i, i_c.xp(2).zp(2)) = (h11_x[i] * h11_z[i] / 4); + weights(i, i_c) = (h00_x[i] * h00_z[i]) + (h11_x[i] * h11_z[i] / 4) + - (h00_x[i] * h11_z[i] / 2) - (h11_x[i] * h00_z[i] / 2); + weights(i, i_c.zm()) = (h11_x[i] * h10_z[i] / 4) - (h00_x[i] * h10_z[i] / 2); + weights(i, i_c.zp()) = (h00_x[i] * h01_z[i]) + (h00_x[i] * h10_z[i] / 2) + - (h11_x[i] * h01_z[i] / 2) - (h11_x[i] * h10_z[i] / 4); + weights(i, i_c.zp(2)) = (h00_x[i] * h11_z[i] / 2) - (h11_x[i] * h11_z[i] / 4); + } + + weights.assemble(); +} + +void PetscXZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, + const BoutMask& mask, const std::string& region) { + setMask(mask); + calcWeights(delta_x, delta_z, region); +} + +/*! + * Return position and weight of points needed to approximate the function value at the + * point that the field line through (i,j,k) meets the (j+1)-plane. For the case where + * only the z-direction is not aligned to grid points, the approximation is: f(i,j+1,k*) = + * h00_z * f(i,j+1,k) + h01_z * f(i,j+1,k+1) + * + h10_z * dfdz(i,j+1,k) + h11_z * dfdz(i,j+1,k+1) + * = h00_z * f(i,j+1,k) + h01_z * f(i,j+1,k+1) + * + h10_z * ( f(i,j+1,k+1) - f(i,j+1,k-1) ) / 2 + * + h11_z * ( f(i,j+1,k+2) - f(i,j+1,k) ) / 2 + * for k* a point between k and k+1. Therefore, this function returns + * position weight + * (i, j+1, k-1) - h10_z / 2 + * (i, j+1, k) h00_z - h11_z / 2 + * (i, j+1, k+1) h01_z + h10_z / 2 + * (i, j+1, k+2) h11_z / 2 + */ +std::vector +PetscXZHermiteSpline::getWeightsForYApproximation(int i, int j, int k, int yoffset) { + const int nz = localmesh->LocalNz; + const int k_mod = k_corner(i, j, k); + const int k_mod_m1 = (k_mod > 0) ? (k_mod - 1) : (nz - 1); + const int k_mod_p1 = (k_mod == nz) ? 0 : k_mod + 1; + const int k_mod_p2 = (k_mod_p1 == nz) ? 0 : k_mod_p1 + 1; + + return {{i, j + yoffset, k_mod_m1, -0.5 * h10_z(i, j, k)}, + {i, j + yoffset, k_mod, h00_z(i, j, k) - 0.5 * h11_z(i, j, k)}, + {i, j + yoffset, k_mod_p1, h01_z(i, j, k) + 0.5 * h10_z(i, j, k)}, + {i, j + yoffset, k_mod_p2, 0.5 * h11_z(i, j, k)}}; +} + +Field3D +PetscXZHermiteSpline::interpolate(const Field3D& f, + [[maybe_unused]] const std::string& region) const { + + ASSERT1(f.getMesh() == localmesh); + const PetscVector f_petsc(f, indexer); + const PetscVector 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(); +} + +Field3D PetscXZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, + const Field3D& delta_z, + const std::string& region) { + calcWeights(delta_x, delta_z, region); + return interpolate(f, region); +} + +Field3D PetscXZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, + const Field3D& delta_z, const BoutMask& mask, + const std::string& region) { + calcWeights(delta_x, delta_z, mask, region); + return interpolate(f, region); +} + +#endif diff --git a/src/mesh/interpolation_xz.cxx b/src/mesh/interpolation_xz.cxx index f7f0b457f2..52c348921c 100644 --- a/src/mesh/interpolation_xz.cxx +++ b/src/mesh/interpolation_xz.cxx @@ -83,13 +83,3 @@ const Field3D interpolate(const Field2D& f, const Field3D& delta_x) { } return result; } - -void XZInterpolationFactory::ensureRegistered() {} - -namespace { -RegisterXZInterpolation registerinterphermitespline{"hermitespline"}; -RegisterXZInterpolation registerinterpmonotonichermitespline{ - "monotonichermitespline"}; -RegisterXZInterpolation registerinterplagrange4pt{"lagrange4pt"}; -RegisterXZInterpolation registerinterpbilinear{"bilinear"}; -} // namespace diff --git a/tests/integrated/test-interpolate/data/BOUT.inp b/tests/integrated/test-interpolate/data/BOUT.inp index 101c63f3c7..804e780bbe 100644 --- a/tests/integrated/test-interpolate/data/BOUT.inp +++ b/tests/integrated/test-interpolate/data/BOUT.inp @@ -4,7 +4,6 @@ # MZ = 4 # Z size -NXPE = 1 ZMAX = 1 MXG = 2 diff --git a/tests/integrated/test-interpolate/runtest b/tests/integrated/test-interpolate/runtest index 08975cfd33..77e8d0b99b 100755 --- a/tests/integrated/test-interpolate/runtest +++ b/tests/integrated/test-interpolate/runtest @@ -6,6 +6,7 @@ from boututils.run_wrapper import build_and_log, shell, launch_safe from boutdata import collect +import boutconfig from numpy import sqrt, max, abs, mean, array, log, polyfit from sys import stdout, exit @@ -15,9 +16,6 @@ show_plot = False # List of NX values to use nxlist = [16, 32, 64, 128] -# Only testing 2D (x, z) slices, so only need one processor -nproc = 1 - # Variables to compare varlist = ["a", "b", "c"] markers = ["bo", "r^", "kx"] @@ -29,6 +27,8 @@ methods = { "bilinear": 2, } +if boutconfig.has["petsc"]: + methods["petschermitespline"] = 3 build_and_log("Interpolation test") @@ -37,7 +37,7 @@ success = True for method in methods: print("------------------------------") - print("Using {} interpolation".format(method)) + print(f"Using {method} interpolation") error_2 = {} error_inf = {} @@ -48,25 +48,20 @@ for method in methods: for nx in nxlist: dx = 1.0 / (nx) - args = ( - " mesh:nx={nx4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method}".format( - nx4=nx + 4, dx=dx, nx=nx, method=method - ) - ) - - cmd = "./test_interpolate" + args + nproc = 2 if method == "petschermitespline" else 1 + cmd = f"./test_interpolate mesh:nx={nx + 4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method} NXPE={nproc}" shell("rm data/BOUT.dmp.*.nc") s, out = launch_safe(cmd, nproc=nproc, pipe=True) - with open("run.log.{}.{}".format(method, nx), "w") as f: + with open(f"run.log.{method}.{nx}", "w") as f: f.write(out) # Collect output data for var in varlist: - interp = collect(var + "_interp", path="data", xguards=False, info=False) + interp = collect(f"{var}_interp", path="data", xguards=False, info=False) solution = collect( - var + "_solution", path="data", xguards=False, info=False + f"{var}_solution", path="data", xguards=False, info=False ) E = interp - solution diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index a090877552..30b0899316 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -11,14 +11,19 @@ #include #include "bout/bout.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" #include "bout/constants.hxx" +#include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" #include "bout/interpolation_xz.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" #include "bout/sys/generator_context.hxx" /// Get a FieldGenerator from the options for a variable -std::shared_ptr getGeneratorFromOptions(const std::string& varname, - std::string& func) { +auto getGeneratorFromOptions(const std::string& varname, std::string& func) { Options* options = Options::getRoot()->getSection(varname); options->get("solution", func, "0.0"); @@ -30,88 +35,89 @@ std::shared_ptr getGeneratorFromOptions(const std::string& varna int main(int argc, char** argv) { BoutInitialise(argc, argv); - - // Random number generator - std::default_random_engine generator; - // Uniform distribution of BoutReals from 0 to 1 - std::uniform_real_distribution distribution{0.0, 1.0}; - - using bout::globals::mesh; - - FieldFactory f(mesh); - - // Set up generators and solutions for three different analtyic functions - std::string a_func; - auto a_gen = getGeneratorFromOptions("a", a_func); - Field3D a = f.create3D(a_func); - Field3D a_solution = 0.0; - Field3D a_interp = 0.0; - - std::string b_func; - auto b_gen = getGeneratorFromOptions("b", b_func); - Field3D b = f.create3D(b_func); - Field3D b_solution = 0.0; - Field3D b_interp = 0.0; - - std::string c_func; - auto c_gen = getGeneratorFromOptions("c", c_func); - Field3D c = f.create3D(c_func); - Field3D c_solution = 0.0; - Field3D c_interp = 0.0; - - // x and z displacements - Field3D deltax = 0.0; - Field3D deltaz = 0.0; - - // Bind the random number generator and distribution into a single function - auto dice = std::bind(distribution, generator); - - for (const auto& index : deltax) { - // Get some random displacements - BoutReal dx = index.x() + dice(); - BoutReal dz = index.z() + dice(); - // For the last point, put the displacement inwards - // Otherwise we try to interpolate in the guard cells, which doesn't work so well - if (index.x() >= mesh->xend) { - dx = index.x() - dice(); + { + // Random number generator + const std::default_random_engine generator; + // Uniform distribution of BoutReals from 0 to 1 + const std::uniform_real_distribution distribution{0.0, 1.0}; + + using bout::globals::mesh; + + const FieldFactory fieldfact(mesh); + + // Set up generators and solutions for three different analtyic functions + std::string a_func; + auto a_gen = getGeneratorFromOptions("a", a_func); + const Field3D a_field = fieldfact.create3D(a_func); + Field3D a_solution = 0.0; + Field3D a_interp = 0.0; + + std::string b_func; + auto b_gen = getGeneratorFromOptions("b", b_func); + const Field3D b_field = fieldfact.create3D(b_func); + Field3D b_solution = 0.0; + Field3D b_interp = 0.0; + + std::string c_func; + auto c_gen = getGeneratorFromOptions("c", c_func); + const Field3D c_field = fieldfact.create3D(c_func); + Field3D c_solution = 0.0; + Field3D c_interp = 0.0; + + // x and z displacements + Field3D deltax = 0.0; + Field3D deltaz = 0.0; + + // Bind the random number generator and distribution into a single function + auto dice = std::bind(distribution, generator); + + for (const auto& index : deltax) { + // Get some random displacements + BoutReal dx = index.x() + dice(); + const BoutReal dz = index.z() + dice(); + // For the last point, put the displacement inwards + // Otherwise we try to interpolate in the guard cells, which doesn't work so well + if (index.x() >= mesh->xend && mesh->getNXPE() - 1 == mesh->getXProcIndex()) { + dx = index.x() - dice(); + } + deltax[index] = dx; + deltaz[index] = dz; + // Get the global indices + bout::generator::Context pos{index, CELL_CENTRE, deltax.getMesh(), 0.0}; + pos.set("x", mesh->GlobalX(dx), "z", TWOPI * dz / mesh->LocalNz); + // Generate the analytic solution at the displacements + a_solution[index] = a_gen->generate(pos); + b_solution[index] = b_gen->generate(pos); + c_solution[index] = c_gen->generate(pos); } - deltax[index] = dx; - deltaz[index] = dz; - // Get the global indices - bout::generator::Context pos{index, CELL_CENTRE, deltax.getMesh(), 0.0}; - pos.set("x", mesh->GlobalX(dx), "z", - TWOPI * static_cast(dz) / static_cast(mesh->LocalNz)); - // Generate the analytic solution at the displacements - a_solution[index] = a_gen->generate(pos); - b_solution[index] = b_gen->generate(pos); - c_solution[index] = c_gen->generate(pos); - } - // Create the interpolation object from the input options - auto interp = XZInterpolationFactory::getInstance().create(); + deltax += (mesh->LocalNx - mesh->xstart * 2) * mesh->getXProcIndex(); + // Create the interpolation object from the input options + auto interp = XZInterpolationFactory::getInstance().create(); - // Interpolate the analytic functions at the displacements - a_interp = interp->interpolate(a, deltax, deltaz); - b_interp = interp->interpolate(b, deltax, deltaz); - c_interp = interp->interpolate(c, deltax, deltaz); + // Interpolate the analytic functions at the displacements + a_interp = interp->interpolate(a_field, deltax, deltaz); + b_interp = interp->interpolate(b_field, deltax, deltaz); + c_interp = interp->interpolate(c_field, deltax, deltaz); - Options dump; + Options dump; - dump["a"] = a; - dump["a_interp"] = a_interp; - dump["a_solution"] = a_solution; + dump["a"] = a_field; + dump["a_interp"] = a_interp; + dump["a_solution"] = a_solution; - dump["b"] = b; - dump["b_interp"] = b_interp; - dump["b_solution"] = b_solution; + dump["b"] = b_field; + dump["b_interp"] = b_interp; + dump["b_solution"] = b_solution; - dump["c"] = c; - dump["c_interp"] = c_interp; - dump["c_solution"] = c_solution; + dump["c"] = c_field; + dump["c_interp"] = c_interp; + dump["c_solution"] = c_solution; - bout::writeDefaultOutputFile(dump); + bout::writeDefaultOutputFile(dump); - bout::checkForUnusedOptions(); + bout::checkForUnusedOptions(); + } BoutFinalise(); return 0;