From 567318b8e9b100fa2b19782873c830eabd5e2e12 Mon Sep 17 00:00:00 2001 From: Hannes Signer Date: Thu, 17 Aug 2023 10:22:56 +0200 Subject: [PATCH] omp functionality and profiling --- examples/profiling.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++ src/BTCSv2.cpp | 37 +++++++++++++++++++++---------- 2 files changed, 76 insertions(+), 11 deletions(-) create mode 100644 examples/profiling.cpp diff --git a/examples/profiling.cpp b/examples/profiling.cpp new file mode 100644 index 0000000..a162705 --- /dev/null +++ b/examples/profiling.cpp @@ -0,0 +1,50 @@ +#include +#include +#include + +int main(int argc, char *argv[]) { + + int n[4] = {10, 20, 50, 100}; + int iterations[1] = {100}; + int repetition = 10; + + ofstream myfile; + myfile.open("btcs_time_measurement_openmp_10.csv"); + + for (int i = 0; i < size(n); i++){ + cout << "Grid size: " << n[i] << " x " << n[i] << endl << endl; + myfile << "Grid size: " << n[i] << " x " << n[i] << endl; + for(int j = 0; j < size(iterations); j++){ + cout << "Iterations: " << iterations[j] << endl; + for (int k = 0; k < repetition; k++){ + cout << "Wiederholung: " << k << endl; + Grid grid = Grid(n[i], n[i]); + grid.setDomain(n[i], n[i]); + + MatrixXd concentrations = MatrixXd::Constant(n[i], n[i], 0); + concentrations(5,5) = 1; + grid.setConcentrations(concentrations); + MatrixXd alpha = MatrixXd::Constant(n[i], n[i], 0.5); + + Boundary bc = Boundary(grid); + + Simulation sim = Simulation(grid, bc, BTCS_APPROACH); + + + sim.setTimestep(0.001); + sim.setIterations(iterations[j]); + sim.setOutputCSV(CSV_OUTPUT_OFF); + + auto begin = std::chrono::high_resolution_clock::now(); + sim.run(); + auto end = std::chrono::high_resolution_clock::now(); + auto milliseconds = std::chrono::duration_cast(end - begin); + myfile << milliseconds.count() << endl; + } + } + cout << endl; + myfile << endl; + + } + myfile.close(); +} \ No newline at end of file diff --git a/src/BTCSv2.cpp b/src/BTCSv2.cpp index 8b7ed21..775520d 100644 --- a/src/BTCSv2.cpp +++ b/src/BTCSv2.cpp @@ -7,6 +7,9 @@ #include "FTCS.cpp" #include +#include + +#define NUM_THREADS_BTCS 1 using namespace Eigen; @@ -62,7 +65,7 @@ static tuple calcRightBoundaryCoeffClosed(MatrixXd &alpha, int & // creates coefficient matrix for next time step from alphas in x-direction -static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector bcLeft, vector bcRight, int numCols, int rowIndex, double sx) { +static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector &bcLeft, vector &bcRight, int &numCols, int &rowIndex, double &sx) { // square matrix of column^2 dimension for the coefficients SparseMatrix cm(numCols, numCols); @@ -84,6 +87,7 @@ static SparseMatrix createCoeffMatrix(MatrixXd &alpha, vector &bcLeft, vector &bcRight, vector &bcTop, vector &bcBottom, - int length, int rowIndex, double sx, double sy) { + int &length, int &rowIndex, double &sx, double &sy) { VectorXd sv(length); int numRows = concentrations.rows(); @@ -198,6 +202,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // inner rows if (rowIndex > 0 && rowIndex < numRows-1) { + #pragma omp parallel for num_threads(NUM_THREADS_BTCS) for (int i = 0; i < length; i++) { sv(i) = sy * calcAlphaIntercell(alphaY(rowIndex,i), alphaY(rowIndex+1,i)) * concentrations(rowIndex+1,i) @@ -215,6 +220,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // first row if (rowIndex == 0) { + #pragma omp parallel for num_threads(NUM_THREADS_BTCS) for (int i = 0; i < length; i++) { type = bcTop[i].getType(); if (type == BC_TYPE_CONSTANT) { @@ -229,6 +235,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // last row if (rowIndex == numRows-1) { + #pragma omp parallel for num_threads(NUM_THREADS_BTCS) for (int i = 0; i < length; i++) { type = bcBottom[i].getType(); if (type == BC_TYPE_CONSTANT) { @@ -258,6 +265,7 @@ static VectorXd createSolutionVector(MatrixXd &concentrations, MatrixXd &alphaX, // solver for linear equation system; A corresponds to coefficient matrix, // b to the solution vector static VectorXd solve(SparseMatrix &A, VectorXd &b) { + SparseLU> solver; solver.analyzePattern(A); solver.factorize(A); @@ -281,7 +289,8 @@ static void BTCS_1D(Grid &grid, Boundary &bc, double ×tep) { vector bcRight = bc.getBoundarySide(BC_SIDE_RIGHT); MatrixXd concentrations = grid.getConcentrations(); - A = createCoeffMatrix(alpha, bcLeft, bcRight, length, 0, sx); // this is exactly same as in 2D + int rowIndex = 0; + A = createCoeffMatrix(alpha, bcLeft, bcRight, length, rowIndex, sx); // this is exactly same as in 2D for (int i = 0; i < length; i++) { b(i) = concentrations(0,i); } @@ -323,33 +332,39 @@ static void BTCS_2D(Grid &grid, Boundary &bc, double ×tep) { vector bcBottom = bc.getBoundarySide(BC_SIDE_BOTTOM); MatrixXd concentrations = grid.getConcentrations(); + #pragma omp parallel for num_threads(NUM_THREADS_BTCS) private(A, b, row_t1) + { for (int i = 0; i < rowMax; i++) { + A = createCoeffMatrix(alphaX, bcLeft, bcRight, colMax, i, sx); b = createSolutionVector(concentrations, alphaX, alphaY, bcLeft, bcRight, bcTop, bcBottom, colMax, i, sx, sy); + + SparseLU> solver; + row_t1 = solve(A, b); - for (int j = 0; j < colMax; j++) { - concentrations_t1(i,j) = row_t1(j); // can potentially be improved by using Eigen method - } - + concentrations_t1.row(i) = row_t1; + } } concentrations_t1.transposeInPlace(); concentrations.transposeInPlace(); alphaX.transposeInPlace(); alphaY.transposeInPlace(); + + #pragma omp parallel for num_threads(NUM_THREADS_BTCS) private(A, b, row_t1) + { for (int i = 0; i < colMax; i++) { - // swap alphas, boundary conditions and sx/sy for column-wise calculation A = createCoeffMatrix(alphaY, bcLeft, bcRight, rowMax, i, sy); b = createSolutionVector(concentrations_t1, alphaY, alphaX, bcTop, bcBottom, bcLeft, bcRight, rowMax, i, sy, sx); + row_t1 = solve(A, b); - for (int j = 0; j < rowMax; j++) { - concentrations(i,j) = row_t1(j); // can potentially be improved by using Eigen method - } + concentrations.row(i) = row_t1; + } } concentrations.transposeInPlace();