WIP: Tidying FCI implementation

CMakeLists.txt

@@ -593,9 +593,6 @@ 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)

include/bout/field.hxx

@@ -527,22 +527,36 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") {
  * and uses checkData() to, if CHECK >= 3, check
  * result for non-finite numbers
  *
+ * If the input field has parallel slices, then those will also be
+ * operated on. To avoid this, use `withoutParallelSlices()` to
+ * discard slices before calling.
+ *
  */
 #ifdef FIELD_FUNC
 #error This macro has already been defined
 #else
-#define FIELD_FUNC(_name, func)                                        \
-  template <typename T, typename = bout::utils::EnableIfField<T>>      \
-  inline T _name(const T& f, const std::string& rgn = "RGN_ALL") {     \
-    AUTO_TRACE();                                                      \
-    /* Check if the input is allocated */                              \
-    checkData(f);                                                      \
-    /* Define and allocate the output result */                        \
-    T result{emptyFrom(f)};                                            \
-    BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); }     \
-    result.name = std::string(#_name "(") + f.name + std::string(")"); \
-    checkData(result);                                                 \
-    return result;                                                     \
+#define FIELD_FUNC(_name, func)                                                          \
+  template <typename T, typename = bout::utils::EnableIfField<T>>                        \
+  inline T _name(const T& f, const std::string& rgn = "RGN_ALL") {                       \
+    AUTO_TRACE();                                                                        \
+    /* Check if the input is allocated */                                                \
+    checkData(f);                                                                        \
+    /* Define and allocate the output result */                                          \
+    T result{emptyFrom(f)};                                                              \
+    BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); }                       \
+    if (f.hasParallelSlices()) {                                                         \
+      /* Operate on parallel slices */                                                   \
+      result.splitParallelSlicesAndAllocate();                                           \
+      for (size_t i{0}; i != f.numberParallelSlices(); ++i) {                            \
+        BOUT_FOR(d, result.getRegion(rgn)) { result.yup(i)[d] = func(f.yup(i)[d]); }     \
+        checkData(result.yup(i));                                                        \
+        BOUT_FOR(d, result.getRegion(rgn)) { result.ydown(i)[d] = func(f.ydown(i)[d]); } \
+        checkData(result.ydown(i));                                                      \
+      }                                                                                  \
+    }                                                                                    \
+    result.name = std::string(#_name "(") + f.name + std::string(")");                   \
+    checkData(result);                                                                   \
+    return result;                                                                       \
   }
 #endif
 
@@ -670,6 +684,8 @@ T copy(const T& f) {
 /// Apply a floor value \p f to a field \p var. Any value lower than
 /// the floor is set to the floor.
 ///
+/// If the field has parallel slices then they will also be floored.
+///
 /// @param[in] var  Variable to apply floor to
 /// @param[in] f    The floor value
 /// @param[in] rgn  The region to calculate the result over
@@ -683,8 +699,8 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") {
       result[d] = f;
     }
   }
-#if BOUT_USE_FCI_AUTOMAGIC
-  if (var.isFci()) {
+
+  if (var.hasParallelSlices()) {
     for (size_t i = 0; i < result.numberParallelSlices(); ++i) {
       BOUT_FOR(d, result.yup(i).getRegion(rgn)) {
         if (result.yup(i)[d] < f) {
@@ -697,10 +713,6 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") {
         }
       }
     }
-  } else
-#endif
-  {
-    result.clearParallelSlices();
   }
   return result;
 }

include/bout/field2d.hxx

@@ -134,12 +134,14 @@ public:
     return *this;
   }
 
-  /// Dummy functions to increase portability
+  /// Dummy functions to replicate Field3D interface
   bool hasParallelSlices() const { return true; }
   void calcParallelSlices() const {}
   void splitParallelSlices() const {}
+  void splitParallelSlicesAndAllocate() const {}
   void clearParallelSlices() const {}
