Allow X splitting with FCI

.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

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()

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,57 +52,52 @@ public:
 protected:
   Mesh* localmesh{nullptr};
 
-  std::string region_name;
-  std::shared_ptr<Region<Ind3D>> region{nullptr};
+  int region_id{-1};
 
 public:
   XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr)
       : y_offset(y_offset),
         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<Ind3D>> 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<Ind3D>& 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<Ind3D>>& region) {
-    this->region_name = "";
-    this->region = region;
+    this->region_id = localmesh->getRegionID(region_name);
   }
+  void setRegion(const std::unique_ptr<Region<Ind3D>> region) { setRegion(*region); }
   void setRegion(const Region<Ind3D>& region) {
-    this->region_name = "";
-    this->region = std::make_shared<Region<Ind3D>>(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<Ind3D> getRegion() const {
-    if (!region_name.empty()) {
-      return localmesh->getRegion(region_name);
-    }
-    ASSERT1(region != nullptr);
-    return *region;
+  const Region<Ind3D>& getRegion() const {
+    ASSERT2(region_id != -1);
+    return localmesh->getRegion(region_id);
   }
-  Region<Ind3D> 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<Ind3D>& 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<int> i_corner; // x-index of bottom-left grid point
-  Tensor<int> k_corner; // z-index of bottom-left grid point
+  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
 
   // 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<Field3D> 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,

include/bout/utils.hxx

@@ -361,6 +361,13 @@ 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);
+    // ASSERT2(i.ny == n2);
+    ASSERT2(0 <= i.ind && i.ind < n1 * n2 * n3);
+    return data[i.ind];
+  }
 
   Tensor& operator=(const T& val) {
     for (auto& i : data) {