diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index e20aef2459..4c530d8f46 100644 --- a/include/bout/difops.hxx +++ b/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); } diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 5a7df5b088..8373cfeefe 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -975,6 +976,16 @@ public: /// @returns A @ref Parallel struct containing all assembled operators. Parallel getParallel() const; + class Cache { + public: + static std::shared_ptr get(Mesh* mesh); + static std::shared_ptr getParallel(Mesh* mesh); + + private: + static std::map> petscOp; + static std::map> parOp; + }; + private: /// @brief Read a Field3D from the mesh, throwing if the field is absent. /// diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 194d561c0b..ab50f684c7 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -40,6 +40,7 @@ #include #include // Delp2 uses same coefficients as inversion code #include +#include #include #include #include @@ -378,6 +379,26 @@ 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(); @@ -385,7 +406,6 @@ Field3D Div_par_K_Grad_par_mod_impl(const Field3DParallel& Kin, 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(Kin, fin, flow_ylow, bndry_flux); case Harmonic: return Div_par_K_Grad_par_mod_impl(Kin, fin, flow_ylow, bndry_flux); + case PetscOperator: + return Div_par_K_Grad_par_mod_impl(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 diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index b05f013023..66d16ad465 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -1,4 +1,5 @@ #include "bout/build_defines.hxx" +#include #if BOUT_HAS_PETSC @@ -559,4 +560,29 @@ PetscCellOperator makeNeumannOperator(const PetscCellMappingPtr& mapping, return PetscCellOperator(mapping, mapping, std::move(out)); } +std::shared_ptr PetscOperators::Cache::get(Mesh* mesh) { + if (auto search = petscOp.find(mesh); search != petscOp.end()) { + if (std::shared_ptr spt = search->second.lock()) { + return spt; + } + } + auto op = std::make_shared(mesh); + petscOp[mesh] = op; + return op; +} + +std::shared_ptr PetscOperators::Cache::getParallel(Mesh* mesh) { + if (auto search = parOp.find(mesh); search != parOp.end()) { + if (std::shared_ptr spt = search->second.lock()) { + return spt; + } + } + auto op = get(mesh); + auto par = std::make_shared(op->getParallel()); + parOp[mesh] = par; + return par; +} + +std::map> PetscOperators::Cache::parOp; +std::map> PetscOperators::Cache::petscOp; #endif diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 75204de221..94cd50aa46 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/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 diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index a1bc20e0a5..9f06b79114 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/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 diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 4c776c7364..cd84cd6805 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/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,