-  int numberParallelSlices() const { return 0; }
+  Field2D withoutParallelSlices() const { return *this; }
+  size_t numberParallelSlices() const { 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 {

include/bout/field3d.hxx

@@ -237,13 +237,23 @@ public:
    */
   void splitParallelSlices();
 
+  /// Separate fields for yup and ydown, with memory allocated. Note:
+  /// After this the parallel slices will be allocated and unique, but
+  /// may contain uninitialised values.
   void splitParallelSlicesAndAllocate();
 
   /*!
    * Clear the parallel slices, yup and ydown
    */
   void clearParallelSlices();
 
+  /// Returns a shallow copy without parallel slices
+  Field3D withoutParallelSlices() const {
+    Field3D result{getMesh(), getLocation(), getDirections(), getRegionID()};
+    result.data = data;
+    return result;
+  }
+
   /// Check if this field has yup and ydown fields
   bool hasParallelSlices() const {
 #if CHECK > 2
@@ -488,10 +498,11 @@ public:
   friend class Vector2D;
 
   Field3D& calcParallelSlices();
-  void allowParallelSlices([[maybe_unused]] bool allow) {
+  Field3D& allowParallelSlices([[maybe_unused]] bool allow) {
 #if CHECK > 0
     allowCalcParallelSlices = allow;
 #endif
+    return *this;
   }
 
   void applyBoundary(bool init = false) override;

include/bout/fieldperp.hxx

@@ -158,11 +158,14 @@ public:
     return *this;
   }
 
-  /// Dummy functions to increase portability
+  /// Dummy functions to match Field3D interface
   bool hasParallelSlices() const { return true; }
+  void splitParallelSlices() const {}
+  void splitParallelSlicesAndAllocate() const {}
   void calcParallelSlices() const {}
-  void clearParallelSlices() {}
-  int numberParallelSlices() { return 0; }
+  void clearParallelSlices() const {}
+  FieldPerp withoutParallelSlices() const { return *this; }
+  size_t numberParallelSlices() const { return 0; }
 
   FieldPerp& yup(std::vector<FieldPerp>::size_type UNUSED(index) = 0) { return *this; }
   const FieldPerp& yup(std::vector<FieldPerp>::size_type UNUSED(index) = 0) const {

include/bout/fv_ops.hxx

@@ -188,16 +188,13 @@ void communicateFluxes(Field3D& f);
 /// @param[in] fixflux     Fix the flux at the boundary to be the value at the
 ///                        midpoint (for boundary conditions)
 ///
-/// NB: Uses to/from FieldAligned coordinates
+/// Handles FCI by calling ::Div_par(f_in, v_in)
 template <typename CellEdges = MC>
 const Field3D Div_par(const Field3D& f_in, const Field3D& v_in,
                       const Field3D& wave_speed_in, bool fixflux = true) {
-
-#if BOUT_USE_FCI_AUTOMAGIC
   if (f_in.isFci()) {
     return ::Div_par(f_in, v_in);
   }
-#endif
 
   ASSERT1_FIELDS_COMPATIBLE(f_in, v_in);
   ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in);

include/bout/index_derivs_interface.hxx

@@ -4,10 +4,9 @@
  * Definition of main derivative kernels
  *
  **************************************************************************
- * Copyright 2018
- *    D.Dickinson
+ * Copyright 2018 - 2025 BOUT++ contributors
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
  * This file is part of BOUT++.
  *
@@ -204,13 +203,9 @@ T DDY(const T& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "D
     ASSERT1(f.getDirectionY() == YDirectionType::Standard);
     T f_tmp = f;
     if (!f.hasParallelSlices()) {
-#if BOUT_USE_FCI_AUTOMAGIC
-      f_tmp.calcParallelSlices();
-#else
       throw BoutException(
           "parallel slices needed for parallel derivatives. Make sure to communicate and "
           "apply parallel boundary conditions before calling derivative");
-#endif
     }
     return standardDerivative<T, DIRECTION::YOrthogonal, DERIV::Standard>(f_tmp, outloc,
                                                                           method, region);

include/bout/parallel_boundary_region.hxx

@@ -150,7 +150,7 @@ public:
     return bndry_position != rhs.bndry_position;
   }
 
-#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i)
+#define ITER() for (int i = 0; i <= localmesh->ystart - abs_offset(); ++i)
   // dirichlet boundary code
   void dirichlet_o1(Field3D& f, BoutReal value) const {
     ITER() { getAt(f, i) = value; }

src/field/field3d.cxx

@@ -93,15 +93,14 @@ Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field3D(localmesh) {
   TRACE("Field3D: Copy constructor from value");
 
   *this = val;
-#if BOUT_USE_FCI_AUTOMAGIC
+
   if (this->isFci()) {
     splitParallelSlices();
     for (size_t i = 0; i < numberParallelSlices(); ++i) {
       yup(i) = val;
       ydown(i) = val;
     }
   }
-#endif
 }
 
 Field3D::Field3D(Array<BoutReal> data_in, Mesh* localmesh, CELL_LOC datalocation,
@@ -361,18 +360,13 @@ Field3D& Field3D::operator=(const BoutReal val) {
   TRACE("Field3D = BoutReal");
   track(val, "operator=");
 
-#if BOUT_USE_FCI_AUTOMAGIC
-  if (isFci() && hasParallelSlices()) {
+  if (hasParallelSlices()) {
     for (size_t i = 0; i < numberParallelSlices(); ++i) {
       yup(i) = val;
       ydown(i) = val;
     }
   }
-#else
-  // Delete existing parallel slices. We don't copy parallel slices, so any
-  // that currently exist will be incorrect.
-  clearParallelSlices();
-#endif
+
   resetRegion();
 
   allocate();
@@ -386,11 +380,13 @@ Field3D& Field3D::operator=(const BoutReal val) {
 Field3D& Field3D::calcParallelSlices() {
   ASSERT2(allowCalcParallelSlices);
   getCoordinates()->getParallelTransform().calcParallelSlices(*this);
-#if BOUT_USE_FCI_AUTOMAGIC
-  if (this->isFci()) {
-    this->applyParallelBoundary("parallel_neumann_o2");
-  }
-#endif
+
+  // This seems like it could work but give wrong results.
+  // Better to insert NaNs into boundary cells
+  // if (this->isFci()) {
+  //   this->applyParallelBoundary("parallel_neumann_o2");
+  // }
+
   return *this;
 }
 
@@ -548,7 +544,12 @@ void Field3D::applyParallelBoundary() {
   TRACE("Field3D::applyParallelBoundary()");
 
   checkData(*this);
-  ASSERT1(hasParallelSlices());
+  if (isFci()) {
+    ASSERT1(hasParallelSlices());
+  }
+  if (!hasParallelSlices()) {
+    return;
+  }
 
   // Apply boundary to this field
   for (const auto& bndry : getBoundaryOpPars()) {
@@ -561,7 +562,12 @@ void Field3D::applyParallelBoundary(BoutReal t) {
   TRACE("Field3D::applyParallelBoundary(t)");
 
   checkData(*this);
-  ASSERT1(hasParallelSlices());
+  if (isFci()) {
+    ASSERT1(hasParallelSlices());
+  }
+  if (!hasParallelSlices()) {
+    return;
+  }
 
   // Apply boundary to this field
   for (const auto& bndry : getBoundaryOpPars()) {
@@ -574,7 +580,12 @@ void Field3D::applyParallelBoundary(const std::string& condition) {
   TRACE("Field3D::applyParallelBoundary(condition)");
 
   checkData(*this);
-  ASSERT1(hasParallelSlices());
+  if (isFci()) {
+    ASSERT1(hasParallelSlices());
+  }
+  if (!hasParallelSlices()) {
+    return;
+  }
 
   /// Get the boundary factory (singleton)
   BoundaryFactory* bfact = BoundaryFactory::getInstance();
@@ -593,7 +604,12 @@ void Field3D::applyParallelBoundary(const std::string& region,
   TRACE("Field3D::applyParallelBoundary(region, condition)");
 
   checkData(*this);
-  ASSERT1(hasParallelSlices());
+  if (isFci()) {
+    ASSERT1(hasParallelSlices());
+  }
+  if (!hasParallelSlices()) {
+    return;
+  }
 
   /// Get the boundary factory (singleton)
   BoundaryFactory* bfact = BoundaryFactory::getInstance();
@@ -615,7 +631,12 @@ void Field3D::applyParallelBoundary(const std::string& region,
   TRACE("Field3D::applyParallelBoundary(region, condition, f)");
 
   checkData(*this);
-  ASSERT1(hasParallelSlices());
+  if (isFci()) {
+    ASSERT1(hasParallelSlices());
+  }
+  if (!hasParallelSlices()) {
+    return;
+  }
 
   /// Get the boundary factory (singleton)
   BoundaryFactory* bfact = BoundaryFactory::getInstance();