Add option to Div_par_K_Grad_par_mod to use petsc support operator method

include/bout/difops.hxx

@@ -217,7 +217,10 @@ Field3D Div_par_K_Grad_par(const Field3DParallel& kY, const Field3DParallel& f,
 /// :math:`K J / (g_{22} dy)`. This better matches a series-resistance
 /// interpretation of two adjacent half-cells and can give noticeably
 /// different results when coefficients or cell sizes vary strongly.
-BOUT_ENUM_CLASS_NS(bout, ConductionMethod, Original, ProductJK, Harmonic);
+///
+/// PetscOperator:
+///  Uses the PETSc support operator method. For FCI only.
+BOUT_ENUM_CLASS_NS(bout, ConductionMethod, Original, ProductJK, Harmonic, PetscOperator);
 
 namespace bout {
 inline std::string format_as(const ConductionMethod& method) { return toString(method); }

include/bout/petsc_operators.hxx

@@ -23,6 +23,7 @@
 #include <petscvec.h>
 
 #include <functional>
+#include <map>
 #include <memory>
 #include <string>
 #include <type_traits>
@@ -975,6 +976,16 @@ public:
   /// @returns A @ref Parallel struct containing all assembled operators.
   Parallel getParallel() const;
 
+  class Cache {
+  public:
+    static std::shared_ptr<PetscOperators> get(Mesh* mesh);
+    static std::shared_ptr<PetscOperators::Parallel> getParallel(Mesh* mesh);
+
+  private:
+    static std::map<Mesh*, std::weak_ptr<PetscOperators>> petscOp;
+    static std::map<Mesh*, std::weak_ptr<PetscOperators::Parallel>> parOp;
+  };
+
 private:
   /// @brief Read a Field3D from the mesh, throwing if the field is absent.
   ///

src/mesh/difops.cxx

@@ -40,6 +40,7 @@
 #include <bout/interpolation.hxx>
 #include <bout/invert_laplace.hxx> // Delp2 uses same coefficients as inversion code
 #include <bout/mesh.hxx>
+#include <bout/petsc_operators.hxx>
 #include <bout/region.hxx>
 #include <bout/solver.hxx>
 #include <bout/unused.hxx>
@@ -378,14 +379,33 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin,
     //       -> Only B, dy and g_22 can be used in yup/ydown
     //          Others (e.g J) may not be averaged between y planes.
 
+    flow_ylow = zeroFrom(fin);
+
+    if constexpr (conduction_method == bout::ConductionMethod::PetscOperator) {
+#if BOUT_HAS_PETSC
+      auto par = PetscOperators::Cache::getParallel(Kin.getMesh());
+      return par->Div_par_K_Grad_par(Kin, fin);
+#else
+      throw BoutException("Cannot use `bout::ConductionMethod::PetscOperator`: BOUT++ "
+                          "was compiled without PETSc.");
+#endif
+    }
+
+    if constexpr (conduction_method != bout::ConductionMethod::Original) {
+      throw BoutException(
+          "Only bout::ConductionMethod::Original and "
+          "\nbout::ConductionMethod::PetscOperator (with PETSc) are supported for FCI"
+          "\nbut found {}.",
+          conduction_method);
+    }
+
     const auto& K_up = Kin.yup();
     const auto& K_down = Kin.ydown();
 
     const auto& f_up = fin.yup();
     const auto& f_down = fin.ydown();
 
     Field3D result{zeroFrom(fin)};
-    flow_ylow = zeroFrom(fin);
 
     const auto yboundary = coord->getYBoundary();
 
@@ -473,6 +493,8 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin,
                 : 0.0;
 
         flux = C_edge * (f[iyp] - f[i]);
+      } else {
+        throw BoutException("Unsupported option {}", conduction_method);
       }
       result[i] += flux / (coord->dy()[i] * coord->J()[i]);
     }
@@ -545,9 +567,12 @@ Field3D Div_par_K_Grad_par_mod(const Field3DParallel& Kin, const Field3DParallel
     return Div_par_K_Grad_par_mod_impl<ProductJK>(Kin, fin, flow_ylow, bndry_flux);
   case Harmonic:
     return Div_par_K_Grad_par_mod_impl<Harmonic>(Kin, fin, flow_ylow, bndry_flux);
+  case PetscOperator:
+    return Div_par_K_Grad_par_mod_impl<PetscOperator>(Kin, fin, flow_ylow, bndry_flux);
   }
