Do things automatically for 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

CMakeLists.txt

@@ -591,6 +591,9 @@ else()
 endif()
 set(BOUT_USE_METRIC_3D ${BOUT_ENABLE_METRIC_3D})
 
+option(BOUT_ENABLE_FCI_AUTOMAGIC "Enable (slow?) automatic features for FCI" ON)
+set(BOUT_USE_FCI_AUTOMAGIC ${BOUT_ENABLE_FCI_AUTOMAGIC})
+
 include(CheckCXXSourceCompiles)
 check_cxx_source_compiles("int main() { const char* name = __PRETTY_FUNCTION__; }"
   HAS_PRETTY_FUNCTION)

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

cmake/SetupBOUTThirdParty.cmake

@@ -288,7 +288,7 @@ if (BOUT_USE_SUNDIALS)
     set(EXAMPLES_ENABLE_C OFF CACHE BOOL "" FORCE)
     set(EXAMPLES_INSTALL OFF CACHE BOOL "" FORCE)
     set(ENABLE_MPI ${BOUT_USE_MPI} CACHE BOOL "" FORCE)
-    set(ENABLE_OPENMP OFF CACHE BOOL "" FORCE)
+    set(ENABLE_OPENMP ${BOUT_USE_OPENMP} CACHE BOOL "" FORCE)
     if (BUILD_SHARED_LIBS)
       set(BUILD_STATIC_LIBS OFF CACHE BOOL "" FORCE)
     else()

cmake_build_defines.hxx.in

@@ -35,6 +35,7 @@
 #cmakedefine BOUT_METRIC_TYPE @BOUT_METRIC_TYPE@
 #cmakedefine01 BOUT_USE_METRIC_3D
 #cmakedefine01 BOUT_USE_MSGSTACK
+#cmakedefine01 BOUT_USE_FCI_AUTOMAGIC
 
 // CMake build does not support legacy interface
 #define BOUT_HAS_LEGACY_NETCDF 0

examples/fci-wave-logn/boundary/BOUT.inp

@@ -40,5 +40,5 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_yup = parallel_dirichlet(+1.0)
-bndry_par_ydown = parallel_dirichlet(-1.0)
+bndry_par_yup = parallel_dirichlet_o2(+1.0)
+bndry_par_ydown = parallel_dirichlet_o2(-1.0)

examples/fci-wave-logn/div-integrate/BOUT.inp

@@ -40,4 +40,4 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2

examples/fci-wave-logn/expanded/BOUT.inp

@@ -40,4 +40,4 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2

examples/fci-wave/div-integrate/BOUT.inp

@@ -41,4 +41,4 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2

examples/fci-wave/div/BOUT.inp

@@ -41,4 +41,4 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2

examples/fci-wave/logn/BOUT.inp

@@ -41,4 +41,4 @@ bndry_par_ydown = parallel_neumann_o2
 
 [nv]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2

include/bout/boundary_iterator.hxx

