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

Pass BlockedMatrix to BdaSolvers

Browse Source
This commit is contained in:
Tong Dong Qiu
2022-02-07 13:39:38 +01:00
parent 7f5322f7d4
commit d963820e48
10 changed files with 119 additions and 164 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -205,9 +205,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unu
result.converged = false;
const int dim = (*bridgeMat)[0][0].N();
const int Nb = bridgeMat->N();
const int N = Nb * dim;
const int nnzb = bridgeMat->nonzeroes();
const int nnz = nnzb * dim * dim;
if (dim != 3) {
OpmLog::warning("BdaSolver only accepts blocksize = 3 at this time, will use Dune for the remainder of the program");
@@ -237,10 +235,9 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unu
SolverStatus status;
if (numJacobiBlocks < 2) {
// assume that underlying data (nonzeroes) from mat (Dune::BCRSMatrix) are contiguous, if this is not the case, the chosen BdaSolver is expected to perform undefined behaviour
status = backend->solve_system(N, nnz, dim, static_cast<double*>(&(((*bridgeMat)[0][0][0][0]))), h_rows.data(), h_cols.data(), static_cast<double*>(&(b[0][0])), wellContribs, result);
status = backend->solve_system(matrix, static_cast<double*>(&(b[0][0])), wellContribs, result);
} else {
const int jacNnzb = (h_jacRows.empty()) ? jacMat->nonzeroes() : h_jacRows.back();
const int jacNnz = jacNnzb * dim * dim;
if (!jacMatrix) {
h_jacRows.reserve(Nb+1);
@@ -256,9 +253,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unu
out << "Checking zeros for jacMat took: " << t_zeros2.stop() << " s, found " << jacNumZeros << " zeros";
OpmLog::info(out.str());
}
status = backend->solve_system2(N, nnz, dim, static_cast<double*>(&(((*bridgeMat)[0][0][0][0]))), h_rows.data(), h_cols.data(), static_cast<double*>(&(b[0][0])),
jacNnz, static_cast<double*>(&(((*jacMat)[0][0][0][0]))), h_jacRows.data(), h_jacCols.data(),
wellContribs, result);
status = backend->solve_system2(matrix, static_cast<double*>(&(b[0][0])), jacMatrix, wellContribs, result);
}
switch(status) {

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

@@ -22,6 +22,7 @@
#include <opm/simulators/linalg/bda/BdaResult.hpp>
#include <opm/simulators/linalg/bda/BlockedMatrix.hpp>
#include <string>
@@ -85,13 +86,11 @@ namespace Accelerator {
virtual ~BdaSolver() {};
/// Define as pure virtual functions, so derivedclass must implement them
virtual SolverStatus solve_system(int N, int nnz, int dim,
double *vals, int *rows, int *cols,
virtual SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix,
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,
virtual SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) = 0;
virtual void get_result(double *x) = 0;

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

@@ -225,16 +225,16 @@ void FpgaSolverBackend<block_size>::get_result(double *x_)
template <unsigned int block_size>
SolverStatus FpgaSolverBackend<block_size>::solve_system(int N_, int nnz_, int dim, double *vals, int *rows, int *cols, double *b, WellContributions&, BdaResult &res)
SolverStatus FpgaSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions&, BdaResult &res)
{
if (initialized == false) {
initialize(N_, nnz_, dim, vals, rows, cols);
initialize(matrix->Nb, matrix->nnzbzs, matrix->nnzValues, matrix->rowPointers, matrix->colIndices);
if (!analyse_matrix()) {
return SolverStatus::BDA_SOLVER_ANALYSIS_FAILED;
}
}
perf_call.emplace_back();
update_system(vals, b);
update_system(matrix->nnzValues, b);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
@@ -251,31 +251,31 @@ 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)
SolverStatus FpgaSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
(void) nnz2; (void) vals2; (void) rows2; (void) cols2;
return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res);
return solve_system(matrix, b, wellContribs, res);
}
template <unsigned int block_size>
void FpgaSolverBackend<block_size>::initialize(int N_, int nnz_, int dim, double *vals, int *rows, int *cols)
void FpgaSolverBackend<block_size>::initialize(int Nb_, int nnzbs, double *vals, int *rows, int *cols)
{
double start = second();
this->N = N_;
this->nnz = nnz_;
this->nnzb = nnz_ / block_size / block_size;
Nb = (N + dim - 1) / dim;
this->Nb = Nb_;
this->N = Nb * block_size;
this->nnzb = nnzbs;
this->nnz = nnzbs * block_size * block_size;
// allocate host memory for matrices and vectors
// actual data for mat points to std::vector.data() in ISTLSolverEbos, so no alloc/free here
mat.reset(new BlockedMatrix(N_ / block_size, nnz_ / block_size / block_size, block_size, vals, cols, rows));
std::ostringstream oss;
oss << "Initializing FPGA data, matrix size: " << this->N << " blocks, nnz: " << this->nnzb << " blocks, " << \
"block size: " << dim << ", total nnz: " << this->nnz << "\n";
oss << "Initializing FPGA data, matrix size: " << this->Nb << " blockrows, nnz: " << this->nnzb << " blocks, " << \
"block size: " << block_size << ", total nnz: " << this->nnz << "\n";
oss << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance;
OpmLog::info(oss.str());

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

@@ -201,13 +201,12 @@ private:
unsigned short rst_settle_cycles = 0;
/// Allocate host memory
/// \param[in] N number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] Nb number of blockrows
/// \param[in] nnzbs number of blocks
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
void initialize(int N, int nnz, int dim, double *vals, int *rows, int *cols);
void initialize(int Nb, int nnzbs, int dim, double *vals, int *rows, int *cols);
/// Reorder the linear system so it corresponds with the coloring
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
@@ -243,21 +242,15 @@ public:
~FpgaSolverBackend();
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] N number of rows, divide by dim to get number of blockrows
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] wellContribs WellContributions, not used in FPGA solver because it requires them already added to matrix A
/// \param[inout] res summary of solver result
/// \return status code
SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override;
SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, 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,
SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) override;
/// Get result after linear solve, and peform postprocessing if necessary

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

