From 9dec7145909d62519a23c75e4c6e894d3f7cdc2e Mon Sep 17 00:00:00 2001 From: Tong Dong Qiu Date: Tue, 22 Feb 2022 17:18:03 +0100 Subject: [PATCH] Allow to use jacobi matrix for ILU with CPR Fix whitespace --- opm/simulators/linalg/bda/BdaBridge.cpp | 3 +- opm/simulators/linalg/bda/BdaSolver.hpp | 8 ++-- .../linalg/bda/FPGASolverBackend.cpp | 6 +-- .../linalg/bda/FPGASolverBackend.hpp | 6 +-- .../linalg/bda/amgclSolverBackend.cpp | 6 +-- .../linalg/bda/amgclSolverBackend.hpp | 6 +-- .../linalg/bda/cuda/cusparseSolverBackend.cu | 6 +-- .../linalg/bda/cuda/cusparseSolverBackend.hpp | 6 +-- opm/simulators/linalg/bda/opencl/BILU0.cpp | 9 +--- opm/simulators/linalg/bda/opencl/CPR.cpp | 41 +++++++++++++++---- .../linalg/bda/opencl/openclSolverBackend.cpp | 1 + .../linalg/bda/opencl/openclSolverBackend.hpp | 2 +- 12 files changed, 61 insertions(+), 39 deletions(-) diff --git a/opm/simulators/linalg/bda/BdaBridge.cpp b/opm/simulators/linalg/bda/BdaBridge.cpp index 58f229606..fb920eadb 100644 --- a/opm/simulators/linalg/bda/BdaBridge.cpp +++ b/opm/simulators/linalg/bda/BdaBridge.cpp @@ -245,7 +245,7 @@ void BdaBridge::solve_system([[maybe_unu static std::vector bm_h_cols; const int bm_nnzb = (bm_h_rows.empty()) ? blockMat->nonzeroes() : bm_h_rows.back(); const int bm_nnz = bm_nnzb * dim * dim; - + if (bm_h_rows.capacity() == 0) { bm_h_rows.reserve(Nb+1); bm_h_cols.reserve(bm_nnzb); @@ -258,7 +258,6 @@ void BdaBridge::solve_system([[maybe_unu out << "Checking zeros took: " << t_zeros.stop() << " s, found " << numZeros << " zeros"; OpmLog::info(out.str()); } - status = backend->solve_system2(N, nnz, dim, static_cast(&(((*mat)[0][0][0][0]))), h_rows.data(), h_cols.data(), static_cast(&(b[0][0])), bm_nnz, static_cast(&(((*blockMat)[0][0][0][0]))), bm_h_rows.data(), bm_h_cols.data(), wellContribs, result); diff --git a/opm/simulators/linalg/bda/BdaSolver.hpp b/opm/simulators/linalg/bda/BdaSolver.hpp index 7a65d458b..799cd0785 100644 --- a/opm/simulators/linalg/bda/BdaSolver.hpp +++ b/opm/simulators/linalg/bda/BdaSolver.hpp @@ -89,10 +89,10 @@ namespace Accelerator { double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) = 0; - virtual SolverStatus solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) = 0; + virtual SolverStatus solve_system2(int N_, int nnz_, int dim, + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) = 0; virtual void get_result(double *x) = 0; diff --git a/opm/simulators/linalg/bda/FPGASolverBackend.cpp b/opm/simulators/linalg/bda/FPGASolverBackend.cpp index a20b406c8..4ab683582 100644 --- a/opm/simulators/linalg/bda/FPGASolverBackend.cpp +++ b/opm/simulators/linalg/bda/FPGASolverBackend.cpp @@ -252,9 +252,9 @@ SolverStatus FpgaSolverBackend::solve_system(int N_, int nnz_, int d template SolverStatus FpgaSolverBackend::solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) { (void) nnz2; (void) vals2; (void) rows2; (void) cols2; return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res); diff --git a/opm/simulators/linalg/bda/FPGASolverBackend.hpp b/opm/simulators/linalg/bda/FPGASolverBackend.hpp index 12293e2e7..ee3e05e25 100644 --- a/opm/simulators/linalg/bda/FPGASolverBackend.hpp +++ b/opm/simulators/linalg/bda/FPGASolverBackend.hpp @@ -256,9 +256,9 @@ public: SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override; SolverStatus solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) override; + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) override; /// Get result after linear solve, and peform postprocessing if necessary /// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array diff --git a/opm/simulators/linalg/bda/amgclSolverBackend.cpp b/opm/simulators/linalg/bda/amgclSolverBackend.cpp index 5bc158572..f05c7b089 100644 --- a/opm/simulators/linalg/bda/amgclSolverBackend.cpp +++ b/opm/simulators/linalg/bda/amgclSolverBackend.cpp @@ -361,9 +361,9 @@ void amgclSolverBackend::get_result(double *x_) { template SolverStatus amgclSolverBackend::solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) { (void) nnz2; (void) vals2; (void) rows2; (void) cols2; return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res); diff --git a/opm/simulators/linalg/bda/amgclSolverBackend.hpp b/opm/simulators/linalg/bda/amgclSolverBackend.hpp index 95f3c928f..9ceb1fe9c 100644 --- a/opm/simulators/linalg/bda/amgclSolverBackend.hpp +++ b/opm/simulators/linalg/bda/amgclSolverBackend.hpp @@ -146,9 +146,9 @@ public: SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override; SolverStatus solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) override; + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) override; /// Get result after linear solve, and peform postprocessing if necessary /// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array diff --git a/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.cu b/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.cu index ab543c3b4..aec915976 100644 --- a/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.cu +++ b/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.cu @@ -502,9 +502,9 @@ SolverStatus cusparseSolverBackend::solve_system(int N, int nnz, int } template SolverStatus cusparseSolverBackend::solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) { (void) nnz2; (void) vals2; (void) rows2; (void) cols2; return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res); diff --git a/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.hpp b/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.hpp index 8b05e1a8a..e215c8d2d 100644 --- a/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.hpp +++ b/opm/simulators/linalg/bda/cuda/cusparseSolverBackend.hpp @@ -139,9 +139,9 @@ public: SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override; SolverStatus solve_system2(int N_, int nnz_, int dim, - double *vals, int *rows, int *cols, double *b, - int nnz2, double *vals2, int *rows2, int *cols2, - WellContributions& wellContribs, BdaResult &res) override; + double *vals, int *rows, int *cols, double *b, + int nnz2, double *vals2, int *rows2, int *cols2, + WellContributions& wellContribs, BdaResult &res) override; /// Get resulting vector x after linear solve, also includes post processing if necessary /// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array diff --git a/opm/simulators/linalg/bda/opencl/BILU0.cpp b/opm/simulators/linalg/bda/opencl/BILU0.cpp index ccbd5df84..a73d3d66e 100644 --- a/opm/simulators/linalg/bda/opencl/BILU0.cpp +++ b/opm/simulators/linalg/bda/opencl/BILU0.cpp @@ -178,7 +178,7 @@ bool BILU0::analyze_matrix(BlockedMatrix *mat, BlockedMatrix *jacMat rmat = std::make_shared(mat->Nb, mat->nnzbs, block_size); rJacMat = std::make_shared(jacMat->Nb, jacMat->nnzbs, block_size); LUmat = std::make_unique(*rJacMat); - + Timer t_convert; csrPatternToCsc(jacMat->colIndices, jacMat->rowPointers, CSCRowIndices.data(), CSCColPointers.data(), jacMat->Nb); if(verbosity >= 3){ @@ -193,20 +193,16 @@ bool BILU0::analyze_matrix(BlockedMatrix *mat, BlockedMatrix *jacMat if (opencl_ilu_reorder == ILUReorder::LEVEL_SCHEDULING) { out << "BILU0 reordering strategy: " << "level_scheduling\n"; findLevelScheduling(jacMat->colIndices, jacMat->rowPointers, CSCRowIndices.data(), CSCColPointers.data(), jacMat->Nb, &numColors, toOrder.data(), fromOrder.data(), rowsPerColor); - for (int iii = 0; iii < numColors; ++iii) { out << "rpc: "<< rowsPerColor[iii] << "\n";} - out << "numColors: "<< numColors << "\n"; } else if (opencl_ilu_reorder == ILUReorder::GRAPH_COLORING) { out << "BILU0 reordering strategy: " << "graph_coloring\n"; findGraphColoring(jacMat->colIndices, jacMat->rowPointers, CSCRowIndices.data(), CSCColPointers.data(), jacMat->Nb, jacMat->Nb, jacMat->Nb, &numColors, toOrder.data(), fromOrder.data(), rowsPerColor); - for (int iii = 0; iii < numColors; ++iii) { out << "rpc: "<< rowsPerColor[iii] << "\n";} - out << "numColors: "<< numColors << "\n"; } else if (opencl_ilu_reorder == ILUReorder::NONE) { out << "BILU0 reordering strategy: none\n"; // numColors = 1; // rowsPerColor.emplace_back(Nb); numColors = Nb; for(int i = 0; i < Nb; ++i){ - rowsPerColor.emplace_back(1); + rowsPerColor.emplace_back(1); } } else { OPM_THROW(std::logic_error, "Error ilu reordering strategy not set correctly\n"); @@ -369,7 +365,6 @@ bool BILU0::create_preconditioner(BlockedMatrix *mat, BlockedMatrix auto *jm = jacMat; if (opencl_ilu_reorder != ILUReorder::NONE) { - jm = rJacMat.get(); Timer t_reorder; reorderBlockedMatrixByPattern(mat, toOrder.data(), fromOrder.data(), rmat.get()); diff --git a/opm/simulators/linalg/bda/opencl/CPR.cpp b/opm/simulators/linalg/bda/opencl/CPR.cpp index 95973207e..9ae847c99 100644 --- a/opm/simulators/linalg/bda/opencl/CPR.cpp +++ b/opm/simulators/linalg/bda/opencl/CPR.cpp @@ -82,8 +82,39 @@ bool CPR::analyze_matrix(BlockedMatrix *mat_) { template bool CPR::analyze_matrix(BlockedMatrix *mat_, BlockedMatrix *jacMat) { - (void) jacMat; - return analyze_matrix(mat_); + this->Nb = mat_->Nb; + this->nnzb = mat_->nnzbs; + this->N = Nb * block_size; + this->nnz = nnzb * block_size * block_size; + + bool success = bilu0->analyze_matrix(mat_, jacMat); + + if (opencl_ilu_reorder == ILUReorder::NONE) { + mat = mat_; + } else { + mat = bilu0->getRMat(); + } + return success; +} + +template +bool CPR::create_preconditioner(BlockedMatrix *mat_, BlockedMatrix *jacMat) { + Dune::Timer t_bilu0; + bool result = bilu0->create_preconditioner(mat_, jacMat); + if (verbosity >= 3) { + std::ostringstream out; + out << "CPR create_preconditioner bilu0(): " << t_bilu0.stop() << " s"; + OpmLog::info(out.str()); + } + + Dune::Timer t_amg; + create_preconditioner_amg(mat); // already points to bilu0::rmat if needed + if (verbosity >= 3) { + std::ostringstream out; + out << "CPR create_preconditioner_amg(): " << t_amg.stop() << " s"; + OpmLog::info(out.str()); + } + return result; } template @@ -106,11 +137,7 @@ bool CPR::create_preconditioner(BlockedMatrix *mat_) { return result; } -template -bool CPR::create_preconditioner(BlockedMatrix *mat_, BlockedMatrix *jacMat) { - (void) jacMat; - return create_preconditioner(mat_); -} + // return the absolute value of the N elements for which the absolute value is highest double get_absmax(const double *data, const int N) { return std::abs(*std::max_element(data, data + N, [](double a, double b){return std::fabs(a) < std::fabs(b);})); diff --git a/opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp b/opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp index fbf6e65e4..22c6b0e2f 100644 --- a/opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp +++ b/opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp @@ -454,6 +454,7 @@ void openclSolverBackend::initialize2(int N_, int nnz_, int dim, dou Nb = (N + dim - 1) / dim; std::ostringstream out; out << "Initializing GPU, matrix size: " << N << " blocks, nnzb: " << nnzb << "\n"; + out << "Blocks in ILU matrix: " << jac_nnzb << "\n"; out << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance << "\n"; out << "PlatformID: " << platformID << ", deviceID: " << deviceID << "\n"; OpmLog::info(out.str()); diff --git a/opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp b/opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp index f32eaeb9d..e04375808 100644 --- a/opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp +++ b/opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp @@ -144,7 +144,7 @@ private: void initialize(int N, int nnz, int dim, double *vals, int *rows, int *cols); void initialize2(int N, int nnz, int dim, double *vals, int *rows, int *cols, - int nnz2, double *vals2, int *rows2, int *cols2); + int nnz2, double *vals2, int *rows2, int *cols2); /// Clean memory void finalize();