OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Allow to use jacobi matrix for ILU with CPR

Browse Source 
Fix whitespace
This commit is contained in:
Tong Dong Qiu 2022-02-22 17:18:03 +01:00
parent 61f693dbaf
commit 9dec714590
12 changed files with 61 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
opm/simulators/linalg/bda/BdaBridge.cpp
View File

@ -245,7 +245,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unu
static std::vector<int> 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<BridgeMatrix, BridgeVector, block_size>::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<double*>(&(((*mat)[0][0][0][0]))), h_rows.data(), h_cols.data(), static_cast<double*>(&(b[0][0])),
bm_nnz, static_cast<double*>(&(((*blockMat)[0][0][0][0]))), bm_h_rows.data(), bm_h_cols.data(),
wellContribs, result);

8
opm/simulators/linalg/bda/BdaSolver.hpp
View File

@ -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;

6
opm/simulators/linalg/bda/FPGASolverBackend.cpp
View File

@ -252,9 +252,9 @@ SolverStatus FpgaSolverBackend<block_size>::solve_system(int N_, int nnz_, int d
template <unsigned int block_size>
SolverStatus FpgaSolverBackend<block_size>::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);

6
opm/simulators/linalg/bda/FPGASolverBackend.hpp
View File

@ -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

6
opm/simulators/linalg/bda/amgclSolverBackend.cpp
View File

@ -361,9 +361,9 @@ void amgclSolverBackend<block_size>::get_result(double *x_) {
template <unsigned int block_size>
SolverStatus amgclSolverBackend<block_size>::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);

6
opm/simulators/linalg/bda/amgclSolverBackend.hpp
View File

@ -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

6
opm/simulators/linalg/bda/cuda/cusparseSolverBackend.cu
View File

@ -502,9 +502,9 @@ SolverStatus cusparseSolverBackend<block_size>::solve_system(int N, int nnz, int
}
template <unsigned int block_size>
SolverStatus cusparseSolverBackend<block_size>::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);

6
opm/simulators/linalg/bda/cuda/cusparseSolverBackend.hpp
View File

@ -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

9
opm/simulators/linalg/bda/opencl/BILU0.cpp
View File

@ -178,7 +178,7 @@ bool BILU0<block_size>::analyze_matrix(BlockedMatrix *mat, BlockedMatrix *jacMat
rmat = std::make_shared<BlockedMatrix>(mat->Nb, mat->nnzbs, block_size);
rJacMat = std::make_shared<BlockedMatrix>(jacMat->Nb, jacMat->nnzbs, block_size);
LUmat = std::make_unique<BlockedMatrix>(*rJacMat);
Timer t_convert;
csrPatternToCsc(jacMat->colIndices, jacMat->rowPointers, CSCRowIndices.data(), CSCColPointers.data(), jacMat->Nb);
if(verbosity >= 3){
@ -193,20 +193,16 @@ bool BILU0<block_size>::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<block_size>(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<block_size>::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());

41
opm/simulators/linalg/bda/opencl/CPR.cpp
View File

@ -82,8 +82,39 @@ bool CPR<block_size>::analyze_matrix(BlockedMatrix *mat_) {
template <unsigned int block_size>
bool CPR<block_size>::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 <unsigned int block_size>
bool CPR<block_size>::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 <unsigned int block_size>
@ -106,11 +137,7 @@ bool CPR<block_size>::create_preconditioner(BlockedMatrix *mat_) {
return result;
}
template <unsigned int block_size>
bool CPR<block_size>::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);}));

1
opm/simulators/linalg/bda/opencl/openclSolverBackend.cpp
View File

@ -454,6 +454,7 @@ void openclSolverBackend<block_size>::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());

2
opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp
View File

@ -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();