@@ -52,14 +52,14 @@ amgclSolverBackend<block_size>::~amgclSolverBackend() {}
template <unsigned int block_size>
void amgclSolverBackend<block_size>::initialize(int N_, int nnz_, int dim) {
this->N = N_;
this->nnz = nnz_;
this->nnzb = nnz_ / block_size / block_size;
void amgclSolverBackend<block_size>::initialize(int Nb_, int nnzbs) {
this->Nb = Nb_;
this->N = Nb * block_size;
this->nnzb = nnzbs;
this->nnz = nnzbs * block_size * block_size;
Nb = (N + dim - 1) / dim;
std::ostringstream out;
out << "Initializing amgclSolverBackend, matrix size: " << N << " blocks, nnzb: " << nnzb << "\n";
out << "Initializing amgclSolverBackend, matrix size: " << Nb << " blockrows, nnzb: " << nnzb << " blocks\n";
out << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance << "\n";
out << "DeviceID: " << deviceID << "\n";
OpmLog::info(out.str());
@@ -360,22 +360,22 @@ 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)
SolverStatus amgclSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
(void) nnz2; (void) vals2; (void) rows2; (void) cols2;
return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res);
return solve_system(matrix, b, wellContribs, res);
}
template <unsigned int block_size>
SolverStatus amgclSolverBackend<block_size>::solve_system(int N_, int nnz_, int dim, double *vals, int *rows, int *cols, double *b, WellContributions&, BdaResult &res) {
SolverStatus amgclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions&, BdaResult &res) {
if (initialized == false) {
initialize(N_, nnz_, dim);
convert_sparsity_pattern(rows, cols);
initialize(matrix->Nb, matrix->nnzbs);
convert_sparsity_pattern(matrix->rowPointers, matrix->colIndices);
}
convert_data(vals, rows);
convert_data(matrix->nnzValues, matrix->rowPointers);
solve_system(b, res);
return SolverStatus::BDA_SOLVER_SUCCESS;
}

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

