From 289ddf2a9c89ec5b39e8f22f329f9c4bcade0703 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 12 Mar 2026 08:50:06 -0700 Subject: [PATCH 01/26] PetscOperators: FCI operators using PETSc Read operators from the mesh in CSR format. PETSc handles communication pattern for MPI decomposition. --- CMakeLists.txt | 2 + include/bout/petsc_interface.hxx | 4 +- include/bout/petsc_operators.hxx | 30 ++++++ src/mesh/petsc_operators.cxx | 124 +++++++++++++++++++++++ tests/unit/CMakeLists.txt | 1 + tests/unit/mesh/test_petsc_operators.cxx | 44 ++++++++ 6 files changed, 203 insertions(+), 2 deletions(-) create mode 100644 include/bout/petsc_operators.hxx create mode 100644 src/mesh/petsc_operators.cxx create mode 100644 tests/unit/mesh/test_petsc_operators.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index 5f054c1c1f..0b14e5c3d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -177,6 +177,7 @@ set(BOUT_SOURCES ./include/bout/parallel_boundary_region.hxx ./include/bout/paralleltransform.hxx ./include/bout/petsc_interface.hxx + ./include/bout/petsc_operators.hxx ./include/bout/petsclib.hxx ./include/bout/physicsmodel.hxx ./include/bout/rajalib.hxx @@ -308,6 +309,7 @@ set(BOUT_SOURCES ./src/mesh/parallel/shiftedmetricinterp.hxx ./src/mesh/parallel_boundary_op.cxx ./src/mesh/parallel_boundary_region.cxx + ./src/mesh/petsc_operators.cxx ./src/mesh/surfaceiter.cxx ./src/physics/gyro_average.cxx ./src/physics/physicsmodel.cxx diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index 5239039f0c..ef17a63e4a 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -248,7 +248,7 @@ private: }; /*! - * A class which wraps PETSc vector objects, allowing them to be + * A class which wraps PETSc matrix objects, allowing them to be * indexed using the BOUT++ scheme. It provides the option of setting * a y-offset that interpolates onto field lines. */ @@ -289,7 +289,7 @@ public: } // Construct a matrix capable of operating on the specified field, - // preallocating memory if requeted and possible. + // preallocating memory if requested and possible. PetscMatrix(IndexerPtr indConverter, bool preallocate = true) : matrix(new Mat()), indexConverter(indConverter), pt(&indConverter->getMesh()->getCoordinates()->getParallelTransform()) { diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx new file mode 100644 index 0000000000..6801418996 --- /dev/null +++ b/include/bout/petsc_operators.hxx @@ -0,0 +1,30 @@ +/// Represent operators using PETSc matrices +/// + +#include "bout/mesh.hxx" +#include + +#include + +class PetscOperators { +public: + PetscOperators(Mesh* mesh); + + // Static functions for unit testing + static Region create_region(Field3D cell_number); + static Region create_region_xin(Field3D cell_number); + static Region create_region_xout(Field3D cell_number); + +private: + PetscLib lib{}; // Initialize and finalize PETSc + + Region evolving_region; ///< Evolving cells + Region xin_region; ///< X boundary cells in inner boundary + Region xout_region; ///< X boundary cells in outer boundary + Region yup_region; ///< Y boundary cells in forward direction + Region ydown_region; ///< Y boundary cells in backward direction + + Mat mat_mesh_to_petsc; + + Mat mat_yforward; +}; diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx new file mode 100644 index 0000000000..6d3f98791e --- /dev/null +++ b/src/mesh/petsc_operators.cxx @@ -0,0 +1,124 @@ +#include "bout/petsc_operators.hxx" +#include "bout/boutexception.hxx" +#include "bout/output_bout_types.hxx" + +Region PetscOperators::create_region(Field3D cell_number) { + Region::RegionIndices indices; + BOUT_FOR_SERIAL(i, cell_number.getRegion("RGN_NOBNDRY")) { + if (cell_number[i] > -1) { + indices.push_back(i); + } + } + return Region(indices); +} + +Region PetscOperators::create_region_xin(Field3D cell_number) { + Region::RegionIndices xin_indices; + const Mesh* mesh = cell_number.getMesh(); + if (mesh->firstX()) { + for (int i = 0; i < mesh->xstart; ++i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = mesh->zstart; k <= mesh->zend; ++k) { + if (cell_number(i, j, k) > 0) { + xin_indices.push_back(Ind3D{i, j, k}); + } + } + } + } + } + return Region(xin_indices); +} + +Region PetscOperators::create_region_xout(Field3D cell_number) { + Region::RegionIndices xout_indices; + const Mesh* mesh = cell_number.getMesh(); + if (mesh->lastX()) { + for (int i = mesh->xend + 1; i < mesh->LocalNx; ++i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = mesh->zstart; k <= mesh->zend; ++k) { + if (cell_number(i, j, k) > 0) { + xout_indices.push_back(Ind3D{i, j, k}); + } + } + } + } + } + return Region(xout_indices); +} + +PetscOperators::PetscOperators(Mesh* mesh) { + + // Get the cell center global numbering + // Note: This is an integer but Field3D only stores BoutReal. + Field3D cell_number; + if (mesh->get(cell_number, "cell_number") != 0) { + throw BoutException("PetscOperators requires field 'cell_number'"); + } + + // Forward Y boundary cell numbers. + // >0 only if a forward boundary cell + Field3D forward_cell_number; + if (mesh->get(forward_cell_number, "forward_cell_number") != 0) { + throw BoutException("PetscOperators requires field 'forward_cell_number'"); + } + + // Backward Y boundary cell numbers. + // >0 only if a backward boundary cell + Field3D backward_cell_number; + if (mesh->get(backward_cell_number, "backward_cell_number") != 0) { + throw BoutException("PetscOperators requires field 'backward_cell_number'"); + } + + // Calculate Regions for evolving variables and boundaries + evolving_region = create_region(cell_number); + xin_region = create_region_xin(cell_number); + xout_region = create_region_xout(cell_number); + yup_region = create_region(forward_cell_number); + ydown_region = create_region(backward_cell_number); + + const std::size_t nlocal = evolving_region.size() + xin_region.size() + + xout_region.size() + yup_region.size() + + ydown_region.size(); + output.write("Local {} : evolving {} xin {} xout {} yup {} ydown {}", nlocal, + evolving_region.size(), xin_region.size(), xout_region.size(), + yup_region.size(), ydown_region.size()); + + int mesh_total_cells; + if (mesh->get(mesh_total_cells, "total_cells") == 0) { + // Check total number of cells + output.write("Total cells in mesh: {}\n", mesh_total_cells); + } + + // Create a PETSc matrix + // Note: Numbering is different from that used in PetscVector / PetscMatrix + // because yup/down boundaries are included. + + MPI_Comm comm = BoutComm::get(); + + // Renumbering matrix. + // Maps PETSc row (or column) indices to the global indices used in + // the mesh. This is needed because the PETSc indices depend on the + // number of processors. + MatCreate(comm, &mat_mesh_to_petsc); + MatSetSizes(mat_mesh_to_petsc, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); + MatSetType(mat_mesh_to_petsc, MATMPIAIJ); + // Each row will have one non-zero entry, which could be in + // either the "diagonal" or "off-diagonal" block. + MatMPIAIJSetPreallocation(mat_mesh_to_petsc, 1, nullptr, 1, nullptr); + + // Get the range of rows owned by this processor + PetscInt row_start, row_end; + MatGetOwnershipRange(mat_mesh_to_petsc, &row_start, &row_end); + output.write("Local row range: {} -> {}\n", row_start, row_end); + + PetscInt row = row_start; + const PetscScalar ONE = 1.0; + BOUT_FOR_SERIAL(i, evolving_region) { + PetscInt col = ROUND(cell_number[i]); + MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); + ++row; + } + + //MatCreate(comm, &mat_yforward); + //MatSetSizes(mat_yforward, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); +} diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index b347de7354..d3d690417b 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -97,6 +97,7 @@ set(serial_tests_source ./mesh/test_invert3x3.cxx ./mesh/test_mesh.cxx ./mesh/test_paralleltransform.cxx + ./mesh/test_petsc_operators.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_solver.cxx diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx new file mode 100644 index 0000000000..949029f19c --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -0,0 +1,44 @@ +#include "gtest/gtest.h" + +#include "bout/output_bout_types.hxx" +#include "bout/petsc_operators.hxx" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +namespace { +Ind3D indexAt(const Field3D& f, int x, int y, int z) { + int ny = f.getNy(); + int nz = f.getNz(); + return Ind3D{(((x * ny) + y) * nz) + z, ny, nz}; +} +} // namespace + +using bout::globals::mesh; + +// Reuse the "standard" fixture for FakeMesh +using PetscOperatorsTest = FakeMeshFixture; + +TEST_F(PetscOperatorsTest, create_region_empty) { + Field3D cell_number{-1.0}; // No cells >= 0 + + auto rgn = PetscOperators::create_region(cell_number); + + ASSERT_EQ(0, rgn.size()); +} + +TEST_F(PetscOperatorsTest, create_region) { + Field3D cell_number{-1.0}; // No cells >= 0 + + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 0; // Corner + cell_number(1, 1, 0) = 1; // In domain + cell_number(0, 1, 1) = 2; // Xin boundary + + auto rgn = PetscOperators::create_region(cell_number); + ASSERT_EQ(1, rgn.size()); + + Ind3D expected_ind = indexAt(cell_number, 1, 1, 0); + output.write("{} : {}\n", expected_ind, *rgn.begin()); + ASSERT_TRUE(expected_ind == *rgn.begin()); +} From c73c29aea4b563faecac53f0bbe0664c1fc50517 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 12 Mar 2026 14:57:43 -0700 Subject: [PATCH 02/26] PetscOperators: Splitting out PetscMapping Handles the mapping between mesh indices and PETSc global indices. --- include/bout/field3d.hxx | 16 +++- include/bout/petsc_operators.hxx | 48 ++++++++-- include/bout/region.hxx | 17 ++-- src/mesh/petsc_operators.cxx | 107 ++++++++++++----------- tests/unit/mesh/test_petsc_operators.cxx | 58 ++++++++---- 5 files changed, 157 insertions(+), 89 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 54d9fb85c9..66c6d54a9e 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -1,8 +1,8 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2026 BOUT++ contributors + * + * Contact: Ben Dudson, dudson2@llnl.gov * - * Contact: Ben Dudson, bd512@york.ac.uk - * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -348,6 +348,14 @@ public: return std::end(getRegion("RGN_ALL")); }; + /// Convert (x, y, z) to an Ind3D + /// Requires mesh size. + Ind3D indexAt(int x, int y, int z) const { + const int ny = this->getNy(); + const int nz = this->getNz(); + return Ind3D{(((x * ny) + y) * nz) + z, ny, nz}; + } + BoutReal& operator[](const Ind3D& d) { return data[d.ind]; } const BoutReal& operator[](const Ind3D& d) const { return data[d.ind]; } @@ -361,7 +369,7 @@ public: * Direct access to the underlying data array * * If CHECK > 2 then bounds checking is performed - * + * * If CHECK <= 2 then no checks are performed, to * allow inlining and optimisation of inner loops */ diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 6801418996..bda4017c5e 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -1,22 +1,45 @@ /// Represent operators using PETSc matrices /// +#include "bout/field3d.hxx" #include "bout/mesh.hxx" -#include +#include "bout/petsclib.hxx" +#include "bout/region.hxx" #include -class PetscOperators { +#include + +/// Handles mapping between global cell index used in the mesh file, +/// and the global indexing used in PETSc that depends on the mesh +/// decomposition across MPI ranks. +/// +/// Each processor reads a chunk of the global arrays into Field3Ds: +/// - cell_number +/// - forward_cell_number +/// - backward_cell_number +/// Non-negative values in these arrays indicate a cell number. +/// +/// A Petsc Mat is constructed that maps between petsc indices and +/// mesh indices. +class PetscMapping { public: - PetscOperators(Mesh* mesh); + PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number); // Static functions for unit testing - static Region create_region(Field3D cell_number); - static Region create_region_xin(Field3D cell_number); - static Region create_region_xout(Field3D cell_number); + static Region create_region(const Field3D& cell_number); + static Region create_region_xin(const Field3D& cell_number); + static Region create_region_xout(const Field3D& cell_number); + + /// Total number of cells, including X and Y boundaries + unsigned int size() const { + return evolving_region.size() + xin_region.size() + xout_region.size() + + yup_region.size() + ydown_region.size(); + } private: - PetscLib lib{}; // Initialize and finalize PETSc + PetscLib lib; // Initialize and finalize PETSc Region evolving_region; ///< Evolving cells Region xin_region; ///< X boundary cells in inner boundary @@ -25,6 +48,15 @@ private: Region ydown_region; ///< Y boundary cells in backward direction Mat mat_mesh_to_petsc; +}; + +class PetscOperators { +public: + PetscOperators(Mesh* mesh); + +private: + /// Read a Field3D from the mesh or throw BoutException if not found. + Field3D meshGetField3D(Mesh* mesh, const std::string& name); - Mat mat_yforward; + PetscMapping mapping; }; diff --git a/include/bout/region.hxx b/include/bout/region.hxx index e00ad6d41d..d423c3ef4c 100644 --- a/include/bout/region.hxx +++ b/include/bout/region.hxx @@ -1,9 +1,9 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2026 BOUT++ contributors * * Region class by D. Dickinson 2018 * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -170,8 +170,8 @@ struct SpecificInd { int ny = -1, nz = -1; ///< Sizes of y and z dimensions SpecificInd() = default; - SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){}; - explicit SpecificInd(int i) : ind(i){}; + SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz) {}; + explicit SpecificInd(int i) : ind(i) {}; /// Allow explicit conversion to an int explicit operator int() const { return ind; } @@ -491,10 +491,9 @@ template class Region { // Following prevents a Region being created with anything other // than Ind2D, Ind3D or IndPerp as template type - static_assert( - std::is_base_of_v< - Ind2D, T> || std::is_base_of_v || std::is_base_of_v, - "Region must be templated with one of IndPerp, Ind2D or Ind3D"); + static_assert(std::is_base_of_v || std::is_base_of_v + || std::is_base_of_v, + "Region must be templated with one of IndPerp, Ind2D or Ind3D"); public: using data_type = T; @@ -570,7 +569,7 @@ public: }; Region(RegionIndices& indices, int maxregionblocksize = MAXREGIONBLOCKSIZE) - : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)){}; + : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)) {}; // We need to first set the blocks, and only after that call getRegionIndices. // Do not put in the member initialisation diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 6d3f98791e..f5d312d8fb 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -1,8 +1,15 @@ #include "bout/petsc_operators.hxx" #include "bout/boutexception.hxx" +#include "bout/field3d.hxx" +#include "bout/output.hxx" #include "bout/output_bout_types.hxx" +#include "bout/region.hxx" -Region PetscOperators::create_region(Field3D cell_number) { +#include +#include +#include + +Region PetscMapping::create_region(const Field3D& cell_number) { Region::RegionIndices indices; BOUT_FOR_SERIAL(i, cell_number.getRegion("RGN_NOBNDRY")) { if (cell_number[i] > -1) { @@ -12,7 +19,7 @@ Region PetscOperators::create_region(Field3D cell_number) { return Region(indices); } -Region PetscOperators::create_region_xin(Field3D cell_number) { +Region PetscMapping::create_region_xin(const Field3D& cell_number) { Region::RegionIndices xin_indices; const Mesh* mesh = cell_number.getMesh(); if (mesh->firstX()) { @@ -20,7 +27,7 @@ Region PetscOperators::create_region_xin(Field3D cell_number) { for (int j = mesh->ystart; j <= mesh->yend; ++j) { for (int k = mesh->zstart; k <= mesh->zend; ++k) { if (cell_number(i, j, k) > 0) { - xin_indices.push_back(Ind3D{i, j, k}); + xin_indices.push_back(cell_number.indexAt(i, j, k)); } } } @@ -29,7 +36,7 @@ Region PetscOperators::create_region_xin(Field3D cell_number) { return Region(xin_indices); } -Region PetscOperators::create_region_xout(Field3D cell_number) { +Region PetscMapping::create_region_xout(const Field3D& cell_number) { Region::RegionIndices xout_indices; const Mesh* mesh = cell_number.getMesh(); if (mesh->lastX()) { @@ -37,7 +44,7 @@ Region PetscOperators::create_region_xout(Field3D cell_number) { for (int j = mesh->ystart; j <= mesh->yend; ++j) { for (int k = mesh->zstart; k <= mesh->zend; ++k) { if (cell_number(i, j, k) > 0) { - xout_indices.push_back(Ind3D{i, j, k}); + xout_indices.push_back(cell_number.indexAt(i, j, k)); } } } @@ -46,60 +53,28 @@ Region PetscOperators::create_region_xout(Field3D cell_number) { return Region(xout_indices); } -PetscOperators::PetscOperators(Mesh* mesh) { - - // Get the cell center global numbering - // Note: This is an integer but Field3D only stores BoutReal. - Field3D cell_number; - if (mesh->get(cell_number, "cell_number") != 0) { - throw BoutException("PetscOperators requires field 'cell_number'"); - } - - // Forward Y boundary cell numbers. - // >0 only if a forward boundary cell - Field3D forward_cell_number; - if (mesh->get(forward_cell_number, "forward_cell_number") != 0) { - throw BoutException("PetscOperators requires field 'forward_cell_number'"); - } - - // Backward Y boundary cell numbers. - // >0 only if a backward boundary cell - Field3D backward_cell_number; - if (mesh->get(backward_cell_number, "backward_cell_number") != 0) { - throw BoutException("PetscOperators requires field 'backward_cell_number'"); - } - - // Calculate Regions for evolving variables and boundaries - evolving_region = create_region(cell_number); - xin_region = create_region_xin(cell_number); - xout_region = create_region_xout(cell_number); - yup_region = create_region(forward_cell_number); - ydown_region = create_region(backward_cell_number); - - const std::size_t nlocal = evolving_region.size() + xin_region.size() - + xout_region.size() + yup_region.size() - + ydown_region.size(); +PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number) + : evolving_region(create_region(cell_number)), + xin_region(create_region_xin(cell_number)), + xout_region(create_region_xout(cell_number)), + yup_region(create_region(forward_cell_number)), + ydown_region(create_region(backward_cell_number)) { + // Calculate size of each region + const unsigned int nlocal = this->size(); output.write("Local {} : evolving {} xin {} xout {} yup {} ydown {}", nlocal, evolving_region.size(), xin_region.size(), xout_region.size(), yup_region.size(), ydown_region.size()); - int mesh_total_cells; - if (mesh->get(mesh_total_cells, "total_cells") == 0) { - // Check total number of cells - output.write("Total cells in mesh: {}\n", mesh_total_cells); - } - // Create a PETSc matrix // Note: Numbering is different from that used in PetscVector / PetscMatrix // because yup/down boundaries are included. - MPI_Comm comm = BoutComm::get(); - // Renumbering matrix. // Maps PETSc row (or column) indices to the global indices used in // the mesh. This is needed because the PETSc indices depend on the // number of processors. - MatCreate(comm, &mat_mesh_to_petsc); + MatCreate(BoutComm::get(), &mat_mesh_to_petsc); MatSetSizes(mat_mesh_to_petsc, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); MatSetType(mat_mesh_to_petsc, MATMPIAIJ); // Each row will have one non-zero entry, which could be in @@ -113,10 +88,40 @@ PetscOperators::PetscOperators(Mesh* mesh) { PetscInt row = row_start; const PetscScalar ONE = 1.0; - BOUT_FOR_SERIAL(i, evolving_region) { - PetscInt col = ROUND(cell_number[i]); - MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); - ++row; + // Iterate through regions in this order + const std::vector>> regions = { + evolving_region, xin_region, xout_region, yup_region, ydown_region}; + for (const auto& region : regions) { + BOUT_FOR_SERIAL(i, region.get()) { + const PetscInt col = ROUND(cell_number[i]); + MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); + ++row; + } + } + MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); +} + +Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { + Field3D result; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires field '{}'", name); + } + return result; +} + +PetscOperators::PetscOperators(Mesh* mesh) + : mapping(PetscMapping(meshGetField3D(mesh, "cell_number"), + meshGetField3D(mesh, "forward_cell_number"), + meshGetField3D(mesh, "backward_cell_number"))) { + + int mesh_total_cells; + if (mesh->get(mesh_total_cells, "total_cells") == 0) { + // Check total number of cells + if (mesh_total_cells != mapping.size()) { + throw BoutException("Total cells in mesh {} doesn't match mapping size {}", + mesh_total_cells, mapping.size()); + } } //MatCreate(comm, &mat_yforward); diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 949029f19c..b33d41a7b1 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -1,33 +1,40 @@ #include "gtest/gtest.h" +#include "bout/output.hxx" #include "bout/output_bout_types.hxx" #include "bout/petsc_operators.hxx" +#include "bout/region.hxx" #include "fake_mesh_fixture.hxx" #include "test_extras.hxx" -namespace { -Ind3D indexAt(const Field3D& f, int x, int y, int z) { - int ny = f.getNy(); - int nz = f.getNz(); - return Ind3D{(((x * ny) + y) * nz) + z, ny, nz}; -} -} // namespace +// Reuse the "standard" fixture for FakeMesh +using PetscMappingTest = FakeMeshFixture; -using bout::globals::mesh; +TEST_F(PetscMappingTest, create_region_empty) { + const Field3D cell_number{-1.0}; // No cells >= 0 -// Reuse the "standard" fixture for FakeMesh -using PetscOperatorsTest = FakeMeshFixture; + auto rgn = PetscMapping::create_region(cell_number); -TEST_F(PetscOperatorsTest, create_region_empty) { + ASSERT_EQ(0, rgn.size()); +} + +TEST_F(PetscMappingTest, create_region) { Field3D cell_number{-1.0}; // No cells >= 0 - auto rgn = PetscOperators::create_region(cell_number); + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 0; // Corner + cell_number(1, 1, 0) = 1; // In domain + cell_number(0, 1, 1) = 2; // Xin boundary - ASSERT_EQ(0, rgn.size()); + auto rgn = PetscMapping::create_region(cell_number); + ASSERT_EQ(1, rgn.size()); + + const Ind3D expected_ind = cell_number.indexAt(1, 1, 0); + ASSERT_TRUE(expected_ind == *rgn.begin()); } -TEST_F(PetscOperatorsTest, create_region) { +TEST_F(PetscMappingTest, create_region_xin) { Field3D cell_number{-1.0}; // No cells >= 0 // Note: 1 boundary cell in X and Y @@ -35,10 +42,27 @@ TEST_F(PetscOperatorsTest, create_region) { cell_number(1, 1, 0) = 1; // In domain cell_number(0, 1, 1) = 2; // Xin boundary - auto rgn = PetscOperators::create_region(cell_number); + auto rgn = PetscMapping::create_region_xin(cell_number); ASSERT_EQ(1, rgn.size()); - Ind3D expected_ind = indexAt(cell_number, 1, 1, 0); - output.write("{} : {}\n", expected_ind, *rgn.begin()); + const Ind3D expected_ind = cell_number.indexAt(0, 1, 1); + output.write("Expecting {} == {}\n", expected_ind, *rgn.begin()); ASSERT_TRUE(expected_ind == *rgn.begin()); } + +TEST_F(PetscMappingTest, mapping) { + Field3D cell_number{-1.0}; // No cells >= 0 + + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 0; // Corner + cell_number(1, 1, 0) = 1; // In domain + cell_number(0, 1, 1) = 2; // Xin boundary + + const Field3D forward_cell_number{-1.0}; + const Field3D backward_cell_number{-1.0}; + + const PetscMapping mapping(cell_number, forward_cell_number, backward_cell_number); + + // Two cells: one evolving and one in xin + ASSERT_EQ(2, mapping.size()); +} From 724572af852b03ab0f34ecbb7486cd2704f9e322 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 12 Mar 2026 20:49:43 -0700 Subject: [PATCH 03/26] Mesh::get Arrays of int and BoutReal Simple wrapper functions around the Options functionality, so that Array types can be read from the mesh file. --- include/bout/griddata.hxx | 20 ++++++++++++++++++-- include/bout/mesh.hxx | 10 +++++++--- src/mesh/data/gridfromfile.cxx | 23 +++++++++++++++++++++-- src/mesh/mesh.cxx | 14 +++++++++++--- 4 files changed, 57 insertions(+), 10 deletions(-) diff --git a/include/bout/griddata.hxx b/include/bout/griddata.hxx index 29a32e5779..e9d9214846 100644 --- a/include/bout/griddata.hxx +++ b/include/bout/griddata.hxx @@ -2,9 +2,9 @@ * Functions to read input variables on a grid/mesh * ************************************************************************** - * 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++. * @@ -29,12 +29,16 @@ class GridDataSource; #define BOUT_GRIDDATA_H #include "mesh.hxx" +#include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/options.hxx" #include #include +#include +#include + /// Interface class to serve grid data /*! * Provides a generic interface for sources of @@ -71,6 +75,15 @@ public: static constexpr Direction Y = Direction::Y; static constexpr Direction Z = Direction::Z; + virtual bool get([[maybe_unused]] Array& var, + [[maybe_unused]] const std::string& name) { + return false; + } + virtual bool get([[maybe_unused]] Array& var, + [[maybe_unused]] const std::string& name) { + return false; + } + virtual bool get(Mesh* m, std::vector& var, const std::string& name, int len, int offset = 0, Direction dir = GridDataSource::X) = 0; virtual bool get(Mesh* m, std::vector& var, const std::string& name, int len, @@ -115,6 +128,9 @@ public: return getField(m, var, name, def, location); } + bool get(Array& var, const std::string& name) override; + bool get(Array& var, const std::string& name) override; + bool get(Mesh* m, std::vector& var, const std::string& name, int len, int offset = 0, GridDataSource::Direction dir = GridDataSource::X) override; bool get(Mesh* m, std::vector& var, const std::string& name, int len, diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index a1ed6a9011..2422ce283d 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -3,7 +3,7 @@ * * Interface for mesh classes. Contains standard variables and useful * routines. - * + * * Changelog * ========= * @@ -11,7 +11,7 @@ * * Removing coordinate system into separate * Coordinates class * * Adding index derivative functions from derivs.cxx - * + * * 2010-06 Ben Dudson, Sean Farley * * Initial version, adapted from GridData class * * Incorporates code from topology.cpp and Communicator @@ -20,7 +20,7 @@ * Copyright 2010-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 @@ -43,6 +43,7 @@ class Mesh; #ifndef BOUT_MESH_H #define BOUT_MESH_H +#include "bout/array.hxx" #include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" #include "bout/coordinates.hxx" // Coordinates class @@ -170,6 +171,9 @@ public: /// @returns zero if successful, non-zero on failure int get(bool& bval, const std::string& name, bool def = false); + int get(Array& var, const std::string& name); + int get(Array& var, const std::string& name); + /// Get a Field2D from the input source /// including communicating guard cells /// diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 10f6f6926e..e2e9850612 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -14,7 +15,9 @@ #include #include +#include #include +#include GridFile::GridFile(std::string gridfilename) : GridDataSource(true), data(bout::OptionsIO::create(gridfilename)->read()), @@ -120,7 +123,7 @@ bool GridFile::get(Mesh* UNUSED(m), BoutReal& rval, const std::string& name, /*! * Reads a 2D, 3D or FieldPerp field variable from a file - * + * * Successfully reads Field2D or FieldPerp if the variable in the file is 0-D or 2-D. * Successfully reads Field3D if the variable in the file is 0-D, 2-D or 3-D. */ @@ -378,7 +381,7 @@ void GridFile::readField(Mesh* m, const std::string& name, int ys, int yd, int n for (int x = xs; x < xs + nx_to_read; ++x) { for (int y = ys; y < ys + ny_to_read; ++y) { - BoutReal const value = full_var(x, y); + const BoutReal value = full_var(x, y); for (int z = 0; z < var.getNz(); z++) { var(x - xs + xd, y - ys + yd, z) = value; } @@ -466,6 +469,22 @@ void GridFile::readField(Mesh* m, const std::string& name, int UNUSED(ys), int U } } +bool GridFile::get(Array& var, const std::string& name) { + if (not data.isSet(name)) { + return false; + } + var = data[name].as>(); + return true; +} + +bool GridFile::get(Array& var, const std::string& name) { + if (not data.isSet(name)) { + return false; + } + var = data[name].as>(); + return true; +} + bool GridFile::get([[maybe_unused]] Mesh* m, [[maybe_unused]] std::vector& var, [[maybe_unused]] const std::string& name, [[maybe_unused]] int len, [[maybe_unused]] int offset, diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 829c150a15..034da6eecd 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -179,13 +179,21 @@ int Mesh::get(bool& bval, const std::string& name, bool def) { if (source == nullptr) { warn_default_used(def, name); bval = def; - return true; + return 1; } int bval_as_int = 0; - bool success = source->get(this, bval_as_int, name, def); + const bool success = source->get(this, bval_as_int, name, int(def)); bval = bool(bval_as_int); - return !success; + return success ? 0 : 1; +} + +int Mesh::get(Array& var, const std::string& name) { + return source->get(var, name) ? 0 : 1; +} + +int Mesh::get(Array& var, const std::string& name) { + return source->get(var, name) ? 0 : 1; } int Mesh::get(Field2D& var, const std::string& name, BoutReal def, bool communicate, From af65816e753ed4cda2c90f78099a558a847c2120 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 12 Mar 2026 21:46:49 -0700 Subject: [PATCH 04/26] PetscOperators: Implement operators from CSR input Reads row, column and weights from Mesh, constructs Petsc Mat objects. Working on composition and operation on Field3D objects. --- include/bout/petsc_operators.hxx | 87 +++++++++++++++++++++- src/mesh/petsc_operators.cxx | 121 +++++++++++++++++++++++++------ 2 files changed, 184 insertions(+), 24 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index bda4017c5e..c19146a105 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -1,14 +1,21 @@ /// Represent operators using PETSc matrices /// +#include "bout/array.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/mesh.hxx" #include "bout/petsclib.hxx" #include "bout/region.hxx" +#include "bout/utils.hxx" #include +#include +#include #include +#include /// Handles mapping between global cell index used in the mesh file, /// and the global indexing used in PETSc that depends on the mesh @@ -38,25 +45,103 @@ public: + yup_region.size() + ydown_region.size(); } + /// Loops over cells and calls the given function with arguments: + /// - PetscInt PETSc row index + /// - Ind3D Index into Field3D variables + /// - int Global cell number in mesh + /// + template + void map_evolving(Function func) const { + PetscInt row = this->row_start; + BOUT_FOR_SERIAL(i, this->evolving_region) { + func(row, i, ROUND(cell_number[i])); + ++row; + } + } + + /// Note: It doesn't make sense to pass the Ind3D because + /// some reference Field3D, some the yup/ydown fields. + template + void map(Function func) const { + const std::vector>> regions = { + evolving_region, xin_region, xout_region, yup_region, ydown_region}; + PetscInt row = row_start; + for (const auto& region : regions) { + BOUT_FOR_SERIAL(i, region.get()) { + func(row, ROUND(cell_number[i])); + ++row; + } + } + } + private: PetscLib lib; // Initialize and finalize PETSc + Field3D cell_number; + Region evolving_region; ///< Evolving cells Region xin_region; ///< X boundary cells in inner boundary Region xout_region; ///< X boundary cells in outer boundary Region yup_region; ///< Y boundary cells in forward direction Region ydown_region; ///< Y boundary cells in backward direction + PetscInt row_start, row_end; ///< Local row indices Mat mat_mesh_to_petsc; }; +/// Shared pointer to const PetscMapping +/// Mapping is const after creation. +using PetscMappingPtr = std::shared_ptr; + +/// Represents an operator on Field3Ds +class PetscOperator { +public: + /// Create using CSR format arrays + PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, + Array weights); + + /// Perform operation + Field3D operator()(const Field3D& rhs) const; + + /// Operator composition + PetscOperator operator*(const PetscOperator& rhs) const; + + /// Operator addition + PetscOperator operator+(const PetscOperator& rhs) const; + + /// Operator subtraction + PetscOperator operator-(const PetscOperator& rhs) const; + +private: + PetscOperator(PetscMappingPtr mapping, Mat mat) + : mapping(std::move(mapping)), mat_operator(mat) {} + PetscMappingPtr mapping; + Mat mat_operator; +}; + class PetscOperators { public: + /// + /// Reads from the mesh: + /// - cell_number : Field3D + /// - forward_cell_number : Field3D + /// - backward_cell_number : Field3D + /// - total_cells : int, optional PetscOperators(Mesh* mesh); private: /// Read a Field3D from the mesh or throw BoutException if not found. Field3D meshGetField3D(Mesh* mesh, const std::string& name); - PetscMapping mapping; + /// Read a 1D Array from the mesh or throw BoutException + template + Array meshGetArray(Mesh* mesh, const std::string& name) { + Array result; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires Array '{}'", name); + } + return result; + } + + PetscMappingPtr mapping; }; diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index f5d312d8fb..ef48d75dd1 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -1,13 +1,14 @@ #include "bout/petsc_operators.hxx" +#include "bout/array.hxx" +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/output.hxx" #include "bout/output_bout_types.hxx" #include "bout/region.hxx" - -#include +#include #include -#include Region PetscMapping::create_region(const Field3D& cell_number) { Region::RegionIndices indices; @@ -55,7 +56,7 @@ Region PetscMapping::create_region_xout(const Field3D& cell_number) { PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, const Field3D& backward_cell_number) - : evolving_region(create_region(cell_number)), + : cell_number(cell_number), evolving_region(create_region(cell_number)), xin_region(create_region_xin(cell_number)), xout_region(create_region_xout(cell_number)), yup_region(create_region(forward_cell_number)), @@ -82,26 +83,97 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce MatMPIAIJSetPreallocation(mat_mesh_to_petsc, 1, nullptr, 1, nullptr); // Get the range of rows owned by this processor - PetscInt row_start, row_end; MatGetOwnershipRange(mat_mesh_to_petsc, &row_start, &row_end); output.write("Local row range: {} -> {}\n", row_start, row_end); - PetscInt row = row_start; - const PetscScalar ONE = 1.0; // Iterate through regions in this order - const std::vector>> regions = { - evolving_region, xin_region, xout_region, yup_region, ydown_region}; - for (const auto& region : regions) { - BOUT_FOR_SERIAL(i, region.get()) { - const PetscInt col = ROUND(cell_number[i]); - MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); - ++row; - } - } + this->map([&](PetscInt row, PetscInt col) { + // `row` is the PETSc index; `col` is the Mesh index + const PetscScalar ONE = 1.0; + MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); + }); MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); } +PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, + Array weights) + : mapping(mapping) { + + MatCreate(BoutComm::get(), &mat_operator); + int nlocal = mapping->size(); + MatSetSizes(mat_operator, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); + + mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { + if (mesh_index >= rows.size()) { + return; // No weights -> skip + } + // Get the range of indices into columns and weights + int start_ind = rows[mesh_index]; + int end_ind = (mesh_index + 1 < rows.size()) ? rows[mesh_index + 1] : cols.size(); + output.write("Mapping {} : {} -> {}\n", row, start_ind, end_ind); + + MatSetValues(mat_operator, 1, &row, end_ind - start_ind, &cols[start_ind], + &weights[start_ind], INSERT_VALUES); + }); + MatAssemblyBegin(mat_operator, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(mat_operator, MAT_FINAL_ASSEMBLY); +} + +/// Perform operation +Field3D PetscOperator::operator()(const Field3D& rhs) const { + Vec rhs_vec; + VecCreate(BoutComm::get(), &rhs_vec); + VecSetSizes(rhs_vec, this->mapping->size(), PETSC_DETERMINE); + + // Fill vec from rhs + PetscScalar* x; + VecGetArray(rhs_vec, &x); + throw BoutException("Not implemented"); + VecRestoreArray(rhs_vec, &x); + + // Perform Mat-Vec muliplication + Vec result_vec; + VecDuplicate(rhs_vec, &result_vec); + MatMult(this->mat_operator, rhs_vec, result_vec); + + // Copy result_vec into a Field3D + Field3D result; + const PetscScalar* r; + VecGetArrayRead(result_vec, &r); + throw BoutException("Not implemented"); + VecRestoreArrayRead(result_vec, &r); + + return result; +} + +/// Operator composition +PetscOperator PetscOperator::operator*(const PetscOperator& rhs) const { + ASSERT0(this->mapping == rhs.mapping); + Mat mat; + MatMatMult(this->mat_operator, rhs.mat_operator, MAT_INITIAL_MATRIX, PETSC_DETERMINE, + &mat); + return PetscOperator(this->mapping, mat); +} + +/// Operator addition +PetscOperator PetscOperator::operator+(const PetscOperator& rhs) const { + ASSERT0(this->mapping == rhs.mapping); + Mat mat; + MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat); + MatAXPY(mat, 1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN); + return PetscOperator(this->mapping, mat); +} + +/// Operator subtraction +PetscOperator PetscOperator::operator-(const PetscOperator& rhs) const { + ASSERT0(this->mapping == rhs.mapping); + Mat mat; + MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat); + MatAXPY(mat, -1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN); + return PetscOperator(this->mapping, mat); +} + Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { Field3D result; if (mesh->get(result, name) != 0) { @@ -111,19 +183,22 @@ Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { } PetscOperators::PetscOperators(Mesh* mesh) - : mapping(PetscMapping(meshGetField3D(mesh, "cell_number"), - meshGetField3D(mesh, "forward_cell_number"), - meshGetField3D(mesh, "backward_cell_number"))) { + : mapping(std::make_shared( + meshGetField3D(mesh, "cell_number"), + meshGetField3D(mesh, "forward_cell_number"), + meshGetField3D(mesh, "backward_cell_number"))) { int mesh_total_cells; if (mesh->get(mesh_total_cells, "total_cells") == 0) { // Check total number of cells - if (mesh_total_cells != mapping.size()) { + if (mesh_total_cells != mapping->size()) { throw BoutException("Total cells in mesh {} doesn't match mapping size {}", - mesh_total_cells, mapping.size()); + mesh_total_cells, mapping->size()); } } - //MatCreate(comm, &mat_yforward); - //MatSetSizes(mat_yforward, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); + // Read forward operator + PetscOperator forward(mapping, this->meshGetArray(mesh, "forward_rows"), + this->meshGetArray(mesh, "forward_columns"), + this->meshGetArray(mesh, "forward_weights")); } From 8a3af0cc8da3809584fe9ed8ec38203a81606c1b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 13 Mar 2026 09:14:31 -0700 Subject: [PATCH 05/26] PetscOperators: Operate on Field3Ds Map from Field3D, including yup/ydown boundaries, to PETSc Vec, perform matrix-vector product, and map back to a Field3D. --- include/bout/petsc_operators.hxx | 37 ++++++++++++++++++++++++++++++++ src/mesh/petsc_operators.cxx | 33 +++++++++++++++++++++------- 2 files changed, 62 insertions(+), 8 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index c19146a105..72f7476121 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -15,6 +15,7 @@ #include #include #include +#include #include /// Handles mapping between global cell index used in the mesh file, @@ -74,6 +75,41 @@ public: } } + template + void map_local_field(Function func) const { + const std::vector>> regions = { + evolving_region, xin_region, xout_region}; + PetscInt row = 0; // Starting from 0, not row_start + for (const auto& region : regions) { + BOUT_FOR_SERIAL(i, region.get()) { + func(row, i); + ++row; + } + } + } + + template + void map_local_yup(Function func) const { + PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); + BOUT_FOR_SERIAL(i, yup_region) { + func(row, i); + ++row; + } + } + + template + void map_local_ydown(Function func) const { + PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() + + yup_region.size(); + BOUT_FOR_SERIAL(i, ydown_region) { + func(row, i); + ++row; + } + } + + /// Return a matrix that reorders vectors from petsc global index to mesh index + Mat getPetscToMesh() const { return mat_petsc_to_mesh; } + private: PetscLib lib; // Initialize and finalize PETSc @@ -87,6 +123,7 @@ private: PetscInt row_start, row_end; ///< Local row indices Mat mat_mesh_to_petsc; + Mat mat_petsc_to_mesh; }; /// Shared pointer to const PetscMapping diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index ef48d75dd1..2d6c7743ec 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -94,15 +94,19 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce }); MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); + + // Transpose to map petsc indices to mesh indices + MatTranspose(mat_mesh_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_mesh); } PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights) : mapping(mapping) { - MatCreate(BoutComm::get(), &mat_operator); + Mat mat; + MatCreate(BoutComm::get(), &mat); int nlocal = mapping->size(); - MatSetSizes(mat_operator, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); + MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { if (mesh_index >= rows.size()) { @@ -113,11 +117,16 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array {}\n", row, start_ind, end_ind); - MatSetValues(mat_operator, 1, &row, end_ind - start_ind, &cols[start_ind], - &weights[start_ind], INSERT_VALUES); + MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], &weights[start_ind], + INSERT_VALUES); }); - MatAssemblyBegin(mat_operator, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(mat_operator, MAT_FINAL_ASSEMBLY); + MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY); + + // Row indices are PETSc indices but columns are mesh indices + // Multiply on the right by PetscToMesh. + MatMatMult(mat, mapping->getPetscToMesh(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, + &this->mat_operator); } /// Perform operation @@ -128,8 +137,16 @@ Field3D PetscOperator::operator()(const Field3D& rhs) const { // Fill vec from rhs PetscScalar* x; + // Copy rows from Field3D cells VecGetArray(rhs_vec, &x); - throw BoutException("Not implemented"); + mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = rhs[i]; }); + + // Copy Yup / Ydown values from boundaries + const Field3D& yup = rhs.yup(); + mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + + const Field3D& ydown = rhs.yup(); + mapping->map_local_ydown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); VecRestoreArray(rhs_vec, &x); // Perform Mat-Vec muliplication @@ -141,7 +158,7 @@ Field3D PetscOperator::operator()(const Field3D& rhs) const { Field3D result; const PetscScalar* r; VecGetArrayRead(result_vec, &r); - throw BoutException("Not implemented"); + mapping->map_local_field([&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); VecRestoreArrayRead(result_vec, &r); return result; From 2b712591e4ef502562336c560e8f3ca1afc48dbe Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 13 Mar 2026 11:30:53 -0700 Subject: [PATCH 06/26] PetscOperators: Add documentation and error checking Docstrings to explain PetscMapping, PetscOperator and PetscOperators classes. Extend `BOUT_DO_PETSC` macro to print error message, file and line number. --- include/bout/petsc_operators.hxx | 288 ++++++++++++++++++++++++++----- include/bout/petsclib.hxx | 31 ++-- src/mesh/petsc_operators.cxx | 110 ++++++------ 3 files changed, 326 insertions(+), 103 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 72f7476121..3f4d0a596e 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -1,4 +1,12 @@ -/// Represent operators using PETSc matrices +/// Represent finite-difference operators using PETSc matrices. +/// +/// This header defines: +/// - `PetscMapping`, which maps between mesh-global numbering, +/// PETSc numbering, and `Field3D` storage; +/// - `PetscOperator`, which wraps a PETSc `Mat` as a linear operator on +/// `Field3D` objects; and +/// - `PetscOperators`, a helper for reading mappings and operators from +/// a `Mesh`. /// #include "bout/array.hxx" @@ -18,39 +26,84 @@ #include #include -/// Handles mapping between global cell index used in the mesh file, -/// and the global indexing used in PETSc that depends on the mesh -/// decomposition across MPI ranks. +/// Map between mesh-global cell numbering and PETSc numbering. +/// +/// The mesh file stores operator stencils and weights using a mesh-global +/// numbering scheme. This numbering includes: +/// - evolving cells in the main domain; +/// - X-boundary cells; and +/// - Y-boundary cells represented by the `yup` and `ydown` fields. /// -/// Each processor reads a chunk of the global arrays into Field3Ds: -/// - cell_number -/// - forward_cell_number -/// - backward_cell_number -/// Non-negative values in these arrays indicate a cell number. +/// PETSc uses a different global numbering based on the MPI distribution +/// of rows and vector entries. `PetscMapping` records the correspondence +/// between these numbering systems and provides helper methods to iterate +/// over the local cells in a consistent order. /// -/// A Petsc Mat is constructed that maps between petsc indices and -/// mesh indices. +/// The mapping order used by this class is: +/// 1. evolving cells; +/// 2. inner X-boundary cells; +/// 3. outer X-boundary cells; +/// 4. `yup` boundary cells; and +/// 5. `ydown` boundary cells. +/// +/// A pair of PETSc matrices is constructed internally to reorder vectors +/// between mesh-global and PETSc-global numbering. class PetscMapping { public: + /// Construct the mapping from mesh numbering fields. + /// + /// @param cell_number Mesh-global cell numbers for evolving cells. + /// @param forward_cell_number Mesh-global cell numbers associated with + /// forward Y-boundary (`yup`) values. + /// @param backward_cell_number Mesh-global cell numbers associated with + /// backward Y-boundary (`ydown`) values. + /// + /// Non-negative values in the numbering fields identify valid cells. + /// The constructor builds the local regions and the PETSc permutation + /// matrices needed to translate between mesh and PETSc numbering. PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, const Field3D& backward_cell_number); - // Static functions for unit testing + /// Destroy the PETSc matrices owned by this mapping. + ~PetscMapping(); + + /// Build a region containing all valid cells in a numbering field. + /// + /// @param cell_number Field whose entries contain mesh-global cell numbers. + /// @return A region containing all indices with non-negative cell numbers. + /// + /// This helper is public for unit testing. static Region create_region(const Field3D& cell_number); + + /// Build a region containing valid cells on the inner X boundary. + /// + /// @param cell_number Field whose entries contain mesh-global cell numbers. + /// @return A region containing valid cells adjacent to the inner X boundary. static Region create_region_xin(const Field3D& cell_number); + + /// Build a region containing valid cells on the outer X boundary. + /// + /// @param cell_number Field whose entries contain mesh-global cell numbers. + /// @return A region containing valid cells adjacent to the outer X boundary. static Region create_region_xout(const Field3D& cell_number); - /// Total number of cells, including X and Y boundaries + /// Return the total number of locally mapped cells. + /// + /// The count includes evolving cells and all boundary cells represented + /// in the mapping on this processor. unsigned int size() const { return evolving_region.size() + xin_region.size() + xout_region.size() + yup_region.size() + ydown_region.size(); } - /// Loops over cells and calls the given function with arguments: - /// - PetscInt PETSc row index - /// - Ind3D Index into Field3D variables - /// - int Global cell number in mesh + /// Iterate over locally owned evolving cells in PETSc row order. + /// + /// @tparam Function Callable with signature + /// `func(PetscInt petsc_row, Ind3D index, int mesh_global_index)`. + /// @param func Function called once for each evolving cell owned locally. /// + /// The PETSc row index starts at the locally owned global row `row_start`. + /// The supplied `Ind3D` index refers to the main `Field3D` storage. template void map_evolving(Function func) const { PetscInt row = this->row_start; @@ -60,8 +113,17 @@ public: } } - /// Note: It doesn't make sense to pass the Ind3D because - /// some reference Field3D, some the yup/ydown fields. + /// Iterate over all locally mapped cells in PETSc row order. + /// + /// @tparam Function Callable with signature + /// `func(PetscInt petsc_row, int mesh_global_index)`. + /// @param func Function called once for each mapped cell owned locally. + /// + /// The iteration order is evolving, inner X boundary, outer X boundary, + /// `yup`, then `ydown`. + /// + /// No `Ind3D` is passed because some entries refer to interior `Field3D` + /// storage while others refer to boundary storage in `yup` or `ydown`. template void map(Function func) const { const std::vector>> regions = { @@ -75,6 +137,15 @@ public: } } + /// Iterate over locally stored entries packed from the main `Field3D`. + /// + /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. + /// @param func Function called once for each evolving or X-boundary cell + /// stored in the local PETSc vector layout. + /// + /// The index passed to `func` is a local vector index beginning at zero. + /// Only the main field storage is traversed here; Y-boundary values are + /// handled separately by `map_local_yup()` and `map_local_ydown()`. template void map_local_field(Function func) const { const std::vector>> regions = { @@ -88,6 +159,12 @@ public: } } + /// Iterate over local vector entries corresponding to `yup` boundary values. + /// + /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. + /// @param func Function called once for each locally stored `yup` entry. + /// + /// The local index is the offset into the packed local PETSc vector. template void map_local_yup(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); @@ -97,6 +174,12 @@ public: } } + /// Iterate over local vector entries corresponding to `ydown` boundary values. + /// + /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. + /// @param func Function called once for each locally stored `ydown` entry. + /// + /// The local index is the offset into the packed local PETSc vector. template void map_local_ydown(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() @@ -107,72 +190,185 @@ public: } } - /// Return a matrix that reorders vectors from petsc global index to mesh index + /// Return the PETSc matrix that maps PETSc ordering to mesh ordering. + /// + /// @return PETSc matrix representing the permutation from PETSc global + /// indices to mesh-global indices. Mat getPetscToMesh() const { return mat_petsc_to_mesh; } private: PetscLib lib; // Initialize and finalize PETSc + /// Mesh-global numbering for cells stored on this rank. Field3D cell_number; - Region evolving_region; ///< Evolving cells - Region xin_region; ///< X boundary cells in inner boundary - Region xout_region; ///< X boundary cells in outer boundary - Region yup_region; ///< Y boundary cells in forward direction - Region ydown_region; ///< Y boundary cells in backward direction + /// Evolving cells in the main domain. + Region evolving_region; - PetscInt row_start, row_end; ///< Local row indices + /// X-boundary cells on the inner radial boundary. + Region xin_region; + + /// X-boundary cells on the outer radial boundary. + Region xout_region; + + /// Y-boundary cells in the forward direction. + Region yup_region; + + /// Y-boundary cells in the backward direction. + Region ydown_region; + + /// First and one-past-last PETSc global rows owned by this rank. + PetscInt row_start, row_end; + + /// Permutation matrix mapping mesh-global ordering to PETSc ordering. Mat mat_mesh_to_petsc; + + /// Permutation matrix mapping PETSc ordering to mesh-global ordering. Mat mat_petsc_to_mesh; }; -/// Shared pointer to const PetscMapping -/// Mapping is const after creation. +/// Shared pointer to an immutable `PetscMapping`. +/// +/// The mapping is constructed once and then shared by operators that use it. using PetscMappingPtr = std::shared_ptr; -/// Represents an operator on Field3Ds +/// Linear operator on `Field3D` backed by a PETSc matrix. +/// +/// `PetscOperator` wraps a PETSc `Mat` whose rows and columns use the PETSc +/// numbering defined by a corresponding `PetscMapping`. Operators are usually +/// constructed from CSR arrays read from the mesh file and can then be: +/// - applied to a `Field3D`; +/// - composed with other operators; and +/// - added or subtracted. class PetscOperator { public: - /// Create using CSR format arrays + /// Construct an operator from CSR data in mesh-global numbering. + /// + /// @param mapping Shared mapping between mesh and PETSc numbering. + /// @param rows CSR row offsets. + /// @param cols CSR column indices in mesh-global numbering. + /// @param weights CSR non-zero values. + /// + /// The CSR arrays define an operator in mesh-global numbering. The + /// constructor converts this representation into a PETSc matrix using + /// the supplied mapping. PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights); - /// Perform operation + /// Destroy the PETSc matrix and working vectors owned by this operator. + ~PetscOperator(); + + /// Apply the operator to a field. + /// + /// @param rhs Input field, including any required `yup` and `ydown` + /// boundary values. + /// @return Result of applying the operator to `rhs`. + /// + /// The input field is packed into PETSc vector storage using the mapping, + /// multiplied by the PETSc matrix, then unpacked back into a `Field3D`. Field3D operator()(const Field3D& rhs) const; - /// Operator composition + /// Compose this operator with another operator. + /// + /// @param rhs Right-hand-side operator. + /// @return The composed operator `(*this) * rhs`. + /// + /// Both operators must use compatible mappings. PetscOperator operator*(const PetscOperator& rhs) const; - /// Operator addition + /// Add two operators. + /// + /// @param rhs Right-hand-side operator. + /// @return The sum of the two operators. + /// + /// Both operators must use compatible mappings. PetscOperator operator+(const PetscOperator& rhs) const; - /// Operator subtraction + /// Subtract one operator from another. + /// + /// @param rhs Right-hand-side operator. + /// @return The difference of the two operators. + /// + /// Both operators must use compatible mappings. PetscOperator operator-(const PetscOperator& rhs) const; private: + /// Construct directly from an existing PETSc matrix. + /// + /// @param mapping Shared mapping between mesh and PETSc numbering. + /// @param mat PETSc matrix implementing the operator. + /// + /// This constructor is used internally when creating operators from + /// PETSc matrix algebra such as composition, addition, or subtraction. PetscOperator(PetscMappingPtr mapping, Mat mat) - : mapping(std::move(mapping)), mat_operator(mat) {} + : mapping(std::move(mapping)), mat_operator(mat) { + // Allocate working vectors + VecCreate(BoutComm::get(), &this->rhs_vec); + VecSetSizes(this->rhs_vec, this->mapping->size(), PETSC_DETERMINE); + VecDuplicate(this->rhs_vec, &this->result_vec); + } + + /// Shared mapping between mesh and PETSc numbering. PetscMappingPtr mapping; + + /// PETSc matrix implementing the operator. Mat mat_operator; + + /// Work vector holding the packed input field. + Vec rhs_vec; + + /// Work vector holding the packed result. + Vec result_vec; }; +/// Collection of PETSc operators read from a mesh. +/// +/// `PetscOperators` owns the `PetscMapping` derived from the mesh numbering +/// fields and provides access to named operators stored in the mesh file +/// as CSR arrays. class PetscOperators { public: + /// Construct PETSc operator support from a mesh. + /// + /// Reads the following mesh variables: + /// - `cell_number` + /// - `forward_cell_number` + /// - `backward_cell_number` /// - /// Reads from the mesh: - /// - cell_number : Field3D - /// - forward_cell_number : Field3D - /// - backward_cell_number : Field3D - /// - total_cells : int, optional + /// and uses them to construct the shared `PetscMapping`. + /// + /// @param mesh Mesh from which numbering fields and operators are read. PetscOperators(Mesh* mesh); + /// Retrieve a named PETSc operator from the mesh. + /// + /// @param name Base name of the operator in the mesh file. + /// @return A `PetscOperator` constructed from the arrays + /// `{name}_rows`, `{name}_columns`, and `{name}_weights`. + /// + /// Throws `BoutException` if any of the required arrays are missing. + PetscOperator get(const std::string& name) const; + private: - /// Read a Field3D from the mesh or throw BoutException if not found. - Field3D meshGetField3D(Mesh* mesh, const std::string& name); + /// Read a `Field3D` from the mesh. + /// + /// @param mesh Mesh from which to read. + /// @param name Name of the field. + /// @return The requested `Field3D`. + /// + /// Throws `BoutException` if the field is not present. + Field3D meshGetField3D(Mesh* mesh, const std::string& name) const; - /// Read a 1D Array from the mesh or throw BoutException + /// Read a one-dimensional `Array` from the mesh. + /// + /// @tparam T Array element type. + /// @param mesh Mesh from which to read. + /// @param name Name of the array. + /// @return The requested array. + /// + /// Throws `BoutException` if the array is not present. template - Array meshGetArray(Mesh* mesh, const std::string& name) { + Array meshGetArray(Mesh* mesh, const std::string& name) const { Array result; if (mesh->get(result, name) != 0) { throw BoutException("PetscOperators requires Array '{}'", name); @@ -180,5 +376,9 @@ private: return result; } + /// Mesh from which operator data are read. + Mesh* mesh; + + /// Shared mapping used by all operators created from this mesh. PetscMappingPtr mapping; }; diff --git a/include/bout/petsclib.hxx b/include/bout/petsclib.hxx index 41d7618fd2..637943bd12 100644 --- a/include/bout/petsclib.hxx +++ b/include/bout/petsclib.hxx @@ -1,31 +1,31 @@ /*!************************************************************************ * Provides access to the PETSc library, handling initialisation and - * finalisation. + * finalisation. * * Usage * ----- * * #include - * + * * class MyClass { * public: - * + * * private: * PetscLib lib; * }; - * + * * * This will then automatically initialise Petsc the first time an object * is created, and finalise it when the last object is destroyed. * * This header tries to workaround some annoying PETSc features, and * so it *must* be included before *any* PETSc header. - * + * ************************************************************************** - * Copyright 2012 - 2025 BOUT++ contributors + * Copyright 2012 - 2026 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 @@ -48,6 +48,8 @@ #include "bout/build_defines.hxx" +#include + class Options; #if BOUT_HAS_PETSC @@ -66,7 +68,7 @@ class Options; #include "bout/boutexception.hxx" // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define BOUT_DO_PETSC(cmd) PetscLib::assertIerr((cmd), #cmd) +#define BOUT_DO_PETSC(cmd) PetscLib::assertIerr((cmd), #cmd, __FILE__, __LINE__) /*! * Handles initialisation and finalisation of PETSc library. @@ -124,10 +126,17 @@ public: static void cleanup(); static inline void assertIerr(PetscErrorCode ierr, - const std::string& petsc_op = "PETSc operation") { - if (ierr != 0) { - throw BoutException("{:s} failed with {:d}", petsc_op, ierr); + const std::string& petsc_op = "PETSc operation", + const char* file = nullptr, int linenr = 0) { + if (PetscLikely(ierr == PETSC_SUCCESS)) { + return; } + + const char* err_msg = nullptr; + PetscErrorMessage(ierr, &err_msg, nullptr); + throw BoutException("{:s} failed at {}:{} with {:d} [{}]", petsc_op, + (file != nullptr) ? file : "unknown_file", linenr, ierr, + (err_msg != nullptr) ? err_msg : ""); } static BoutException SNESFailure(SNES& snes); diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 2d6c7743ec..8df2f4eeb7 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -6,9 +6,11 @@ #include "bout/field3d.hxx" #include "bout/output.hxx" #include "bout/output_bout_types.hxx" +#include "bout/petsclib.hxx" #include "bout/region.hxx" #include #include +#include Region PetscMapping::create_region(const Field3D& cell_number) { Region::RegionIndices indices; @@ -27,7 +29,7 @@ Region PetscMapping::create_region_xin(const Field3D& cell_number) { for (int i = 0; i < mesh->xstart; ++i) { for (int j = mesh->ystart; j <= mesh->yend; ++j) { for (int k = mesh->zstart; k <= mesh->zend; ++k) { - if (cell_number(i, j, k) > 0) { + if (cell_number(i, j, k) > -1) { xin_indices.push_back(cell_number.indexAt(i, j, k)); } } @@ -44,7 +46,7 @@ Region PetscMapping::create_region_xout(const Field3D& cell_number) { for (int i = mesh->xend + 1; i < mesh->LocalNx; ++i) { for (int j = mesh->ystart; j <= mesh->yend; ++j) { for (int k = mesh->zstart; k <= mesh->zend; ++k) { - if (cell_number(i, j, k) > 0) { + if (cell_number(i, j, k) > -1) { xout_indices.push_back(cell_number.indexAt(i, j, k)); } } @@ -75,38 +77,45 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce // Maps PETSc row (or column) indices to the global indices used in // the mesh. This is needed because the PETSc indices depend on the // number of processors. - MatCreate(BoutComm::get(), &mat_mesh_to_petsc); - MatSetSizes(mat_mesh_to_petsc, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); - MatSetType(mat_mesh_to_petsc, MATMPIAIJ); + BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat_mesh_to_petsc)); + BOUT_DO_PETSC( + MatSetSizes(mat_mesh_to_petsc, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); + BOUT_DO_PETSC(MatSetType(mat_mesh_to_petsc, MATMPIAIJ)); // Each row will have one non-zero entry, which could be in // either the "diagonal" or "off-diagonal" block. - MatMPIAIJSetPreallocation(mat_mesh_to_petsc, 1, nullptr, 1, nullptr); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_mesh_to_petsc, 1, nullptr, 1, nullptr)); // Get the range of rows owned by this processor - MatGetOwnershipRange(mat_mesh_to_petsc, &row_start, &row_end); + BOUT_DO_PETSC(MatGetOwnershipRange(mat_mesh_to_petsc, &row_start, &row_end)); output.write("Local row range: {} -> {}\n", row_start, row_end); // Iterate through regions in this order this->map([&](PetscInt row, PetscInt col) { // `row` is the PETSc index; `col` is the Mesh index const PetscScalar ONE = 1.0; - MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES); + BOUT_DO_PETSC(MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES)); }); - MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY); + BOUT_DO_PETSC(MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY)); // Transpose to map petsc indices to mesh indices MatTranspose(mat_mesh_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_mesh); } +PetscMapping::~PetscMapping() { + MatDestroy(&this->mat_mesh_to_petsc); + MatDestroy(&this->mat_petsc_to_mesh); +} + PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights) - : mapping(mapping) { + : mapping(std::move(mapping)) { Mat mat; - MatCreate(BoutComm::get(), &mat); + BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); + BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); int nlocal = mapping->size(); - MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE); + BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { if (mesh_index >= rows.size()) { @@ -115,51 +124,55 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array {}\n", row, start_ind, end_ind); - MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], &weights[start_ind], - INSERT_VALUES); + BOUT_DO_PETSC(MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], + &weights[start_ind], INSERT_VALUES)); }); - MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY); + BOUT_DO_PETSC(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); // Row indices are PETSc indices but columns are mesh indices // Multiply on the right by PetscToMesh. - MatMatMult(mat, mapping->getPetscToMesh(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, - &this->mat_operator); + BOUT_DO_PETSC(MatMatMult(mat, mapping->getPetscToMesh(), MAT_INITIAL_MATRIX, + PETSC_DETERMINE, &this->mat_operator)); + + // Allocate vectors for operator input and output + BOUT_DO_PETSC(VecCreate(BoutComm::get(), &this->rhs_vec)); + BOUT_DO_PETSC(VecSetSizes(this->rhs_vec, this->mapping->size(), PETSC_DETERMINE)); + BOUT_DO_PETSC(VecDuplicate(this->rhs_vec, &this->result_vec)); +} + +PetscOperator::~PetscOperator() { + MatDestroy(&this->mat_operator); + VecDestroy(&this->rhs_vec); + VecDestroy(&this->result_vec); } /// Perform operation Field3D PetscOperator::operator()(const Field3D& rhs) const { - Vec rhs_vec; - VecCreate(BoutComm::get(), &rhs_vec); - VecSetSizes(rhs_vec, this->mapping->size(), PETSC_DETERMINE); - // Fill vec from rhs PetscScalar* x; // Copy rows from Field3D cells - VecGetArray(rhs_vec, &x); + BOUT_DO_PETSC(VecGetArray(this->rhs_vec, &x)); mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = rhs[i]; }); // Copy Yup / Ydown values from boundaries const Field3D& yup = rhs.yup(); mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); - const Field3D& ydown = rhs.yup(); + const Field3D& ydown = rhs.ydown(); mapping->map_local_ydown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); - VecRestoreArray(rhs_vec, &x); + BOUT_DO_PETSC(VecRestoreArray(this->rhs_vec, &x)); // Perform Mat-Vec muliplication - Vec result_vec; - VecDuplicate(rhs_vec, &result_vec); - MatMult(this->mat_operator, rhs_vec, result_vec); + BOUT_DO_PETSC(MatMult(this->mat_operator, rhs_vec, result_vec)); // Copy result_vec into a Field3D - Field3D result; + Field3D result{emptyFrom(rhs)}; // This allocates memory const PetscScalar* r; - VecGetArrayRead(result_vec, &r); + BOUT_DO_PETSC(VecGetArrayRead(result_vec, &r)); mapping->map_local_field([&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); - VecRestoreArrayRead(result_vec, &r); + BOUT_DO_PETSC(VecRestoreArrayRead(result_vec, &r)); return result; } @@ -168,8 +181,8 @@ Field3D PetscOperator::operator()(const Field3D& rhs) const { PetscOperator PetscOperator::operator*(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); Mat mat; - MatMatMult(this->mat_operator, rhs.mat_operator, MAT_INITIAL_MATRIX, PETSC_DETERMINE, - &mat); + BOUT_DO_PETSC(MatMatMult(this->mat_operator, rhs.mat_operator, MAT_INITIAL_MATRIX, + PETSC_DETERMINE, &mat)); return PetscOperator(this->mapping, mat); } @@ -177,8 +190,8 @@ PetscOperator PetscOperator::operator*(const PetscOperator& rhs) const { PetscOperator PetscOperator::operator+(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); Mat mat; - MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat); - MatAXPY(mat, 1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN); + BOUT_DO_PETSC(MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat)); + BOUT_DO_PETSC(MatAXPY(mat, 1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); return PetscOperator(this->mapping, mat); } @@ -186,12 +199,12 @@ PetscOperator PetscOperator::operator+(const PetscOperator& rhs) const { PetscOperator PetscOperator::operator-(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); Mat mat; - MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat); - MatAXPY(mat, -1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN); + BOUT_DO_PETSC(MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat)); + BOUT_DO_PETSC(MatAXPY(mat, -1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); return PetscOperator(this->mapping, mat); } -Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { +Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const { Field3D result; if (mesh->get(result, name) != 0) { throw BoutException("PetscOperators requires field '{}'", name); @@ -200,10 +213,10 @@ Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { } PetscOperators::PetscOperators(Mesh* mesh) - : mapping(std::make_shared( - meshGetField3D(mesh, "cell_number"), - meshGetField3D(mesh, "forward_cell_number"), - meshGetField3D(mesh, "backward_cell_number"))) { + : mesh(mesh), mapping(std::make_shared( + meshGetField3D(mesh, "cell_number"), + meshGetField3D(mesh, "forward_cell_number"), + meshGetField3D(mesh, "backward_cell_number"))) { int mesh_total_cells; if (mesh->get(mesh_total_cells, "total_cells") == 0) { @@ -213,9 +226,10 @@ PetscOperators::PetscOperators(Mesh* mesh) mesh_total_cells, mapping->size()); } } +} - // Read forward operator - PetscOperator forward(mapping, this->meshGetArray(mesh, "forward_rows"), - this->meshGetArray(mesh, "forward_columns"), - this->meshGetArray(mesh, "forward_weights")); +PetscOperator PetscOperators::get(const std::string& name) const { + return PetscOperator(this->mapping, this->meshGetArray(this->mesh, name + "_rows"), + this->meshGetArray(this->mesh, name + "_columns"), + this->meshGetArray(this->mesh, name + "_weights")); } From e8f562f5ff31fa4a79b12cfb7525ca32bd6740a1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 13 Mar 2026 22:51:32 -0700 Subject: [PATCH 07/26] Array: Add static fromValues(initializer_list) method Alternative to initializer_list constructor that avoids confusion with e.g. `Array{2}` : Is it an array of size 2, or an array of size 1 containing value 2? `std::vector` uses `vector(2)` means size 2; `vector{2}` means value 2. With this commit, both `Array{2}` and `Array(2)` mean size 2; `Array::fromValues({2})` means size 1, value 2. It's more verbose but less prone to typos. --- include/bout/array.hxx | 20 ++++++++++++-- src/invert/fft_fftw.cxx | 15 ++++++---- src/solver/impls/imex-bdf2/imex-bdf2.cxx | 25 ++++++++--------- tests/unit/include/bout/test_array.cxx | 29 ++++++++++++-------- tests/unit/include/test_cyclic_reduction.cxx | 13 +++++---- tests/unit/invert/test_fft.cxx | 8 +++--- 6 files changed, 66 insertions(+), 44 deletions(-) diff --git a/include/bout/array.hxx b/include/bout/array.hxx index 82677fd38a..d8ce207dac 100644 --- a/include/bout/array.hxx +++ b/include/bout/array.hxx @@ -27,6 +27,7 @@ #define BOUT_ARRAY_H #include +#include #include #include #include @@ -187,9 +188,22 @@ public: Array() noexcept : ptr(nullptr) {} /*! - * Create an array of given length + * Create an array of given length. */ - Array(size_type len) { ptr = get(len); } + Array(size_type len) : ptr(get(len)) {} + + /*! + * Create an array with initializer list. + * This is explicit to avoid confusion with (size_type) constructor. + */ + static Array fromValues(std::initializer_list init) { + if (init.size() == 0) { + return Array(); + } + Array array(init.size()); + std::copy(init.begin(), init.end(), array.begin()); + return array; + } /*! * Destructor. Releases the underlying dataBlock @@ -199,7 +213,7 @@ public: /*! * Copy constructor */ - Array(const Array& other) noexcept { ptr = other.ptr; } + Array(const Array& other) noexcept : ptr(other.ptr) {} /*! * Assignment operator diff --git a/src/invert/fft_fftw.cxx b/src/invert/fft_fftw.cxx index 2c633da6d5..2ada5bc98a 100644 --- a/src/invert/fft_fftw.cxx +++ b/src/invert/fft_fftw.cxx @@ -4,9 +4,9 @@ * \brief FFT routines using external libraries * ************************************************************************** - * 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++. * @@ -25,9 +25,12 @@ * **************************************************************************/ -#include "bout/build_defines.hxx" - +#include +#include +#include +#include #include +#include #include #include #include @@ -518,7 +521,7 @@ Array rfft(const Array& in) { ASSERT1(!in.empty()); int size{in.size()}; - Array out{(size / 2) + 1}; + Array out((size / 2) + 1); bout::fft::rfft(in.begin(), size, out.begin()); return out; @@ -528,7 +531,7 @@ Array irfft(const Array& in, int length) { ASSERT1(!in.empty()); ASSERT1(in.size() == (length / 2) + 1); - Array out{length}; + Array out(length); bout::fft::irfft(in.begin(), length, out.begin()); return out; diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 4ccf27981b..56d21c91d0 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -4,11 +4,15 @@ #include "imex-bdf2.hxx" +#include #include +#include #include #include #include +#include #include +#include #include #include @@ -119,12 +123,9 @@ static PetscErrorCode FormFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, } static PetscErrorCode imexbdf2PCapply(PC pc, Vec x, Vec y) { - int ierr; - // Get the context IMEXBDF2* s; - ierr = PCShellGetContext(pc, reinterpret_cast(&s)); - CHKERRQ(ierr); + PetscCall(PCShellGetContext(pc, reinterpret_cast(&s))); PetscFunctionReturn(s->precon(x, y)); } @@ -197,8 +198,8 @@ int IMEXBDF2::init() { // Allocate memory and initialise structures u.reallocate(nlocal); for (int i = 0; i < maxOrder; i++) { - uV.emplace_back(Array{nlocal}); - fV.emplace_back(Array{nlocal}); + uV.emplace_back(Array(nlocal)); + fV.emplace_back(Array(nlocal)); timesteps.push_back(timestep); uFac.push_back(0.0); fFac.push_back(0.0); @@ -220,16 +221,12 @@ int IMEXBDF2::init() { } // Initialise PETSc components - int ierr; // Vectors - ierr = VecCreate(BoutComm::get(), &snes_x); - CHKERRQ(ierr); - ierr = VecSetSizes(snes_x, nlocal, PETSC_DECIDE); - CHKERRQ(ierr); - ierr = VecSetFromOptions(snes_x); - CHKERRQ(ierr); - VecDuplicate(snes_x, &snes_f); + PetscCall(VecCreate(BoutComm::get(), &snes_x)); + PetscCall(VecSetSizes(snes_x, nlocal, PETSC_DECIDE)); + PetscCall(VecSetFromOptions(snes_x)); + PetscCall(VecDuplicate(snes_x, &snes_f)); // The SNES solver object(s) constructSNES(&snes); diff --git a/tests/unit/include/bout/test_array.cxx b/tests/unit/include/bout/test_array.cxx index 21be6fe8a5..4efd968f57 100644 --- a/tests/unit/include/bout/test_array.cxx +++ b/tests/unit/include/bout/test_array.cxx @@ -1,9 +1,7 @@ #include "gtest/gtest.h" #include "bout/array.hxx" -#include "bout/boutexception.hxx" -#include #include // In order to keep these tests independent, they need to use @@ -19,7 +17,7 @@ class ArrayTest : public ::testing::Test { }; TEST_F(ArrayTest, ArraySize) { - Array a(5); + const Array a(5); ASSERT_FALSE(a.empty()); EXPECT_EQ(a.size(), 5); @@ -36,7 +34,7 @@ TEST_F(ArrayTest, ArrayValues) { } TEST_F(ArrayTest, CopyArrayConstructor) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -53,7 +51,7 @@ TEST_F(ArrayTest, CopyArrayConstructor) { } TEST_F(ArrayTest, CopyArrayOperator) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -71,7 +69,7 @@ TEST_F(ArrayTest, CopyArrayOperator) { } TEST_F(ArrayTest, CopyArrayNonMemberFunction) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -86,11 +84,11 @@ TEST_F(ArrayTest, CopyArrayNonMemberFunction) { } TEST_F(ArrayTest, SwapArray) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); - Array b{10}; + Array b(10); std::iota(b.begin(), b.end(), 5); @@ -108,7 +106,7 @@ TEST_F(ArrayTest, SwapArray) { } TEST_F(ArrayTest, MoveArrayConstructor) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -123,7 +121,7 @@ TEST_F(ArrayTest, MoveArrayConstructor) { } TEST_F(ArrayTest, Reallocate) { - Array a{}; + Array a; ASSERT_TRUE(a.empty()); @@ -235,7 +233,7 @@ TEST_F(ArrayTest, RetrieveData) { TEST_F(ArrayTest, Assignment) { Array a(35); - Array b(35); + const Array b(35); // Assign a = b; @@ -251,3 +249,12 @@ TEST_F(ArrayTest, OutOfBoundsThrow) { EXPECT_THROW(a[-1] = 1.0, BoutException); } #endif + +TEST_F(ArrayTest, FromValues) { + const auto a = Array::fromValues({3, 1, 4, 2}); + ASSERT_EQ(4, a.size()); + EXPECT_EQ(3, a[0]); + EXPECT_EQ(1, a[1]); + EXPECT_EQ(4, a[2]); + EXPECT_EQ(2, a[3]); +} diff --git a/tests/unit/include/test_cyclic_reduction.cxx b/tests/unit/include/test_cyclic_reduction.cxx index e057c41bef..5eb9acff0a 100644 --- a/tests/unit/include/test_cyclic_reduction.cxx +++ b/tests/unit/include/test_cyclic_reduction.cxx @@ -2,6 +2,7 @@ #include "test_extras.hxx" #include "bout/array.hxx" +#include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" #include "bout/cyclic_reduction.hxx" @@ -19,7 +20,7 @@ Array makeArrayFromVector(const std::vector& values) { using std::begin; using std::end; - Array array{static_cast::size_type>(values.size())}; + Array array(static_cast::size_type>(values.size())); std::copy(begin(values), end(values), begin(array)); return array; } @@ -28,9 +29,9 @@ Matrix makeMatrixFromVector(const std::vector>& using std::begin; using std::end; - Matrix matrix{static_cast::size_type>(values.size()), - static_cast::size_type>(values[0].size())}; - auto start = begin(matrix); + Matrix matrix(static_cast::size_type>(values.size()), + static_cast::size_type>(values[0].size())); + auto* start = begin(matrix); for (const auto& sub_values : values) { start = std::copy(begin(sub_values), end(sub_values), start); } @@ -48,7 +49,7 @@ TEST(CyclicReduction, SerialSolveSingleArray) { reduce.setCoefs(a, b, c); auto rhs = makeArrayFromVector({0., 1., 2., 2., 3.}); - Array x{reduction_size}; + Array x(reduction_size); reduce.solve(rhs, x); @@ -70,7 +71,7 @@ TEST(CyclicReduction, SerialSolveSingleMatrix) { reduce.setCoefs(a, b, c); auto rhs = makeMatrixFromVector({{0., 1., 2., 2., 3.}}); - Matrix x{1, reduction_size}; + Matrix x(1, reduction_size); reduce.solve(rhs, x); diff --git a/tests/unit/invert/test_fft.cxx b/tests/unit/invert/test_fft.cxx index 1fb30c5eb7..36b344dbb1 100644 --- a/tests/unit/invert/test_fft.cxx +++ b/tests/unit/invert/test_fft.cxx @@ -4,12 +4,12 @@ #include "test_extras.hxx" #include "bout/array.hxx" +#include "bout/bout_types.hxx" #include "bout/constants.hxx" #include "bout/dcomplex.hxx" #include "bout/fft.hxx" #include -#include #include #if BOUT_HAS_FFTW @@ -19,7 +19,7 @@ class FFTTest : public ::testing::TestWithParam { : size(GetParam()), nmodes((size / 2) + 1), real_signal(size), fft_signal(nmodes) { // Make grid indices from [0, size - 1] - Array indices{size}; + Array indices(size); std::iota(indices.begin(), indices.end(), 0.0); // Calculate sin(x) + cos(2x) on [0, 2pi] @@ -48,7 +48,7 @@ INSTANTIATE_TEST_SUITE_P(FFTEvenAndOddSamples, FFTTest, ::testing::Values(8, 9)) TEST_P(FFTTest, rfft) { - Array output{nmodes}; + Array output(nmodes); // Compute forward real FFT rfft(real_signal.begin(), size, output.begin()); @@ -63,7 +63,7 @@ TEST_P(FFTTest, rfft) { TEST_P(FFTTest, irfft) { - Array output{size}; + Array output(size); // Compute inverse real FFT irfft(fft_signal.begin(), size, output.begin()); From 8ba9503807e2f526912356404621281b97b259cb Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 14 Mar 2026 13:26:30 -0700 Subject: [PATCH 08/26] PetscOperators: mapping fixes and unit tests Some care needed with yup/down and yp/ym indices. --- include/bout/petsc_operators.hxx | 30 ++++++--- src/mesh/petsc_operators.cxx | 60 ++++++++++++----- tests/unit/mesh/test_petsc_operators.cxx | 85 +++++++++++++++++++++++- 3 files changed, 147 insertions(+), 28 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 3f4d0a596e..0b5c2851cb 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -10,6 +10,7 @@ /// #include "bout/array.hxx" +#include "bout/assert.hxx" #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/field3d.hxx" @@ -127,7 +128,7 @@ public: template void map(Function func) const { const std::vector>> regions = { - evolving_region, xin_region, xout_region, yup_region, ydown_region}; + evolving_region, xin_region, xout_region}; PetscInt row = row_start; for (const auto& region : regions) { BOUT_FOR_SERIAL(i, region.get()) { @@ -135,6 +136,15 @@ public: ++row; } } + BOUT_FOR_SERIAL(i, yup_region) { + func(row, ROUND(forward_cell_number[i])); + ++row; + } + BOUT_FOR_SERIAL(i, ydown_region) { + func(row, ROUND(backward_cell_number[i])); + ++row; + } + ASSERT0(row == row_end); } /// Iterate over locally stored entries packed from the main `Field3D`. @@ -169,7 +179,7 @@ public: void map_local_yup(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); BOUT_FOR_SERIAL(i, yup_region) { - func(row, i); + func(row, i.yp()); ++row; } } @@ -185,7 +195,7 @@ public: PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() + yup_region.size(); BOUT_FOR_SERIAL(i, ydown_region) { - func(row, i); + func(row, i.ym()); ++row; } } @@ -201,6 +211,8 @@ private: /// Mesh-global numbering for cells stored on this rank. Field3D cell_number; + Field3D forward_cell_number; + Field3D backward_cell_number; /// Evolving cells in the main domain. Region evolving_region; @@ -301,12 +313,9 @@ private: /// This constructor is used internally when creating operators from /// PETSc matrix algebra such as composition, addition, or subtraction. PetscOperator(PetscMappingPtr mapping, Mat mat) - : mapping(std::move(mapping)), mat_operator(mat) { - // Allocate working vectors - VecCreate(BoutComm::get(), &this->rhs_vec); - VecSetSizes(this->rhs_vec, this->mapping->size(), PETSC_DETERMINE); - VecDuplicate(this->rhs_vec, &this->result_vec); - } + : mapping(std::move(mapping)), mat_operator(mat), + rhs_vec(createVec(this->mapping->size())), + result_vec(createVec(this->mapping->size())) {} /// Shared mapping between mesh and PETSc numbering. PetscMappingPtr mapping; @@ -319,6 +328,9 @@ private: /// Work vector holding the packed result. Vec result_vec; + + // Allocate MPI vector + static Vec createVec(PetscInt local_size); }; /// Collection of PETSc operators read from a mesh. diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 8df2f4eeb7..2cc66b87f8 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -58,7 +58,9 @@ Region PetscMapping::create_region_xout(const Field3D& cell_number) { PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, const Field3D& backward_cell_number) - : cell_number(cell_number), evolving_region(create_region(cell_number)), + : cell_number(cell_number), forward_cell_number(forward_cell_number), + backward_cell_number(backward_cell_number), + evolving_region(create_region(cell_number)), xin_region(create_region_xin(cell_number)), xout_region(create_region_xout(cell_number)), yup_region(create_region(forward_cell_number)), @@ -109,22 +111,36 @@ PetscMapping::~PetscMapping() { PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights) - : mapping(std::move(mapping)) { + : mapping(std::move(mapping)), rhs_vec(createVec(this->mapping->size())), + result_vec(createVec(this->mapping->size())) { - Mat mat; + output.write("PetscOperator Array sizes {} {} {}\n", rows.size(), cols.size(), + weights.size()); + + Mat mat{nullptr}; BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); - int nlocal = mapping->size(); + int nlocal = this->mapping->size(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { + this->mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { if (mesh_index >= rows.size()) { return; // No weights -> skip } // Get the range of indices into columns and weights int start_ind = rows[mesh_index]; - int end_ind = (mesh_index + 1 < rows.size()) ? rows[mesh_index + 1] : cols.size(); - + if (start_ind < 0) { + return; // No entries + } + int end_ind = cols.size(); // End of the columns array + for (int i = mesh_index + 1; i < rows.size(); ++i) { + // rows[i] can be -1 if no weights + if (rows[i] > -1) { + // This is the next entry in the columns / weights array + end_ind = rows[i]; + break; + } + } BOUT_DO_PETSC(MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], &weights[start_ind], INSERT_VALUES)); }); @@ -133,13 +149,18 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, ArraygetPetscToMesh(), MAT_INITIAL_MATRIX, + BOUT_DO_PETSC(MatMatMult(mat, this->mapping->getPetscToMesh(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, &this->mat_operator)); + // Destroy temporary matrix + BOUT_DO_PETSC(MatDestroy(&mat)); +} - // Allocate vectors for operator input and output - BOUT_DO_PETSC(VecCreate(BoutComm::get(), &this->rhs_vec)); - BOUT_DO_PETSC(VecSetSizes(this->rhs_vec, this->mapping->size(), PETSC_DETERMINE)); - BOUT_DO_PETSC(VecDuplicate(this->rhs_vec, &this->result_vec)); +Vec PetscOperator::createVec(PetscInt local_size) { + Vec vec{nullptr}; + BOUT_DO_PETSC(VecCreate(BoutComm::get(), &vec)); + BOUT_DO_PETSC(VecSetType(vec, VECMPI)); + BOUT_DO_PETSC(VecSetSizes(vec, local_size, PETSC_DETERMINE)); + return vec; } PetscOperator::~PetscOperator() { @@ -150,18 +171,24 @@ PetscOperator::~PetscOperator() { /// Perform operation Field3D PetscOperator::operator()(const Field3D& rhs) const { + ASSERT1(rhs.hasParallelSlices()); + ASSERT1(rhs.yup().isAllocated()); + ASSERT1(rhs.ydown().isAllocated()); + // Fill vec from rhs PetscScalar* x; // Copy rows from Field3D cells BOUT_DO_PETSC(VecGetArray(this->rhs_vec, &x)); - mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = rhs[i]; }); + this->mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = rhs[i]; }); // Copy Yup / Ydown values from boundaries const Field3D& yup = rhs.yup(); - mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + this->mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); const Field3D& ydown = rhs.ydown(); - mapping->map_local_ydown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + this->mapping->map_local_ydown( + [&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + BOUT_DO_PETSC(VecRestoreArray(this->rhs_vec, &x)); // Perform Mat-Vec muliplication @@ -171,7 +198,8 @@ Field3D PetscOperator::operator()(const Field3D& rhs) const { Field3D result{emptyFrom(rhs)}; // This allocates memory const PetscScalar* r; BOUT_DO_PETSC(VecGetArrayRead(result_vec, &r)); - mapping->map_local_field([&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); + this->mapping->map_local_field( + [&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); BOUT_DO_PETSC(VecRestoreArrayRead(result_vec, &r)); return result; diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index b33d41a7b1..5528e8b5e1 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -1,5 +1,7 @@ #include "gtest/gtest.h" +#include "bout/array.hxx" +#include "bout/bout_types.hxx" #include "bout/output.hxx" #include "bout/output_bout_types.hxx" #include "bout/petsc_operators.hxx" @@ -8,6 +10,8 @@ #include "fake_mesh_fixture.hxx" #include "test_extras.hxx" +#include + // Reuse the "standard" fixture for FakeMesh using PetscMappingTest = FakeMeshFixture; @@ -54,9 +58,9 @@ TEST_F(PetscMappingTest, mapping) { Field3D cell_number{-1.0}; // No cells >= 0 // Note: 1 boundary cell in X and Y - cell_number(0, 0, 0) = 0; // Corner - cell_number(1, 1, 0) = 1; // In domain - cell_number(0, 1, 1) = 2; // Xin boundary + cell_number(0, 0, 0) = 2; // Corner + cell_number(1, 1, 0) = 0; // In domain + cell_number(0, 1, 1) = 1; // Xin boundary const Field3D forward_cell_number{-1.0}; const Field3D backward_cell_number{-1.0}; @@ -66,3 +70,78 @@ TEST_F(PetscMappingTest, mapping) { // Two cells: one evolving and one in xin ASSERT_EQ(2, mapping.size()); } + +using PetscOperatorTest = FakeMeshFixture; + +TEST_F(PetscOperatorTest, identity) { + Field3D cell_number{-1.0}; + const Field3D forward_cell_number{-1.0}; + const Field3D backward_cell_number{-1.0}; + + // Three cells active + cell_number(1, 1, 0) = 0; + cell_number(1, 2, 0) = 1; + cell_number(1, 2, 1) = 2; + + auto mapping = std::make_shared(cell_number, forward_cell_number, + backward_cell_number); + ASSERT_EQ(3, mapping->size()); + + const auto rows = Array::fromValues({0, 1, 2}); + const auto cols = Array::fromValues({0, 1, 2}); + const auto weights = Array::fromValues({1.0, 1.0, 1.0}); + + const PetscOperator identity(mapping, rows, cols, weights); + + Field3D value{0.0}; + value(1, 1, 0) = 10; + value(1, 2, 0) = 21; + value(1, 2, 1) = 32; + value.splitParallelSlices(); + value.yup() = -1.0; + value.ydown() = -1.0; + + const Field3D result = identity(value); + + EXPECT_EQ(10, result(1, 1, 0)); + EXPECT_EQ(21, result(1, 2, 0)); + EXPECT_EQ(32, result(1, 2, 1)); +} + +TEST_F(PetscOperatorTest, identity_yupdown) { + Field3D cell_number{-1.0}; + Field3D forward_cell_number{-1.0}; + Field3D backward_cell_number{-1.0}; + + // Three cells active + cell_number(1, 1, 0) = 0; + cell_number(1, 2, 0) = 1; + cell_number(1, 2, 1) = 2; + + forward_cell_number(1, 2, 0) = 3; + backward_cell_number(1, 1, 2) = 4; + + auto mapping = std::make_shared(cell_number, forward_cell_number, + backward_cell_number); + ASSERT_EQ(5, mapping->size()); + + const auto rows = Array::fromValues({0, 1, 2}); + const auto cols = Array::fromValues({0, 1, 2}); + const auto weights = Array::fromValues({1.0, 1.0, 1.0}); + + const PetscOperator identity(mapping, rows, cols, weights); + + Field3D value{0.0}; + value(1, 1, 0) = 10; + value(1, 2, 0) = 21; + value(1, 2, 1) = 32; + value.splitParallelSlices(); + value.yup() = -1.0; + value.ydown() = -1.0; + + const Field3D result = identity(value); + + EXPECT_EQ(10, result(1, 1, 0)); + EXPECT_EQ(21, result(1, 2, 0)); + EXPECT_EQ(32, result(1, 2, 1)); +} From 16e2768061dd48700381433d5902af810b6d4a11 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 15 Mar 2026 15:02:48 -0700 Subject: [PATCH 09/26] PetscOperator: Add diagonal and transpose methods `diagonal` creates diagonal matrices from Field3D. `transpose` creates a new matrix that is a transpose of a PetscOperator. This is useful for converting gradients into divergences. --- include/bout/petsc_operators.hxx | 15 +++++++ src/mesh/petsc_operators.cxx | 70 +++++++++++++++++++++++--------- 2 files changed, 65 insertions(+), 20 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 0b5c2851cb..046057bc5f 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -270,6 +270,9 @@ public: /// Destroy the PETSc matrix and working vectors owned by this operator. ~PetscOperator(); + /// Create a diagonal operator + static PetscOperator diagonal(PetscMappingPtr mapping, const Field3D& f); + /// Apply the operator to a field. /// /// @param rhs Input field, including any required `yup` and `ydown` @@ -304,6 +307,9 @@ public: /// Both operators must use compatible mappings. PetscOperator operator-(const PetscOperator& rhs) const; + /// Calculate transpose as a new matrix + PetscOperator transpose() const; + private: /// Construct directly from an existing PETSc matrix. /// @@ -329,6 +335,10 @@ private: /// Work vector holding the packed result. Vec result_vec; + /// Copy contents of f, f.yup() and f.ydown() into vec. + /// Assumes that vec has already been allocated. + static void copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec); + // Allocate MPI vector static Vec createVec(PetscInt local_size); }; @@ -361,6 +371,11 @@ public: /// Throws `BoutException` if any of the required arrays are missing. PetscOperator get(const std::string& name) const; + /// Create a diagonal operator + PetscOperator diagonal(const Field3D& f) const { + return PetscOperator::diagonal(mapping, f); + } + private: /// Read a `Field3D` from the mesh. /// diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 2cc66b87f8..89cf1a9503 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -120,7 +120,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymapping->size(); + const int nlocal = this->mapping->size(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); this->mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { @@ -128,7 +128,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array skip } // Get the range of indices into columns and weights - int start_ind = rows[mesh_index]; + const int start_ind = rows[mesh_index]; if (start_ind < 0) { return; // No entries } @@ -155,6 +155,26 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymap_local_field([&](PetscInt row, const Ind3D& i) { x[row] = f[i]; }); + + // Copy Yup / Ydown values from boundaries + const Field3D& yup = f.yup(); + mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + + const Field3D& ydown = f.ydown(); + mapping->map_local_ydown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + + BOUT_DO_PETSC(VecRestoreArray(vec, &x)); +} + Vec PetscOperator::createVec(PetscInt local_size) { Vec vec{nullptr}; BOUT_DO_PETSC(VecCreate(BoutComm::get(), &vec)); @@ -169,34 +189,38 @@ PetscOperator::~PetscOperator() { VecDestroy(&this->result_vec); } -/// Perform operation -Field3D PetscOperator::operator()(const Field3D& rhs) const { - ASSERT1(rhs.hasParallelSlices()); - ASSERT1(rhs.yup().isAllocated()); - ASSERT1(rhs.ydown().isAllocated()); +PetscOperator PetscOperator::diagonal(PetscMappingPtr mapping, const Field3D& f) { + Vec diag{createVec(mapping->size())}; + copyToVec(mapping, f, diag); - // Fill vec from rhs - PetscScalar* x; - // Copy rows from Field3D cells - BOUT_DO_PETSC(VecGetArray(this->rhs_vec, &x)); - this->mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = rhs[i]; }); + Mat mat{nullptr}; + // Note: MatMatMult with diagonal and mpiaij not supported + // -> Create MPIAIJ + BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); + BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); + const int nlocal = mapping->size(); + BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat, 1, nullptr, 0, nullptr)); + BOUT_DO_PETSC(MatDiagonalSet(mat, diag, INSERT_VALUES)); + BOUT_DO_PETSC(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); - // Copy Yup / Ydown values from boundaries - const Field3D& yup = rhs.yup(); - this->mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + BOUT_DO_PETSC(VecDestroy(&diag)); - const Field3D& ydown = rhs.ydown(); - this->mapping->map_local_ydown( - [&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + return PetscOperator(std::move(mapping), mat); +} - BOUT_DO_PETSC(VecRestoreArray(this->rhs_vec, &x)); +/// Perform operation +Field3D PetscOperator::operator()(const Field3D& rhs) const { + // Fill vec from rhs + copyToVec(this->mapping, rhs, this->rhs_vec); // Perform Mat-Vec muliplication BOUT_DO_PETSC(MatMult(this->mat_operator, rhs_vec, result_vec)); // Copy result_vec into a Field3D Field3D result{emptyFrom(rhs)}; // This allocates memory - const PetscScalar* r; + const PetscScalar* r{nullptr}; BOUT_DO_PETSC(VecGetArrayRead(result_vec, &r)); this->mapping->map_local_field( [&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); @@ -232,6 +256,12 @@ PetscOperator PetscOperator::operator-(const PetscOperator& rhs) const { return PetscOperator(this->mapping, mat); } +PetscOperator PetscOperator::transpose() const { + Mat mat{nullptr}; + BOUT_DO_PETSC(MatTranspose(this->mat_operator, MAT_INITIAL_MATRIX, &mat)); + return PetscOperator(this->mapping, mat); +} + Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const { Field3D result; if (mesh->get(result, name) != 0) { From 5182ac40159a532a5a47f7865af48f672e06d29f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 16 Mar 2026 13:15:14 -0700 Subject: [PATCH 10/26] PetscOperators: map matrix rows in boundary if present Allows input operators to contain entries in boundary rows. This is important so that the gradient and divergence operators depend on boundary cell values. --- include/bout/parallel_boundary_region.hxx | 2 ++ include/bout/petsc_operators.hxx | 2 ++ src/mesh/petsc_operators.cxx | 4 +++- 3 files changed, 7 insertions(+), 1 deletion(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 4d5278d00f..0833c8d3bc 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -108,6 +108,8 @@ public: void next() final { ++bndry_position; } bool isDone() final { return (bndry_position == end(bndry_points)); } + std::size_t size() { return this->bndry_points.size(); } + // getter Ind3D ind() const { return bndry_position->index; } BoutReal s_x() const { return bndry_position->intersection.s_x; } diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 046057bc5f..d3b0e89db9 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -310,6 +310,8 @@ public: /// Calculate transpose as a new matrix PetscOperator transpose() const; + void view() const { MatView(this->mat_operator, PETSC_VIEWER_STDOUT_WORLD); } + private: /// Construct directly from an existing PETSc matrix. /// diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 89cf1a9503..dab858bd20 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -94,6 +94,8 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce // Iterate through regions in this order this->map([&](PetscInt row, PetscInt col) { // `row` is the PETSc index; `col` is the Mesh index + ASSERT1(row >= 0); + ASSERT1(col >= 0); const PetscScalar ONE = 1.0; BOUT_DO_PETSC(MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES)); }); @@ -123,7 +125,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymapping->size(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - this->mapping->map_evolving([&](PetscInt row, Ind3D, PetscInt mesh_index) { + this->mapping->map([&](PetscInt row, PetscInt mesh_index) { if (mesh_index >= rows.size()) { return; // No weights -> skip } From 68428b3c952171aae7cd010413947ddfa5ba76d4 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 16 Mar 2026 14:52:44 -0700 Subject: [PATCH 11/26] PetscOperators: Add BOUT_HAS_PETSC guards Not available unless configured with PETSc. --- include/bout/petsc_operators.hxx | 10 ++++++++++ src/mesh/petsc_operators.cxx | 8 +++++++- tests/unit/mesh/test_petsc_operators.cxx | 6 ++++++ 3 files changed, 23 insertions(+), 1 deletion(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index d3b0e89db9..a745ff9e0b 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -9,6 +9,10 @@ /// a `Mesh`. /// +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + #include "bout/array.hxx" #include "bout/assert.hxx" #include "bout/bout_types.hxx" @@ -411,3 +415,9 @@ private: /// Shared mapping used by all operators created from this mesh. PetscMappingPtr mapping; }; + +#else // BOUT_HAS_PETSC + +#warning PETSc not available. No PetscOperators. + +#endif // BOUT_HAS_PETSC diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index dab858bd20..d2bf7b76c0 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -1,4 +1,7 @@ -#include "bout/petsc_operators.hxx" +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + #include "bout/array.hxx" #include "bout/assert.hxx" #include "bout/bout_types.hxx" @@ -6,6 +9,7 @@ #include "bout/field3d.hxx" #include "bout/output.hxx" #include "bout/output_bout_types.hxx" +#include "bout/petsc_operators.hxx" #include "bout/petsclib.hxx" #include "bout/region.hxx" #include @@ -293,3 +297,5 @@ PetscOperator PetscOperators::get(const std::string& name) const { this->meshGetArray(this->mesh, name + "_columns"), this->meshGetArray(this->mesh, name + "_weights")); } + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 5528e8b5e1..bc0a81cc76 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -1,3 +1,7 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + #include "gtest/gtest.h" #include "bout/array.hxx" @@ -145,3 +149,5 @@ TEST_F(PetscOperatorTest, identity_yupdown) { EXPECT_EQ(21, result(1, 2, 0)); EXPECT_EQ(32, result(1, 2, 1)); } + +#endif // BOUT_HAS_PETSC From 2c45584d26d86f020581deb9c40aa40cb97f8996 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 16 Mar 2026 16:46:53 -0700 Subject: [PATCH 12/26] PetscMapping: Add localSize and globalSize Check number of cells against global mapping size, not local size. --- include/bout/petsc_operators.hxx | 15 ++++++++++++--- src/mesh/petsc_operators.cxx | 18 +++++++++--------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index a745ff9e0b..e95ccb6fc2 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -96,11 +96,20 @@ public: /// /// The count includes evolving cells and all boundary cells represented /// in the mapping on this processor. - unsigned int size() const { + unsigned int localSize() const { return evolving_region.size() + xin_region.size() + xout_region.size() + yup_region.size() + ydown_region.size(); } + unsigned int globalSize() const { + // Get global size + PetscInt globalRows{0}; + PetscInt globalCols{0}; + BOUT_DO_PETSC(MatGetSize(this->mat_mesh_to_petsc, &globalRows, &globalCols)); + ASSERT1(globalRows == globalCols); + return globalRows; + } + /// Iterate over locally owned evolving cells in PETSc row order. /// /// @tparam Function Callable with signature @@ -326,8 +335,8 @@ private: /// PETSc matrix algebra such as composition, addition, or subtraction. PetscOperator(PetscMappingPtr mapping, Mat mat) : mapping(std::move(mapping)), mat_operator(mat), - rhs_vec(createVec(this->mapping->size())), - result_vec(createVec(this->mapping->size())) {} + rhs_vec(createVec(this->mapping->localSize())), + result_vec(createVec(this->mapping->localSize())) {} /// Shared mapping between mesh and PETSc numbering. PetscMappingPtr mapping; diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index d2bf7b76c0..c0cc5d0012 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -70,7 +70,7 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce yup_region(create_region(forward_cell_number)), ydown_region(create_region(backward_cell_number)) { // Calculate size of each region - const unsigned int nlocal = this->size(); + const unsigned int nlocal = this->localSize(); output.write("Local {} : evolving {} xin {} xout {} yup {} ydown {}", nlocal, evolving_region.size(), xin_region.size(), xout_region.size(), yup_region.size(), ydown_region.size()); @@ -117,8 +117,8 @@ PetscMapping::~PetscMapping() { PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights) - : mapping(std::move(mapping)), rhs_vec(createVec(this->mapping->size())), - result_vec(createVec(this->mapping->size())) { + : mapping(std::move(mapping)), rhs_vec(createVec(this->mapping->localSize())), + result_vec(createVec(this->mapping->localSize())) { output.write("PetscOperator Array sizes {} {} {}\n", rows.size(), cols.size(), weights.size()); @@ -126,7 +126,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymapping->size(); + const int nlocal = this->mapping->localSize(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); this->mapping->map([&](PetscInt row, PetscInt mesh_index) { @@ -196,7 +196,7 @@ PetscOperator::~PetscOperator() { } PetscOperator PetscOperator::diagonal(PetscMappingPtr mapping, const Field3D& f) { - Vec diag{createVec(mapping->size())}; + Vec diag{createVec(mapping->localSize())}; copyToVec(mapping, f, diag); Mat mat{nullptr}; @@ -204,7 +204,7 @@ PetscOperator PetscOperator::diagonal(PetscMappingPtr mapping, const Field3D& f) // -> Create MPIAIJ BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); - const int nlocal = mapping->size(); + const int nlocal = mapping->localSize(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat, 1, nullptr, 0, nullptr)); BOUT_DO_PETSC(MatDiagonalSet(mat, diag, INSERT_VALUES)); @@ -285,9 +285,9 @@ PetscOperators::PetscOperators(Mesh* mesh) int mesh_total_cells; if (mesh->get(mesh_total_cells, "total_cells") == 0) { // Check total number of cells - if (mesh_total_cells != mapping->size()) { - throw BoutException("Total cells in mesh {} doesn't match mapping size {}", - mesh_total_cells, mapping->size()); + if (mesh_total_cells != mapping->globalSize()) { + throw BoutException("Total cells in mesh {} doesn't match global mapping size {}", + mesh_total_cells, mapping->globalSize()); } } } From d8597124024f2f3a923f5fd35cbeadb817bb8b45 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 16 Mar 2026 21:20:23 -0700 Subject: [PATCH 13/26] PetscOperator: Handle standard CSR and variations Handle non-standard CSR inputs: - Negative indices for missing rows - Missing final entry in rows array --- src/mesh/petsc_operators.cxx | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index c0cc5d0012..5b64b15775 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -129,6 +129,8 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymapping->localSize(); BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); + // This reads CSR format but is defensive in handling negative indices and missing + // final 'row' entry. this->mapping->map([&](PetscInt row, PetscInt mesh_index) { if (mesh_index >= rows.size()) { return; // No weights -> skip @@ -136,17 +138,22 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array -1) { // This is the next entry in the columns / weights array end_ind = rows[i]; break; } } + if (end_ind == start_ind) { + // Empty row in CSR format + return; + } BOUT_DO_PETSC(MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], &weights[start_ind], INSERT_VALUES)); }); From 350eba58c76e705b84f9e24fd5a662b0a35812b1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 16 Mar 2026 21:41:58 -0700 Subject: [PATCH 14/26] PetscOperators: Fix unit tests Changes PetscMapping::size() to PetscMapping::localSize() --- tests/unit/mesh/test_petsc_operators.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index bc0a81cc76..219a667bac 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -72,7 +72,7 @@ TEST_F(PetscMappingTest, mapping) { const PetscMapping mapping(cell_number, forward_cell_number, backward_cell_number); // Two cells: one evolving and one in xin - ASSERT_EQ(2, mapping.size()); + ASSERT_EQ(2, mapping.localSize()); } using PetscOperatorTest = FakeMeshFixture; @@ -89,7 +89,7 @@ TEST_F(PetscOperatorTest, identity) { auto mapping = std::make_shared(cell_number, forward_cell_number, backward_cell_number); - ASSERT_EQ(3, mapping->size()); + ASSERT_EQ(3, mapping->localSize()); const auto rows = Array::fromValues({0, 1, 2}); const auto cols = Array::fromValues({0, 1, 2}); @@ -127,7 +127,7 @@ TEST_F(PetscOperatorTest, identity_yupdown) { auto mapping = std::make_shared(cell_number, forward_cell_number, backward_cell_number); - ASSERT_EQ(5, mapping->size()); + ASSERT_EQ(5, mapping->localSize()); const auto rows = Array::fromValues({0, 1, 2}); const auto cols = Array::fromValues({0, 1, 2}); From fb34b2a662cf15e9035f150914a9d307682a81b0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 17 Mar 2026 21:30:54 -0700 Subject: [PATCH 15/26] PetscOperator: Use bout::petsc::UniqueVec and UniqueMat Defined in `bout/petsc_interface.hxx`, `bout::petsc::UniqueVec and `UniqueMat` are `std::unique_ptr` wrappers around PETSc Vec and Mat pointers, with custom deleters. PetscOperator copy constructor and assignment deleted, move assignment set to default. This enables some optimizations to operate in-place on temporaries. --- include/bout/petsc_interface.hxx | 62 +++++++++++++-------- include/bout/petsc_operators.hxx | 46 ++++++++++++---- src/mesh/petsc_operators.cxx | 92 +++++++++++++++----------------- 3 files changed, 119 insertions(+), 81 deletions(-) diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index ab77743163..cfe27597ff 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -62,6 +62,36 @@ #include +namespace bout { +namespace petsc { +struct VecDeleter { + void operator()(Vec* vec_ptr) const { + VecDestroy(vec_ptr); + delete vec_ptr; + } +}; + +/// Unique pointer wrapper around PETSc Vec +/// Assumes that the Vec has been allocated on the heap. +/// +/// e.g. UniqueVec(new Vec{}); +using UniqueVec = std::unique_ptr; + +struct MatDeleter { + void operator()(Mat* mat_ptr) const { + MatDestroy(mat_ptr); + delete mat_ptr; + } +}; + +/// Unique pointer wrapper around PETSc Mat +/// Assumes that the Mat has been allocated on the heap. +/// +/// e.g. UniqueMat(new Mat{}); +using UniqueMat = std::unique_ptr; +} // namespace petsc +} // namespace bout + /*! * A class which wraps PETSc vector objects, allowing them to be * indexed using the BOUT++ scheme. Note that boundaries are only @@ -70,7 +100,7 @@ template class PetscVector; template -void swap(PetscVector& first, PetscVector& second); +void swap(PetscVector& first, PetscVector& second) noexcept; template inline MPI_Comm getComm([[maybe_unused]] const T& field) { @@ -88,13 +118,6 @@ public: static_assert(bout::utils::is_Field_v, "PetscVector only works with Fields"); using ind_type = typename T::ind_type; - struct VectorDeleter { - void operator()(Vec* vec) const { - VecDestroy(vec); - delete vec; - } - }; - PetscVector() : vector(new Vec{}) {} ~PetscVector() = default; @@ -160,7 +183,7 @@ public: return *this; } - friend void swap(PetscVector& first, PetscVector& second); + friend void swap(PetscVector& first, PetscVector& second) noexcept; BoutReal& operator()(const ind_type& index) { #if CHECKLEVEL >= 1 @@ -239,12 +262,12 @@ public: const Vec* get() const { return vector.get(); } private: - PetscLib lib{}; - std::unique_ptr vector = nullptr; + PetscLib lib; + bout::petsc::UniqueVec vector = nullptr; IndexerPtr indexConverter{}; CELL_LOC location = CELL_LOC::deflt; bool initialised = false; - Array vector_values{}; + Array vector_values; }; /*! @@ -255,7 +278,7 @@ private: template class PetscMatrix; template -void swap(PetscMatrix& first, PetscMatrix& second); +void swap(PetscMatrix& first, PetscMatrix& second) noexcept; template class PetscMatrix { @@ -263,13 +286,6 @@ public: static_assert(bout::utils::is_Field_v, "PetscMatrix only works with Fields"); using ind_type = typename T::ind_type; - struct MatrixDeleter { - void operator()(Mat* mat) const { - MatDestroy(mat); - delete mat; - } - }; - /// Default constructor does nothing PetscMatrix() : matrix(new Mat()) {} ~PetscMatrix() = default; @@ -327,7 +343,7 @@ public: rhs.initialised = false; return *this; } - friend void swap(PetscMatrix& first, PetscMatrix& second); + friend void swap(PetscMatrix& first, PetscMatrix& second) noexcept; /*! * A class which is used to assign to a particular element of a PETSc @@ -540,7 +556,7 @@ private: * vice versa. */ template -void swap(PetscVector& first, PetscVector& second) { +void swap(PetscVector& first, PetscVector& second) noexcept { std::swap(first.vector, second.vector); std::swap(first.indexConverter, second.indexConverter); std::swap(first.location, second.location); @@ -552,7 +568,7 @@ void swap(PetscVector& first, PetscVector& second) { * vice versa. */ template -void swap(PetscMatrix& first, PetscMatrix& second) { +void swap(PetscMatrix& first, PetscMatrix& second) noexcept { std::swap(first.matrix, second.matrix); std::swap(first.indexConverter, second.indexConverter); std::swap(first.pt, second.pt); diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index e95ccb6fc2..9c7e53fa01 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -19,6 +19,7 @@ #include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/mesh.hxx" +#include "bout/petsc_interface.hxx" #include "bout/petsclib.hxx" #include "bout/region.hxx" #include "bout/utils.hxx" @@ -266,6 +267,9 @@ using PetscMappingPtr = std::shared_ptr; /// - composed with other operators; and /// - added or subtracted. class PetscOperator { + using UniqueVec = bout::petsc::UniqueVec; // unique_ptr to Vec + using UniqueMat = bout::petsc::UniqueMat; // unique_ptr to Mat + public: /// Construct an operator from CSR data in mesh-global numbering. /// @@ -280,8 +284,16 @@ public: PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights); + // No copying + PetscOperator(const PetscOperator&) = delete; + PetscOperator& operator=(const PetscOperator&) = delete; + + // Allow Moving + PetscOperator(PetscOperator&&) noexcept = default; + PetscOperator& operator=(PetscOperator&&) noexcept = default; + /// Destroy the PETSc matrix and working vectors owned by this operator. - ~PetscOperator(); + ~PetscOperator() = default; /// Create a diagonal operator static PetscOperator diagonal(PetscMappingPtr mapping, const Field3D& f); @@ -323,7 +335,23 @@ public: /// Calculate transpose as a new matrix PetscOperator transpose() const; - void view() const { MatView(this->mat_operator, PETSC_VIEWER_STDOUT_WORLD); } + void view() const { MatView(*this->mat_operator, PETSC_VIEWER_STDOUT_WORLD); } + + /// Multiply by a scalar + /// This version allocates a new Mat and so is for lvalues. + PetscOperator operator*(BoutReal scalar) const& { + UniqueMat mat{new Mat{}}; + BOUT_DO_PETSC(MatDuplicate(*mat_operator, MAT_COPY_VALUES, mat.get())); + BOUT_DO_PETSC(MatScale(*mat, scalar)); + return PetscOperator(this->mapping, std::move(mat)); + } + + /// Multiply by a scalar. + /// This version is for rvalue temporaries and modifies in-place. + PetscOperator operator*(BoutReal scalar) && { + BOUT_DO_PETSC(MatScale(*this->mat_operator, scalar)); + return std::move(*this); + } private: /// Construct directly from an existing PETSc matrix. @@ -333,8 +361,8 @@ private: /// /// This constructor is used internally when creating operators from /// PETSc matrix algebra such as composition, addition, or subtraction. - PetscOperator(PetscMappingPtr mapping, Mat mat) - : mapping(std::move(mapping)), mat_operator(mat), + PetscOperator(PetscMappingPtr mapping, UniqueMat mat) + : mapping(std::move(mapping)), mat_operator(std::move(mat)), rhs_vec(createVec(this->mapping->localSize())), result_vec(createVec(this->mapping->localSize())) {} @@ -342,20 +370,20 @@ private: PetscMappingPtr mapping; /// PETSc matrix implementing the operator. - Mat mat_operator; + UniqueMat mat_operator; /// Work vector holding the packed input field. - Vec rhs_vec; + UniqueVec rhs_vec; /// Work vector holding the packed result. - Vec result_vec; + UniqueVec result_vec; /// Copy contents of f, f.yup() and f.ydown() into vec. /// Assumes that vec has already been allocated. static void copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec); // Allocate MPI vector - static Vec createVec(PetscInt local_size); + static UniqueVec createVec(PetscInt local_size); }; /// Collection of PETSc operators read from a mesh. @@ -399,7 +427,7 @@ private: /// @return The requested `Field3D`. /// /// Throws `BoutException` if the field is not present. - Field3D meshGetField3D(Mesh* mesh, const std::string& name) const; + static Field3D meshGetField3D(Mesh* mesh, const std::string& name); /// Read a one-dimensional `Array` from the mesh. /// diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 5b64b15775..ae682781e9 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -9,6 +9,7 @@ #include "bout/field3d.hxx" #include "bout/output.hxx" #include "bout/output_bout_types.hxx" +#include "bout/petsc_interface.hxx" #include "bout/petsc_operators.hxx" #include "bout/petsclib.hxx" #include "bout/region.hxx" @@ -117,7 +118,8 @@ PetscMapping::~PetscMapping() { PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, Array weights) - : mapping(std::move(mapping)), rhs_vec(createVec(this->mapping->localSize())), + : mapping(std::move(mapping)), mat_operator(new Mat{nullptr}), + rhs_vec(createVec(this->mapping->localSize())), result_vec(createVec(this->mapping->localSize())) { output.write("PetscOperator Array sizes {} {} {}\n", rows.size(), cols.size(), @@ -163,7 +165,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraymapping->getPetscToMesh(), MAT_INITIAL_MATRIX, - PETSC_DETERMINE, &this->mat_operator)); + PETSC_DETERMINE, this->mat_operator.get())); // Destroy temporary matrix BOUT_DO_PETSC(MatDestroy(&mat)); } @@ -188,56 +190,48 @@ void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec BOUT_DO_PETSC(VecRestoreArray(vec, &x)); } -Vec PetscOperator::createVec(PetscInt local_size) { - Vec vec{nullptr}; - BOUT_DO_PETSC(VecCreate(BoutComm::get(), &vec)); - BOUT_DO_PETSC(VecSetType(vec, VECMPI)); - BOUT_DO_PETSC(VecSetSizes(vec, local_size, PETSC_DETERMINE)); +PetscOperator::UniqueVec PetscOperator::createVec(PetscInt local_size) { + UniqueVec vec(new Vec{nullptr}); + BOUT_DO_PETSC(VecCreate(BoutComm::get(), vec.get())); + BOUT_DO_PETSC(VecSetType(*vec, VECMPI)); + BOUT_DO_PETSC(VecSetSizes(*vec, local_size, PETSC_DETERMINE)); return vec; } -PetscOperator::~PetscOperator() { - MatDestroy(&this->mat_operator); - VecDestroy(&this->rhs_vec); - VecDestroy(&this->result_vec); -} - PetscOperator PetscOperator::diagonal(PetscMappingPtr mapping, const Field3D& f) { - Vec diag{createVec(mapping->localSize())}; - copyToVec(mapping, f, diag); + const UniqueVec diag{createVec(mapping->localSize())}; + copyToVec(mapping, f, *diag); - Mat mat{nullptr}; + UniqueMat mat(new Mat{nullptr}); // Note: MatMatMult with diagonal and mpiaij not supported // -> Create MPIAIJ - BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); - BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); - const int nlocal = mapping->localSize(); - BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat, 1, nullptr, 0, nullptr)); - BOUT_DO_PETSC(MatDiagonalSet(mat, diag, INSERT_VALUES)); - BOUT_DO_PETSC(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); - BOUT_DO_PETSC(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); - - BOUT_DO_PETSC(VecDestroy(&diag)); - - return PetscOperator(std::move(mapping), mat); + BOUT_DO_PETSC(MatCreate(BoutComm::get(), mat.get())); + BOUT_DO_PETSC(MatSetType(*mat, MATMPIAIJ)); + const PetscInt nlocal = mapping->localSize(); + BOUT_DO_PETSC(MatSetSizes(*mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*mat, 1, nullptr, 0, nullptr)); + BOUT_DO_PETSC(MatDiagonalSet(*mat, *diag, INSERT_VALUES)); + BOUT_DO_PETSC(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); + + return PetscOperator(std::move(mapping), std::move(mat)); } /// Perform operation Field3D PetscOperator::operator()(const Field3D& rhs) const { // Fill vec from rhs - copyToVec(this->mapping, rhs, this->rhs_vec); + copyToVec(this->mapping, rhs, *this->rhs_vec); // Perform Mat-Vec muliplication - BOUT_DO_PETSC(MatMult(this->mat_operator, rhs_vec, result_vec)); + BOUT_DO_PETSC(MatMult(*this->mat_operator, *rhs_vec, *result_vec)); // Copy result_vec into a Field3D Field3D result{emptyFrom(rhs)}; // This allocates memory const PetscScalar* r{nullptr}; - BOUT_DO_PETSC(VecGetArrayRead(result_vec, &r)); + BOUT_DO_PETSC(VecGetArrayRead(*result_vec, &r)); this->mapping->map_local_field( [&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); - BOUT_DO_PETSC(VecRestoreArrayRead(result_vec, &r)); + BOUT_DO_PETSC(VecRestoreArrayRead(*result_vec, &r)); return result; } @@ -245,37 +239,37 @@ Field3D PetscOperator::operator()(const Field3D& rhs) const { /// Operator composition PetscOperator PetscOperator::operator*(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); - Mat mat; - BOUT_DO_PETSC(MatMatMult(this->mat_operator, rhs.mat_operator, MAT_INITIAL_MATRIX, - PETSC_DETERMINE, &mat)); - return PetscOperator(this->mapping, mat); + UniqueMat mat(new Mat{nullptr}); + BOUT_DO_PETSC(MatMatMult(*this->mat_operator, *rhs.mat_operator, MAT_INITIAL_MATRIX, + PETSC_DETERMINE, mat.get())); + return PetscOperator(this->mapping, std::move(mat)); } /// Operator addition PetscOperator PetscOperator::operator+(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); - Mat mat; - BOUT_DO_PETSC(MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat)); - BOUT_DO_PETSC(MatAXPY(mat, 1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); - return PetscOperator(this->mapping, mat); + UniqueMat mat(new Mat{nullptr}); + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, mat.get())); + BOUT_DO_PETSC(MatAXPY(*mat, 1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(this->mapping, std::move(mat)); } /// Operator subtraction PetscOperator PetscOperator::operator-(const PetscOperator& rhs) const { ASSERT0(this->mapping == rhs.mapping); - Mat mat; - BOUT_DO_PETSC(MatDuplicate(mat_operator, MAT_COPY_VALUES, &mat)); - BOUT_DO_PETSC(MatAXPY(mat, -1.0, rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); - return PetscOperator(this->mapping, mat); + UniqueMat mat(new Mat{nullptr}); + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, mat.get())); + BOUT_DO_PETSC(MatAXPY(*mat, -1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(this->mapping, std::move(mat)); } PetscOperator PetscOperator::transpose() const { - Mat mat{nullptr}; - BOUT_DO_PETSC(MatTranspose(this->mat_operator, MAT_INITIAL_MATRIX, &mat)); - return PetscOperator(this->mapping, mat); + UniqueMat mat(new Mat{nullptr}); + BOUT_DO_PETSC(MatTranspose(*this->mat_operator, MAT_INITIAL_MATRIX, mat.get())); + return PetscOperator(this->mapping, std::move(mat)); } -Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const { +Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { Field3D result; if (mesh->get(result, name) != 0) { throw BoutException("PetscOperators requires field '{}'", name); @@ -289,7 +283,7 @@ PetscOperators::PetscOperators(Mesh* mesh) meshGetField3D(mesh, "forward_cell_number"), meshGetField3D(mesh, "backward_cell_number"))) { - int mesh_total_cells; + int mesh_total_cells{0}; if (mesh->get(mesh_total_cells, "total_cells") == 0) { // Check total number of cells if (mesh_total_cells != mapping->globalSize()) { From 1c48e719bcd898e8f63d053d41395234662cf911 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 18 Mar 2026 09:00:34 -0700 Subject: [PATCH 16/26] examples/fci-operators Shows how to use the PetscOperators, and compares against the yup/ydown operators. Includes a script to generate the mesh, which currently requires a branch of Zoidberg to run. The code to calculate Grad_par, Div_par and Div_par_Grad_par from forward and backward operators is now in `PetscOperators::getParallel()`. That calculates and returns a collection of `Parallel` operators. This example could be turned into a test: - `Grad_par` should be close between the two methods - Volume integrals of `Div_par` and `Div_par_Grad_par` results should be near zero (flux conserved). --- examples/CMakeLists.txt | 1 + examples/fci-operators/CMakeLists.txt | 14 +++ examples/fci-operators/data/BOUT.inp | 5 + .../fci-operators/fci_operators_example.cxx | 48 ++++++++ examples/fci-operators/makeplots.py | 81 +++++++++++++ examples/fci-operators/rotating-ellipse.py | 91 +++++++++++++++ include/bout/petsc_operators.hxx | 15 +++ src/mesh/petsc_operators.cxx | 106 ++++++++++++++++++ 8 files changed, 361 insertions(+) create mode 100644 examples/fci-operators/CMakeLists.txt create mode 100644 examples/fci-operators/data/BOUT.inp create mode 100644 examples/fci-operators/fci_operators_example.cxx create mode 100644 examples/fci-operators/makeplots.py create mode 100755 examples/fci-operators/rotating-ellipse.py diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index f2444e6bf6..5eabaf045d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -18,6 +18,7 @@ add_subdirectory(dalf3) add_subdirectory(eigen-box) add_subdirectory(elm-pb) add_subdirectory(elm-pb-outerloop) +add_subdirectory(fci-operators) add_subdirectory(fci-wave) add_subdirectory(finite-volume/diffusion) add_subdirectory(finite-volume/fluid) diff --git a/examples/fci-operators/CMakeLists.txt b/examples/fci-operators/CMakeLists.txt new file mode 100644 index 0000000000..e010f45b50 --- /dev/null +++ b/examples/fci-operators/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 3.13) + +project(fci-operators LANGUAGES CXX) + +if(NOT TARGET bout++::bout++) + find_package(bout++ REQUIRED) +endif() + +bout_add_example( + fci-operators + SOURCES fci_operators_example.cxx + DATA_DIRS data + EXTRA_FILES rotating-ellipse-20x8x20.fci.nc makeplots.py +) diff --git a/examples/fci-operators/data/BOUT.inp b/examples/fci-operators/data/BOUT.inp new file mode 100644 index 0000000000..182d49cb96 --- /dev/null +++ b/examples/fci-operators/data/BOUT.inp @@ -0,0 +1,5 @@ +[mesh] +file = "rotating-ellipse-20x8x20.fci.nc" + +[mesh:paralleltransform] +type = fci diff --git a/examples/fci-operators/fci_operators_example.cxx b/examples/fci-operators/fci_operators_example.cxx new file mode 100644 index 0000000000..682d7c287c --- /dev/null +++ b/examples/fci-operators/fci_operators_example.cxx @@ -0,0 +1,48 @@ +#include "bout/bout.hxx" +#include "bout/difops.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/petsc_operators.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + const PetscOperators ops(bout::globals::mesh); + + // Parallel operators + const auto parallel = ops.getParallel(); + + // Test field + Field3D f = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f); + f.applyParallelBoundary("parallel_dirichlet_o2"); + + Field3D f_neumann = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f_neumann); + f_neumann.applyParallelBoundary("parallel_neumann_o1"); + + Options dump; + dump["f"] = f; + dump["dV"] = parallel.dV; + + dump["grad_par_op"] = parallel.Grad_par(f); + dump["grad_par_yud"] = Grad_par(f); + + dump["div_par_op"] = parallel.Div_par(f); + + auto* coords = bout::globals::mesh->getCoordinates(); + bout::globals::mesh->communicate(coords->Bxy); + coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); + dump["div_par_yud"] = Div_par(f); + + dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); + dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); + return 0; +} diff --git a/examples/fci-operators/makeplots.py b/examples/fci-operators/makeplots.py new file mode 100644 index 0000000000..b10a5829f3 --- /dev/null +++ b/examples/fci-operators/makeplots.py @@ -0,0 +1,81 @@ +from boutdata import collect +import matplotlib.pyplot as plt +import numpy as np + +path = "data" + +grad_par_op = collect("grad_par_op", path=path) +grad_par_yud = collect("grad_par_yud", path=path) + +div_par_op = collect("div_par_op", path=path) +div_par_yud = collect("div_par_yud", path=path) + +dV = collect("dV", path=path) # Cell volume + +# Integrate divergence +op_div_int = np.sum((dV * div_par_op)[2:-2, :, :]) +yud_div_int = np.sum((dV * div_par_yud)[2:-2, :, :]) + +print(f"Div_par volume integral Op: {op_div_int} Yup/down: {yud_div_int}") + +div_par_grad_par_op = collect("div_par_grad_par_op", path=path) +div_par_grad_par_yud = collect("div_par_grad_par_yud", path=path) + +# Integrate divergence +op_div2_int = np.sum((dV * div_par_grad_par_op)[2:-2, :, :]) +yud_div2_int = np.sum((dV * div_par_grad_par_yud)[2:-2, :, :]) + +print(f"Div_par_Grad_par volume integral Op: {op_div2_int} Yup/down: {yud_div2_int}") + + +def plot_comparison(a, alabel, b, blabel): + fig, axs = plt.subplots(1, 3, figsize=(10, 3)) + vmax = max([np.amax(a), np.amax(b)]) + vmin = min([np.amin(a), np.amin(b)]) + + cs = axs[0].contourf(a, vmin=vmin, vmax=vmax) + fig.colorbar(cs, ax=axs[0]) + axs[0].set_title(alabel) + + axs[1].contourf(b, vmin=vmin, vmax=vmax) + axs[1].set_title(blabel) + + cs = axs[2].contourf(a - b) + axs[2].set_title("Difference") + fig.colorbar(cs, ax=axs[2]) + + return fig + + +fig = plot_comparison( + grad_par_op[2:-2, 2, :], "Operator", grad_par_yud[2:-2, 2, :], "Yup/down" +) +fig.suptitle("Parallel Gradient") +fig.tight_layout() +fig.savefig("gradient.pdf") +fig.savefig("gradient.png") + +fig = plot_comparison( + div_par_op[2:-2, 2, :], + f"Operator [integral {op_div_int:.2e}]", + div_par_yud[2:-2, 2, :], + f"Yup/down [integral {yud_div_int:.2e}]", +) +fig.suptitle("Parallel Divergence") +fig.tight_layout() +fig.savefig("divergence.pdf") +fig.savefig("divergence.png") + +fig = plot_comparison( + div_par_grad_par_op[2:-2, 2, :], + f"Operator [integral {op_div2_int:.2e}]", + div_par_grad_par_yud[2:-2, 2, :], + f"Yup/down [integral {yud_div2_int:.2e}]", +) +fig.suptitle("Parallel diffusion") +fig.tight_layout() +fig.savefig("diffusion.pdf") +fig.savefig("diffusion.png") + + +plt.show() diff --git a/examples/fci-operators/rotating-ellipse.py b/examples/fci-operators/rotating-ellipse.py new file mode 100755 index 0000000000..79e95a9347 --- /dev/null +++ b/examples/fci-operators/rotating-ellipse.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python + +# Create grids for a rotating ellipse stellarator, based on curvilinear grids + +import numpy as np + +import zoidberg + +############################################################################# +# Define the magnetic field + +# Length in y after which the coils return to their starting (R,Z) locations +yperiod = 2.0 * np.pi / 5 + +xcentre = 2.0 + +magnetic_field = zoidberg.field.RotatingEllipse( + xcentre=xcentre, # Major radius of axis [m] + radius=0.8, # Minor radius of the coils [m] + I_coil=0.4, # Coil current + yperiod=yperiod, + Btor=1.0, # Toroidal magnetic field +) + +############################################################################# +# Create the inner flux surface, starting at a point at phi=0 +# To do this we need to define the y locations of the poloidal points +# where we will construct grids + +start_r = xcentre + 0.4 +start_z = 0.0 + +nslices = 8 # Number of poloidal slices +ycoords = np.linspace(0, yperiod, nslices) +npoints = 20 # Points per poloidal slice + +rzcoord, ycoords = zoidberg.fieldtracer.trace_poincare( + magnetic_field, start_r, start_z, yperiod, y_slices=ycoords, revs=npoints, nover=1 +) + +inner_lines = [] +for i in range(nslices): + r = rzcoord[:, i, 0, 0] + z = rzcoord[:, i, 0, 1] + line = zoidberg.rzline.line_from_points(r, z) + # Re-map the points so they're approximately uniform in distance along the surface + # Note that this results in some motion of the line + line = line.equallySpaced() + inner_lines.append(line) + +# Now have a list of y coordinates (ycoords) and inner lines (inner_lines) + +############################################################################# +# Generate a fixed circle for the outer boundary + +outer_line = zoidberg.rzline.circle(R0=xcentre, r=0.6) + +############################################################################# +# Now have inner and outer boundaries for each poloidal grid +# Generate a grid on each poloidal slice using the elliptic grid generator + +nx = 20 +nz = 20 + +pol_grids = [ + zoidberg.poloidal_grid.grid_elliptic(inner_line, outer_line, nx, nz) + for inner_line in inner_lines +] + +############################################################################# +# Create a grid, then calculate forward and backward maps + +grid = zoidberg.grid.Grid(pol_grids, ycoords, yperiod, yperiodic=True) + +maps = zoidberg.make_maps(grid, magnetic_field) + +############################################################################# +# Interpolation weights + +weight_vars = zoidberg.weights.calculate_parallel_map_operators(maps) +maps.update(weight_vars) + +############################################################################# +# Write grid file + +filename = f"rotating-ellipse-{nx}x{nslices}x{nz}.fci.nc" + +print(f"Writing to grid file '{filename}'") +zoidberg.write_maps( + grid, magnetic_field, maps, gridfile=filename, new_names=False, metric2d=False +) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 9c7e53fa01..cf87048420 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -419,6 +419,21 @@ public: return PetscOperator::diagonal(mapping, f); } + /// A standard collection of operators + /// parallel to the magnetic field + struct Parallel { + PetscOperator Grad_par; ///< Gradient + PetscOperator Div_par; ///< Divergence + PetscOperator Div_par_Grad_par; ///< Diffusion + Field3D dV; ///< Cell volume + }; + + /// Calculate a set of parallel operators + /// + /// This assumes that "forward" and "backward" operators + /// are stored in the mesh. + Parallel getParallel() const; + private: /// Read a `Field3D` from the mesh. /// diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index ae682781e9..88dfa94750 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -299,4 +299,110 @@ PetscOperator PetscOperators::get(const std::string& name) const { this->meshGetArray(this->mesh, name + "_weights")); } +PetscOperators::Parallel PetscOperators::getParallel() const { + // Read maps from the mesh + auto forward = this->get("forward"); + auto backward = this->get("backward"); + + // ---- Construct Grad_par ---- + // + // Create a parallel gradient operator by combining the parallel + // length dl = dy * sqrt(g_22) with forward & backward operators + auto* coords = this->mesh->getCoordinates(); + Field3D dl = coords->dy * sqrt(coords->g_22); + dl.splitParallelSlices(); + dl.yup() = 0.0; + dl.ydown() = 0.0; + dl.applyParallelBoundary("parallel_neumann_o1"); + + auto inv_2dl = 0.5 / dl; + inv_2dl.splitParallelSlices(); + inv_2dl.yup() = 0.0; + inv_2dl.ydown() = 0.0; + inv_2dl.applyParallelBoundary("parallel_neumann_o1"); + + inv_2dl.yup() *= 0.5; + inv_2dl.ydown() *= 0.5; + + auto inv_2dl_op = this->diagonal(inv_2dl); + + auto Grad_par = inv_2dl_op * (forward - backward); + + // ---- Construct Div_par ---- + // + // Use the Support Operator Method (SOM) to calculate + // Div_par from Grad_par and cell volumes. + + // Cell volume + auto dV = coords->J * coords->dx * coords->dy * coords->dz; + dV.splitParallelSlices(); + dV.yup() = 0.0; + dV.ydown() = 0.0; + dV.applyParallelBoundary("parallel_neumann_o1"); + auto dV_op = this->diagonal(dV); + + Field3D neg_inv_dV = -1. / dV; + neg_inv_dV.splitParallelSlices(); + neg_inv_dV.yup() = 0.0; + neg_inv_dV.ydown() = 0.0; + neg_inv_dV.applyParallelBoundary("parallel_neumann_o1"); + auto neg_inv_dV_op = this->diagonal(neg_inv_dV); + + auto Div_par = neg_inv_dV_op * Grad_par.transpose() * dV_op; + + // ---- Construct Div_par_Grad_par ---- + // + // Requires gradients between planes, at +1/2 and -1/2, and interpolation + // operators to calculate quantities between cells. + + // Identity operator + Field3D one{1.0}; + one.splitParallelSlices(); + one.yup() = 1.0; + one.ydown() = 1.0; + const auto identity = this->diagonal(one); + + // Interpolate at + 1/2 + const auto interp_plus_op = (identity + forward) * 0.5; + + // dl averaged at +1/2 + const Field3D dl_plus = interp_plus_op(dl); + Field3D inv_dl_plus = 1. / dl_plus; + inv_dl_plus.splitParallelSlices(); + inv_dl_plus.yup() = 0.0; + inv_dl_plus.ydown() = 0.0; + inv_dl_plus.applyParallelBoundary("parallel_neumann_o1"); + const auto inv_dl_plus_op = this->diagonal(inv_dl_plus); + + // Gradient at + 1/2 + const auto Grad_plus = inv_dl_plus_op * (forward - identity); + + // Divergence at -1/2 + const auto Div_minus = neg_inv_dV_op * Grad_plus.transpose() * dV_op; + + // Interpolate at - 1/2 + auto interp_minus_op = (identity + backward) * 0.5; + + // dl averaged at -1/2 + const Field3D dl_minus = interp_minus_op(dl); + Field3D inv_dl_minus = 1. / dl_minus; + inv_dl_minus.splitParallelSlices(); + inv_dl_minus.yup() = 0.0; + inv_dl_minus.ydown() = 0.0; + inv_dl_minus.applyParallelBoundary("parallel_neumann_o1"); + const auto inv_dl_minus_op = this->diagonal(inv_dl_minus); + + // Gradient at - 1/2 + auto Grad_minus = inv_dl_minus_op * (identity - backward); + + // Divergence at +1/2 + const auto Div_plus = neg_inv_dV_op * Grad_minus.transpose() * dV_op; + + // Div(Grad_par()) operator + auto Div_par_Grad_par = ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; + + return Parallel{std::move(Grad_par), std::move(Div_par), std::move(Div_par_Grad_par), + std::move(dV)}; +} + #endif // BOUT_HAS_PETSC From d5e1218b8c6d7bec3e9b3877c23ca129df9f2006 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 18 Mar 2026 11:13:56 -0700 Subject: [PATCH 17/26] PetscOperators::Parallel::Div_par_K_Grad_par This is a function rather than a PetscOperator because it requires at least two matrix-vector products. This implementation requires six Mat-Vec products that could probably be optimized. --- include/bout/petsc_operators.hxx | 20 ++++++++++++++++---- src/mesh/petsc_operators.cxx | 32 +++++++++++++++++++++++++------- 2 files changed, 41 insertions(+), 11 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index cf87048420..8ab95d0f17 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -244,13 +244,13 @@ private: Region ydown_region; /// First and one-past-last PETSc global rows owned by this rank. - PetscInt row_start, row_end; + PetscInt row_start{}, row_end{}; /// Permutation matrix mapping mesh-global ordering to PETSc ordering. - Mat mat_mesh_to_petsc; + Mat mat_mesh_to_petsc{nullptr}; /// Permutation matrix mapping PETSc ordering to mesh-global ordering. - Mat mat_petsc_to_mesh; + Mat mat_petsc_to_mesh{nullptr}; }; /// Shared pointer to an immutable `PetscMapping`. @@ -380,7 +380,7 @@ private: /// Copy contents of f, f.yup() and f.ydown() into vec. /// Assumes that vec has already been allocated. - static void copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec); + static void copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec); // Allocate MPI vector static UniqueVec createVec(PetscInt local_size); @@ -426,6 +426,18 @@ public: PetscOperator Div_par; ///< Divergence PetscOperator Div_par_Grad_par; ///< Diffusion Field3D dV; ///< Cell volume + + PetscOperator Grad_minus; + PetscOperator Grad_plus; + PetscOperator Div_minus; + PetscOperator Div_plus; + + PetscOperator Interp_minus; + PetscOperator Interp_plus; + + /// This is a bilinear operator that requires multiple + /// matrix-vector products. + Field3D Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const; }; /// Calculate a set of parallel operators diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 88dfa94750..98db8eea08 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -170,7 +170,7 @@ PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Arraydiagonal(one); // Interpolate at + 1/2 - const auto interp_plus_op = (identity + forward) * 0.5; + auto interp_plus_op = (identity + forward) * 0.5; // dl averaged at +1/2 const Field3D dl_plus = interp_plus_op(dl); @@ -375,10 +375,10 @@ PetscOperators::Parallel PetscOperators::getParallel() const { const auto inv_dl_plus_op = this->diagonal(inv_dl_plus); // Gradient at + 1/2 - const auto Grad_plus = inv_dl_plus_op * (forward - identity); + auto Grad_plus = inv_dl_plus_op * (forward - identity); // Divergence at -1/2 - const auto Div_minus = neg_inv_dV_op * Grad_plus.transpose() * dV_op; + auto Div_minus = neg_inv_dV_op * Grad_plus.transpose() * dV_op; // Interpolate at - 1/2 auto interp_minus_op = (identity + backward) * 0.5; @@ -396,13 +396,31 @@ PetscOperators::Parallel PetscOperators::getParallel() const { auto Grad_minus = inv_dl_minus_op * (identity - backward); // Divergence at +1/2 - const auto Div_plus = neg_inv_dV_op * Grad_minus.transpose() * dV_op; + auto Div_plus = neg_inv_dV_op * Grad_minus.transpose() * dV_op; // Div(Grad_par()) operator auto Div_par_Grad_par = ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; - return Parallel{std::move(Grad_par), std::move(Div_par), std::move(Div_par_Grad_par), - std::move(dV)}; + return Parallel{std::move(Grad_par), std::move(Div_par), + std::move(Div_par_Grad_par), std::move(dV), + std::move(Grad_minus), std::move(Grad_plus), + std::move(Div_minus), std::move(Div_plus), + std::move(interp_minus_op), std::move(interp_plus_op)}; +} + +Field3D PetscOperators::Parallel::Div_par_K_Grad_par(const Field3D& K, + const Field3D& f) const { + // Calculate gradients in + and - directions + const Field3D grad_plus = this->Grad_plus(f); + const Field3D grad_minus = this->Grad_minus(f); + + // Interpolate K to + and - locations + const Field3D K_plus = this->Interp_plus(K); + const Field3D K_minus = this->Interp_minus(K); + + // Calculate divergence + return (this->Div_minus(K_plus * grad_plus) + this->Div_plus(K_minus * grad_minus)) + * 0.5; } #endif // BOUT_HAS_PETSC From 9a44ff8a565a75fe6dcd7c3f1bbdd279b254e472 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 19 Mar 2026 20:37:51 -0700 Subject: [PATCH 18/26] WIP: Sub-cell sampling Trace multiple field lines per cell, to construct smoother maps. A complication is cells that partially intersect a boundary. --- examples/fci-operators/rotating-ellipse.py | 15 +++++++++-- src/mesh/petsc_operators.cxx | 29 +++++++++++++++++++--- 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/examples/fci-operators/rotating-ellipse.py b/examples/fci-operators/rotating-ellipse.py index 79e95a9347..3aa485c093 100755 --- a/examples/fci-operators/rotating-ellipse.py +++ b/examples/fci-operators/rotating-ellipse.py @@ -72,18 +72,29 @@ grid = zoidberg.grid.Grid(pol_grids, ycoords, yperiod, yperiodic=True) -maps = zoidberg.make_maps(grid, magnetic_field) +samples_per_dim = 5 # Sub-sampling in each cell + +maps = zoidberg.make_maps(grid, magnetic_field, samples_per_dim=samples_per_dim) ############################################################################# # Interpolation weights weight_vars = zoidberg.weights.calculate_parallel_map_operators(maps) maps.update(weight_vars) +# Remove sub-sampled maps from output +for var in [ + "forward_xt_primes", + "forward_zt_primes", + "backward_xt_primes", + "backward_zt_primes", + "subcell_weights", +]: + maps.pop(var, None) ############################################################################# # Write grid file -filename = f"rotating-ellipse-{nx}x{nslices}x{nz}.fci.nc" +filename = f"rotating-ellipse-{nx}x{nslices}x{nz}-s{samples_per_dim}.fci.nc" print(f"Writing to grid file '{filename}'") zoidberg.write_maps( diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 98db8eea08..2b11c3aa7c 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -321,9 +321,6 @@ PetscOperators::Parallel PetscOperators::getParallel() const { inv_2dl.ydown() = 0.0; inv_2dl.applyParallelBoundary("parallel_neumann_o1"); - inv_2dl.yup() *= 0.5; - inv_2dl.ydown() *= 0.5; - auto inv_2dl_op = this->diagonal(inv_2dl); auto Grad_par = inv_2dl_op * (forward - backward); @@ -339,7 +336,19 @@ PetscOperators::Parallel PetscOperators::getParallel() const { dV.yup() = 0.0; dV.ydown() = 0.0; dV.applyParallelBoundary("parallel_neumann_o1"); - auto dV_op = this->diagonal(dV); + + // Fractional boundary cells + const auto forward_boundary_fraction = + meshGetField3D(this->mesh, "forward_boundary_fraction"); + const auto backward_boundary_fraction = + meshGetField3D(this->mesh, "backward_boundary_fraction"); + + BOUT_FOR(i, dV.getRegion("RGN_NOBNDRY")) { + dV.yup()[i.yp()] *= forward_boundary_fraction[i]; + dV.ydown()[i.ym()] *= backward_boundary_fraction[i]; + } + + const auto dV_op = this->diagonal(dV); Field3D neg_inv_dV = -1. / dV; neg_inv_dV.splitParallelSlices(); @@ -410,6 +419,18 @@ PetscOperators::Parallel PetscOperators::getParallel() const { Field3D PetscOperators::Parallel::Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const { + + // There are 4 matrix-vector products that could be performed in parallel. + // The best way to optimize this is probably to pack f and K into one Vec, + // and assemble a single sparse matrix that contains the four blocks: + // + // (Grad_plus 0 ) (grad_plus ) + // (Grad_minus 0 ) (f) -> (grad_minus ) + // ( 0 Interp_plus ) (K) (K_plus ) + // ( 0 Interp_minus) (K_minus ) + // + // For now we perform each operation in sequence. + // Calculate gradients in + and - directions const Field3D grad_plus = this->Grad_plus(f); const Field3D grad_minus = this->Grad_minus(f); From 3992b8c7be99c66e6733600e09f8e3261f03d708 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 19 Mar 2026 20:46:53 -0700 Subject: [PATCH 19/26] fci-operators example: Don't include mesh file Mesh not available in repo, only generation script. --- examples/fci-operators/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fci-operators/CMakeLists.txt b/examples/fci-operators/CMakeLists.txt index e010f45b50..3a19ea9dcc 100644 --- a/examples/fci-operators/CMakeLists.txt +++ b/examples/fci-operators/CMakeLists.txt @@ -10,5 +10,5 @@ bout_add_example( fci-operators SOURCES fci_operators_example.cxx DATA_DIRS data - EXTRA_FILES rotating-ellipse-20x8x20.fci.nc makeplots.py + EXTRA_FILES rotating-ellipse.py makeplots.py ) From d2f26bb31dac5d8165b5978ce94d06d83f4f57c7 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 19 Mar 2026 21:05:51 -0700 Subject: [PATCH 20/26] PetscOperators: Specify dV is Field3D Fix for when metrics are 2D. --- src/mesh/petsc_operators.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 2b11c3aa7c..6767ebe716 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -331,7 +331,7 @@ PetscOperators::Parallel PetscOperators::getParallel() const { // Div_par from Grad_par and cell volumes. // Cell volume - auto dV = coords->J * coords->dx * coords->dy * coords->dz; + Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; dV.splitParallelSlices(); dV.yup() = 0.0; dV.ydown() = 0.0; From ed56609095b1aba1b8ac50d168f1ba566b3e8453 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 23 Mar 2026 15:18:06 -0700 Subject: [PATCH 21/26] PetscOperators: Change to use C, L+ and L- spaces More consistent treatment of where values are defined: On cells (C), forward (L+) and backward (L-) legs. Uses typed vectors and operators to ensure that quantities are used consistently. Coauthor: ChatGPT. --- include/bout/petsc_operators.hxx | 744 ++++++++++++++++--------------- src/mesh/petsc_operators.cxx | 605 +++++++++++-------------- 2 files changed, 643 insertions(+), 706 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 8ab95d0f17..b9d35a9793 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -1,14 +1,3 @@ -/// Represent finite-difference operators using PETSc matrices. -/// -/// This header defines: -/// - `PetscMapping`, which maps between mesh-global numbering, -/// PETSc numbering, and `Field3D` storage; -/// - `PetscOperator`, which wraps a PETSc `Mat` as a linear operator on -/// `Field3D` objects; and -/// - `PetscOperators`, a helper for reading mappings and operators from -/// a `Mesh`. -/// - #include "bout/build_defines.hxx" #if BOUT_HAS_PETSC @@ -25,156 +14,103 @@ #include "bout/utils.hxx" #include +#include #include #include #include +#include +#include #include #include -/// Map between mesh-global cell numbering and PETSc numbering. -/// -/// The mesh file stores operator stencils and weights using a mesh-global -/// numbering scheme. This numbering includes: -/// - evolving cells in the main domain; -/// - X-boundary cells; and -/// - Y-boundary cells represented by the `yup` and `ydown` fields. -/// -/// PETSc uses a different global numbering based on the MPI distribution -/// of rows and vector entries. `PetscMapping` records the correspondence -/// between these numbering systems and provides helper methods to iterate -/// over the local cells in a consistent order. -/// -/// The mapping order used by this class is: -/// 1. evolving cells; -/// 2. inner X-boundary cells; -/// 3. outer X-boundary cells; -/// 4. `yup` boundary cells; and -/// 5. `ydown` boundary cells. +/// Tag type for the cell space C. +struct CellSpaceTag {}; + +/// Tag type for the forward leg space L+. +struct ForwardLegSpaceTag {}; + +/// Tag type for the backward leg space L-. +struct BackwardLegSpaceTag {}; + +/// Mapping between stored indices in a mesh file and PETSc row ownership. /// -/// A pair of PETSc matrices is constructed internally to reorder vectors -/// between mesh-global and PETSc-global numbering. -class PetscMapping { +/// Every space is defined by a global stored numbering written to the mesh file. +/// PETSc distributes rows contiguously across MPI ranks, so this class stores the +/// local owned stored indices in row order and a permutation matrix that maps +/// from stored numbering to PETSc numbering. +class PetscIndexMapping { public: - /// Construct the mapping from mesh numbering fields. - /// - /// @param cell_number Mesh-global cell numbers for evolving cells. - /// @param forward_cell_number Mesh-global cell numbers associated with - /// forward Y-boundary (`yup`) values. - /// @param backward_cell_number Mesh-global cell numbers associated with - /// backward Y-boundary (`ydown`) values. - /// - /// Non-negative values in the numbering fields identify valid cells. - /// The constructor builds the local regions and the PETSc permutation - /// matrices needed to translate between mesh and PETSc numbering. - PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, - const Field3D& backward_cell_number); - - /// Destroy the PETSc matrices owned by this mapping. - ~PetscMapping(); - - /// Build a region containing all valid cells in a numbering field. - /// - /// @param cell_number Field whose entries contain mesh-global cell numbers. - /// @return A region containing all indices with non-negative cell numbers. - /// - /// This helper is public for unit testing. - static Region create_region(const Field3D& cell_number); + PetscIndexMapping() = default; + virtual ~PetscIndexMapping(); - /// Build a region containing valid cells on the inner X boundary. - /// - /// @param cell_number Field whose entries contain mesh-global cell numbers. - /// @return A region containing valid cells adjacent to the inner X boundary. - static Region create_region_xin(const Field3D& cell_number); + PetscIndexMapping(const PetscIndexMapping&) = delete; + PetscIndexMapping& operator=(const PetscIndexMapping&) = delete; + PetscIndexMapping(PetscIndexMapping&&) = delete; + PetscIndexMapping& operator=(PetscIndexMapping&&) = delete; - /// Build a region containing valid cells on the outer X boundary. - /// - /// @param cell_number Field whose entries contain mesh-global cell numbers. - /// @return A region containing valid cells adjacent to the outer X boundary. - static Region create_region_xout(const Field3D& cell_number); - - /// Return the total number of locally mapped cells. - /// - /// The count includes evolving cells and all boundary cells represented - /// in the mapping on this processor. - unsigned int localSize() const { - return evolving_region.size() + xin_region.size() + xout_region.size() - + yup_region.size() + ydown_region.size(); + /// Number of rows owned locally by this MPI rank. + PetscInt localSize() const { + return static_cast(local_stored_indices.size()); } - unsigned int globalSize() const { - // Get global size - PetscInt globalRows{0}; - PetscInt globalCols{0}; - BOUT_DO_PETSC(MatGetSize(this->mat_mesh_to_petsc, &globalRows, &globalCols)); - ASSERT1(globalRows == globalCols); - return globalRows; - } + /// Global number of rows in this space. + PetscInt globalSize() const { return global_size; } - /// Iterate over locally owned evolving cells in PETSc row order. - /// - /// @tparam Function Callable with signature - /// `func(PetscInt petsc_row, Ind3D index, int mesh_global_index)`. - /// @param func Function called once for each evolving cell owned locally. - /// - /// The PETSc row index starts at the locally owned global row `row_start`. - /// The supplied `Ind3D` index refers to the main `Field3D` storage. - template - void map_evolving(Function func) const { - PetscInt row = this->row_start; - BOUT_FOR_SERIAL(i, this->evolving_region) { - func(row, i, ROUND(cell_number[i])); - ++row; - } - } + /// First global PETSc row owned by this MPI rank. + PetscInt rowStart() const { return row_start; } + + /// One-past-last global PETSc row owned by this MPI rank. + PetscInt rowEnd() const { return row_end; } + + /// Stored indices owned locally in PETSc row order. + const std::vector& localStoredIndices() const { return local_stored_indices; } + + /// Map a locally owned stored index to a PETSc global row. + PetscInt storedToPetsc(int stored_index) const; + + /// Permutation matrix mapping PETSc ordering to stored numbering. + Mat getPetscToStored() const { return mat_petsc_to_stored; } + +protected: + PetscLib lib; + + void buildPermutation(PetscInt nlocal, PetscInt nglobal, + const std::vector& stored_indices); + + PetscInt global_size{0}; + PetscInt row_start{0}; + PetscInt row_end{0}; + std::vector local_stored_indices; + std::unordered_map local_stored_to_petsc; + Mat mat_stored_to_petsc{nullptr}; + Mat mat_petsc_to_stored{nullptr}; +}; + +/// Cell-space mapping between stored cell numbers, PETSc numbering, and Field3D. +class PetscCellMapping : public PetscIndexMapping { +public: + PetscCellMapping(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells); + + static Region create_region(const Field3D& cell_number); + static Region create_region_xin(const Field3D& cell_number); + static Region create_region_xout(const Field3D& cell_number); - /// Iterate over all locally mapped cells in PETSc row order. - /// - /// @tparam Function Callable with signature - /// `func(PetscInt petsc_row, int mesh_global_index)`. - /// @param func Function called once for each mapped cell owned locally. - /// - /// The iteration order is evolving, inner X boundary, outer X boundary, - /// `yup`, then `ydown`. - /// - /// No `Ind3D` is passed because some entries refer to interior `Field3D` - /// storage while others refer to boundary storage in `yup` or `ydown`. template - void map(Function func) const { - const std::vector>> regions = { - evolving_region, xin_region, xout_region}; - PetscInt row = row_start; - for (const auto& region : regions) { - BOUT_FOR_SERIAL(i, region.get()) { - func(row, ROUND(cell_number[i])); - ++row; - } - } - BOUT_FOR_SERIAL(i, yup_region) { - func(row, ROUND(forward_cell_number[i])); - ++row; - } - BOUT_FOR_SERIAL(i, ydown_region) { - func(row, ROUND(backward_cell_number[i])); + void mapOwnedInteriorCells(Function func) const { + PetscInt row = rowStart(); + BOUT_FOR_SERIAL(i, evolving_region) { + func(row, i, ROUND(cell_number[i])); ++row; } - ASSERT0(row == row_end); } - /// Iterate over locally stored entries packed from the main `Field3D`. - /// - /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. - /// @param func Function called once for each evolving or X-boundary cell - /// stored in the local PETSc vector layout. - /// - /// The index passed to `func` is a local vector index beginning at zero. - /// Only the main field storage is traversed here; Y-boundary values are - /// handled separately by `map_local_yup()` and `map_local_ydown()`. template - void map_local_field(Function func) const { + void mapLocalField(Function func) const { + PetscInt row = 0; const std::vector>> regions = { evolving_region, xin_region, xout_region}; - PetscInt row = 0; // Starting from 0, not row_start for (const auto& region : regions) { BOUT_FOR_SERIAL(i, region.get()) { func(row, i); @@ -183,14 +119,8 @@ public: } } - /// Iterate over local vector entries corresponding to `yup` boundary values. - /// - /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. - /// @param func Function called once for each locally stored `yup` entry. - /// - /// The local index is the offset into the packed local PETSc vector. template - void map_local_yup(Function func) const { + void mapLocalYup(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); BOUT_FOR_SERIAL(i, yup_region) { func(row, i.yp()); @@ -198,14 +128,8 @@ public: } } - /// Iterate over local vector entries corresponding to `ydown` boundary values. - /// - /// @tparam Function Callable with signature `func(PetscInt local_index, Ind3D index)`. - /// @param func Function called once for each locally stored `ydown` entry. - /// - /// The local index is the offset into the packed local PETSc vector. template - void map_local_ydown(Function func) const { + void mapLocalYdown(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() + yup_region.size(); BOUT_FOR_SERIAL(i, ydown_region) { @@ -214,274 +138,372 @@ public: } } - /// Return the PETSc matrix that maps PETSc ordering to mesh ordering. - /// - /// @return PETSc matrix representing the permutation from PETSc global - /// indices to mesh-global indices. - Mat getPetscToMesh() const { return mat_petsc_to_mesh; } - private: - PetscLib lib; // Initialize and finalize PETSc - - /// Mesh-global numbering for cells stored on this rank. Field3D cell_number; Field3D forward_cell_number; Field3D backward_cell_number; - - /// Evolving cells in the main domain. Region evolving_region; - - /// X-boundary cells on the inner radial boundary. Region xin_region; - - /// X-boundary cells on the outer radial boundary. Region xout_region; - - /// Y-boundary cells in the forward direction. Region yup_region; - - /// Y-boundary cells in the backward direction. Region ydown_region; +}; + +/// Mapping for forward or backward leg spaces. +class PetscLegMapping : public PetscIndexMapping { +public: + PetscLegMapping(int total_legs, std::vector local_leg_indices); +}; + +using PetscCellMappingPtr = std::shared_ptr; +using PetscLegMappingPtr = std::shared_ptr; + +/// Wrapper around a PETSc Vec with a compile-time space tag. +template +class PetscSpaceVector { + using UniqueVec = bout::petsc::UniqueVec; + +public: + PetscSpaceVector() = default; + + explicit PetscSpaceVector(std::shared_ptr mapping) + : mapping(std::move(mapping)), + vec(createVec(this->mapping->localSize(), this->mapping->globalSize())) {} + + PetscSpaceVector(std::shared_ptr mapping, UniqueVec vec) + : mapping(std::move(mapping)), vec(std::move(vec)) {} + + PetscSpaceVector(const PetscSpaceVector&) = delete; + PetscSpaceVector& operator=(const PetscSpaceVector&) = delete; + PetscSpaceVector(PetscSpaceVector&&) noexcept = default; + PetscSpaceVector& operator=(PetscSpaceVector&&) noexcept = default; + + Vec raw() const { return *this->vec; } + const std::shared_ptr& getMapping() const { return mapping; } + + PetscSpaceVector duplicate() const { + UniqueVec out{new Vec{nullptr}}; + BOUT_DO_PETSC(VecDuplicate(*this->vec, out.get())); + return PetscSpaceVector(this->mapping, std::move(out)); + } + + PetscSpaceVector operator*(BoutReal scalar) const { + auto out = this->duplicate(); + BOUT_DO_PETSC(VecScale(out.raw(), scalar)); + return out; + } + + PetscSpaceVector operator+(const PetscSpaceVector& rhs) const { + ASSERT0(mapping == rhs.mapping); + auto out = this->duplicate(); + BOUT_DO_PETSC(VecAXPY(out.raw(), 1.0, rhs.raw())); + return out; + } + + PetscSpaceVector operator-(const PetscSpaceVector& rhs) const { + ASSERT0(mapping == rhs.mapping); + auto out = this->duplicate(); + BOUT_DO_PETSC(VecAXPY(out.raw(), -1.0, rhs.raw())); + return out; + } + + static PetscSpaceVector pointwiseMultiply(const PetscSpaceVector& lhs, + const PetscSpaceVector& rhs) { + ASSERT0(lhs.mapping == rhs.mapping); + auto out = lhs.duplicate(); + BOUT_DO_PETSC(VecPointwiseMult(out.raw(), lhs.raw(), rhs.raw())); + return out; + } - /// First and one-past-last PETSc global rows owned by this rank. - PetscInt row_start{}, row_end{}; + static PetscSpaceVector reciprocal(const PetscSpaceVector& in) { + auto out = in.duplicate(); + BOUT_DO_PETSC(VecReciprocal(out.raw())); + return out; + } - /// Permutation matrix mapping mesh-global ordering to PETSc ordering. - Mat mat_mesh_to_petsc{nullptr}; +private: + static UniqueVec createVec(PetscInt local_size, PetscInt global_size) { + UniqueVec out{new Vec{nullptr}}; + BOUT_DO_PETSC(VecCreate(BoutComm::get(), out.get())); + BOUT_DO_PETSC(VecSetType(*out, VECMPI)); + BOUT_DO_PETSC(VecSetSizes(*out, local_size, global_size)); + return out; + } - /// Permutation matrix mapping PETSc ordering to mesh-global ordering. - Mat mat_petsc_to_mesh{nullptr}; + std::shared_ptr mapping; + UniqueVec vec; }; -/// Shared pointer to an immutable `PetscMapping`. -/// -/// The mapping is constructed once and then shared by operators that use it. -using PetscMappingPtr = std::shared_ptr; +// Tagged vectors provide compile-time check of compatible operations +using PetscCellVector = PetscSpaceVector; +using PetscForwardLegVector = PetscSpaceVector; +using PetscBackwardLegVector = PetscSpaceVector; -/// Linear operator on `Field3D` backed by a PETSc matrix. -/// -/// `PetscOperator` wraps a PETSc `Mat` whose rows and columns use the PETSc -/// numbering defined by a corresponding `PetscMapping`. Operators are usually -/// constructed from CSR arrays read from the mesh file and can then be: -/// - applied to a `Field3D`; -/// - composed with other operators; and -/// - added or subtracted. +PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, + const Field3D& field); +Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype); + +/// Typed linear operator between two spaces. +template class PetscOperator { - using UniqueVec = bout::petsc::UniqueVec; // unique_ptr to Vec - using UniqueMat = bout::petsc::UniqueMat; // unique_ptr to Mat + using UniqueMat = bout::petsc::UniqueMat; public: - /// Construct an operator from CSR data in mesh-global numbering. - /// - /// @param mapping Shared mapping between mesh and PETSc numbering. - /// @param rows CSR row offsets. - /// @param cols CSR column indices in mesh-global numbering. - /// @param weights CSR non-zero values. - /// - /// The CSR arrays define an operator in mesh-global numbering. The - /// constructor converts this representation into a PETSc matrix using - /// the supplied mapping. - PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, - Array weights); - - // No copying + PetscOperator() = default; + ~PetscOperator() = default; + + PetscOperator(std::shared_ptr out_mapping, + std::shared_ptr in_mapping, Array rows, + Array cols, Array weights) + : out_mapping(std::move(out_mapping)), in_mapping(std::move(in_mapping)) { + assembleFromCSR(rows, cols, weights); + } + + PetscOperator(std::shared_ptr out_mapping, + std::shared_ptr in_mapping, UniqueMat mat) + : out_mapping(std::move(out_mapping)), in_mapping(std::move(in_mapping)), + mat_operator(std::move(mat)) {} + + // Cannot be copied but can be moved PetscOperator(const PetscOperator&) = delete; PetscOperator& operator=(const PetscOperator&) = delete; - - // Allow Moving PetscOperator(PetscOperator&&) noexcept = default; PetscOperator& operator=(PetscOperator&&) noexcept = default; - /// Destroy the PETSc matrix and working vectors owned by this operator. - ~PetscOperator() = default; + PetscSpaceVector + operator()(const PetscSpaceVector& rhs) const { + ASSERT0(in_mapping == rhs.getMapping()); + PetscSpaceVector result(out_mapping); + BOUT_DO_PETSC(MatMult(*this->mat_operator, rhs.raw(), result.raw())); + return result; + } + + template >> + PetscSpaceVector operator()(const Field3D& rhs) const { + auto in = makePetscCellVector( + std::static_pointer_cast(in_mapping), rhs); + return (*this)(in); + } + + template >> + Field3D applyToField(const PetscSpaceVector& rhs, + const Field3D& prototype) const { + auto out = (*this)(rhs); + return toField3D(out, prototype); + } + + template + && std::is_same_v>> + Field3D operator()(const Field3D& rhs, const Field3D& prototype) const { + return applyToField((*this)(rhs), prototype); + } + + PetscOperator transpose() const { + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatTranspose(*mat_operator, MAT_INITIAL_MATRIX, out.get())); + return PetscOperator(in_mapping, out_mapping, + std::move(out)); + } + + PetscOperator operator+(const PetscOperator& rhs) const { + ASSERT0(out_mapping == rhs.out_mapping); + ASSERT0(in_mapping == rhs.in_mapping); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatAXPY(*out, 1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); + } + + PetscOperator operator-(const PetscOperator& rhs) const { + ASSERT0(out_mapping == rhs.out_mapping); + ASSERT0(in_mapping == rhs.in_mapping); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatAXPY(*out, -1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); + } - /// Create a diagonal operator - static PetscOperator diagonal(PetscMappingPtr mapping, const Field3D& f); - - /// Apply the operator to a field. - /// - /// @param rhs Input field, including any required `yup` and `ydown` - /// boundary values. - /// @return Result of applying the operator to `rhs`. - /// - /// The input field is packed into PETSc vector storage using the mapping, - /// multiplied by the PETSc matrix, then unpacked back into a `Field3D`. - Field3D operator()(const Field3D& rhs) const; - - /// Compose this operator with another operator. - /// - /// @param rhs Right-hand-side operator. - /// @return The composed operator `(*this) * rhs`. - /// - /// Both operators must use compatible mappings. - PetscOperator operator*(const PetscOperator& rhs) const; - - /// Add two operators. - /// - /// @param rhs Right-hand-side operator. - /// @return The sum of the two operators. - /// - /// Both operators must use compatible mappings. - PetscOperator operator+(const PetscOperator& rhs) const; - - /// Subtract one operator from another. - /// - /// @param rhs Right-hand-side operator. - /// @return The difference of the two operators. - /// - /// Both operators must use compatible mappings. - PetscOperator operator-(const PetscOperator& rhs) const; - - /// Calculate transpose as a new matrix - PetscOperator transpose() const; - - void view() const { MatView(*this->mat_operator, PETSC_VIEWER_STDOUT_WORLD); } - - /// Multiply by a scalar - /// This version allocates a new Mat and so is for lvalues. PetscOperator operator*(BoutReal scalar) const& { - UniqueMat mat{new Mat{}}; - BOUT_DO_PETSC(MatDuplicate(*mat_operator, MAT_COPY_VALUES, mat.get())); - BOUT_DO_PETSC(MatScale(*mat, scalar)); - return PetscOperator(this->mapping, std::move(mat)); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatScale(*out, scalar)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); } - /// Multiply by a scalar. - /// This version is for rvalue temporaries and modifies in-place. PetscOperator operator*(BoutReal scalar) && { BOUT_DO_PETSC(MatScale(*this->mat_operator, scalar)); return std::move(*this); } + Mat raw() const { return *this->mat_operator; } + + static PetscOperator diagonal(std::shared_ptr mapping, + const PetscSpaceVector& diag) { + static_assert(std::is_same_v, + "Diagonal operators require the same input and output space"); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatCreate(BoutComm::get(), out.get())); + BOUT_DO_PETSC(MatSetType(*out, MATMPIAIJ)); + BOUT_DO_PETSC(MatSetSizes(*out, mapping->localSize(), mapping->localSize(), + mapping->globalSize(), mapping->globalSize())); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*out, 1, nullptr, 0, nullptr)); + BOUT_DO_PETSC(MatDiagonalSet(*out, diag.raw(), INSERT_VALUES)); + BOUT_DO_PETSC(MatAssemblyBegin(*out, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*out, MAT_FINAL_ASSEMBLY)); + return PetscOperator(mapping, mapping, std::move(out)); + } + private: - /// Construct directly from an existing PETSc matrix. - /// - /// @param mapping Shared mapping between mesh and PETSc numbering. - /// @param mat PETSc matrix implementing the operator. - /// - /// This constructor is used internally when creating operators from - /// PETSc matrix algebra such as composition, addition, or subtraction. - PetscOperator(PetscMappingPtr mapping, UniqueMat mat) - : mapping(std::move(mapping)), mat_operator(std::move(mat)), - rhs_vec(createVec(this->mapping->localSize())), - result_vec(createVec(this->mapping->localSize())) {} - - /// Shared mapping between mesh and PETSc numbering. - PetscMappingPtr mapping; - - /// PETSc matrix implementing the operator. - UniqueMat mat_operator; - - /// Work vector holding the packed input field. - UniqueVec rhs_vec; - - /// Work vector holding the packed result. - UniqueVec result_vec; - - /// Copy contents of f, f.yup() and f.ydown() into vec. - /// Assumes that vec has already been allocated. - static void copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec); - - // Allocate MPI vector - static UniqueVec createVec(PetscInt local_size); + void assembleFromCSR(const Array& rows, const Array& cols, + const Array& weights) { + UniqueMat temp{new Mat{nullptr}}; + BOUT_DO_PETSC(MatCreate(BoutComm::get(), temp.get())); + BOUT_DO_PETSC(MatSetType(*temp, MATMPIAIJ)); + BOUT_DO_PETSC(MatSetSizes(*temp, out_mapping->localSize(), PETSC_DECIDE, + out_mapping->globalSize(), in_mapping->globalSize())); + + const auto& local_rows = out_mapping->localStoredIndices(); + for (PetscInt local_row = 0; local_row < static_cast(local_rows.size()); + ++local_row) { + const int stored_row = local_rows[local_row]; + if (stored_row < 0 || stored_row + 1 >= static_cast(rows.size())) { + continue; + } + const int start = rows[stored_row]; + const int end = rows[stored_row + 1]; + if (end <= start) { + continue; + } + const PetscInt petsc_row = out_mapping->rowStart() + local_row; + BOUT_DO_PETSC(MatSetValues(*temp, 1, &petsc_row, end - start, &cols[start], + &weights[start], INSERT_VALUES)); + } + BOUT_DO_PETSC(MatAssemblyBegin(*temp, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*temp, MAT_FINAL_ASSEMBLY)); + + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatMatMult(*temp, in_mapping->getPetscToStored(), MAT_INITIAL_MATRIX, + PETSC_DETERMINE, out.get())); + mat_operator = std::move(out); + } + + std::shared_ptr out_mapping; + std::shared_ptr in_mapping; + UniqueMat mat_operator{new Mat{nullptr}}; + + template + friend PetscOperator operator*(const PetscOperator& lhs, + const PetscOperator& rhs); }; -/// Collection of PETSc operators read from a mesh. -/// -/// `PetscOperators` owns the `PetscMapping` derived from the mesh numbering -/// fields and provides access to named operators stored in the mesh file -/// as CSR arrays. +template +PetscOperator +operator*(const PetscOperator& lhs, + const PetscOperator& rhs) { + ASSERT0(lhs.in_mapping == rhs.out_mapping); + typename PetscOperator::UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC( + MatMatMult(lhs.raw(), rhs.raw(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, out.get())); + return PetscOperator(lhs.out_mapping, rhs.in_mapping, + std::move(out)); +} + +using PetscForwardOperator = PetscOperator; +using PetscBackwardOperator = PetscOperator; +using PetscForwardToCellOperator = PetscOperator; +using PetscBackwardToCellOperator = PetscOperator; +using PetscCellOperator = PetscOperator; +using PetscForwardLegOperator = PetscOperator; +using PetscBackwardLegOperator = PetscOperator; + +/// Factory for operators defined on cell and leg spaces. class PetscOperators { public: - /// Construct PETSc operator support from a mesh. - /// - /// Reads the following mesh variables: - /// - `cell_number` - /// - `forward_cell_number` - /// - `backward_cell_number` - /// - /// and uses them to construct the shared `PetscMapping`. - /// - /// @param mesh Mesh from which numbering fields and operators are read. - PetscOperators(Mesh* mesh); - - /// Retrieve a named PETSc operator from the mesh. - /// - /// @param name Base name of the operator in the mesh file. - /// @return A `PetscOperator` constructed from the arrays - /// `{name}_rows`, `{name}_columns`, and `{name}_weights`. - /// - /// Throws `BoutException` if any of the required arrays are missing. - PetscOperator get(const std::string& name) const; - - /// Create a diagonal operator - PetscOperator diagonal(const Field3D& f) const { - return PetscOperator::diagonal(mapping, f); - } + explicit PetscOperators(Mesh* mesh); - /// A standard collection of operators - /// parallel to the magnetic field - struct Parallel { - PetscOperator Grad_par; ///< Gradient - PetscOperator Div_par; ///< Divergence - PetscOperator Div_par_Grad_par; ///< Diffusion - Field3D dV; ///< Cell volume + PetscForwardOperator forward() const; + PetscBackwardOperator backward() const; - PetscOperator Grad_minus; - PetscOperator Grad_plus; - PetscOperator Div_minus; - PetscOperator Div_plus; + PetscForwardOperator injectForward() const; + PetscBackwardOperator injectBackward() const; - PetscOperator Interp_minus; - PetscOperator Interp_plus; + PetscForwardLegVector forwardLegWeights() const; + PetscBackwardLegVector backwardLegWeights() const; + + PetscCellOperator diagonal(const Field3D& f) const; + PetscForwardLegOperator diagonal(const PetscForwardLegVector& v) const; + PetscBackwardLegOperator diagonal(const PetscBackwardLegVector& v) const; + + struct Parallel { + PetscForwardOperator Forward; + PetscBackwardOperator Backward; + PetscForwardOperator Inject_plus; + PetscBackwardOperator Inject_minus; + PetscForwardOperator Interp_plus; + PetscBackwardOperator Interp_minus; + PetscForwardOperator Grad_plus; + PetscBackwardOperator Grad_minus; + PetscForwardToCellOperator Div_minus; + PetscBackwardToCellOperator Div_plus; + PetscCellOperator Div_par_Grad_par; + Field3D dV; - /// This is a bilinear operator that requires multiple - /// matrix-vector products. Field3D Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const; }; - /// Calculate a set of parallel operators - /// - /// This assumes that "forward" and "backward" operators - /// are stored in the mesh. Parallel getParallel() const; private: - /// Read a `Field3D` from the mesh. - /// - /// @param mesh Mesh from which to read. - /// @param name Name of the field. - /// @return The requested `Field3D`. - /// - /// Throws `BoutException` if the field is not present. static Field3D meshGetField3D(Mesh* mesh, const std::string& name); - /// Read a one-dimensional `Array` from the mesh. - /// - /// @tparam T Array element type. - /// @param mesh Mesh from which to read. - /// @param name Name of the array. - /// @return The requested array. - /// - /// Throws `BoutException` if the array is not present. template - Array meshGetArray(Mesh* mesh, const std::string& name) const { + Array meshGetArray(const std::string& name) const { Array result; if (mesh->get(result, name) != 0) { - throw BoutException("PetscOperators requires Array '{}'", name); + throw BoutException("PetscOperators requires array '{}'", name); } return result; } - /// Mesh from which operator data are read. - Mesh* mesh; + int meshGetInt(const std::string& name) const; + + template + PetscOperator + buildInjection(const Field3D& interior_leg_number, const Field3D& boundary_leg_number, + std::shared_ptr leg_mapping) const; + + template + LegVector loadLegWeights(const std::string& name, + std::shared_ptr mapping) const { + auto values = meshGetArray(name); + LegVector out(mapping); + PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArray(out.raw(), &x)); + const auto& local = mapping->localStoredIndices(); + for (PetscInt i = 0; i < static_cast(local.size()); ++i) { + const int stored = local[i]; + x[i] = (stored >= 0 && stored < static_cast(values.size())) ? values[stored] + : 0.0; + } + BOUT_DO_PETSC(VecRestoreArray(out.raw(), &x)); + return out; + } - /// Shared mapping used by all operators created from this mesh. - PetscMappingPtr mapping; + Mesh* mesh{nullptr}; + PetscCellMappingPtr cell_mapping; + PetscLegMappingPtr forward_leg_mapping; + PetscLegMappingPtr backward_leg_mapping; + Field3D forward_leg_interior_number; + Field3D forward_leg_boundary_number; + Field3D backward_leg_interior_number; + Field3D backward_leg_boundary_number; }; -#else // BOUT_HAS_PETSC +#else #warning PETSc not available. No PetscOperators. -#endif // BOUT_HAS_PETSC +#endif diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 6767ebe716..d3682d2b00 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -2,67 +2,99 @@ #if BOUT_HAS_PETSC -#include "bout/array.hxx" #include "bout/assert.hxx" #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/field3d.hxx" -#include "bout/output.hxx" -#include "bout/output_bout_types.hxx" -#include "bout/petsc_interface.hxx" #include "bout/petsc_operators.hxx" -#include "bout/petsclib.hxx" #include "bout/region.hxx" +#include "bout/utils.hxx" + +#include #include #include #include +#include + +PetscIndexMapping::~PetscIndexMapping() { + if (mat_stored_to_petsc != nullptr) { + MatDestroy(&mat_stored_to_petsc); + } + if (mat_petsc_to_stored != nullptr) { + MatDestroy(&mat_petsc_to_stored); + } +} -Region PetscMapping::create_region(const Field3D& cell_number) { +PetscInt PetscIndexMapping::storedToPetsc(int stored_index) const { + auto it = local_stored_to_petsc.find(stored_index); + if (it == local_stored_to_petsc.end()) { + throw BoutException("Stored index {} is not owned locally", stored_index); + } + return it->second; +} + +void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, + const std::vector& stored_indices) { + global_size = nglobal; + local_stored_indices = stored_indices; + + BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat_stored_to_petsc)); + BOUT_DO_PETSC(MatSetSizes(mat_stored_to_petsc, nlocal, nlocal, nglobal, nglobal)); + BOUT_DO_PETSC(MatSetType(mat_stored_to_petsc, MATMPIAIJ)); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_stored_to_petsc, 1, nullptr, 1, nullptr)); + BOUT_DO_PETSC(MatGetOwnershipRange(mat_stored_to_petsc, &row_start, &row_end)); + ASSERT1(row_end - row_start == nlocal); + + const PetscScalar one = 1.0; + for (PetscInt i = 0; i < nlocal; ++i) { + const PetscInt petsc_row = row_start + i; + const PetscInt stored_col = stored_indices[i]; + ASSERT1(stored_col >= 0); + local_stored_to_petsc[static_cast(stored_col)] = petsc_row; + BOUT_DO_PETSC(MatSetValues(mat_stored_to_petsc, 1, &petsc_row, 1, &stored_col, &one, + INSERT_VALUES)); + } + BOUT_DO_PETSC(MatAssemblyBegin(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC( + MatTranspose(mat_stored_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_stored)); +} + +Region PetscCellMapping::create_region(const Field3D& cell_number) { Region::RegionIndices indices; - BOUT_FOR_SERIAL(i, cell_number.getRegion("RGN_NOBNDRY")) { - if (cell_number[i] > -1) { + BOUT_FOR(i, cell_number.getRegion("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { indices.push_back(i); } } return Region(indices); } -Region PetscMapping::create_region_xin(const Field3D& cell_number) { - Region::RegionIndices xin_indices; - const Mesh* mesh = cell_number.getMesh(); - if (mesh->firstX()) { - for (int i = 0; i < mesh->xstart; ++i) { - for (int j = mesh->ystart; j <= mesh->yend; ++j) { - for (int k = mesh->zstart; k <= mesh->zend; ++k) { - if (cell_number(i, j, k) > -1) { - xin_indices.push_back(cell_number.indexAt(i, j, k)); - } - } - } +Region PetscCellMapping::create_region_xin(const Field3D& cell_number) { + Region::RegionIndices indices; + const auto& region = cell_number.getRegion("RGN_INNER_X"); + BOUT_FOR(i, region) { + if (cell_number[i] >= 0) { + indices.push_back(i); } } - return Region(xin_indices); + return Region(indices); } -Region PetscMapping::create_region_xout(const Field3D& cell_number) { - Region::RegionIndices xout_indices; - const Mesh* mesh = cell_number.getMesh(); - if (mesh->lastX()) { - for (int i = mesh->xend + 1; i < mesh->LocalNx; ++i) { - for (int j = mesh->ystart; j <= mesh->yend; ++j) { - for (int k = mesh->zstart; k <= mesh->zend; ++k) { - if (cell_number(i, j, k) > -1) { - xout_indices.push_back(cell_number.indexAt(i, j, k)); - } - } - } +Region PetscCellMapping::create_region_xout(const Field3D& cell_number) { + Region::RegionIndices indices; + const auto& region = cell_number.getRegion("RGN_OUTER_X"); + BOUT_FOR(i, region) { + if (cell_number[i] >= 0) { + indices.push_back(i); } } - return Region(xout_indices); + return Region(indices); } -PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, - const Field3D& backward_cell_number) +PetscCellMapping::PetscCellMapping(const Field3D& cell_number, + const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells) : cell_number(cell_number), forward_cell_number(forward_cell_number), backward_cell_number(backward_cell_number), evolving_region(create_region(cell_number)), @@ -70,378 +102,261 @@ PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_ce xout_region(create_region_xout(cell_number)), yup_region(create_region(forward_cell_number)), ydown_region(create_region(backward_cell_number)) { - // Calculate size of each region - const unsigned int nlocal = this->localSize(); - output.write("Local {} : evolving {} xin {} xout {} yup {} ydown {}", nlocal, - evolving_region.size(), xin_region.size(), xout_region.size(), - yup_region.size(), ydown_region.size()); - - // Create a PETSc matrix - // Note: Numbering is different from that used in PetscVector / PetscMatrix - // because yup/down boundaries are included. - - // Renumbering matrix. - // Maps PETSc row (or column) indices to the global indices used in - // the mesh. This is needed because the PETSc indices depend on the - // number of processors. - BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat_mesh_to_petsc)); - BOUT_DO_PETSC( - MatSetSizes(mat_mesh_to_petsc, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - BOUT_DO_PETSC(MatSetType(mat_mesh_to_petsc, MATMPIAIJ)); - // Each row will have one non-zero entry, which could be in - // either the "diagonal" or "off-diagonal" block. - BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_mesh_to_petsc, 1, nullptr, 1, nullptr)); - - // Get the range of rows owned by this processor - BOUT_DO_PETSC(MatGetOwnershipRange(mat_mesh_to_petsc, &row_start, &row_end)); - output.write("Local row range: {} -> {}\n", row_start, row_end); - - // Iterate through regions in this order - this->map([&](PetscInt row, PetscInt col) { - // `row` is the PETSc index; `col` is the Mesh index - ASSERT1(row >= 0); - ASSERT1(col >= 0); - const PetscScalar ONE = 1.0; - BOUT_DO_PETSC(MatSetValues(mat_mesh_to_petsc, 1, &row, 1, &col, &ONE, INSERT_VALUES)); - }); - BOUT_DO_PETSC(MatAssemblyBegin(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY)); - BOUT_DO_PETSC(MatAssemblyEnd(mat_mesh_to_petsc, MAT_FINAL_ASSEMBLY)); + std::vector local_indices; + const auto push_region = [&](const Region& region, const Field3D& values) { + BOUT_FOR_SERIAL(i, region) { local_indices.push_back(ROUND(values[i])); } + }; + push_region(evolving_region, this->cell_number); + push_region(xin_region, this->cell_number); + push_region(xout_region, this->cell_number); + push_region(yup_region, this->forward_cell_number); + push_region(ydown_region, this->backward_cell_number); + buildPermutation(static_cast(local_indices.size()), total_cells, + local_indices); +} + +PetscLegMapping::PetscLegMapping(int total_legs, std::vector local_leg_indices) { + std::sort(local_leg_indices.begin(), local_leg_indices.end()); + local_leg_indices.erase(std::unique(local_leg_indices.begin(), local_leg_indices.end()), + local_leg_indices.end()); + buildPermutation(static_cast(local_leg_indices.size()), total_legs, + local_leg_indices); +} + +PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, + const Field3D& field) { + ASSERT1(field.hasParallelSlices()); + ASSERT1(field.yup().isAllocated()); + ASSERT1(field.ydown().isAllocated()); + PetscCellVector out(mapping); + PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArray(out.raw(), &x)); + mapping->mapLocalField([&](PetscInt row, const Ind3D& i) { x[row] = field[i]; }); + const auto& yup = field.yup(); + mapping->mapLocalYup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + const auto& ydown = field.ydown(); + mapping->mapLocalYdown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + BOUT_DO_PETSC(VecRestoreArray(out.raw(), &x)); + return out; +} - // Transpose to map petsc indices to mesh indices - MatTranspose(mat_mesh_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_mesh); +Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype) { + auto mapping = std::static_pointer_cast(vec.getMapping()); + Field3D result{emptyFrom(prototype)}; + const PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArrayRead(vec.raw(), &x)); + mapping->mapLocalField([&](PetscInt row, const Ind3D& i) { result[i] = x[row]; }); + mapping->mapLocalYup([&](PetscInt row, const Ind3D& i) { result.yup()[i] = x[row]; }); + mapping->mapLocalYdown( + [&](PetscInt row, const Ind3D& i) { result.ydown()[i] = x[row]; }); + BOUT_DO_PETSC(VecRestoreArrayRead(vec.raw(), &x)); + return result; } -PetscMapping::~PetscMapping() { - MatDestroy(&this->mat_mesh_to_petsc); - MatDestroy(&this->mat_petsc_to_mesh); +Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { + Field3D result; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires field '{}'", name); + } + return result; +} + +int PetscOperators::meshGetInt(const std::string& name) const { + int result{0}; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires int '{}'", name); + } + return result; } -PetscOperator::PetscOperator(PetscMappingPtr mapping, Array rows, Array cols, - Array weights) - : mapping(std::move(mapping)), mat_operator(new Mat{nullptr}), - rhs_vec(createVec(this->mapping->localSize())), - result_vec(createVec(this->mapping->localSize())) { - - output.write("PetscOperator Array sizes {} {} {}\n", rows.size(), cols.size(), - weights.size()); - - Mat mat{nullptr}; - BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat)); - BOUT_DO_PETSC(MatSetType(mat, MATMPIAIJ)); - const int nlocal = this->mapping->localSize(); - BOUT_DO_PETSC(MatSetSizes(mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - - // This reads CSR format but is defensive in handling negative indices and missing - // final 'row' entry. - this->mapping->map([&](PetscInt row, PetscInt mesh_index) { - if (mesh_index >= rows.size()) { - return; // No weights -> skip +PetscOperators::PetscOperators(Mesh* mesh) + : mesh(mesh), + cell_mapping(std::make_shared( + meshGetField3D(mesh, "cell_number"), + meshGetField3D(mesh, "forward_cell_number"), + meshGetField3D(mesh, "backward_cell_number"), meshGetInt("total_cells"))), + forward_leg_interior_number(meshGetField3D(mesh, "forward_leg_interior_number")), + forward_leg_boundary_number(meshGetField3D(mesh, "forward_leg_boundary_number")), + backward_leg_interior_number(meshGetField3D(mesh, "backward_leg_interior_number")), + backward_leg_boundary_number(meshGetField3D(mesh, "backward_leg_boundary_number")) { + std::vector local_forward_legs; + std::vector local_backward_legs; + cell_mapping->mapOwnedInteriorCells([&](PetscInt, const Ind3D& i, int) { + const int f_int = ROUND(forward_leg_interior_number[i]); + const int f_bnd = ROUND(forward_leg_boundary_number[i]); + const int b_int = ROUND(backward_leg_interior_number[i]); + const int b_bnd = ROUND(backward_leg_boundary_number[i]); + if (f_int >= 0) { + local_forward_legs.push_back(f_int); } - // Get the range of indices into columns and weights - const int start_ind = rows[mesh_index]; - if (start_ind < 0) { - return; // No entries (non-standard CSR) + if (f_bnd >= 0) { + local_forward_legs.push_back(f_bnd); } - int end_ind = - cols.size(); // End of the columns array (should be last element of rows) - for (int i = mesh_index + 1; i < rows.size(); ++i) { - // rows[i] can be -1 if no weights (non-standard CSR) - if (rows[i] > -1) { - // This is the next entry in the columns / weights array - end_ind = rows[i]; - break; - } + if (b_int >= 0) { + local_backward_legs.push_back(b_int); } - if (end_ind == start_ind) { - // Empty row in CSR format - return; + if (b_bnd >= 0) { + local_backward_legs.push_back(b_bnd); } - BOUT_DO_PETSC(MatSetValues(mat, 1, &row, end_ind - start_ind, &cols[start_ind], - &weights[start_ind], INSERT_VALUES)); }); - BOUT_DO_PETSC(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); - BOUT_DO_PETSC(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); - - // Row indices are PETSc indices but columns are mesh indices - // Multiply on the right by PetscToMesh. - BOUT_DO_PETSC(MatMatMult(mat, this->mapping->getPetscToMesh(), MAT_INITIAL_MATRIX, - PETSC_DETERMINE, this->mat_operator.get())); - // Destroy temporary matrix - BOUT_DO_PETSC(MatDestroy(&mat)); -} - -void PetscOperator::copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec) { - ASSERT1(f.hasParallelSlices()); - ASSERT1(f.yup().isAllocated()); - ASSERT1(f.ydown().isAllocated()); - - PetscScalar* x{nullptr}; - // Copy rows from Field3D cells - BOUT_DO_PETSC(VecGetArray(vec, &x)); - mapping->map_local_field([&](PetscInt row, const Ind3D& i) { x[row] = f[i]; }); - // Copy Yup / Ydown values from boundaries - const Field3D& yup = f.yup(); - mapping->map_local_yup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); - - const Field3D& ydown = f.ydown(); - mapping->map_local_ydown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); - - BOUT_DO_PETSC(VecRestoreArray(vec, &x)); + forward_leg_mapping = std::make_shared(meshGetInt("n_forward_legs"), + std::move(local_forward_legs)); + backward_leg_mapping = std::make_shared( + meshGetInt("n_backward_legs"), std::move(local_backward_legs)); } -PetscOperator::UniqueVec PetscOperator::createVec(PetscInt local_size) { - UniqueVec vec(new Vec{nullptr}); - BOUT_DO_PETSC(VecCreate(BoutComm::get(), vec.get())); - BOUT_DO_PETSC(VecSetType(*vec, VECMPI)); - BOUT_DO_PETSC(VecSetSizes(*vec, local_size, PETSC_DETERMINE)); - return vec; +PetscForwardOperator PetscOperators::forward() const { + return PetscForwardOperator( + forward_leg_mapping, cell_mapping, meshGetArray("forward_rows"), + meshGetArray("forward_columns"), meshGetArray("forward_weights")); } -PetscOperator PetscOperator::diagonal(PetscMappingPtr mapping, const Field3D& f) { - const UniqueVec diag{createVec(mapping->localSize())}; - copyToVec(mapping, f, *diag); +PetscBackwardOperator PetscOperators::backward() const { + return PetscBackwardOperator( + backward_leg_mapping, cell_mapping, meshGetArray("backward_rows"), + meshGetArray("backward_columns"), meshGetArray("backward_weights")); +} - UniqueMat mat(new Mat{nullptr}); - // Note: MatMatMult with diagonal and mpiaij not supported - // -> Create MPIAIJ +template +PetscOperator +PetscOperators::buildInjection(const Field3D& interior_leg_number, + const Field3D& boundary_leg_number, + std::shared_ptr leg_mapping) const { + using Op = PetscOperator; + bout::petsc::UniqueMat mat{new Mat{nullptr}}; BOUT_DO_PETSC(MatCreate(BoutComm::get(), mat.get())); BOUT_DO_PETSC(MatSetType(*mat, MATMPIAIJ)); - const PetscInt nlocal = mapping->localSize(); - BOUT_DO_PETSC(MatSetSizes(*mat, nlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE)); - BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*mat, 1, nullptr, 0, nullptr)); - BOUT_DO_PETSC(MatDiagonalSet(*mat, *diag, INSERT_VALUES)); + BOUT_DO_PETSC(MatSetSizes(*mat, leg_mapping->localSize(), cell_mapping->localSize(), + leg_mapping->globalSize(), cell_mapping->globalSize())); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*mat, 1, nullptr, 1, nullptr)); + const PetscScalar one = 1.0; + cell_mapping->mapOwnedInteriorCells([&](PetscInt cell_row, const Ind3D& i, int) { + const int interior = ROUND(interior_leg_number[i]); + const int boundary = ROUND(boundary_leg_number[i]); + if (interior >= 0) { + const PetscInt leg_row = leg_mapping->storedToPetsc(interior); + BOUT_DO_PETSC(MatSetValues(*mat, 1, &leg_row, 1, &cell_row, &one, INSERT_VALUES)); + } + if (boundary >= 0) { + const PetscInt leg_row = leg_mapping->storedToPetsc(boundary); + BOUT_DO_PETSC(MatSetValues(*mat, 1, &leg_row, 1, &cell_row, &one, INSERT_VALUES)); + } + }); BOUT_DO_PETSC(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); BOUT_DO_PETSC(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); - - return PetscOperator(std::move(mapping), std::move(mat)); + return Op(leg_mapping, cell_mapping, std::move(mat)); } -/// Perform operation -Field3D PetscOperator::operator()(const Field3D& rhs) const { - // Fill vec from rhs - copyToVec(this->mapping, rhs, *this->rhs_vec); - - // Perform Mat-Vec muliplication - BOUT_DO_PETSC(MatMult(*this->mat_operator, *rhs_vec, *result_vec)); +template PetscForwardOperator PetscOperators::buildInjection( + const Field3D&, const Field3D&, std::shared_ptr) const; +template PetscBackwardOperator PetscOperators::buildInjection( + const Field3D&, const Field3D&, std::shared_ptr) const; - // Copy result_vec into a Field3D - Field3D result{emptyFrom(rhs)}; // This allocates memory - const PetscScalar* r{nullptr}; - BOUT_DO_PETSC(VecGetArrayRead(*result_vec, &r)); - this->mapping->map_local_field( - [&](PetscInt row, const Ind3D& i) { result[i] = r[row]; }); - BOUT_DO_PETSC(VecRestoreArrayRead(*result_vec, &r)); - - return result; +PetscForwardOperator PetscOperators::injectForward() const { + return buildInjection( + forward_leg_interior_number, forward_leg_boundary_number, forward_leg_mapping); } -/// Operator composition -PetscOperator PetscOperator::operator*(const PetscOperator& rhs) const { - ASSERT0(this->mapping == rhs.mapping); - UniqueMat mat(new Mat{nullptr}); - BOUT_DO_PETSC(MatMatMult(*this->mat_operator, *rhs.mat_operator, MAT_INITIAL_MATRIX, - PETSC_DETERMINE, mat.get())); - return PetscOperator(this->mapping, std::move(mat)); +PetscBackwardOperator PetscOperators::injectBackward() const { + return buildInjection( + backward_leg_interior_number, backward_leg_boundary_number, backward_leg_mapping); } -/// Operator addition -PetscOperator PetscOperator::operator+(const PetscOperator& rhs) const { - ASSERT0(this->mapping == rhs.mapping); - UniqueMat mat(new Mat{nullptr}); - BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, mat.get())); - BOUT_DO_PETSC(MatAXPY(*mat, 1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); - return PetscOperator(this->mapping, std::move(mat)); +PetscForwardLegVector PetscOperators::forwardLegWeights() const { + return loadLegWeights("forward_leg_weights", + forward_leg_mapping); } -/// Operator subtraction -PetscOperator PetscOperator::operator-(const PetscOperator& rhs) const { - ASSERT0(this->mapping == rhs.mapping); - UniqueMat mat(new Mat{nullptr}); - BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, mat.get())); - BOUT_DO_PETSC(MatAXPY(*mat, -1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); - return PetscOperator(this->mapping, std::move(mat)); +PetscBackwardLegVector PetscOperators::backwardLegWeights() const { + return loadLegWeights("backward_leg_weights", + backward_leg_mapping); } -PetscOperator PetscOperator::transpose() const { - UniqueMat mat(new Mat{nullptr}); - BOUT_DO_PETSC(MatTranspose(*this->mat_operator, MAT_INITIAL_MATRIX, mat.get())); - return PetscOperator(this->mapping, std::move(mat)); +PetscCellOperator PetscOperators::diagonal(const Field3D& f) const { + return PetscCellOperator::diagonal(cell_mapping, makePetscCellVector(cell_mapping, f)); } -Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { - Field3D result; - if (mesh->get(result, name) != 0) { - throw BoutException("PetscOperators requires field '{}'", name); - } - return result; +PetscForwardLegOperator PetscOperators::diagonal(const PetscForwardLegVector& v) const { + return PetscForwardLegOperator::diagonal(forward_leg_mapping, v); } -PetscOperators::PetscOperators(Mesh* mesh) - : mesh(mesh), mapping(std::make_shared( - meshGetField3D(mesh, "cell_number"), - meshGetField3D(mesh, "forward_cell_number"), - meshGetField3D(mesh, "backward_cell_number"))) { - - int mesh_total_cells{0}; - if (mesh->get(mesh_total_cells, "total_cells") == 0) { - // Check total number of cells - if (mesh_total_cells != mapping->globalSize()) { - throw BoutException("Total cells in mesh {} doesn't match global mapping size {}", - mesh_total_cells, mapping->globalSize()); - } - } -} - -PetscOperator PetscOperators::get(const std::string& name) const { - return PetscOperator(this->mapping, this->meshGetArray(this->mesh, name + "_rows"), - this->meshGetArray(this->mesh, name + "_columns"), - this->meshGetArray(this->mesh, name + "_weights")); +PetscBackwardLegOperator PetscOperators::diagonal(const PetscBackwardLegVector& v) const { + return PetscBackwardLegOperator::diagonal(backward_leg_mapping, v); } PetscOperators::Parallel PetscOperators::getParallel() const { - // Read maps from the mesh - auto forward = this->get("forward"); - auto backward = this->get("backward"); - - // ---- Construct Grad_par ---- - // - // Create a parallel gradient operator by combining the parallel - // length dl = dy * sqrt(g_22) with forward & backward operators - auto* coords = this->mesh->getCoordinates(); + auto Forward = forward(); + auto Backward = backward(); + auto Inject_plus = injectForward(); + auto Inject_minus = injectBackward(); + + auto W_plus = diagonal(forwardLegWeights()); + auto W_minus = diagonal(backwardLegWeights()); + + auto* coords = mesh->getCoordinates(); Field3D dl = coords->dy * sqrt(coords->g_22); dl.splitParallelSlices(); dl.yup() = 0.0; dl.ydown() = 0.0; dl.applyParallelBoundary("parallel_neumann_o1"); - auto inv_2dl = 0.5 / dl; - inv_2dl.splitParallelSlices(); - inv_2dl.yup() = 0.0; - inv_2dl.ydown() = 0.0; - inv_2dl.applyParallelBoundary("parallel_neumann_o1"); + auto Interp_plus = W_plus * ((Forward + Inject_plus) * 0.5); + auto Interp_minus = W_minus * ((Backward + Inject_minus) * 0.5); - auto inv_2dl_op = this->diagonal(inv_2dl); + const auto dl_plus = Interp_plus(dl); + const auto dl_minus = Interp_minus(dl); + const auto inv_dl_plus = PetscForwardLegVector::reciprocal(dl_plus); + const auto inv_dl_minus = PetscBackwardLegVector::reciprocal(dl_minus); + const auto Inv_dl_plus = diagonal(inv_dl_plus); + const auto Inv_dl_minus = diagonal(inv_dl_minus); - auto Grad_par = inv_2dl_op * (forward - backward); + auto Grad_plus = Inv_dl_plus * (Forward - Inject_plus); + auto Grad_minus = Inv_dl_minus * (Inject_minus - Backward); - // ---- Construct Div_par ---- - // - // Use the Support Operator Method (SOM) to calculate - // Div_par from Grad_par and cell volumes. - - // Cell volume Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; dV.splitParallelSlices(); dV.yup() = 0.0; dV.ydown() = 0.0; dV.applyParallelBoundary("parallel_neumann_o1"); - // Fractional boundary cells - const auto forward_boundary_fraction = - meshGetField3D(this->mesh, "forward_boundary_fraction"); - const auto backward_boundary_fraction = - meshGetField3D(this->mesh, "backward_boundary_fraction"); - - BOUT_FOR(i, dV.getRegion("RGN_NOBNDRY")) { - dV.yup()[i.yp()] *= forward_boundary_fraction[i]; - dV.ydown()[i.ym()] *= backward_boundary_fraction[i]; - } - - const auto dV_op = this->diagonal(dV); - - Field3D neg_inv_dV = -1. / dV; + Field3D neg_inv_dV = -1.0 / dV; neg_inv_dV.splitParallelSlices(); neg_inv_dV.yup() = 0.0; neg_inv_dV.ydown() = 0.0; neg_inv_dV.applyParallelBoundary("parallel_neumann_o1"); - auto neg_inv_dV_op = this->diagonal(neg_inv_dV); - - auto Div_par = neg_inv_dV_op * Grad_par.transpose() * dV_op; - - // ---- Construct Div_par_Grad_par ---- - // - // Requires gradients between planes, at +1/2 and -1/2, and interpolation - // operators to calculate quantities between cells. - - // Identity operator - Field3D one{1.0}; - one.splitParallelSlices(); - one.yup() = 1.0; - one.ydown() = 1.0; - const auto identity = this->diagonal(one); - - // Interpolate at + 1/2 - auto interp_plus_op = (identity + forward) * 0.5; - - // dl averaged at +1/2 - const Field3D dl_plus = interp_plus_op(dl); - Field3D inv_dl_plus = 1. / dl_plus; - inv_dl_plus.splitParallelSlices(); - inv_dl_plus.yup() = 0.0; - inv_dl_plus.ydown() = 0.0; - inv_dl_plus.applyParallelBoundary("parallel_neumann_o1"); - const auto inv_dl_plus_op = this->diagonal(inv_dl_plus); - - // Gradient at + 1/2 - auto Grad_plus = inv_dl_plus_op * (forward - identity); - - // Divergence at -1/2 - auto Div_minus = neg_inv_dV_op * Grad_plus.transpose() * dV_op; - - // Interpolate at - 1/2 - auto interp_minus_op = (identity + backward) * 0.5; - - // dl averaged at -1/2 - const Field3D dl_minus = interp_minus_op(dl); - Field3D inv_dl_minus = 1. / dl_minus; - inv_dl_minus.splitParallelSlices(); - inv_dl_minus.yup() = 0.0; - inv_dl_minus.ydown() = 0.0; - inv_dl_minus.applyParallelBoundary("parallel_neumann_o1"); - const auto inv_dl_minus_op = this->diagonal(inv_dl_minus); - - // Gradient at - 1/2 - auto Grad_minus = inv_dl_minus_op * (identity - backward); - - // Divergence at +1/2 - auto Div_plus = neg_inv_dV_op * Grad_minus.transpose() * dV_op; - - // Div(Grad_par()) operator + const auto Neg_inv_dV = diagonal(neg_inv_dV); + + const auto dV_plus = Interp_plus(dV); + const auto dV_minus = Interp_minus(dV); + const auto DV_plus = diagonal(dV_plus); + const auto DV_minus = diagonal(dV_minus); + + auto Div_minus = Neg_inv_dV * Grad_plus.transpose() * DV_plus; + auto Div_plus = Neg_inv_dV * Grad_minus.transpose() * DV_minus; + auto Div_par_Grad_par = ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; - return Parallel{std::move(Grad_par), std::move(Div_par), - std::move(Div_par_Grad_par), std::move(dV), - std::move(Grad_minus), std::move(Grad_plus), - std::move(Div_minus), std::move(Div_plus), - std::move(interp_minus_op), std::move(interp_plus_op)}; + return Parallel{ + std::move(Forward), std::move(Backward), std::move(Inject_plus), + std::move(Inject_minus), std::move(Interp_plus), std::move(Interp_minus), + std::move(Grad_plus), std::move(Grad_minus), std::move(Div_minus), + std::move(Div_plus), std::move(Div_par_Grad_par), std::move(dV)}; } Field3D PetscOperators::Parallel::Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const { + const auto grad_plus = Grad_plus(f); + const auto grad_minus = Grad_minus(f); + const auto K_plus = Interp_plus(K); + const auto K_minus = Interp_minus(K); + + const auto flux_plus = PetscForwardLegVector::pointwiseMultiply(K_plus, grad_plus); + const auto flux_minus = PetscBackwardLegVector::pointwiseMultiply(K_minus, grad_minus); - // There are 4 matrix-vector products that could be performed in parallel. - // The best way to optimize this is probably to pack f and K into one Vec, - // and assemble a single sparse matrix that contains the four blocks: - // - // (Grad_plus 0 ) (grad_plus ) - // (Grad_minus 0 ) (f) -> (grad_minus ) - // ( 0 Interp_plus ) (K) (K_plus ) - // ( 0 Interp_minus) (K_minus ) - // - // For now we perform each operation in sequence. - - // Calculate gradients in + and - directions - const Field3D grad_plus = this->Grad_plus(f); - const Field3D grad_minus = this->Grad_minus(f); - - // Interpolate K to + and - locations - const Field3D K_plus = this->Interp_plus(K); - const Field3D K_minus = this->Interp_minus(K); - - // Calculate divergence - return (this->Div_minus(K_plus * grad_plus) + this->Div_plus(K_minus * grad_minus)) - * 0.5; + const auto div_minus = Div_minus(flux_plus); + const auto div_plus = Div_plus(flux_minus); + return (toField3D(div_minus, K) + toField3D(div_plus, K)) * 0.5; } -#endif // BOUT_HAS_PETSC +#endif From 6676ac9f19c87d9716d0b6e1fc4cb3548ef221b4 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 23 Mar 2026 21:23:49 -0700 Subject: [PATCH 22/26] PetscSpaceVector: Fix duplicate -> copy Use VecCopy to make a copy of the contents. Rename function from `duplicate` to `copy`. --- include/bout/petsc_operators.hxx | 38 ++++++---- src/mesh/petsc_operators.cxx | 124 +++++++++++++++++++++++-------- 2 files changed, 115 insertions(+), 47 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index b9d35a9793..ce69049542 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -181,28 +181,30 @@ public: Vec raw() const { return *this->vec; } const std::shared_ptr& getMapping() const { return mapping; } - PetscSpaceVector duplicate() const { + /// Make a copy of the vector + PetscSpaceVector copy() const { UniqueVec out{new Vec{nullptr}}; BOUT_DO_PETSC(VecDuplicate(*this->vec, out.get())); + BOUT_DO_PETSC(VecCopy(*this->vec, *out.get())); return PetscSpaceVector(this->mapping, std::move(out)); } PetscSpaceVector operator*(BoutReal scalar) const { - auto out = this->duplicate(); + auto out = this->copy(); BOUT_DO_PETSC(VecScale(out.raw(), scalar)); return out; } PetscSpaceVector operator+(const PetscSpaceVector& rhs) const { ASSERT0(mapping == rhs.mapping); - auto out = this->duplicate(); + auto out = this->copy(); BOUT_DO_PETSC(VecAXPY(out.raw(), 1.0, rhs.raw())); return out; } PetscSpaceVector operator-(const PetscSpaceVector& rhs) const { ASSERT0(mapping == rhs.mapping); - auto out = this->duplicate(); + auto out = this->copy(); BOUT_DO_PETSC(VecAXPY(out.raw(), -1.0, rhs.raw())); return out; } @@ -210,13 +212,13 @@ public: static PetscSpaceVector pointwiseMultiply(const PetscSpaceVector& lhs, const PetscSpaceVector& rhs) { ASSERT0(lhs.mapping == rhs.mapping); - auto out = lhs.duplicate(); + auto out = lhs.copy(); BOUT_DO_PETSC(VecPointwiseMult(out.raw(), lhs.raw(), rhs.raw())); return out; } static PetscSpaceVector reciprocal(const PetscSpaceVector& in) { - auto out = in.duplicate(); + auto out = in.copy(); BOUT_DO_PETSC(VecReciprocal(out.raw())); return out; } @@ -438,18 +440,24 @@ public: PetscBackwardLegOperator diagonal(const PetscBackwardLegVector& v) const; struct Parallel { - PetscForwardOperator Forward; - PetscBackwardOperator Backward; - PetscForwardOperator Inject_plus; - PetscBackwardOperator Inject_minus; - PetscForwardOperator Interp_plus; - PetscBackwardOperator Interp_minus; - PetscForwardOperator Grad_plus; + PetscCellOperator Grad_par; + PetscCellOperator Div_par; + PetscCellOperator Div_par_Grad_par; + Field3D dV; + PetscBackwardOperator Grad_minus; + PetscForwardOperator Grad_plus; PetscForwardToCellOperator Div_minus; PetscBackwardToCellOperator Div_plus; - PetscCellOperator Div_par_Grad_par; - Field3D dV; + + PetscBackwardOperator Inject_minus; + PetscForwardOperator Inject_plus; + + PetscBackwardOperator Interp_minus; + PetscForwardOperator Interp_plus; + + PetscForwardToCellOperator Restrict_minus; + PetscBackwardToCellOperator Restrict_plus; Field3D Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const; }; diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index d3682d2b00..e48643153c 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -6,6 +6,8 @@ #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/field3d.hxx" +#include "bout/output.hxx" +#include "bout/output_bout_types.hxx" #include "bout/petsc_operators.hxx" #include "bout/region.hxx" #include "bout/utils.hxx" @@ -38,10 +40,19 @@ void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, global_size = nglobal; local_stored_indices = stored_indices; + // Renumbering matrix. + // Maps PETSc row (or column) indices to the global indices stored in + // the mesh. This is needed because the PETSc indices depend on the + // number of processors. BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat_stored_to_petsc)); BOUT_DO_PETSC(MatSetSizes(mat_stored_to_petsc, nlocal, nlocal, nglobal, nglobal)); BOUT_DO_PETSC(MatSetType(mat_stored_to_petsc, MATMPIAIJ)); + + // Each row will have one non-zero entry, which could be in + // either the "diagonal" or "off-diagonal" block. BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_stored_to_petsc, 1, nullptr, 1, nullptr)); + + // Get the range of rows owned by this processor BOUT_DO_PETSC(MatGetOwnershipRange(mat_stored_to_petsc, &row_start, &row_end)); ASSERT1(row_end - row_start == nlocal); @@ -56,6 +67,7 @@ void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, } BOUT_DO_PETSC(MatAssemblyBegin(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); BOUT_DO_PETSC(MatAssemblyEnd(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); + // Transpose to map petsc indices to stored mesh indices BOUT_DO_PETSC( MatTranspose(mat_stored_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_stored)); } @@ -143,6 +155,10 @@ PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype) { auto mapping = std::static_pointer_cast(vec.getMapping()); Field3D result{emptyFrom(prototype)}; + result.splitParallelSlices(); + result.yup().allocate(); + result.ydown().allocate(); + const PetscScalar* x{nullptr}; BOUT_DO_PETSC(VecGetArrayRead(vec.raw(), &x)); mapping->mapLocalField([&](PetscInt row, const Ind3D& i) { result[i] = x[row]; }); @@ -150,6 +166,7 @@ Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype) { mapping->mapLocalYdown( [&](PetscInt row, const Ind3D& i) { result.ydown()[i] = x[row]; }); BOUT_DO_PETSC(VecRestoreArrayRead(vec.raw(), &x)); + return result; } @@ -286,34 +303,32 @@ PetscBackwardLegOperator PetscOperators::diagonal(const PetscBackwardLegVector& } PetscOperators::Parallel PetscOperators::getParallel() const { - auto Forward = forward(); - auto Backward = backward(); - auto Inject_plus = injectForward(); - auto Inject_minus = injectBackward(); - - auto W_plus = diagonal(forwardLegWeights()); - auto W_minus = diagonal(backwardLegWeights()); + // Primitive operators from the mesh / metadata + auto Forward = this->forward(); // L+ <- C + auto Backward = this->backward(); // L- <- C + auto Inject_plus = this->injectForward(); // L+ <- C + auto Inject_minus = this->injectBackward(); // L- <- C + + // Leg weights + const auto w_plus_vec = forwardLegWeights(); + const auto w_minus_vec = backwardLegWeights(); + const auto W_plus = diagonal(w_plus_vec); // L+ <- L+ + const auto W_minus = diagonal(w_minus_vec); // L- <- L- + + // Weighted adjoints: C <- L+ and C <- L- + auto Restrict_minus = Inject_plus.transpose() * W_plus; + auto Restrict_plus = Inject_minus.transpose() * W_minus; auto* coords = mesh->getCoordinates(); + + // Parallel spacing in cell space Field3D dl = coords->dy * sqrt(coords->g_22); dl.splitParallelSlices(); dl.yup() = 0.0; dl.ydown() = 0.0; dl.applyParallelBoundary("parallel_neumann_o1"); - auto Interp_plus = W_plus * ((Forward + Inject_plus) * 0.5); - auto Interp_minus = W_minus * ((Backward + Inject_minus) * 0.5); - - const auto dl_plus = Interp_plus(dl); - const auto dl_minus = Interp_minus(dl); - const auto inv_dl_plus = PetscForwardLegVector::reciprocal(dl_plus); - const auto inv_dl_minus = PetscBackwardLegVector::reciprocal(dl_minus); - const auto Inv_dl_plus = diagonal(inv_dl_plus); - const auto Inv_dl_minus = diagonal(inv_dl_minus); - - auto Grad_plus = Inv_dl_plus * (Forward - Inject_plus); - auto Grad_minus = Inv_dl_minus * (Inject_minus - Backward); - + // Cell volume Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; dV.splitParallelSlices(); dV.yup() = 0.0; @@ -325,23 +340,68 @@ PetscOperators::Parallel PetscOperators::getParallel() const { neg_inv_dV.yup() = 0.0; neg_inv_dV.ydown() = 0.0; neg_inv_dV.applyParallelBoundary("parallel_neumann_o1"); + + const auto DV = diagonal(dV); const auto Neg_inv_dV = diagonal(neg_inv_dV); - const auto dV_plus = Interp_plus(dV); - const auto dV_minus = Interp_minus(dV); - const auto DV_plus = diagonal(dV_plus); - const auto DV_minus = diagonal(dV_minus); + // Unweighted interpolation to +1/2 and -1/2 legs + auto Interp_plus = (Forward + Inject_plus) * 0.5; // L+ <- C + auto Interp_minus = (Backward + Inject_minus) * 0.5; // L- <- C + + // Leg-centered dl + const auto dl_plus = Interp_plus(dl); + const auto dl_minus = Interp_minus(dl); + + const auto inv_dl_plus = PetscForwardLegVector::reciprocal(dl_plus); + const auto inv_dl_minus = PetscBackwardLegVector::reciprocal(dl_minus); + + const auto Inv_dl_plus = diagonal(inv_dl_plus); // L+ <- L+ + const auto Inv_dl_minus = diagonal(inv_dl_minus); // L- <- L- + + // Half-step gradients + auto Grad_plus = Inv_dl_plus * (Forward - Inject_plus); // L+ <- C + auto Grad_minus = Inv_dl_minus * (Inject_minus - Backward); // L- <- C + + // Cell-centered central gradient + auto Grad_par = + ((Restrict_minus * Grad_plus) + (Restrict_plus * Grad_minus)) * 0.5; // C <- C + + // Leg-centered volumes from unweighted interpolation + const auto dV_plus = Interp_plus(dV); // L+ + const auto dV_minus = Interp_minus(dV); // L- + + // Physical leg volumes include the leg weights + const auto dV_leg_plus = PetscForwardLegVector::pointwiseMultiply(dV_plus, w_plus_vec); + const auto dV_leg_minus = + PetscBackwardLegVector::pointwiseMultiply(dV_minus, w_minus_vec); + + const auto DV_leg_plus = diagonal(dV_leg_plus); // L+ <- L+ + const auto DV_leg_minus = diagonal(dV_leg_minus); // L- <- L- + + // Support-operator divergence from half-step gradients + auto Div_minus = Neg_inv_dV * Grad_plus.transpose() * DV_leg_plus; // C <- L+ + auto Div_plus = Neg_inv_dV * Grad_minus.transpose() * DV_leg_minus; // C <- L- + + // Cell-centered divergence paired with the cell-centered gradient + auto Div_par = Neg_inv_dV * Grad_par.transpose() * DV; // C <- C + + // Diffusion operator assembled from the half-step pieces + auto Div_par_Grad_par = + ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; // C <- C + + return Parallel{// Cell <- Cell operators + std::move(Grad_par), std::move(Div_par), std::move(Div_par_Grad_par), + std::move(dV), + + // Leg <-> Cell operators + std::move(Grad_minus), std::move(Grad_plus), std::move(Div_minus), + std::move(Div_plus), - auto Div_minus = Neg_inv_dV * Grad_plus.transpose() * DV_plus; - auto Div_plus = Neg_inv_dV * Grad_minus.transpose() * DV_minus; + std::move(Inject_minus), std::move(Inject_plus), - auto Div_par_Grad_par = ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; + std::move(Interp_minus), std::move(Interp_plus), - return Parallel{ - std::move(Forward), std::move(Backward), std::move(Inject_plus), - std::move(Inject_minus), std::move(Interp_plus), std::move(Interp_minus), - std::move(Grad_plus), std::move(Grad_minus), std::move(Div_minus), - std::move(Div_plus), std::move(Div_par_Grad_par), std::move(dV)}; + std::move(Restrict_minus), std::move(Restrict_plus)}; } Field3D PetscOperators::Parallel::Div_par_K_Grad_par(const Field3D& K, From b1b0dc0868459624f3db4e4a8444358296d03c12 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 8 Apr 2026 15:33:51 -0700 Subject: [PATCH 23/26] PetscOperators: Fix tests, add outputs Still debugging --- .../fci-operators/fci_operators_example.cxx | 17 +++++++++-- examples/fci-operators/makeplots.py | 9 ++++++ include/bout/petsc_operators.hxx | 6 ++++ src/mesh/petsc_operators.cxx | 3 +- tests/unit/mesh/test_petsc_operators.cxx | 30 ++++++++++--------- 5 files changed, 47 insertions(+), 18 deletions(-) diff --git a/examples/fci-operators/fci_operators_example.cxx b/examples/fci-operators/fci_operators_example.cxx index 682d7c287c..66c28d3b6b 100644 --- a/examples/fci-operators/fci_operators_example.cxx +++ b/examples/fci-operators/fci_operators_example.cxx @@ -28,17 +28,28 @@ int main(int argc, char** argv) { dump["f"] = f; dump["dV"] = parallel.dV; - dump["grad_par_op"] = parallel.Grad_par(f); + // Test1 : (parallel.Restrict_minus * parallel.Inject_plus) is the identity in the interior (not X boundaries) + + dump["grad_par_op"] = parallel.Grad_par(f, f); + Field3D forward{0.0}; + BOUT_FOR(i, forward.getRegion("RGN_NOBNDRY")) { forward[i] = f.yup()[i.yp()]; } + + dump["forward"] = forward; + + dump["forward_op"] = (parallel.Restrict_minus * parallel.Forward)(f, f); + dump["grad_par_yud"] = Grad_par(f); - dump["div_par_op"] = parallel.Div_par(f); + // Test2 : Sum of dV * Div_par over domain is zero + + dump["div_par_op"] = parallel.Div_par(f, f); auto* coords = bout::globals::mesh->getCoordinates(); bout::globals::mesh->communicate(coords->Bxy); coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); dump["div_par_yud"] = Div_par(f); - dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); + dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann, f); dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); bout::writeDefaultOutputFile(dump); diff --git a/examples/fci-operators/makeplots.py b/examples/fci-operators/makeplots.py index b10a5829f3..53c3d38f05 100644 --- a/examples/fci-operators/makeplots.py +++ b/examples/fci-operators/makeplots.py @@ -4,6 +4,9 @@ path = "data" +forward_op = collect("forward_op", path=path) +forward = collect("forward", path=path) + grad_par_op = collect("grad_par_op", path=path) grad_par_yud = collect("grad_par_yud", path=path) @@ -47,6 +50,12 @@ def plot_comparison(a, alabel, b, blabel): return fig +fig = plot_comparison( + forward_op[2:-2, 2, :], "Operator", forward[2:-2, 2, :], "Yup/down" +) +fig.suptitle("Forward") +fig.tight_layout() + fig = plot_comparison( grad_par_op[2:-2, 2, :], "Operator", grad_par_yud[2:-2, 2, :], "Yup/down" ) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index ce69049542..03befb4d9d 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -8,6 +8,8 @@ #include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/mesh.hxx" +#include "bout/output.hxx" +#include "bout/output_bout_types.hxx" #include "bout/petsc_interface.hxx" #include "bout/petsclib.hxx" #include "bout/region.hxx" @@ -189,6 +191,7 @@ public: return PetscSpaceVector(this->mapping, std::move(out)); } + /// Multiply by a scalar, returning a new PetscSpaceVector PetscSpaceVector operator*(BoutReal scalar) const { auto out = this->copy(); BOUT_DO_PETSC(VecScale(out.raw(), scalar)); @@ -450,6 +453,9 @@ public: PetscForwardToCellOperator Div_minus; PetscBackwardToCellOperator Div_plus; + PetscBackwardOperator Backward; + PetscForwardOperator Forward; + PetscBackwardOperator Inject_minus; PetscForwardOperator Inject_plus; diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index e48643153c..a63407ec4a 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -397,7 +397,8 @@ PetscOperators::Parallel PetscOperators::getParallel() const { std::move(Grad_minus), std::move(Grad_plus), std::move(Div_minus), std::move(Div_plus), - std::move(Inject_minus), std::move(Inject_plus), + std::move(Backward), std::move(Forward), std::move(Inject_minus), + std::move(Inject_plus), std::move(Interp_minus), std::move(Interp_plus), diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 219a667bac..923b17ea84 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -22,7 +22,7 @@ using PetscMappingTest = FakeMeshFixture; TEST_F(PetscMappingTest, create_region_empty) { const Field3D cell_number{-1.0}; // No cells >= 0 - auto rgn = PetscMapping::create_region(cell_number); + auto rgn = PetscCellMapping::create_region(cell_number); ASSERT_EQ(0, rgn.size()); } @@ -35,7 +35,7 @@ TEST_F(PetscMappingTest, create_region) { cell_number(1, 1, 0) = 1; // In domain cell_number(0, 1, 1) = 2; // Xin boundary - auto rgn = PetscMapping::create_region(cell_number); + auto rgn = PetscCellMapping::create_region(cell_number); ASSERT_EQ(1, rgn.size()); const Ind3D expected_ind = cell_number.indexAt(1, 1, 0); @@ -50,7 +50,7 @@ TEST_F(PetscMappingTest, create_region_xin) { cell_number(1, 1, 0) = 1; // In domain cell_number(0, 1, 1) = 2; // Xin boundary - auto rgn = PetscMapping::create_region_xin(cell_number); + auto rgn = PetscCellMapping::create_region_xin(cell_number); ASSERT_EQ(1, rgn.size()); const Ind3D expected_ind = cell_number.indexAt(0, 1, 1); @@ -69,10 +69,12 @@ TEST_F(PetscMappingTest, mapping) { const Field3D forward_cell_number{-1.0}; const Field3D backward_cell_number{-1.0}; - const PetscMapping mapping(cell_number, forward_cell_number, backward_cell_number); + const PetscCellMapping mapping(cell_number, forward_cell_number, backward_cell_number, + 2); // Two cells: one evolving and one in xin ASSERT_EQ(2, mapping.localSize()); + ASSERT_EQ(2, mapping.globalSize()); } using PetscOperatorTest = FakeMeshFixture; @@ -87,15 +89,15 @@ TEST_F(PetscOperatorTest, identity) { cell_number(1, 2, 0) = 1; cell_number(1, 2, 1) = 2; - auto mapping = std::make_shared(cell_number, forward_cell_number, - backward_cell_number); + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, 3); ASSERT_EQ(3, mapping->localSize()); - const auto rows = Array::fromValues({0, 1, 2}); + const auto rows = Array::fromValues({0, 1, 2, 3}); const auto cols = Array::fromValues({0, 1, 2}); const auto weights = Array::fromValues({1.0, 1.0, 1.0}); - const PetscOperator identity(mapping, rows, cols, weights); + const PetscCellOperator identity(mapping, mapping, rows, cols, weights); Field3D value{0.0}; value(1, 1, 0) = 10; @@ -105,7 +107,7 @@ TEST_F(PetscOperatorTest, identity) { value.yup() = -1.0; value.ydown() = -1.0; - const Field3D result = identity(value); + const Field3D result = identity(value, value); EXPECT_EQ(10, result(1, 1, 0)); EXPECT_EQ(21, result(1, 2, 0)); @@ -125,15 +127,15 @@ TEST_F(PetscOperatorTest, identity_yupdown) { forward_cell_number(1, 2, 0) = 3; backward_cell_number(1, 1, 2) = 4; - auto mapping = std::make_shared(cell_number, forward_cell_number, - backward_cell_number); + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, 5); ASSERT_EQ(5, mapping->localSize()); - const auto rows = Array::fromValues({0, 1, 2}); + const auto rows = Array::fromValues({0, 1, 2, 3}); const auto cols = Array::fromValues({0, 1, 2}); const auto weights = Array::fromValues({1.0, 1.0, 1.0}); - const PetscOperator identity(mapping, rows, cols, weights); + const PetscCellOperator identity(mapping, mapping, rows, cols, weights); Field3D value{0.0}; value(1, 1, 0) = 10; @@ -143,7 +145,7 @@ TEST_F(PetscOperatorTest, identity_yupdown) { value.yup() = -1.0; value.ydown() = -1.0; - const Field3D result = identity(value); + const Field3D result = identity(value, value); EXPECT_EQ(10, result(1, 1, 0)); EXPECT_EQ(21, result(1, 2, 0)); From 7f3b584122693da82902f2ca97fee91653e0539d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 8 Apr 2026 22:50:41 -0700 Subject: [PATCH 24/26] PetscOperators: Fix operator() and applyToField issues `applyToField` was ambiguous in meaning apply and convert to Field, and `operator()(Field3D, Field3D)` ended up applying the operator twice. Simplified the application to vectors and Field3D objects. --- .../fci-operators/fci_operators_example.cxx | 10 +++++--- include/bout/petsc_operators.hxx | 25 ++++++++++--------- src/mesh/petsc_operators.cxx | 2 +- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/examples/fci-operators/fci_operators_example.cxx b/examples/fci-operators/fci_operators_example.cxx index 66c28d3b6b..3366622ff0 100644 --- a/examples/fci-operators/fci_operators_example.cxx +++ b/examples/fci-operators/fci_operators_example.cxx @@ -24,32 +24,34 @@ int main(int argc, char** argv) { bout::globals::mesh->communicate(f_neumann); f_neumann.applyParallelBoundary("parallel_neumann_o1"); + auto forward_op = parallel.Restrict_minus * parallel.Forward; + Options dump; dump["f"] = f; dump["dV"] = parallel.dV; // Test1 : (parallel.Restrict_minus * parallel.Inject_plus) is the identity in the interior (not X boundaries) - dump["grad_par_op"] = parallel.Grad_par(f, f); + dump["grad_par_op"] = parallel.Grad_par(f); Field3D forward{0.0}; BOUT_FOR(i, forward.getRegion("RGN_NOBNDRY")) { forward[i] = f.yup()[i.yp()]; } dump["forward"] = forward; - dump["forward_op"] = (parallel.Restrict_minus * parallel.Forward)(f, f); + dump["forward_op"] = forward_op(f); dump["grad_par_yud"] = Grad_par(f); // Test2 : Sum of dV * Div_par over domain is zero - dump["div_par_op"] = parallel.Div_par(f, f); + dump["div_par_op"] = parallel.Div_par(f); auto* coords = bout::globals::mesh->getCoordinates(); bout::globals::mesh->communicate(coords->Bxy); coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); dump["div_par_yud"] = Div_par(f); - dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann, f); + dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); bout::writeDefaultOutputFile(dump); diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 03befb4d9d..e2502dd723 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -275,6 +275,8 @@ public: PetscOperator(PetscOperator&&) noexcept = default; PetscOperator& operator=(PetscOperator&&) noexcept = default; + /// Apply operator to a vector in the input space, + /// returning a vector in the output space. PetscSpaceVector operator()(const PetscSpaceVector& rhs) const { ASSERT0(in_mapping == rhs.getMapping()); @@ -283,27 +285,26 @@ public: return result; } - template >> + /// Apply operator to a Field3D, returning a vector. + /// Enable if input is Cell Space and output is not. + template + && !std::is_same_v>> PetscSpaceVector operator()(const Field3D& rhs) const { auto in = makePetscCellVector( std::static_pointer_cast(in_mapping), rhs); return (*this)(in); } - template >> - Field3D applyToField(const PetscSpaceVector& rhs, - const Field3D& prototype) const { - auto out = (*this)(rhs); - return toField3D(out, prototype); - } - + /// Apply operator to a Field3D, returning a Field3D. + /// Enable if both input and output mapppings are in Cell Space template && std::is_same_v>> - Field3D operator()(const Field3D& rhs, const Field3D& prototype) const { - return applyToField((*this)(rhs), prototype); + Field3D operator()(const Field3D& rhs) const { + auto in = makePetscCellVector( + std::static_pointer_cast(in_mapping), rhs); + return toField3D((*this)(in), rhs); } PetscOperator transpose() const { diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index a63407ec4a..0fe564bcd3 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -52,7 +52,7 @@ void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, // either the "diagonal" or "off-diagonal" block. BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_stored_to_petsc, 1, nullptr, 1, nullptr)); - // Get the range of rows owned by this processor + // Get the range of global rows owned by this processor BOUT_DO_PETSC(MatGetOwnershipRange(mat_stored_to_petsc, &row_start, &row_end)); ASSERT1(row_end - row_start == nlocal); From af4a5b33d26c13c2825693d691fa574000245e23 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 8 Apr 2026 23:30:10 -0700 Subject: [PATCH 25/26] PetscOperators: Fix unit tests, add operators PetscSpaceVector operators on temporaries can modify in-place. --- include/bout/petsc_operators.hxx | 25 +++++++++++++++++++++--- tests/unit/mesh/test_petsc_operators.cxx | 4 ++-- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index e2502dd723..bc19b73d18 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -192,26 +192,44 @@ public: } /// Multiply by a scalar, returning a new PetscSpaceVector - PetscSpaceVector operator*(BoutReal scalar) const { + PetscSpaceVector operator*(BoutReal scalar) const& { auto out = this->copy(); BOUT_DO_PETSC(VecScale(out.raw(), scalar)); return out; } - PetscSpaceVector operator+(const PetscSpaceVector& rhs) const { + /// Multiply by a scalar, modifying a temporary in-place + PetscSpaceVector operator*(BoutReal scalar) && { + BOUT_DO_PETSC(VecScale(this->raw(), scalar)); + return std::move(*this); + } + + PetscSpaceVector operator+(const PetscSpaceVector& rhs) const& { ASSERT0(mapping == rhs.mapping); auto out = this->copy(); BOUT_DO_PETSC(VecAXPY(out.raw(), 1.0, rhs.raw())); return out; } - PetscSpaceVector operator-(const PetscSpaceVector& rhs) const { + PetscSpaceVector operator+(const PetscSpaceVector& rhs) && { + ASSERT0(mapping == rhs.mapping); + BOUT_DO_PETSC(VecAXPY(this->raw(), 1.0, rhs.raw())); + return std::move(*this); + } + + PetscSpaceVector operator-(const PetscSpaceVector& rhs) const& { ASSERT0(mapping == rhs.mapping); auto out = this->copy(); BOUT_DO_PETSC(VecAXPY(out.raw(), -1.0, rhs.raw())); return out; } + PetscSpaceVector operator-(const PetscSpaceVector& rhs) && { + ASSERT0(mapping == rhs.mapping); + BOUT_DO_PETSC(VecAXPY(this->raw(), -1.0, rhs.raw())); + return std::move(*this); + } + static PetscSpaceVector pointwiseMultiply(const PetscSpaceVector& lhs, const PetscSpaceVector& rhs) { ASSERT0(lhs.mapping == rhs.mapping); @@ -355,6 +373,7 @@ public: BOUT_DO_PETSC(MatSetType(*out, MATMPIAIJ)); BOUT_DO_PETSC(MatSetSizes(*out, mapping->localSize(), mapping->localSize(), mapping->globalSize(), mapping->globalSize())); + // Note: Off-diagonal blocks are empty because this is a square matrix BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*out, 1, nullptr, 0, nullptr)); BOUT_DO_PETSC(MatDiagonalSet(*out, diag.raw(), INSERT_VALUES)); BOUT_DO_PETSC(MatAssemblyBegin(*out, MAT_FINAL_ASSEMBLY)); diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 923b17ea84..2ba63f4528 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -107,7 +107,7 @@ TEST_F(PetscOperatorTest, identity) { value.yup() = -1.0; value.ydown() = -1.0; - const Field3D result = identity(value, value); + const Field3D result = identity(value); EXPECT_EQ(10, result(1, 1, 0)); EXPECT_EQ(21, result(1, 2, 0)); @@ -145,7 +145,7 @@ TEST_F(PetscOperatorTest, identity_yupdown) { value.yup() = -1.0; value.ydown() = -1.0; - const Field3D result = identity(value, value); + const Field3D result = identity(value); EXPECT_EQ(10, result(1, 1, 0)); EXPECT_EQ(21, result(1, 2, 0)); From c949d7987109fb5aff9235c60fd9b8050c4a83cd Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 00:03:35 -0700 Subject: [PATCH 26/26] PetscOperators: Add docstrings Coauthor: Claude --- include/bout/petsc_operators.hxx | 593 ++++++++++++++++++++++++++++--- 1 file changed, 537 insertions(+), 56 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index bc19b73d18..0d9c6c6912 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -26,21 +26,42 @@ #include #include -/// Tag type for the cell space C. +/// @brief Tag type identifying the cell space C. +/// +/// Used as a compile-time template parameter to distinguish vectors and operators +/// that live in the cell space from those in the forward or backward leg spaces. +/// Mismatched space tags produce a compile error rather than a runtime failure. struct CellSpaceTag {}; -/// Tag type for the forward leg space L+. +/// @brief Tag type identifying the forward leg space L+. +/// +/// Operators and vectors tagged with this type map to or from the forward (y+1) +/// leg space. See @ref BackwardLegSpaceTag and @ref CellSpaceTag. struct ForwardLegSpaceTag {}; -/// Tag type for the backward leg space L-. +/// @brief Tag type identifying the backward leg space L-. +/// +/// Operators and vectors tagged with this type map to or from the backward (y-1) +/// leg space. See @ref ForwardLegSpaceTag and @ref CellSpaceTag. struct BackwardLegSpaceTag {}; -/// Mapping between stored indices in a mesh file and PETSc row ownership. +/// @brief Bidirectional index mapping between mesh-file stored numbering and PETSc +/// row ownership. +/// +/// Every vector space (cell space, forward leg space, backward leg space) is defined +/// by a global stored numbering written to the mesh file by the Python weights module. +/// PETSc distributes rows contiguously across MPI ranks, so the local PETSc row range +/// [rowStart(), rowEnd()) does not in general coincide with the stored indices. +/// +/// This class maintains: +/// - the list of stored indices owned by the local rank, in PETSc row order; +/// - a hash map for O(1) stored → PETSc lookup; +/// - a permutation matrix @c mat_stored_to_petsc (PETSc ordering → stored numbering) +/// and its transpose @c mat_petsc_to_stored (stored numbering → PETSc ordering), +/// used to translate operator column indices during matrix assembly. /// -/// Every space is defined by a global stored numbering written to the mesh file. -/// PETSc distributes rows contiguously across MPI ranks, so this class stores the -/// local owned stored indices in row order and a permutation matrix that maps -/// from stored numbering to PETSc numbering. +/// Subclasses specialise the constructor to populate these structures for a particular +/// space from mesh metadata. class PetscIndexMapping { public: PetscIndexMapping() = default; @@ -51,54 +72,137 @@ public: PetscIndexMapping(PetscIndexMapping&&) = delete; PetscIndexMapping& operator=(PetscIndexMapping&&) = delete; - /// Number of rows owned locally by this MPI rank. + /// @brief Number of rows owned locally by this MPI rank. PetscInt localSize() const { return static_cast(local_stored_indices.size()); } - /// Global number of rows in this space. + /// @brief Total number of rows in this space across all MPI ranks. PetscInt globalSize() const { return global_size; } - /// First global PETSc row owned by this MPI rank. + /// @brief Global PETSc index of the first row owned by this MPI rank. PetscInt rowStart() const { return row_start; } - /// One-past-last global PETSc row owned by this MPI rank. + /// @brief Global PETSc index one past the last row owned by this MPI rank. + /// + /// The locally owned rows are the half-open interval [rowStart(), rowEnd()). PetscInt rowEnd() const { return row_end; } - /// Stored indices owned locally in PETSc row order. + /// @brief Stored mesh-file indices owned locally, listed in PETSc row order. + /// + /// Entry @c i corresponds to global PETSc row rowStart() + i. const std::vector& localStoredIndices() const { return local_stored_indices; } - /// Map a locally owned stored index to a PETSc global row. + /// @brief Convert a locally owned stored index to its global PETSc row index. + /// + /// @param stored_index A mesh-file stored index that must be owned by this rank. + /// @returns The corresponding global PETSc row index. + /// @throws BoutException if @p stored_index is not owned locally. PetscInt storedToPetsc(int stored_index) const; - /// Permutation matrix mapping PETSc ordering to stored numbering. + /// @brief Return the permutation matrix that maps PETSc row ordering to stored + /// mesh-file numbering. + /// + /// This is the transpose of @c mat_stored_to_petsc and is used in + /// @ref PetscOperator::assembleFromCSR to translate stored column indices into + /// PETSc column indices via a post-multiplication. Mat getPetscToStored() const { return mat_petsc_to_stored; } protected: + /// @brief BOUT++ PETSc library handle; ensures PETSc is initialised. PetscLib lib; + /// @brief Construct the permutation matrices and populate internal lookup structures. + /// + /// Assigns contiguous PETSc rows starting from the rank-local offset determined by + /// PETSc's ownership partitioning. Each stored index in @p stored_indices maps to + /// one PETSc row, in the order given. + /// + /// @param nlocal Number of rows owned by this rank. + /// @param nglobal Total number of rows across all ranks. + /// @param stored_indices Mesh-file stored indices for the locally owned rows, in + /// the order they should appear in the PETSc vector. void buildPermutation(PetscInt nlocal, PetscInt nglobal, const std::vector& stored_indices); - PetscInt global_size{0}; - PetscInt row_start{0}; - PetscInt row_end{0}; - std::vector local_stored_indices; - std::unordered_map local_stored_to_petsc; - Mat mat_stored_to_petsc{nullptr}; - Mat mat_petsc_to_stored{nullptr}; + PetscInt global_size{0}; ///< Total global row count. + PetscInt row_start{0}; ///< First global PETSc row owned locally. + PetscInt row_end{0}; ///< One-past-last global PETSc row owned locally. + std::vector local_stored_indices; ///< Stored indices in local PETSc row order. + std::unordered_map local_stored_to_petsc; ///< Stored → PETSc lookup. + Mat mat_stored_to_petsc{nullptr}; ///< Permutation: stored numbering → PETSc ordering. + Mat mat_petsc_to_stored{nullptr}; ///< Permutation: PETSc ordering → stored numbering. }; -/// Cell-space mapping between stored cell numbers, PETSc numbering, and Field3D. +/// @brief Cell-space index mapping that bridges stored cell numbers, PETSc row +/// ownership, and BOUT++ Field3D index arithmetic. +/// +/// The cell space C contains five disjoint subsets of cells, appended in order: +/// 1. Evolving interior cells (the primary unknowns). +/// 2. Inner radial X-boundary cells. +/// 3. Outer radial X-boundary cells. +/// 4. Forward parallel boundary cells (one virtual cell per interior source cell +/// whose forward map intersects the radial boundary). +/// 5. Backward parallel boundary cells (analogous for the backward map). +/// +/// The local PETSc vector stores these in the same order. The @c mapLocal* methods +/// iterate over subsets 1–3 (interior + X-boundary) while @c mapLocalYup and +/// @c mapLocalYdown cover subsets 4 and 5 respectively. +/// +/// @c mapOwnedInteriorCells iterates only the evolving interior cells (subset 1) and +/// provides global PETSc row indices, making it suitable for assembling matrices. class PetscCellMapping : public PetscIndexMapping { public: + /// @brief Construct the cell-space mapping from mesh metadata fields. + /// + /// Populates all five cell subsets and calls @ref buildPermutation with the + /// concatenated list of stored indices. + /// + /// @param cell_number Field3D of stored cell numbers for interior and + /// X-boundary cells. Negative entries are skipped. + /// @param forward_cell_number Field3D of stored cell numbers for forward parallel + /// boundary virtual cells. Negative entries are skipped. + /// @param backward_cell_number Field3D of stored cell numbers for backward parallel + /// boundary virtual cells. Negative entries are skipped. + /// @param total_cells Total number of cells in the global cell space, + /// as computed by the Python weights module. PetscCellMapping(const Field3D& cell_number, const Field3D& forward_cell_number, const Field3D& backward_cell_number, int total_cells); + /// @brief Build the region of evolving interior cells from a cell-number field. + /// + /// Iterates RGN_NOBNDRY and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of all valid interior cells. static Region create_region(const Field3D& cell_number); + + /// @brief Build the region of inner radial X-boundary cells. + /// + /// Iterates RGN_INNER_X and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of valid inner-X boundary cells. static Region create_region_xin(const Field3D& cell_number); + + /// @brief Build the region of outer radial X-boundary cells. + /// + /// Iterates RGN_OUTER_X and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of valid outer-X boundary cells. static Region create_region_xout(const Field3D& cell_number); + /// @brief Iterate over locally owned evolving interior cells, providing global + /// PETSc row indices. + /// + /// Calls @p func(petsc_row, field_index, stored_cell_number) for each cell in + /// the evolving region that is owned by this MPI rank. @p petsc_row is the global + /// PETSc index and is suitable for use as a row or column argument to MatSetValues. + /// + /// @tparam Function Callable with signature + /// (PetscInt row, Ind3D i, int stored) -> void. + /// @param func Callback invoked for each owned interior cell. template void mapOwnedInteriorCells(Function func) const { PetscInt row = rowStart(); @@ -108,6 +212,17 @@ public: } } + /// @brief Iterate over all locally stored field cells (interior + X-boundary), + /// providing local PETSc vector indices. + /// + /// Covers the evolving interior region, inner-X boundary, and outer-X boundary + /// in that order. The local index starts at zero on each rank and increases + /// contiguously. Use this to read from or write to a local PETSc Vec array + /// obtained via VecGetArray. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i) -> void. + /// @param func Callback invoked for each locally stored field cell. template void mapLocalField(Function func) const { PetscInt row = 0; @@ -121,6 +236,16 @@ public: } } + /// @brief Iterate over locally stored forward parallel boundary cells, providing + /// local PETSc vector indices. + /// + /// The local index starts immediately after the field cells covered by + /// @ref mapLocalField. The Field3D index passed to @p func is the y+1 neighbour + /// of the source cell (i.yp()), matching how field.yup() is indexed. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i_yp) -> void. + /// @param func Callback invoked for each locally stored forward boundary cell. template void mapLocalYup(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); @@ -130,6 +255,16 @@ public: } } + /// @brief Iterate over locally stored backward parallel boundary cells, providing + /// local PETSc vector indices. + /// + /// The local index starts immediately after the forward boundary cells covered by + /// @ref mapLocalYup. The Field3D index passed to @p func is the y-1 neighbour + /// of the source cell (i.ym()), matching how field.ydown() is indexed. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i_ym) -> void. + /// @param func Callback invoked for each locally stored backward boundary cell. template void mapLocalYdown(Function func) const { PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() @@ -141,37 +276,78 @@ public: } private: - Field3D cell_number; - Field3D forward_cell_number; - Field3D backward_cell_number; - Region evolving_region; - Region xin_region; - Region xout_region; - Region yup_region; - Region ydown_region; + Field3D cell_number; ///< Stored cell numbers for interior/X-boundary cells. + Field3D + forward_cell_number; ///< Stored cell numbers for forward boundary virtual cells. + Field3D + backward_cell_number; ///< Stored cell numbers for backward boundary virtual cells. + Region evolving_region; ///< Interior evolving cells. + Region xin_region; ///< Inner radial X-boundary cells. + Region xout_region; ///< Outer radial X-boundary cells. + Region yup_region; ///< Forward parallel boundary virtual cells. + Region ydown_region; ///< Backward parallel boundary virtual cells. }; -/// Mapping for forward or backward leg spaces. +/// @brief Index mapping for forward (L+) or backward (L-) leg spaces. +/// +/// Each interior source cell contributes one or two leg rows depending on whether +/// its mapped sub-samples are entirely interior, entirely boundary, or mixed. +/// The leg indices are assigned by the Python weights module and stored in the mesh +/// file as Field3D arrays (e.g. forward_leg_interior_number). +/// +/// On construction, duplicate indices are removed and the remaining indices are +/// sorted before being passed to @ref buildPermutation. class PetscLegMapping : public PetscIndexMapping { public: + /// @brief Construct a leg-space mapping. + /// + /// @param total_legs Global number of leg rows in this space. + /// @param local_leg_indices Stored leg indices owned by this MPI rank. Duplicates + /// are removed internally before building the permutation. PetscLegMapping(int total_legs, std::vector local_leg_indices); }; +/// @brief Shared-pointer alias for a const PetscCellMapping. using PetscCellMappingPtr = std::shared_ptr; + +/// @brief Shared-pointer alias for a const PetscLegMapping. using PetscLegMappingPtr = std::shared_ptr; -/// Wrapper around a PETSc Vec with a compile-time space tag. +/// @brief Type-safe wrapper around a PETSc Vec belonging to a particular vector space. +/// +/// The compile-time @p SpaceTag (one of @ref CellSpaceTag, @ref ForwardLegSpaceTag, +/// @ref BackwardLegSpaceTag) prevents accidental arithmetic between vectors from +/// incompatible spaces. All binary operations assert that both operands share the +/// same mapping pointer, catching dimension mismatches at runtime. +/// +/// Vectors are non-copyable but movable, matching the ownership semantics of the +/// underlying PETSc Vec. +/// +/// @tparam SpaceTag One of CellSpaceTag, ForwardLegSpaceTag, BackwardLegSpaceTag. template class PetscSpaceVector { using UniqueVec = bout::petsc::UniqueVec; public: + /// @brief Default constructor; produces an empty (unusable) vector. PetscSpaceVector() = default; + /// @brief Construct a zero-initialised vector for the given space. + /// + /// Allocates a PETSc MPI vector sized according to @p mapping->localSize() and + /// @p mapping->globalSize(). + /// + /// @param mapping Shared ownership of the space mapping. explicit PetscSpaceVector(std::shared_ptr mapping) : mapping(std::move(mapping)), vec(createVec(this->mapping->localSize(), this->mapping->globalSize())) {} + /// @brief Construct from an existing mapping and a pre-built PETSc Vec. + /// + /// Takes ownership of @p vec. + /// + /// @param mapping Shared ownership of the space mapping. + /// @param vec An already-assembled PETSc Vec. Ownership is transferred. PetscSpaceVector(std::shared_ptr mapping, UniqueVec vec) : mapping(std::move(mapping)), vec(std::move(vec)) {} @@ -180,10 +356,15 @@ public: PetscSpaceVector(PetscSpaceVector&&) noexcept = default; PetscSpaceVector& operator=(PetscSpaceVector&&) noexcept = default; + /// @brief Return the underlying PETSc Vec handle (non-owning). Vec raw() const { return *this->vec; } + + /// @brief Return the shared space mapping. const std::shared_ptr& getMapping() const { return mapping; } - /// Make a copy of the vector + /// @brief Return a deep copy of this vector. + /// + /// Allocates a new PETSc Vec with the same size and mapping, then copies values. PetscSpaceVector copy() const { UniqueVec out{new Vec{nullptr}}; BOUT_DO_PETSC(VecDuplicate(*this->vec, out.get())); @@ -191,7 +372,9 @@ public: return PetscSpaceVector(this->mapping, std::move(out)); } - /// Multiply by a scalar, returning a new PetscSpaceVector + /// @brief Return a new vector equal to @c (*this) * @p scalar. + /// + /// @param scalar Scalar multiplier. PetscSpaceVector operator*(BoutReal scalar) const& { auto out = this->copy(); BOUT_DO_PETSC(VecScale(out.raw(), scalar)); @@ -204,6 +387,9 @@ public: return std::move(*this); } + /// @brief Return a new vector equal to @c (*this) + @p rhs. + /// + /// @param rhs Must share the same mapping as @c *this. PetscSpaceVector operator+(const PetscSpaceVector& rhs) const& { ASSERT0(mapping == rhs.mapping); auto out = this->copy(); @@ -217,6 +403,9 @@ public: return std::move(*this); } + /// @brief Return a new vector equal to @c (*this) - @p rhs. + /// + /// @param rhs Must share the same mapping as @c *this. PetscSpaceVector operator-(const PetscSpaceVector& rhs) const& { ASSERT0(mapping == rhs.mapping); auto out = this->copy(); @@ -230,6 +419,11 @@ public: return std::move(*this); } + /// @brief Return a new vector whose entries are the element-wise product + /// @c lhs[i] * rhs[i]. + /// + /// @param lhs Left operand; must share a mapping with @p rhs. + /// @param rhs Right operand. static PetscSpaceVector pointwiseMultiply(const PetscSpaceVector& lhs, const PetscSpaceVector& rhs) { ASSERT0(lhs.mapping == rhs.mapping); @@ -238,6 +432,9 @@ public: return out; } + /// @brief Return a new vector whose entries are @c 1 / in[i]. + /// + /// @param in Input vector. No entry may be zero. static PetscSpaceVector reciprocal(const PetscSpaceVector& in) { auto out = in.copy(); BOUT_DO_PETSC(VecReciprocal(out.raw())); @@ -245,6 +442,10 @@ public: } private: + /// @brief Allocate a new zeroed PETSc MPI vector of the given size. + /// + /// @param local_size Number of locally owned entries. + /// @param global_size Total number of entries across all ranks. static UniqueVec createVec(PetscInt local_size, PetscInt global_size) { UniqueVec out{new Vec{nullptr}}; BOUT_DO_PETSC(VecCreate(BoutComm::get(), out.get())); @@ -253,28 +454,77 @@ private: return out; } - std::shared_ptr mapping; - UniqueVec vec; + std::shared_ptr mapping; ///< Space mapping (shared ownership). + UniqueVec vec; ///< Underlying PETSc Vec (owned). }; -// Tagged vectors provide compile-time check of compatible operations +// Tagged vector aliases — these provide compile-time checks of compatible operations. +/// @brief PETSc vector in the cell space C. using PetscCellVector = PetscSpaceVector; +/// @brief PETSc vector in the forward leg space L+. using PetscForwardLegVector = PetscSpaceVector; +/// @brief PETSc vector in the backward leg space L-. using PetscBackwardLegVector = PetscSpaceVector; +/// @brief Populate a PetscCellVector from a Field3D and its parallel slices. +/// +/// The cell vector is filled in the order defined by @ref PetscCellMapping: +/// evolving interior cells (from @p field), inner/outer X-boundary cells +/// (from @p field), forward boundary virtual cells (from @p field.yup()), and +/// backward boundary virtual cells (from @p field.ydown()). +/// +/// @param mapping Shared cell-space mapping. +/// @param field Source field; must have parallel slices split and allocated. +/// @returns A fully assembled PetscCellVector. PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, const Field3D& field); + +/// @brief Extract a Field3D (with parallel slices) from a PetscCellVector. +/// +/// Reverses @ref makePetscCellVector: interior and X-boundary values are written to +/// @c result, forward boundary values to @c result.yup(), and backward boundary +/// values to @c result.ydown(). +/// +/// @param vec Source vector; must use a PetscCellMapping. +/// @param prototype Field3D whose mesh, coordinates, and metadata are copied to +/// initialise the result. +/// @returns A newly allocated Field3D with split parallel slices. Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype); -/// Typed linear operator between two spaces. +/// @brief Type-safe sparse linear operator mapping from one space to another. +/// +/// Wraps a PETSc Mat and enforces compatible space pairings at compile time via +/// the @p OutSpaceTag and @p InSpaceTag template parameters. The underlying matrix +/// stores values in PETSc's internal ordering; translation from mesh-file stored +/// column indices to PETSc column indices is performed during construction via a +/// post-multiplication by @ref PetscIndexMapping::getPetscToStored. +/// +/// Operators are non-copyable but movable. Binary operations (addition, subtraction, +/// scalar multiplication, composition) return new operators without modifying +/// their operands, unless the rvalue-qualified @c operator*(BoutReal) overload is +/// selected, which scales in place. +/// +/// @tparam OutSpaceTag Tag for the output (row) space. +/// @tparam InSpaceTag Tag for the input (column) space. template class PetscOperator { using UniqueMat = bout::petsc::UniqueMat; public: + /// @brief Default constructor; produces an empty (unusable) operator. PetscOperator() = default; ~PetscOperator() = default; + /// @brief Construct an operator from a CSR-format sparse matrix stored in the mesh. + /// + /// Assembles the PETSc matrix from stored (non-PETSc) row and column indices by + /// calling @ref assembleFromCSR. + /// + /// @param out_mapping Shared mapping for the output (row) space. + /// @param in_mapping Shared mapping for the input (column) space. + /// @param rows CSR row-pointer array (size = number_of_rows + 1). + /// @param cols CSR column-index array, using stored (mesh-file) indices. + /// @param weights CSR nonzero-value array. PetscOperator(std::shared_ptr out_mapping, std::shared_ptr in_mapping, Array rows, Array cols, Array weights) @@ -282,19 +532,29 @@ public: assembleFromCSR(rows, cols, weights); } + /// @brief Construct an operator from an already-assembled PETSc Mat. + /// + /// Takes ownership of @p mat. + /// + /// @param out_mapping Shared mapping for the output (row) space. + /// @param in_mapping Shared mapping for the input (column) space. + /// @param mat Pre-assembled PETSc Mat. Ownership is transferred. PetscOperator(std::shared_ptr out_mapping, std::shared_ptr in_mapping, UniqueMat mat) : out_mapping(std::move(out_mapping)), in_mapping(std::move(in_mapping)), mat_operator(std::move(mat)) {} - // Cannot be copied but can be moved + // Operators are non-copyable (PETSc matrices are expensive to copy) but movable. PetscOperator(const PetscOperator&) = delete; PetscOperator& operator=(const PetscOperator&) = delete; PetscOperator(PetscOperator&&) noexcept = default; PetscOperator& operator=(PetscOperator&&) noexcept = default; - /// Apply operator to a vector in the input space, - /// returning a vector in the output space. + /// @brief Apply the operator to a typed vector: @c result = A * rhs. + /// + /// @param rhs Input vector; must share the same mapping as this operator's + /// input space. + /// @returns A new output-space vector holding the result. PetscSpaceVector operator()(const PetscSpaceVector& rhs) const { ASSERT0(in_mapping == rhs.getMapping()); @@ -303,8 +563,14 @@ public: return result; } - /// Apply operator to a Field3D, returning a vector. - /// Enable if input is Cell Space and output is not. + /// @brief Apply the operator to a Field3D: converts to a cell vector then applies. + /// + /// Only available when @p InSpaceTag is @ref CellSpaceTag and @p OutSpaceTag is + /// not @ref CellSpaceTag. The Field3D is packed into a PetscCellVector via + /// @ref makePetscCellVector before applying the matrix. + /// + /// @param rhs Source field; must have parallel slices split and allocated. + /// @returns A new output-space vector holding the result. template && !std::is_same_v>> @@ -325,6 +591,11 @@ public: return toField3D((*this)(in), rhs); } + /// @brief Return the transpose operator mapping from OutSpaceTag to InSpaceTag. + /// + /// Allocates a new PETSc Mat via MatTranspose (MAT_INITIAL_MATRIX). The mappings + /// are swapped: the transposed operator's input mapping is this operator's output + /// mapping, and vice versa. PetscOperator transpose() const { UniqueMat out{new Mat{nullptr}}; BOUT_DO_PETSC(MatTranspose(*mat_operator, MAT_INITIAL_MATRIX, out.get())); @@ -332,6 +603,13 @@ public: std::move(out)); } + /// @brief Return a new operator equal to @c (*this) + @p rhs. + /// + /// Both operands must share the same input and output mappings. Uses + /// UNKNOWN_NONZERO_PATTERN, so the sparsity structure of the result may differ + /// from either operand. + /// + /// @param rhs Operator to add; must have the same in/out mappings. PetscOperator operator+(const PetscOperator& rhs) const { ASSERT0(out_mapping == rhs.out_mapping); ASSERT0(in_mapping == rhs.in_mapping); @@ -341,6 +619,12 @@ public: return PetscOperator(out_mapping, in_mapping, std::move(out)); } + /// @brief Return a new operator equal to @c (*this) - @p rhs. + /// + /// Both operands must share the same input and output mappings. Uses + /// UNKNOWN_NONZERO_PATTERN; see @ref operator+. + /// + /// @param rhs Operator to subtract; must have the same in/out mappings. PetscOperator operator-(const PetscOperator& rhs) const { ASSERT0(out_mapping == rhs.out_mapping); ASSERT0(in_mapping == rhs.in_mapping); @@ -350,6 +634,11 @@ public: return PetscOperator(out_mapping, in_mapping, std::move(out)); } + /// @brief Return a new operator equal to @c (*this) * @p scalar (lvalue overload). + /// + /// Duplicates the matrix before scaling, leaving @c *this unchanged. + /// + /// @param scalar Scalar multiplier. PetscOperator operator*(BoutReal scalar) const& { UniqueMat out{new Mat{nullptr}}; BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); @@ -357,13 +646,32 @@ public: return PetscOperator(out_mapping, in_mapping, std::move(out)); } + /// @brief Scale this operator in place and return it (rvalue overload). + /// + /// More efficient than the lvalue overload when @c *this is a temporary, + /// as no matrix duplication is needed. + /// + /// @param scalar Scalar multiplier. PetscOperator operator*(BoutReal scalar) && { BOUT_DO_PETSC(MatScale(*this->mat_operator, scalar)); return std::move(*this); } + /// @brief Return the underlying PETSc Mat handle (non-owning). + /// + /// Intended for use by the free @c operator* (matrix composition) and other + /// low-level callers that need direct PETSc access. Mat raw() const { return *this->mat_operator; } + /// @brief Construct a diagonal operator from a vector of diagonal entries. + /// + /// Only available when @p OutSpaceTag == @p InSpaceTag (square, same-space + /// operator). Builds an MPIAIJ matrix with one nonzero per row via + /// MatDiagonalSet. + /// + /// @param mapping Shared mapping for both the row and column space. + /// @param diag Vector of diagonal values; must share @p mapping. + /// @returns A new diagonal operator. static PetscOperator diagonal(std::shared_ptr mapping, const PetscSpaceVector& diag) { static_assert(std::is_same_v, @@ -382,6 +690,18 @@ public: } private: + /// @brief Assemble the operator matrix from a mesh-file CSR representation. + /// + /// Iterates over the locally owned output rows (identified via + /// @c out_mapping->localStoredIndices()), looks up each row's nonzeros in the + /// CSR arrays, and inserts them using global PETSc row indices. Column indices + /// are given in stored (mesh-file) numbering; the final matrix is obtained by + /// post-multiplying the assembled temporary by + /// @c in_mapping->getPetscToStored() to convert to PETSc column ordering. + /// + /// @param rows CSR row-pointer array (size = total_stored_rows + 1). + /// @param cols CSR column-index array in stored (mesh-file) numbering. + /// @param weights CSR nonzero-value array. void assembleFromCSR(const Array& rows, const Array& cols, const Array& weights) { UniqueMat temp{new Mat{nullptr}}; @@ -409,21 +729,34 @@ private: BOUT_DO_PETSC(MatAssemblyBegin(*temp, MAT_FINAL_ASSEMBLY)); BOUT_DO_PETSC(MatAssemblyEnd(*temp, MAT_FINAL_ASSEMBLY)); + // Post-multiply by the stored→PETSc permutation to convert column numbering. UniqueMat out{new Mat{nullptr}}; BOUT_DO_PETSC(MatMatMult(*temp, in_mapping->getPetscToStored(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, out.get())); mat_operator = std::move(out); } - std::shared_ptr out_mapping; - std::shared_ptr in_mapping; - UniqueMat mat_operator{new Mat{nullptr}}; + std::shared_ptr out_mapping; ///< Row-space mapping. + std::shared_ptr in_mapping; ///< Column-space mapping. + UniqueMat mat_operator{new Mat{nullptr}}; ///< Underlying PETSc Mat. template friend PetscOperator operator*(const PetscOperator& lhs, const PetscOperator& rhs); }; +/// @brief Compose two compatible operators: @c result = lhs * rhs. +/// +/// The output space of @p rhs must match the input space of @p lhs, enforced at +/// compile time by the @p MidSpaceTag parameter. The resulting operator maps from +/// @p rhs's input space to @p lhs's output space. +/// +/// @tparam OutSpaceTag Output space of @p lhs and of the result. +/// @tparam MidSpaceTag Input space of @p lhs and output space of @p rhs. +/// @tparam InSpaceTag Input space of @p rhs and of the result. +/// @param lhs Left (outer) operator. +/// @param rhs Right (inner) operator. +/// @returns A new operator representing the composition. template PetscOperator operator*(const PetscOperator& lhs, @@ -436,63 +769,185 @@ operator*(const PetscOperator& lhs, std::move(out)); } +/// @brief Operator mapping from cell space C to forward leg space L+. using PetscForwardOperator = PetscOperator; +/// @brief Operator mapping from cell space C to backward leg space L-. using PetscBackwardOperator = PetscOperator; +/// @brief Operator mapping from forward leg space L+ to cell space C. using PetscForwardToCellOperator = PetscOperator; +/// @brief Operator mapping from backward leg space L- to cell space C. using PetscBackwardToCellOperator = PetscOperator; +/// @brief Operator mapping within the cell space C. using PetscCellOperator = PetscOperator; +/// @brief Operator mapping within the forward leg space L+. using PetscForwardLegOperator = PetscOperator; +/// @brief Operator mapping within the backward leg space L-. using PetscBackwardLegOperator = PetscOperator; -/// Factory for operators defined on cell and leg spaces. +/// @brief Factory that constructs cell-space and leg-space operators from mesh metadata. +/// +/// On construction, reads all required fields and integer scalars from the mesh and +/// sets up the @ref PetscCellMapping, forward @ref PetscLegMapping, and backward +/// @ref PetscLegMapping. Individual primitive operators (interpolation, injection, +/// leg weights) are constructed on demand. The high-level set of parallel operators +/// needed for Support-Operator-Method diffusion is assembled by @ref getParallel. class PetscOperators { public: + /// @brief Construct the operator factory from a BOUT++ mesh. + /// + /// Reads the following fields from @p mesh: + /// - @c cell_number, @c forward_cell_number, @c backward_cell_number (Field3D) + /// - @c total_cells, @c n_forward_legs, @c n_backward_legs (int) + /// - @c forward_leg_interior_number, @c forward_leg_boundary_number (Field3D) + /// - @c backward_leg_interior_number, @c backward_leg_boundary_number (Field3D) + /// + /// @param mesh Non-owning pointer to the BOUT++ mesh. Must outlive this object. + /// @throws BoutException if any required field or scalar is missing from the mesh. explicit PetscOperators(Mesh* mesh); + /// @brief Build the forward interpolation operator F: L+ ← C. + /// + /// Maps cell-space values to forward leg space using the CSR arrays stored in the + /// mesh (@c forward_rows, @c forward_columns, @c forward_weights). PetscForwardOperator forward() const; + + /// @brief Build the backward interpolation operator B: L- ← C. + /// + /// Maps cell-space values to backward leg space using the CSR arrays stored in the + /// mesh (@c backward_rows, @c backward_columns, @c backward_weights). PetscBackwardOperator backward() const; + /// @brief Build the forward injection operator I+: L+ ← C. + /// + /// Places the value of each interior source cell into its corresponding forward leg + /// row(s), with weight 1. Boundary leg rows point to the virtual boundary cell in C. PetscForwardOperator injectForward() const; + + /// @brief Build the backward injection operator I-: L- ← C. + /// + /// Analogous to @ref injectForward for the backward leg space. PetscBackwardOperator injectBackward() const; + /// @brief Load the forward leg-weight vector w+: L+. + /// + /// Each entry holds the quadrature weight of the corresponding forward leg row, + /// equal to the interior or boundary fraction of that leg's sub-sample integral. PetscForwardLegVector forwardLegWeights() const; + + /// @brief Load the backward leg-weight vector w-: L-. + /// + /// Analogous to @ref forwardLegWeights for the backward leg space. PetscBackwardLegVector backwardLegWeights() const; + /// @brief Build a diagonal cell-space operator from a Field3D of values. + /// + /// @param f Field3D of diagonal entries; must have parallel slices allocated. + /// @returns A diagonal PetscCellOperator. PetscCellOperator diagonal(const Field3D& f) const; + + /// @brief Build a diagonal forward-leg-space operator from a leg vector. + /// + /// @param v Forward leg vector of diagonal entries. + /// @returns A diagonal PetscForwardLegOperator. PetscForwardLegOperator diagonal(const PetscForwardLegVector& v) const; + + /// @brief Build a diagonal backward-leg-space operator from a leg vector. + /// + /// @param v Backward leg vector of diagonal entries. + /// @returns A diagonal PetscBackwardLegOperator. PetscBackwardLegOperator diagonal(const PetscBackwardLegVector& v) const; + /// @brief Collection of parallel differential operators assembled for the + /// Support Operator Method (SOM). + /// + /// All operators are constructed together in @ref getParallel so that intermediate + /// quantities (leg weights, interpolated dl, dV) are shared correctly. + /// + /// Naming convention: the subscript on Grad/Div/Restrict refers to which side of + /// the cell face the operator acts on, not the leg space it lives in. + /// - @c Grad_plus / @c Div_plus: forward-side half-step operators (in L-). + /// - @c Grad_minus / @c Div_minus: backward-side half-step operators (in L+). + /// - @c Restrict_minus = I+^T * W+: weighted restriction from L+ back to C, + /// paired with the minus-side gradient. + /// - @c Restrict_plus = I-^T * W-: weighted restriction from L- back to C, + /// paired with the plus-side gradient. struct Parallel { + /// Central parallel gradient: C ← C. PetscCellOperator Grad_par; + /// Support-operator parallel divergence paired with Grad_par: C ← C. PetscCellOperator Div_par; + /// Combined SOM parallel Laplacian (Div_par ∘ Grad_par): C ← C. PetscCellOperator Div_par_Grad_par; + /// Cell volume J * dx * dy * dz. Field3D dV; + /// Half-step gradient on the backward (y-1) side: L- ← C. PetscBackwardOperator Grad_minus; + /// Half-step gradient on the forward (y+1) side: L+ ← C. PetscForwardOperator Grad_plus; + /// SOM divergence paired with Grad_plus: C ← L+. PetscForwardToCellOperator Div_minus; + /// SOM divergence paired with Grad_minus: C ← L-. PetscBackwardToCellOperator Div_plus; + /// Raw backward interpolation operator B: L- ← C. PetscBackwardOperator Backward; + /// Raw forward interpolation operator F: L+ ← C. PetscForwardOperator Forward; + /// Injection operator into backward leg space I-: L- ← C. PetscBackwardOperator Inject_minus; + /// Injection operator into forward leg space I+: L+ ← C. PetscForwardOperator Inject_plus; + /// Average interpolation to backward leg positions: L- ← C. PetscBackwardOperator Interp_minus; + /// Average interpolation to forward leg positions: L+ ← C. PetscForwardOperator Interp_plus; + /// Weighted restriction from L+ to C: I+^T * W+. PetscForwardToCellOperator Restrict_minus; + /// Weighted restriction from L- to C: I-^T * W-. PetscBackwardToCellOperator Restrict_plus; + /// @brief Compute the SOM parallel diffusion Div_par(K * Grad_par(f)). + /// + /// Evaluates the half-step fluxes separately on each side, multiplies by the + /// interpolated diffusion coefficient K, and averages the two divergence + /// contributions. This is the primary user-facing SOM operator. + /// + /// @param K Diffusion coefficient field; must have parallel slices allocated. + /// @param f Field to differentiate; must have parallel slices allocated. + /// @returns Field3D containing Div_par(K * Grad_par(f)) at interior cells. Field3D Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const; }; + /// @brief Assemble and return the full set of parallel SOM operators. + /// + /// Constructs all primitive operators (Forward, Backward, Inject_plus, + /// Inject_minus, leg weights), derives the interpolation, gradient, restriction, + /// and divergence operators, and returns them packaged in a @ref Parallel struct. + /// Intermediate Field3D quantities (dl, dV) are computed from the mesh + /// coordinates and boundary conditions are applied before use. + /// + /// @returns A @ref Parallel struct containing all assembled operators. Parallel getParallel() const; private: + /// @brief Read a Field3D from the mesh, throwing if the field is absent. + /// + /// @param mesh Non-owning mesh pointer. + /// @param name Field name as stored in the mesh file. + /// @returns The requested Field3D. + /// @throws BoutException if the field is not found. static Field3D meshGetField3D(Mesh* mesh, const std::string& name); + /// @brief Read an Array from the mesh, throwing if absent. + /// + /// @tparam T Element type of the array. + /// @param name Array name as stored in the mesh file. + /// @returns The requested Array. + /// @throws BoutException if the array is not found. template Array meshGetArray(const std::string& name) const { Array result; @@ -502,13 +957,39 @@ private: return result; } + /// @brief Read a scalar integer from the mesh, throwing if absent. + /// + /// @param name Scalar name as stored in the mesh file. + /// @returns The requested integer value. + /// @throws BoutException if the scalar is not found. int meshGetInt(const std::string& name) const; + /// @brief Build a leg-space injection operator for the given direction. + /// + /// Iterates owned interior cells and inserts a unit entry for each interior or + /// boundary leg row, mapping to the corresponding cell-space column. Column + /// indices come from @p interior_leg_number and @p boundary_leg_number. + /// + /// @tparam LegTag @ref ForwardLegSpaceTag or @ref BackwardLegSpaceTag. + /// @param interior_leg_number Field3D of stored interior leg row numbers. + /// @param boundary_leg_number Field3D of stored boundary leg row numbers. + /// @param leg_mapping Shared mapping for the leg space. + /// @returns The assembled injection operator LegTag ← CellSpaceTag. template PetscOperator buildInjection(const Field3D& interior_leg_number, const Field3D& boundary_leg_number, std::shared_ptr leg_mapping) const; + /// @brief Load a leg-weight vector from a named mesh array. + /// + /// Reads the flat array @p name from the mesh, then scatters values into the + /// local PETSc vector using the mapping's stored index list. Stored indices that + /// are out of range are assigned weight 0. + /// + /// @tparam LegVector @ref PetscForwardLegVector or @ref PetscBackwardLegVector. + /// @param name Array name in the mesh file (e.g. "forward_leg_weights"). + /// @param mapping Shared mapping for the leg space. + /// @returns A fully populated leg-weight vector. template LegVector loadLegWeights(const std::string& name, std::shared_ptr mapping) const { @@ -526,14 +1007,14 @@ private: return out; } - Mesh* mesh{nullptr}; - PetscCellMappingPtr cell_mapping; - PetscLegMappingPtr forward_leg_mapping; - PetscLegMappingPtr backward_leg_mapping; - Field3D forward_leg_interior_number; - Field3D forward_leg_boundary_number; - Field3D backward_leg_interior_number; - Field3D backward_leg_boundary_number; + Mesh* mesh{nullptr}; ///< Non-owning pointer to the BOUT++ mesh. + PetscCellMappingPtr cell_mapping; ///< Cell-space index mapping. + PetscLegMappingPtr forward_leg_mapping; ///< Forward leg-space index mapping. + PetscLegMappingPtr backward_leg_mapping; ///< Backward leg-space index mapping. + Field3D forward_leg_interior_number; ///< Stored interior leg numbers for L+. + Field3D forward_leg_boundary_number; ///< Stored boundary leg numbers for L+. + Field3D backward_leg_interior_number; ///< Stored interior leg numbers for L-. + Field3D backward_leg_boundary_number; ///< Stored boundary leg numbers for L-. }; #else