From 5ce6cae6ea94a76417a4c058fdd848db384d36d1 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Tue, 4 Mar 2025 17:04:44 -0800 Subject: [PATCH 01/12] First version of code to eliminate boundary equations to improve HYPRE solves --- include/bout/hypre_interface.hxx | 193 +++++++++++++- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 6 + .../laplace/impls/hypre3d/hypre_boundary.c | 252 ++++++++++++++++++ 3 files changed, 450 insertions(+), 1 deletion(-) create mode 100644 src/invert/laplace/impls/hypre3d/hypre_boundary.c diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index ae392de4f3..72eec916e2 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,6 +20,56 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" +HYPRE_Int +AdjustBCMatrixEquations( + 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 + HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) + HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) + HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) + HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) + HYPRE_Int *na_ptr, // number of interior equations to adjust + HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) + HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) + +HYPRE_Int +AdjustBCRightHandSideEquations( + HYPRE_Complex *rhs, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex **brhs_array_ptr, + HYPRE_Int na, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array); + +HYPRE_Int +AdjustBCSolutionEquations( + HYPRE_Complex *solution, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array); + +HYPRE_Int +AdjustBCEquationsFree( + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array); + #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -174,6 +224,24 @@ public: HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); } + // Data for eliminating boundary equations (TODO: convert this to a structure) + bool elimBErhs = false; + bool elimBEsol = false; + HYPRE_Complex *rhs; + HYPRE_Int nb; + HYPRE_Int *binum_array; + HYPRE_Int *bjnum_array; + HYPRE_Complex *bii_array; + HYPRE_Complex *bij_array; + HYPRE_Complex *brhs_array; // This is not created from the HypreMatrix class + HYPRE_Int na; + HYPRE_Int *aknum_array; + HYPRE_Complex *aki_array; + + void syncElimBErhs(HypreVector& rhs) { + brhs_array = rhs.brhs_array; + } + void assemble() { CALI_CXX_MARK_FUNCTION; writeCacheToHypre(); @@ -183,11 +251,20 @@ public: } void writeCacheToHypre() { + if (elimBErhs) + { + AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, + na, aknum_array, aki_array); + } checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); } void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); + if (elimBEsol) + { + AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); + } } T toField() { @@ -667,6 +744,43 @@ public: return Element(*this, global_row, global_column, positions, weights); } + // Data for eliminating boundary equations (TODO: convert this to a structure) + bool elimBE = false; + 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; + + void setElimBE() { + elimBE = true; + } + + void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { + sol.elimBEsol = elimBE; + sol.nb = nb; + sol.binum_array = binum_array; + sol.bjnum_array = bjnum_array; + sol.bii_array = bii_array; + sol.bij_array = bij_array; + sol.na = na; + sol.aknum_array = aknum_array; + sol.aki_array = aki_array; + + rhs.elimBErhs = elimBE; + rhs.nb = nb; + rhs.binum_array = binum_array; + rhs.bjnum_array = bjnum_array; + rhs.bii_array = bii_array; + rhs.bij_array = bij_array; + rhs.na = na; + rhs.aknum_array = aknum_array; + rhs.aki_array = aki_array; + } + void assemble() { CALI_CXX_MARK_FUNCTION; @@ -695,8 +809,85 @@ public: entry++; } } - checkHypreError( + + // 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 + 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++; + } + AdjustBCMatrixEquations(num_rows, num_cols, rawI, &row_indexes, cols, vals, + nb, bi_array, + &binum_array, &bjnum_array, &bii_array, &bij_array, + &na, &aknum_array, &aki_array); + HypreFree(bi_array); +// { +// FILE *file; +// int i; +// file = fopen("zbout.setvalues.out2", "w"); +// fprintf(file, "nrows %d\n", num_rows); +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "ncols %d: %d\n", i, num_cols[i]); +// } +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "rows %d: %d\n", i, rawI[i]); +// } +// for (i = 0; i < num_entries; i++) +// { +// fprintf(file, "cols %d: %d\n", i, cols[i]); +// } +// for (i = 0; i < num_entries; i++) +// { +// fprintf(file, "vals %d: %f\n", i, vals[i]); +// } +// fclose(file); +// +// file = fopen("zbout.setvalues.out3", "w"); +// for (i = 0; i < num_rows; i++) +// { +// fprintf(file, "row_indexes %d: %d\n", i, row_indexes[i]); +// } +// fprintf(file, "nb %d\n", nb); +// for (i = 0; i < nb; i++) +// { +// fprintf(file, "bijnum_array %d: %d %d\n", i, binum_array[i], bjnum_array[i]); +// } +// for (i = 0; i < nb; i++) +// { +// fprintf(file, "biiij_array %d: %f %f\n", i, bii_array[i], bij_array[i]); +// } +// fprintf(file, "na %d\n", na); +// for (i = 0; i < na; i++) +// { +// fprintf(file, "aknum_array %d: %d\n", i, aknum_array[i]); +// } +// for (i = 0; i < na; i++) +// { +// fprintf(file, "bki_array %d: %f\n", i, aki_array[i]); +// } +// fclose(file); +// } + 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))); diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index c50be1db85..def54b36b5 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -28,6 +28,7 @@ #if BOUT_HAS_HYPRE #include "hypre3d_laplace.hxx" +#include "hypre_boundary.c" #include #include @@ -218,11 +219,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 +416,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/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c new file mode 100644 index 0000000000..c4e7393c3f --- /dev/null +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -0,0 +1,252 @@ + +/* + * 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. + * + * Returned info arrays can be freed with the function AdjustEquationsFree(). + * + * 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. + */ + +HYPRE_Int +AdjustBCMatrixEquations( + 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 + HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) + HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) + HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) + HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) + HYPRE_Int *na_ptr, // number of interior equations to adjust + HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) + HYPRE_Complex **aki_array_ptr) // coefficient a_ki (for each interior equation) +{ + HYPRE_Int *row_indexes; + HYPRE_Int na, *binum_array, *bjnum_array, *aknum_array; + HYPRE_Complex *bii_array, *bij_array, *aki_array; + HYPRE_Int i, j, k, m, mkj, anum, bnum, acoeffnum, bcoeffnum; + HYPRE_Int binum, aknum; + HYPRE_Complex bii, bij, aki; + + /* Create the row_indexes array */ + row_indexes = (HYPRE_Int *)malloc(sizeof(HYPRE_Int) * nrows); + row_indexes[0] = 0; + for (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 return 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); + + binum = 0; + aknum = 0; + for (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) */ + i = bi_array[bnum]; + for (; binum < nrows; binum++) + { + if (i == rows[binum]) + { + break; // Found row i in rows array + } + } + bcoeffnum = row_indexes[binum]; + for (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) */ + k = j; // Assume equation k = j + for (; aknum < nrows; aknum++) + { + if (k == rows[aknum]) + { + break; // Found row k in rows array + } + } + acoeffnum = row_indexes[aknum]; + for (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 return arrays */ + 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; + *binum_array_ptr = binum_array; + *bjnum_array_ptr = bjnum_array; + *bii_array_ptr = bii_array; + *bij_array_ptr = bij_array; + *na_ptr = na; + *aknum_array_ptr = aknum_array; + *aki_array_ptr = aki_array; + + return 0; +} + +HYPRE_Int +AdjustBCRightHandSideEquations( + HYPRE_Complex *rhs, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex **brhs_array_ptr, + HYPRE_Int na, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array) +{ + HYPRE_Complex *brhs_array; + HYPRE_Int anum, bnum, binum, aknum; + + HypreMalloc(brhs_array, sizeof(HYPRE_Complex) * nb); + + for (bnum = 0; bnum < nb; bnum++) + { + binum = binum_array[bnum]; + brhs_array[bnum] = rhs[binum]; + } + + for (anum = 0; anum < na; anum++) + { + bnum = anum; // Assume only one interior equation per boundary equation + aknum = aknum_array[anum]; + rhs[aknum] -= aki_array[anum] * bij_array[bnum] / bii_array[bnum]; + } + + *brhs_array_ptr = brhs_array; + + return 0; +} + +HYPRE_Int +AdjustBCSolutionEquations( + HYPRE_Complex *solution, + HYPRE_Int nb, + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array) +{ + HYPRE_Int bnum, binum, bjnum; + + for (bnum = 0; bnum < nb; bnum++) + { + binum = binum_array[bnum]; + bjnum = bjnum_array[bnum]; + solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; + } + + return 0; +} + +HYPRE_Int +AdjustBCEquationsFree( + HYPRE_Int *binum_array, + HYPRE_Int *bjnum_array, + HYPRE_Complex *bii_array, + HYPRE_Complex *bij_array, + HYPRE_Complex *brhs_array, + HYPRE_Int *aknum_array, + HYPRE_Complex *aki_array) +{ + HypreFree(binum_array); + HypreFree(bjnum_array); + HypreFree(bii_array); + HypreFree(bij_array); + HypreFree(brhs_array); + HypreFree(aknum_array); + HypreFree(aki_array); + + return 0; +} + From e3f4a98ac6235b0bef6ffc4d90ec2da2f1ecb668 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 08:35:17 -0800 Subject: [PATCH 02/12] Added HYPRE GMRES SetKDim call. Need to make this an input file parameter. --- include/bout/hypre_interface.hxx | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 72eec916e2..efb620241e 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1067,6 +1067,22 @@ public: setMaxIter( options["maxits"].doc("Maximum iterations for Hypre solver").withDefault(10000)); + switch (solver_type) { + case HYPRE_SOLVER_TYPE::gmres: { + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + 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 From ef069675e39a6e92355036c45f3f20183a67e628 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 09:53:21 -0800 Subject: [PATCH 03/12] Added code to free up data. The HypreMatrix destroy needs to be updated. --- examples/elm-pb/elm_pb.cxx | 134 ++++++++++++++---- include/bout/hypre_interface.hxx | 41 +++--- .../laplace/impls/hypre3d/hypre_boundary.c | 35 +---- 3 files changed, 129 insertions(+), 81 deletions(-) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index f830f3d98a..c332edd223 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -46,6 +46,7 @@ class ELMpb : public PhysicsModel { Coordinates::FieldMetric U0; // 0th vorticity of equilibrium flow, // radial flux coordinate, normalized radial flux coordinate + bool laplace_perp; // Use Laplace_perp or Delp2? bool constn0; // the total height, average width and center of profile of N0 BoutReal n0_height, n0_ave, n0_width, n0_center, n0_bottom_x, Nbar, Tibar, Tebar; @@ -350,7 +351,8 @@ class ELMpb : public PhysicsModel { ////////////////////////////////////////////////////////////// auto& globalOptions = Options::root(); auto& options = globalOptions["highbeta"]; - + laplace_perp = options["laplace_perp"].withDefault(false); + // Use Laplace_perp rather than Delp2 constn0 = options["constn0"].withDefault(true); // use the hyperbolic profile of n0. If both n0_fake_prof and // T0_fake_prof are false, use the profiles from grid file @@ -1315,7 +1317,11 @@ class ELMpb : public PhysicsModel { Ctmp.applyBoundary(); Ctmp -= phi; // Now contains error in the boundary - C_phi = Delp2(phi) - U; // Error in the bulk + if (laplace_perp) { + C_phi = Laplace_perp(phi) - U; // Error in the bulk + } else { + C_phi = Delp2(phi) - U; // Error in the bulk + } C_phi.setBoundaryTo(Ctmp); } else { @@ -1347,8 +1353,12 @@ class ELMpb : public PhysicsModel { } else { ubyn = U / N0; if (diamag) { - ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); - mesh->communicate(ubyn); + if (laplace_perp) { + ubyn -= 0.5 * dnorm / (N0 * B0) * Laplace_perp(P); + } else { + ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); + } + mesh->communicate(ubyn); } // Invert laplacian for phi phiSolver->setCoefC(N0); @@ -1361,9 +1371,18 @@ class ELMpb : public PhysicsModel { if (!evolve_jpar) { // Get J from Psi - Jpar = Delp2(Psi); + if (laplace_perp) { + Jpar = Laplace_perp(Psi); + } else { + Jpar = Delp2(Psi); + } + if (include_rmp) { - Jpar += Delp2(rmp_Psi); + if (laplace_perp) { + Jpar += Laplace_perp(rmp_Psi); + } else { + Jpar += Delp2(rmp_Psi); + } } Jpar.applyBoundary(); @@ -1397,8 +1416,11 @@ class ELMpb : public PhysicsModel { } // Get Delp2(J) from J - Jpar2 = Delp2(Jpar); - + if (laplace_perp) { + Jpar2 = Laplace_perp(Jpar); + } else { + Jpar2 = Delp2(Jpar); + } Jpar2.applyBoundary(); mesh->communicate(Jpar2); @@ -1494,7 +1516,11 @@ class ELMpb : public PhysicsModel { // Jpar Field3D B0U = B0 * U; mesh->communicate(B0U); - ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + if (laplace_perp) { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Laplace_perp(Jpar); + }else { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + } if (relax_j_vac) { // Make ddt(Jpar) relax to zero. @@ -1524,11 +1550,19 @@ class ELMpb : public PhysicsModel { } if (hyperresist > 0.0) { // Hyper-resistivity - ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + if (laplace_perp) { + ddt(Psi) -= eta * hyperresist * Laplace_perp(Jpar); + } else { + ddt(Psi) -= eta * hyperresist * Delp2(Jpar); + } } if (ehyperviscos > 0.0) { // electron Hyper-viscosity coefficient - ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + if (laplace_perp) { + ddt(Psi) -= eta * ehyperviscos * Laplace_perp(Jpar2); + } else { + ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + } } // Parallel hyper-viscous diffusion for vector potential @@ -1599,7 +1633,11 @@ class ELMpb : public PhysicsModel { } if (viscos_perp > 0.0) { - ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + if (laplace_perp) { + ddt(U) += viscos_perp * Laplace_perp(U); // Perpendicular viscosity + } else { + ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + } } // Hyper-viscosity @@ -1626,21 +1664,40 @@ class ELMpb : public PhysicsModel { Pi = 0.5 * P; Pi0 = 0.5 * P0; - Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); - Dperp2Phi0.applyBoundary(); - mesh->communicate(Dperp2Phi0); + if (laplace_perp) { + Dperp2Phi0 = Field3D(Laplace_perp(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); - Dperp2Phi = Delp2(B0 * phi); - Dperp2Phi.applyBoundary(); - mesh->communicate(Dperp2Phi); + Dperp2Phi = Laplace_perp(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); - Dperp2Pi0 = Field3D(Delp2(Pi0)); - Dperp2Pi0.applyBoundary(); - mesh->communicate(Dperp2Pi0); + Dperp2Pi0 = Field3D(Laplace_perp(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); - Dperp2Pi = Delp2(Pi); - Dperp2Pi.applyBoundary(); - mesh->communicate(Dperp2Pi); + Dperp2Pi = Laplace_perp(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } else { + Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); + + Dperp2Phi = Delp2(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); + + Dperp2Pi0 = Field3D(Delp2(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); + + Dperp2Pi = Delp2(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + + } bracketPhi0P = bracket(B0 * phi0, Pi, bm_exb); bracketPhi0P.applyBoundary(); @@ -1658,8 +1715,13 @@ class ELMpb : public PhysicsModel { mesh->communicate(B0phi0); ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP0) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + } if (nonlinear) { Field3D B0phi = B0 * phi; @@ -1670,7 +1732,11 @@ class ELMpb : public PhysicsModel { ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / B0; ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + } } } @@ -1808,7 +1874,12 @@ class ELMpb : public PhysicsModel { int precon(BoutReal UNUSED(t), BoutReal gamma, BoutReal UNUSED(delta)) { // First matrix, applying L mesh->communicate(ddt(Psi)); - Field3D Jrhs = Delp2(ddt(Psi)); + Field3D Jrhs; + if (laplace_perp) { + Jrhs = Laplace_perp(ddt(Psi)); + } else { + Jrhs = Delp2(ddt(Psi)); + } Jrhs.applyBoundary("neumann"); if (jpar_bndry_width > 0) { @@ -1880,8 +1951,11 @@ class ELMpb : public PhysicsModel { phi = phiSolver->solve(ddt(U)); - Jpar = Delp2(ddt(Psi)); - + if (laplace_perp) { + Jpar = Laplace_perp(ddt(Psi)); + } else { + Jpar = Delp2(ddt(Psi)); + } mesh->communicate(phi, Jpar); Field3D JP = -b0xGrad_dot_Grad(phi, P0); diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index efb620241e..0025ed062b 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,7 +20,7 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" -HYPRE_Int +void AdjustBCMatrixEquations( HYPRE_Int nrows, HYPRE_Int *ncols, @@ -38,7 +38,7 @@ AdjustBCMatrixEquations( HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) -HYPRE_Int +void AdjustBCRightHandSideEquations( HYPRE_Complex *rhs, HYPRE_Int nb, @@ -50,7 +50,7 @@ AdjustBCRightHandSideEquations( HYPRE_Int *aknum_array, HYPRE_Complex *aki_array); -HYPRE_Int +void AdjustBCSolutionEquations( HYPRE_Complex *solution, HYPRE_Int nb, @@ -60,16 +60,6 @@ AdjustBCSolutionEquations( HYPRE_Complex *bij_array, HYPRE_Complex *brhs_array); -HYPRE_Int -AdjustBCEquationsFree( - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array); - #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -136,6 +126,9 @@ public: checkHypreError(HYPRE_IJVectorDestroy(hypre_vector)); HypreFree(I); HypreFree(V); + if (elimBErhs) { + HypreFree(brhs_array); + } } // Disable copy, at least for now: not clear that HYPRE_IJVector is @@ -251,8 +244,7 @@ public: } void writeCacheToHypre() { - if (elimBErhs) - { + if (elimBErhs) { AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, na, aknum_array, aki_array); } @@ -261,8 +253,7 @@ public: void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); - if (elimBEsol) - { + if (elimBEsol) { AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); } } @@ -419,6 +410,19 @@ public: using ind_type = typename T::ind_type; HypreMatrix() = default; +// The 'if' block below needs to be called when the HypreMatrix is destroyed, +// but I don't know where to put it. The HypreMatrix class is handled +// differently from HypreVector. +// ~HypreVector() { +// if (elimBE) { +// HypreFree(binum_array); +// HypreFree(bjnum_array); +// HypreFree(bii_array); +// HypreFree(bij_array); +// HypreFree(aknum_array); +// HypreFree(aki_array); +// } +// } HypreMatrix(const HypreMatrix&) = delete; HypreMatrix(HypreMatrix&& other) : comm(other.comm), ilower(other.ilower), iupper(other.iupper), @@ -811,8 +815,7 @@ public: } // Eliminate boundary condition equations in hypre SetValues input arguments - if (elimBE) - { + if (elimBE) { HYPRE_Int *bi_array; HYPRE_Int *row_indexes; // There must be an easier way to get nb diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c index c4e7393c3f..69e9a2d1b0 100644 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -10,8 +10,6 @@ * returns info needed to adjust the right-hand-side and solution vector through * the functions AdjustRightHandSideEquations and AdjustSolutionEquations. * - * Returned info arrays can be freed with the function AdjustEquationsFree(). - * * NOTE: It may make sense from an organizational standpoint to collect many of * these arguments in a structure of some sort. * @@ -49,7 +47,7 @@ * the coefficient for b_ij. */ -HYPRE_Int +void AdjustBCMatrixEquations( HYPRE_Int nrows, HYPRE_Int *ncols, @@ -168,11 +166,9 @@ AdjustBCMatrixEquations( *na_ptr = na; *aknum_array_ptr = aknum_array; *aki_array_ptr = aki_array; - - return 0; } -HYPRE_Int +void AdjustBCRightHandSideEquations( HYPRE_Complex *rhs, HYPRE_Int nb, @@ -203,11 +199,9 @@ AdjustBCRightHandSideEquations( } *brhs_array_ptr = brhs_array; - - return 0; } -HYPRE_Int +void AdjustBCSolutionEquations( HYPRE_Complex *solution, HYPRE_Int nb, @@ -225,28 +219,5 @@ AdjustBCSolutionEquations( bjnum = bjnum_array[bnum]; solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; } - - return 0; -} - -HYPRE_Int -AdjustBCEquationsFree( - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array) -{ - HypreFree(binum_array); - HypreFree(bjnum_array); - HypreFree(bii_array); - HypreFree(bij_array); - HypreFree(brhs_array); - HypreFree(aknum_array); - HypreFree(aki_array); - - return 0; } From 67660069a1bd2174bdc36a0750f8a9d79b0a88a9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 5 Mar 2025 12:11:45 -0800 Subject: [PATCH 04/12] hypre_interface: Wrap BC arrays in structs & shared_ptr Use structs to group together arrays for boundary row elimination. Coefficients from HypreMatrix are stored in one struct, using a shared_ptr to share values with vectors. Eliminated boundary row values need to be passed from RHS to solution vector. This now uses a shared pointer that will manage the HypreFree call. --- include/bout/hypre_interface.hxx | 202 ++++++++++++++----------------- 1 file changed, 93 insertions(+), 109 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 0025ed062b..7ad50a008d 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -98,6 +98,82 @@ 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() { + 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; + +/// Contains information needed to eliminate boundary equations from +/// RHS vectors, and restore boundary values in solution vectors. +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 + : nb(nb) { + AdjustBCMatrixEquations(nrows, ncols, rows, row_indexes_ptr, cols, values, + nb, bi_array, + // Outputs + &binum_array, &bjnum_array, &bii_array, &bij_array, + &na, &aknum_array, &aki_array); + } + + ~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) { + BCValuesPtr brhs = std::make_shared(); + AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, bij_array, &brhs->data, + na, aknum_array, aki_array); + return brhs; + } + + /// Apply boundary conditions to the solution. + /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() + void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex *solution) { + AdjustBCSolutionEquations(solution, nb, + binum_array, bjnum_array, bii_array, bij_array, brhs->data); + } +}; + +/// A shared pointer to a BCMatrixEquations object +using BCMatrixPtr = std::shared_ptr; + + template class HypreVector { MPI_Comm comm; @@ -126,9 +202,6 @@ public: checkHypreError(HYPRE_IJVectorDestroy(hypre_vector)); HypreFree(I); HypreFree(V); - if (elimBErhs) { - HypreFree(brhs_array); - } } // Disable copy, at least for now: not clear that HYPRE_IJVector is @@ -217,22 +290,14 @@ public: HypreMalloc(V, vsize * sizeof(HYPRE_Complex)); } - // Data for eliminating boundary equations (TODO: convert this to a structure) + // Data for eliminating boundary equation bool elimBErhs = false; bool elimBEsol = false; - HYPRE_Complex *rhs; - HYPRE_Int nb; - HYPRE_Int *binum_array; - HYPRE_Int *bjnum_array; - HYPRE_Complex *bii_array; - HYPRE_Complex *bij_array; - HYPRE_Complex *brhs_array; // This is not created from the HypreMatrix class - HYPRE_Int na; - HYPRE_Int *aknum_array; - HYPRE_Complex *aki_array; + BCMatrixPtr bcmatrix; + BCValuesPtr bcvalues; /// Stores rhs values of BC rows void syncElimBErhs(HypreVector& rhs) { - brhs_array = rhs.brhs_array; + bcvalues = rhs.bcvalues; } void assemble() { @@ -245,8 +310,7 @@ public: void writeCacheToHypre() { if (elimBErhs) { - AdjustBCRightHandSideEquations(V, nb, binum_array, bii_array, bij_array, &brhs_array, - na, aknum_array, aki_array); + bcvalues = bcmatrix->adjustBCRightHandSideEquations(V); } checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V)); } @@ -254,7 +318,7 @@ public: void readCacheFromHypre() { checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V)); if (elimBEsol) { - AdjustBCSolutionEquations(V, nb, binum_array, bjnum_array, bii_array, bij_array, brhs_array); + bcmatrix->adjustBCSolutionEquations(bcvalues, V); } } @@ -410,19 +474,6 @@ public: using ind_type = typename T::ind_type; HypreMatrix() = default; -// The 'if' block below needs to be called when the HypreMatrix is destroyed, -// but I don't know where to put it. The HypreMatrix class is handled -// differently from HypreVector. -// ~HypreVector() { -// if (elimBE) { -// HypreFree(binum_array); -// HypreFree(bjnum_array); -// HypreFree(bii_array); -// HypreFree(bij_array); -// HypreFree(aknum_array); -// HypreFree(aki_array); -// } -// } HypreMatrix(const HypreMatrix&) = delete; HypreMatrix(HypreMatrix&& other) : comm(other.comm), ilower(other.ilower), iupper(other.iupper), @@ -748,16 +799,9 @@ public: return Element(*this, global_row, global_column, positions, weights); } - // Data for eliminating boundary equations (TODO: convert this to a structure) + // Data for eliminating boundary equations bool elimBE = false; - 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; + BCMatrixPtr bcmatrix; // Shared pointer void setElimBE() { elimBE = true; @@ -765,24 +809,10 @@ public: void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { sol.elimBEsol = elimBE; - sol.nb = nb; - sol.binum_array = binum_array; - sol.bjnum_array = bjnum_array; - sol.bii_array = bii_array; - sol.bij_array = bij_array; - sol.na = na; - sol.aknum_array = aknum_array; - sol.aki_array = aki_array; + sol.bcmatrix = bcmatrix; rhs.elimBErhs = elimBE; - rhs.nb = nb; - rhs.binum_array = binum_array; - rhs.bjnum_array = bjnum_array; - rhs.bii_array = bii_array; - rhs.bij_array = bij_array; - rhs.na = na; - rhs.aknum_array = aknum_array; - rhs.aki_array = aki_array; + rhs.bcmatrix = bcmatrix; } void assemble() { @@ -819,7 +849,7 @@ public: HYPRE_Int *bi_array; HYPRE_Int *row_indexes; // There must be an easier way to get nb - nb = 0; + int nb = 0; BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; } @@ -829,59 +859,13 @@ public: bi_array[nb] = index_converter->getGlobal(i); nb++; } - AdjustBCMatrixEquations(num_rows, num_cols, rawI, &row_indexes, cols, vals, - nb, bi_array, - &binum_array, &bjnum_array, &bii_array, &bij_array, - &na, &aknum_array, &aki_array); + + bcmatrix = std::make_shared(num_rows, num_cols, + rawI, &row_indexes, + cols, vals, + nb, bi_array); HypreFree(bi_array); -// { -// FILE *file; -// int i; -// file = fopen("zbout.setvalues.out2", "w"); -// fprintf(file, "nrows %d\n", num_rows); -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "ncols %d: %d\n", i, num_cols[i]); -// } -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "rows %d: %d\n", i, rawI[i]); -// } -// for (i = 0; i < num_entries; i++) -// { -// fprintf(file, "cols %d: %d\n", i, cols[i]); -// } -// for (i = 0; i < num_entries; i++) -// { -// fprintf(file, "vals %d: %f\n", i, vals[i]); -// } -// fclose(file); -// -// file = fopen("zbout.setvalues.out3", "w"); -// for (i = 0; i < num_rows; i++) -// { -// fprintf(file, "row_indexes %d: %d\n", i, row_indexes[i]); -// } -// fprintf(file, "nb %d\n", nb); -// for (i = 0; i < nb; i++) -// { -// fprintf(file, "bijnum_array %d: %d %d\n", i, binum_array[i], bjnum_array[i]); -// } -// for (i = 0; i < nb; i++) -// { -// fprintf(file, "biiij_array %d: %f %f\n", i, bii_array[i], bij_array[i]); -// } -// fprintf(file, "na %d\n", na); -// for (i = 0; i < na; i++) -// { -// fprintf(file, "aknum_array %d: %d\n", i, aknum_array[i]); -// } -// for (i = 0; i < na; i++) -// { -// fprintf(file, "bki_array %d: %f\n", i, aki_array[i]); -// } -// fclose(file); -// } + checkHypreError( HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, row_indexes, cols, vals)); HypreFree(row_indexes); From 29b65c592a229b122da3660fe98a40a23d1f3ed5 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 5 Mar 2025 12:40:04 -0800 Subject: [PATCH 05/12] Fixed a bug in rhs update for boundary elimination code --- include/bout/hypre_interface.hxx | 3 +-- src/invert/laplace/impls/hypre3d/hypre_boundary.c | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 7ad50a008d..c8377f7f5a 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -44,7 +44,6 @@ AdjustBCRightHandSideEquations( HYPRE_Int nb, HYPRE_Int *binum_array, HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, HYPRE_Complex **brhs_array_ptr, HYPRE_Int na, HYPRE_Int *aknum_array, @@ -157,7 +156,7 @@ struct BCMatrixEquations { /// boundary conditions to a solution vector. BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex *rhs) { BCValuesPtr brhs = std::make_shared(); - AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, bij_array, &brhs->data, + AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, &brhs->data, na, aknum_array, aki_array); return brhs; } diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c index 69e9a2d1b0..30f0968765 100644 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ b/src/invert/laplace/impls/hypre3d/hypre_boundary.c @@ -174,7 +174,6 @@ AdjustBCRightHandSideEquations( HYPRE_Int nb, HYPRE_Int *binum_array, HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, HYPRE_Complex **brhs_array_ptr, HYPRE_Int na, HYPRE_Int *aknum_array, @@ -195,7 +194,7 @@ AdjustBCRightHandSideEquations( { bnum = anum; // Assume only one interior equation per boundary equation aknum = aknum_array[anum]; - rhs[aknum] -= aki_array[anum] * bij_array[bnum] / bii_array[bnum]; + rhs[aknum] -= aki_array[anum] * brhs_array[bnum] / bii_array[bnum]; } *brhs_array_ptr = brhs_array; From 08b1261f86e38eb45864265a4d48ec644545256b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 6 Mar 2025 09:43:05 -0800 Subject: [PATCH 06/12] hypre_interface: Moved BC elimination to src/sys/hypre_interface.cxx Changed C free functions to be members of BCMatrixEquations. Clang formatted, minor style changes. Test results unchanged. --- CMakeLists.txt | 1 + include/bout/hypre_interface.hxx | 187 +++++++-------- .../laplace/impls/hypre3d/hypre3d_laplace.cxx | 1 - .../laplace/impls/hypre3d/hypre_boundary.c | 222 ------------------ src/sys/hypre_interface.cxx | 131 +++++++++++ 5 files changed, 215 insertions(+), 327 deletions(-) delete mode 100644 src/invert/laplace/impls/hypre3d/hypre_boundary.c create mode 100644 src/sys/hypre_interface.cxx diff --git a/CMakeLists.txt b/CMakeLists.txt index c45fca3b72..af628a7db5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -335,6 +335,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 c8377f7f5a..33e6218c8f 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -20,45 +20,6 @@ #include "HYPRE_utilities.h" #include "_hypre_utilities.h" -void -AdjustBCMatrixEquations( - 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 - HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) - HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) - HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) - HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) - HYPRE_Int *na_ptr, // number of interior equations to adjust - HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) - HYPRE_Complex **aki_array_ptr); // coefficient a_ki (for each interior equation) - -void -AdjustBCRightHandSideEquations( - HYPRE_Complex *rhs, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex **brhs_array_ptr, - HYPRE_Int na, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array); - -void -AdjustBCSolutionEquations( - HYPRE_Complex *solution, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array); - #include // BOUT_ENUM_CLASS does not work inside namespaces @@ -99,46 +60,83 @@ int checkHypreError(int error) { /// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed. struct HypreComplexArray { - HYPRE_Complex *data; + HYPRE_Complex* data; - ~HypreComplexArray() { - HypreFree(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; -/// Contains information needed to eliminate boundary equations from -/// RHS vectors, and restore boundary values in solution vectors. +/*! + * 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; + 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 - : nb(nb) { - AdjustBCMatrixEquations(nrows, ncols, rows, row_indexes_ptr, cols, values, - nb, bi_array, - // Outputs - &binum_array, &bjnum_array, &bii_array, &bij_array, - &na, &aknum_array, &aki_array); - } + 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 @@ -154,25 +152,16 @@ struct BCMatrixEquations { /// /// Returns an array of boundary values that can be used to apply /// boundary conditions to a solution vector. - BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex *rhs) { - BCValuesPtr brhs = std::make_shared(); - AdjustBCRightHandSideEquations(rhs, nb, binum_array, bii_array, &brhs->data, - na, aknum_array, aki_array); - return brhs; - } + BCValuesPtr adjustBCRightHandSideEquations(HYPRE_Complex* rhs); /// Apply boundary conditions to the solution. /// Uses the BCValuesPtr returned from adjustBCRightHandSideEquations() - void adjustBCSolutionEquations(BCValuesPtr brhs, HYPRE_Complex *solution) { - AdjustBCSolutionEquations(solution, nb, - binum_array, bjnum_array, bii_array, bij_array, brhs->data); - } + 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; @@ -293,11 +282,9 @@ public: bool elimBErhs = false; bool elimBEsol = false; BCMatrixPtr bcmatrix; - BCValuesPtr bcvalues; /// Stores rhs values of BC rows + BCValuesPtr bcvalues; /// Stores rhs values of BC rows - void syncElimBErhs(HypreVector& rhs) { - bcvalues = rhs.bcvalues; - } + void syncElimBErhs(HypreVector& rhs) { bcvalues = rhs.bcvalues; } void assemble() { CALI_CXX_MARK_FUNCTION; @@ -802,9 +789,7 @@ public: bool elimBE = false; BCMatrixPtr bcmatrix; // Shared pointer - void setElimBE() { - elimBE = true; - } + void setElimBE() { elimBE = true; } void setElimBEVectors(HypreVector& sol, HypreVector& rhs) { sol.elimBEsol = elimBE; @@ -845,34 +830,28 @@ public: // Eliminate boundary condition equations in hypre SetValues input arguments if (elimBE) { - HYPRE_Int *bi_array; - HYPRE_Int *row_indexes; + 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++; - } + 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++; + bi_array[nb] = index_converter->getGlobal(i); + nb++; } - bcmatrix = std::make_shared(num_rows, num_cols, - rawI, &row_indexes, - cols, vals, - nb, bi_array); + 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)); + checkHypreError(HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI, + row_indexes, cols, vals)); HypreFree(row_indexes); - } - else - { + } else { checkHypreError( - HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); + HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals)); } checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix)); checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix, @@ -1055,7 +1034,7 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter break; } case HYPRE_SOLVER_TYPE::bicgstab: { diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index def54b36b5..beb83a216d 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -28,7 +28,6 @@ #if BOUT_HAS_HYPRE #include "hypre3d_laplace.hxx" -#include "hypre_boundary.c" #include #include diff --git a/src/invert/laplace/impls/hypre3d/hypre_boundary.c b/src/invert/laplace/impls/hypre3d/hypre_boundary.c deleted file mode 100644 index 30f0968765..0000000000 --- a/src/invert/laplace/impls/hypre3d/hypre_boundary.c +++ /dev/null @@ -1,222 +0,0 @@ - -/* - * 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. - */ - -void -AdjustBCMatrixEquations( - 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 - HYPRE_Int **binum_array_ptr, // data index for row i (for each boundary equation) - HYPRE_Int **bjnum_array_ptr, // data index for col j (for each boundary equation) - HYPRE_Complex **bii_array_ptr, // coefficient b_ii (for each boundary equation) - HYPRE_Complex **bij_array_ptr, // coefficient b_ij (for each boundary equation) - HYPRE_Int *na_ptr, // number of interior equations to adjust - HYPRE_Int **aknum_array_ptr, // data index for row k (for each interior equation) - HYPRE_Complex **aki_array_ptr) // coefficient a_ki (for each interior equation) -{ - HYPRE_Int *row_indexes; - HYPRE_Int na, *binum_array, *bjnum_array, *aknum_array; - HYPRE_Complex *bii_array, *bij_array, *aki_array; - HYPRE_Int i, j, k, m, mkj, anum, bnum, acoeffnum, bcoeffnum; - HYPRE_Int binum, aknum; - HYPRE_Complex bii, bij, aki; - - /* Create the row_indexes array */ - row_indexes = (HYPRE_Int *)malloc(sizeof(HYPRE_Int) * nrows); - row_indexes[0] = 0; - for (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 return 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); - - binum = 0; - aknum = 0; - for (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) */ - i = bi_array[bnum]; - for (; binum < nrows; binum++) - { - if (i == rows[binum]) - { - break; // Found row i in rows array - } - } - bcoeffnum = row_indexes[binum]; - for (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) */ - k = j; // Assume equation k = j - for (; aknum < nrows; aknum++) - { - if (k == rows[aknum]) - { - break; // Found row k in rows array - } - } - acoeffnum = row_indexes[aknum]; - for (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 return arrays */ - 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; - *binum_array_ptr = binum_array; - *bjnum_array_ptr = bjnum_array; - *bii_array_ptr = bii_array; - *bij_array_ptr = bij_array; - *na_ptr = na; - *aknum_array_ptr = aknum_array; - *aki_array_ptr = aki_array; -} - -void -AdjustBCRightHandSideEquations( - HYPRE_Complex *rhs, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex **brhs_array_ptr, - HYPRE_Int na, - HYPRE_Int *aknum_array, - HYPRE_Complex *aki_array) -{ - HYPRE_Complex *brhs_array; - HYPRE_Int anum, bnum, binum, aknum; - - HypreMalloc(brhs_array, sizeof(HYPRE_Complex) * nb); - - for (bnum = 0; bnum < nb; bnum++) - { - binum = binum_array[bnum]; - brhs_array[bnum] = rhs[binum]; - } - - for (anum = 0; anum < na; anum++) - { - bnum = anum; // Assume only one interior equation per boundary equation - aknum = aknum_array[anum]; - rhs[aknum] -= aki_array[anum] * brhs_array[bnum] / bii_array[bnum]; - } - - *brhs_array_ptr = brhs_array; -} - -void -AdjustBCSolutionEquations( - HYPRE_Complex *solution, - HYPRE_Int nb, - HYPRE_Int *binum_array, - HYPRE_Int *bjnum_array, - HYPRE_Complex *bii_array, - HYPRE_Complex *bij_array, - HYPRE_Complex *brhs_array) -{ - HYPRE_Int bnum, binum, bjnum; - - for (bnum = 0; bnum < nb; bnum++) - { - binum = binum_array[bnum]; - bjnum = bjnum_array[bnum]; - solution[binum] = (brhs_array[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum]; - } -} - 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 eea57fc2c1a798a80eca0d7687fa8219ffe09f9a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 6 Mar 2025 09:44:53 -0800 Subject: [PATCH 07/12] test-laplace-hypre3d: Suppress BOUT++ log outputs Difficult to see the test results amongst the log outputs. Now pipe the logs and save to files. --- tests/integrated/test-laplace-hypre3d/runtest | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/tests/integrated/test-laplace-hypre3d/runtest b/tests/integrated/test-laplace-hypre3d/runtest index b50c5993b7..f1f5950547 100755 --- a/tests/integrated/test-laplace-hypre3d/runtest +++ b/tests/integrated/test-laplace-hypre3d/runtest @@ -19,20 +19,24 @@ build_and_log("Laplace 3D with Hypre") 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) + s, out = launch_safe(command, nproc=nproc, pipe=True) + # Save output to log file + with open("run.log." + directory, "w") as f: + f.write(out) error_max = collect("error_max", path=directory, info=False) if error_max > tolerance: - print(directory + " failed with maximum error {}".format(error_max)) + print(" => " + directory + " failed with maximum error {}".format(error_max)) success = False else: - print(directory + " passed with maximum error {}".format(error_max)) + print(" => " + directory + " passed with maximum error {}".format(error_max)) if success: - print("All passed") + print("=> All passed") exit(0) else: + print("=> Some tests failed") exit(1) From ef9fbed3e68ff6ea45294e95cb174457cc5c939b Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Thu, 6 Mar 2025 17:46:47 +0000 Subject: [PATCH 08/12] Apply black changes --- tests/integrated/test_suite | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrated/test_suite b/tests/integrated/test_suite index 307a8d84b3..77ad7882c4 100755 --- a/tests/integrated/test_suite +++ b/tests/integrated/test_suite @@ -188,7 +188,7 @@ class Test(threading.Thread): self.output += "\n(It is likely that a timeout occured)" else: # ❌ Failed - print("\u274C", end="") # No newline + print("\u274c", end="") # No newline print(" %7.3f s" % (time.time() - self.local.start_time), flush=True) def _cost(self): From 9d1b46908fab0c61950905eefa73ad735b4a4c27 Mon Sep 17 00:00:00 2001 From: Rob Falgout Date: Wed, 19 Mar 2025 14:39:19 -0700 Subject: [PATCH 09/12] Turn off test for false convergence in hypre GMRES solver --- include/bout/hypre_interface.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 33e6218c8f..1837d5e275 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1034,7 +1034,8 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter + HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); // TODO: Make this an input file parameter break; } case HYPRE_SOLVER_TYPE::bicgstab: { From f969c011249b60b337817432c026b1bd312291cc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 19 Mar 2025 14:43:14 -0700 Subject: [PATCH 10/12] hypre_interface: Add kdim and skip_real_residual_check options These options control the GMRES solver. --- include/bout/hypre_interface.hxx | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 1837d5e275..15ca97cb8a 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1034,8 +1034,16 @@ public: switch (solver_type) { case HYPRE_SOLVER_TYPE::gmres: { - HYPRE_ParCSRGMRESSetKDim(solver, 30); // TODO: Make this an input file parameter - HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); // TODO: Make this an input file parameter + 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: { From 3853d3dca954e42001dd1231eaae6631c35a2426 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 27 May 2025 16:38:01 -0700 Subject: [PATCH 11/12] Use python f-strings Cleans up test code Co-authored-by: David Bold --- tests/integrated/test-laplace-hypre3d/runtest | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/integrated/test-laplace-hypre3d/runtest b/tests/integrated/test-laplace-hypre3d/runtest index f1f5950547..89acf79364 100755 --- a/tests/integrated/test-laplace-hypre3d/runtest +++ b/tests/integrated/test-laplace-hypre3d/runtest @@ -29,10 +29,10 @@ for directory, nproc in test_directories: error_max = collect("error_max", path=directory, info=False) if error_max > tolerance: - print(" => " + directory + " failed with maximum error {}".format(error_max)) + print(f" => {directory} failed with maximum error {error_max}") success = False else: - print(" => " + directory + " passed with maximum error {}".format(error_max)) + print(f" => {directory} passed with maximum error {error_max}") if success: print("=> All passed") From f880dcb62c6dba8df336bcb2e7d1d100651c1807 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Sat, 21 Jun 2025 04:21:03 +0000 Subject: [PATCH 12/12] Apply clang-format changes --- include/bout/hypre_interface.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 15ca97cb8a..365dd715e6 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -1036,12 +1036,12 @@ public: case HYPRE_SOLVER_TYPE::gmres: { HYPRE_ParCSRGMRESSetKDim(solver, options["kdim"] - .doc("Set the maximum size of the Krylov space") - .withDefault(30)); + .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)) { + .doc("Skip the evaluation and the check of the actual residual?") + .withDefault(false)) { HYPRE_GMRESSetSkipRealResidualCheck(solver, 1); } break;