@@ -96,10 +96,9 @@ private:
std::once_flag vexcl_initialize;
#endif
/// Initialize host memory and determine amgcl parameters
/// \param[in] N number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
void initialize(int N, int nnz, int dim);
/// \param[in] Nb number of blockrows
/// \param[in] nnzbs number of blocks
void initialize(int Nb, int nnzbs);
/// Convert the BCSR sparsity pattern to a CSR one
/// \param[in] rows array of rowPointers, contains N/dim+1 values
@@ -129,25 +128,15 @@ public:
~amgclSolverBackend();
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] N number of rows, divide by dim to get number of blockrows
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz_prec number of nonzeroes of matrix for ILU0, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
/// \param[in] vals_prec array of nonzeroes for preconditioner
/// \param[in] rows_prec array of rowPointers for preconditioner
/// \param[in] cols_prec array of columnIndices for preconditioner
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] wellContribs WellContributions, to apply them separately, instead of adding them to matrix A
/// \param[inout] res summary of solver result
/// \return status code
SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override;
SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, 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,
SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) override;
/// Get result after linear solve, and peform postprocessing if necessary

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

@@ -188,17 +188,15 @@ void cusparseSolverBackend<block_size>::gpu_pbicgstab(WellContributions& wellCon
template <unsigned int block_size>
void cusparseSolverBackend<block_size>::initialize(int N, int nnz, int dim) {
this->N = N;
this->nnz = nnz;
this->nnzb = nnz / block_size / block_size;
Nb = (N + dim - 1) / dim;
void cusparseSolverBackend<block_size>::initialize(int Nb_, int nnzb) {
this->Nb = Nb_;
this->N = Nb * block_size;
this->nnzb = nnzb;
this->nnz = nnzb * block_size * block_size;
std::ostringstream out;
out << "Initializing GPU, matrix size: " << N << " blocks, nnz: " << nnzb << " blocks";
OpmLog::info(out.str());
out.str("");
out.clear();
out << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance;
out << "Initializing GPU, matrix size: " << Nb << " blockrows, nnz: " << nnzb << " blocks\n";
out << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance << "\n";
OpmLog::info(out.str());
cudaSetDevice(deviceID);
@@ -480,12 +478,12 @@ void cusparseSolverBackend<block_size>::get_result(double *x) {
template <unsigned int block_size>
SolverStatus cusparseSolverBackend<block_size>::solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) {
SolverStatus cusparseSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions& wellContribs, BdaResult &res) {
if (initialized == false) {
initialize(N, nnz, dim);
copy_system_to_gpu(vals, rows, cols, b);
initialize(matrix->Nb, matrix->nnzbs);
copy_system_to_gpu(matrix->nnzValues, matrix->rowPointers, matrix->colIndices, b);
} else {
update_system_on_gpu(vals, rows, b);
update_system_on_gpu(matrix->nnzValues, matrix->rowPointers, b);
}
if (analysis_done == false) {
if (!analyse_matrix()) {
@@ -500,14 +498,15 @@ SolverStatus cusparseSolverBackend<block_size>::solve_system(int N, int nnz, int
}
return SolverStatus::BDA_SOLVER_SUCCESS;
}
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)
SolverStatus cusparseSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
(void) nnz2; (void) vals2; (void) rows2; (void) cols2;
return solve_system(N_, nnz_, dim, vals, rows, cols, b, wellContribs, res);
return solve_system(matrix, b, wellContribs, res);
}
#define INSTANTIATE_BDA_FUNCTIONS(n) \

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

@@ -75,10 +75,9 @@ private:
void gpu_pbicgstab(WellContributions& wellContribs, BdaResult& res);
/// Initialize GPU and allocate memory
/// \param[in] N number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
void initialize(int N, int nnz, int dim);
/// \param[in] Nb number of blockrows
/// \param[in] nnzbs number of blocks
void initialize(int Nb, int nnzbs);
/// Clean memory
void finalize();
@@ -126,21 +125,15 @@ public:
~cusparseSolverBackend();
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] N number of rows, divide by dim to get number of blockrows
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] wellContribs contains all WellContributions, to apply them separately, instead of adding them to matrix A
/// \param[inout] res summary of solver result
/// \return status code
SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override;
SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, 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,
SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) override;
/// Get resulting vector x after linear solve, also includes post processing if necessary

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

