From c8d1b08e284c2e05053f0efa04dfb801f591a976 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 18:33:20 +0200 Subject: [PATCH 1/8] add diagonal optimization approach --- include/tug/Core/Numeric/BTCS.hpp | 117 +++++++++++++++--------------- 1 file changed, 60 insertions(+), 57 deletions(-) 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; From 06b890fe81dc2a5b71a5a7d93944660475e057cf Mon Sep 17 00:00:00 2001 From: Hannes Martin Signer Date: Tue, 14 Oct 2025 18:49:02 +0200 Subject: [PATCH 2/8] add EigenLUSolver test case --- test/testDiffusion.cpp | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 15fa18b..e89c8c2 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -114,6 +114,31 @@ DIFFUSION_TEST(EqualityBTCS) { EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01)); } +DIFFUSION_TEST(EqualityEigenLU) { + // set string from the header file + string test_path = testSimulationCSVDir; + RowMajMat reference = CSV2Eigen(test_path); + cout << "BTCS Test: " << endl; + + RowMajMat concentrations = MatrixXd::Constant(row, col, 0); + + Diffusion sim = + setupSimulation(concentrations, timestep, + iterations); // Boundary + + // Boundary bc = Boundary(grid); + + // Simulation + // Diffusion sim(grid, bc); + // sim.setOutputConsole(CONSOLE_OUTPUT_ON); + // sim.setTimestep(timestep); + // sim.setIterations(iterations); + sim.run(); + + cout << endl; + EXPECT_TRUE(checkSimilarityV2(reference, sim.getConcentrationMatrix(), 0.01)); +} + DIFFUSION_TEST(InitializeEnvironment) { int rc = 12; RowMajMat concentrations(rc, rc); From 91ae02cbf1d8304e9b2377fc1298c5a6ce51c73f Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 18:54:18 +0200 Subject: [PATCH 3/8] fix error for missing matching function --- include/tug/Core/Numeric/BTCS.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/tug/Core/Numeric/BTCS.hpp b/include/tug/Core/Numeric/BTCS.hpp index 3a875ec..0d2375a 100644 --- a/include/tug/Core/Numeric/BTCS.hpp +++ b/include/tug/Core/Numeric/BTCS.hpp @@ -464,7 +464,7 @@ template void BTCS_LU(SimulationInput &input, int numThreads) { if (input.dim == 1) { BTCS_1D(input, EigenLUAlgorithm); } else { - BTCS_2D(input.dim, EigenLUAlgorithm, numThreads); + BTCS_2D(input, EigenLUAlgorithm, numThreads); } } From 30d2099d8a3e236811d513c9b88fb02c2731e1a1 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 19:01:32 +0200 Subject: [PATCH 4/8] add solver to template --- test/testDiffusion.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index e89c8c2..d1c8225 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -20,7 +20,7 @@ using namespace tug; constexpr int row = 11; constexpr int col = 11; -template +template Diffusion setupSimulation(RowMajMat &concentrations, double timestep, int iterations) { int domain_row = 10; @@ -30,7 +30,7 @@ Diffusion setupSimulation(RowMajMat &concentrations, // RowMajMat concentrations = MatrixXd::Constant(row, col, 0); concentrations(5, 5) = 1; - Diffusion diffusiongrid(concentrations); + Diffusion diffusiongrid(concentrations); diffusiongrid.getConcentrationMatrix() = concentrations; diffusiongrid.setDomain(domain_row, domain_col); @@ -72,7 +72,7 @@ DIFFUSION_TEST(EqualityFTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); - Diffusion sim = + Diffusion sim = setupSimulation(concentrations, timestep, iterations); // Boundary bc = Boundary(grid); @@ -97,7 +97,7 @@ DIFFUSION_TEST(EqualityBTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); - Diffusion sim = + Diffusion sim = setupSimulation(concentrations, timestep, iterations); // Boundary From 3701dea34ec008d89d49d900cf00fcef0e5fa7e0 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 19:05:25 +0200 Subject: [PATCH 5/8] fix test case --- test/testDiffusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index d1c8225..edb9749 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -72,7 +72,7 @@ DIFFUSION_TEST(EqualityFTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); - Diffusion sim = + Diffusion sim = setupSimulation(concentrations, timestep, iterations); // Boundary bc = Boundary(grid); From a20d5d53e64c5ea5504aa4b84e8f7affdd2fe997 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 19:08:16 +0200 Subject: [PATCH 6/8] fix test case --- test/testDiffusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index edb9749..d1c8225 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -72,7 +72,7 @@ DIFFUSION_TEST(EqualityFTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); - Diffusion sim = + Diffusion sim = setupSimulation(concentrations, timestep, iterations); // Boundary bc = Boundary(grid); From 3e270ee01c7b7a896fc7f4ad53bb051b0a4fdf28 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 19:10:03 +0200 Subject: [PATCH 7/8] add solver --- test/testDiffusion.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index d1c8225..9530795 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -21,8 +21,9 @@ constexpr int row = 11; constexpr int col = 11; template -Diffusion setupSimulation(RowMajMat &concentrations, - double timestep, int iterations) { +Diffusion +setupSimulation(RowMajMat &concentrations, double timestep, + int iterations) { int domain_row = 10; int domain_col = 10; From 5b144aea3adc09ab0d742dacdaeea51e98016738 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Tue, 14 Oct 2025 19:12:02 +0200 Subject: [PATCH 8/8] add solver --- test/testDiffusion.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/testDiffusion.cpp b/test/testDiffusion.cpp index 9530795..1a7c11a 100644 --- a/test/testDiffusion.cpp +++ b/test/testDiffusion.cpp @@ -74,7 +74,8 @@ DIFFUSION_TEST(EqualityFTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); Diffusion sim = - setupSimulation(concentrations, timestep, iterations); + setupSimulation( + concentrations, timestep, iterations); // Boundary bc = Boundary(grid); @@ -99,8 +100,9 @@ DIFFUSION_TEST(EqualityBTCS) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); Diffusion sim = - setupSimulation(concentrations, timestep, - iterations); // Boundary + setupSimulation( + concentrations, timestep, + iterations); // Boundary // Boundary bc = Boundary(grid); @@ -124,8 +126,9 @@ DIFFUSION_TEST(EqualityEigenLU) { RowMajMat concentrations = MatrixXd::Constant(row, col, 0); Diffusion sim = - setupSimulation(concentrations, timestep, - iterations); // Boundary + setupSimulation( + concentrations, timestep, + iterations); // Boundary // Boundary bc = Boundary(grid);