From 325019aeacbae942d680dfb7671a98cfe72eff3d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Apr 2026 14:29:21 -0700 Subject: [PATCH 01/18] 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 02/18] 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 03/18] 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 04/18] 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 05/18] 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 06/18] 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 07/18] 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 08/18] 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 09/18] 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 10/18] 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 From 6441eadd81ef2c7915857f1949876802954c3a72 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 29 Apr 2026 10:33:46 -0700 Subject: [PATCH 11/18] PetscOperators: Add debugging print Test passes locally but fails in CI --- src/mesh/petsc_operators.cxx | 1 + tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index d5bfe87bd5..aadad200c4 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -60,6 +60,7 @@ void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, const PetscInt petsc_row = row_start + i; const PetscInt stored_col = stored_indices[i]; ASSERT1(stored_col >= 0); + output.write("{}: {}, {} | {} {}\n", i, petsc_row, stored_col, nlocal, nglobal); local_stored_to_petsc[static_cast(stored_col)] = petsc_row; BOUT_DO_PETSC(MatSetValues(mat_stored_to_petsc, 1, &petsc_row, 1, &stored_col, &one, INSERT_VALUES)); diff --git a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx index c3962a1b3f..19420b5dfa 100644 --- a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx +++ b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx @@ -184,8 +184,7 @@ int main(int argc, char** argv) { BoutComm::get()); // ── Report ──────────────────────────────────────────────────────────────── - int my_rank; - MPI_Comm_rank(BoutComm::get(), &my_rank); + const int my_rank = BoutComm::rank(); // Each rank writes its own diagnostics; rank 0 also writes the summary. for (const auto& m : mismatches) { From cafee4bdac645f9d1c3cd7f26dd3b2af75d71d53 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 29 Apr 2026 22:18:53 -0700 Subject: [PATCH 12/18] test-petsc-ordering: Fix cell numbering Use global numbering --- .../fci-operators/fci_operators_example.cxx | 63 ++-- .../test_petsc_ordering.cxx | 285 +++++++++--------- 2 files changed, 178 insertions(+), 170 deletions(-) diff --git a/examples/fci-operators/fci_operators_example.cxx b/examples/fci-operators/fci_operators_example.cxx index 3366622ff0..80625db9fc 100644 --- a/examples/fci-operators/fci_operators_example.cxx +++ b/examples/fci-operators/fci_operators_example.cxx @@ -10,52 +10,53 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); - const PetscOperators ops(bout::globals::mesh); + { + const PetscOperators ops(bout::globals::mesh); - // Parallel operators - const auto parallel = ops.getParallel(); + // Parallel operators + const auto parallel = ops.getParallel(); - // Test field - Field3D f = FieldFactory::get()->create3D("sin(y) + cos(z)"); - bout::globals::mesh->communicate(f); - f.applyParallelBoundary("parallel_dirichlet_o2"); + // Test field + Field3D f = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f); + f.applyParallelBoundary("parallel_dirichlet_o2"); - Field3D f_neumann = FieldFactory::get()->create3D("sin(y) + cos(z)"); - bout::globals::mesh->communicate(f_neumann); - f_neumann.applyParallelBoundary("parallel_neumann_o1"); + Field3D f_neumann = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f_neumann); + f_neumann.applyParallelBoundary("parallel_neumann_o1"); - auto forward_op = parallel.Restrict_minus * parallel.Forward; + auto forward_op = parallel.Restrict_minus * parallel.Forward; - Options dump; - dump["f"] = f; - dump["dV"] = parallel.dV; + Options dump; + dump["f"] = f; + dump["dV"] = parallel.dV; - // Test1 : (parallel.Restrict_minus * parallel.Inject_plus) is the identity in the interior (not X boundaries) + // Test1 : (parallel.Restrict_minus * parallel.Inject_plus) is the identity in the interior (not X boundaries) - dump["grad_par_op"] = parallel.Grad_par(f); - Field3D forward{0.0}; - BOUT_FOR(i, forward.getRegion("RGN_NOBNDRY")) { forward[i] = f.yup()[i.yp()]; } + dump["grad_par_op"] = parallel.Grad_par(f); + Field3D forward{0.0}; + BOUT_FOR(i, forward.getRegion("RGN_NOBNDRY")) { forward[i] = f.yup()[i.yp()]; } - dump["forward"] = forward; + dump["forward"] = forward; - dump["forward_op"] = forward_op(f); + dump["forward_op"] = forward_op(f); - dump["grad_par_yud"] = Grad_par(f); + dump["grad_par_yud"] = Grad_par(f); - // Test2 : Sum of dV * Div_par over domain is zero + // Test2 : Sum of dV * Div_par over domain is zero - dump["div_par_op"] = parallel.Div_par(f); + dump["div_par_op"] = parallel.Div_par(f); - auto* coords = bout::globals::mesh->getCoordinates(); - bout::globals::mesh->communicate(coords->Bxy); - coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); - dump["div_par_yud"] = Div_par(f); + auto* coords = bout::globals::mesh->getCoordinates(); + bout::globals::mesh->communicate(coords->Bxy); + coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); + dump["div_par_yud"] = Div_par(f); - dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); - dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); - - bout::writeDefaultOutputFile(dump); + dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); + dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); + bout::writeDefaultOutputFile(dump); + } BoutFinalise(); return 0; } diff --git a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx index 19420b5dfa..6a0fa25d0b 100644 --- a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx +++ b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx @@ -1,4 +1,4 @@ -/// Integrated test for Unit 3: ordering equivalence between PetscCellMapping +/// Integrated test for ordering equivalence between PetscCellMapping /// and the SNES solver's globalIndex traversal. /// /// For every locally owned evolving interior cell this program checks: @@ -42,170 +42,177 @@ int main() { return EXIT_SUCCESS; } #include int main(int argc, char** argv) { + int total_mismatches = 0; + int total_shift_fail = 0; + // 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++; + 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}; + { + const int ngy = mesh->GlobalNyNoBoundaries; + const int ngz = mesh->GlobalNzNoBoundaries; + + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + const int global_x = + mesh->getGlobalXIndexNoBoundaries(x); // Exclude boundary cells + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + const int global_y = mesh->getGlobalYIndexNoBoundaries(y); + for (int z = 0; z < mesh->LocalNz; ++z) { + const int global_z = mesh->getGlobalZIndexNoBoundaries(z); + cell_number(x, y, z) = ((global_x * ngy) + global_y) * ngz + global_z; + } } } + // Communicate so guard cells are filled (needed by some mesh operations) + mesh->communicate(cell_number); } - // 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}; + 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; + // 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) { - snes_index[i] = local_idx++; + ++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; - }); - } + // 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; + // 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; - } + 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; - } + 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 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; - } + // 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; - } + ++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()); + // ── Gather results across ranks ─────────────────────────────────────────── + const int rank_mismatches = static_cast(mismatches.size()); + 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()); + int rank_shift_fail = shift_constant ? 0 : 1; + MPI_Allreduce(&rank_shift_fail, &total_shift_fail, 1, MPI_INT, MPI_SUM, + BoutComm::get()); - // ── Report ──────────────────────────────────────────────────────────────── - const int my_rank = BoutComm::rank(); + // ── Report ──────────────────────────────────────────────────────────────── + const int my_rank = BoutComm::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); - } + // 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"); + 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; From fdf757af26af8036fdfc75d7525c10d0c75807e9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 24 Jun 2026 17:26:07 -0700 Subject: [PATCH 13/18] Fix 2D metric build with BinaryExpr Field2D expressions need to be evaluated before assignment to Field3D. --- src/mesh/petsc_operators.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 640a831156..9827512a5f 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -359,14 +359,14 @@ PetscOperators::Parallel PetscOperators::getParallel() const { auto* coords = mesh->getCoordinates(); // Parallel spacing in cell space - Field3D dl = coords->dy * sqrt(coords->g_22); + Field3D dl = Coordinates::FieldMetric{coords->dy * sqrt(coords->g_22)}; dl.splitParallelSlices(); dl.yup() = 0.0; dl.ydown() = 0.0; dl.applyParallelBoundary("parallel_neumann_o1"); // Cell volume - Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; + Field3D dV = Coordinates::FieldMetric{coords->J * coords->dx * coords->dy * coords->dz}; dV.splitParallelSlices(); dV.yup() = 0.0; dV.ydown() = 0.0; From fee5e520c4c046f014e02209fbd39de7d4aa593c Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 24 Jun 2026 17:32:24 -0700 Subject: [PATCH 14/18] Mesh::get(Array) source guards Handle the case where `source` is `nullptr`. --- src/mesh/mesh.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index c44a62aa7e..33368e6d6a 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -189,10 +189,16 @@ int Mesh::get(bool& bval, const std::string& name, bool def) { } int Mesh::get(Array& var, const std::string& name) { + if (source == nullptr) { + return 1; + } return source->get(var, name) ? 0 : 1; } int Mesh::get(Array& var, const std::string& name) { + if (source == nullptr) { + return 1; + } return source->get(var, name) ? 0 : 1; } From c8fc9b9bbe3c4bad99356223a5e8d6a16230073e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 24 Jun 2026 17:37:54 -0700 Subject: [PATCH 15/18] addOperatorSparsity: Document assumptions & limitations Assumes no evolving Field2D or evolving boundary cells. --- include/bout/petsc_jacobian.hxx | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/include/bout/petsc_jacobian.hxx b/include/bout/petsc_jacobian.hxx index 73d4d2fc71..120279c683 100644 --- a/include/bout/petsc_jacobian.hxx +++ b/include/bout/petsc_jacobian.hxx @@ -17,6 +17,13 @@ /// @p sub produces an entry at (r * nvars + out_var, c * nvars + in_var) /// in @p Jfd. /// +/// Limitation: this helper assumes that @p Jfd uses a uniform per-cell +/// interleaving of the form global_row = cell * nvars + var over the +/// evolving cells only. It is therefore not valid for solver layouts that mix +/// Field2D and Field3D variables, include evolving boundary variables, or use +/// any other non-uniform cell-to-row mapping. Callers must check that the +/// solver Jacobian ordering matches this assumption before using this helper. +/// /// @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. @@ -34,6 +41,10 @@ void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var); /// 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. /// +/// The same layout restriction applies as for the block-specific overload: +/// @p Jfd must use the uniform cell * nvars + var ordering over +/// evolving cells only. +/// /// @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); From 50de5cb0d24f75c23ac9be3eab37a8a715913fa0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 1 Jul 2026 15:47:07 -0700 Subject: [PATCH 16/18] Solver::addJacobianPattern(PetscOperator) Allows the sparsity pattern in PetscPreconditioner for SNES, CVODE and PETSc interfaces to be set using a PetscOperator. This should enable preconditioning of 3D FCI simulations. --- include/bout/petsc_operators.hxx | 10 +++ include/bout/solver.hxx | 67 ++++++++++++++++++ src/solver/impls/cvode/cvode.cxx | 1 + src/solver/impls/cvode/cvode.hxx | 7 ++ src/solver/impls/petsc/petsc.cxx | 1 + src/solver/impls/petsc/petsc.hxx | 5 ++ src/solver/impls/snes/snes.cxx | 1 + src/solver/impls/snes/snes.hxx | 7 ++ src/solver/solver.cxx | 109 ++++++++++++++++++++++++++++++ tests/unit/solver/test_solver.cxx | 33 +++++++++ 10 files changed, 241 insertions(+) diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index e17880ed3c..83069b9f60 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -709,6 +709,16 @@ public: /// low-level callers that need direct PETSc access. Mat raw() const { return *this->mat_operator; } + /// @brief Return the shared mapping for the operator's output (row) space. + const std::shared_ptr& getOutMapping() const { + return out_mapping; + } + + /// @brief Return the shared mapping for the operator's input (column) space. + const std::shared_ptr& getInMapping() const { + return in_mapping; + } + /// @brief Construct a diagonal operator from a vector of diagonal entries. /// /// Only available when @p OutSpaceTag == @p InSpaceTag (square, same-space diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 09ede6a32b..0a5a53ca81 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -72,6 +72,10 @@ using TimestepMonitorFunc = int (*)(Solver* solver, BoutReal simtime, BoutReal l #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/generic_factory.hxx" +#if BOUT_HAS_PETSC +#include "bout/petsc_interface.hxx" +#include "bout/petsc_operators.hxx" +#endif #include "bout/vector2d.hxx" #include "bout/vector3d.hxx" @@ -81,6 +85,7 @@ using TimestepMonitorFunc = int (*)(Solver* solver, BoutReal simtime, BoutReal l #include #include +#include using SolverType = std::string; constexpr auto SOLVERCVODE = "cvode"; @@ -210,6 +215,29 @@ using RegisterUnavailableSolver = SolverFactory::RegisterUnavailableInFactory; */ class Solver { public: + class VarRef { + public: + static constexpr int AllValue = -1; + static constexpr int InvalidValue = -2; + + constexpr VarRef() = default; + + static constexpr VarRef All() { return VarRef(AllValue); } + static constexpr VarRef Invalid() { return VarRef(InvalidValue); } + + constexpr bool isAll() const { return value == AllValue; } + constexpr bool isInvalid() const { return value == InvalidValue; } + constexpr bool isConcrete() const { return value >= 0; } + constexpr bool isValid() const { return isAll() || isConcrete(); } + constexpr int index() const { return value; } + + private: + explicit constexpr VarRef(int value) : value(value) {} + + int value{InvalidValue}; + friend class Solver; + }; + Solver(Options* opts = nullptr); virtual ~Solver() = default; @@ -265,6 +293,21 @@ public: virtual void add(Vector3D& v, const std::string& name, const std::string& description = ""); + /// Get the solver-variable reference for a scalar field/component by name. + VarRef getVarRef(std::string_view name) const; + +#if BOUT_HAS_PETSC + /// Register a Jacobian sparsity contribution for all variable blocks. + bool addJacobianPattern(const PetscCellOperator& op); + + /// Register a Jacobian sparsity contribution for one or more variable blocks. + /// + /// Either variable reference may be VarRef::All(), which expands when the + /// Jacobian matrix is created during solver initialisation. + virtual bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var); +#endif + /// Returns true if constraints available virtual bool constraints() { return has_constraints; } @@ -688,6 +731,26 @@ private: /// Physics model being evolved PhysicsModel* model{nullptr}; +protected: +#if BOUT_HAS_PETSC + struct DeferredJacobianPattern { + bout::petsc::UniqueMat submatrix{new Mat{nullptr}}; + VarRef out_var{}; + VarRef in_var{}; + }; + + /// Queue a Jacobian-pattern contribution for PETSc-preconditioner-based solvers. + bool queueJacobianPattern(const PetscCellOperator& op, VarRef out_var, VarRef in_var); + + /// Insert any queued Jacobian-pattern contributions into an existing Jacobian matrix. + void applyQueuedJacobianPatterns(Mat Jfd) const; + + /// Check whether the current solver variable layout is compatible with + /// addOperatorSparsity(). + bool canApplyQueuedJacobianPatterns() const; +#endif + +private: /// Should non-split physics models be treated as diffusive? bool is_nonsplit_model_diffusive{true}; @@ -704,6 +767,10 @@ private: /// List of timestep monitor functions std::list timestep_monitors; +#if BOUT_HAS_PETSC + std::vector deferred_jacobian_patterns; +#endif + /// Should be run before user RHS is called void pre_rhs(BoutReal t); /// Should be run after user RHS is called diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index 2cdeeed9c3..2edac8e292 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -491,6 +491,7 @@ int CvodeSolver::init() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, local_N, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); PetscCall( petsc_preconditioner.updateColoring(CvodeSolver::petscFormFunction, this)); PetscCall(MatFDColoringSetF(petsc_preconditioner.coloring(), petsc_f)); diff --git a/src/solver/impls/cvode/cvode.hxx b/src/solver/impls/cvode/cvode.hxx index c242e4128a..3e8c116718 100644 --- a/src/solver/impls/cvode/cvode.hxx +++ b/src/solver/impls/cvode/cvode.hxx @@ -79,12 +79,19 @@ class CvodeSolver : public Solver { public: explicit CvodeSolver(Options* opts = nullptr); ~CvodeSolver() override; + using Solver::addJacobianPattern; BoutReal getCurrentTimestep() override { return hcur; } int init() override; int run() override; BoutReal run(BoutReal tout); +#if BOUT_HAS_PETSC + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } +#endif void resetInternalFields() override; diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index d123b56508..e6fbe95249 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -519,6 +519,7 @@ int PetscSolver::init() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, nlocal, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); output_progress.write("Creating Jacobian coloring\n"); updateColoring(); } else { diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index ab8e5bf9e0..2438f7de1c 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -63,9 +63,14 @@ class PetscSolver : public Solver { public: PetscSolver(Options* opts = nullptr); ~PetscSolver(); + using Solver::addJacobianPattern; int init() override; int run() override; + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } // These functions used internally (but need to be public) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 07d16652f2..415eb77a09 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -78,6 +78,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, nlocal, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); output_progress.write("Creating Jacobian coloring\n"); PetscCall(petsc_preconditioner.updateColoring(FormFunctionForColoring, this)); diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 6cdaac082c..21f1d72f01 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -78,9 +78,16 @@ class SNESSolver : public Solver { public: explicit SNESSolver(Options* opts = nullptr); ~SNESSolver() override = default; + using Solver::addJacobianPattern; int init() override; int run() override; +#if BOUT_HAS_PETSC + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } +#endif /// Nonlinear function. This is called by PETSc SNES object /// via a static C-style function. For implicit diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 1f7b654427..f6c032268e 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -47,6 +47,10 @@ #include "bout/vector2d.hxx" #include "bout/vector3d.hxx" +#if BOUT_HAS_PETSC +#include "bout/petsc_jacobian.hxx" +#endif + #include #include @@ -350,6 +354,111 @@ void Solver::add(Vector3D& v, const std::string& name, const std::string& descri v3d.emplace_back(std::move(d)); } +Solver::VarRef Solver::getVarRef(std::string_view name) const { + int index = 0; + for (const auto& field : f2d) { + if (field.name == name) { + return VarRef(index); + } + ++index; + } + for (const auto& field : f3d) { + if (field.name == name) { + return VarRef(index); + } + ++index; + } + return VarRef::Invalid(); +} + +#if BOUT_HAS_PETSC +bool Solver::addJacobianPattern(const PetscCellOperator& op) { + return addJacobianPattern(op, VarRef::All(), VarRef::All()); +} + +bool Solver::addJacobianPattern(const PetscCellOperator& UNUSED(op), + VarRef UNUSED(out_var), VarRef UNUSED(in_var)) { + return false; +} + +bool Solver::queueJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) { + if (initialised || !out_var.isValid() || !in_var.isValid()) { + return false; + } + + auto out_mapping = + std::dynamic_pointer_cast(op.getOutMapping()); + auto in_mapping = std::dynamic_pointer_cast(op.getInMapping()); + if (!out_mapping || !in_mapping || (out_mapping != in_mapping)) { + return false; + } + + DeferredJacobianPattern pattern; + *pattern.submatrix = out_mapping->extractEvolvingSubmatrix(op); + pattern.out_var = out_var; + pattern.in_var = in_var; + deferred_jacobian_patterns.emplace_back(std::move(pattern)); + return true; +} + +bool Solver::canApplyQueuedJacobianPatterns() const { + if (deferred_jacobian_patterns.empty()) { + return true; + } + + if (n2Dvars() != 0 || n3Dvars() == 0) { + return false; + } + + const auto evolve_boundary = [](const auto& field) { return field.evolve_bndry; }; + return std::none_of(f2d.begin(), f2d.end(), evolve_boundary) + && std::none_of(f3d.begin(), f3d.end(), evolve_boundary); +} + +void Solver::applyQueuedJacobianPatterns(Mat Jfd) const { + if (deferred_jacobian_patterns.empty()) { + return; + } + + if (!canApplyQueuedJacobianPatterns()) { + throw BoutException( + "Queued Jacobian patterns require a uniform evolving-cell layout with only " + "Field3D variables and no evolving boundary cells"); + } + + const int nvars = n2Dvars() + n3Dvars(); + ASSERT1(nvars > 0); + + BOUT_DO_PETSC(MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); + + for (const auto& pattern : deferred_jacobian_patterns) { + if (pattern.out_var.isAll() && pattern.in_var.isAll()) { + addOperatorSparsity(Jfd, *pattern.submatrix); + continue; + } + if (pattern.out_var.isAll()) { + for (int out_var = 0; out_var < nvars; ++out_var) { + addOperatorSparsity(Jfd, *pattern.submatrix, out_var, pattern.in_var.index()); + } + continue; + } + if (pattern.in_var.isAll()) { + for (int in_var = 0; in_var < nvars; ++in_var) { + addOperatorSparsity(Jfd, *pattern.submatrix, pattern.out_var.index(), in_var); + } + continue; + } + addOperatorSparsity(Jfd, *pattern.submatrix, pattern.out_var.index(), + pattern.in_var.index()); + } + + BOUT_DO_PETSC(MatAssemblyBegin(Jfd, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(Jfd, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); +} +#endif + /************************************************************************** * Constraints **************************************************************************/ diff --git a/tests/unit/solver/test_solver.cxx b/tests/unit/solver/test_solver.cxx index 525d52767f..9f8367b684 100644 --- a/tests/unit/solver/test_solver.cxx +++ b/tests/unit/solver/test_solver.cxx @@ -381,6 +381,39 @@ TEST_F(SolverTest, AddVector3D) { EXPECT_EQ(solver.listVector3DNames(), expected_names); } +TEST_F(SolverTest, GetVarRef) { + Options options; + FakeSolver solver{&options}; + + Field2D field2d{}; + Field3D field3d{}; + Vector3D vector{}; + + solver.add(field2d, "phi"); + solver.add(field3d, "n"); + solver.add(vector, "v"); + + const auto phi = solver.getVarRef("phi"); + const auto n = solver.getVarRef("n"); + const auto vx = solver.getVarRef("v_x"); + const auto vy = solver.getVarRef("v_y"); + const auto vz = solver.getVarRef("v_z"); + const auto missing = solver.getVarRef("missing"); + + EXPECT_TRUE(phi.isConcrete()); + EXPECT_EQ(phi.index(), 0); + EXPECT_TRUE(n.isConcrete()); + EXPECT_EQ(n.index(), 1); + EXPECT_TRUE(vx.isConcrete()); + EXPECT_EQ(vx.index(), 2); + EXPECT_TRUE(vy.isConcrete()); + EXPECT_EQ(vy.index(), 3); + EXPECT_TRUE(vz.isConcrete()); + EXPECT_EQ(vz.index(), 4); + EXPECT_TRUE(missing.isInvalid()); + EXPECT_TRUE(Solver::VarRef::All().isAll()); +} + TEST_F(SolverTest, ConstraintField2D) { Options options; FakeSolver solver{&options}; From a48cad0d00062c3e048d52023e3e5bdbf547dcf5 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 1 Jul 2026 16:33:42 -0700 Subject: [PATCH 17/18] Documentation of solvers and preconditioning Some reorganisation, removing duplicate and obsolete sections. --- manual/sphinx/index.rst | 2 +- .../sphinx/user_docs/parallel_operators.rst | 82 +++ manual/sphinx/user_docs/preconditioning.rst | 467 ++++++++---------- manual/sphinx/user_docs/time_integration.rst | 392 ++++----------- 4 files changed, 376 insertions(+), 567 deletions(-) diff --git a/manual/sphinx/index.rst b/manual/sphinx/index.rst index 1d4aaeab10..b63d5c1236 100644 --- a/manual/sphinx/index.rst +++ b/manual/sphinx/index.rst @@ -67,6 +67,7 @@ The documentation is divided into the following sections: :name: bout-interfaces user_docs/time_integration + user_docs/preconditioning user_docs/parallel-transforms user_docs/differential_operators user_docs/parallel_operators @@ -92,7 +93,6 @@ The documentation is divided into the following sections: :name: coordinates user_docs/coordinates - user_docs/preconditioning user_docs/BOUT_Gradperp_op .. toctree:: diff --git a/manual/sphinx/user_docs/parallel_operators.rst b/manual/sphinx/user_docs/parallel_operators.rst index b44e60ad70..23f4c73417 100644 --- a/manual/sphinx/user_docs/parallel_operators.rst +++ b/manual/sphinx/user_docs/parallel_operators.rst @@ -292,6 +292,86 @@ Current limitations include: - PETSc and 3D metrics are required - The grid metadata must be generated ahead of time +Using operator stencils in solver Jacobians +------------------------------------------- + +The ``PetscCellOperator`` objects can also be used to augment the Jacobian +pattern of PETSc-backed implicit solvers. This is useful when the default +solver stencil is too small for the operators used in the model, but the +coupling still has a sparse cell-to-cell structure. + +At a high level: + +- Build or reuse a ``PetscCellOperator`` +- Get solver variable references with ``solver->getVarRef(...)`` +- Register the operator with ``solver->addJacobianPattern(...)`` before + ``solver->init()`` + +For example: + +.. code-block:: C++ + + #if BOUT_HAS_PETSC + PetscOperators ops(mesh); + auto parallel = ops.getParallel(); + + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); + + // Apply the operator stencil to every Jacobian block + solver->addJacobianPattern(parallel.Div_par_Grad_par); + + // Apply the same cell stencil only to dF_n / dT + solver->addJacobianPattern(parallel.Grad_par, n, T); + #endif + +The one-argument overload is equivalent to: + +.. code-block:: C++ + + solver->addJacobianPattern(op, Solver::VarRef::All(), Solver::VarRef::All()); + +The ``VarRef::All()`` sentinel expands when the solver Jacobian is created: + +- ``(All, All)`` inserts the stencil into every variable block +- ``(out, All)`` inserts it into one output-variable row block against all inputs +- ``(All, in)`` inserts it into all output-variable row blocks against one input +- ``(out, in)`` inserts it into a single block + +This interface is only available when BOUT++ is built with PETSc. Even then, +``addJacobianPattern(...)`` may return ``false`` if the chosen solver does not +use the PETSc preconditioner/Jacobian path. See :ref:`sec-time-integration` and +:ref:`sec-preconditioning` for the solver-side behavior. + + +How solver registration works +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +``PetscOperator`` tracks its row and column mappings, and exposes them through +the operator API. For Jacobian-pattern registration, BOUT++ uses the operator's +cell mapping to restrict the full cell-space matrix to the evolving subset of +cells that actually appear in the solver state vector. + +The resulting evolving-cell submatrix is then inserted into one or more +Jacobian variable blocks using the helper documented in +``bout/petsc_jacobian.hxx``. This means the operator contributes only sparsity +information here; the Jacobian entries are still computed later by the solver's +finite-difference machinery. + +Current limitations +~~~~~~~~~~~~~~~~~~~ + +This solver-Jacobian path currently assumes the Jacobian uses a uniform +per-cell interleaving of the evolving variables. In practice that means it is +currently intended for the same cases documented in +``bout/petsc_jacobian.hxx``: + +- Evolving ``Field3D`` variables only +- No evolving boundary cells in the solver state +- No solver layouts with mixed or non-uniform cell-to-row mappings + +If those assumptions do not hold, use the standard solver coloring stencil or a +different Jacobian/preconditioning strategy instead. See also -------- @@ -299,4 +379,6 @@ See also - :ref:`sec-diffops` - :ref:`sec-fci` - :ref:`sec-parallel-bc-fci` +- :ref:`sec-time-integration` +- :ref:`sec-preconditioning` - :ref:`sec-petsc` diff --git a/manual/sphinx/user_docs/preconditioning.rst b/manual/sphinx/user_docs/preconditioning.rst index 79d8e050ae..78a25814bd 100644 --- a/manual/sphinx/user_docs/preconditioning.rst +++ b/manual/sphinx/user_docs/preconditioning.rst @@ -1,351 +1,284 @@ .. default-role:: math -====================== -BOUT++ preconditioning -====================== +.. _sec-preconditioning: -:Author: B.Dudson, University of York +===================================== +BOUT++ preconditioning and Jacobians +===================================== -Introduction -============ +Implicit solvers repeatedly solve linearised systems of the form +`(I - \gamma J) x = b`. BOUT++ supports several ways to help those solves: -This manual describes some of the ways BOUT++ could (and in some cases -does) support preconditioning, Jacobian calculations and other methods -to speed up simulations. This manual assumes that you’re familiar with -how BOUT++ works internally. +- a user-supplied preconditioner registered with ``solver->setPrecon(...)`` +- a user-supplied Jacobian-vector product registered with + ``solver->setJacobian(...)`` +- PETSc finite-difference Jacobians, usually accelerated with coloring +- extra Jacobian sparsity supplied from ``PetscCellOperator`` stencils with + ``solver->addJacobianPattern(...)`` -Some notation: The ODE being solved is of the form +This page starts with the user-facing workflow and then moves into the current +implementation details and limitations. -.. math:: {\frac{\partial {\mathbf{f}}}{\partial t}} = {\mathbf{F}}\left({\mathbf{f}}\right) -Here the state vector `f = \left(f_0, f_1, f_2, \ldots\right)^T` -is a vector containing the evolving (3D) variables -`f_i\left(x,y,z\right)`. +Choosing an approach +-------------------- -The Jacobian of this system is then +The most useful question is usually not "which feature exists?" but "which +approximation is easiest to provide for this model?". -.. math:: {\mathbb{J}}= {\frac{\partial {\mathbf{F}}}{\partial {\mathbf{f}}}} +- Use matrix-free methods when assembling a Jacobian is too expensive or the + coupling is not well represented by a sparse local stencil. +- Add a user preconditioner when you know an approximate inverse for the stiff + part of the physics. +- Use PETSc finite-difference Jacobians when the coupling is sparse and local + enough that a structured sparsity pattern is a good approximation. +- Add ``PetscCellOperator`` sparsity when the solver's default stencil is too + small, but the coupling is still naturally represented by a sparse + cell-to-cell operator. -The order of the elements in the vector `{\mathbf{f}}` -is determined in the solver code and SUNDIALS, so here just assume that -there exists a map `\mathbb{I}` between a global index `k` -and (variable, position) i.e. `\left(i,x,y,z\right)` +The solver options that choose between these modes are described in +:ref:`sec-time-integration`. PETSc option prefixes are described in +:ref:`sec-petsc`. -.. math:: \mathbf{I} : \left(i,x,y,z\right) \mapsto k -and its inverse +User-supplied preconditioners +----------------------------- -.. math:: \mathbf{I}^{-1} : k \mapsto \left(i,x,y,z\right) +A user preconditioner is a callback that approximately solves +`(I - \gamma J) x = b` for the solver's current linearisation. The callback +signature is: -Some problem-specific operations which can be used to speed up the -timestepping +.. code-block:: C++ -#. Jacobian-vector multiply: Given a vector, multiply it by - `{\mathbb{J}}` + int precon(BoutReal t, BoutReal gamma, BoutReal delta); -#. Preconditioner multiply: Given a vector, multiply by an approximate - inverse of `\mathbb{M} = \mathbb{I} - \gamma\mathbb{J}` +Register it during ``PhysicsModel::init``: -#. Calculate the stencils i.e. non-zero elements in - `{\mathbb{J}}` +.. code-block:: C++ -#. Calculate the non-zero elements of `{\mathbb{J}}` + int init(bool restarting) override { + solver->setPrecon(precon); + return 0; + } -Physics problems -================ +Inside the callback: -Some interesting physics problems of increasing difficulty +- the current state is stored in the usual evolving variables +- the vector to be preconditioned is stored in the corresponding + ``ddt(variable)`` fields +- the callback should overwrite those ``ddt(...)`` fields with the + preconditioned result -Resistive drift-interchange instability ---------------------------------------- +For CVODE, enable a user preconditioner with for example: -A “simple” test problem of 2 fields, which results in non-trivial -turbulent results. Supports resistive drift wave and interchange -instabilities. +.. code-block:: cfg -.. math:: - - \begin{aligned} - {\frac{\partial N_i}{\partial t}} + {{\mathbf{v}}_E}\cdot\nabla N_i &=& 0 \\ - {\frac{\partial \omega}{\partial t}} + {{\mathbf{v}}_E}\cdot\nabla\omega &=& 2\omega_{ci}{\mathbf{b}}\times\kappa\cdot\nabla P + N_iZ_i e\frac{4\pi V_A^2}{c^2}\nabla_{||}j_{||} \\ - \nabla_\perp^2\omega / N_i &=& \phi \\ - 0.51\nu_{ei}j_{||} &=& \frac{e}{m_e}\partial_{||}\phi + \frac{T_e}{N_i m_e}\partial_{||} N_i\end{aligned} - -Reduced 3-field MHD -------------------- - -This is a 3-field system of pressure `P`, magnetic flux -`\psi` and vorticity `U`: + [solver] + type = cvode + cvode_precon_method = user -.. math:: - - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U - \end{array}\right) +The ``examples/test-precon`` case shows the standard pattern. For the simple +wave system .. math:: - \begin{aligned} - {\frac{\partial \psi}{\partial t}} &=& -\frac{1}{B_0}\nabla_{||}\phi \\ - &=& -\frac{1}{B_0}\left[{\mathbf{b}}_0 - \left({\mathbf{b}}_0\times\nabla\psi\right)\right]\cdot\nabla\phi \\ - &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla\phi - \frac{1}{B_0}\left({\mathbf{b}}_0\times\nabla\phi\right)\cdot\nabla\psi \\ - \Rightarrow \frac{d \psi}{dt} &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla \phi\end{aligned} + \frac{\partial u}{\partial t} = \partial_{||} v + \qquad + \frac{\partial v}{\partial t} = \partial_{||} u -The coupled set of equations to be solved are therefore +one useful approximation factors `(I - \gamma J)^{-1}` into inexpensive pieces: .. math:: - \begin{aligned} - \frac{1}{B_0}\nabla_\perp^2\phi &=& U \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)\psi &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla\phi \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)P &=& 0 \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)U &=& \frac{1}{\rho}B_0^2\left[{\mathbf{b}}_0 - \left({\mathbf{b}}_0\times\nabla\psi\right)\right]\cdot\left(\frac{J_{||0}}{B_0} - \frac{1}{\mu_0}\nabla_\perp^2\psi\right) \nonumber \\ - &+& \frac{1}{\rho}{\mathbf{b}}_0\times{\mathbf{\kappa}}_0\cdot\nabla P \\ - {\mathbf{v}}_E &=& \frac{1}{B_0}{\mathbf{b}}_0\times\nabla\phi\end{aligned} + \left(\begin{array}{cc} + 1 & -\gamma \partial_{||} \\ + -\gamma \partial_{||} & 1 + \end{array}\right)^{-1} + = + \left(\begin{array}{cc} + 1 & \gamma \partial_{||} \\ + 0 & 1 + \end{array}\right) + \left(\begin{array}{cc} + 1 & 0 \\ + 0 & (1 - \gamma^2 \partial_{||}^2)^{-1} + \end{array}\right) + \left(\begin{array}{cc} + 1 & 0 \\ + \gamma \partial_{||} & 1 + \end{array}\right) -The Jacobian of this system is therefore: +In BOUT++ that becomes a sequence of field operations and an ``InvertPar`` solve: -.. math:: - :label: eq:mhdjacobian +.. code-block:: C++ - \mathbb{J} = - \left[ \begin{array}{c|c|c} - \color{blue}{-{\mathbf{v}}_E\cdot\nabla} & 0 & \left[{\mathbf{b}}_0\times\nabla\left(P_0 + \color{blue}{P}\right)\cdot\nabla\right]\nabla_\perp^{-2} \\ - \hline - 0 & \color{blue}{-{\mathbf{v}}_E\cdot\nabla} & \left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^{-2} \\ - \hline - 2{\mathbf{b}}_0\times{\mathbf{\kappa}}_0\cdot\nabla& -\frac{B_0^2}{\mu_0\rho}\left({\mathbf{b}}_0 \color{blue}{-{\mathbf{b}}_0\times\nabla\psi}\right)\cdot\nabla\nabla_\perp^2& \color{blue}{-{\mathbf{v}}_E\cdot\nabla} \\ - & + \frac{B_0^2}{\rho}\left[{\mathbf{b}}_0\times\nabla\left(\frac{J_{||0}}{B_0}\right)\right]\cdot\nabla & \\ - & + \color{blue}{\frac{B_0^2}{\mu_0\rho}\nabla\left(\nabla_\perp^2\psi\right)\cdot\left({\mathbf{b}}_0\times\nabla\right)} & - \end{array}\right] + int precon(BoutReal t, BoutReal gamma, BoutReal delta) { + mesh->communicate(ddt(u)); + ddt(v) = gamma * Grad_par(ddt(u)) + ddt(v); -Where the blue terms are only included in nonlinear simulations. + inv->setCoefB(-SQ(gamma)); + ddt(v) = inv->solve(ddt(v)); -This Jacobian has large dense blocks because of the Laplacian inversion -terms (involving `\nabla_\perp^{-2}` which couples together all -points in an X-Z plane. The way to make `{\mathbb{J}}` -sparse is to solve `\phi` as a constraint (using e.g. the IDA -solver) which moves the Laplacian inversion to the preconditioner. + mesh->communicate(ddt(v)); + ddt(u) = ddt(u) + gamma * Grad_par(ddt(v)); -Solving `\phi` as a constraint ------------------------------------- + ddt(u).applyBoundary("dirichlet"); + ddt(v).applyBoundary("dirichlet"); + return 0; + } -The evolving state vector becomes +This is the main design pattern for analytic preconditioners: keep only the +stiff physics, approximate the rest, and return something that is cheap but +captures the dominant fast scales. -.. math:: - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U \\ - \phi - \end{array}\right) +Jacobian-vector products +------------------------ -UEDGE equations ---------------- +If the model can apply the Jacobian to a vector directly, register that +operation with ``solver->setJacobian(...)``. This is primarily useful for +matrix-free methods, where the solver needs Jacobian-vector products but does +not assemble a full sparse Jacobian matrix. -The UEDGE benchmark is a 4-field model with the following equations: +For problems where a sparse matrix is more useful than a matrix-free operator, +PETSc-backed solvers can instead assemble the Jacobian by finite differences. -.. math:: - \begin{aligned} - {\frac{\partial N_i}{\partial t}} + {V_{||}}\partial_{||}N_i &=& -N_i\nabla_{||}{V_{||}}+\nabla_\psi\left(D_\perp \partial_\psi N_i\right) \\ - {\frac{\partial \left(N_i{V_{||}}\right)}{\partial t}} + {V_{||}}\partial_{||}\left(N_i{V_{||}}\right) &=& -\partial_{||}P + \nabla_\psi\left(N_i\mu_\perp\partial_\psi{V_{||}}\right) \\ - \frac{3}{2}{\frac{\partial }{\partial t}}\left(N_iT_e\right) &=& \nabla_{||}\left(\kappa_e\partial_{||}T_e\right) + \nabla_\psi\left(N_i\chi_\perp\partial_\perp T_e\right) \\ - \frac{3}{2}{\frac{\partial }{\partial t}}\left(N_iT_i\right) &=& \nabla_{||}\left(\kappa_i\partial_{||}T_i\right) + \nabla_\psi\left(N_i\chi_\perp\partial_\perp T_i\right)\end{aligned} +PETSc finite-difference Jacobians +--------------------------------- -This set of equations is good in that there is no inversion needed, and -so the Jacobian is sparse everywhere. The state vector is +When a PETSc-backed implicit solver runs with ``matrix_free = false``, BOUT++ +can assemble an explicit sparse Jacobian by finite differences. In most cases +``use_coloring = true`` is the right choice: -.. math:: +.. code-block:: cfg - {\mathbf{f}} = \left(\begin{array}{c} - N_i \\ - {V_{||}}\\ - T_e \\ - T_i \\ - \end{array}\right) + [solver] + matrix_free = false + use_coloring = true + lag_jacobian = 5 -The Jacobian is: + stencil:taxi = 2 + stencil:square = 0 + stencil:cross = 0 -.. math:: +The default pattern comes from a solver-side stencil and is then passed to PETSc +coloring, which groups independent columns so the Jacobian can be built with far +fewer RHS evaluations than a brute-force finite-difference calculation. - \mathbb{J} = - \left( \begin{array}{c|c|c|c} - -{V_{||}}\partial_{||} - \nabla_{||}{V_{||}}+ \nabla_\psi D_\perp\partial_\psi & -\partial_{||}N_i - N_i\nabla_{||} & 0 & 0 \\ - -\frac{1}{N_i}{\frac{\partial {V_{||}}}{\partial t}} - \frac{1}{N_i}{V_{||}}{\mathbb{J}}_{N_iN_i} & & & - \end{array}\right) +The most important current assumption is that the Jacobian is sparse and local +in the solver ordering. If the RHS includes long-range couplings such as Fourier +transforms, matrix inversions, or other nonlocal operators, the default stencil +can miss real dependencies and the implicit solve may fail to converge. In those +cases the usual alternatives are: -If instead the state vector is +- stay matrix-free +- reformulate the problem so the long-range solve becomes a constraint +- supply a better sparsity pattern when one exists -.. math:: - - {\mathbf{f}} = \left(\begin{array}{c} - N_i \\ - N_i{V_{||}}\\ - N_iT_e \\ - N_iT_i \\ - \end{array}\right) - -then the Jacobian is - -.. Result is missing! - -2-fluid turbulence ------------------- - -Jacobian-vector multiply -======================== - -This is currently implemented into the CVODE (SUNDIALS) solver. - -Preconditioner-vector multiply -============================== +The main tuning options are: -.. _reduced-3-field-mhd-1: +- ``use_coloring = false`` to fall back to a brute-force finite-difference + Jacobian +- ``stencil:taxi``, ``stencil:square``, and ``stencil:cross`` to widen the + default solver stencil +- ``force_symmetric_coloring = true`` to make the coloring pattern symmetric +- ``lag_jacobian`` to reuse the assembled Jacobian across successive nonlinear + iterations -Reduced 3-field MHD -------------------- -The matrix `\mathbb{M}` to be inverted can therefore be written +Augmenting Jacobian structure with PetscOperator +------------------------------------------------ -.. math:: +If the coupling you want to expose to the solver is already available as a +``PetscCellOperator``, you can register its sparsity pattern directly with the +solver before ``solver->init()``. - \mathbb{M} = - \left[ \begin{array}{ccc} - \mathbb{D} & 0 & \mathbb{U}_P \\ - 0 & \mathbb{D} & \mathbb{U}_\psi \\ - \mathbb{L}_P & \mathbb{L}_\psi & \mathbb{D} - \end{array}\right] +For example: -where +.. code-block:: C++ -.. math:: \mathbb{D} = \mathbb{I} \color{blue}{+ \gamma{\mathbf{v}}_E\cdot\nabla} + #if BOUT_HAS_PETSC + PetscOperators ops(mesh); + auto parallel = ops.getParallel(); -For small flow velocities, the inverse of `\mathbb{D}` can be -approximated using the Binomial theorem: + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); -.. math:: - :label: eq:dapprox + // Insert the operator stencil into every Jacobian block + solver->addJacobianPattern(parallel.Div_par_Grad_par); - \mathbb{D}^{-1} \simeq \mathbb{I} \color{blue}{- \gamma{\mathbf{v}}_E\cdot\nabla} + // Insert the operator stencil only into dF_n / dT + solver->addJacobianPattern(parallel.Grad_par, n, T); + #endif -Following [chacon-2008]_, [chacon-2002]_, `\mathbb{M}` can be -re-written as +The one-argument form is shorthand for: -.. math:: +.. code-block:: C++ - \mathbb{M} = - \left[ \begin{array}{cc} - \mathbb{E} & \mathbb{U} \\ - \mathbb{L} & \mathbb{D} - \end{array}\right] \qquad \mathbb{E} = - \left[ \begin{array}{cc} - \mathbb{D} & 0 \\ - 0 & \mathbb{D} - \end{array}\right] \qquad \mathbb{U} = - \left(\begin{array}{c} - \mathbb{U}_P \\ - \mathbb{U}_\psi - \end{array}\right) \qquad \mathbb{L} = \left(\mathbb{L}_P \quad \mathbb{L}_\psi\right) - -The Schur factorization of `\mathbb{M}` yields ([chacon-2008]_) + solver->addJacobianPattern(op, Solver::VarRef::All(), Solver::VarRef::All()); -.. math:: +``Solver::VarRef::All()`` expands when the solver Jacobian is created: - \mathbb{M}^{-1} = - \left[ \begin{array}{cc} - \mathbb{E} & \mathbb{U} \\ - \mathbb{L} & \mathbb{D} - \end{array}\right]^{-1} = - \left[ \begin{array}{cc} - \mathbb{I} & -\mathbb{E}^{-1}\mathbb{U} \\ - 0 & \mathbb{I} - \end{array}\right] - \left[ \begin{array}{cc} - \mathbb{E}^{-1} & 0 \\ - 0 & \mathbb{P}_{Schur}^{-1} - \end{array}\right] - \left[ \begin{array}{cc} - \mathbb{I} & 0 \\ - -\mathbb{L}\mathbb{E}^{-1} & \mathbb{I} - \end{array}\right] - -Where -`\mathbb{P}_{Schur} = \mathbb{D} - \mathbb{L}\mathbb{E}^{-1}\mathbb{U}` -is the Schur complement. Note that this inversion is exact so far. Since -`\mathbb{E}` is block-diagonal, and `\mathbb{D}` can be -easily approximated using equation :eq:`eq:dapprox`, this -simplifies the problem to inverting `\mathbb{P}_{Schur}`, which is -much smaller than `\mathbb{M}`. - -A possible approximation to `\mathbb{P}_{Schur}` is to neglect: - -- All drive terms - - - the curvature term `\mathbb{L}_P` - - - the `J_{||0}` term in `\mathbb{L}_\psi` - -- All nonlinear terms (blue terms in equation :eq:`eq:mhdjacobian`), - including perpendicular terms (so `\mathbb{D} = \mathbb{I}`) - -This gives +- ``(All, All)`` inserts the stencil into every variable block +- ``(out, All)`` inserts it into one output-variable row block against all inputs +- ``(All, in)`` inserts it into all output-variable row blocks against one input +- ``(out, in)`` inserts it into a single block -.. math:: +This interface is guarded by ``BOUT_HAS_PETSC``. Even in PETSc-enabled builds, +``addJacobianPattern(...)`` may return ``false`` if the chosen solver does not +use the PETSc preconditioner/Jacobian path. - \begin{aligned} - \mathbb{P}_{Schur} &\simeq& \mathbb{I} + \gamma^2 \frac{B_0^2}{\mu_0\rho}\left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^2\left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^{-2} \nonumber \\ - &\simeq& \mathbb{I} + \gamma^2 V_A^2 \left({\mathbf{b}}_0\cdot\nabla\right)^2\end{aligned} -Where the commutation of parallel and perpendicular derivatives is also -an approximation. This remaining term is just the shear Alfvén wave -propagating along field-lines, the fastest wave supported by these -equations. +How operator-supplied sparsity is applied +----------------------------------------- -Stencils -======== +The operator itself lives on the full PETSc cell space, including boundary and +virtual-cell entries used by the operator construction. The solver Jacobian does +not: it contains only the evolving degrees of freedom. -Jacobian calculation -==================== +To bridge those two layouts, BOUT++ uses the operator's stored cell mapping to +extract the evolving-cell submatrix first. That restricted submatrix is then +inserted into one or more solver Jacobian blocks using the helper documented in +``bout/petsc_jacobian.hxx``. -The (sparse) Jacobian matrix elements can be calculated automatically -from the physics code by keeping track of the (linearised) operations -going through the RHS function. +The important consequence is that ``addJacobianPattern(...)`` contributes only +the nonzero structure. The Jacobian entries themselves are still computed later +by the solver's finite-difference machinery. -For each point, keep the value (as usual), plus the non-zero elements in -that row of `{\mathbb{J}}` and the constant: result = -Ax + b Keep track of elements using product rule. +At the moment the sequence is: -:: +1. Model or component code calls ``solver->addJacobianPattern(...)`` during + setup. +2. The solver stores the restricted operator submatrix and the requested output + and input variable references. +3. During ``solver->init()``, a supported PETSc-backed solver creates its + default Jacobian pattern. +4. The queued operator patterns are replayed into that Jacobian. +5. PETSc coloring and finite-difference Jacobian assembly use the augmented + sparsity pattern. - class Field3D { - data[ngx][ngy][ngz]; // The data as now - - int JacIndex; // Variable index in Jacobian - SparseMatrix *jac; // Set of rows for indices (JacIndex,*,*,*) - }; -JacIndex is set by the solver, so for the system +Current limitations +------------------- -.. math:: +The helper used to insert operator sparsity into the solver Jacobian currently +assumes a uniform per-cell variable ordering of the form +``global_row = cell * nvars + var`` over evolving cells only. That means this +path is currently intended for: - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U - \end{array}\right) +- evolving ``Field3D`` variables +- no evolving boundary cells +- solver layouts that use a uniform per-cell interleaving of variables -``P.JacIndex = 0``, ``psi.JacIndex = 1`` and ``U.JacIndex = 2``. All -other fields are given ``JacIndex = -1``. +It is not currently valid for mixed ``Field2D`` and ``Field3D`` solver layouts, +for layouts that include evolving boundary variables, or for any other +non-uniform cell-to-row mapping. -SparseMatrix stores the non-zero Jacobian components for the set of rows -corresponding to this variable. Evolving variables do not have an -associated ``SparseMatrix`` object, but any fields which result from -operations on evolving fields will have one. -.. [chacon-2008] L. Chacón, An optimal, parallel, fully implicit Newton-Krylov solver for three-dimensional viscoresistive magnetohydrodynamics, POP, 2008, 15, 056103 +See also +-------- -.. [chacon-2002] L. Chacón, D.A. Knoll, and J.M. Finn, An Implicit, Nonlinear Reduced Resistive MHD Solver, JCP, 2002, 178, 15-36 +- :ref:`sec-time-integration` +- :ref:`sec-parallel-operators-petsc-fci` +- :ref:`sec-petsc` diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 097b62eeda..58ba7c4694 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -721,6 +721,14 @@ The default and recommended approach for most problems: The coloring algorithm exploits the sparse structure of the Jacobian to reduce the number of function evaluations needed for finite differencing. +If the default solver stencil is too narrow for part of the physics, PETSc +cell-space operators can add extra sparsity information before coloring is set +up. Register those contributions during model setup with +``solver->addJacobianPattern(...)``. This is particularly useful for +``PetscCellOperator`` stencils assembled from the FCI support-operator tools in +:ref:`sec-parallel-operators-petsc-fci`. See :ref:`sec-preconditioning` for the +full solver-side behavior and current limitations. + Jacobian coloring stencil ^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -997,263 +1005,42 @@ Finally, delete the model and solver when finished:: **Note:** If an ODE needs to be solved multiple times, at the moment it is recommended to delete the solver, and create a new one each time. -.. _sec-preconditioning: - Preconditioning --------------- -At every time step, an implicit scheme such as BDF has to solve a -non-linear problem to find the next solution. This is usually done using -Newton’s method, each step of which involves solving a linear (matrix) -problem. For :math:`N` evolving variables is an :math:`N\times N` matrix -and so can be very large. By default matrix-free methods are used, in -which the Jacobian :math:`\mathcal{J}` is approximated by finite -differences (see next subsection), and so this matrix never needs to be -explicitly calculated. Finding a solution to this matrix can still be -difficult, particularly as :math:`\delta t` gets large compared with -some time-scales in the system (i.e. a stiff problem). - -A preconditioner is a function which quickly finds an approximate -solution to this matrix, speeding up convergence to a solution. A -preconditioner does not need to include all the terms in the problem -being solved, as the preconditioner only affects the convergence rate -and not the final solution. A good preconditioner can therefore -concentrate on solving the parts of the problem with the fastest -time-scales. - -A simple example [1]_ is a coupled wave equation, solved in the -``test-precon`` example code: - -.. math:: - - \frac{\partial u}{\partial t} = \partial_{||}v \qquad \frac{\partial - v}{\partial t} = \partial_{||} u - -First, calculate the Jacobian of this set of equations by taking -partial derivatives of the time-derivatives with respect to each of the -evolving variables - -.. math:: - - \mathcal{J} = (\begin{array}{cc} - \frac{\partial}{\partial u}\frac{\partial u}{\partial t} & - \frac{\partial}{\partial v}\frac{\partial u}{\partial t}\\ - \frac{\partial}{\partial u}\frac{\partial v}{\partial t} & - \frac{\partial}{\partial v}\frac{\partial v}{\partial t} - \end{array} - ) = (\begin{array}{cc} - 0 & \partial_{||} \\ - \partial_{||} & 0 - \end{array} - ) - -In this case :math:`\frac{\partial u}{\partial t}` doesn’t depend on -:math:`u` nor :math:`\frac{\partial v}{\partial t}` on :math:`v`, so the -diagonal is empty. Since the equations are linear, the Jacobian doesn’t -depend on :math:`u` or :math:`v` and so - -.. math:: - - \frac{\partial}{\partial t}(\begin{array}{c} u \\ - v \end{array}) = \mathcal{J} (\begin{array}{c} u \\ - v \end{array} ) - -In general for non-linear functions :math:`\mathcal{J}` gives the -change in time-derivatives in response to changes in the state variables -:math:`u` and :math:`v`. - -In implicit time stepping, the preconditioner needs to solve an equation - -.. math:: - - \mathcal{I} - \gamma \mathcal{J} - -where :math:`\mathcal{I}` is the identity matrix, and :math:`\gamma` -depends on the time step and method (e.g. :math:`\gamma = \delta t` for -backwards Euler method). For the simple wave equation problem, this is - -.. math:: - - \mathcal{I} - \gamma \mathcal{J} = (\begin{array}{cc} - 1 & -\gamma\partial_{||} \\ - -\gamma\partial_{||} & 1 - \end{array} - ) - -This matrix can be block inverted using Schur factorisation [2]_ - -.. math:: +Implicit solvers repeatedly solve linearised systems of the form +:math:`(I - \gamma J)\,x = b`. Preconditioning improves the convergence of +those inner solves by applying a cheap approximate inverse, while +matrix-free or finite-difference Jacobian information controls how the solver +models the couplings in :math:`J`. - (\begin{array}{cc} - {\mathbf{E}} & {\mathbf{U}} \\ - {\mathbf{L}} & {\mathbf{D}} - \end{array})^{-1} - = (\begin{array}{cc} - {\mathbf{I}} & -{\mathbf{E}}^{-1}{\mathbf{U}} \\ - 0 & {\mathbf{I}} - \end{array} - )(\begin{array}{cc} - {\mathbf{E}}^{-1} & 0 \\ - 0 & {\mathbf{P}}_{Schur}^{-1} - \end{array} - )(\begin{array}{cc} - {\mathbf{I}} & 0 \\ - -{\mathbf{L}}{\mathbf{E}}^{-1} & {\mathbf{I}} - \end{array} - ) - -where -:math:`{\mathbf{P}}_{Schur} = {\mathbf{D}} - {\mathbf{L}}{\mathbf{E}}^{-1}{\mathbf{U}}` -Using this, the wave problem becomes: +For most models, the relevant entry points are: -.. math:: - :label: precon - - (\begin{array}{cc} 1 & -\gamma\partial_{||} \\ - -\gamma\partial_{||} & 1 \end{array})^{-1} = (\begin{array}{cc} 1 & \gamma\partial_{||}\\ - 0 & 1 \end{array} )(\begin{array}{cc} 1 & 0 \\ - 0 & (1 -\gamma^2\partial^2_{||})^{-1} \end{array} )(\begin{array}{cc} 1 & 0\\ - \gamma\partial_{||} & 1 \end{array} ) - -The preconditioner is implemented by defining a function of the form - -:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - ... - } - -which takes as input the current time, the :math:`\gamma` factor -appearing above, and :math:`\delta` which is only important for -constrained problems (not discussed here... yet). The current state of -the system is stored in the state variables (here ``u`` and ``v`` ), -whilst the vector to be preconditioned is stored in the time derivatives -(here ``ddt(u)`` and ``ddt(v)`` ). At the end of the preconditioner the -result should be in the time derivatives. A preconditioner which is just -the identity matrix and so does nothing is therefore:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - } - -To implement the preconditioner in equation :eq:`precon`, first apply the -rightmost matrix to the given vector: - -.. math:: +- ``solver->setPrecon(...)`` for a user-supplied preconditioner +- ``solver->setJacobian(...)`` for a Jacobian-vector product +- ``matrix_free = false`` and ``use_coloring = true`` for PETSc-backed + finite-difference Jacobians +- ``solver->addJacobianPattern(...)`` when a ``PetscCellOperator`` provides a + better sparsity pattern than the default solver stencil - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) = (\begin{array}{cc} - 1 & 0 \\ - \gamma\partial_{||} & 1 - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) +The detailed discussion, worked example, and implementation notes are in +:ref:`sec-preconditioning`. This section keeps the focus on solver options and +setup. -:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - mesh->communicate(ddt(u)); - //ddt(u) = ddt(u); - ddt(v) = gamma*Grad_par(ddt(u)) + ddt(v); - -note that since the preconditioner is linear, it doesn’t depend on -:math:`u` or :math:`v`. As in the RHS function, since we are taking a -differential of ``ddt(u)``, it first needs to be communicated to -exchange guard cell values. - -The second matrix - -.. math:: - - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) \rightarrow (\begin{array}{cc} - 1 & 0 \\ - 0 & (1 - \gamma^2\partial^2_{||})^{-1} - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) - -doesn’t alter :math:`u`, but solves a parabolic equation in the -parallel direction. There is a solver class to do this called -`InvertPar` which solves the equation :math:`(A + B\partial_{||}^2)x = -b` where :math:`A` and :math:`B` are `Field2D` or constants [3]_. In -`PhysicsModel::init` we create one of these solvers:: - - InvertPar *inv; // Parallel inversion class - int init(bool restarting) { - ... - inv = InvertPar::Create(); - inv->setCoefA(1.0); - ... - } - -In the preconditioner we then use this solver to update :math:`v`:: - - inv->setCoefB(-SQ(gamma)); - ddt(v) = inv->solve(ddt(v)); - -which solves -:math:`ddt(v) \rightarrow (1 - \gamma^2\partial_{||}^2)^{-1} ddt(v)`. -The final matrix just updates :math:`u` using this new solution for -:math:`v` - -.. math:: - - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) \rightarrow (\begin{array}{cc} - 1 & \gamma\partial_{||} \\ - 0 & 1 - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) - -:: - - mesh->communicate(ddt(v)); - ddt(u) = ddt(u) + gamma*Grad_par(ddt(v)); - -Finally, boundary conditions need to be imposed, which should be -consistent with the conditions used in the RHS:: - - ddt(u).applyBoundary("dirichlet"); - ddt(v).applyBoundary("dirichlet"); - -To use the preconditioner, pass the function to the solver in -`PhysicsModel::init`:: - - int init(bool restarting) { - solver->setPrecon(precon); - ... - } - -then in the ``BOUT.inp`` settings file switch on the preconditioner - -.. code-block:: bash +Jacobian-vector function +------------------------ - [solver] - type = cvode # Need CVODE or PETSc - cvode_precon_method = user # Use user-supplied preconditioner - rightprec = false # Use Right preconditioner (default left) +If the physics model can apply the Jacobian to a vector directly, it can provide +that routine through ``solver->setJacobian(...)`` during ``PhysicsModel::init``. +This is mainly useful in matrix-free workflows, where the solver needs +Jacobian-vector products but does not assemble a full sparse matrix. -Jacobian function ------------------ +When a PETSc-backed implicit solver is configured with +``matrix_free = false``, BOUT++ can instead assemble an explicit sparse +Jacobian by finite differences. In that case, the coloring options described +below determine how the sparsity pattern is used, and +``solver->addJacobianPattern(...)`` can add extra block structure from +``PetscCellOperator`` stencils before the colored finite differencing is set up. DAE constraint equations ------------------------ @@ -1458,73 +1245,80 @@ here:\ https://computation.llnl.gov/casc/sundials/support/notes.html). This may in some cases be less efficient. -Implementation internals ------------------------- +Solver interface +---------------- -The solver is the interface between BOUT++ and the time-integration -code such as SUNDIALS. All solvers implement the `Solver` -class interface (see ``src/solver/generic_solver.hxx``). +The solver is the interface between BOUT++ and the time-integration library. +All time integrators implement the ``Solver`` interface, and most user-facing +setup happens during ``PhysicsModel::init`` before ``solver->init()`` is called. -First all the fields which are to be evolved need to be added to the -solver. These are always done in pairs, the first specifying the field, -and the second the time-derivative:: +The main setup methods are: - void add(Field2D &v, Field2D &F_v, const char* name); +.. code-block:: C++ -This is normally called in the `PhysicsModel::init` initialisation routine. -Some solvers (e.g. IDA) can support constraints, which need to be added -in the same way as evolving fields:: + void add(Field2D& v, const std::string& name, + const std::string& description = ""); + void add(Field3D& v, const std::string& name, + const std::string& description = ""); + void add(Vector2D& v, const std::string& name, + const std::string& description = ""); + void add(Vector3D& v, const std::string& name, + const std::string& description = ""); - bool constraints(); - void constraint(Field2D &v, Field2D &C_v, const char* name); + bool constraints(); + void constraint(Field3D& v, Field3D& residual, std::string name); -The ``constraints()`` function tests whether or not the current solver -supports constraints. The format of ``constraint(...)`` is the same as -``add``, except that now the solver will attempt to make ``C_v`` zero. -If ``constraint`` is called when the solver doesn’t support them then an -error should occur. + typedef int (*PhysicsPrecon)(BoutReal t, BoutReal gamma, BoutReal delta); + void setPrecon(PhysicsPrecon f); -If the physics model implements a preconditioner or Jacobian-vector -multiplication routine, these can be passed to the solver during -initialisation:: + typedef int (*Jacobian)(BoutReal t); + void setJacobian(Jacobian j); - typedef int (*PhysicsPrecon)(BoutReal t, BoutReal gamma, BoutReal delta); - void setPrecon(PhysicsPrecon f); // Specify a preconditioner - typedef int (*Jacobian)(BoutReal t); - void setJacobian(Jacobian j); // Specify a Jacobian +#if BOUT_HAS_PETSC + Solver::VarRef getVarRef(std::string_view name) const; + bool addJacobianPattern(const PetscCellOperator& op); + bool addJacobianPattern(const PetscCellOperator& op, Solver::VarRef out_var, + Solver::VarRef in_var); +#endif -If the solver doesn’t support these functions then the calls will just -be ignored. +``add(...)`` registers evolving variables. ``constraint(...)`` registers DAE +constraints for solvers that support them. ``setPrecon(...)`` and +``setJacobian(...)`` register optional callbacks used by implicit solvers. -Once the problem to be solved has been specified, the solver can be -initialised using:: +``getVarRef(...)`` resolves the solver's internal variable numbering from the +name used when the variable or component was added. Those references can then be +used with ``addJacobianPattern(...)`` to target one Jacobian block, one row or +column of blocks, or all blocks via ``Solver::VarRef::All()``. - int init(); +For example: -which returns an error code (0 on success). This is currently called in -:doc:`bout++.cxx<../_breathe_autogen/file/bout_09_09_8cxx>`:: +.. code-block:: C++ - if (solver.init()) { - output.write("Failed to initialise solver. Aborting\n"); - return(1); - } + #if BOUT_HAS_PETSC + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); -which passes the (physics module) RHS function `PhysicsModel::rhs` to the -solver along with the number and size of the output steps. + solver->addJacobianPattern(op); + solver->addJacobianPattern(op, n, T); + #endif -:: +These Jacobian-pattern registrations are designed to be called during model or +component setup, before the solver knows the final matrix size. BOUT++ therefore +queues them and applies them later when a supported PETSc-backed solver creates +its Jacobian structure. + +Behavior depends on the solver: - typedef int (*MonitorFunc)(BoutReal simtime, int iter, int NOUT); - int run(MonitorFunc f); +- ``setPrecon(...)`` and ``setJacobian(...)`` are ignored by solvers that do + not use those callbacks +- ``addJacobianPattern(...)`` is only available in PETSc-enabled builds +- even with PETSc, ``addJacobianPattern(...)`` may return ``false`` if the + chosen solver does not use the PETSc preconditioner/Jacobian path -.. [1] - Taken from a talk by L.Chacon available here - https://bout2011.llnl.gov/pdf/talks/Chacon_bout2011.pdf +Once the problem has been specified, the solver is initialised with:: + + int init(); -.. [2] - See paper https://arxiv.org/abs/1209.2054 for an application to - 2-fluid equations +and then run in the usual way:: -.. [3] This `InvertPar` class can handle cases with closed - field-lines and twist-shift boundary conditions for tokamak - simulations + int run(); From 2749c0e7d8c2feb18b0c11282c34ea21ee97a028 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 3 Jul 2026 08:22:35 -0700 Subject: [PATCH 18/18] Clang-tidy suggestions --- include/bout/solver.hxx | 23 +++++++---------------- src/mesh/petsc_operators.cxx | 1 + src/solver/impls/cvode/cvode.hxx | 2 +- src/solver/impls/petsc/petsc.hxx | 6 ++++-- src/solver/impls/snes/snes.hxx | 3 +-- src/solver/solver.cxx | 3 +++ 6 files changed, 17 insertions(+), 21 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 0a5a53ca81..bfaf73b7d8 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -1,20 +1,10 @@ /************************************************************************** * Base class for all solvers. Specifies required interface functions * - * Changelog: - * - * 2009-08 Ben Dudson, Sean Farley - * * Major overhaul, and changed API. Trying to make consistent - * interface to PETSc and SUNDIALS solvers - * - * 2013-08 Ben Dudson - * * Added OO-style API, to allow multiple physics models to coexist - * For now both APIs are supported - * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2026 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -45,7 +35,6 @@ #include "bout/monitor.hxx" #include "bout/options.hxx" #include "bout/region.hxx" -#include "bout/unused.hxx" #include #include @@ -86,6 +75,7 @@ using TimestepMonitorFunc = int (*)(Solver* solver, BoutReal simtime, BoutReal l #include #include #include +#include using SolverType = std::string; constexpr auto SOLVERCVODE = "cvode"; @@ -215,6 +205,7 @@ using RegisterUnavailableSolver = SolverFactory::RegisterUnavailableInFactory; */ class Solver { public: + /// Variable reference handle class VarRef { public: static constexpr int AllValue = -1; @@ -437,7 +428,7 @@ protected: bool covariant{false}; /// For vectors bool evolve_bndry{false}; /// Are the boundary regions being evolved? std::string name; /// Name of the variable - std::string description{""}; /// Description of what the variable is + std::string description; /// Description of what the variable is }; /// A structure for iterating over fields @@ -735,8 +726,8 @@ protected: #if BOUT_HAS_PETSC struct DeferredJacobianPattern { bout::petsc::UniqueMat submatrix{new Mat{nullptr}}; - VarRef out_var{}; - VarRef in_var{}; + VarRef out_var; + VarRef in_var; }; /// Queue a Jacobian-pattern contribution for PETSc-preconditioner-based solvers. diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index be583b675c..581e02eb24 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -17,6 +17,7 @@ #include #include +#include #include #include #include diff --git a/src/solver/impls/cvode/cvode.hxx b/src/solver/impls/cvode/cvode.hxx index 3e8c116718..bdf601cfd2 100644 --- a/src/solver/impls/cvode/cvode.hxx +++ b/src/solver/impls/cvode/cvode.hxx @@ -44,10 +44,10 @@ RegisterUnavailableSolver #include "../../sundials_nvector_interface.hxx" #include "bout/bout_types.hxx" -#include "bout/region.hxx" #include "bout/sundials_backports.hxx" #if BOUT_HAS_PETSC +#include "bout/petsc_operators.hxx" #include "bout/petsc_preconditioner.hxx" #include "bout/petsclib.hxx" diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index 2438f7de1c..fea8a03411 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -2,7 +2,7 @@ * Interface to PETSc solver * ************************************************************************** - * Copyright 2010 - 2025 BOUT++ contributors + * Copyright 2010 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -27,6 +27,7 @@ #define BOUT_PETSC_SOLVER_H #include "bout/build_defines.hxx" +#include "bout/petsc_operators.hxx" #include "bout/solver.hxx" #if not BOUT_HAS_PETSC @@ -40,6 +41,7 @@ RegisterUnavailableSolver class PetscSolver; +#include #include #include #include @@ -53,7 +55,7 @@ class PetscSolver; #include #include -#include +#include namespace { RegisterSolver registersolverpetsc("petsc"); diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 21f1d72f01..4c749f75e0 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -82,12 +82,11 @@ public: int init() override; int run() override; -#if BOUT_HAS_PETSC + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, VarRef in_var) override { return queueJacobianPattern(op, out_var, in_var); } -#endif /// Nonlinear function. This is called by PETSc SNES object /// via a static C-style function. For implicit diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index f6c032268e..b18f30c045 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -49,10 +49,13 @@ #if BOUT_HAS_PETSC #include "bout/petsc_jacobian.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/petsclib.hxx" #endif #include +#include #include #include #include