Field3DParallel

.ci_fedora.sh

@@ -47,6 +47,8 @@ test . != ".$1" && mpi="$1" || mpi=openmpi
     export OMP_NUM_THREADS=1
     cd
     cd BOUT-dev
+    python3 -m ensurepip
+    python3 -m pip install -r requirements.txt
     echo "starting configure"
     time cmake -S . -B build -DBOUT_USE_PETSC=ON \
 	 -DBOUT_UPDATE_GIT_SUBMODULE=OFF \

.github/workflows/tests.yml

@@ -187,7 +187,7 @@ jobs:
         run: BUILD_PETSC=${{ matrix.config.build_petsc }} ./.build_petsc_for_ci.sh ${{ matrix.config.build_petsc_branch }}
 
       - name: Build ADIOS2
-        run: BUILD_ADIOS2=${{ matrix.config.build_adios2 }} ./.build_adios2_for_ci.sh 
+        run: BUILD_ADIOS2=${{ matrix.config.build_adios2 }} ./.build_adios2_for_ci.sh
 
       - name: Build BOUT++
         run: UNIT_ONLY=${{ matrix.config.unit_only }} ./.ci_with_cmake.sh ${{ matrix.config.cmake_options }}

include/bout/coordinates.hxx

@@ -165,7 +165,7 @@ public:
               const std::string& method = "DEFAULT",
               const std::string& region = "RGN_NOBNDRY");
 
-  Field3D DDY(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+  Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
               const std::string& method = "DEFAULT",
               const std::string& region = "RGN_NOBNDRY") const;
 
@@ -177,30 +177,30 @@ public:
   FieldMetric Grad_par(const Field2D& var, CELL_LOC outloc = CELL_DEFAULT,
                        const std::string& method = "DEFAULT");
 
