diff --git a/.ci_fedora.sh b/.ci_fedora.sh index b8805abb15..583bb79d7a 100755 --- a/.ci_fedora.sh +++ b/.ci_fedora.sh @@ -41,7 +41,7 @@ then # Ignore weak depencies echo "install_weak_deps=False" >> /etc/dnf/dnf.conf time dnf -y install dnf5 - time dnf5 -y install dnf5-plugins cmake python3-zoidberg python3-natsort + time dnf5 -y install dnf5-plugins cmake python3-zoidberg python3-natsort wget # Allow to override packages - see #2073 time dnf5 copr enable -y davidsch/fixes4bout || : time dnf5 -y upgrade diff --git a/cmake/BOUT++functions.cmake b/cmake/BOUT++functions.cmake index 40e45f99be..77279dfd4b 100644 --- a/cmake/BOUT++functions.cmake +++ b/cmake/BOUT++functions.cmake @@ -162,7 +162,7 @@ endfunction() # function(bout_add_integrated_or_mms_test BUILD_CHECK_TARGET TESTNAME) set(options USE_RUNTEST USE_DATA_BOUT_INP) - set(oneValueArgs EXECUTABLE_NAME PROCESSORS) + set(oneValueArgs EXECUTABLE_NAME PROCESSORS DOWNLOAD DOWNLOAD_NAME) set(multiValueArgs SOURCES EXTRA_FILES REQUIRES CONFLICTS TESTARGS EXTRA_DEPENDS) cmake_parse_arguments(BOUT_TEST_OPTIONS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) @@ -202,6 +202,20 @@ function(bout_add_integrated_or_mms_test BUILD_CHECK_TARGET TESTNAME) add_custom_target(${TESTNAME}) endif() + if (BOUT_TEST_OPTIONS_DOWNLOAD) + if (NOT BOUT_TEST_OPTIONS_DOWNLOAD_NAME) + message(FATAL_ERROR "We need DOWNLOAD_NAME if we should DOWNLOAD!") + endif() + set(output ) + add_custom_command(OUTPUT ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME} + COMMAND wget ${BOUT_TEST_OPTIONS_DOWNLOAD} -O ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMENT "Downloading ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME}" + ) + add_custom_target(download_test_data DEPENDS ${BOUT_TEST_OPTIONS_DOWNLOAD_NAME}) + add_dependencies(${TESTNAME} download_test_data) + endif() + if (BOUT_TEST_OPTIONS_EXTRA_DEPENDS) add_dependencies(${TESTNAME} ${BOUT_TEST_OPTIONS_EXTRA_DEPENDS}) endif() diff --git a/include/bout/interpolation_xz.hxx b/include/bout/interpolation_xz.hxx index 52dc38f174..def8a60a3e 100644 --- a/include/bout/interpolation_xz.hxx +++ b/include/bout/interpolation_xz.hxx @@ -26,6 +26,15 @@ #include "bout/mask.hxx" +#define USE_NEW_WEIGHTS 1 +#if BOUT_HAS_PETSC +#define HS_USE_PETSC 1 +#endif + +#ifdef HS_USE_PETSC +#include "bout/petsclib.hxx" +#endif + class Options; /// Interpolate a field onto a perturbed set of points @@ -43,8 +52,7 @@ public: protected: Mesh* localmesh{nullptr}; - std::string region_name; - std::shared_ptr> region{nullptr}; + int region_id{-1}; public: XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) @@ -52,48 +60,44 @@ public: localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn) {} XZInterpolation(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZInterpolation(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setMask(mask); } XZInterpolation(const std::string& region_name, int y_offset = 0, Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh), region_name(region_name) {} - XZInterpolation(std::shared_ptr> region, int y_offset = 0, - Mesh* mesh = nullptr) - : y_offset(y_offset), localmesh(mesh), region(std::move(region)) {} + : y_offset(y_offset), localmesh(mesh), + region_id(localmesh->getRegionID(region_name)) {} + XZInterpolation(const Region& region, int y_offset = 0, Mesh* mesh = nullptr) + : y_offset(y_offset), localmesh(mesh) { + setRegion(region); + } virtual ~XZInterpolation() = default; - void setMask(const BoutMask& mask) { - region = regionFromMask(mask, localmesh); - region_name = ""; - } + void setMask(const BoutMask& mask) { setRegion(regionFromMask(mask, localmesh)); } void setRegion(const std::string& region_name) { - this->region_name = region_name; - this->region = nullptr; - } - void setRegion(const std::shared_ptr>& region) { - this->region_name = ""; - this->region = region; + this->region_id = localmesh->getRegionID(region_name); } + void setRegion(const std::unique_ptr> region) { setRegion(*region); } void setRegion(const Region& region) { - this->region_name = ""; - this->region = std::make_shared>(region); + std::string name; + int i = 0; + do { + name = fmt::format("unsec_reg_xz_interp_{:d}", i++); + } while (localmesh->hasRegion3D(name)); + localmesh->addRegion(name, region); + this->region_id = localmesh->getRegionID(name); } - Region getRegion() const { - if (!region_name.empty()) { - return localmesh->getRegion(region_name); - } - ASSERT1(region != nullptr); - return *region; + const Region& getRegion() const { + ASSERT2(region_id != -1); + return localmesh->getRegion(region_id); } - Region getRegion(const std::string& region) const { - const bool has_region = !region_name.empty() or this->region != nullptr; - if (!region.empty() and region != "RGN_ALL") { - if (has_region) { - return intersection(localmesh->getRegion(region), getRegion()); - } + const Region& getRegion(const std::string& region) const { + if (region_id == -1) { return localmesh->getRegion(region); } - ASSERT1(has_region); - return getRegion(); + if (region == "" or region == "RGN_ALL") { + return getRegion(); + } + return localmesh->getRegion( + localmesh->getCommonRegion(localmesh->getRegionID(region), region_id)); } virtual void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region = "RGN_NOBNDRY") = 0; @@ -134,8 +138,8 @@ protected: /// This is protected rather than private so that it can be /// extended and used by HermiteSplineMonotonic - Tensor i_corner; // x-index of bottom-left grid point - Tensor k_corner; // z-index of bottom-left grid point + Tensor> i_corner; // 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 @@ -152,12 +156,32 @@ protected: Field3D h10_z; Field3D h11_z; + std::vector newWeights; + +#if HS_USE_PETSC + PetscLib* petsclib; + bool isInit{false}; + Mat petscWeights; + Vec rhs, result; +#endif + public: XZHermiteSpline(Mesh* mesh = 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) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); + } + ~XZHermiteSpline() { +#if HS_USE_PETSC + if (isInit) { + MatDestroy(&petscWeights); + VecDestroy(&rhs); + VecDestroy(&result); + isInit = false; + delete petsclib; + } +#endif } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -188,11 +212,23 @@ public: /// problems most obviously occur. class XZMonotonicHermiteSpline : public XZHermiteSpline { public: - XZMonotonicHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) {} + XZMonotonicHermiteSpline(Mesh* mesh = nullptr) : XZHermiteSpline(0, mesh) { + 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) {} + : 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) {} + : 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. @@ -213,7 +249,7 @@ public: XZLagrange4pt(int y_offset = 0, Mesh* mesh = nullptr); XZLagrange4pt(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZLagrange4pt(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -246,7 +282,7 @@ public: XZBilinear(int y_offset = 0, Mesh* mesh = nullptr); XZBilinear(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZBilinear(y_offset, mesh) { - region = regionFromMask(mask, localmesh); + setRegion(regionFromMask(mask, localmesh)); } void calcWeights(const Field3D& delta_x, const Field3D& delta_z, diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx index b45152fbcc..e965c5d433 100644 --- a/include/bout/utils.hxx +++ b/include/bout/utils.hxx @@ -361,7 +361,6 @@ public: ASSERT2(0 <= i.ind && i.ind < n1 * n2 * n3); return data[i.ind]; } - T& operator[](Ind3D i) { // ny and nz are private :-( // ASSERT2(i.nz == n3); diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 57c6658891..a3cdf0078b 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -147,8 +147,10 @@ shifted``, see :ref:`sec-shifted-metric`), the recommended method is to apply boundary conditions directly to the ``yup`` and ``ydown`` parallel slices. This can be done by setting ``bndry_par_yup`` and ``bndry_par_ydown``, or ``bndry_par_all`` to set both at once. The -possible values are ``parallel_dirichlet``, ``parallel_dirichlet_O3`` -and ``parallel_neumann``. The stencils used are the same as for the +possible values are ``parallel_dirichlet_o1``, +``parallel_dirichlet_o2``, ``parallel_dirichlet_o3`` +and ``parallel_neumann_o1``, ``parallel_neumann_o2``, +``parallel_neumann_o3``. The stencils used are the same as for the standard boundary conditions without the ``parallel_`` prefix, but are applied directly to parallel slices. The boundary condition can only be applied after the parallel slices are calculated, which is usually @@ -168,7 +170,7 @@ For example, for an evolving variable ``f``, put a section in the [f] bndry_xin = dirichlet bndry_xout = dirichlet - bndry_par_all = parallel_neumann + bndry_par_all = parallel_neumann_o2 bndry_ydown = none bndry_yup = none @@ -278,7 +280,7 @@ cells of the base variable. For example, for an evolving variable [f] bndry_xin = dirichlet bndry_xout = dirichlet - bndry_par_all = parallel_dirichlet + bndry_par_all = parallel_dirichlet_o2 bndry_ydown = none bndry_yup = none @@ -289,7 +291,7 @@ communication, while the perpendicular ones before: f.applyBoundary(); mesh->communicate(f); - f.applyParallelBoundary("parallel_neumann"); + f.applyParallelBoundary("parallel_neumann_o2"); Note that during grid generation care has to be taken to ensure that there are no "short" connection lengths. Otherwise it can happen that for a point on a diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index a8e5df7cdf..c0040d096e 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -20,13 +20,87 @@ * **************************************************************************/ +#include "../impls/bout/boutmesh.hxx" #include "bout/globals.hxx" #include "bout/index_derivs_interface.hxx" #include "bout/interpolation_xz.hxx" -#include "bout/mesh.hxx" #include +class IndConverter { +public: + IndConverter(Mesh* mesh) + : mesh(dynamic_cast(mesh)), nxpe(mesh->getNXPE()), nype(mesh->getNYPE()), + xstart(mesh->xstart), ystart(mesh->ystart), zstart(0), + lnx(mesh->LocalNx - 2 * xstart), lny(mesh->LocalNy - 2 * ystart), + lnz(mesh->LocalNz - 2 * zstart) {} + // ix and iy are global indices + // iy is local + int fromMeshToGlobal(int ix, int iy, int iz) { + const int xstart = mesh->xstart; + const int lnx = mesh->LocalNx - xstart * 2; + // x-proc-id + int pex = divToNeg(ix - xstart, lnx); + if (pex < 0) { + pex = 0; + } + if (pex >= nxpe) { + pex = nxpe - 1; + } + const int zstart = 0; + const int lnz = mesh->LocalNz - zstart * 2; + // z-proc-id + // pez only for wrapping around ; later needs similar treatment than pey + const int pez = divToNeg(iz - zstart, lnz); + // y proc-id - y is already local + const int ystart = mesh->ystart; + const int lny = mesh->LocalNy - ystart * 2; + const int pey_offset = divToNeg(iy - ystart, lny); + int pey = pey_offset + mesh->getYProcIndex(); + while (pey < 0) { + pey += nype; + } + while (pey >= nype) { + pey -= nype; + } + ASSERT2(pex >= 0); + ASSERT2(pex < nxpe); + ASSERT2(pey >= 0); + ASSERT2(pey < nype); + return fromLocalToGlobal(ix - pex * lnx, iy - pey_offset * lny, iz - pez * lnz, pex, + pey, 0); + } + int fromLocalToGlobal(const int ilocalx, const int ilocaly, const int ilocalz) { + return fromLocalToGlobal(ilocalx, ilocaly, ilocalz, mesh->getXProcIndex(), + mesh->getYProcIndex(), 0); + } + int fromLocalToGlobal(const int ilocalx, const int ilocaly, const int ilocalz, + const int pex, const int pey, const int pez) { + ASSERT3(ilocalx >= 0); + ASSERT3(ilocaly >= 0); + ASSERT3(ilocalz >= 0); + const int ilocal = ((ilocalx * mesh->LocalNy) + ilocaly) * mesh->LocalNz + ilocalz; + const int ret = ilocal + + mesh->LocalNx * mesh->LocalNy * mesh->LocalNz + * ((pey * nxpe + pex) * nzpe + pez); + ASSERT3(ret >= 0); + ASSERT3(ret < nxpe * nype * mesh->LocalNx * mesh->LocalNy * mesh->LocalNz); + return ret; + } + +private: + // number of procs + BoutMesh* mesh; + const int nxpe; + const int nype; + const int nzpe{1}; + const int xstart, ystart, zstart; + const int lnx, lny, lnz; + static int divToNeg(const int n, const int d) { + return (n < 0) ? ((n - d + 1) / d) : (n / d); + } +}; + XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* mesh) : XZInterpolation(y_offset, mesh), h00_x(localmesh), h01_x(localmesh), h10_x(localmesh), h11_x(localmesh), h00_z(localmesh), h01_z(localmesh), @@ -38,7 +112,6 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* mesh) // Initialise in order to avoid 'uninitialized value' errors from Valgrind when using // guard-cell values - i_corner = -1; k_corner = -1; // Allocate Field3D members @@ -50,45 +123,77 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* mesh) h01_z.allocate(); h10_z.allocate(); h11_z.allocate(); + +#if USE_NEW_WEIGHTS + newWeights.reserve(16); + for (int w = 0; w < 16; ++w) { + newWeights.emplace_back(localmesh); + newWeights[w].allocate(); + } +#ifdef HS_USE_PETSC + petsclib = new PetscLib( + &Options::root()["mesh:paralleltransform:xzinterpolation:hermitespline"]); + // MatCreate(MPI_Comm comm,Mat *A) + // MatCreate(MPI_COMM_WORLD, &petscWeights); + // MatSetSizes(petscWeights, m, m, M, M); + // PetscErrorCode MatCreateAIJ(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, + // PetscInt N, PetscInt d_nz, const PetscInt d_nnz[], + // PetscInt o_nz, const PetscInt o_nnz[], Mat *A) + // MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) + const int m = localmesh->LocalNx * localmesh->LocalNy * localmesh->LocalNz; + const int M = m * localmesh->getNXPE() * localmesh->getNYPE(); + MatCreateAIJ(MPI_COMM_WORLD, m, m, M, M, 16, nullptr, 16, nullptr, &petscWeights); +#endif +#endif +#ifndef HS_USE_PETSC + if (localmesh->getNXPE() > 1) { + throw BoutException("Require PETSc for MPI splitting in X"); + } +#endif } void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, const std::string& region) { - const int ncz = localmesh->LocalNz; - const auto curregion{getRegion(region)}; - BOUT_FOR(i, curregion) { + const int ny = localmesh->LocalNy; + const int nz = localmesh->LocalNz; + const int xend = (localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE() + + localmesh->xstart - 1; +#ifdef HS_USE_PETSC + IndConverter conv{localmesh}; +#endif + BOUT_FOR(i, getRegion(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 - i_corner(x, y, z) = static_cast(floor(delta_x(x, y, z))); + int i_corn = static_cast(floor(delta_x(x, y, z))); k_corner(x, y, z) = static_cast(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_corner(x, y, z)); + BoutReal t_x = delta_x(x, y, z) - static_cast(i_corn); 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_corner(x, y, z) >= localmesh->xend) { - i_corner(x, y, z) = localmesh->xend - 1; + if (i_corn >= xend) { + i_corn = xend - 1; t_x = 1.0; } - if (i_corner(x, y, z) < localmesh->xstart) { - i_corner(x, y, z) = localmesh->xstart; + if (i_corn < localmesh->xstart) { + i_corn = localmesh->xstart; t_x = 0.0; } - k_corner(x, y, z) = ((k_corner(x, y, z) % ncz) + ncz) % ncz; + 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_corner={:d})", t_x, - x, y, z, delta_x(x, y, z), i_corner(x, y, z)); + "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)) { @@ -97,18 +202,107 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z x, y, z, delta_z(x, y, z), k_corner(x, y, z)); } - h00_x(x, y, z) = (2. * t_x * t_x * t_x) - (3. * t_x * t_x) + 1.; - h00_z(x, y, z) = (2. * t_z * t_z * t_z) - (3. * t_z * t_z) + 1.; + i_corner[i] = SpecificInd( + (((i_corn * ny) + (y + y_offset)) * nz + k_corner(x, y, z)), ny, nz); + + 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); - h01_x(x, y, z) = (-2. * t_x * t_x * t_x) + (3. * t_x * t_x); - h01_z(x, y, z) = (-2. * t_z * t_z * t_z) + (3. * t_z * 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); - h10_x(x, y, z) = t_x * (1. - t_x) * (1. - t_x); - h10_z(x, y, z) = t_z * (1. - t_z) * (1. - t_z); +#if USE_NEW_WEIGHTS - h11_x(x, y, z) = (t_x * t_x * t_x) - (t_x * t_x); - h11_z(x, y, z) = (t_z * t_z * t_z) - (t_z * t_z); + for (int w = 0; w < 16; ++w) { + newWeights[w][i] = 0; + } + // The distribution of our weights: + // 0 4 8 12 + // 1 5 9 13 + // 2 6 10 14 + // 3 7 11 15 + // e.g. 1 == ic.xm(); 4 == ic.zm(); 5 == ic; 7 == ic.zp(2); + + // f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; + newWeights[5][i] += h00_x[i] * h00_z[i]; + newWeights[9][i] += h01_x[i] * h00_z[i]; + newWeights[9][i] += h10_x[i] * h00_z[i] / 2; + newWeights[1][i] -= h10_x[i] * h00_z[i] / 2; + newWeights[13][i] += h11_x[i] * h00_z[i] / 2; + newWeights[5][i] -= h11_x[i] * h00_z[i] / 2; + + // f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + + // fx[iczp] * h10_x[i] + fx[icxpzp] * h11_x[i]; + newWeights[6][i] += h00_x[i] * h01_z[i]; + newWeights[10][i] += h01_x[i] * h01_z[i]; + newWeights[10][i] += h10_x[i] * h01_z[i] / 2; + newWeights[2][i] -= h10_x[i] * h01_z[i] / 2; + newWeights[14][i] += h11_x[i] * h01_z[i] / 2; + newWeights[6][i] -= h11_x[i] * h01_z[i] / 2; + + // fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + + // fxz[ic] * h10_x[i]+ fxz[icxp] * h11_x[i]; + newWeights[6][i] += h00_x[i] * h10_z[i] / 2; + newWeights[4][i] -= h00_x[i] * h10_z[i] / 2; + newWeights[10][i] += h01_x[i] * h10_z[i] / 2; + newWeights[8][i] -= h01_x[i] * h10_z[i] / 2; + newWeights[10][i] += h10_x[i] * h10_z[i] / 4; + newWeights[8][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[2][i] -= h10_x[i] * h10_z[i] / 4; + newWeights[0][i] += h10_x[i] * h10_z[i] / 4; + newWeights[14][i] += h11_x[i] * h10_z[i] / 4; + newWeights[12][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[6][i] -= h11_x[i] * h10_z[i] / 4; + newWeights[4][i] += h11_x[i] * h10_z[i] / 4; + + // fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + // fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; + newWeights[7][i] += h00_x[i] * h11_z[i] / 2; + newWeights[5][i] -= h00_x[i] * h11_z[i] / 2; + newWeights[11][i] += h01_x[i] * h11_z[i] / 2; + newWeights[9][i] -= h01_x[i] * h11_z[i] / 2; + newWeights[11][i] += h10_x[i] * h11_z[i] / 4; + newWeights[9][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[3][i] -= h10_x[i] * h11_z[i] / 4; + newWeights[1][i] += h10_x[i] * h11_z[i] / 4; + newWeights[15][i] += h11_x[i] * h11_z[i] / 4; + newWeights[13][i] -= h11_x[i] * h11_z[i] / 4; + newWeights[7][i] -= h11_x[i] * h11_z[i] / 4; + newWeights[5][i] += h11_x[i] * h11_z[i] / 4; +#ifdef HS_USE_PETSC + PetscInt idxn[1] = {conv.fromLocalToGlobal(x, y + y_offset, z)}; + // output.write("debug: {:d} -> {:d}: {:d}:{:d} -> {:d}:{:d}\n", + // conv.fromLocalToGlobal(x, y + y_offset, z), + // conv.fromMeshToGlobal(i_corn, y + y_offset, k_corner(x, y, z)), + // x, z, i_corn, k_corner(x, y, z)); + // ixstep = mesh->LocalNx * mesh->LocalNz; + for (int j = 0; j < 4; ++j) { + PetscInt idxm[4]; + PetscScalar vals[4]; + for (int k = 0; k < 4; ++k) { + idxm[k] = conv.fromMeshToGlobal(i_corn - 1 + j, y + y_offset, + k_corner(x, y, z) - 1 + k); + vals[k] = newWeights[j * 4 + k][i]; + } + MatSetValues(petscWeights, 1, idxn, 4, idxm, vals, INSERT_VALUES); + } +#endif +#endif } +#ifdef HS_USE_PETSC + MatAssemblyBegin(petscWeights, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(petscWeights, MAT_FINAL_ASSEMBLY); + if (!isInit) { + MatCreateVecs(petscWeights, &rhs, &result); + } + isInit = true; +#endif } void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z, @@ -135,11 +329,11 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z */ std::vector XZHermiteSpline::getWeightsForYApproximation(int i, int j, int k, int yoffset) { - const int ncz = localmesh->LocalNz; + 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) : (ncz - 1); - const int k_mod_p1 = (k_mod + 1) % ncz; - const int k_mod_p2 = (k_mod + 2) % ncz; + 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)}, @@ -152,76 +346,80 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region ASSERT1(f.getMesh() == localmesh); Field3D f_interp{emptyFrom(f)}; - // Derivatives are used for tension and need to be on dimensionless - // coordinates - Field3D fx = bout::derivatives::index::DDX(f, CELL_DEFAULT, "DEFAULT"); - localmesh->communicateXZ(fx); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fx); - localmesh->wait(h); - } - Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", "RGN_ALL"); - localmesh->communicateXZ(fz); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fz); - localmesh->wait(h); +#if USE_NEW_WEIGHTS +#ifdef HS_USE_PETSC + BoutReal* ptr; + const BoutReal* cptr; + VecGetArray(rhs, &ptr); + BOUT_FOR(i, f.getRegion("RGN_NOY")) { ptr[int(i)] = f[i]; } + VecRestoreArray(rhs, &ptr); + MatMult(petscWeights, rhs, result); + VecGetArrayRead(result, &cptr); + const auto region2 = y_offset == 0 ? region : fmt::format("RGN_YPAR_{:+d}", y_offset); + BOUT_FOR(i, f.getRegion(region2)) { + f_interp[i] = cptr[int(i)]; + ASSERT2(std::isfinite(cptr[int(i)])); } - Field3D fxz = bout::derivatives::index::DDX(fz, CELL_DEFAULT, "DEFAULT"); - localmesh->communicateXZ(fxz); - // communicate in y, but do not calculate parallel slices - { - auto h = localmesh->sendY(fxz); - localmesh->wait(h); + VecRestoreArrayRead(result, &cptr); +#else + BOUT_FOR(i, getRegion(region)) { + auto ic = i_corner[i]; + auto iyp = i.yp(y_offset); + + f_interp[iyp] = 0; + for (int w = 0; w < 4; ++w) { + f_interp[iyp] += newWeights[w * 4 + 0][i] * f[ic.zm().xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 1][i] * f[ic.xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 2][i] * f[ic.zp().xp(w - 1)]; + f_interp[iyp] += newWeights[w * 4 + 3][i] * f[ic.zp(2).xp(w - 1)]; + } } - - const auto curregion{getRegion(region)}; - BOUT_FOR(i, curregion) { - const int x = i.x(); - const int y = i.y(); - const int z = i.z(); - - // Due to lack of guard cells in z-direction, we need to ensure z-index - // wraps around - const int z_mod = k_corner(x, y, z); - const int z_mod_p1 = (z_mod + 1) % localmesh->LocalNz; - - const int y_next = y + y_offset; +#endif + return f_interp; +#else + // Derivatives are used for tension and need to be on dimensionless + // coordinates + const auto region2 = fmt::format("RGN_YPAR_{:+d}", y_offset); + // f has been communcated, and thus we can assume that the x-boundaries are + // also valid in the y-boundary. Thus the differentiated field needs no + // extra comms. + Field3D fx = bout::derivatives::index::DDX(f, CELL_DEFAULT, "DEFAULT", region2); + Field3D fz = bout::derivatives::index::DDZ(f, CELL_DEFAULT, "DEFAULT", region2); + Field3D fxz = bout::derivatives::index::DDZ(fx, CELL_DEFAULT, "DEFAULT", region2); + + BOUT_FOR(i, getRegion(region)) { + const auto iyp = i.yp(y_offset); + + const auto ic = i_corner[i]; + const auto iczp = ic.zp(); + const auto icxp = ic.xp(); + const auto icxpzp = iczp.xp(); // Interpolate f in X at Z - const BoutReal f_z = f(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal f_z = + f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; // Interpolate f in X at Z+1 - const BoutReal f_zp1 = f(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal f_zp1 = f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + fx[iczp] * h10_x[i] + + fx[icxpzp] * h11_x[i]; // Interpolate fz in X at Z - const BoutReal fz_z = fz(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal fz_z = fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + fxz[ic] * h10_x[i] + + fxz[icxp] * h11_x[i]; // Interpolate fz in X at Z+1 - const BoutReal fz_zp1 = - fz(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal fz_zp1 = fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; // Interpolate in Z - f_interp(x, y_next, z) = +f_z * h00_z(x, y, z) + f_zp1 * h01_z(x, y, z) - + fz_z * h10_z(x, y, z) + fz_zp1 * h11_z(x, y, z); + f_interp[iyp] = + +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; - ASSERT2(std::isfinite(f_interp(x, y_next, z)) || x < localmesh->xstart - || x > localmesh->xend); + ASSERT2(std::isfinite(f_interp[iyp]) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); } return f_interp; +#endif } Field3D XZHermiteSpline::interpolate(const Field3D& f, const Field3D& delta_x, diff --git a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx index abedb27733..4b84bcd265 100644 --- a/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx +++ b/src/mesh/interpolation/monotonic_hermite_spline_xz.cxx @@ -59,48 +59,35 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, const auto curregion{getRegion(region)}; BOUT_FOR(i, curregion) { - const int x = i.x(); - const int y = i.y(); - const int z = i.z(); + const auto iyp = i.yp(y_offset); - // Due to lack of guard cells in z-direction, we need to ensure z-index - // wraps around - const int ncz = localmesh->LocalNz; - const int z_mod = ((k_corner(x, y, z) % ncz) + ncz) % ncz; - const int z_mod_p1 = (z_mod + 1) % ncz; - - const int y_next = y + y_offset; + const auto ic = i_corner[i]; + const auto iczp = ic.zp(); + const auto icxp = ic.xp(); + const auto icxpzp = iczp.xp(); // Interpolate f in X at Z - const BoutReal f_z = f(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal f_z = + f[ic] * h00_x[i] + f[icxp] * h01_x[i] + fx[ic] * h10_x[i] + fx[icxp] * h11_x[i]; // Interpolate f in X at Z+1 - const BoutReal f_zp1 = f(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + f(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fx(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fx(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal f_zp1 = f[iczp] * h00_x[i] + f[icxpzp] * h01_x[i] + fx[iczp] * h10_x[i] + + fx[icxpzp] * h11_x[i]; // Interpolate fz in X at Z - const BoutReal fz_z = fz(i_corner(x, y, z), y_next, z_mod) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod) * h11_x(x, y, z); + const BoutReal fz_z = fz[ic] * h00_x[i] + fz[icxp] * h01_x[i] + fxz[ic] * h10_x[i] + + fxz[icxp] * h11_x[i]; // Interpolate fz in X at Z+1 - const BoutReal fz_zp1 = - fz(i_corner(x, y, z), y_next, z_mod_p1) * h00_x(x, y, z) - + fz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h01_x(x, y, z) - + fxz(i_corner(x, y, z), y_next, z_mod_p1) * h10_x(x, y, z) - + fxz(i_corner(x, y, z) + 1, y_next, z_mod_p1) * h11_x(x, y, z); + const BoutReal fz_zp1 = fz[iczp] * h00_x[i] + fz[icxpzp] * h01_x[i] + + fxz[iczp] * h10_x[i] + fxz[icxpzp] * h11_x[i]; // Interpolate in Z - BoutReal result = +f_z * h00_z(x, y, z) + f_zp1 * h01_z(x, y, z) - + fz_z * h10_z(x, y, z) + fz_zp1 * h11_z(x, y, z); + BoutReal result = + +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; - ASSERT2(std::isfinite(result) || x < localmesh->xstart || x > localmesh->xend); + ASSERT2(std::isfinite(result) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); // Monotonicity // Force the interpolated result to be in the range of the @@ -108,18 +95,14 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, // but also degrades accuracy near maxima and minima. // Perhaps should only impose near boundaries, since that is where // problems most obviously occur. - const BoutReal localmax = BOUTMAX(f(i_corner(x, y, z), y_next, z_mod), - f(i_corner(x, y, z) + 1, y_next, z_mod), - f(i_corner(x, y, z), y_next, z_mod_p1), - f(i_corner(x, y, z) + 1, y_next, z_mod_p1)); + const BoutReal localmax = BOUTMAX(f[ic], f[icxp], f[iczp], f[icxpzp]); - const BoutReal localmin = BOUTMIN(f(i_corner(x, y, z), y_next, z_mod), - f(i_corner(x, y, z) + 1, y_next, z_mod), - f(i_corner(x, y, z), y_next, z_mod_p1), - f(i_corner(x, y, z) + 1, y_next, z_mod_p1)); + const BoutReal localmin = BOUTMIN(f[ic], f[icxp], f[iczp], f[icxpzp]); - ASSERT2(std::isfinite(localmax) || x < localmesh->xstart || x > localmesh->xend); - ASSERT2(std::isfinite(localmin) || x < localmesh->xstart || x > localmesh->xend); + ASSERT2(std::isfinite(localmax) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); + ASSERT2(std::isfinite(localmin) || i.x() < localmesh->xstart + || i.x() > localmesh->xend); if (result > localmax) { result = localmax; @@ -128,7 +111,7 @@ Field3D XZMonotonicHermiteSpline::interpolate(const Field3D& f, result = localmin; } - f_interp(x, y_next, z) = result; + f_interp[iyp] = result; } return f_interp; } diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index cb8c19bbd7..3363d331e1 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -160,6 +160,8 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& const int ncz = map_mesh.LocalNz; BoutMask to_remove(map_mesh); + const int xend = + map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1; // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { @@ -173,7 +175,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& } } - if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= map_mesh.xend)) { + if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= xend)) { // Not a boundary continue; } diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 802f70f99b..b4825c6207 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -21,4 +21,4 @@ y_periodic = true z_periodic = true [mesh:paralleltransform:xzinterpolation] -type = lagrange4pt +type = hermitespline diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 204a9cc271..b51c311a50 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -21,7 +21,7 @@ from sys import stdout import zoidberg as zb -nx = 3 # Not changed for these tests +nx = 4 # Not changed for these tests # Resolution in y and z nlist = [8, 16, 32, 64, 128] @@ -49,92 +49,115 @@ failures = [] build_and_log("FCI MMS test") for nslice in nslices: - error_2[nslice] = [] - error_inf[nslice] = [] - - # Which central difference scheme to use and its expected order - order = nslice * 2 - method_orders[nslice] = {"name": "C{}".format(order), "order": order} - - for n in nlist: - # Define the magnetic field using new poloidal gridding method - # Note that the Bz and Bzprime parameters here must be the same as in mms.py - field = zb.field.Slab(Bz=0.05, Bzprime=0.1) - # Create rectangular poloidal grids - poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid(nx, n, 0.1, 1.0) - # Set the ylength and y locations - ylength = 10.0 - - if yperiodic: - ycoords = linspace(0.0, ylength, n, endpoint=False) + for method in [ + "hermitespline", + "lagrange4pt", + "bilinear", + # "monotonichermitespline", + ]: + error_2[nslice] = [] + error_inf[nslice] = [] + + # Which central difference scheme to use and its expected order + order = nslice * 2 + method_orders[nslice] = {"name": "C{}".format(order), "order": order} + + for n in nlist: + # Define the magnetic field using new poloidal gridding method + # Note that the Bz and Bzprime parameters here must be the same as in mms.py + field = zb.field.Slab(Bz=0.05, Bzprime=0.1) + # Create rectangular poloidal grids + poloidal_grid = zb.poloidal_grid.RectangularPoloidalGrid( + nx, n, 0.1, 1.0, MXG=1 + ) + # Set the ylength and y locations + ylength = 10.0 + + if yperiodic: + ycoords = linspace(0.0, ylength, n, endpoint=False) + else: + # Doesn't include the end points + ycoords = (arange(n) + 0.5) * ylength / float(n) + + # Create the grid + grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) + # Make and write maps + maps = zb.make_maps(grid, field, nslice=nslice, quiet=True, MXG=1) + zb.write_maps( + grid, + field, + maps, + new_names=False, + metric2d=conf.isMetric2D(), + quiet=True, + ) + + args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={} NXPE={}".format( + n, + nslice, + yperiodic, + method_orders[nslice]["name"], + 2 if conf.has["petsc"] and method == "hermitespline" else 1, + ) + args += f" mesh:paralleltransform:xzinterpolation:type={method}" + + # Command to run + cmd = "./fci_mms " + args + + print("Running command: " + cmd) + + # Launch using MPI + s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) + + # Save output to log file + with open("run.log." + str(n), "w") as f: + f.write(out) + + if s: + print("Run failed!\nOutput was:\n") + print(out) + exit(s) + + # Collect data + l_2 = collect( + "l_2", + tind=[1, 1], + info=False, + path=directory, + xguards=False, + yguards=False, + ) + l_inf = collect( + "l_inf", + tind=[1, 1], + info=False, + path=directory, + xguards=False, + yguards=False, + ) + + error_2[nslice].append(l_2) + error_inf[nslice].append(l_inf) + + print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) + + dx = 1.0 / array(nlist) + + # Calculate convergence order + fit = polyfit(log(dx), log(error_2[nslice]), 1) + order = fit[0] + stdout.write("Convergence order = {:f} (fit)".format(order)) + + order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) + stdout.write(", {:f} (small spacing)".format(order)) + + # Should be close to the expected order + if order > method_orders[nslice]["order"] * 0.95: + print("............ PASS\n") else: - # Doesn't include the end points - ycoords = (arange(n) + 0.5) * ylength / float(n) - - # Create the grid - grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=yperiodic) - # Make and write maps - maps = zb.make_maps(grid, field, nslice=nslice, quiet=True) - zb.write_maps( - grid, field, maps, new_names=False, metric2d=conf.isMetric2D(), quiet=True - ) - - args = " MZ={} MYG={} mesh:paralleltransform:y_periodic={} mesh:ddy:first={}".format( - n, nslice, yperiodic, method_orders[nslice]["name"] - ) - - # Command to run - cmd = "./fci_mms " + args - - print("Running command: " + cmd) - - # Launch using MPI - s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) - - # Save output to log file - with open("run.log." + str(n), "w") as f: - f.write(out) - - if s: - print("Run failed!\nOutput was:\n") - print(out) - exit(s) - - # Collect data - l_2 = collect( - "l_2", tind=[1, 1], info=False, path=directory, xguards=False, yguards=False - ) - l_inf = collect( - "l_inf", - tind=[1, 1], - info=False, - path=directory, - xguards=False, - yguards=False, - ) - - error_2[nslice].append(l_2) - error_inf[nslice].append(l_inf) - - print("Errors : l-2 {:f} l-inf {:f}".format(l_2, l_inf)) - - dx = 1.0 / array(nlist) - - # Calculate convergence order - fit = polyfit(log(dx), log(error_2[nslice]), 1) - order = fit[0] - stdout.write("Convergence order = {:f} (fit)".format(order)) - - order = log(error_2[nslice][-2] / error_2[nslice][-1]) / log(dx[-2] / dx[-1]) - stdout.write(", {:f} (small spacing)".format(order)) - - # Should be close to the expected order - if order > method_orders[nslice]["order"] * 0.95: - print("............ PASS\n") - else: - print("............ FAIL\n") - success = False - failures.append(method_orders[nslice]["name"]) + print("............ FAIL\n") + success = False + failures.append(method_orders[nslice]["name"]) with open("fci_mms.pkl", "wb") as output: diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index ef173db7df..e11403efb9 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -10,6 +10,7 @@ add_subdirectory(test-delp2) add_subdirectory(test-datafilefacade) add_subdirectory(test-drift-instability) add_subdirectory(test-drift-instability-staggered) +add_subdirectory(test-fci-mpi) add_subdirectory(test-fieldgroupComm) add_subdirectory(test-fci-boundary) add_subdirectory(test-griddata) diff --git a/tests/integrated/test-fci-mpi/CMakeLists.txt b/tests/integrated/test-fci-mpi/CMakeLists.txt new file mode 100644 index 0000000000..0dd38487a3 --- /dev/null +++ b/tests/integrated/test-fci-mpi/CMakeLists.txt @@ -0,0 +1,9 @@ +bout_add_mms_test(test-fci-mpi + SOURCES fci_mpi.cxx + USE_RUNTEST + USE_DATA_BOUT_INP + PROCESSORS 6 + DOWNLOAD https://zenodo.org/record/7614499/files/W7X-conf4-36x8x128.fci.nc?download=1 + DOWNLOAD_NAME grid.fci.nc + REQUIRES BOUT_HAS_PETSC +) diff --git a/tests/integrated/test-fci-mpi/data/BOUT.inp b/tests/integrated/test-fci-mpi/data/BOUT.inp new file mode 100644 index 0000000000..47272dab61 --- /dev/null +++ b/tests/integrated/test-fci-mpi/data/BOUT.inp @@ -0,0 +1,28 @@ +grid = grid.fci.nc + +[mesh] +symmetricglobalx = true + +[mesh:ddy] +first = C2 +second = C2 + +[mesh:paralleltransform] +type = fci +y_periodic = true +z_periodic = true + +[mesh:paralleltransform:xzinterpolation] +type = hermitespline + +[input_0] +function = sin(z) + +[input_1] +function = cos(y) + +[input_2] +function = sin(x) + +[input_3] +function = sin(x) * sin(z) * cos(y) diff --git a/tests/integrated/test-fci-mpi/fci_mpi.cxx b/tests/integrated/test-fci-mpi/fci_mpi.cxx new file mode 100644 index 0000000000..94520dd4a6 --- /dev/null +++ b/tests/integrated/test-fci-mpi/fci_mpi.cxx @@ -0,0 +1,38 @@ +#include "bout/bout.hxx" +#include "bout/derivs.hxx" +#include "bout/field_factory.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + { + using bout::globals::mesh; + Options* options = Options::getRoot(); + int i = 0; + const std::string default_str{"not_set"}; + Options dump; + while (true) { + std::string temp_str; + options->get(fmt::format("input_{:d}:function", i), temp_str, default_str); + if (temp_str == default_str) { + break; + } + Field3D input{FieldFactory::get()->create3D(fmt::format("input_{:d}:function", i), + Options::getRoot(), mesh)}; + // options->get(fmt::format("input_{:d}:boundary_perp", i), temp_str, s"free_o3"); + mesh->communicate(input); + input.applyParallelBoundary("parallel_neumann_o2"); + for (int slice = -mesh->ystart; slice <= mesh->ystart; ++slice) { + if (slice != 0) { + Field3D tmp{0.}; + BOUT_FOR(i, tmp.getRegion("RGN_NOBNDRY")) { + tmp[i] = input.ynext(slice)[i.yp(slice)]; + } + dump[fmt::format("output_{:d}_{:+d}", i, slice)] = tmp; + } + } + ++i; + } + bout::writeDefaultOutputFile(dump); + } + BoutFinalise(); +} diff --git a/tests/integrated/test-fci-mpi/runtest b/tests/integrated/test-fci-mpi/runtest new file mode 100755 index 0000000000..6676f8f7a5 --- /dev/null +++ b/tests/integrated/test-fci-mpi/runtest @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +# Cores: 8 +# requires: metric_3d + +from boututils.run_wrapper import build_and_log, launch_safe, shell_safe +from boutdata.collect import collect +import boutconfig as conf +import itertools + +import numpy as np + +# Resolution in x and y +nlist = [1, 2, 4] + +maxcores = 8 + +nslices = [1] + +success = True + +build_and_log("FCI MMS test") + +for nslice in nslices: + for NXPE, NYPE in itertools.product(nlist, nlist): + if NXPE * NYPE > maxcores: + continue + + args = f"NXPE={NXPE} NYPE={NYPE}" + # Command to run + cmd = f"./fci_mpi {args}" + + print(f"Running command: {cmd}") + + mthread = maxcores // (NXPE * NYPE) + # Launch using MPI + _, out = launch_safe(cmd, nproc=NXPE * NYPE, mthread=mthread, pipe=True) + + # Save output to log file + with open(f"run.log.{NXPE}.{NYPE}.{nslice}.log", "w") as f: + f.write(out) + + collect_kw = dict(info=False, xguards=False, yguards=False, path="data") + if NXPE == NYPE == 1: + # reference data! + ref = {} + for i in range(4): + for yp in range(1, nslice + 1): + for y in [-yp, yp]: + name = f"output_{i}_{y:+d}" + ref[name] = collect(name, **collect_kw) + else: + for name, val in ref.items(): + assert np.allclose(val, collect(name, **collect_kw)) 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..f5460aff2a 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"] @@ -48,11 +46,9 @@ 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 - ) - ) + args = f" mesh:nx={nx + 4} mesh:dx={dx} MZ={nx} xzinterpolation:type={method}" + nproc = 2 if method == "hermitespline" and boutconfig.has["petsc"] else 1 + args += f" NXPE={nproc}" cmd = "./test_interpolate" + args @@ -71,6 +67,17 @@ for method in methods: E = interp - solution + if False: + import matplotlib.pyplot as plt + + def myplot(f, lbl=None): + plt.plot(f[:, 0, 6], label=lbl) + + myplot(interp, "interp") + myplot(solution, "sol") + plt.legend() + plt.show() + l2 = float(sqrt(mean(E**2))) linf = float(max(abs(E))) diff --git a/tests/integrated/test-interpolate/test_interpolate.cxx b/tests/integrated/test-interpolate/test_interpolate.cxx index a090877552..33963dbb9e 100644 --- a/tests/integrated/test-interpolate/test_interpolate.cxx +++ b/tests/integrated/test-interpolate/test_interpolate.cxx @@ -30,88 +30,90 @@ 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 = 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 = 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 = 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(); + 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 * 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); } - 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, deltax, deltaz); + b_interp = interp->interpolate(b, deltax, deltaz); + c_interp = interp->interpolate(c, deltax, deltaz); - Options dump; + Options dump; - dump["a"] = a; - dump["a_interp"] = a_interp; - dump["a_solution"] = a_solution; + dump["a"] = a; + 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; + 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; + dump["c_interp"] = c_interp; + dump["c_solution"] = c_solution; - bout::writeDefaultOutputFile(dump); + bout::writeDefaultOutputFile(dump); - bout::checkForUnusedOptions(); + bout::checkForUnusedOptions(); + } BoutFinalise(); return 0;