-  throw BoutException(
-      "Unknown method `{}` - choose from `Original`, `ProductJK` or `Harmonic`.", method);
+  throw BoutException("Unknown method `{}` - choose from `Original`, `ProductJK`, "
+                      "`Harmonic` or `PetscOperator`.",
+                      method);
 }
 /*******************************************************************************
 * Delp2

src/mesh/petsc_operators.cxx

@@ -1,4 +1,5 @@
 #include "bout/build_defines.hxx"
+#include <map>
 
 #if BOUT_HAS_PETSC
 
@@ -559,4 +560,29 @@ PetscCellOperator makeNeumannOperator(const PetscCellMappingPtr& mapping,
   return PetscCellOperator(mapping, mapping, std::move(out));
 }
 
+std::shared_ptr<PetscOperators> PetscOperators::Cache::get(Mesh* mesh) {
+  if (auto search = petscOp.find(mesh); search != petscOp.end()) {
+    if (std::shared_ptr<PetscOperators> spt = search->second.lock()) {
+      return spt;
+    }
+  }
+  auto op = std::make_shared<PetscOperators>(mesh);
+  petscOp[mesh] = op;
+  return op;
+}
+
+std::shared_ptr<PetscOperators::Parallel> PetscOperators::Cache::getParallel(Mesh* mesh) {
+  if (auto search = parOp.find(mesh); search != parOp.end()) {
+    if (std::shared_ptr<PetscOperators::Parallel> spt = search->second.lock()) {
+      return spt;
+    }
+  }
+  auto op = get(mesh);
+  auto par = std::make_shared<PetscOperators::Parallel>(op->getParallel());
+  parOp[mesh] = par;
+  return par;
+}
+
+std::map<Mesh*, std::weak_ptr<PetscOperators::Parallel>> PetscOperators::Cache::parOp;
+std::map<Mesh*, std::weak_ptr<PetscOperators>> PetscOperators::Cache::petscOp;
 #endif

tests/MMS/spatial/fci/data/BOUT.inp

@@ -9,8 +9,7 @@ FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((
 FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0)
 div_par_K_grad_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0)
 div_par_K_grad_par_mod_Original_solution = div_par_K_grad_par_mod_solution
-div_par_K_grad_par_mod_ProductJK_solution = div_par_K_grad_par_mod_solution
-div_par_K_grad_par_mod_Harmonic_solution = div_par_K_grad_par_mod_solution
+div_par_K_grad_par_mod_PetscOperator_solution = div_par_K_grad_par_mod_solution
 
 [mesh]
 symmetricglobalx = true

tests/MMS/spatial/fci/fci_mms.cxx

@@ -66,11 +66,8 @@ int main(int argc, char** argv) {
       "div_par_K_grad_par_mod_Original", dump, input,
       Div_par_K_Grad_par_mod(K, input, flow_ylow, true, ConductionMethod::Original));
   fci_op_test(
-      "div_par_K_grad_par_mod_ProductJK", dump, input,
-      Div_par_K_Grad_par_mod(K, input, flow_ylow, true, ConductionMethod::ProductJK));
-  fci_op_test(
-      "div_par_K_grad_par_mod_Harmonic", dump, input,
-      Div_par_K_Grad_par_mod(K, input, flow_ylow, true, ConductionMethod::Harmonic));
+      "div_par_K_grad_par_mod_PetscOperator", dump, input,
+      Div_par_K_Grad_par_mod(K, input, flow_ylow, true, ConductionMethod::PetscOperator));
   fci_op_test("laplace_par", dump, input, Laplace_par(input));
 
   // Finite volume methods

tests/MMS/spatial/fci/runtest

@@ -25,8 +25,7 @@ OPERATORS = (
     "div_par",
     "div_par_K_grad_par",
     "div_par_K_grad_par_mod_Original",
-    "div_par_K_grad_par_mod_ProductJK",
-    "div_par_K_grad_par_mod_Harmonic",
+    "div_par_K_grad_par_mod_PetscOperator",
     "laplace_par",
     "FV_div_par_mod",
     "FV_div_par_fvv",
@@ -102,6 +101,7 @@ def make_grid(nz: int):
     # Create the grid
     grid = zb.grid.Grid(poloidal_grid, ycoords, ylength, yperiodic=True)
     maps = zb.make_maps(grid, field, nslice=1, quiet=True, MXG=1)
+    zb.weights.modify_maps(maps, partial_boundaries=False)
     zb.write_maps(
         grid,
         field,