@@ -0,0 +1,117 @@
+#pragma once
+
+#include "bout/mesh.hxx"
+#include "bout/sys/parallel_stencils.hxx"
+#include "bout/sys/range.hxx"
+
+class BoundaryRegionIter {
+public:
+  BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh)
+      : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) {
+    ASSERT3(bx * by == 0);
+  }
+  bool operator!=(const BoundaryRegionIter& rhs) { return ind() != rhs.ind(); }
+
+  Ind3D ind() const { return xyz2ind(x, y, z); }
+  BoundaryRegionIter& operator++() {
+    ASSERT3(z < nz());
+    z++;
+    if (z == nz()) {
+      z = 0;
+      _next();
+    }
+    return *this;
+  }
+  virtual void _next() = 0;
+  BoundaryRegionIter& operator*() { return *this; }
+
+  void dirichlet_o2(Field3D& f, BoutReal value) const {
+    ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value);
+  }
+
+  BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); }
+
+  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] * 3 - yprev(f)) * 0.5;
+  }
+
+  BoutReal extrapolate_next_o2(const Field3D& f) const { return 2 * f[ind()] - yprev(f); }
+
+  BoutReal
+  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return 2 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx));
+  }
+
+  BoutReal interpolate_sheath(const Field3D& f) const {
+    return (f[ind()] + ynext(f)) * 0.5;
+  }
+
+  BoutReal& ynext(Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+
+  const int dir;
+
+protected:
+  int z{0};
+  int x;
+  int y;
+  const int bx;
+  const int by;
+
+private:
+  Mesh* localmesh;
+  int nx() const { return localmesh->LocalNx; }
+  int ny() const { return localmesh->LocalNy; }
+  int nz() const { return localmesh->LocalNz; }
+
+  Ind3D xyz2ind(int x, int y, int z) const {
+    return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()};
+  }
+};
+
+class BoundaryRegionIterY : public BoundaryRegionIter {
+public:
+  BoundaryRegionIterY(RangeIterator r, int y, int dir, bool is_end, Mesh* mesh)
+      : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {}
+
+  bool operator!=(const BoundaryRegionIterY& rhs) {
+    ASSERT2(y == rhs.y);
+    if (is_end) {
+      if (rhs.is_end) {
+        return false;
+      }
+      return !rhs.r.isDone();
+    }
+    if (rhs.is_end) {
+      return !r.isDone();
+    }
+    return x != rhs.x;
+  }
+
+  virtual void _next() override {
+    ++r;
+    x = r.ind;
+  }
+
+private:
+  RangeIterator r;
+  bool is_end;
+};
+
+class NewBoundaryRegionY {
+public:
+  NewBoundaryRegionY(Mesh* mesh, bool lower, RangeIterator r)
+      : mesh(mesh), lower(lower), r(std::move(r)) {}
+  BoundaryRegionIterY begin(bool begin = true) {
+    return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1,
+                               !begin, mesh);
+  }
+  BoundaryRegionIterY end() { return begin(false); }
+
+private:
+  Mesh* mesh;
+  bool lower;
+  RangeIterator r;
+};

include/bout/boundary_region.hxx

@@ -4,6 +4,9 @@ class BoundaryRegion;
 #ifndef BOUT_BNDRY_REGION_H
 #define BOUT_BNDRY_REGION_H
 
+#include "bout/mesh.hxx"
+#include "bout/region.hxx"
+#include "bout/sys/parallel_stencils.hxx"
 #include <string>
 #include <utility>
 

include/bout/coordinates.hxx

@@ -133,9 +133,10 @@ public:
     transform = std::move(pt);
   }
 
+  bool hasParallelTransform() const { return transform != nullptr; }
   /// Return the parallel transform
-  ParallelTransform& getParallelTransform() {
-    ASSERT1(transform != nullptr);
+  ParallelTransform& getParallelTransform() const {
+    ASSERT1(hasParallelTransform());
     return *transform;
   }
 

include/bout/field.hxx

@@ -86,6 +86,8 @@ public:
 
   std::string name;
 
+  bool isFci() const;
+
 #if CHECK > 0
   // Routines to test guard/boundary cells set
 
@@ -677,7 +679,27 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") {
       result[d] = f;
     }
   }
-
+  if constexpr (bout::utils::is_Field3D<T>()) {
+#if BOUT_USE_FCI_AUTOMAGIC
+    if (var.isFci()) {
+      for (size_t i = 0; i < result.numberParallelSlices(); ++i) {
+        BOUT_FOR(d, result.yup(i).getRegion(rgn)) {
+          if (result.yup(i)[d] < f) {
+            result.yup(i)[d] = f;
+          }
+        }
+        BOUT_FOR(d, result.ydown(i).getRegion(rgn)) {
+          if (result.ydown(i)[d] < f) {
+            result.ydown(i)[d] = f;
+          }
+        }
+      }
+    } else
+#endif
+    {
+      result.clearParallelSlices();
+    }
+  }
   return result;
 }
 

