diff --git a/Code/Source/linear_solver/precond.cpp b/Code/Source/linear_solver/precond.cpp index 37be32c03..ae1f3a4d8 100644 --- a/Code/Source/linear_solver/precond.cpp +++ b/Code/Source/linear_solver/precond.cpp @@ -73,7 +73,7 @@ void pos_mul(const Array& rowPtr, const Vector& colPtr, const int nNo, int a = colPtr(i); for (int b = 0; b < dof; b++) { int j = dof*(dof-1) + b; - for (int k = b; k < j; k += dof) { + for (int k = b; k <= j; k += dof) { Val(k,i) = Val(k,i)*W(b,a); } } @@ -303,52 +303,16 @@ void precond_rcs(fsi_linear_solver::FSILS_lhsType& lhs, const Array& rowPtr pos_mul(rowPtr, colPtr, lhs.nNo, lhs.nnz, dof, Val, Wr); - // Set diagonal term to one + // Preserve free diagonals: subtracting and adding one loses small values. // - switch (dof) { - case 1: - for (int Ac = 0; Ac < nNo; Ac++) { - int d = diagPtr(Ac); - Val(0,d) = Wr(0,Ac) * (Val(0,d) - 1.0) + 1.0; - } - break; - - case 2: - for (int Ac = 0; Ac < nNo; Ac++) { - int d = diagPtr(Ac); - Val(0,d) = Wr(0,Ac)*(Val(0,d)-1.0) + 1.0; - Val(3,d) = Wr(1,Ac)*(Val(3,d)-1.0) + 1.0; - } - break; - - case 3: - for (int Ac = 0; Ac < nNo; Ac++) { - int d = diagPtr(Ac); - Val(0,d) = Wr(0,Ac)*(Val(0,d)-1.0) + 1.0; - Val(4,d) = Wr(1,Ac)*(Val(4,d)-1.0) + 1.0; - Val(8,d) = Wr(2,Ac)*(Val(8,d)-1.0) + 1.0; + for (int Ac = 0; Ac < nNo; ++Ac) { + const int d = diagPtr(Ac); + for (int i = 0; i < dof; ++i) { + if (Wr(i, Ac) == 0.0) { + Val(i * dof + i, d) = 1.0; } - break; - - case 4: - for (int Ac = 0; Ac < nNo; Ac++) { - int d = diagPtr(Ac); - Val(0 ,d) = Wr(0,Ac)*(Val(0 ,d)-1.0) + 1.0; - Val(5 ,d) = Wr(1,Ac)*(Val(5 ,d)-1.0) + 1.0; - Val(10,d) = Wr(2,Ac)*(Val(10,d)-1.0) + 1.0; - Val(15,d) = Wr(3,Ac)*(Val(15,d)-1.0) + 1.0; - } - break; - - default: - for (int Ac = 0; Ac < nNo; Ac++) { - int d = diagPtr(Ac); - for (int i = 0; i < dof; i++) { - Val(i*dof+i,d) = Wr(i,Ac)*(Val(i*dof+i,d) - 1.0) + 1.0; - } - } - break; - } + } + } //***************************************************** // Row and column scaling @@ -475,7 +439,7 @@ void precond_rcs(fsi_linear_solver::FSILS_lhsType& lhs, const Array& rowPtr int b = rowPtr(1,Ac); for (int i = 0; i < dof; i++) { - int j = i*dof + 1; + int j = i*dof; int k = (i+1)*dof - 1; auto vals = Val.values({j,k}, {a,b}); Wr(i,Ac) = fabs(*std::max_element(vals.begin(), vals.end(), max_func)); diff --git a/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp b/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp new file mode 100644 index 000000000..a049b2a8d --- /dev/null +++ b/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp @@ -0,0 +1,136 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#include "precond.h" +#include "FE/Math/DenseLinearAlgebra.h" + +#include + +#include +#include +#include + +namespace { + +// Boundary modes: excluded face, first component fixed, all components fixed. +using RcsParameters = std::tuple; +class RcsPreconditioner : public testing::TestWithParam {}; + +TEST_P(RcsPreconditioner, PreservesConstrainedSystem) { + const auto [dof, scale, boundary] = GetParam(); + const int n = 3 * dof; + const auto fixed = [=](int row) { + return row < dof && boundary != 0 && (boundary == 2 || row == 0); + }; + + // Three-node tridiagonal stencil with coupled, positive-definite blocks. + std::vector original(n * n, 0.0), exact(n), original_rhs(n, 0.0); + for (int row = 0; row < n; ++row) { + exact[row] = fixed(row) ? 0.0 : row + 1.0; + for (int col = 0; col < n; ++col) { + const int a = row / dof; + const int b = col / dof; + if (std::abs(a - b) <= 1) { + original[row * n + col] = scale * (a == b ? 2.0 : -1.0) + * (row % dof == col % dof ? 4.0 : 0.25); + } + } + } + for (int row = 0; row < n; ++row) { + for (int col = 0; col < n; ++col) { + original_rhs[row] += original[row * n + col] * exact[col]; + } + } + + fsi_linear_solver::FSILS_lhsType lhs{}; + lhs.nNo = lhs.gnNo = lhs.mynNo = 3; + lhs.nnz = 7; + lhs.commu.nTasks = 1; + lhs.rowPtr.resize(2, 3); + lhs.colPtr = Vector{0, 1, 0, 1, 2, 1, 2}; + lhs.diagPtr = Vector{0, 3, 6}; + const int starts[] = {0, 2, 5}; + const int ends[] = {1, 4, 6}; + for (int a = 0; a < 3; ++a) { + lhs.rowPtr(0, a) = starts[a]; + lhs.rowPtr(1, a) = ends[a]; + } + lhs.nFaces = 1; + lhs.face.resize(1); + auto& face = lhs.face[0]; + face.incFlag = boundary != 0; + face.bGrp = fsi_linear_solver::BcType::BC_TYPE_Dir; + face.nNo = 1; + face.dof = dof; + face.glob = Vector{0}; + face.val.resize(dof, 1); + for (int i = 0; i < dof; ++i) { + face.val(i, 0) = (boundary == 2 || i == 0) ? 0.0 : 1.0; + } + + Array val(dof * dof, lhs.nnz), rhs(dof, 3); + Array left(dof, 3), right(dof, 3); + for (int a = 0; a < 3; ++a) { + for (int i = 0; i < dof; ++i) { + rhs(i, a) = original_rhs[a * dof + i]; + } + for (int k = starts[a]; k <= ends[a]; ++k) { + for (int i = 0; i < dof; ++i) { + for (int j = 0; j < dof; ++j) { + val(i * dof + j, k) = original[(a * dof + i) * n + lhs.colPtr(k) * dof + j]; + } + } + } + } + + precond::precond_rcs(lhs, lhs.rowPtr, lhs.colPtr, lhs.diagPtr, + dof, val, rhs, left, right); + + std::vector matrix(n * n, 0.0), solution(n); + for (int a = 0; a < 3; ++a) { + for (int i = 0; i < dof; ++i) { + ASSERT_TRUE(std::isfinite(left(i, a))); + ASSERT_TRUE(std::isfinite(right(i, a))); + ASSERT_GT(left(i, a), 0.0); + ASSERT_GT(right(i, a), 0.0); + const int row = a * dof + i; + const double expected_rhs = fixed(row) ? 0.0 : original_rhs[row]; + ASSERT_NEAR(rhs(i, a) / left(i, a) / scale, expected_rhs / scale, 1e-12); + solution[row] = rhs(i, a); + for (int k = starts[a]; k <= ends[a]; ++k) { + for (int j = 0; j < dof; ++j) { + const int col = lhs.colPtr(k) * dof + j; + const bool constrained = fixed(row) || fixed(col); + const double expected = constrained ? (row == col ? 1.0 : 0.0) + : original[row * n + col]; + const double restored = val(i * dof + j, k) / left(i, a) / right(j, lhs.colPtr(k)); + // Normalize by coefficient scale so erased small entries cannot pass. + const double reference_scale = constrained && row == col ? 1.0 : scale; + ASSERT_NEAR(restored / reference_scale, expected / reference_scale, 1e-12) + << "row=" << row << ", col=" << col; + matrix[row * n + col] = val(i * dof + j, k); + } + } + } + } + + svmp::FE::math::factor_dense_matrix(matrix, n, "RCS regression").solve_in_place(solution); + for (int row = 0; row < n; ++row) { + solution[row] *= right(row % dof, row / dof); + EXPECT_NEAR(solution[row], exact[row], 1e-11); + } + for (int row = 0; row < n; ++row) { + if (fixed(row)) continue; + double residual = -original_rhs[row]; + for (int col = 0; col < n; ++col) { + residual += original[row * n + col] * solution[col]; + } + EXPECT_NEAR(residual / scale, 0.0, 1e-11); + } +} + +INSTANTIATE_TEST_SUITE_P(Blocks, RcsPreconditioner, + testing::Combine(testing::Values(1, 2, 3, 4, 5), testing::Values(1.0, 1e-11, 1e-15, 1e-17, 1e-20), + testing::Values(0, 1, 2))); + +} // namespace