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..73d4d2fc71 --- /dev/null +++ b/include/bout/petsc_jacobian.hxx @@ -0,0 +1,43 @@ +#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); + +/// @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); + +#endif // BOUT_HAS_PETSC + +#endif // BOUT_PETSC_JACOBIAN_H diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 0d9c6c6912..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 @@ -8,7 +13,6 @@ #include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/mesh.hxx" -#include "bout/output.hxx" #include "bout/output_bout_types.hxx" #include "bout/petsc_interface.hxx" #include "bout/petsclib.hxx" @@ -45,6 +49,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. /// @@ -275,6 +283,37 @@ 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; + + /// @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 @@ -865,8 +904,8 @@ public: /// /// Naming convention: the subscript on Grad/Div/Restrict refers to which side of /// the cell face the operator acts on, not the leg space it lives in. - /// - @c Grad_plus / @c Div_plus: forward-side half-step operators (in L-). - /// - @c Grad_minus / @c Div_minus: backward-side half-step operators (in L+). + /// - @c Grad_plus / @c Div_plus: forward-side half-step operators (C to L+). + /// - @c Grad_minus / @c Div_minus: backward-side half-step operators (C to L-). /// - @c Restrict_minus = I+^T * W+: weighted restriction from L+ back to C, /// paired with the minus-side gradient. /// - @c Restrict_plus = I-^T * W-: weighted restriction from L- back to C, @@ -914,7 +953,7 @@ public: /// /// Evaluates the half-step fluxes separately on each side, multiplies by the /// interpolated diffusion coefficient K, and averages the two divergence - /// contributions. This is the primary user-facing SOM operator. + /// contributions. /// /// @param K Diffusion coefficient field; must have parallel slices allocated. /// @param f Field to differentiate; must have parallel slices allocated. @@ -1021,4 +1060,6 @@ private: #warning PETSc not available. No PetscOperators. -#endif +#endif // BOUT_HAS_PETSC + +#endif //BOUT_PETSC_OPERATORS_H diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx new file mode 100644 index 0000000000..457c667de5 --- /dev/null +++ b/src/mesh/petsc_jacobian.cxx @@ -0,0 +1,65 @@ +#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) { + 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 diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 0fe564bcd3..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" @@ -127,6 +126,33 @@ 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; +} + +Mat PetscCellMapping::extractEvolvingSubmatrix( + const PetscOperator& op) const { + 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) { 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/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 diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index d3d690417b..aca8b833f1 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -98,6 +98,9 @@ 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 + ./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..ac2e8e5d08 --- /dev/null +++ b/tests/unit/mesh/test_petsc_jacobian.cxx @@ -0,0 +1,638 @@ +// ============================================================================ +// 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 +#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); +} + +// ============================================================================ +// 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 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 new file mode 100644 index 0000000000..4ef74318f2 --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx @@ -0,0 +1,489 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "bout/array.hxx" +#include "bout/bout_types.hxx" +#include "bout/globals.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 new file mode 100644 index 0000000000..9dbaed584e --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx @@ -0,0 +1,484 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#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