include/bout/field2d.hxx

@@ -135,8 +135,11 @@ public:
     return *this;
   }
 
-  /// Check if this field has yup and ydown fields
+  /// Dummy functions to increase portability
   bool hasParallelSlices() const { return true; }
+  void calcParallelSlices() const {}
+  void clearParallelSlices() {}
+  int numberParallelSlices() { return 0; }
 
   Field2D& yup(std::vector<Field2D>::size_type UNUSED(index) = 0) { return *this; }
   const Field2D& yup(std::vector<Field2D>::size_type UNUSED(index) = 0) const {
@@ -280,7 +283,7 @@ public:
 
   friend void swap(Field2D& first, Field2D& second) noexcept;
 
-  int size() const override { return nx * ny; };
+  int size() const override { return nx * ny; }
 
 private:
   /// Internal data array. Handles allocation/freeing of memory

include/bout/field3d.hxx

@@ -261,27 +261,38 @@ public:
 #endif
   }
 
+  /// get number of parallel slices
+  size_t numberParallelSlices() const {
+    // Do checks
+    hasParallelSlices();
+    return yup_fields.size();
+  }
+
   /// Check if this field has yup and ydown fields
   /// Return reference to yup field
   Field3D& yup(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(allow_parallel_slices);
     return yup_fields[index];
   }
   /// Return const reference to yup field
   const Field3D& yup(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < yup_fields.size());
+    ASSERT2(allow_parallel_slices);
     return yup_fields[index];
   }
 
   /// Return reference to ydown field
   Field3D& ydown(std::vector<Field3D>::size_type index = 0) {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(allow_parallel_slices);
     return ydown_fields[index];
   }
 
   /// Return const reference to ydown field
   const Field3D& ydown(std::vector<Field3D>::size_type index = 0) const {
     ASSERT2(index < ydown_fields.size());
+    ASSERT2(allow_parallel_slices);
     return ydown_fields[index];
   }
 
@@ -473,6 +484,11 @@ public:
   friend class Vector2D;
 
   Field3D& calcParallelSlices();
+  void allowParallelSlices([[maybe_unused]] bool allow) {
+#if CHECK > 0
+    allow_parallel_slices = allow;
+#endif
+  }
 
   void applyBoundary(bool init = false) override;
   void applyBoundary(BoutReal t);
@@ -487,6 +503,7 @@ public:
   /// Note: does not just copy values in boundary region.
   void setBoundaryTo(const Field3D& f3d);
 
+  using FieldData::applyParallelBoundary;
   void applyParallelBoundary() override;
   void applyParallelBoundary(BoutReal t) override;
   void applyParallelBoundary(const std::string& condition) override;
@@ -514,6 +531,8 @@ private:
 
   /// RegionID over which the field is valid
   std::optional<size_t> regionID;
+
+  bool allow_parallel_slices{true};
 };
 
 // Non-member overloaded operators
@@ -656,4 +675,14 @@ bool operator==(const Field3D& a, const Field3D& b);
 /// Output a string describing a Field3D to a stream
 std::ostream& operator<<(std::ostream& out, const Field3D& value);
 
+inline Field3D copy(const Field3D& f) {
+  Field3D result{f};
+  result.allocate();
+  for (size_t i = 0; i < result.numberParallelSlices(); ++i) {
+    result.yup(i).allocate();
+    result.ydown(i).allocate();
+  }
+  return result;
+}
+
 #endif /* BOUT_FIELD3D_H */

include/bout/field_data.hxx

@@ -43,7 +43,8 @@ class BoundaryOpPar;
 class Coordinates;
 class Mesh;
 
-#include "bout/boundary_region.hxx"
+class BoundaryRegion;
+//#include "bout/boundary_region.hxx"
 class BoundaryRegionPar;
 enum class BndryLoc;