Skip to content
Snippets Groups Projects
Verified Commit 5736e065 authored by Dorian Stoll's avatar Dorian Stoll
Browse files

tug: cpp: Remove alternative solver

parent 697ba15b
No related branches found
No related tags found
No related merge requests found
...@@ -123,8 +123,6 @@ createCoeffMatrix(const RowMajMat<T> &alpha, ...@@ -123,8 +123,6 @@ createCoeffMatrix(const RowMajMat<T> &alpha,
} }
} }
cm.makeCompressed(); // important for Eigen solver
return cm; return cm;
} }
...@@ -244,20 +242,6 @@ createSolutionVector(const RowMajMat<T> &concentrations, ...@@ -244,20 +242,6 @@ createSolutionVector(const RowMajMat<T> &concentrations,
return sv; return sv;
} }
// solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector
// use of EigenLU solver
template <class T>
static Eigen::VectorX<T> EigenLUAlgorithm(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b) {
Eigen::SparseLU<Eigen::SparseMatrix<T>> solver;
solver.analyzePattern(A);
solver.factorize(A);
return solver.solve(b);
}
// solver for linear equation system; A corresponds to coefficient matrix, // solver for linear equation system; A corresponds to coefficient matrix,
// b to the solution vector // b to the solution vector
// implementation of Thomas Algorithm // implementation of Thomas Algorithm
...@@ -329,9 +313,7 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A, ...@@ -329,9 +313,7 @@ static Eigen::VectorX<T> ThomasAlgorithm(Eigen::SparseMatrix<T> &A,
// BTCS solution for 2D grid // BTCS solution for 2D grid
template <class T> template <class T>
static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep, static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep) {
Eigen::VectorX<T> (*solverFunc)(Eigen::SparseMatrix<T> &A,
Eigen::VectorX<T> &b)) {
int rowMax = grid.getRow(); int rowMax = grid.getRow();
int colMax = grid.getCol(); int colMax = grid.getCol();
T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol()); T sx = timestep / (2 * grid.getDeltaCol() * grid.getDeltaCol());
...@@ -358,7 +340,7 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep, ...@@ -358,7 +340,7 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight,
bcTop, bcBottom, colMax, i, sx, sy); bcTop, bcBottom, colMax, i, sx, sy);
concentrations_t1.row(i) = solverFunc(A, b); concentrations_t1.row(i) = ThomasAlgorithm(A, b);
} }
concentrations_t1.transposeInPlace(); concentrations_t1.transposeInPlace();
...@@ -372,18 +354,7 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep, ...@@ -372,18 +354,7 @@ static void BTCS_2D(Grid<T> &grid, Boundary<T> &bc, T timestep,
b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom,
bcLeft, bcRight, rowMax, i, sy, sx); bcLeft, bcRight, rowMax, i, sy, sx);
concentrations.col(i) = solverFunc(A, b); concentrations.col(i) = ThomasAlgorithm(A, b);
}
}
// entry point for EigenLU solver; differentiate between 1D and 2D grid
template <class T>
void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep) {
if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, EigenLUAlgorithm);
} else {
throw_invalid_argument(
"Error: Only 2-dimensional grids are defined!");
} }
} }
...@@ -391,10 +362,10 @@ void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep) { ...@@ -391,10 +362,10 @@ void BTCS_LU(Grid<T> &grid, Boundary<T> &bc, T timestep) {
template <class T> template <class T>
void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep) { void BTCS_Thomas(Grid<T> &grid, Boundary<T> &bc, T timestep) {
if (grid.getDim() == 2) { if (grid.getDim() == 2) {
BTCS_2D(grid, bc, timestep, ThomasAlgorithm); BTCS_2D(grid, bc, timestep);
} else { } else {
throw_invalid_argument( throw_invalid_argument(
"Error: Only 1- and 2-dimensional grids are defined!"); "Error: Only 2-dimensional grids are defined!");
} }
} }
} // namespace tug } // namespace tug
......
...@@ -30,16 +30,6 @@ ...@@ -30,16 +30,6 @@
namespace tug { namespace tug {
/**
* @brief Enum defining the Linear Equation solvers
*
*/
enum SOLVER {
EIGEN_LU_SOLVER, /*!< EigenLU solver */
THOMAS_ALGORITHM_SOLVER /*!< Thomas Algorithm solver; more efficient for
tridiagonal matrices */
};
/** /**
* @brief Enum holding different options for .csv output. * @brief Enum holding different options for .csv output.
* *
...@@ -69,11 +59,8 @@ enum CONSOLE_OUTPUT { ...@@ -69,11 +59,8 @@ enum CONSOLE_OUTPUT {
* *
* @tparam T the type of the internal data structures for grid, boundary * @tparam T the type of the internal data structures for grid, boundary
* condition and timestep * condition and timestep
* @tparam approach Set the SLE scheme to be used
* @tparam solver Set the solver to be used
*/ */
template <class T, APPROACH approach = BTCS_APPROACH, template <class T>
SOLVER solver = THOMAS_ALGORITHM_SOLVER>
class Simulation { class Simulation {
public: public:
/** /**
...@@ -86,7 +73,6 @@ public: ...@@ -86,7 +73,6 @@ public:
* *
* @param grid Valid grid object * @param grid Valid grid object
* @param bc Valid boundary condition object * @param bc Valid boundary condition object
* @param approach Approach to solving the problem. Either FTCS or BTCS.
*/ */
Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){}; Simulation(Grid<T> &_grid, Boundary<T> &_bc) : grid(_grid), bc(_bc){};
...@@ -283,18 +269,7 @@ public: ...@@ -283,18 +269,7 @@ public:
} }
{ {
if constexpr (solver == EIGEN_LU_SOLVER) { {
for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole();
}
if (csv_output >= CSV_OUTPUT_VERBOSE) {
printConcentrationsCSV(filename);
}
BTCS_LU(this->grid, this->bc, this->timestep);
}
} else if constexpr (solver == THOMAS_ALGORITHM_SOLVER) {
for (int i = 0; i < iterations; i++) { for (int i = 0; i < iterations; i++) {
if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) { if (console_output == CONSOLE_OUTPUT_VERBOSE && i > 0) {
printConcentrationsConsole(); printConcentrationsConsole();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment