diff --git a/include/tug/Core/Numeric/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp index 5e561c5..3a875ec 100644 --- a/include/tug/Core/Numeric/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -28,6 +28,18 @@ namespace tug { +// optimization to remove Eigen sparse matrix +template class Diagonals { +public: + Diagonals() : left(), center(), right() {}; + Diagonals(std::size_t size) : left(size), center(size), right(size) {}; + +public: + std::vector left; + std::vector center; + std::vector right; +}; + // calculates coefficient for boundary in constant case template constexpr std::pair calcBoundaryCoeffConstant(T alpha_center, @@ -52,7 +64,7 @@ constexpr std::pair calcBoundaryCoeffClosed(T alpha_center, T alpha_side, // creates coefficient matrix for next time step from alphas in x-direction template -static Eigen::SparseMatrix +static Diagonals createCoeffMatrix(const RowMajMat &alpha, const std::vector> &bcLeft, const std::vector> &bcRight, @@ -60,26 +72,25 @@ createCoeffMatrix(const RowMajMat &alpha, int rowIndex, T sx) { // square matrix of column^2 dimension for the coefficients - Eigen::SparseMatrix cm(numCols, numCols); - cm.reserve(Eigen::VectorXi::Constant(numCols, 3)); + Diagonals cm(numCols); // left column if (inner_bc[0].first) { - cm.insert(0, 0) = 1; + cm.center[0] = 1; } else { switch (bcLeft[rowIndex].getType()) { case BC_TYPE_CONSTANT: { auto [centerCoeffTop, rightCoeffTop] = calcBoundaryCoeffConstant(alpha(rowIndex, 0), alpha(rowIndex, 1), sx); - cm.insert(0, 0) = centerCoeffTop; - cm.insert(0, 1) = rightCoeffTop; + cm.center[0] = centerCoeffTop; + cm.right[0] = rightCoeffTop; break; } case BC_TYPE_CLOSED: { auto [centerCoeffTop, rightCoeffTop] = calcBoundaryCoeffClosed(alpha(rowIndex, 0), alpha(rowIndex, 1), sx); - cm.insert(0, 0) = centerCoeffTop; - cm.insert(0, 1) = rightCoeffTop; + cm.center[0] = centerCoeffTop; + cm.right[0] = rightCoeffTop; break; } default: { @@ -93,36 +104,36 @@ createCoeffMatrix(const RowMajMat &alpha, int n = numCols - 1; for (int i = 1; i < n; i++) { if (inner_bc[i].first) { - cm.insert(i, i) = 1; + cm.center[i] = 1; continue; } - cm.insert(i, i - 1) = + cm.left[i] = -sx * calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i)); - cm.insert(i, i) = + cm.center[i] = 1 + sx * (calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)) + calcAlphaIntercell(alpha(rowIndex, i - 1), alpha(rowIndex, i))); - cm.insert(i, i + 1) = + cm.right[i] = -sx * calcAlphaIntercell(alpha(rowIndex, i), alpha(rowIndex, i + 1)); } // right column if (inner_bc[n].first) { - cm.insert(n, n) = 1; + cm.center[n] = 1; } else { switch (bcRight[rowIndex].getType()) { case BC_TYPE_CONSTANT: { auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffConstant( alpha(rowIndex, n), alpha(rowIndex, n - 1), sx); - cm.insert(n, n - 1) = leftCoeffBottom; - cm.insert(n, n) = centerCoeffBottom; + cm.left[n] = leftCoeffBottom; + cm.center[n] = centerCoeffBottom; break; } case BC_TYPE_CLOSED: { auto [centerCoeffBottom, leftCoeffBottom] = calcBoundaryCoeffClosed( alpha(rowIndex, n), alpha(rowIndex, n - 1), sx); - cm.insert(n, n - 1) = leftCoeffBottom; - cm.insert(n, n) = centerCoeffBottom; + cm.left[n] = leftCoeffBottom; + cm.center[n] = centerCoeffBottom; break; } default: { @@ -132,8 +143,6 @@ createCoeffMatrix(const RowMajMat &alpha, } } - cm.makeCompressed(); // important for Eigen solver - return cm; } @@ -271,12 +280,28 @@ createSolutionVector(const EigenType &concentrations, // b to the solution vector // use of EigenLU solver template -static Eigen::VectorX EigenLUAlgorithm(Eigen::SparseMatrix &A, +static Eigen::VectorX EigenLUAlgorithm(Diagonals &A, Eigen::VectorX &b) { + // convert A to Eigen sparse matrix + size_t dim_A = A.center.size(); + Eigen::SparseMatrix A_sparse(dim_A, dim_A); + + A_sparse.insert(0, 0) = A.center[0]; + A_sparse.insert(0, 1) = A.right[0]; + + for (size_t i = 1; i < dim_A - 1; i++) { + A_sparse.insert(i, i - 1) = A.left[i]; + A_sparse.insert(i, i) = A.center[i]; + A_sparse.insert(i, i + 1) = A.right[i]; + } + + A_sparse.insert(dim_A - 1, dim_A - 2) = A.left[dim_A - 1]; + A_sparse.insert(dim_A - 1, dim_A - 1) = A.center[dim_A - 1]; + Eigen::SparseLU> solver; - solver.analyzePattern(A); - solver.factorize(A); + solver.analyzePattern(A_sparse); + solver.factorize(A_sparse); return solver.solve(b); } @@ -285,27 +310,11 @@ static Eigen::VectorX EigenLUAlgorithm(Eigen::SparseMatrix &A, // b to the solution vector // implementation of Thomas Algorithm template -static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, +static Eigen::VectorX ThomasAlgorithm(Diagonals &A, Eigen::VectorX &b) { Eigen::Index n = b.size(); - - Eigen::VectorX a_diag(n); - Eigen::VectorX b_diag(n); - Eigen::VectorX c_diag(n); Eigen::VectorX x_vec = b; - // Fill diagonals vectors - b_diag[0] = A.coeff(0, 0); - c_diag[0] = A.coeff(0, 1); - - for (Eigen::Index i = 1; i < n - 1; i++) { - a_diag[i] = A.coeff(i, i - 1); - b_diag[i] = A.coeff(i, i); - c_diag[i] = A.coeff(i, i + 1); - } - a_diag[n - 1] = A.coeff(n - 1, n - 2); - b_diag[n - 1] = A.coeff(n - 1, n - 1); - // HACK: write CSV to file #ifdef WRITE_THOMAS_CSV #include @@ -331,20 +340,20 @@ static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, // start solving - c_diag and x_vec are overwritten n--; - c_diag[0] /= b_diag[0]; - x_vec[0] /= b_diag[0]; + A.right[0] /= A.center[0]; + x_vec[0] /= A.center[0]; for (Eigen::Index i = 1; i < n; i++) { - c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i - 1]; - x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i - 1]) / - (b_diag[i] - a_diag[i] * c_diag[i - 1]); + A.right[i] /= A.center[i] - A.left[i] * A.right[i - 1]; + x_vec[i] = (x_vec[i] - A.left[i] * x_vec[i - 1]) / + (A.center[i] - A.left[i] * A.right[i - 1]); } - x_vec[n] = (x_vec[n] - a_diag[n] * x_vec[n - 1]) / - (b_diag[n] - a_diag[n] * c_diag[n - 1]); + x_vec[n] = (x_vec[n] - A.left[n] * x_vec[n - 1]) / + (A.center[n] - A.left[n] * A.right[n - 1]); for (Eigen::Index i = n; i-- > 0;) { - x_vec[i] -= c_diag[i] * x_vec[i + 1]; + x_vec[i] -= A.right[i] * x_vec[i + 1]; } return x_vec; @@ -353,14 +362,14 @@ static Eigen::VectorX ThomasAlgorithm(Eigen::SparseMatrix &A, // BTCS solution for 1D grid template static void BTCS_1D(SimulationInput &input, - Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX (*solverFunc)(Diagonals &A, Eigen::VectorX &b)) { const std::size_t &length = input.colMax; T sx = input.timestep / (input.deltaCol * input.deltaCol); Eigen::VectorX concentrations_t1(length); - Eigen::SparseMatrix A; + Diagonals A; Eigen::VectorX b(length); const auto &alpha = input.alphaX; @@ -389,18 +398,12 @@ static void BTCS_1D(SimulationInput &input, } concentrations = solverFunc(A, b); - - // for (int j = 0; j < length; j++) { - // concentrations(0, j) = concentrations_t1(j); - // } - - // grid.setConcentrations(concentrations); } // BTCS solution for 2D grid template static void BTCS_2D(SimulationInput &input, - Eigen::VectorX (*solverFunc)(Eigen::SparseMatrix &A, + Eigen::VectorX (*solverFunc)(Diagonals &A, Eigen::VectorX &b), int numThreads) { const std::size_t &rowMax = input.rowMax; @@ -410,7 +413,7 @@ static void BTCS_2D(SimulationInput &input, RowMajMat concentrations_t1(rowMax, colMax); - Eigen::SparseMatrix A; + Diagonals A; Eigen::VectorX b; const RowMajMat &alphaX = input.alphaX;