-  Field3D Grad_par(const Field3D& var, CELL_LOC outloc = CELL_DEFAULT,
+  Field3D Grad_par(const Field3DParallel& var, CELL_LOC outloc = CELL_DEFAULT,
                    const std::string& method = "DEFAULT");
 
   /// Advection along magnetic field V*b.Grad(f)
   FieldMetric Vpar_Grad_par(const Field2D& v, const Field2D& f,
                             CELL_LOC outloc = CELL_DEFAULT,
                             const std::string& method = "DEFAULT");
 
-  Field3D Vpar_Grad_par(const Field3D& v, const Field3D& f,
+  Field3D Vpar_Grad_par(const Field3D& v, const Field3DParallel& f,
                         CELL_LOC outloc = CELL_DEFAULT,
                         const std::string& method = "DEFAULT");
 
   /// Divergence along magnetic field  Div(b*f) = B.Grad(f/B)
   FieldMetric Div_par(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                       const std::string& method = "DEFAULT");
 
-  Field3D Div_par(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+  Field3D Div_par(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
                   const std::string& method = "DEFAULT");
 
   // Second derivative along magnetic field
   FieldMetric Grad2_par2(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                          const std::string& method = "DEFAULT");
 
-  Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+  Field3D Grad2_par2(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
                      const std::string& method = "DEFAULT");
   // Perpendicular Laplacian operator, using only X-Z derivatives
   // NOTE: This might be better bundled with the Laplacian inversion code
@@ -212,13 +212,13 @@ public:
   // Full parallel Laplacian operator on scalar field
   // Laplace_par(f) = Div( b (b dot Grad(f)) )
   FieldMetric Laplace_par(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT);
-  Field3D Laplace_par(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT);
+  Field3D Laplace_par(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT);
 
   // Full Laplacian operator on scalar field
   FieldMetric Laplace(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                       const std::string& dfdy_boundary_conditions = "free_o3",
                       const std::string& dfdy_dy_region = "");
-  Field3D Laplace(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+  Field3D Laplace(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
                   const std::string& dfdy_boundary_conditions = "free_o3",
                   const std::string& dfdy_dy_region = "");
 

include/bout/deriv_store.hxx

@@ -33,8 +33,10 @@
 #include <functional>
 #include <map>
 #include <set>
+#include <type_traits>
 #include <unordered_map>
 
+#include "bout/field3d.hxx"
 #include <bout/scorepwrapper.hxx>
 
 #include <bout/bout_types.hxx>
@@ -522,4 +524,13 @@ private:
   }
 };
 
+template <typename T>
+auto& getStore() {
+  if constexpr (std::is_same<T, Field3DParallel>::value) {
+    return DerivativeStore<Field3D>::getInstance();
+  } else {
+    return DerivativeStore<T>::getInstance();
+  }
+}
+
 #endif

include/bout/derivs.hxx

@@ -82,7 +82,7 @@ Coordinates::FieldMetric DDX(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
 ///                    If not given, defaults to DIFF_DEFAULT
 /// @param[in] region  What region is expected to be calculated
 ///                    If not given, defaults to RGN_NOBNDRY
-Field3D DDY(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
             const std::string& method = "DEFAULT",
             const std::string& region = "RGN_NOBNDRY");
 
@@ -410,7 +410,7 @@ Coordinates::FieldMetric VDDX(const Field2D& v, const Field2D& f,
 ///                    If not given, defaults to DIFF_DEFAULT
 /// @param[in] region  What region is expected to be calculated
 ///                    If not given, defaults to RGN_NOBNDRY
-Field3D VDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+Field3D VDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
              const std::string& method = "DEFAULT",
              const std::string& region = "RGN_NOBNDRY");
 
@@ -533,7 +533,7 @@ Coordinates::FieldMetric FDDX(const Field2D& v, const Field2D& f,
 ///                    If not given, defaults to DIFF_DEFAULT
 /// @param[in] region  What region is expected to be calculated
 ///                    If not given, defaults to RGN_NOBNDRY
-Field3D FDDY(const Field3D& v, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
+Field3D FDDY(const Field3D& v, const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
              const std::string& method = "DEFAULT",
              const std::string& region = "RGN_NOBNDRY");
 

include/bout/field.hxx

@@ -137,6 +137,7 @@ public:
   virtual void calcParallelSlices() {}
   virtual void splitParallelSlices() {}
   virtual void clearParallelSlices() {}
+  virtual size_t numberParallelSlices() const { return 0; }
 
 private:
   /// Labels for the type of coordinate system this field is defined over
@@ -189,8 +190,8 @@ inline bool areFieldsCompatible(const Field& field1, const Field& field2) {
 template <typename T>
 inline T emptyFrom(const T& f) {
   static_assert(bout::utils::is_Field_v<T>, "emptyFrom only works on Fields");
-  return T(f.getMesh(), f.getLocation(), {f.getDirectionY(), f.getDirectionZ()},
-           f.getRegionID())
+  return T(f.getMesh(), f.getLocation(),
+           DirectionTypes{f.getDirectionY(), f.getDirectionZ()}, f.getRegionID())
       .allocate();
 }
 
@@ -524,20 +525,27 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") {
  * result for non-finite numbers
  *
  */
+class Field3DParallel;
 #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") {  \
-                                                                   \
-    /* 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]); } \
-    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") {     \
+    /* 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 constexpr (std::is_base_of_v<Field3DParallel, T>) {             \
+      for (int i = 0; i < f.numberParallelSlices(); ++i) {             \
+        result.yup(i) = func(f.yup(i));                                \
+        result.ydown(i) = func(f.ydown(i));                            \
+      }                                                                \
+    }                                                                  \
+    result.name = std::string(#_name "(") + f.name + std::string(")"); \
+    checkData(result);                                                 \
+    return result;                                                     \
   }
 #endif
 
@@ -661,6 +669,8 @@ T copy(const T& f) {
   return result;
 }
 
+class Field3DParallel;
+
 /// Apply a floor value \p f to a field \p var. Any value lower than
 /// the floor is set to the floor.
 ///
@@ -677,7 +687,36 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") {
       result[d] = f;
     }
   }
-
+  if constexpr (std::is_same_v<T, Field3DParallel>) {
+    if (var.hasParallelSlices()) {
+      for (size_t i = 0; i < result.numberParallelSlices(); ++i) {
+        if (result.yup(i).isAllocated()) {
+          BOUT_FOR(d, result.yup(i).getRegion(rgn)) {
+            if (result.yup(i)[d] < f) {
+              result.yup(i)[d] = f;
+            }
+          }
+        } else {
+          if (result.isFci()) {
+            throw BoutException("Expected parallel slice to be allocated");
+          }
+        }
+        if (result.ydown(i).isAllocated()) {
+          BOUT_FOR(d, result.ydown(i).getRegion(rgn)) {
+            if (result.ydown(i)[d] < f) {
+              result.ydown(i)[d] = f;
+            }
+          }
+        } else {
+          if (result.isFci()) {
+            throw BoutException("Expected parallel slice to be allocated");
+          }
+        }
+      }
+    }
+  } else {
+    result.clearParallelSlices();
+  }
   return result;
 }
 

include/bout/field2d.hxx

@@ -4,9 +4,9 @@
  * \brief Definition of 2D scalar field class
  *
  **************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010 - 2026 BOUT++ contributors
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
  * This file is part of BOUT++.
  *
@@ -33,13 +33,17 @@ class Field2D;
 #define BOUT_FIELD2D_H
 
 #include "bout/array.hxx"
-#include "bout/build_config.hxx"
+#include "bout/assert.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/build_defines.hxx"
 #include "bout/field.hxx"
 #include "bout/field_data.hxx"
 #include "bout/fieldperp.hxx"
 #include "bout/region.hxx"
 
 #include <cstddef>
+#include <ostream>
+#include <string>
 
 #if BOUT_HAS_RAJA
 #include "RAJA/RAJA.hpp" // using RAJA lib
@@ -70,7 +74,7 @@ public:
   Field2D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE,
           DirectionTypes directions_in = {YDirectionType::Standard,
                                           ZDirectionType::Average},
-          std::optional<size_t> region = {});
+          std::optional<size_t> regionID = {});
 
   /*!
    * Copy constructor. After this both fields
@@ -206,7 +210,7 @@ public:
     }
 #endif
 
-    return data[jx * ny + jy];
+    return data[(jx * ny) + jy];
   }
   inline const BoutReal& operator()(int jx, int jy) const {
 #if CHECK > 2 && !BOUT_HAS_CUDA
@@ -220,7 +224,7 @@ public:
     }
 #endif
 
-    return data[jx * ny + jy];
+    return data[(jx * ny) + jy];
   }
 
   /*!
@@ -274,6 +278,11 @@ public:
 
   int size() const override { return nx * ny; }
 
+  /// Return a field that can be used to perform parallel (Y) derivatives.
+  /// This is used to write generic FA/FCI code.
+  Field2D& asField3DParallel() { return *this; }
+  const Field2D& asField3DParallel() const { return *this; }
+
 private:
   /// Internal data array. Handles allocation/freeing of memory
   Array<BoutReal> data;
@@ -286,6 +295,10 @@ private:
 };
 
 // Non-member overloaded operators
+FieldPerp operator+(const Field2D& lhs, const FieldPerp& rhs);
+FieldPerp operator-(const Field2D& lhs, const FieldPerp& rhs);
+FieldPerp operator*(const Field2D& lhs, const FieldPerp& rhs);
+FieldPerp operator/(const Field2D& lhs, const FieldPerp& rhs);
 
 Field2D operator+(const Field2D& lhs, const Field2D& rhs);
 Field2D operator-(const Field2D& lhs, const Field2D& rhs);