From 87f6feb331adeed9fcc171a295ee363a509e52d9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 16 Jul 2026 12:03:17 -0700 Subject: [PATCH 1/7] Reapply "Merge pull request #3082 from boutproject/next-elimBE" This reverts commit 81befc23a51139f22c8b2835cacebd16dfdd5457. --- CMakeLists.txt | 1 + include/bout/hypre_interface.hxx | 185 +++++++++++++++++- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 5 + src/sys/hypre_interface.cxx | 131 +++++++++++++ 4 files changed, 320 insertions(+), 2 deletions(-) create mode 100644 src/sys/hypre_interface.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index e6e424b55c..6b1c4290db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -371,6 +371,7 @@ set(BOUT_SOURCES ./src/sys/generator_context.cxx ./include/bout/hyprelib.hxx ./src/sys/hyprelib.cxx + ./src/sys/hypre_interface.cxx ./src/sys/msg_stack.cxx ./src/sys/options.cxx ./src/sys/options/optionparser.hxx diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 189e91d159..f9f97b5bb3 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -58,6 +58,110 @@ int checkHypreError(int error) { // TODO: set sizes // TODO: set contiguous blocks at once +/// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed. +struct HypreComplexArray { + HYPRE_Complex* data; + + HypreComplexArray(int n) { HypreMalloc(data, sizeof(HYPRE_Complex) * n); } + + ~HypreComplexArray() { HypreFree(data); } +}; + +/// Shared pointter to a HypreComplexArray. When the last copy is destroyed +/// the HYPRE_Complex array inside will be free'd. +using BCValuesPtr = std::shared_ptr; + +/*! + * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to + * eliminate the boundary condition equations (see below for details on how the + * equations are adjusted). It modifies the arrays ncols, rows, cols, and + * values. It also returns a row_indexes array. This can then be passed to the + * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre. + * + * The arguments nb and bi_array indicate the boundary equations. The routine + * returns info needed to adjust the right-hand-side and solution vector through + * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. + * + * NOTE: It may make sense from an organizational standpoint to collect many of + * these arguments in a structure of some sort. + * + * Notation, assumptions, and other details: + * + * - Boundary equation i is assumed to have two coefficients + * + * b_ii * u_i + b_ij * u_j = rhs_i + * + * - We also assume that each boundary equation has only one interior equation k + * coupled to it (such that k = j) with coupling coefficient a_ki + * + * a_ki * u_i + a_kj * u_j + ... = rhs_k + * + * - Each equation k is adjusted as follows: + * + * a_kj = a_kj - a_ki * b_ij / b_ii + * a_ki = 0 + * + * - Boundary equations are adjusted to be identity equations in the matrix, but + * the boundary coefficients (b_ii, b_ij) are returned for use later + * + * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as + * follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii + * + * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as + * follows: u_i = (rhs_i - b_ij * u_j) / b_ii + * + * - Naming conventions: Arrays starting with 'b' are boundary equation arrays + * indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays + * (interior matrix equations) indexed by 'anum'. When 'num' is prefixed with + * a row or column number 'i', 'j', or 'k', the array holds the corresponding + * local data index for that row or column (e.g., an index into the local + * solution vector). Matrix coefficients are named as above, e.g., 'bij' is + * the coefficient for b_ij. + * + * NOTE: Implementation in src/sys/hypre_interface.cxx + */ +struct BCMatrixEquations { + HYPRE_Int nb; + HYPRE_Int* binum_array; + HYPRE_Int* bjnum_array; + HYPRE_Complex* bii_array; + HYPRE_Complex* bij_array; + HYPRE_Int na; + HYPRE_Int* aknum_array; + HYPRE_Complex* aki_array; + + BCMatrixEquations() = delete; + + BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows, + HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols, + HYPRE_Complex* values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int* bi_array); // row i for each boundary equation + + ~BCMatrixEquations() { + // Free arrays + HypreFree(binum_array); + HypreFree(bjnum_array); + HypreFree(bii_array); + HypreFree(bij_array); + HypreFree(aknum_array); + HypreFree(aki_array); + } + + /// Applies in-place modification of the rhs array. + /// + /// Returns an array of boundary values that can be used to apply + /// boundary conditions to a solution vector. + BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex* rhs); + + /// Apply boundary conditions to the solution. + /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() + void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex* solution); +}; + +/// A shared pointer to a BCMatrixEquations object +using BCMatrixPtr = std::shared_ptr; + template class HypreVector { MPI_Comm comm; @@ -174,6 +278,14 @@ public: HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); } + // Data for eliminating boundary equation + bool elimBErhs = false; + bool elimBEsol = false; + BCMatrixPtr bcmatrix; + BCValuesPtr bcvalues; /// Stores rhs values of BC rows + + void syncElimBErhs(HypreVector& rhs) { bcvalues = rhs.bcvalues; } + void assemble() { CALI_CXX_MARK_FUNCTION; writeCacheToHypre(); @@ -183,11 +295,17 @@ public: } void writeCacheToHypre() { + if (elimBErhs) { + bcvalues = bcmatrix->adjustBCRightHandSideEquations(V); + } checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); } void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); + if (elimBEsol) { + bcmatrix->adjustBCSolutionEquations(bcvalues, V); + } } T toField() { @@ -667,6 +785,20 @@ public: return Element(*this, global_row, global_column, positions, weights); } + // Data for eliminating boundary equations + bool elimBE = false; + BCMatrixPtr bcmatrix; // Shared pointer + + void setElimBE() { elimBE = true; } + + void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { + sol.elimBEsol = elimBE; + sol.bcmatrix = bcmatrix; + + rhs.elimBErhs = elimBE; + rhs.bcmatrix = bcmatrix; + } + void assemble() { CALI_CXX_MARK_FUNCTION; @@ -695,8 +827,32 @@ public: entry++; } } - checkHypreError( - HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); + + // Eliminate boundary condition equations in hypre SetValues input arguments + if (elimBE) { + HYPRE_Int* bi_array; + HYPRE_Int* row_indexes; + // There must be an easier way to get nb + int nb = 0; + BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; } + HypreMalloc(bi_array, nb * sizeof(HYPRE_Int)); + nb = 0; + BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { + bi_array[nb] = index_converter->getGlobal(i); + nb++; + } + + bcmatrix = std::make_shared( + num_rows, num_cols, rawI, &row_indexes, cols, vals, nb, bi_array); + HypreFree(bi_array); + + checkHypreError(HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, + row_indexes, cols, vals)); + HypreFree(row_indexes); + } else { + checkHypreError( + HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); + } checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix)); checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix, reinterpret_cast(¶llel_matrix))); @@ -876,6 +1032,31 @@ public: setMaxIter( options["maxits"].doc("Maximum iterations for Hypre solver").withDefault(10000)); + switch (solver_type) { + case HYPRE_SOLVER_TYPE::gmres: { + HYPRE_ParCSRGMRESSetKDim(solver, + options["kdim"] + .doc("Set the maximum size of the Krylov space") + .withDefault(30)); + + if (options["skip_real_residual_check"] + .doc("Skip the evaluation and the check of the actual residual?") + .withDefault(false)) { + HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); + } + break; + } + case HYPRE_SOLVER_TYPE::bicgstab: { + break; + } + case HYPRE_SOLVER_TYPE::pcg: { + break; + } + default: { + throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type)); + } + } + HYPRE_BoomerAMGCreate(&precon); HYPRE_BoomerAMGSetOldDefault(precon); #if BOUT_HAS_CUDA diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index c50be1db85..beb83a216d 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -218,11 +218,15 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { CALI_MARK_BEGIN("LaplaceHypre3d_solve:vectorAssemble"); + operator3D.setElimBEVectors(solution, rhs); + rhs.importValuesFromField(b); solution.importValuesFromField(x0); rhs.assemble(); solution.assemble(); + solution.syncElimBErhs(rhs); + CALI_MARK_END("LaplaceHypre3d_solve:vectorAssemble"); CALI_MARK_BEGIN("LaplaceHypre3d_solve:solve"); @@ -411,6 +415,7 @@ void LaplaceHypre3d::updateMatrix3D() { operator3D.ydown(ydown)(l, l.ym().zp()) += -C_d2f_dydz; operator3D.ydown(ydown)(l, l.ym().zm()) += C_d2f_dydz; } + operator3D.setElimBE(); operator3D.assemble(); if (print_matrix) { diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx new file mode 100644 index 0000000000..1838bd5158 --- /dev/null +++ b/src/sys/hypre_interface.cxx @@ -0,0 +1,131 @@ + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_HYPRE + +#include "bout/hypre_interface.hxx" + +namespace bout { + +BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, + HYPRE_BigInt* rows, HYPRE_Int** row_indexes_ptr, + HYPRE_BigInt* cols, HYPRE_Complex* values, + HYPRE_Int nb, HYPRE_Int* bi_array) + : nb(nb) { + HYPRE_Int* row_indexes; + + // Create the row_indexes array + row_indexes = (HYPRE_Int*)malloc(sizeof(HYPRE_Int) * nrows); + row_indexes[0] = 0; + for (HYPRE_Int i = 1; i < nrows; i++) { + row_indexes[i] = row_indexes[i - 1] + ncols[i - 1]; + } + + // Assume just one interior equation coupled to each boundary equation + na = nb; + + // Allocate arrays + HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); + HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); + HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); + + HYPRE_Int binum = 0; + HYPRE_Int aknum = 0; + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + // Get boundary equation information and adjust boundary equations + // Find row i in rows array (assume i increases and rows is sorted) + HYPRE_Int i = bi_array[bnum]; + for (; binum < nrows; binum++) { + if (i == rows[binum]) { + break; // Found row i in rows array + } + } + HYPRE_Int bcoeffnum = row_indexes[binum]; + HYPRE_Complex bii{0.0}, bij{0.0}; + HYPRE_Int j = 0; + + for (HYPRE_Int m = 0; m < 2; m++) { // Assume only two boundary equation coefficients + if (cols[bcoeffnum + m] == i) { + bii = values[bcoeffnum + m]; + values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix) + } else { + j = cols[bcoeffnum + m]; + bij = values[bcoeffnum + m]; + values[bcoeffnum + m] = 0.0; // Identity equation + } + } + ncols[binum] = 1; // Identity equation + + /* Get interior equation information and adjust interior equations */ + /* Find row k in rows array (assume k increases and rows is sorted) */ + HYPRE_Int k = j; // Assume equation k = j + for (; aknum < nrows; aknum++) { + if (k == rows[aknum]) { + break; // Found row k in rows array + } + } + HYPRE_Int acoeffnum = row_indexes[aknum]; + + HYPRE_Int mkj = 0; + HYPRE_Complex aki{0.0}; + for (HYPRE_Int m = 0; m < ncols[aknum]; m++) { + if (cols[acoeffnum + m] == j) { + mkj = m; // Save for update of akj value below + } + if (cols[acoeffnum + m] == i) { + aki = values[acoeffnum + m]; + values[acoeffnum + m] = 0.0; // Eliminate coupling to boundary equation + } + } + values[acoeffnum + mkj] -= aki * bij / bii; // Update akj value + + // Update arrays + HYPRE_Int anum = bnum; // Assume only one interior equation k + binum_array[bnum] = binum; + bjnum_array[bnum] = aknum; // Assume only one interior equation k + bii_array[bnum] = bii; + bij_array[bnum] = bij; + aknum_array[anum] = aknum; + aki_array[anum] = aki; + } + + // Set return arguments + *row_indexes_ptr = row_indexes; +} + +BCValuesPtr BCMatrixEquations::adjustBCRightHandSideEquations(HYPRE_Complex* rhs) { + + // Allocate array to store boundary row values + BCValuesPtr brhs = std::make_shared(nb); + + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + HYPRE_Int binum = binum_array[bnum]; + brhs->data[bnum] = rhs[binum]; + } + + for (HYPRE_Int anum = 0; anum < na; anum++) { + HYPRE_Int bnum = anum; // Assume only one interior equation per boundary equation + HYPRE_Int aknum = aknum_array[anum]; + rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum]; + } + + return brhs; +} + +void BCMatrixEquations::adjustBCSolutionEquations(BCValuesPtr brhs, + HYPRE_Complex* solution) { + + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + HYPRE_Int binum = binum_array[bnum]; + HYPRE_Int bjnum = bjnum_array[bnum]; + solution[binum] = + (brhs->data[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; + } +} + +} // namespace bout + +#endif // BOUT_HAS_HYPRE From 4ed7b90933732f679a473fc1d31b7588fd00b752 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 16 Jul 2026 13:49:59 -0700 Subject: [PATCH 2/7] hypre_interface: Use std::isfinite --- include/bout/hypre_interface.hxx | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index f9f97b5bb3..f016ee96d8 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -5,8 +5,11 @@ #if BOUT_HAS_HYPRE +#include "bout/assert.hxx" #include "bout/bout_enum_class.hxx" +#include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" +#include "bout/boutexception.hxx" #include "bout/caliper_wrapper.hxx" #include "bout/field.hxx" #include "bout/globalindexer.hxx" @@ -383,19 +386,19 @@ public: } Element& operator=(const Element& other) { - ASSERT3(finite(static_cast(other))); + ASSERT3(std::isfinite(static_cast(other))); return *this = static_cast(other); } Element& operator=(BoutReal value_) { - ASSERT3(finite(value_)); + ASSERT3(std::isfinite(value_)); value = value_; vector->V[vec_i] = value_; return *this; } Element& operator+=(BoutReal value_) { - ASSERT3(finite(value_)); + ASSERT3(std::isfinite(value_)); value += value_; - ASSERT3(finite(value)); + ASSERT3(std::isfinite(value)); vector->V[vec_i] += value_; return *this; } @@ -550,7 +553,7 @@ public: weights(weights_) { #if CHECK > 2 for (const auto val : weights) { - ASSERT3(finite(val)); + ASSERT3(std::isfinite(val)); } #endif ASSERT2(positions.size() == weights.size()); @@ -563,19 +566,19 @@ public: } Element& operator=(const Element& other) { - ASSERT3(finite(static_cast(other))); + ASSERT3(std::isfinite(static_cast(other))); return *this = static_cast(other); } Element& operator=(BoutReal value_) { - ASSERT3(finite(value_)); + ASSERT3(std::isfinite(value_)); value = value_; setValues(value); return *this; } Element& operator+=(BoutReal value_) { - ASSERT3(finite(value_)); + ASSERT3(std::isfinite(value_)); auto column_position = std::find(cbegin(positions), cend(positions), column); if (column_position != cend(positions)) { const auto i = std::distance(cbegin(positions), column_position); From f6cbdb02a2b03053419e3861cc4f3704ddb964e5 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 16 Jul 2026 16:09:15 -0700 Subject: [PATCH 3/7] HypreVector::importValuesFromField don't call assemble The vector should be assembled once after importing values. Assembling twice with elimBE over-corrects the RHS. --- include/bout/hypre_interface.hxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index f016ee96d8..c595d4cca8 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -331,6 +331,7 @@ public: return result; } + // Loads values into I and V arrays but does not assemble void importValuesFromField(const T& f) { CALI_CXX_MARK_FUNCTION; @@ -345,8 +346,6 @@ public: } ASSERT2(vec_i == vsize); - // writeCacheToHypre(); // redundant assemble already performs writeCacheToHypre - assemble(); have_indices = true; } From 9aae2b7d969f0e9e0f076a16217de36ef4951de7 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 16 Jul 2026 16:19:41 -0700 Subject: [PATCH 4/7] Hypre BCMatrixEquations: Use HypreMalloc row_indexes is free'd using HypreFree, so should be allocated using HypreMalloc. --- src/sys/hypre_interface.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx index 1838bd5158..d90392cf1e 100644 --- a/src/sys/hypre_interface.cxx +++ b/src/sys/hypre_interface.cxx @@ -15,7 +15,7 @@ BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_Int* row_indexes; // Create the row_indexes array - row_indexes = (HYPRE_Int*)malloc(sizeof(HYPRE_Int) * nrows); + HypreMalloc(row_indexes, sizeof(HYPRE_Int) * nrows); row_indexes[0] = 0; for (HYPRE_Int i = 1; i < nrows; i++) { row_indexes[i] = row_indexes[i - 1] + ncols[i - 1]; From 24d1742f1466922a082afa8a5596a6ff127638af Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Aug 2026 06:14:44 -0700 Subject: [PATCH 5/7] Clang tidy fixes --- include/bout/hypre_interface.hxx | 26 ++++++++++++++++--- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 18 +++++++++---- .../invert/laplace/test_laplace_hypre3d.cxx | 18 ++++++++----- 3 files changed, 46 insertions(+), 16 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index c595d4cca8..a118ff4787 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -12,8 +12,13 @@ #include "bout/boutexception.hxx" #include "bout/caliper_wrapper.hxx" #include "bout/field.hxx" +#include "bout/field2d.hxx" #include "bout/globalindexer.hxx" #include "bout/hyprelib.hxx" +#include "bout/options.hxx" +#include "bout/output.hxx" +#include "bout/paralleltransform.hxx" +#include "bout/region.hxx" #include "bout/utils.hxx" #include "HYPRE.h" @@ -23,7 +28,10 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" +#include +#include #include +#include // BOUT_ENUM_CLASS does not work inside namespaces BOUT_ENUM_CLASS(HYPRE_SOLVER_TYPE, gmres, bicgstab, pcg); @@ -213,7 +221,7 @@ public: other.V = nullptr; } - HypreVector& operator=(HypreVector&& other) { + HypreVector& operator=(HypreVector&& other) noexcept { comm = other.comm; jlower = other.jlower; jupper = other.jupper; @@ -291,6 +299,10 @@ public: void assemble() { CALI_CXX_MARK_FUNCTION; + + // Should not already be assembled + ASSERT1(parallel_vector == nullptr); + writeCacheToHypre(); checkHypreError(HYPRE_IJVectorAssemble(hypre_vector)); checkHypreError(HYPRE_IJVectorGetObject(hypre_vector, @@ -352,8 +364,14 @@ public: HYPRE_IJVector get() { return hypre_vector; } const HYPRE_IJVector& get() const { return hypre_vector; } - HYPRE_ParVector getParallel() { return parallel_vector; } - const HYPRE_ParVector& getParallel() const { return parallel_vector; } + HYPRE_ParVector getParallel() { + ASSERT1(parallel_vector != nullptr); + return parallel_vector; + } + const HYPRE_ParVector& getParallel() const { + ASSERT1(parallel_vector != nullptr); + return parallel_vector; + } class Element { HypreVector* vector; @@ -777,7 +795,7 @@ public: pw.begin(), pw.end(), std::back_inserter(positions), [this, ny, nz](ParallelTransform::PositionsAndWeights p) -> HYPRE_Int { return this->index_converter->getGlobal( - ind_type(p.i * ny * nz + p.j * nz + p.k, ny, nz)); + ind_type((((p.i * ny) + p.j) * nz) + p.k, ny, nz)); }); std::transform(pw.begin(), pw.end(), std::back_inserter(weights), [](ParallelTransform::PositionsAndWeights p) -> HYPRE_Complex { diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index beb83a216d..05e2ddd83f 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -3,7 +3,7 @@ * Using Hypre Solvers * ************************************************************************** - * Copyright 2021 - 2025 BOUT++ contributors + * Copyright 2021 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -30,17 +30,25 @@ #include "hypre3d_laplace.hxx" #include +#include #include +#include #include #include +#include +#include +#include #include #include #include +#include #include #include #include #include +#include +#include LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Solver*) : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), @@ -232,7 +240,7 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { CALI_MARK_BEGIN("LaplaceHypre3d_solve:solve"); // Invoke solver { - Timer timer("hypresolve"); + const Timer timer("hypresolve"); linearSystem.solve(); } @@ -309,7 +317,7 @@ void LaplaceHypre3d::updateMatrix3D() { } BoutReal C_d2f_dx2 = coords->g11[l]; - BoutReal C_d2f_dy2 = (coords->g22[l] - 1.0 / coords->g_22[l]); + BoutReal C_d2f_dy2 = (coords->g22[l] - (1.0 / coords->g_22[l])); BoutReal C_d2f_dz2 = coords->g33[l]; if (issetD) { C_d2f_dx2 *= D[l]; @@ -364,7 +372,7 @@ void LaplaceHypre3d::updateMatrix3D() { // Must add these (rather than assign) so that elements used in // interpolation don't overwrite each other. BOUT_FOR_SERIAL(l, indexer->getRegionNobndry()) { - BoutReal C_df_dy = (coords->G2[l] - dJ_dy[l] / coords->J[l]); + BoutReal C_df_dy = (coords->G2[l] - (dJ_dy[l] / coords->J[l])); if (issetD) { C_df_dy *= D[l]; } @@ -375,7 +383,7 @@ void LaplaceHypre3d::updateMatrix3D() { / C1[l]; } - BoutReal C_d2f_dy2 = (coords->g22[l] - 1.0 / coords->g_22[l]); + BoutReal C_d2f_dy2 = (coords->g22[l] - (1.0 / coords->g_22[l])); if (issetD) { C_d2f_dy2 *= D[l]; } diff --git a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx index a721b96833..92d7b3194b 100644 --- a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx +++ b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx @@ -10,14 +10,18 @@ #include "bout/invert_laplace.hxx" #include "gtest/gtest.h" +#include "bout/bout_types.hxx" +#include "bout/coordinates.hxx" #include "bout/derivs.hxx" #include "bout/difops.hxx" #include "bout/field2d.hxx" #include "bout/field3d.hxx" +#include "bout/globals.hxx" #include "bout/griddata.hxx" #include "bout/hypre_interface.hxx" #include "bout/mesh.hxx" #include "bout/options.hxx" +#include "bout/region.hxx" #include "bout/vecops.hxx" #include "fake_mesh_fixture.hxx" @@ -38,7 +42,7 @@ class ForwardOperator { upper_y_neumann = yup_neumann; } - const Field3D operator()(Field3D& f) { + Field3D operator()(Field3D& f) { Field3D result = d * Laplace_perp(f, CELL_DEFAULT, "free", "RGN_NOY") + (Grad(f) * Grad(c2) - DDY(c2) * DDY(f) / coords->g_22) / c1 + a * f + ex * DDX(f) + ez * DDZ(f); @@ -53,7 +57,7 @@ class ForwardOperator { bool inner_x_neumann, outer_x_neumann, // If false then use Dirichlet conditions lower_y_neumann, upper_y_neumann; - void applyBoundaries(Field3D& newF, Field3D& f) { + void applyBoundaries(Field3D& newF, Field3D& f) const { BOUT_FOR(i, f.getMesh()->getRegion3D("RGN_INNER_X")) { if (inner_x_neumann) { newF[i] = (f[i.xp()] - f[i]) / coords->dx[i] / sqrt(coords->g_11[i]); @@ -104,14 +108,14 @@ class LaplaceHypre3dTest coef3.allocate(); BOUT_FOR(i, mesh->getRegion2D("RGN_ALL")) { - BoutReal x = i.x() / (BoutReal)nx - 0.5; - BoutReal y = i.y() / (BoutReal)ny - 0.5; + const BoutReal x = i.x() / (BoutReal)nx - 0.5; + const BoutReal y = i.y() / (BoutReal)ny - 0.5; coef2[i] = x + y; } BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) { - BoutReal x = i.x() / (BoutReal)nx - 0.5; - BoutReal y = i.y() / (BoutReal)ny - 0.5; - BoutReal z = i.z() / (BoutReal)nz - 0.5; + const BoutReal x = i.x() / (BoutReal)nx - 0.5; + const BoutReal y = i.y() / (BoutReal)ny - 0.5; + const BoutReal z = i.z() / (BoutReal)nz - 0.5; f3[i] = 1e3 * exp(-0.5 * sqrt(x * x + y * y + z * z) / sigmasq); coef3[i] = x + y + sin(2 * 3.14159265358979323846 * z); } From 524ceab35d9a90739f3a3942622f665ecaa954b9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 14 Aug 2026 10:07:52 -0700 Subject: [PATCH 6/7] Refactored Hypre boundary row elimination The row elimination is now an implementation detail that doesn't require changes to the HypreMatrix and HypreVector interface. The elimination is applied to temporary working vectors when assembling the HypreVector, so the input values are not mutated. New tests check the boundary elimination. --- include/bout/hypre_interface.hxx | 230 +++++++++++++----- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 6 - .../laplace/impls/hypre3d/hypre3d_laplace.hxx | 2 +- src/sys/hypre_interface.cxx | 128 +++++++--- .../include/bout/test_hypre_interface.cxx | 229 +++++++++++++++-- .../invert/laplace/test_laplace_hypre3d.cxx | 11 + 6 files changed, 489 insertions(+), 117 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index a118ff4787..d71659ca3a 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -28,6 +28,7 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" +#include #include #include #include @@ -131,47 +132,70 @@ using BCValuesPtr = std::shared_ptr; * * NOTE: Implementation in src/sys/hypre_interface.cxx */ -struct BCMatrixEquations { +struct BoundaryElimination { HYPRE_Int nb; HYPRE_Int* binum_array; HYPRE_Int* bjnum_array; + HYPRE_Int* bdep_array; HYPRE_Complex* bii_array; HYPRE_Complex* bij_array; HYPRE_Int na; HYPRE_Int* aknum_array; HYPRE_Complex* aki_array; + std::vector reduction_order; + std::vector expansion_order; - BCMatrixEquations() = delete; + BoundaryElimination() = delete; - BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows, - HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols, - HYPRE_Complex* values, - HYPRE_Int nb, // number of boundary equations - HYPRE_Int* bi_array); // row i for each boundary equation + BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows, + HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols, + HYPRE_Complex* values, + HYPRE_Int nb, // number of boundary equations + HYPRE_Int* bi_array); // row i for each boundary equation - ~BCMatrixEquations() { + ~BoundaryElimination() { // Free arrays HypreFree(binum_array); HypreFree(bjnum_array); + HypreFree(bdep_array); HypreFree(bii_array); HypreFree(bij_array); HypreFree(aknum_array); HypreFree(aki_array); } + HYPRE_Int size() const { return nb; } + + /// Copy boundary-row entries from an existing full-space vector. + BCValuesPtr copyBoundaryRowValues(const HYPRE_Complex* values) const; + + /// Evaluate the original boundary equations using the supplied full-space vector. + BCValuesPtr evaluateBoundaryEquations(const HYPRE_Complex* values) const; + /// Applies in-place modification of the rhs array. /// - /// Returns an array of boundary values that can be used to apply - /// boundary conditions to a solution vector. - BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex* rhs); + /// Returns boundary values needed to reconstruct the full-space solution. + BCValuesPtr reduceRightHandSideInPlace(HYPRE_Complex* rhs) const; /// Apply boundary conditions to the solution. - /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() - void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex* solution); + /// Uses the BCValuesPtr returned from reduceRightHandSideInPlace() + void expandSolutionInPlace(BCValuesPtr brhs, HYPRE_Complex* solution) const; + + /// Reconstruct the action of the full operator from the reduced matrix-vector product. + void expandMatvecResultInPlace(BCValuesPtr boundary_operator_values, + BCValuesPtr full_boundary_values, + HYPRE_Complex* result) const; +}; + +/// A shared pointer to a BoundaryElimination object +using BoundaryEliminationPtr = std::shared_ptr; + +struct BoundaryEliminationState { + BCValuesPtr interior_values; + BCValuesPtr boundary_values; }; -/// A shared pointer to a BCMatrixEquations object -using BCMatrixPtr = std::shared_ptr; +enum class HypreVectorReadMode { standard, solution, matvec }; template class HypreVector { @@ -185,6 +209,8 @@ class HypreVector { bool have_indices{false}; HYPRE_BigInt* I{nullptr}; HYPRE_Complex* V{nullptr}; + HYPRE_Complex* workV{nullptr}; + bool cache_current{false}; HypreLib hyprelib{}; public: @@ -201,6 +227,7 @@ public: checkHypreError(HYPRE_IJVectorDestroy(hypre_vector)); HypreFree(I); HypreFree(V); + HypreFree(workV); } // Disable copy, at least for now: not clear that HYPRE_IJVector is @@ -212,13 +239,15 @@ public: : comm(other.comm), jlower(other.jlower), jupper(other.jupper), vsize(other.vsize), indexConverter(other.indexConverter), location(other.location), initialised(other.initialised), have_indices(other.have_indices), I(other.I), - V(other.V) { + V(other.V), workV(other.workV), cache_current(other.cache_current) { std::swap(hypre_vector, other.hypre_vector); std::swap(parallel_vector, other.parallel_vector); other.initialised = false; other.have_indices = false; other.I = nullptr; other.V = nullptr; + other.workV = nullptr; + other.cache_current = false; } HypreVector& operator=(HypreVector&& other) noexcept { @@ -236,8 +265,12 @@ public: other.have_indices = false; I = other.I; V = other.V; + workV = other.workV; + cache_current = other.cache_current; other.I = nullptr; other.V = nullptr; + other.workV = nullptr; + other.cache_current = false; return *this; } @@ -267,7 +300,9 @@ public: initialised = true; HypreMalloc(I, vsize * sizeof(HYPRE_BigInt)); HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); + HypreMalloc(workV, vsize * sizeof(HYPRE_Complex)); importValuesFromField(f); + assemble(); } /// Construct a vector with given index set, but don't set any values @@ -287,40 +322,70 @@ public: location = CELL_LOC::centre; HypreMalloc(I, vsize * sizeof(HYPRE_BigInt)); HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); + HypreMalloc(workV, vsize * sizeof(HYPRE_Complex)); } - // Data for eliminating boundary equation - bool elimBErhs = false; - bool elimBEsol = false; - BCMatrixPtr bcmatrix; - BCValuesPtr bcvalues; /// Stores rhs values of BC rows - - void syncElimBErhs(HypreVector& rhs) { bcvalues = rhs.bcvalues; } - - void assemble() { + void assemble(const BoundaryElimination* boundary_elimination = nullptr, + BoundaryEliminationState* elimination_state = nullptr) { CALI_CXX_MARK_FUNCTION; - // Should not already be assembled - ASSERT1(parallel_vector == nullptr); - - writeCacheToHypre(); + parallel_vector = nullptr; + checkHypreError(HYPRE_IJVectorInitialize(hypre_vector)); + writeCacheToHypre(boundary_elimination, elimination_state); checkHypreError(HYPRE_IJVectorAssemble(hypre_vector)); checkHypreError(HYPRE_IJVectorGetObject(hypre_vector, reinterpret_cast(¶llel_vector))); + cache_current = true; } - void writeCacheToHypre() { - if (elimBErhs) { - bcvalues = bcmatrix->adjustBCRightHandSideEquations(V); + void writeCacheToHypre(const BoundaryElimination* boundary_elimination = nullptr, + BoundaryEliminationState* elimination_state = nullptr) { + HYPRE_Complex* values = V; + if (boundary_elimination != nullptr) { + ASSERT1(elimination_state != nullptr); + std::copy(V, V + vsize, workV); + values = workV; + elimination_state->interior_values = + boundary_elimination->reduceRightHandSideInPlace(values); + elimination_state->boundary_values = elimination_state->interior_values; + } else if (elimination_state != nullptr) { + elimination_state->interior_values = nullptr; + elimination_state->boundary_values = nullptr; } - checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); + checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, values)); } - void readCacheFromHypre() { - checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); - if (elimBEsol) { - bcmatrix->adjustBCSolutionEquations(bcvalues, V); + void readCacheFromHypre(const BoundaryElimination* boundary_elimination = nullptr, + const BoundaryEliminationState* elimination_state = nullptr, + HypreVectorReadMode mode = HypreVectorReadMode::standard) { + HYPRE_Complex* values = V; + if ((boundary_elimination != nullptr) and (mode != HypreVectorReadMode::standard)) { + values = workV; + } + checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, values)); + if (boundary_elimination != nullptr) { + switch (mode) { + case HypreVectorReadMode::standard: + break; + case HypreVectorReadMode::solution: + ASSERT1(elimination_state != nullptr); + ASSERT1(elimination_state->interior_values != nullptr); + boundary_elimination->expandSolutionInPlace(elimination_state->interior_values, + values); + std::copy(values, values + vsize, V); + break; + case HypreVectorReadMode::matvec: + ASSERT1(elimination_state != nullptr); + ASSERT1(elimination_state->interior_values != nullptr); + ASSERT1(elimination_state->boundary_values != nullptr); + boundary_elimination->expandMatvecResultInPlace( + elimination_state->interior_values, elimination_state->boundary_values, + values); + std::copy(values, values + vsize, V); + break; + } } + cache_current = true; } T toField() { @@ -329,7 +394,9 @@ public: T result(indexConverter->getMesh()); result.allocate().setLocation(location); - readCacheFromHypre(); + if ((!cache_current) and (parallel_vector != nullptr)) { + readCacheFromHypre(); + } // Note that this only populates boundaries to a depth of 1 int count = 0; BOUT_FOR_SERIAL(i, indexConverter->getRegionAll()) { @@ -359,10 +426,12 @@ public: ASSERT2(vec_i == vsize); have_indices = true; + cache_current = true; } HYPRE_IJVector get() { return hypre_vector; } const HYPRE_IJVector& get() const { return hypre_vector; } + const HYPRE_Complex* getValues() const { return V; } HYPRE_ParVector getParallel() { ASSERT1(parallel_vector != nullptr); @@ -410,6 +479,7 @@ public: ASSERT3(std::isfinite(value_)); value = value_; vector->V[vec_i] = value_; + vector->cache_current = true; return *this; } Element& operator+=(BoutReal value_) { @@ -417,6 +487,7 @@ public: value += value_; ASSERT3(std::isfinite(value)); vector->V[vec_i] += value_; + vector->cache_current = true; return *this; } operator BoutReal() const { return value; } @@ -439,11 +510,20 @@ public: friend void swap(HypreVector& lhs, HypreVector& rhs) { using std::swap; + swap(lhs.comm, rhs.comm); + swap(lhs.jlower, rhs.jlower); + swap(lhs.jupper, rhs.jupper); + swap(lhs.vsize, rhs.vsize); swap(lhs.hypre_vector, rhs.hypre_vector); swap(lhs.parallel_vector, rhs.parallel_vector); swap(lhs.indexConverter, rhs.indexConverter); swap(lhs.location, rhs.location); swap(lhs.initialised, rhs.initialised); + swap(lhs.have_indices, rhs.have_indices); + swap(lhs.I, rhs.I); + swap(lhs.V, rhs.V); + swap(lhs.workV, rhs.workV); + swap(lhs.cache_current, rhs.cache_current); } }; @@ -463,6 +543,8 @@ class HypreMatrix { std::vector* I; std::vector>* J; std::vector>* V; + bool elimBE{false}; + BoundaryEliminationPtr boundary_elimination; HypreLib hyprelib{}; // todo also take care of I,J,V @@ -486,7 +568,8 @@ public: index_converter(other.index_converter), location(other.location), initialised(other.initialised), yoffset(other.yoffset), parallel_transform(other.parallel_transform), assembled(other.assembled), - num_rows(other.num_rows), I(other.I), J(other.J), V(other.V) { + num_rows(other.num_rows), I(other.I), J(other.J), V(other.V), + elimBE(other.elimBE), boundary_elimination(other.boundary_elimination) { std::swap(hypre_matrix, other.hypre_matrix); std::swap(parallel_matrix, other.parallel_matrix); } @@ -507,6 +590,8 @@ public: I = other.I; J = other.J; V = other.V; + elimBE = other.elimBE; + boundary_elimination = other.boundary_elimination; return *this; } @@ -805,33 +890,21 @@ public: return Element(*this, global_row, global_column, positions, weights); } - // Data for eliminating boundary equations - bool elimBE = false; - BCMatrixPtr bcmatrix; // Shared pointer - void setElimBE() { elimBE = true; } - void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { - sol.elimBEsol = elimBE; - sol.bcmatrix = bcmatrix; - - rhs.elimBErhs = elimBE; - rhs.bcmatrix = bcmatrix; - } - void assemble() { CALI_CXX_MARK_FUNCTION; HYPRE_BigInt num_entries = 0; - HYPRE_BigInt* num_cols; + HYPRE_Int* num_cols; HYPRE_BigInt* cols; HYPRE_BigInt* rawI; HYPRE_Complex* vals; - HypreMalloc(num_cols, num_rows * sizeof(HYPRE_BigInt)); + HypreMalloc(num_cols, num_rows * sizeof(HYPRE_Int)); for (HYPRE_BigInt i = 0; i < num_rows; ++i) { - num_cols[i] = (*J)[i].size(); - num_entries += (*J)[i].size(); + num_cols[i] = static_cast((*J)[i].size()); + num_entries += num_cols[i]; } HypreMalloc(rawI, num_rows * sizeof(HYPRE_BigInt)); @@ -862,7 +935,7 @@ public: nb++; } - bcmatrix = std::make_shared( + boundary_elimination = std::make_shared( num_rows, num_cols, rawI, &row_indexes, cols, vals, nb, bi_array); HypreFree(bi_array); @@ -870,6 +943,7 @@ public: row_indexes, cols, vals)); HypreFree(row_indexes); } else { + boundary_elimination = nullptr; checkHypreError( HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); } @@ -917,12 +991,17 @@ public: result.I = I; // We want the pointer to transfer so this works like a view result.J = J; result.V = V; + result.elimBE = elimBE; + result.boundary_elimination = boundary_elimination; return result; } HYPRE_IJMatrix get() { return *hypre_matrix; } const HYPRE_IJMatrix& get() const { return *hypre_matrix; } + const BoundaryElimination* getBoundaryElimination() const { + return boundary_elimination.get(); + } HYPRE_ParCSRMatrix getParallel() { return parallel_matrix; } const HYPRE_ParCSRMatrix& getParallel() const { return parallel_matrix; } @@ -932,8 +1011,33 @@ public: void computeAxpby(double alpha, HypreVector& x, double beta, HypreVector& y) { CALI_CXX_MARK_FUNCTION; + BoundaryEliminationState elimination_state; + if (boundary_elimination != nullptr) { + elimination_state.interior_values = + boundary_elimination->evaluateBoundaryEquations(x.getValues()); + elimination_state.boundary_values = + std::make_shared(boundary_elimination->size()); + BCValuesPtr y_boundary_values; + if (beta != 0.0) { + y_boundary_values = boundary_elimination->copyBoundaryRowValues(y.getValues()); + } + for (HYPRE_Int i = 0; i < boundary_elimination->size(); ++i) { + const HYPRE_Complex operator_value = elimination_state.interior_values->data[i]; + elimination_state.interior_values->data[i] = alpha * operator_value; + elimination_state.boundary_values->data[i] = + alpha * operator_value + + (y_boundary_values != nullptr ? beta * y_boundary_values->data[i] : 0.0); + } + } + + x.assemble(); + y.assemble(); checkHypreError(HYPRE_ParCSRMatrixMatvec(alpha, parallel_matrix, x.getParallel(), beta, y.getParallel())); + y.readCacheFromHypre(boundary_elimination.get(), + boundary_elimination != nullptr ? &elimination_state : nullptr, + boundary_elimination != nullptr ? HypreVectorReadMode::matvec + : HypreVectorReadMode::standard); } // y = A*x @@ -941,8 +1045,7 @@ public: void computeAx(HypreVector& x, HypreVector& y) { CALI_CXX_MARK_FUNCTION; - checkHypreError(HYPRE_ParCSRMatrixMatvec(1.0, parallel_matrix, x.getParallel(), 0.0, - y.getParallel())); + computeAxpby(1.0, x, 0.0, y); } }; @@ -1161,6 +1264,11 @@ public: ASSERT2(A != nullptr); ASSERT2(x != nullptr); ASSERT2(b != nullptr); + BoundaryEliminationState elimination_state; + const auto* boundary_elimination = A->getBoundaryElimination(); + b->assemble(boundary_elimination, + boundary_elimination != nullptr ? &elimination_state : nullptr); + x->assemble(); if (not solver_setup) { checkHypreError( solverSetup(solver, A->getParallel(), b->getParallel(), x->getParallel())); @@ -1170,6 +1278,12 @@ public: solve_err = checkHypreError( solverSolve(solver, A->getParallel(), b->getParallel(), x->getParallel())); + x->readCacheFromHypre(boundary_elimination, + boundary_elimination != nullptr ? &elimination_state : nullptr, + boundary_elimination != nullptr + ? HypreVectorReadMode::solution + : HypreVectorReadMode::standard); + return solve_err; } }; // class HypreSystem diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index 05e2ddd83f..f7ae40ed39 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -226,14 +226,8 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { CALI_MARK_BEGIN("LaplaceHypre3d_solve:vectorAssemble"); - operator3D.setElimBEVectors(solution, rhs); - rhs.importValuesFromField(b); solution.importValuesFromField(x0); - rhs.assemble(); - solution.assemble(); - - solution.syncElimBErhs(rhs); CALI_MARK_END("LaplaceHypre3d_solve:vectorAssemble"); diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx index d58ef6f688..a5ecaab45e 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx @@ -145,7 +145,7 @@ public: } // Return a reference to the matrix objects representing the Laplace - // operator. These will be (re)construct if necessary. + // operator. These will be (re)constructed if necessary. bout::HypreMatrix& getMatrix3D(); IndexerPtr getIndexer() { return indexer; } diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx index d90392cf1e..a70848f56a 100644 --- a/src/sys/hypre_interface.cxx +++ b/src/sys/hypre_interface.cxx @@ -5,14 +5,25 @@ #include "bout/hypre_interface.hxx" +#include + namespace bout { -BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, - HYPRE_BigInt* rows, HYPRE_Int** row_indexes_ptr, - HYPRE_BigInt* cols, HYPRE_Complex* values, - HYPRE_Int nb, HYPRE_Int* bi_array) +BoundaryElimination::BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols, + HYPRE_BigInt* rows, HYPRE_Int** row_indexes_ptr, + HYPRE_BigInt* cols, HYPRE_Complex* values, + HYPRE_Int nb, HYPRE_Int* bi_array) : nb(nb) { HYPRE_Int* row_indexes; + const auto find_local_row = [nrows, rows](HYPRE_BigInt row) -> HYPRE_Int { + auto row_position = std::lower_bound(rows, rows + nrows, row); + if ((row_position == rows + nrows) || (*row_position != row)) { + throw BoutException("Could not find local row {} while constructing boundary " + "elimination data", + row); + } + return static_cast(std::distance(rows, row_position)); + }; // Create the row_indexes array HypreMalloc(row_indexes, sizeof(HYPRE_Int) * nrows); @@ -27,22 +38,21 @@ BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, // Allocate arrays HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); + HypreMalloc(bdep_array, sizeof(HYPRE_Int) * nb); HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); - HYPRE_Int binum = 0; - HYPRE_Int aknum = 0; + std::vector boundary_at_row(nrows, -1); + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + boundary_at_row[find_local_row(bi_array[bnum])] = bnum; + } + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { // Get boundary equation information and adjust boundary equations - // Find row i in rows array (assume i increases and rows is sorted) HYPRE_Int i = bi_array[bnum]; - for (; binum < nrows; binum++) { - if (i == rows[binum]) { - break; // Found row i in rows array - } - } + const HYPRE_Int binum = find_local_row(i); HYPRE_Int bcoeffnum = row_indexes[binum]; HYPRE_Complex bii{0.0}, bij{0.0}; HYPRE_Int j = 0; @@ -60,13 +70,8 @@ BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, ncols[binum] = 1; // Identity equation /* Get interior equation information and adjust interior equations */ - /* Find row k in rows array (assume k increases and rows is sorted) */ HYPRE_Int k = j; // Assume equation k = j - for (; aknum < nrows; aknum++) { - if (k == rows[aknum]) { - break; // Found row k in rows array - } - } + const HYPRE_Int aknum = find_local_row(k); HYPRE_Int acoeffnum = row_indexes[aknum]; HYPRE_Int mkj = 0; @@ -86,39 +91,91 @@ BCMatrixEquations::BCMatrixEquations(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_Int anum = bnum; // Assume only one interior equation k binum_array[bnum] = binum; bjnum_array[bnum] = aknum; // Assume only one interior equation k + bdep_array[bnum] = boundary_at_row[aknum]; bii_array[bnum] = bii; bij_array[bnum] = bij; aknum_array[anum] = aknum; aki_array[anum] = aki; } + std::vector dependency_depth(nb, -1); + const auto get_depth = [&dependency_depth, this](const auto& self, + HYPRE_Int bnum) -> HYPRE_Int { + if (dependency_depth[bnum] >= 0) { + return dependency_depth[bnum]; + } + const HYPRE_Int dep = bdep_array[bnum]; + if (dep < 0) { + dependency_depth[bnum] = 0; + } else { + dependency_depth[bnum] = 1 + self(self, dep); + } + return dependency_depth[bnum]; + }; + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + get_depth(get_depth, bnum); + } + + expansion_order.resize(nb); + std::iota(expansion_order.begin(), expansion_order.end(), 0); + std::sort(expansion_order.begin(), expansion_order.end(), + [&dependency_depth](HYPRE_Int lhs, HYPRE_Int rhs) { + if (dependency_depth[lhs] != dependency_depth[rhs]) { + return dependency_depth[lhs] < dependency_depth[rhs]; + } + return lhs < rhs; + }); + reduction_order = expansion_order; + std::reverse(reduction_order.begin(), reduction_order.end()); + // Set return arguments *row_indexes_ptr = row_indexes; } -BCValuesPtr BCMatrixEquations::adjustBCRightHandSideEquations(HYPRE_Complex* rhs) { +BCValuesPtr +BoundaryElimination::copyBoundaryRowValues(const HYPRE_Complex* values) const { + BCValuesPtr boundary_values = std::make_shared(nb); - // Allocate array to store boundary row values - BCValuesPtr brhs = std::make_shared(nb); + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + boundary_values->data[bnum] = values[binum_array[bnum]]; + } + + return boundary_values; +} + +BCValuesPtr +BoundaryElimination::evaluateBoundaryEquations(const HYPRE_Complex* values) const { + BCValuesPtr boundary_values = std::make_shared(nb); for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { - HYPRE_Int binum = binum_array[bnum]; - brhs->data[bnum] = rhs[binum]; + boundary_values->data[bnum] = (bii_array[bnum] * values[binum_array[bnum]]) + + (bij_array[bnum] * values[bjnum_array[bnum]]); } - for (HYPRE_Int anum = 0; anum < na; anum++) { - HYPRE_Int bnum = anum; // Assume only one interior equation per boundary equation + return boundary_values; +} + +BCValuesPtr BoundaryElimination::reduceRightHandSideInPlace(HYPRE_Complex* rhs) const { + + // Allocate array to store boundary row values + BCValuesPtr brhs = copyBoundaryRowValues(rhs); + + for (HYPRE_Int bnum : reduction_order) { + HYPRE_Int anum = bnum; // Assume only one interior equation per boundary equation HYPRE_Int aknum = aknum_array[anum]; rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum]; + if (bdep_array[bnum] >= 0) { + brhs->data[bdep_array[bnum]] = rhs[aknum]; + } } return brhs; } -void BCMatrixEquations::adjustBCSolutionEquations(BCValuesPtr brhs, - HYPRE_Complex* solution) { +void BoundaryElimination::expandSolutionInPlace(BCValuesPtr brhs, + HYPRE_Complex* solution) const { - for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + for (HYPRE_Int bnum : expansion_order) { HYPRE_Int binum = binum_array[bnum]; HYPRE_Int bjnum = bjnum_array[bnum]; solution[binum] = @@ -126,6 +183,21 @@ void BCMatrixEquations::adjustBCSolutionEquations(BCValuesPtr brhs, } } +void BoundaryElimination::expandMatvecResultInPlace(BCValuesPtr boundary_operator_values, + BCValuesPtr full_boundary_values, + HYPRE_Complex* result) const { + for (HYPRE_Int anum = 0; anum < na; anum++) { + HYPRE_Int bnum = anum; // Assume only one interior equation per boundary equation + HYPRE_Int aknum = aknum_array[anum]; + result[aknum] += + aki_array[anum] * boundary_operator_values->data[bnum] / bii_array[bnum]; + } + + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + result[binum_array[bnum]] = full_boundary_values->data[bnum]; + } +} + } // namespace bout #endif // BOUT_HAS_HYPRE diff --git a/tests/unit/include/bout/test_hypre_interface.cxx b/tests/unit/include/bout/test_hypre_interface.cxx index ec22140647..1d83ae2734 100644 --- a/tests/unit/include/bout/test_hypre_interface.cxx +++ b/tests/unit/include/bout/test_hypre_interface.cxx @@ -1,15 +1,19 @@ +#include "bout/build_defines.hxx" + #if BOUT_HAS_HYPRE #include "HYPRE.h" #include "HYPRE_IJ_mv.h" #include "HYPRE_parcsr_ls.h" +#include #include +#include +#include #include "fake_mesh_fixture.hxx" #include "test_extras.hxx" -#include "gmock/gmock.h" #include "gtest/gtest.h" #include "bout/field3d.hxx" @@ -154,13 +158,30 @@ TYPED_TEST(HypreVectorTest, Swap) { ////////////////////////////////////////////////// // HypreMatrix tests -class MockTransform : public ParallelTransformIdentity { +class TestParallelTransform : public ParallelTransformIdentity { public: - explicit MockTransform(Mesh& mesh_in) : ParallelTransformIdentity(mesh_in){}; - MOCK_METHOD(std::vector, getWeightsForYUpApproximation, - (int i, int j, int k), (override)); - MOCK_METHOD(std::vector, getWeightsForYDownApproximation, - (int i, int j, int k), (override)); + explicit TestParallelTransform(Mesh& mesh_in) : ParallelTransformIdentity(mesh_in) {} + + std::vector getWeightsForYUpApproximation(int i, int j, + int k) override { + last_up = {i, j, k}; + up_calls += 1; + return y_up_weights; + } + + std::vector getWeightsForYDownApproximation(int i, int j, + int k) override { + last_down = {i, j, k}; + down_calls += 1; + return y_down_weights; + } + + std::vector y_up_weights; + std::vector y_down_weights; + std::array last_up{-1, -1, -1}; + std::array last_down{-1, -1, -1}; + int up_calls{0}; + int down_calls{0}; }; template @@ -169,7 +190,7 @@ class HypreMatrixTest : public FakeMeshFixture { WithQuietOutput all{output}; T field; IndexerPtr indexer; - MockTransform* pt{nullptr}; + TestParallelTransform* pt{nullptr}; std::vector yUpWeights, yDownWeights; using ind_type = typename T::ind_type; ind_type indexA, indexB, iWU0, iWU1, iWU2, iWD0, iWD1, iWD2; @@ -199,7 +220,8 @@ class HypreMatrixTest : public FakeMeshFixture { iWU1 = indexB; iWU2 = indexB.xp(); - auto transform = bout::utils::make_unique(*bout::globals::mesh); + auto transform = + bout::utils::make_unique(*bout::globals::mesh); ParallelTransform::PositionsAndWeights wUp0 = {iWU0.x(), iWU0.y(), iWU0.z(), 0.5}, wUp1 = {iWU1.x(), iWU1.y(), iWU1.z(), 1.0}, wUp2 = {iWU2.x(), iWU2.y(), iWU2.z(), 0.5}, @@ -208,8 +230,14 @@ class HypreMatrixTest : public FakeMeshFixture { wDown2 = {iWD2.x(), iWD2.y(), iWD2.z(), 0.5}; yUpWeights = {wUp0, wUp1, wUp2}; yDownWeights = {wDown0, wDown1, wDown2}; - pt = transform.get(); - field.getCoordinates()->setParallelTransform(std::move(transform)); + transform->y_up_weights = yUpWeights; + transform->y_down_weights = yDownWeights; + bout::globals::mesh->getCoordinates()->setParallelTransform(std::move(transform)); + pt = dynamic_cast( + &bout::globals::mesh->getCoordinates()->getParallelTransform()); + if (pt == nullptr) { + throw BoutException("Failed to install TestParallelTransform in HypreMatrixTest"); + } } virtual ~HypreMatrixTest() = default; }; @@ -307,7 +335,9 @@ TYPED_TEST(HypreMatrixTest, SetElements) { HYPRE_Int ncolumns{1}; HYPRE_Complex value; BOUT_OMP_SAFE(critical) - { HYPRE_IJMatrixGetValues(raw_matrix, 1, &ncolumns, &i_index, &j_index, &value); } + { + HYPRE_IJMatrixGetValues(raw_matrix, 1, &ncolumns, &i_index, &j_index, &value); + } if (i == j) { EXPECT_EQ(static_cast(value), static_cast(this->indexer->getGlobal(i))); @@ -382,8 +412,6 @@ auto IsHypreMatrixEqual(const HypreMatrix& matrix, const HypreMatrix& refe } TYPED_TEST(HypreMatrixTest, YUp) { - using namespace ::testing; - HypreMatrix matrix(this->indexer); if constexpr (std::is_same_v) { @@ -392,15 +420,11 @@ TYPED_TEST(HypreMatrixTest, YUp) { } HypreMatrix expected(this->indexer); - MockTransform* transform = this->pt; const BoutReal value = 42.0; if constexpr (std::is_same_v) { expected(this->indexA, this->indexB) = value; } else { - EXPECT_CALL(*transform, getWeightsForYUpApproximation( - this->indexB.x(), this->indexA.y(), this->indexB.z())) - .WillOnce(Return(this->yUpWeights)); expected(this->indexA, this->iWU0) = this->yUpWeights[0].weight * value; expected(this->indexA, this->iWU1) = this->yUpWeights[1].weight * value; expected(this->indexA, this->iWU2) = this->yUpWeights[2].weight * value; @@ -408,6 +432,13 @@ TYPED_TEST(HypreMatrixTest, YUp) { matrix.yup()(this->indexA, this->indexB) = value; + if constexpr (std::is_same_v) { + EXPECT_EQ(this->pt->up_calls, 1); + EXPECT_EQ(this->pt->last_up[0], this->indexB.x()); + EXPECT_EQ(this->pt->last_up[1], this->indexA.y()); + EXPECT_EQ(this->pt->last_up[2], this->indexB.z()); + } + expected.assemble(); matrix.assemble(); @@ -415,8 +446,6 @@ TYPED_TEST(HypreMatrixTest, YUp) { } TYPED_TEST(HypreMatrixTest, YDown) { - using namespace ::testing; - HypreMatrix matrix(this->indexer); if constexpr (std::is_same_v) { @@ -425,15 +454,11 @@ TYPED_TEST(HypreMatrixTest, YDown) { } HypreMatrix expected(this->indexer); - MockTransform* transform = this->pt; const BoutReal value = 42.0; if constexpr (std::is_same_v) { expected(this->indexB, this->indexA) = value; } else { - EXPECT_CALL(*transform, getWeightsForYDownApproximation( - this->indexA.x(), this->indexB.y(), this->indexA.z())) - .WillOnce(Return(this->yDownWeights)); expected(this->indexB, this->iWD0) = this->yDownWeights[0].weight * value; expected(this->indexB, this->iWD1) = this->yDownWeights[1].weight * value; expected(this->indexB, this->iWD2) = this->yDownWeights[2].weight * value; @@ -441,6 +466,13 @@ TYPED_TEST(HypreMatrixTest, YDown) { matrix.ydown()(this->indexB, this->indexA) = value; + if constexpr (std::is_same_v) { + EXPECT_EQ(this->pt->down_calls, 1); + EXPECT_EQ(this->pt->last_down[0], this->indexA.x()); + EXPECT_EQ(this->pt->last_down[1], this->indexB.y()); + EXPECT_EQ(this->pt->last_down[2], this->indexA.z()); + } + expected.assemble(); matrix.assemble(); @@ -462,4 +494,153 @@ TYPED_TEST(HypreMatrixTest, YNext0) { EXPECT_TRUE(IsHypreMatrixEqual(matrix, expected)); } +namespace { + +struct RawBoundaryEliminationSystem { + std::vector ncols; + std::vector rows; + std::vector cols; + std::vector values; + std::vector boundary_rows; + HYPRE_Int* row_indexes{nullptr}; + std::unique_ptr elimination; + + RawBoundaryEliminationSystem(std::initializer_list ncols_in, + std::initializer_list rows_in, + std::initializer_list cols_in, + std::initializer_list values_in, + std::initializer_list boundary_rows_in) + : ncols(ncols_in), rows(rows_in), cols(cols_in), values(values_in), + boundary_rows(boundary_rows_in) { + elimination = std::make_unique( + static_cast(rows.size()), ncols.data(), rows.data(), &row_indexes, + cols.data(), values.data(), static_cast(boundary_rows.size()), + boundary_rows.data()); + } + + ~RawBoundaryEliminationSystem() { + if (row_indexes != nullptr) { + HypreFree(row_indexes); + } + } +}; + +RawBoundaryEliminationSystem makeSingleBoundarySystem() { + return {{2, 3, 2}, + {0, 1, 2}, + {0, 1, 0, 1, 2, 1, 2}, + {2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0}, + {0}}; +} + +RawBoundaryEliminationSystem makeBackwardCoupledRowSystem() { + return {{2, 3, 3, 2}, + {0, 1, 2, 3}, + {0, 2, 1, 2, 3, 0, 2, 3, 1, 3}, + {2.0, 3.0, 31.0, 37.0, 29.0, 5.0, 7.0, 11.0, 23.0, 19.0}, + {0, 3}}; +} + +} // namespace + +TEST(BoundaryEliminationTest, TransformsSingleBoundarySystem) { + auto system = makeSingleBoundarySystem(); + + EXPECT_EQ(system.ncols[0], 1); + EXPECT_EQ(system.values[0], -1.0); + EXPECT_EQ(system.values[1], 0.0); + EXPECT_EQ(system.values[2], 0.0); + EXPECT_DOUBLE_EQ(system.values[3], -0.5); + EXPECT_EQ(system.values[4], 11.0); + EXPECT_EQ(system.values[5], 13.0); + EXPECT_EQ(system.values[6], 17.0); + + ASSERT_NE(system.elimination, nullptr); + EXPECT_EQ(system.elimination->binum_array[0], 0); + EXPECT_EQ(system.elimination->bjnum_array[0], 1); + EXPECT_EQ(system.elimination->bdep_array[0], -1); + EXPECT_EQ(system.elimination->bii_array[0], 2.0); + EXPECT_EQ(system.elimination->bij_array[0], 3.0); + EXPECT_EQ(system.elimination->aknum_array[0], 1); + EXPECT_EQ(system.elimination->aki_array[0], 5.0); +} + +TEST(BoundaryEliminationTest, ReducesRightHandSideForSingleBoundarySystem) { + auto system = makeSingleBoundarySystem(); + std::array rhs{{19.0, 23.0, 29.0}}; + + auto brhs = system.elimination->reduceRightHandSideInPlace(rhs.data()); + + ASSERT_NE(brhs, nullptr); + EXPECT_EQ(brhs->data[0], 19.0); + EXPECT_EQ(rhs[0], 19.0); + EXPECT_DOUBLE_EQ(rhs[1], -24.5); + EXPECT_EQ(rhs[2], 29.0); +} + +TEST(BoundaryEliminationTest, ExpandsSolutionForSingleBoundarySystem) { + auto system = makeSingleBoundarySystem(); + auto brhs = std::make_shared(1); + brhs->data[0] = 19.0; + std::array solution{{0.0, 2.0, 3.0}}; + + system.elimination->expandSolutionInPlace(brhs, solution.data()); + + EXPECT_DOUBLE_EQ(solution[0], 6.5); + EXPECT_EQ(solution[1], 2.0); + EXPECT_EQ(solution[2], 3.0); +} + +TEST(BoundaryEliminationTest, ReconstructsMatvecForSingleBoundarySystem) { + auto system = makeSingleBoundarySystem(); + std::array x{{1.0, 2.0, 3.0}}; + std::array reduced_result{{-1.0, 32.0, 77.0}}; + auto boundary_values = system.elimination->evaluateBoundaryEquations(x.data()); + + system.elimination->expandMatvecResultInPlace(boundary_values, boundary_values, + reduced_result.data()); + + EXPECT_EQ(reduced_result[0], 8.0); + EXPECT_EQ(reduced_result[1], 52.0); + EXPECT_EQ(reduced_result[2], 77.0); +} + +TEST(BoundaryEliminationTest, HandlesBackwardCoupledRows) { + auto system = makeBackwardCoupledRowSystem(); + + EXPECT_EQ(system.ncols[0], 1); + EXPECT_EQ(system.ncols[3], 1); + EXPECT_EQ(system.values[0], -1.0); + EXPECT_EQ(system.values[1], 0.0); + EXPECT_DOUBLE_EQ(system.values[2], 31.0 - (29.0 * 23.0 / 19.0)); + EXPECT_EQ(system.values[4], 0.0); + EXPECT_EQ(system.values[5], 0.0); + EXPECT_DOUBLE_EQ(system.values[6], -0.5); + EXPECT_EQ(system.values[8], 0.0); + EXPECT_EQ(system.values[9], -1.0); + + EXPECT_EQ(system.elimination->aknum_array[0], 2); + EXPECT_EQ(system.elimination->aknum_array[1], 1); +} + +TEST(BoundaryEliminationTest, DifferentEliminationObjectsProduceDifferentReducedSystems) { + auto system_a = makeSingleBoundarySystem(); + RawBoundaryEliminationSystem system_b{{2, 3, 2}, + {0, 1, 2}, + {0, 1, 0, 1, 2, 1, 2}, + {4.0, -1.0, 6.0, 9.0, 10.0, 13.0, 17.0}, + {0}}; + std::array rhs_a{{19.0, 23.0, 29.0}}; + std::array rhs_b = rhs_a; + + auto brhs_a = system_a.elimination->reduceRightHandSideInPlace(rhs_a.data()); + auto brhs_b = system_b.elimination->reduceRightHandSideInPlace(rhs_b.data()); + + EXPECT_EQ(brhs_a->data[0], 19.0); + EXPECT_EQ(brhs_b->data[0], 19.0); + EXPECT_DOUBLE_EQ(rhs_a[1], -24.5); + EXPECT_DOUBLE_EQ(rhs_b[1], -5.5); + EXPECT_NE(rhs_a[1], rhs_b[1]); +} + #endif // BOUT_HAS_HYPRE diff --git a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx index 92d7b3194b..3924245bad 100644 --- a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx +++ b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx @@ -351,6 +351,17 @@ TEST_P(LaplaceHypre3dTest, TestSolve3DGuess) { BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) { EXPECT_NEAR(expected[i], actual[i], tol); } } +TEST_P(LaplaceHypre3dTest, TestSolve3DRepeated) { + const Field3D rhs = forward(f3); + const Field3D first = solver.solve(rhs); + const Field3D second = solver.solve(rhs); + + BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) { + EXPECT_NEAR(f3[i], first[i], tol); + EXPECT_NEAR(f3[i], second[i], tol); + } +} + TEST_P(LaplaceHypre3dTest, TestSolvePerp) { FieldPerp f(1.0); EXPECT_THROW(solver.solve(f), BoutException); From 7913512c4959115fd4c509cab4101e104e51a2a6 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 17 Aug 2026 15:36:59 -0700 Subject: [PATCH 7/7] Generalize boundary elimination code Fixes the integrated Laplace-Hypre3D tests. GPT says: The boundary-elimination logic in [src/sys/hypre_interface.cxx (line 48)](/Users/dudson2/code/BOUT-next/src/sys/hypre_interface.cxx:48) was only removing each boundary variable from one coupled row; in the hypre3d stencil, the same Y-boundary value can appear in multiple nearby rows, especially in the SOL case. I changed BoundaryElimination to collect all representable row couplings for each boundary equation, store them with per-boundary offsets, and apply them during matrix reduction, RHS reduction, and matvec reconstruction. The interface metadata was updated in [include/bout/hypre_interface.hxx (line 103)](/Users/dudson2/code/BOUT-next/include/bout/hypre_interface.hxx:103), and I added multi-coupling unit coverage in [tests/unit/include/bout/test_hypre_interface.cxx (line 618)](/Users/dudson2/code/BOUT-next/tests/unit/include/bout/test_hypre_interface.cxx:618). --- include/bout/hypre_interface.hxx | 12 +- src/sys/hypre_interface.cxx | 136 +++++++++++++----- .../test_laplace_hypre3d.py | 2 +- .../include/bout/test_hypre_interface.cxx | 66 +++++++++ 4 files changed, 173 insertions(+), 43 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index d71659ca3a..6a53c27f0c 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -103,12 +103,16 @@ using BCValuesPtr = std::shared_ptr; * * b_ii * u_i + b_ij * u_j = rhs_i * - * - We also assume that each boundary equation has only one interior equation k - * coupled to it (such that k = j) with coupling coefficient a_ki + * - We also assume that each boundary equation has only one retained unknown j + * on the boundary row + * + * b_ii * u_i + b_ij * u_j = rhs_i + * + * - Any number of other equations k may couple to u_i with coefficient a_ki * * a_ki * u_i + a_kj * u_j + ... = rhs_k * - * - Each equation k is adjusted as follows: + * - Each such equation k is adjusted as follows: * * a_kj = a_kj - a_ki * b_ij / b_ii * a_ki = 0 @@ -140,6 +144,7 @@ struct BoundaryElimination { HYPRE_Complex* bii_array; HYPRE_Complex* bij_array; HYPRE_Int na; + HYPRE_Int* aoffset_array; HYPRE_Int* aknum_array; HYPRE_Complex* aki_array; std::vector reduction_order; @@ -160,6 +165,7 @@ struct BoundaryElimination { HypreFree(bdep_array); HypreFree(bii_array); HypreFree(bij_array); + HypreFree(aoffset_array); HypreFree(aknum_array); HypreFree(aki_array); } diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx index a70848f56a..101484913a 100644 --- a/src/sys/hypre_interface.cxx +++ b/src/sys/hypre_interface.cxx @@ -6,6 +6,7 @@ #include "bout/hypre_interface.hxx" #include +#include namespace bout { @@ -32,32 +33,49 @@ BoundaryElimination::BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols, row_indexes[i] = row_indexes[i - 1] + ncols[i - 1]; } - // Assume just one interior equation coupled to each boundary equation - na = nb; - // Allocate arrays HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb); HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb); HypreMalloc(bdep_array, sizeof(HYPRE_Int) * nb); HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb); HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb); - HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); - HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); std::vector boundary_at_row(nrows, -1); for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { boundary_at_row[find_local_row(bi_array[bnum])] = bnum; } + struct BoundaryOccurrence { + HYPRE_Int rownum; + HYPRE_Int row_offset; + }; + std::unordered_map boundary_number_for_row; + boundary_number_for_row.reserve(nb); + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + boundary_number_for_row.emplace(bi_array[bnum], bnum); + } + + std::vector> occurrences_by_boundary(nb); + for (HYPRE_Int rownum = 0; rownum < nrows; rownum++) { + const HYPRE_Int row_start = row_indexes[rownum]; + for (HYPRE_Int m = 0; m < ncols[rownum]; m++) { + const auto boundary_position = boundary_number_for_row.find(cols[row_start + m]); + if (boundary_position != boundary_number_for_row.end()) { + occurrences_by_boundary[boundary_position->second].push_back({rownum, m}); + } + } + } + + std::vector bjrow_array(nb, -1); for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { // Get boundary equation information and adjust boundary equations HYPRE_Int i = bi_array[bnum]; const HYPRE_Int binum = find_local_row(i); HYPRE_Int bcoeffnum = row_indexes[binum]; HYPRE_Complex bii{0.0}, bij{0.0}; - HYPRE_Int j = 0; + HYPRE_BigInt j = -1; - for (HYPRE_Int m = 0; m < 2; m++) { // Assume only two boundary equation coefficients + for (HYPRE_Int m = 0; m < ncols[binum]; m++) { if (cols[bcoeffnum + m] == i) { bii = values[bcoeffnum + m]; values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix) @@ -67,37 +85,75 @@ BoundaryElimination::BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols, values[bcoeffnum + m] = 0.0; // Identity equation } } - ncols[binum] = 1; // Identity equation - - /* Get interior equation information and adjust interior equations */ - HYPRE_Int k = j; // Assume equation k = j - const HYPRE_Int aknum = find_local_row(k); - HYPRE_Int acoeffnum = row_indexes[aknum]; - - HYPRE_Int mkj = 0; - HYPRE_Complex aki{0.0}; - for (HYPRE_Int m = 0; m < ncols[aknum]; m++) { - if (cols[acoeffnum + m] == j) { - mkj = m; // Save for update of akj value below - } - if (cols[acoeffnum + m] == i) { - aki = values[acoeffnum + m]; - values[acoeffnum + m] = 0.0; // Eliminate coupling to boundary equation - } + if (j < 0) { + throw BoutException("Boundary row {} does not contain a retained neighbour", i); + } + if (bii == 0.0) { + throw BoutException("Boundary row {} has zero diagonal coefficient", i); } - values[acoeffnum + mkj] -= aki * bij / bii; // Update akj value + ncols[binum] = 1; // Identity equation // Update arrays - HYPRE_Int anum = bnum; // Assume only one interior equation k binum_array[bnum] = binum; - bjnum_array[bnum] = aknum; // Assume only one interior equation k - bdep_array[bnum] = boundary_at_row[aknum]; + bjnum_array[bnum] = find_local_row(j); + bdep_array[bnum] = boundary_at_row[bjnum_array[bnum]]; bii_array[bnum] = bii; bij_array[bnum] = bij; - aknum_array[anum] = aknum; - aki_array[anum] = aki; + bjrow_array[bnum] = j; } + std::vector aoffsets(nb + 1, 0); + std::vector aknums; + std::vector akis; + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + const HYPRE_Int binum = binum_array[bnum]; + const HYPRE_Int bjnum = bjnum_array[bnum]; + const HYPRE_BigInt j = bjrow_array[bnum]; + + for (const auto& occurrence : occurrences_by_boundary[bnum]) { + const HYPRE_Int aknum = occurrence.rownum; + if (aknum == binum) { + continue; + } + if ((boundary_at_row[aknum] >= 0) && (aknum != bjnum)) { + continue; + } + + const HYPRE_Int acoeffnum = row_indexes[aknum]; + const HYPRE_Int aki_position = acoeffnum + occurrence.row_offset; + const HYPRE_Complex aki = values[aki_position]; + if (aki == 0.0) { + continue; + } + + HYPRE_Int j_position = -1; + for (HYPRE_Int m = 0; m < ncols[aknum]; m++) { + if (cols[acoeffnum + m] == j) { + j_position = acoeffnum + m; + break; + } + } + if (j_position < 0) { + continue; + } + + values[aki_position] = 0.0; // Eliminate coupling to boundary equation + values[j_position] -= aki * bij_array[bnum] / bii_array[bnum]; + + aknums.push_back(aknum); + akis.push_back(aki); + } + aoffsets[bnum + 1] = static_cast(aknums.size()); + } + + na = static_cast(aknums.size()); + HypreMalloc(aoffset_array, sizeof(HYPRE_Int) * (nb + 1)); + std::copy(aoffsets.begin(), aoffsets.end(), aoffset_array); + HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na); + std::copy(aknums.begin(), aknums.end(), aknum_array); + HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na); + std::copy(akis.begin(), akis.end(), aki_array); + std::vector dependency_depth(nb, -1); const auto get_depth = [&dependency_depth, this](const auto& self, HYPRE_Int bnum) -> HYPRE_Int { @@ -161,11 +217,12 @@ BCValuesPtr BoundaryElimination::reduceRightHandSideInPlace(HYPRE_Complex* rhs) BCValuesPtr brhs = copyBoundaryRowValues(rhs); for (HYPRE_Int bnum : reduction_order) { - HYPRE_Int anum = bnum; // Assume only one interior equation per boundary equation - HYPRE_Int aknum = aknum_array[anum]; - rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum]; + for (HYPRE_Int anum = aoffset_array[bnum]; anum < aoffset_array[bnum + 1]; anum++) { + HYPRE_Int aknum = aknum_array[anum]; + rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum]; + } if (bdep_array[bnum] >= 0) { - brhs->data[bdep_array[bnum]] = rhs[aknum]; + brhs->data[bdep_array[bnum]] = rhs[bjnum_array[bnum]]; } } @@ -186,11 +243,12 @@ void BoundaryElimination::expandSolutionInPlace(BCValuesPtr brhs, void BoundaryElimination::expandMatvecResultInPlace(BCValuesPtr boundary_operator_values, BCValuesPtr full_boundary_values, HYPRE_Complex* result) const { - for (HYPRE_Int anum = 0; anum < na; anum++) { - HYPRE_Int bnum = anum; // Assume only one interior equation per boundary equation - HYPRE_Int aknum = aknum_array[anum]; - result[aknum] += - aki_array[anum] * boundary_operator_values->data[bnum] / bii_array[bnum]; + for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { + for (HYPRE_Int anum = aoffset_array[bnum]; anum < aoffset_array[bnum + 1]; anum++) { + HYPRE_Int aknum = aknum_array[anum]; + result[aknum] += + aki_array[anum] * boundary_operator_values->data[bnum] / bii_array[bnum]; + } } for (HYPRE_Int bnum = 0; bnum < nb; bnum++) { diff --git a/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py b/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py index e2fa11c13a..014b91ef00 100755 --- a/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py +++ b/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py @@ -17,7 +17,7 @@ def test_laplace_hypre3d(): success = True for directory, nproc in test_directories: - command = "test-laplace3d -d " + directory + command = "./test-laplace3d -d " + directory print("running on", nproc, "processors:", command) launch_safe(command, nproc=nproc) diff --git a/tests/unit/include/bout/test_hypre_interface.cxx b/tests/unit/include/bout/test_hypre_interface.cxx index 1d83ae2734..2be393f1b6 100644 --- a/tests/unit/include/bout/test_hypre_interface.cxx +++ b/tests/unit/include/bout/test_hypre_interface.cxx @@ -533,6 +533,14 @@ RawBoundaryEliminationSystem makeSingleBoundarySystem() { {0}}; } +RawBoundaryEliminationSystem makeMultiCoupledBoundarySystem() { + return {{2, 3, 3, 2}, + {0, 1, 2, 3}, + {0, 1, 0, 1, 3, 0, 1, 2, 2, 3}, + {2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0, 19.0, 23.0, 29.0}, + {0}}; +} + RawBoundaryEliminationSystem makeBackwardCoupledRowSystem() { return {{2, 3, 3, 2}, {0, 1, 2, 3}, @@ -561,6 +569,8 @@ TEST(BoundaryEliminationTest, TransformsSingleBoundarySystem) { EXPECT_EQ(system.elimination->bdep_array[0], -1); EXPECT_EQ(system.elimination->bii_array[0], 2.0); EXPECT_EQ(system.elimination->bij_array[0], 3.0); + EXPECT_EQ(system.elimination->aoffset_array[0], 0); + EXPECT_EQ(system.elimination->aoffset_array[1], 1); EXPECT_EQ(system.elimination->aknum_array[0], 1); EXPECT_EQ(system.elimination->aki_array[0], 5.0); } @@ -605,6 +615,59 @@ TEST(BoundaryEliminationTest, ReconstructsMatvecForSingleBoundarySystem) { EXPECT_EQ(reduced_result[2], 77.0); } +TEST(BoundaryEliminationTest, EliminatesBoundaryCouplingsFromAllAffectedRows) { + auto system = makeMultiCoupledBoundarySystem(); + + EXPECT_EQ(system.ncols[0], 1); + EXPECT_EQ(system.values[0], -1.0); + EXPECT_EQ(system.values[1], 0.0); + EXPECT_EQ(system.values[2], 0.0); + EXPECT_DOUBLE_EQ(system.values[3], -0.5); + EXPECT_EQ(system.values[4], 11.0); + EXPECT_EQ(system.values[5], 0.0); + EXPECT_DOUBLE_EQ(system.values[6], -2.5); + EXPECT_EQ(system.values[7], 19.0); + EXPECT_EQ(system.values[8], 23.0); + EXPECT_EQ(system.values[9], 29.0); + + EXPECT_EQ(system.elimination->na, 2); + EXPECT_EQ(system.elimination->aoffset_array[0], 0); + EXPECT_EQ(system.elimination->aoffset_array[1], 2); + EXPECT_EQ(system.elimination->aknum_array[0], 1); + EXPECT_EQ(system.elimination->aknum_array[1], 2); + EXPECT_EQ(system.elimination->aki_array[0], 5.0); + EXPECT_EQ(system.elimination->aki_array[1], 13.0); +} + +TEST(BoundaryEliminationTest, ReducesRightHandSideForMultipleAffectedRows) { + auto system = makeMultiCoupledBoundarySystem(); + std::array rhs{{19.0, 23.0, 31.0, 37.0}}; + + auto brhs = system.elimination->reduceRightHandSideInPlace(rhs.data()); + + ASSERT_NE(brhs, nullptr); + EXPECT_EQ(brhs->data[0], 19.0); + EXPECT_EQ(rhs[0], 19.0); + EXPECT_DOUBLE_EQ(rhs[1], -24.5); + EXPECT_DOUBLE_EQ(rhs[2], -92.5); + EXPECT_EQ(rhs[3], 37.0); +} + +TEST(BoundaryEliminationTest, ReconstructsMatvecForMultipleAffectedRows) { + auto system = makeMultiCoupledBoundarySystem(); + std::array x{{1.0, 2.0, 3.0, 4.0}}; + std::array reduced_result{{-1.0, 32.0, 52.0, 185.0}}; + auto boundary_values = system.elimination->evaluateBoundaryEquations(x.data()); + + system.elimination->expandMatvecResultInPlace(boundary_values, boundary_values, + reduced_result.data()); + + EXPECT_EQ(reduced_result[0], 8.0); + EXPECT_EQ(reduced_result[1], 52.0); + EXPECT_EQ(reduced_result[2], 104.0); + EXPECT_EQ(reduced_result[3], 185.0); +} + TEST(BoundaryEliminationTest, HandlesBackwardCoupledRows) { auto system = makeBackwardCoupledRowSystem(); @@ -619,6 +682,9 @@ TEST(BoundaryEliminationTest, HandlesBackwardCoupledRows) { EXPECT_EQ(system.values[8], 0.0); EXPECT_EQ(system.values[9], -1.0); + EXPECT_EQ(system.elimination->aoffset_array[0], 0); + EXPECT_EQ(system.elimination->aoffset_array[1], 1); + EXPECT_EQ(system.elimination->aoffset_array[2], 2); EXPECT_EQ(system.elimination->aknum_array[0], 2); EXPECT_EQ(system.elimination->aknum_array[1], 1); }