From b48c40b1dfb1dbd269ee61cc3fdb76ebdbf1af86 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 14:15:58 -0700 Subject: [PATCH 01/11] PetscOperators: Minor docstring fixes --- include/bout/petsc_operators.hxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 0d9c6c6912..8c36f461d0 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -8,7 +8,6 @@ #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" @@ -865,8 +864,8 @@ public: /// /// 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 Grad_plus / @c Div_plus: forward-side half-step operators (C to L+). + /// - @c Grad_minus / @c Div_minus: backward-side half-step operators (C to 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, @@ -914,7 +913,7 @@ public: /// /// 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. + /// contributions. /// /// @param K Diffusion coefficient field; must have parallel slices allocated. /// @param f Field to differentiate; must have parallel slices allocated. From 325019aeacbae942d680dfb7671a98cfe72eff3d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 14:29:21 -0700 Subject: [PATCH 02/11] PetscIndexMapping::makeEvolvingIS unit tests This function should create an index set of the evolving cells. It will be used to extract the part of the PetscOperator that represents coupling between evolving cells. Implementation to follow. --- include/bout/petsc_operators.hxx | 8 + tests/unit/CMakeLists.txt | 1 + .../test_petsc_operators_make_evolving_is.cxx | 487 ++++++++++++++++++ 3 files changed, 496 insertions(+) create mode 100644 tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 8c36f461d0..3a452d87ac 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -107,6 +107,14 @@ public: /// PETSc column indices via a post-multiplication. Mat getPetscToStored() const { return mat_petsc_to_stored; } + /// Dummy implementation + IS makeEvolvingIS() const { + IS is; + const int idx[1] = {0}; + ISCreateGeneral(MPI_COMM_WORLD, 1, &idx[0], PETSC_USE_POINTER, &is); + return is; + } + protected: /// @brief BOUT++ PETSc library handle; ensures PETSc is initialised. PetscLib lib; diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index d3d690417b..86360b1d35 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -98,6 +98,7 @@ set(serial_tests_source ./mesh/test_mesh.cxx ./mesh/test_paralleltransform.cxx ./mesh/test_petsc_operators.cxx + ./mesh/test_petsc_operators_make_evolving_is.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_solver.cxx diff --git a/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx new file mode 100644 index 0000000000..c28ef78937 --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx @@ -0,0 +1,487 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#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" +#include "bout/region.hxx" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include +#include +#include + +#include + +using bout::globals::mesh; + +// ============================================================================ +// Helper functions shared across unit tests +// ============================================================================ + +namespace { + +/// Build a PetscCellMapping whose evolving region has exactly the cells listed +/// in cell_number with non-negative values, plus whatever boundary cells are in +/// forward_cell_number / backward_cell_number. total_cells must equal the +/// number of distinct non-negative entries across all three fields. +std::shared_ptr +makeMappingFromFields(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells) { + return std::make_shared(cell_number, forward_cell_number, + backward_cell_number, total_cells); +} + +/// Collect the global PETSc indices that mapOwnedInteriorCells visits, in order. +std::vector collectOwnedInteriorIndices(const PetscCellMapping& mapping) { + std::vector indices; + mapping.mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D& /*i*/, int /*stored*/) { indices.push_back(row); }); + return indices; +} + +/// Unwrap an IS into a sorted std::vector. The IS is not destroyed. +std::vector isToSortedVector(IS is) { + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + std::vector v(idx, idx + n); + ISRestoreIndices(is, &idx); + std::sort(v.begin(), v.end()); + return v; +} +} // Namespace + +// ============================================================================ +// Fixture +// ============================================================================ + +// FakeMeshFixture provides a mesh with: +// LocalNx = 5 (xstart=1, xend=3, so mxg=1, 3 interior x-points) +// LocalNy = 4 (ystart=1, yend=2, so mgy=1, 2 interior y-points) +// LocalNz = 2 +// Interior cells: x in [1,3], y in [1,2], z in [0,1] => 3*2*2 = 12 evolving cells. +using PetscEvolvingISTest = FakeMeshFixture; + +// ============================================================================ +// Helper: build the standard mapping used by most tests in this suite. +// +// Stored cell numbers are assigned to interior cells only (indices 0..11). +// No forward or backward boundary cells. +// ============================================================================ +namespace { +std::shared_ptr buildStandardMapping(Mesh* mesh) { + Field3D cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + return makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, + next); +} +} // Namespace + +// ============================================================================ +// IS_size_equals_n_evolving +// +// The IS produced by makeEvolvingIS() must contain exactly one index per +// evolving interior cell. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_size_equals_n_evolving) { + auto mapping = buildStandardMapping(mesh); + + // FakeMesh: 3 x-interior * 2 y-interior * 2 z = 12 evolving cells + const PetscInt expected_n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + + IS is = mapping->makeEvolvingIS(); + + PetscInt local_size; + ISGetLocalSize(is, &local_size); + + EXPECT_EQ(expected_n_evolving, local_size); + + ISDestroy(&is); +} + +// ============================================================================ +// IS_indices_match_mapOwnedInteriorCells_order +// +// The IS indices must appear in exactly the same order as the global PETSc +// row indices produced by mapOwnedInteriorCells. Both traversals must visit +// the same cells in the same sequence. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_match_mapOwnedInteriorCells_order) { + auto mapping = buildStandardMapping(mesh); + + const std::vector from_map = collectOwnedInteriorIndices(*mapping); + + IS is = mapping->makeEvolvingIS(); + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + const std::vector from_is(idx, idx + n); + ISRestoreIndices(is, &idx); + ISDestroy(&is); + + ASSERT_EQ(from_map.size(), from_is.size()); + for (std::size_t i = 0; i < from_map.size(); ++i) { + EXPECT_EQ(from_map[i], from_is[i]) << "Index mismatch at position " << i; + } +} + +// ============================================================================ +// IS_indices_are_in_global_petsc_range +// +// On a single MPI rank the global PETSc row range is [rowStart(), rowEnd()). +// Every index in the IS must lie in that range. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_are_in_global_petsc_range) { + auto mapping = buildStandardMapping(mesh); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + const PetscInt row_start = mapping->rowStart(); + const PetscInt row_end = mapping->rowEnd(); + + for (PetscInt idx : v) { + EXPECT_GE(idx, row_start) << "Index " << idx << " is below rowStart()"; + EXPECT_LT(idx, row_end) << "Index " << idx << " is at or above rowEnd()"; + } +} + +// ============================================================================ +// IS_excludes_xin_boundary_cells +// +// Cells in the inner-X boundary region have stored indices assigned by the +// mapping but must NOT appear in the evolving IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_xin_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + // Assign evolving cells + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // Assign one inner-X boundary cell a stored index + const int xin_stored = next++; + cell_number(0, mesh->ystart, 0) = xin_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt xin_petsc = mapping->storedToPetsc(xin_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), xin_petsc)) + << "Inner-X boundary cell (PETSc index " << xin_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_xout_boundary_cells +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_xout_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // Outer-X boundary: x == LocalNx - 1 + const int xout_stored = next++; + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = xout_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt xout_petsc = mapping->storedToPetsc(xout_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), xout_petsc)) + << "Outer-X boundary cell (PETSc index " << xout_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_forward_boundary_cells +// +// Virtual forward-parallel boundary cells (yup) live in forward_cell_number, +// not cell_number. They must not appear in the evolving IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_forward_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // One virtual forward boundary cell + const int fwd_stored = next++; + forward_cell_number(mesh->xstart, mesh->yend, 0) = fwd_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt fwd_petsc = mapping->storedToPetsc(fwd_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), fwd_petsc)) + << "Forward boundary cell (PETSc index " << fwd_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_backward_boundary_cells +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_backward_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + const int bwd_stored = next++; + backward_cell_number(mesh->xstart, mesh->ystart, 0) = bwd_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt bwd_petsc = mapping->storedToPetsc(bwd_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), bwd_petsc)) + << "Backward boundary cell (PETSc index " << bwd_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_is_empty_when_no_evolving_cells +// +// If cell_number is entirely -1, there are no evolving cells and the IS must +// have size zero. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_is_empty_when_no_evolving_cells) { + const Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + // total_cells = 0: no cells at all + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, 0); + + IS is = mapping->makeEvolvingIS(); + PetscInt local_size; + ISGetLocalSize(is, &local_size); + ISDestroy(&is); + + EXPECT_EQ(0, local_size); +} + +// ============================================================================ +// IS_covers_all_evolving_stored_indices +// +// The set of stored indices reachable from the IS (via storedToPetsc inverse) +// must exactly equal the set of stored indices of interior cells as visited by +// mapOwnedInteriorCells. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_covers_all_evolving_stored_indices) { + auto mapping = buildStandardMapping(mesh); + + // Collect stored indices from mapOwnedInteriorCells + std::vector stored_from_map; + mapping->mapOwnedInteriorCells([&](PetscInt /*row*/, const Ind3D& /*i*/, int stored) { + stored_from_map.push_back(stored); + }); + std::sort(stored_from_map.begin(), stored_from_map.end()); + + // Collect PETSc indices from IS, then map back to stored via localStoredIndices + IS is = mapping->makeEvolvingIS(); + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + + // localStoredIndices[i] is the stored index for PETSc row rowStart()+i + const auto& local_stored = mapping->localStoredIndices(); + const PetscInt rstart = mapping->rowStart(); + std::vector stored_from_is; + for (PetscInt k = 0; k < n; ++k) { + const PetscInt local_row = idx[k] - rstart; + ASSERT_GE(local_row, 0); + ASSERT_LT(local_row, static_cast(local_stored.size())); + stored_from_is.push_back(local_stored[local_row]); + } + ISRestoreIndices(is, &idx); + ISDestroy(&is); + + std::sort(stored_from_is.begin(), stored_from_is.end()); + + EXPECT_EQ(stored_from_map, stored_from_is); +} + +// ============================================================================ +// IS_indices_are_unique +// +// No PETSc index should appear more than once in the IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_are_unique) { + auto mapping = buildStandardMapping(mesh); + + IS is = mapping->makeEvolvingIS(); + const auto sorted = isToSortedVector(is); + ISDestroy(&is); + + for (std::size_t i = 1; i < sorted.size(); ++i) { + EXPECT_LT(sorted[i - 1], sorted[i]) + << "Duplicate index " << sorted[i] << " found at positions " << i - 1 << " and " + << i; + } +} + +// ============================================================================ +// IS_with_mixed_boundary_cells_has_correct_size +// +// When forward, backward, and X-boundary cells are all present alongside +// evolving cells, the IS size equals only the evolving cell count. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_with_mixed_boundary_cells_has_correct_size) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + // Evolving interior cells + const int n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // One of each boundary type + cell_number(0, mesh->ystart, 0) = next++; // xin + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = next++; // xout + forward_cell_number(mesh->xstart, mesh->yend, 0) = next++; // yup + backward_cell_number(mesh->xstart, mesh->ystart, 0) = next++; // ydown + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + IS is = mapping->makeEvolvingIS(); + PetscInt local_size; + ISGetLocalSize(is, &local_size); + ISDestroy(&is); + + EXPECT_EQ(n_evolving, local_size); +} + +// ============================================================================ +// IS_can_be_used_to_create_submatrix +// +// Regression guard: the IS must be acceptable to MatCreateSubMatrix without +// error. Uses a simple identity-pattern matrix over the full cell space. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_can_be_used_to_create_submatrix) { + auto mapping = buildStandardMapping(mesh); + + // Build a diagonal matrix over the full cell space + const PetscInt n = mapping->globalSize(); + const PetscInt nlocal = mapping->localSize(); + Mat A; + MatCreate(BoutComm::get(), &A); + MatSetSizes(A, nlocal, nlocal, n, n); + MatSetType(A, MATMPIAIJ); + MatMPIAIJSetPreallocation(A, 1, nullptr, 0, nullptr); + const PetscInt rstart = mapping->rowStart(); + for (PetscInt i = 0; i < nlocal; ++i) { + const PetscInt global = rstart + i; + const PetscScalar val = 1.0; + MatSetValues(A, 1, &global, 1, &global, &val, INSERT_VALUES); + } + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + + IS is = mapping->makeEvolvingIS(); + + Mat sub = nullptr; + // This must not throw or return a PETSc error + const PetscErrorCode ierr = MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &sub); + + EXPECT_EQ(PETSC_SUCCESS, ierr); + + // The submatrix size should equal n_evolving × n_evolving + if (sub != nullptr) { + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + const PetscInt n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + MatDestroy(&sub); + } + + ISDestroy(&is); + MatDestroy(&A); +} + +#endif // BOUT_HAS_PETSC From b433faf5b9e5494ae631ec0fd2901c0dbeaea532 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 14:45:55 -0700 Subject: [PATCH 03/11] PetscCellMapping::makeEvolvingIS implementation Creates an index set of global PETSc rows that correspond to evolving cells. Passes unit tests. --- include/bout/petsc_operators.hxx | 24 ++++++++++++++++-------- src/mesh/petsc_operators.cxx | 15 +++++++++++++++ 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 3a452d87ac..98e706234b 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -107,14 +107,6 @@ public: /// PETSc column indices via a post-multiplication. Mat getPetscToStored() const { return mat_petsc_to_stored; } - /// Dummy implementation - IS makeEvolvingIS() const { - IS is; - const int idx[1] = {0}; - ISCreateGeneral(MPI_COMM_WORLD, 1, &idx[0], PETSC_USE_POINTER, &is); - return is; - } - protected: /// @brief BOUT++ PETSc library handle; ensures PETSc is initialised. PetscLib lib; @@ -282,6 +274,22 @@ public: } } + /// @brief Create a PETSc IS containing the global PETSc indices of all + /// locally owned evolving interior cells, in the order that + /// mapOwnedInteriorCells visits them. + /// + /// The returned IS selects the evolving subset of the full cell space C, + /// excluding inner/outer X-boundary cells and forward/backward parallel + /// boundary virtual cells. It is the correct IS to pass to + /// MatCreateSubMatrix when restricting a PetscCellOperator to the degrees + /// of freedom that the SNES solver actually integrates. + /// + /// The caller owns the returned IS and must call ISDestroy when finished. + /// + /// @returns A PETSc IS of local size equal to the number of locally owned + /// evolving cells, holding global PETSc row indices. + IS makeEvolvingIS() const; + private: Field3D cell_number; ///< Stored cell numbers for interior/X-boundary cells. Field3D diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 0fe564bcd3..f72a433220 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -127,6 +127,21 @@ PetscCellMapping::PetscCellMapping(const Field3D& cell_number, local_indices); } +IS PetscCellMapping::makeEvolvingIS() const { + // Collect global PETSc indices in mapOwnedInteriorCells order. + // Reserve the known count up front to avoid reallocation. + std::vector indices; + indices.reserve(static_cast(evolving_region.size())); + + mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D& /*i*/, int /*stored*/) { indices.push_back(row); }); + + IS is; + BOUT_DO_PETSC(ISCreateGeneral(BoutComm::get(), static_cast(indices.size()), + indices.data(), PETSC_COPY_VALUES, &is)); + return is; +} + 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()), From 65185d019dfe4ff8b692baf0859b4f68ba8e0f80 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 17:09:45 -0700 Subject: [PATCH 04/11] PetscCellMapping::extractEvolvingSubmatrix unit tests Method that will extract the part of a PetscCellOperator that couples evolving cells to other evolving cells. Dummy implementation with failing tests. --- include/bout/petsc_operators.hxx | 19 + src/mesh/petsc_operators.cxx | 8 + tests/unit/CMakeLists.txt | 1 + ...c_operators_extract_evolving_submatrix.cxx | 488 ++++++++++++++++++ .../test_petsc_operators_make_evolving_is.cxx | 3 - 5 files changed, 516 insertions(+), 3 deletions(-) create mode 100644 tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 98e706234b..14d071fdc6 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -44,6 +44,10 @@ struct ForwardLegSpaceTag {}; /// leg space. See @ref ForwardLegSpaceTag and @ref CellSpaceTag. struct BackwardLegSpaceTag {}; +// Forward declare +template +class PetscOperator; + /// @brief Bidirectional index mapping between mesh-file stored numbering and PETSc /// row ownership. /// @@ -290,6 +294,21 @@ public: /// evolving cells, holding global PETSc row indices. IS makeEvolvingIS() const; + /// @brief Extract the evolving-cell submatrix from a cell-to-cell operator. + /// + /// Restricts @p op to the rows and columns that correspond to evolving + /// interior cells, discarding any rows or columns that belong to inner/outer + /// X-boundary cells or forward/backward parallel boundary virtual cells. + /// + /// The returned Mat is an independent copy (MAT_INITIAL_MATRIX): subsequent + /// modifications to @p op do not affect it. The caller owns the returned + /// Mat and must call MatDestroy when finished. + /// + /// @param op A cell-to-cell operator whose row and column space is the full + /// cell space C managed by this mapping. + /// @returns A square Mat of global size n_evolving × n_evolving. + Mat extractEvolvingSubmatrix(const PetscOperator& op) const; + private: Field3D cell_number; ///< Stored cell numbers for interior/X-boundary cells. Field3D diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index f72a433220..b257113011 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -142,6 +142,14 @@ IS PetscCellMapping::makeEvolvingIS() const { return is; } +// Dummy implementation. Note: Result is owned by the caller. +Mat PetscCellMapping::extractEvolvingSubmatrix( + const PetscOperator& op) const { + Mat new_mat; + BOUT_DO_PETSC(MatDuplicate(op.raw(), MAT_COPY_VALUES, &new_mat)); + return new_mat; +} + 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()), diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 86360b1d35..6ea6f92a6b 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -99,6 +99,7 @@ set(serial_tests_source ./mesh/test_paralleltransform.cxx ./mesh/test_petsc_operators.cxx ./mesh/test_petsc_operators_make_evolving_is.cxx + ./mesh/test_petsc_operators_extract_evolving_submatrix.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_solver.cxx diff --git a/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx new file mode 100644 index 0000000000..ea67202863 --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx @@ -0,0 +1,488 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "bout/array.hxx" +#include "bout/bout_types.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/region.hxx" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include +#include +#include +#include + +#include + +using bout::globals::mesh; + +// ============================================================================ +// PetscCellMapping::extractEvolvingSubmatrix() +// +// Concept: given a PetscCellOperator with a known nonzero pattern, the +// submatrix returned by extractEvolvingSubmatrix() contains exactly the rows +// and columns that correspond to evolving interior cells, with values +// preserved, and no entries from boundary cells. +// +// The fixture is FakeMeshFixture (LocalNx=5, LocalNy=4, LocalNz=2, mxg=1). +// Interior cells: x in [1,3], y in [1,2], z in [0,1] => 12 evolving cells. +// ============================================================================ + +// ── Helpers ───────────────────────────────────────────────────────────────── + +namespace { + +// Build a mapping where evolving cells have stored indices [0, n_evolving), +// with optional extra boundary cells appended after. +struct MappingWithCells { + std::shared_ptr mapping; + int n_evolving; + // Global PETSc index of each evolving cell, in mapOwnedInteriorCells order + std::vector evolving_petsc; + // Stored indices of boundary cells (xin, xout, yup, ydown) + std::vector boundary_stored; +}; + +MappingWithCells buildMappingWithBoundaries(Mesh* mesh, bool add_xin = false, + bool add_xout = false, bool add_yup = false, + bool add_ydown = false) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + const int n_evolving = next; + + std::vector boundary_stored; + if (add_xin) { + cell_number(0, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + if (add_xout) { + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + if (add_yup) { + forward_cell_number(mesh->xstart, mesh->yend, 0) = next; + boundary_stored.push_back(next++); + } + if (add_ydown) { + backward_cell_number(mesh->xstart, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, next); + + std::vector evolving_petsc; + mapping->mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D&, int) { evolving_petsc.push_back(row); }); + + return {std::move(mapping), n_evolving, std::move(evolving_petsc), + std::move(boundary_stored)}; +} + +// Build a PetscCellOperator from explicit (row, col, value) triples over a +// given mapping. Rows and cols are stored indices (not PETSc global). +// The operator is assembled from a hand-built CSR representation. +PetscCellOperator +buildOperatorFromTriples(std::shared_ptr mapping, + const std::vector>& triples) { + + const int n = mapping->globalSize(); + // Build dense CSR: rows array size n+1 + std::vector row_counts(n, 0); + for (const auto& [r, c, v] : triples) { + row_counts[r]++; + } + // rows (CSR row pointer) + Array rows(n + 1); + rows[0] = 0; + for (int i = 0; i < n; ++i) { + rows[i + 1] = rows[i] + row_counts[i]; + } + // cols and weights + Array cols(static_cast(triples.size())); + Array weights(static_cast(triples.size())); + std::vector fill(n, 0); + for (const auto& [r, c, v] : triples) { + const int pos = rows[r] + fill[r]++; + cols[pos] = c; + weights[pos] = v; + } + + return PetscCellOperator(mapping, mapping, rows, cols, weights); +} + +// Return a flat vector of (global_row, global_col, value) for every nonzero +// in a matrix, using MatGetRow. Rows are iterated over the local ownership +// range only. +struct Nonzero { + PetscInt row, col; + PetscScalar val; +}; + +std::vector matNonzeros(Mat A) { + PetscInt rstart, rend; + MatGetOwnershipRange(A, &rstart, &rend); + std::vector nzs; + for (PetscInt row = rstart; row < rend; ++row) { + PetscInt ncols; + const PetscInt* col_idx; + const PetscScalar* vals; + MatGetRow(A, row, &ncols, &col_idx, &vals); + for (PetscInt k = 0; k < ncols; ++k) { + nzs.push_back({row, col_idx[k], vals[k]}); + } + MatRestoreRow(A, row, &ncols, &col_idx, &vals); + } + return nzs; +} + +} // namespace + +// ── Fixture ───────────────────────────────────────────────────────────────── + +using PetscExtractEvolvingTest = FakeMeshFixture; + +// ============================================================================ +// submatrix_global_size_equals_n_evolving_squared +// +// MatGetSize on the result must report n_evolving × n_evolving. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_global_size_equals_n_evolving_squared) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true); + + // Diagonal operator over the full cell space (safe: every cell maps to itself) + const auto& local = mapping->localStoredIndices(); + std::vector> triples; + for (int s : local) { + if (s >= 0) { + triples.emplace_back(s, s, 1.0); + } + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_local_size_equals_n_evolving_locally +// +// MatGetLocalSize must report the local evolving count on each rank. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_local_size_equals_n_evolving_locally) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, true, true, true, true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt local_rows, local_cols; + MatGetLocalSize(sub, &local_rows, &local_cols); + + // On a single rank local == global + EXPECT_EQ(n_evolving, local_rows); + EXPECT_EQ(n_evolving, local_cols); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_values_preserved_for_evolving_block +// +// For entries (r, c) where both r and c are evolving cells, the value in the +// submatrix must equal the value in the original operator. +// +// Strategy: build an operator with known values on evolving-to-evolving +// entries, extract the submatrix, and check each value via MatGetValues. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_values_preserved_for_evolving_block) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true, + /*yup=*/true, /*ydown=*/true); + + // Insert entries with distinct values between pairs of evolving cells. + // Use stored indices directly (0..n_evolving-1). + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + // Diagonal entry + triples.emplace_back(i, i, static_cast(10 + i)); + // One off-diagonal entry per row (wrap around) + const int j = (i + 1) % n_evolving; + triples.emplace_back(i, j, static_cast(100 + i)); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + const auto nzs = matNonzeros(sub); + + // Build expected map: submatrix uses 0-based indices within the evolving + // block. The IS preserves mapOwnedInteriorCells order, which matches + // stored index order (0..n_evolving-1) on a single rank. + // Sub row/col i corresponds to evolving stored index i. + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + + for (const auto& nz : nzs) { + // nz.row and nz.col are global indices in the submatrix space. + // On a single rank these equal the stored indices of the evolving cells. + const int stored_row = static_cast(nz.row - sub_rstart); + const int stored_col = static_cast(nz.col); + + // Find the expected value from our triples + BoutReal expected = 0.0; + for (const auto& [r, c, v] : triples) { + if (r == stored_row && c == stored_col) { + expected = v; + break; + } + } + EXPECT_DOUBLE_EQ(expected, nz.val) + << "Wrong value at submatrix (" << nz.row << ", " << nz.col << ")"; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_nonzero_count_matches_evolving_entries +// +// The total nonzero count in the submatrix must equal the number of entries +// in the original operator whose row AND column are both evolving cells. +// Entries coupling to boundary cells must be absent. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_nonzero_count_matches_evolving_entries) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/false, + /*yup=*/true, /*ydown=*/false); + + std::vector> triples; + // Entries entirely within the evolving block + int expected_in_sub = 0; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + ++expected_in_sub; + } + // Entries from an evolving row into a boundary column — must be excluded + for (int bnd : boundary_stored) { + triples.emplace_back(0, bnd, 99.0); // row 0 (evolving) -> boundary col + } + // Entries from a boundary row into an evolving column — must be excluded + for (int bnd : boundary_stored) { + triples.emplace_back(bnd, 0, 88.0); // boundary row -> col 0 (evolving) + } + + const auto op = buildOperatorFromTriples(mapping, triples); + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt nnz; + MatInfo info; + MatGetInfo(sub, MAT_GLOBAL_SUM, &info); + nnz = static_cast(info.nz_used); + + EXPECT_EQ(expected_in_sub, nnz); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_excludes_boundary_columns +// +// No column index in the submatrix should correspond to a boundary cell's +// PETSc index in the original cell space. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_excludes_boundary_columns) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true, + /*yup=*/true, /*ydown=*/true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + // Add evolving-row → boundary-column coupling + for (int bnd : boundary_stored) { + triples.emplace_back(0, bnd, 7.0); + } + + const auto op = buildOperatorFromTriples(mapping, triples); + + // Collect the PETSc indices of the boundary cells so we can check against them. + // The submatrix uses a re-indexed column space [0, n_evolving), so boundary + // columns from the full space simply will not appear. + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + + // The boundary-coupling entries were given the sentinel value 7.0. + // If any survive into the submatrix, the extraction has failed to exclude them. + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_NE(7.0, nz.val) + << "Boundary-column entry (sentinel value 7.0) survived into the " + "submatrix at (" + << nz.row << ", " << nz.col << ")"; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_excludes_boundary_rows +// +// The submatrix must have no rows for boundary cells. Verified by checking +// that MatGetLocalSize equals n_evolving (not n_evolving + n_boundary). +// Also checked by confirming no row index outside [sub_rstart, sub_rend) +// appears in the nonzeros. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_excludes_boundary_rows) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/false, + /*yup=*/false, /*ydown=*/true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + // boundary-row → evolving-column entries + for (int bnd : boundary_stored) { + triples.emplace_back(bnd, 0, 5.0); + } + + const auto op = buildOperatorFromTriples(mapping, triples); + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt local_rows, local_cols; + MatGetLocalSize(sub, &local_rows, &local_cols); + EXPECT_EQ(n_evolving, local_rows); + + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_GE(nz.row, sub_rstart); + EXPECT_LT(nz.row, sub_rend); + } + + MatDestroy(&sub); +} + +// ============================================================================ +// zero_operator_gives_zero_submatrix +// +// A zero operator (no nonzeros) must produce an assembled zero submatrix of +// the correct size. No crash, no uninitialized memory. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, zero_operator_gives_zero_submatrix) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + // Empty triples: zero operator + const auto op = buildOperatorFromTriples(mapping, {}); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + + MatInfo info; + MatGetInfo(sub, MAT_GLOBAL_SUM, &info); + EXPECT_EQ(0, static_cast(info.nz_used)); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_of_identity_is_identity +// +// An identity operator over the evolving cells only (no boundary entries) +// should produce an identity submatrix: diagonal entries of 1.0, nothing else. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_of_identity_is_identity) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + const auto nzs = matNonzeros(sub); + + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + + ASSERT_EQ(n_evolving, static_cast(nzs.size())); + for (const auto& nz : nzs) { + EXPECT_EQ(nz.row, nz.col) << "Off-diagonal entry in identity submatrix"; + EXPECT_DOUBLE_EQ(1.0, nz.val); + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_is_independent_of_original_operator +// +// Modifying the original operator's underlying matrix after extraction must +// not affect the already-returned submatrix. MatCreateSubMatrix with +// MAT_INITIAL_MATRIX produces an independent copy. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_is_independent_of_original_operator) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 2.0); + } + auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + // Scale the original operator's matrix in place + MatScale(op.raw(), 0.0); + + // The submatrix must still hold the original values + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_DOUBLE_EQ(2.0, nz.val) + << "Submatrix was affected by modification of the source operator"; + } + + MatDestroy(&sub); +} + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx index c28ef78937..9dbaed584e 100644 --- a/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx +++ b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx @@ -4,9 +4,6 @@ #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" #include "bout/region.hxx" From 4920c8616e8c0fefec972d0002235bfd9cda77a4 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 17:14:17 -0700 Subject: [PATCH 05/11] PetscCellMapping::extractEvolvingSubmatrix implementation Passes unit tests --- src/mesh/petsc_operators.cxx | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index b257113011..9700001d23 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -142,12 +142,16 @@ IS PetscCellMapping::makeEvolvingIS() const { return is; } -// Dummy implementation. Note: Result is owned by the caller. Mat PetscCellMapping::extractEvolvingSubmatrix( const PetscOperator& op) const { - Mat new_mat; - BOUT_DO_PETSC(MatDuplicate(op.raw(), MAT_COPY_VALUES, &new_mat)); - return new_mat; + IS is = makeEvolvingIS(); + + Mat sub; + BOUT_DO_PETSC(MatCreateSubMatrix(op.raw(), is, is, MAT_INITIAL_MATRIX, &sub)); + + BOUT_DO_PETSC(ISDestroy(&is)); + + return sub; } PetscLegMapping::PetscLegMapping(int total_legs, std::vector local_leg_indices) { From 9feb2932384e3ec5c62b1bd158339b1c1781ef7e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 22:01:53 -0700 Subject: [PATCH 06/11] test-petsc-ordering integrated test Checks the global ordering of cells for consistency as number of processors and nxpe is varied. --- tests/integrated/CMakeLists.txt | 1 + .../test-petsc-ordering/CMakeLists.txt | 7 + .../test-petsc-ordering/data/BOUT.inp | 14 ++ tests/integrated/test-petsc-ordering/runtest | 154 +++++++++++++ .../test_petsc_ordering.cxx | 215 ++++++++++++++++++ 5 files changed, 391 insertions(+) create mode 100644 tests/integrated/test-petsc-ordering/CMakeLists.txt create mode 100644 tests/integrated/test-petsc-ordering/data/BOUT.inp create mode 100755 tests/integrated/test-petsc-ordering/runtest create mode 100644 tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index 3696ef0372..809777fd58 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -34,6 +34,7 @@ add_subdirectory(test-naulin-laplace) add_subdirectory(test-options-netcdf) add_subdirectory(test-petsc_laplace) add_subdirectory(test-petsc_laplace_MAST-grid) +add_subdirectory(test-petsc-ordering) add_subdirectory(test-restart-io) add_subdirectory(test-restarting) add_subdirectory(test-slepc-solver) diff --git a/tests/integrated/test-petsc-ordering/CMakeLists.txt b/tests/integrated/test-petsc-ordering/CMakeLists.txt new file mode 100644 index 0000000000..4366c38657 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/CMakeLists.txt @@ -0,0 +1,7 @@ +bout_add_integrated_test( + test-petsc-ordering + SOURCES test_petsc_ordering.cxx + REQUIRES BOUT_HAS_PETSC + USE_RUNTEST USE_DATA_BOUT_INP + PROCESSORS 4 +) diff --git a/tests/integrated/test-petsc-ordering/data/BOUT.inp b/tests/integrated/test-petsc-ordering/data/BOUT.inp new file mode 100644 index 0000000000..9b51b9f296 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/data/BOUT.inp @@ -0,0 +1,14 @@ +# Input file for test_petsc_ordering integrated test. +# Uses a small grid large enough to exercise MPI decomposition on 4 ranks +# (nx=10 gives 2 interior x-points per rank when NXPE=4; ny and nz similarly). + +[mesh] +nx = 10 # 8 interior + 2 guard (mxg=1 each side) +ny = 8 +nz = 4 + +ixseps1 = nx +ixseps2 = nx + +[output] +enabled = false # Suppress data file output; we only need stdout diff --git a/tests/integrated/test-petsc-ordering/runtest b/tests/integrated/test-petsc-ordering/runtest new file mode 100755 index 0000000000..26540d5338 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/runtest @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +# requires: petsc +# cores: 4 +""" +Integrated test: Ordering equivalence between PetscCellMapping +and the SNES solver's globalIndex traversal. + +Runs test_petsc_ordering with 1, 2, and 4 MPI ranks and checks that: + - The executable exits with status 0. + - The output contains "ordering_check=PASS". + - No "MISMATCH" or "SHIFT_NOT_CONSTANT" lines appear. +""" + +import argparse +import re +import subprocess +import sys + +from boututils.run_wrapper import build_and_log + + +def parse_args(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--executable", + default="./test_petsc_ordering", + help="Path to the test executable", + ) + parser.add_argument( + "--mpirun", default="mpirun", help="MPI launcher (mpirun, srun, ...)" + ) + parser.add_argument( + "--nprocs", + type=int, + nargs="+", + default=[1, 2, 4], + help="List of MPI rank counts to test", + ) + parser.add_argument( + "--timeout", type=int, default=120, help="Per-run timeout in seconds" + ) + return parser.parse_args() + + +def run_case(mpirun, executable, nproc, nxpe, timeout): + """Launch one run and return (stdout+stderr, returncode).""" + cmd = [mpirun, "-n", str(nproc), executable, f"nxpe={nxpe}"] + try: + result = subprocess.run( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + timeout=timeout, + text=True, + ) + return result.stdout, result.returncode + except subprocess.TimeoutExpired: + return f"TIMED OUT after {timeout}s\n", -1 + except FileNotFoundError as e: + return f"Could not launch: {e}\n", -1 + + +def check_output(output, nproc): + """ + Inspect the combined output for pass/fail markers. + Returns a list of failure strings (empty means pass). + """ + failures = [] + + # Must contain the summary pass marker + if "ordering_check=PASS" not in output: + if "ordering_check=FAIL" in output: + failures.append("ordering_check=FAIL found in output") + else: + failures.append("ordering_check marker not found in output") + + # Must not contain any mismatch lines + mismatch_lines = [ + line for line in output.splitlines() if line.startswith("MISMATCH") + ] + if mismatch_lines: + failures.append( + f"{len(mismatch_lines)} MISMATCH line(s) found:\n" + + "\n".join(f" {line}" for line in mismatch_lines) + ) + + # Must not contain any shift-not-constant lines + shift_lines = [ + line for line in output.splitlines() if line.startswith("SHIFT_NOT_CONSTANT") + ] + if shift_lines: + failures.append( + f"{len(shift_lines)} SHIFT_NOT_CONSTANT line(s) found:\n" + + "\n".join(f" {line}" for line in shift_lines) + ) + + # Extract and report the numeric summary for informational purposes + for marker in ("total_mismatches", "total_shift_failures"): + m = re.search(rf"{marker}=(\d+)", output) + if m: + count = int(m.group(1)) + if count != 0: + failures.append(f"{marker}={count} (expected 0)") + + return failures + + +def main(): + args = parse_args() + + build_and_log("PETSc ordering test") + + overall_pass = True + + for nproc in args.nprocs: + for nxpe in [1, 2, 4]: + if nxpe > nproc: + break + print(f"\n{'=' * 60}") + print(f"Running with {nproc} MPI rank(s), nxpe={nxpe}") + print(f"{'=' * 60}") + + output, returncode = run_case( + args.mpirun, args.executable, nproc, nxpe, args.timeout + ) + + case_failures = [] + + if returncode != 0: + # Note: MPI task can exit non-zero for reasons not connected to test + print(f"Warning: Non-zero exit code: {returncode}") + + case_failures.extend(check_output(output, nproc)) + + if case_failures: + print(output) # Only print output on failure + overall_pass = False + print(f"FAIL (nproc={nproc}, nxpe={nxpe}):") + for f in case_failures: + print(f" - {f}") + else: + print(f"PASS (nproc={nproc}, nxpe={nxpe})") + + print(f"\n{'=' * 60}") + if overall_pass: + print("ALL CASES PASSED") + return 0 + else: + print("ONE OR MORE CASES FAILED") + return 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx new file mode 100644 index 0000000000..c3962a1b3f --- /dev/null +++ b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx @@ -0,0 +1,215 @@ +/// Integrated test for Unit 3: ordering equivalence between PetscCellMapping +/// and the SNES solver's globalIndex traversal. +/// +/// For every locally owned evolving interior cell this program checks: +/// +/// 1. The local offset into mapOwnedInteriorCells equals the local SNES +/// index from globalIndex (ordering equivalence). +/// +/// 2. The shift (rowStart() - Istart) is the same constant for every cell +/// on this rank (scalar-shift property). +/// +/// Results are written to the BOUT++ output so the runtest script can parse +/// them. The program exits with a non-zero status if any check fails. + +#include + +#if !BOUT_HAS_PETSC +// Without PETSc there is nothing to test. Exit cleanly so CTest skips rather +// than fails. +#include +int main() { return EXIT_SUCCESS; } +#else + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +int main(int argc, char** argv) { + // Initialise BOUT++ (reads BOUT.inp, sets up mesh, MPI, PETSc, etc.) + BoutInitialise(argc, argv); + + Mesh* mesh = bout::globals::mesh; + + // ── Build cell_number field ──────────────────────────────────────────────── + // Assign stored indices to interior cells exactly as the Python weights + // module does: consecutive integers over the evolving region. + Field3D cell_number{-1.0, mesh}; + { + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // Communicate so guard cells are filled (needed by some mesh operations) + mesh->communicate(cell_number); + } + + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + // Count total evolving cells across all ranks + int n_local_evolving = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { + ++n_local_evolving; + } + } + int n_global_evolving = 0; + MPI_Allreduce(&n_local_evolving, &n_global_evolving, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + // ── Build PetscCellMapping ───────────────────────────────────────────────── + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, n_global_evolving); + + // ── Build globalIndex field (replicates FDJinitialise logic) ────────────── + // globalIndex(0) gives 0-based local indices over RGN_NOBNDRY. + // After shifting by Istart these become the SNES global indices. + Field3D snes_index{-1.0, mesh}; + { + int local_idx = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { + snes_index[i] = local_idx++; + } + } + } + + // Determine SNES global offset for this rank by constructing a dummy PETSc + // matrix of size n_global_evolving and reading its ownership range. + Mat dummy{nullptr}; + MatCreate(BoutComm::get(), &dummy); + MatSetSizes(dummy, n_local_evolving, n_local_evolving, PETSC_DETERMINE, + PETSC_DETERMINE); + MatSetType(dummy, MATMPIAIJ); + MatSetUp(dummy); + PetscInt Istart, Iend; + MatGetOwnershipRange(dummy, &Istart, &Iend); + MatDestroy(&dummy); + + // ── Compare traversal orderings ─────────────────────────────────────────── + // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that + // the local offset produced by each is identical for every cell. + + struct Mismatch { + int x, y, z; + PetscInt mapping_local; // local offset from mapOwnedInteriorCells + int snes_local; // local offset from globalIndex + }; + std::vector mismatches; + + // Collect (Ind3D -> mapping local offset) from mapOwnedInteriorCells + // so we can look up by field index. + // mapping_local[Ind3D linear] = local offset (0-based on this rank) + const PetscInt rstart = mapping->rowStart(); + std::vector mapping_local_by_ind( + static_cast(mesh->LocalNx * mesh->LocalNy * mesh->LocalNz), -1); + { + PetscInt local_offset = 0; + mapping->mapOwnedInteriorCells( + [&](PetscInt petsc_row, const Ind3D& i, int /*stored*/) { + mapping_local_by_ind[i.ind] = static_cast(petsc_row - rstart); + ++local_offset; + }); + } + + // Walk RGN_NOBNDRY and compare + bool shift_initialised = false; + PetscInt expected_shift = 0; + bool shift_constant = true; + + int snes_local_idx = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] < 0) { + continue; + } + + const int map_local = mapping_local_by_ind[i.ind]; + if (map_local < 0) { + // Cell is in RGN_NOBNDRY but not visited by mapOwnedInteriorCells + mismatches.push_back({i.x(), i.y(), i.z(), -1, snes_local_idx}); + ++snes_local_idx; + continue; + } + + // Check ordering equivalence: local offsets must match + if (map_local != snes_local_idx) { + mismatches.push_back( + {i.x(), i.y(), i.z(), static_cast(map_local), snes_local_idx}); + } + + // Check scalar-shift property: petsc_global - snes_global must be constant + const PetscInt petsc_global = rstart + map_local; + const PetscInt snes_global = Istart + snes_local_idx; + const PetscInt shift = petsc_global - snes_global; + if (!shift_initialised) { + expected_shift = shift; + shift_initialised = true; + } else if (shift != expected_shift) { + shift_constant = false; + } + + ++snes_local_idx; + } + + // ── Gather results across ranks ─────────────────────────────────────────── + const int rank_mismatches = static_cast(mismatches.size()); + int total_mismatches = 0; + MPI_Allreduce(&rank_mismatches, &total_mismatches, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + int rank_shift_fail = shift_constant ? 0 : 1; + int total_shift_fail = 0; + MPI_Allreduce(&rank_shift_fail, &total_shift_fail, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + // ── Report ──────────────────────────────────────────────────────────────── + int my_rank; + MPI_Comm_rank(BoutComm::get(), &my_rank); + + // Each rank writes its own diagnostics; rank 0 also writes the summary. + for (const auto& m : mismatches) { + output.write("MISMATCH rank={} cell=({},{},{}) mapping_local={} snes_local={}\n", + my_rank, m.x, m.y, m.z, m.mapping_local, m.snes_local); + } + if (!shift_constant) { + output.write("SHIFT_NOT_CONSTANT rank={} expected_shift={}\n", my_rank, + expected_shift); + } + + if (my_rank == 0) { + output.write("total_mismatches={}\n", total_mismatches); + output.write("total_shift_failures={}\n", total_shift_fail); + if (total_mismatches == 0 && total_shift_fail == 0) { + output.write("ordering_check=PASS\n"); + } else { + output.write("ordering_check=FAIL\n"); + } + } + + BoutFinalise(); + + return (total_mismatches == 0 && total_shift_fail == 0) ? EXIT_SUCCESS : EXIT_FAILURE; +} + +#endif // BOUT_HAS_PETSC From 7fadf065941e1be85319cd4eff1e7a65b295ea91 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 11 Apr 2026 10:56:49 -0700 Subject: [PATCH 07/11] Minor tidying Include guards, header includes. --- include/bout/petsc_operators.hxx | 9 ++++++++- tests/unit/mesh/test_petsc_operators.cxx | 1 - .../test_petsc_operators_extract_evolving_submatrix.cxx | 1 + 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 14d071fdc6..a4328194a0 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -1,3 +1,8 @@ +#pragma once + +#ifndef BOUT_PETSC_OPERATORS_H +#define BOUT_PETSC_OPERATORS_H + #include "bout/build_defines.hxx" #if BOUT_HAS_PETSC @@ -1055,4 +1060,6 @@ private: #warning PETSc not available. No PetscOperators. -#endif +#endif // BOUT_HAS_PETSC + +#endif //BOUT_PETSC_OPERATORS_H diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 2ba63f4528..853578ba15 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -7,7 +7,6 @@ #include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/output.hxx" -#include "bout/output_bout_types.hxx" #include "bout/petsc_operators.hxx" #include "bout/region.hxx" diff --git a/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx index ea67202863..4ef74318f2 100644 --- a/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx +++ b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx @@ -6,6 +6,7 @@ #include "bout/array.hxx" #include "bout/bout_types.hxx" +#include "bout/globals.hxx" #include "bout/petsc_operators.hxx" #include "bout/region.hxx" From 729f86e8e07d510f803f44c93e54c9a9ad340ba8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 11 Apr 2026 13:56:36 -0700 Subject: [PATCH 08/11] addOperatorSparsity: Add unit tests This function will add non-zeros to a Jacobian matrix, based on a submatrix from a PetscCellOperator. Dummy implementation so unit tests fail. --- CMakeLists.txt | 2 + include/bout/petsc_jacobian.hxx | 32 ++ src/mesh/petsc_jacobian.cxx | 3 + tests/unit/CMakeLists.txt | 1 + tests/unit/mesh/test_petsc_jacobian.cxx | 449 ++++++++++++++++++++++++ 5 files changed, 487 insertions(+) create mode 100644 include/bout/petsc_jacobian.hxx create mode 100644 src/mesh/petsc_jacobian.cxx create mode 100644 tests/unit/mesh/test_petsc_jacobian.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b14e5c3d5..fc73be3c4d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -178,6 +178,7 @@ set(BOUT_SOURCES ./include/bout/paralleltransform.hxx ./include/bout/petsc_interface.hxx ./include/bout/petsc_operators.hxx + ./include/bout/petsc_jacobian.hxx ./include/bout/petsclib.hxx ./include/bout/physicsmodel.hxx ./include/bout/rajalib.hxx @@ -310,6 +311,7 @@ set(BOUT_SOURCES ./src/mesh/parallel_boundary_op.cxx ./src/mesh/parallel_boundary_region.cxx ./src/mesh/petsc_operators.cxx + ./src/mesh/petsc_jacobian.cxx ./src/mesh/surfaceiter.cxx ./src/physics/gyro_average.cxx ./src/physics/physicsmodel.cxx diff --git a/include/bout/petsc_jacobian.hxx b/include/bout/petsc_jacobian.hxx new file mode 100644 index 0000000000..eb9848df3d --- /dev/null +++ b/include/bout/petsc_jacobian.hxx @@ -0,0 +1,32 @@ +#pragma once + +#ifndef BOUT_PETSC_JACOBIAN_H +#define BOUT_PETSC_JACOBIAN_H + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include + +/// Insert the nonzero pattern of @p sub into the variable block +/// (@p out_var, @p in_var) of the Jacobian @p Jfd. +/// +/// @p Jfd is a square matrix of size (n_evolving * nvars) where nvars is +/// inferred as Jfd_global_size / sub_global_size. Each nonzero (r, c) in +/// @p sub produces an entry at (r * nvars + out_var, c * nvars + in_var) +/// in @p Jfd. +/// +/// @p Jfd must already be preallocated. Entries are inserted with +/// INSERT_VALUES; MatAssemblyBegin/End must be called by the caller after +/// all insertions are complete. +/// +/// @param Jfd The Jacobian matrix to populate. Must be preallocated. +/// @param sub Evolving-cell submatrix providing the nonzero pattern. +/// @param out_var Row variable index in [0, nvars). +/// @param in_var Column variable index in [0, nvars). +void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var); + +#endif // BOUT_HAS_PETSC + +#endif // BOUT_PETSC_JACOBIAN_H diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx new file mode 100644 index 0000000000..b820b1a451 --- /dev/null +++ b/src/mesh/petsc_jacobian.cxx @@ -0,0 +1,3 @@ +#include "bout/petsc_jacobian.hxx" + +void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) {} diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 6ea6f92a6b..aca8b833f1 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -100,6 +100,7 @@ set(serial_tests_source ./mesh/test_petsc_operators.cxx ./mesh/test_petsc_operators_make_evolving_is.cxx ./mesh/test_petsc_operators_extract_evolving_submatrix.cxx + ./mesh/test_petsc_jacobian.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_solver.cxx diff --git a/tests/unit/mesh/test_petsc_jacobian.cxx b/tests/unit/mesh/test_petsc_jacobian.cxx new file mode 100644 index 0000000000..e53f6933a4 --- /dev/null +++ b/tests/unit/mesh/test_petsc_jacobian.cxx @@ -0,0 +1,449 @@ +// ============================================================================ +// addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) +// +// Concept: the nonzero pattern of `sub` is scattered into the variable block +// (out_var, in_var) of the pre-allocated Jacobian `Jfd`, with the correct +// index stride. +// +// `sub` is a square matrix of size n_evolving x n_evolving. +// `Jfd` is a square matrix of size (n_evolving * nvars) x (n_evolving * nvars). +// nvars is inferred internally as Jfd_global / sub_global. +// +// A nonzero at (r, c) in sub produces an entry at +// (r * nvars + out_var, c * nvars + in_var) +// in Jfd. +// +// ============================================================================ + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include "bout/petsc_jacobian.hxx" + +#include + +#include +#include +#include + +// ── Helpers ────────────────────────────────────────────────────────────────── + +namespace { + +// Build and assemble a square MPIAIJ matrix of global size n x n, +// with the given (global_row, global_col) nonzero positions (value = 1.0). +// nlocal is the number of locally owned rows on this rank. +// rstart is the first globally owned row. +Mat buildPatternMat(PetscInt n, PetscInt nlocal, PetscInt rstart, + const std::vector>& entries) { + Mat A; + MatCreate(BoutComm::get(), &A); + MatSetSizes(A, nlocal, nlocal, n, n); + MatSetType(A, MATMPIAIJ); + // Generous preallocation: allow up to n nonzeros per row + MatMPIAIJSetPreallocation(A, n, nullptr, n, nullptr); + MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + const PetscScalar one = 1.0; + for (const auto& [r, c] : entries) { + MatSetValues(A, 1, &r, 1, &c, &one, INSERT_VALUES); + } + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + return A; +} + +// Collect all (global_row, global_col) nonzero positions in a matrix, +// iterating only over locally owned rows. +std::vector> matNonzeroPositions(Mat A) { + PetscInt rstart, rend; + MatGetOwnershipRange(A, &rstart, &rend); + std::vector> result; + for (PetscInt row = rstart; row < rend; ++row) { + PetscInt ncols; + const PetscInt* cols; + const PetscScalar* vals; + MatGetRow(A, row, &ncols, &cols, &vals); + for (PetscInt k = 0; k < ncols; ++k) { + result.emplace_back(row, cols[k]); + } + MatRestoreRow(A, row, &ncols, &cols, &vals); + } + return result; +} + +// Build a Jfd of size (n_sub * nvars) x (n_sub * nvars), with nlocal_sub +// rows per variable per rank, and rstart_sub the sub-matrix row offset. +// Returns the matrix plus the Jfd local size and row start. +struct JfdInfo { + Mat jfd; + PetscInt nlocal; + PetscInt rstart; +}; + +JfdInfo buildEmptyJfd(PetscInt n_sub, PetscInt nlocal_sub, PetscInt rstart_sub, + int nvars) { + const PetscInt n = n_sub * nvars; + const PetscInt nlocal = nlocal_sub * nvars; + + Mat jfd; + MatCreate(BoutComm::get(), &jfd); + MatSetSizes(jfd, nlocal, nlocal, n, n); + MatSetType(jfd, MATMPIAIJ); + MatMPIAIJSetPreallocation(jfd, n, nullptr, n, nullptr); + MatSetOption(jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + + PetscInt jfd_rstart, jfd_rend; + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + MatGetOwnershipRange(jfd, &jfd_rstart, &jfd_rend); + return {jfd, nlocal, jfd_rstart}; +} + +// Call addOperatorSparsity and assemble Jfd. +void applyAndAssemble(Mat jfd, Mat sub, int out_var, int in_var) { + addOperatorSparsity(jfd, sub, out_var, in_var); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); +} + +} // namespace + +// ── Fixture ────────────────────────────────────────────────────────────────── + +using PetscAddSparsityTest = FakeMeshFixture; + +// ============================================================================ +// single_variable_diagonal_sub_produces_diagonal_jfd +// +// With nvars=1 a diagonal sub must produce a diagonal Jfd. +// The simplest case: no variable interleaving. +// ============================================================================ +TEST_F(PetscAddSparsityTest, single_variable_diagonal_sub_produces_diagonal_jfd) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; // single rank + const PetscInt rstart_sub = 0; + const int nvars = 1; + + // Diagonal sub: (0,0), (1,1), (2,2) + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, rstart_sub, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = + buildEmptyJfd(n_sub, nlocal_sub, rstart_sub, nvars); + + applyAndAssemble(jfd, sub, /*out_var=*/0, /*in_var=*/0); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + for (const auto& [r, c] : nzs) { + EXPECT_EQ(r, c) << "Expected diagonal entry, got (" << r << "," << c << ")"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// correct_block_is_filled_diagonal_coupling +// +// With nvars=3, out_var=in_var=1 (self-coupling of variable 1), +// a diagonal sub must populate only the (1,1) variable block of Jfd. +// ============================================================================ +TEST_F(PetscAddSparsityTest, correct_block_is_filled_diagonal_coupling) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 1; + const int in_var = 1; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + + for (const auto& [r, c] : nzs) { + // Row must be in the out_var=1 slot: r % nvars == 1 + EXPECT_EQ(out_var, r % nvars) + << "Entry (" << r << "," << c << ") is in wrong row variable"; + // Col must be in the in_var=1 slot: c % nvars == 1 + EXPECT_EQ(in_var, c % nvars) + << "Entry (" << r << "," << c << ") is in wrong col variable"; + // Cell index must match: r / nvars == c / nvars for diagonal sub + EXPECT_EQ(r / nvars, c / nvars) + << "Entry (" << r << "," << c << ") cell indices don't match"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// correct_block_is_filled_off_diagonal_coupling +// +// With nvars=3, out_var=0, in_var=2 (cross-variable coupling), +// all Jfd entries must be in the (0,2) variable block. +// ============================================================================ +TEST_F(PetscAddSparsityTest, correct_block_is_filled_off_diagonal_coupling) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 0; + const int in_var = 2; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + + for (const auto& [r, c] : nzs) { + EXPECT_EQ(out_var, r % nvars) + << "Entry (" << r << "," << c << ") is in wrong row variable"; + EXPECT_EQ(in_var, c % nvars) + << "Entry (" << r << "," << c << ") is in wrong col variable"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// stride_is_correct_for_nvars_2 +// +// With nvars=2 and a sub containing entries (0,0) and (0,1), the Jfd +// entry positions must follow the stride formula: +// Jfd row = sub_row * nvars + out_var +// Jfd col = sub_col * nvars + in_var +// ============================================================================ +TEST_F(PetscAddSparsityTest, stride_is_correct_for_nvars_2) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + const int out_var = 0; + const int in_var = 1; + + // Sub entries: (0,0), (1,2), (2,1) — arbitrary non-diagonal pattern + const std::vector> sub_entries = {{0, 0}, {1, 2}, {2, 1}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + // Build expected Jfd positions using the stride formula + std::vector> expected; + for (const auto& [sr, sc] : sub_entries) { + expected.emplace_back((sr * nvars) + out_var, (sc * nvars) + in_var); + } + std::sort(expected.begin(), expected.end()); + + auto nzs = matNonzeroPositions(jfd); + std::sort(nzs.begin(), nzs.end()); + + ASSERT_EQ(expected.size(), nzs.size()); + for (std::size_t k = 0; k < expected.size(); ++k) { + EXPECT_EQ(expected[k].first, nzs[k].first) << "Row mismatch at entry " << k; + EXPECT_EQ(expected[k].second, nzs[k].second) << "Col mismatch at entry " << k; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// no_other_blocks_written +// +// With nvars=3 and only (out_var=1, in_var=1) populated, the other eight +// variable blocks must remain empty. +// ============================================================================ +TEST_F(PetscAddSparsityTest, no_other_blocks_written) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 1; + const int in_var = 1; + + const std::vector> sub_entries = { + {0, 0}, {1, 1}, {2, 2}, {3, 3}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "Wrong number of entries in Jfd"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(out_var, r % nvars) << "Unexpected entry in row-variable " << r % nvars + << " at (" << r << "," << c << ")"; + EXPECT_EQ(in_var, c % nvars) << "Unexpected entry in col-variable " << c % nvars + << " at (" << r << "," << c << ")"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// multiple_calls_union_patterns +// +// Two calls with different (out_var, in_var) pairs must union their patterns: +// the total nonzero count equals the sum of each operator's count, and the +// entries from each call are in the correct block. +// ============================================================================ +TEST_F(PetscAddSparsityTest, multiple_calls_union_patterns) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + // Sub A: diagonal + const std::vector> entries_a = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub_a = buildPatternMat(n_sub, nlocal_sub, 0, entries_a); + + // Sub B: one off-diagonal entry per row + const std::vector> entries_b = {{0, 1}, {1, 2}, {2, 0}}; + Mat sub_b = buildPatternMat(n_sub, nlocal_sub, 0, entries_b); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + // Insert both without assembling between calls + addOperatorSparsity(jfd, sub_a, /*out_var=*/0, /*in_var=*/0); + addOperatorSparsity(jfd, sub_b, /*out_var=*/1, /*in_var=*/1); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + + // Total count: 3 from sub_a + 3 from sub_b = 6 + EXPECT_EQ(6U, nzs.size()); + + // Check each entry is in the correct variable block + for (const auto& [r, c] : nzs) { + const int rv = static_cast(r % nvars); + const int cv = static_cast(c % nvars); + const bool in_a_block = (rv == 0 && cv == 0); + const bool in_b_block = (rv == 1 && cv == 1); + EXPECT_TRUE(in_a_block || in_b_block) + << "Entry (" << r << "," << c << ") is in unexpected variable block (" << rv + << "," << cv << ")"; + } + + MatDestroy(&sub_a); + MatDestroy(&sub_b); + MatDestroy(&jfd); +} + +// ============================================================================ +// empty_sub_adds_no_entries +// +// A zero-pattern sub must not modify Jfd. +// ============================================================================ +TEST_F(PetscAddSparsityTest, empty_sub_adds_no_entries) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, {}); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, 0, 0); + + const auto nzs = matNonzeroPositions(jfd); + EXPECT_EQ(0U, nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// nonzero_count_equals_sub_nonzero_count +// +// For a single call the total Jfd nonzero count must equal the sub nonzero +// count exactly — no extras, no missing. +// ============================================================================ +TEST_F(PetscAddSparsityTest, nonzero_count_equals_sub_nonzero_count) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + // Arbitrary non-trivial pattern + const std::vector> sub_entries = { + {0, 0}, {0, 1}, {1, 1}, {1, 3}, {2, 2}, {3, 0}, {3, 3}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, 0, 0); + + MatInfo info; + MatGetInfo(jfd, MAT_GLOBAL_SUM, &info); + EXPECT_EQ(static_cast(sub_entries.size()), + static_cast(info.nz_used)); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// first_and_last_variable_blocks_correct +// +// Explicitly check out_var=0,in_var=0 and out_var=nvars-1,in_var=nvars-1 +// to confirm the stride formula is correct at the boundaries. +// ============================================================================ +TEST_F(PetscAddSparsityTest, first_and_last_variable_blocks_correct) { + const PetscInt n_sub = 2; + const PetscInt nlocal_sub = n_sub; + const int nvars = 4; + + const std::vector> sub_entries = { + {0, 0}, {0, 1}, {1, 0}, {1, 1}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Test first block (0,0) + { + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + applyAndAssemble(jfd, sub, 0, 0); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "First block: wrong number of entries"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(0, r % nvars) << "First block row variable wrong"; + EXPECT_EQ(0, c % nvars) << "First block col variable wrong"; + } + MatDestroy(&jfd); + } + + // Test last block (nvars-1, nvars-1) + { + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + applyAndAssemble(jfd, sub, nvars - 1, nvars - 1); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "Last block: wrong number of entries"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(nvars - 1, r % nvars) << "Last block row variable wrong"; + EXPECT_EQ(nvars - 1, c % nvars) << "Last block col variable wrong"; + } + MatDestroy(&jfd); + } + + MatDestroy(&sub); +} + +#endif // BOUT_HAS_PETSC From ba6996e22e0e1897d7638f02defcb3660f918a54 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 11 Apr 2026 17:11:14 -0700 Subject: [PATCH 09/11] addOperatorSparsity implementation Sets a sparsity pattern in a given Jacobian matrix, based on an operator matrix and the variable indices. Passes unit tests. --- src/mesh/petsc_jacobian.cxx | 47 +++++++++++++++++++++++++++++++++++- src/mesh/petsc_operators.cxx | 1 - 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx index b820b1a451..b0a4eec50c 100644 --- a/src/mesh/petsc_jacobian.cxx +++ b/src/mesh/petsc_jacobian.cxx @@ -1,3 +1,48 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "bout/assert.hxx" #include "bout/petsc_jacobian.hxx" +#include "bout/petsclib.hxx" + +#include + +void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) { + // Infer nvars from global sizes + PetscInt jfd_global{0}, sub_global{0}; + BOUT_DO_PETSC(MatGetSize(Jfd, &jfd_global, nullptr)); + BOUT_DO_PETSC(MatGetSize(sub, &sub_global, nullptr)); + + ASSERT1(sub_global > 0); + ASSERT1(jfd_global % sub_global == 0); + + const PetscInt nvars = jfd_global / sub_global; + + ASSERT1(out_var >= 0 && out_var < static_cast(nvars)); + ASSERT1(in_var >= 0 && in_var < static_cast(nvars)); + + // Iterate over locally owned rows of sub and insert into Jfd + PetscInt rstart{0}, rend{0}; + BOUT_DO_PETSC(MatGetOwnershipRange(sub, &rstart, &rend)); + + const PetscScalar one = 1.0; + + for (PetscInt sub_row = rstart; sub_row < rend; ++sub_row) { + PetscInt ncols{0}; + const PetscInt* sub_cols{nullptr}; + const PetscScalar* vals{nullptr}; + BOUT_DO_PETSC(MatGetRow(sub, sub_row, &ncols, &sub_cols, &vals)); + + const PetscInt jfd_row = (sub_row * nvars) + out_var; + + for (PetscInt k = 0; k < ncols; ++k) { + const PetscInt jfd_col = (sub_cols[k] * nvars) + in_var; + BOUT_DO_PETSC(MatSetValues(Jfd, 1, &jfd_row, 1, &jfd_col, &one, INSERT_VALUES)); + } + + BOUT_DO_PETSC(MatRestoreRow(sub, sub_row, &ncols, &sub_cols, &vals)); + } +} -void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) {} +#endif // BOUT_HAS_PETSC diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 9700001d23..d5bfe87bd5 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -6,7 +6,6 @@ #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" From 9de5da3155c7269220f2ca16b60cb2ec27478968 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 11 Apr 2026 17:45:27 -0700 Subject: [PATCH 10/11] addOperatorSparsity(J, mat) unit tests Should apply the connectivity pattern to all pairs of variables in the system Jacobian. --- include/bout/petsc_jacobian.hxx | 4 + src/mesh/petsc_jacobian.cxx | 2 + tests/unit/mesh/test_petsc_jacobian.cxx | 189 ++++++++++++++++++++++++ 3 files changed, 195 insertions(+) diff --git a/include/bout/petsc_jacobian.hxx b/include/bout/petsc_jacobian.hxx index eb9848df3d..19588a8e47 100644 --- a/include/bout/petsc_jacobian.hxx +++ b/include/bout/petsc_jacobian.hxx @@ -27,6 +27,10 @@ /// @param in_var Column variable index in [0, nvars). void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var); +/// @param Jfd The Jacobian matrix to populate. Must be preallocated. +/// @param sub Evolving-cell submatrix providing the nonzero pattern. +void addOperatorSparsity(Mat Jfd, Mat sub); + #endif // BOUT_HAS_PETSC #endif // BOUT_PETSC_JACOBIAN_H diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx index b0a4eec50c..e2648da6a4 100644 --- a/src/mesh/petsc_jacobian.cxx +++ b/src/mesh/petsc_jacobian.cxx @@ -45,4 +45,6 @@ void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) { } } +void addOperatorSparsity(Mat Jfd, Mat sub) {} + #endif // BOUT_HAS_PETSC diff --git a/tests/unit/mesh/test_petsc_jacobian.cxx b/tests/unit/mesh/test_petsc_jacobian.cxx index e53f6933a4..ac2e8e5d08 100644 --- a/tests/unit/mesh/test_petsc_jacobian.cxx +++ b/tests/unit/mesh/test_petsc_jacobian.cxx @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -446,4 +447,192 @@ TEST_F(PetscAddSparsityTest, first_and_last_variable_blocks_correct) { MatDestroy(&sub); } +// ============================================================================ +// all_pairs_single_variable +// +// With nvars=1 the all-pairs overload is identical to addOperatorSparsity +// with out_var=0, in_var=0. +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_single_variable) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 1; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_fills_every_variable_block +// +// With nvars=3 every one of the 9 variable blocks must contain exactly the +// same pattern as sub. The total nonzero count must equal +// nvars * nvars * nnz(sub). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_fills_every_variable_block) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + const std::size_t expected_total = + static_cast(nvars * nvars) * sub_entries.size(); + ASSERT_EQ(expected_total, nzs.size()); + + // Every (out_var, in_var) pair must appear exactly once per sub entry. + // Count how many Jfd entries land in each variable block. + std::vector> block_counts(nvars, std::vector(nvars, 0)); + for (const auto& [r, c] : nzs) { + block_counts[r % nvars][c % nvars]++; + } + for (int ov = 0; ov < nvars; ++ov) { + for (int iv = 0; iv < nvars; ++iv) { + EXPECT_EQ(static_cast(sub_entries.size()), block_counts[ov][iv]) + << "Block (" << ov << "," << iv << ") has wrong entry count"; + } + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_entries_match_single_pair_calls +// +// The all-pairs overload must produce exactly the same Jfd nonzero positions +// as calling the single-pair overload for every (out_var, in_var). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_entries_match_single_pair_calls) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + // Non-trivial pattern so the stride mapping is exercised + const std::vector> sub_entries = { + {0, 0}, {0, 2}, {1, 1}, {2, 0}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Build expected positions by calling the single-pair version for each block + std::vector> expected; + { + auto [jfd_ref, nlocal_ref, rstart_ref] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + for (int ov = 0; ov < nvars; ++ov) { + for (int iv = 0; iv < nvars; ++iv) { + addOperatorSparsity(jfd_ref, sub, ov, iv); + } + } + MatAssemblyBegin(jfd_ref, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd_ref, MAT_FINAL_ASSEMBLY); + expected = matNonzeroPositions(jfd_ref); + MatDestroy(&jfd_ref); + } + + // Build actual positions using the all-pairs overload + std::vector> actual; + { + auto [jfd_act, nlocal_act, rstart_act] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + addOperatorSparsity(jfd_act, sub); + MatAssemblyBegin(jfd_act, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd_act, MAT_FINAL_ASSEMBLY); + actual = matNonzeroPositions(jfd_act); + MatDestroy(&jfd_act); + } + + std::sort(expected.begin(), expected.end()); + std::sort(actual.begin(), actual.end()); + + ASSERT_EQ(expected.size(), actual.size()); + for (std::size_t k = 0; k < expected.size(); ++k) { + EXPECT_EQ(expected[k], actual[k]) << "Mismatch at position " << k; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// all_pairs_empty_sub_fills_nothing +// +// An empty sub must leave Jfd empty regardless of nvars. +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_empty_sub_fills_nothing) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, {}); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + EXPECT_EQ(0U, nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_no_entries_outside_expected_positions +// +// Every entry in Jfd must lie at a position consistent with the stride formula +// applied to some (out_var, in_var) pair. Specifically, for each entry +// (r, c) in Jfd there must exist a (sr, sc) in sub such that +// r == sr * nvars + (r % nvars) and c == sc * nvars + (c % nvars). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_no_entries_outside_expected_positions) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + std::vector> sub_entries = { + {0, 1}, {1, 0}, {1, 2}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Pre-compute sub cell pairs for fast lookup + const std::set> sub_set(sub_entries.begin(), + sub_entries.end()); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + for (const auto& [r, c] : nzs) { + const PetscInt sr = r / nvars; + const PetscInt sc = c / nvars; + EXPECT_TRUE(sub_set.count({sr, sc}) > 0) + << "Entry (" << r << "," << c << ") maps to sub cell (" << sr << "," << sc + << ") which is not in sub"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + #endif // BOUT_HAS_PETSC From 631dcddd7deb73769b92f9ff0c224020fbdc6466 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 12 Apr 2026 21:55:16 -0700 Subject: [PATCH 11/11] addOperatorSparsity(J, mat) implementation Applies the same connectivity in `mat` to all pairs of variables in `J`. Passes unit tests. --- include/bout/petsc_jacobian.hxx | 7 +++++++ src/mesh/petsc_jacobian.cxx | 17 ++++++++++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/include/bout/petsc_jacobian.hxx b/include/bout/petsc_jacobian.hxx index 19588a8e47..73d4d2fc71 100644 --- a/include/bout/petsc_jacobian.hxx +++ b/include/bout/petsc_jacobian.hxx @@ -27,6 +27,13 @@ /// @param in_var Column variable index in [0, nvars). void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var); +/// @brief Insert the nonzero pattern of @p sub into every variable block of +/// @p Jfd. +/// +/// Equivalent to calling addOperatorSparsity(Jfd, sub, out_var, in_var) for +/// every combination of @p out_var and @p in_var in [0, nvars), where +/// @c nvars is inferred as @c Jfd_global / @c sub_global. +/// /// @param Jfd The Jacobian matrix to populate. Must be preallocated. /// @param sub Evolving-cell submatrix providing the nonzero pattern. void addOperatorSparsity(Mat Jfd, Mat sub); diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx index e2648da6a4..457c667de5 100644 --- a/src/mesh/petsc_jacobian.cxx +++ b/src/mesh/petsc_jacobian.cxx @@ -45,6 +45,21 @@ void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) { } } -void addOperatorSparsity(Mat Jfd, Mat sub) {} +void addOperatorSparsity(Mat Jfd, Mat sub) { + PetscInt jfd_global{0}, sub_global{0}; + MatGetSize(Jfd, &jfd_global, nullptr); + MatGetSize(sub, &sub_global, nullptr); + + ASSERT1(sub_global > 0); + ASSERT1(jfd_global % sub_global == 0); + + const int nvars = static_cast(jfd_global / sub_global); + + for (int out_var = 0; out_var < nvars; ++out_var) { + for (int in_var = 0; in_var < nvars; ++in_var) { + addOperatorSparsity(Jfd, sub, out_var, in_var); + } + } +} #endif // BOUT_HAS_PETSC