@@ -383,14 +383,14 @@ void openclSolverBackend<block_size>::gpu_pbicgstab(WellContributions& wellContr
template <unsigned int block_size>
void openclSolverBackend<block_size>::initialize(int N_, int nnz_, int dim, double *vals, int *rows, int *cols) {
this->N = N_;
this->nnz = nnz_;
this->nnzb = nnz_ / block_size / block_size;
void openclSolverBackend<block_size>::initialize(std::shared_ptr<BlockedMatrix> matrix) {
this->Nb = matrix->Nb;
this->N = Nb * block_size;
this->nnzb = matrix->nnzbs;
this->nnz = nnzb * block_size * block_size;
Nb = (N + dim - 1) / dim;
std::ostringstream out;
out << "Initializing GPU, matrix size: " << N << " blocks, nnzb: " << nnzb << "\n";
out << "Initializing GPU, matrix size: " << Nb << " blockrows, nnzb: " << nnzb << "\n";
out << "Maxit: " << maxit << std::scientific << ", tolerance: " << tolerance << "\n";
out << "PlatformID: " << platformID << ", deviceID: " << deviceID << "\n";
OpmLog::info(out.str());
@@ -403,7 +403,8 @@ void openclSolverBackend<block_size>::initialize(int N_, int nnz_, int dim, doub
#if COPY_ROW_BY_ROW
vals_contiguous.resize(nnz);
#endif
mat.reset(new BlockedMatrix(Nb, nnzb, block_size, vals, cols, rows));
// mat.reset(new BlockedMatrix(Nb, nnzb, block_size, vals, cols, rows));
mat = matrix;
d_x = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_b = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
@@ -442,16 +443,15 @@ void openclSolverBackend<block_size>::initialize(int N_, int nnz_, int dim, doub
} // end initialize()
template <unsigned int block_size>
void openclSolverBackend<block_size>::initialize2(int N_, int nnz_, int dim, double *vals, int *rows, int *cols,
int nnz2, double *vals2, int *rows2, int *cols2) {
this->N = N_;
this->nnz = nnz_;
this->nnzb = nnz_ / block_size / block_size;
void openclSolverBackend<block_size>::initialize2(std::shared_ptr<BlockedMatrix> matrix, std::shared_ptr<BlockedMatrix> jacMatrix) {
this->Nb = matrix->Nb;
this->N = Nb * block_size;
this->nnzb = matrix->nnzbs;
this->nnz = nnzb * block_size * block_size;
this->jac_nnz = jacMatrix->nnzbs;
this->jac_nnzb = jac_nnz * block_size * block_size;
this->jac_nnz = nnz2;
this->jac_nnzb = nnz2 / block_size / block_size;
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";
@@ -467,8 +467,10 @@ void openclSolverBackend<block_size>::initialize2(int N_, int nnz_, int dim, dou
#if COPY_ROW_BY_ROW
vals_contiguous = new double[N];
#endif
mat.reset(new BlockedMatrix(Nb, nnzb, block_size, vals, cols, rows));
jacMat.reset(new BlockedMatrix(Nb, jac_nnzb, block_size, vals2, cols2, rows2));
// mat.reset(new BlockedMatrix(Nb, nnzb, block_size, vals, cols, rows));
// jacMat.reset(new BlockedMatrix(Nb, jac_nnzb, block_size, vals2, cols2, rows2));
mat = matrix;
jacMat = jacMatrix;
d_x = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_b = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
@@ -710,21 +712,21 @@ void openclSolverBackend<block_size>::get_result(double *x) {
template <unsigned int block_size>
SolverStatus openclSolverBackend<block_size>::solve_system(int N_, int nnz_, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) {
SolverStatus openclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions& wellContribs, BdaResult &res) {
if (initialized == false) {
initialize(N_, nnz_, dim, vals, rows, cols);
initialize(matrix);
if (analysis_done == false) {
if (!analyze_matrix()) {
return SolverStatus::BDA_SOLVER_ANALYSIS_FAILED;
}
}
update_system(vals, b, wellContribs);
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
copy_system_to_gpu();
} else {
update_system(vals, b, wellContribs);
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
@@ -734,25 +736,25 @@ SolverStatus openclSolverBackend<block_size>::solve_system(int N_, int nnz_, int
return SolverStatus::BDA_SOLVER_SUCCESS;
}
template <unsigned int block_size>
SolverStatus openclSolverBackend<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)
SolverStatus openclSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res)
{
if (initialized == false) {
blockJacVersion = true;
initialize2(N_, nnz_, dim, vals, rows, cols, nnz2, vals2, rows2, cols2);
initialize2(matrix, jacMatrix);
if (analysis_done == false) {
if (!analyze_matrix()) {
return SolverStatus::BDA_SOLVER_ANALYSIS_FAILED;
}
}
update_system(vals, b, wellContribs);
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
copy_system_to_gpu();
} else {
update_system(vals, b, wellContribs);
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}

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

@@ -72,8 +72,8 @@ private:
bool is_root; // allow for nested solvers, the root solver is called by BdaBridge
int *toOrder = nullptr, *fromOrder = nullptr; // BILU0 reorders rows of the matrix via these mappings
bool analysis_done = false;
std::unique_ptr<BlockedMatrix> mat = nullptr; // original matrix
std::unique_ptr<BlockedMatrix> jacMat = nullptr; // matrix for preconditioner
std::shared_ptr<BlockedMatrix> mat = nullptr; // original matrix
std::shared_ptr<BlockedMatrix> jacMat = nullptr; // matrix for preconditioner
BlockedMatrix *rmat = nullptr; // reordered matrix (or original if no reordering), used for spmv
ILUReorder opencl_ilu_reorder; // reordering strategy
std::vector<cl::Event> events;
@@ -135,16 +135,10 @@ private:
void gpu_pbicgstab(WellContributions& wellContribs, BdaResult& res);
/// Initialize GPU and allocate memory
/// \param[in] N number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
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);
/// \param[in] matrix matrix A
/// \param[in] jacMatrix matrix for preconditioner
void initialize(std::shared_ptr<BlockedMatrix> matrix);
void initialize2(std::shared_ptr<BlockedMatrix> matrix, std::shared_ptr<BlockedMatrix> jacMatrix);
/// Clean memory
void finalize();
@@ -197,25 +191,16 @@ public:
~openclSolverBackend();
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] N number of rows, divide by dim to get number of blockrows
/// \param[in] nnz number of nonzeroes, divide by dim*dim to get number of blocks
/// \param[in] nnz_prec number of nonzeroes of matrix for ILU0, divide by dim*dim to get number of blocks
/// \param[in] dim size of block
/// \param[in] vals array of nonzeroes, each block is stored row-wise and contiguous, contains nnz values
/// \param[in] rows array of rowPointers, contains N/dim+1 values
/// \param[in] cols array of columnIndices, contains nnz values
/// \param[in] vals_prec array of nonzeroes for preconditioner
/// \param[in] rows_prec array of rowPointers for preconditioner
/// \param[in] cols_prec array of columnIndices for preconditioner
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] wellContribs WellContributions, to apply them separately, instead of adding them to matrix A
/// \param[inout] res summary of solver result
/// \return status code
SolverStatus solve_system(int N, int nnz, int dim, double *vals, int *rows, int *cols, double *b, WellContributions& wellContribs, BdaResult &res) override;
SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, 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;
SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) override;
/// Solve scalar linear system, for example a coarse system of an AMG preconditioner
/// Data is already on the GPU