Port some Hermes-3 finite volume operators (rebased)

.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 git+https://github.com/boutproject/zoidberg@better-metric
     echo "starting configure"
     time cmake -S . -B build -DBOUT_USE_PETSC=ON \
 	 -DBOUT_UPDATE_GIT_SUBMODULE=OFF \

.github/workflows/tests.yml

@@ -166,6 +166,7 @@ jobs:
         run: |
           python -m pip install --upgrade pip setuptools
           python -m pip install -r requirements.txt
+          python -m pip install git+https://github.com/boutproject/zoidberg@better-metric
 
       - name: Cache Zenodo test data
         uses: actions/cache@v5
@@ -187,7 +188,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

@@ -3,16 +3,16 @@
  *
  * ChangeLog
  * =========
- * 
+ *
  * 2014-11-10 Ben Dudson <bd512@york.ac.uk>
  *    * Created by separating metric from Mesh
  *
- * 
+ *
  **************************************************************************
  * Copyright 2014-2025 BOUT++ contributors
  *
  * Contact: Ben Dudson, dudson2@llnl.gov
- * 
+ *
  * This file is part of BOUT++.
  *
  * BOUT++ is free software: you can redistribute it and/or modify
@@ -160,7 +160,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;
 
@@ -172,30 +172,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
@@ -207,13 +207,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 = "");
 
@@ -256,10 +256,10 @@ private:
 class TokamakCoordinates : public Coordinates {
 public:
   TokamakCoordinates(Mesh *mesh) : Coordinates(mesh) {
-    
+
   }
 private:
-  
+
 };
 */
 

include/bout/deriv_store.hxx

@@ -522,4 +522,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/difops.hxx

@@ -1,11 +1,11 @@
 /*!******************************************************************************
  * \file difops.hxx
- * 
+ *
  * Differential operators
  *
  * Changelog:
  *
- * 2009-01 Ben Dudson <bd512@york.ac.uk> 
+ * 2009-01 Ben Dudson <bd512@york.ac.uk>
  *    * Added two optional parameters which can be put in any order
  *      These determine the method to use (DIFF_METHOD)
  *      and CELL_LOC location of the result.
@@ -15,7 +15,7 @@
  * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
  *
  * Contact: Ben Dudson, bd512@york.ac.uk
- * 
+ *
  * This file is part of BOUT++.
  *
  * BOUT++ is free software: you can redistribute it and/or modify
@@ -30,7 +30,7 @@
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
- * 
+ *
  *******************************************************************************/
 
 #ifndef BOUT_DIFOPS_H
@@ -81,7 +81,7 @@ Field3D Grad_parP(const Field3D& apar, const Field3D& f);
  * \f[
  *    v\mathbf{b}_0 \cdot \nabla f
  * \f]
- * 
+ *
  *
  * @param[in] v  The velocity in y direction
  * @param[in] f  The scalar field to be differentiated
@@ -175,7 +175,7 @@ inline Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method)
 
 /*!
  * Parallel divergence of diffusive flux, K*Grad_par
- * 
+ *
  * \f[
  *    \nabla \cdot ( \mathbf{b}_0 kY (\mathbf{b}_0 \cdot \nabla) f )
  * \f]
@@ -195,12 +195,16 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field2D& f,
 Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f,
                            CELL_LOC outloc = CELL_DEFAULT);
 
+/// Version with energy flow diagnostic
+Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow,
+                               bool bndry_flux = true);
+
 /*!
  * Perpendicular Laplacian operator
  *
  * This version only includes terms in X and Z, dropping
- * derivatives in Y. This is the inverse operation to 
- * the Laplacian inversion class. 
+ * derivatives in Y. This is the inverse operation to
+ * the Laplacian inversion class.
  *
  * For the full perpendicular Laplacian, use Laplace_perp
  */
@@ -212,7 +216,7 @@ FieldPerp Delp2(const FieldPerp& f, CELL_LOC outloc = CELL_DEFAULT, bool useFFT
 /*!
  * Perpendicular Laplacian, keeping y derivatives
  *
- * 
+ *
  */
 Coordinates::FieldMetric
 Laplace_perp(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
@@ -246,18 +250,18 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f);
 
 /*!
  * Terms of form b0 x Grad(phi) dot Grad(A)
- * 
+ *
  */
 Coordinates::FieldMetric b0xGrad_dot_Grad(const Field2D& phi, const Field2D& A,
                                           CELL_LOC outloc = CELL_DEFAULT);
 
 /*!
- * Terms of form 
+ * Terms of form
  *
  * \f[
  *   \mathbf{b}_0 \times \nabla \phi \cdot \nabla A
  * \f]
- * 
+ *
  * @param[in] phi The scalar potential
  * @param[in] A   The field being advected
  * @param[in] outloc  The cell location where the result is defined. By default the same as A.
@@ -291,13 +295,13 @@ constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu;
  * \f[
  *   [f, g] = (1/B) \mathbf{b}_0 \times \nabla f  \cdot \nabla g
  * \f]
- * 
+ *
  * @param[in] f  The potential
  * @param[in] g  The field being advected
  * @param[in] method   The method to use
  * @param[in] outloc   The cell location where the result is defined. Default is the same as g
  * @param[in] solver   Pointer to the time integration solver
- * 
+ *
  */
 Coordinates::FieldMetric bracket(const Field2D& f, const Field2D& g,
                                  BRACKET_METHOD method = BRACKET_STD,

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;
 }