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

Merge some duplicate functions

Browse Source
This commit is contained in:
Tong Dong Qiu 2022-02-23 16:28:37 +01:00
parent d963820e48
commit 3797b7297d
10 changed files with 60 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -192,14 +192,13 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::copySparsityPatternFromI
template <class BridgeMatrix, class BridgeVector, int block_size>
void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unused]] BridgeMatrix* bridgeMat,
[[maybe_unused]] BridgeMatrix* jacMat,
[[maybe_unused]] int numJacobiBlocks,
[[maybe_unused]] BridgeVector& b,
[[maybe_unused]] WellContributions& wellContribs,
[[maybe_unused]] InverseOperatorResult& res)
void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system(BridgeMatrix* bridgeMat,
BridgeMatrix* jacMat,
int numJacobiBlocks,
BridgeVector& b,
WellContributions& wellContribs,
InverseOperatorResult& res)
{
if (use_gpu || use_fpga) {
BdaResult result;
result.converged = false;
@ -229,14 +228,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system([[maybe_unu
OpmLog::info(out.str());
}
/////////////////////////
// actually solve
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(matrix, static_cast<double*>(&(b[0][0])), wellContribs, result);
} else {
if (numJacobiBlocks >= 2) {
const int jacNnzb = (h_jacRows.empty()) ? jacMat->nonzeroes() : h_jacRows.back();
if (!jacMatrix) {
@ -253,9 +245,13 @@ 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(matrix, static_cast<double*>(&(b[0][0])), jacMatrix, wellContribs, result);
}
/////////////////////////
// actually solve
// assume that underlying data (nonzeroes) from b (Dune::BlockVector) are contiguous, if this is not the case, the chosen BdaSolver is expected to perform undefined behaviour
SolverStatus status = backend->solve_system(matrix, static_cast<double*>(&(b[0][0])), jacMatrix, wellContribs, result);
switch(status) {
case SolverStatus::BDA_SOLVER_SUCCESS:
//OpmLog::info("BdaSolver converged");

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

@ -86,12 +86,8 @@ namespace Accelerator {
virtual ~BdaSolver() {};
/// Define as pure virtual functions, so derivedclass must implement them
virtual SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix,
double *b, WellContributions& wellContribs, BdaResult &res) = 0;
virtual SolverStatus solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res) = 0;
virtual SolverStatus solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix, WellContributions& wellContribs, BdaResult &res) = 0;
virtual void get_result(double *x) = 0;

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

@ -225,7 +225,11 @@ void FpgaSolverBackend<block_size>::get_result(double *x_)
template <unsigned int block_size>
SolverStatus FpgaSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions&, BdaResult &res)
SolverStatus FpgaSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions&,
BdaResult &res)
{
if (initialized == false) {
initialize(matrix->Nb, matrix->nnzbzs, matrix->nnzValues, matrix->rowPointers, matrix->colIndices);
@ -250,15 +254,6 @@ SolverStatus FpgaSolverBackend<block_size>::solve_system(std::shared_ptr<Blocked
return SolverStatus::BDA_SOLVER_SUCCESS;
}
template <unsigned int block_size>
SolverStatus FpgaSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
return solve_system(matrix, b, wellContribs, res);
}
template <unsigned int block_size>
void FpgaSolverBackend<block_size>::initialize(int Nb_, int nnzbs, double *vals, int *rows, int *cols)

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

@ -244,14 +244,12 @@ public:
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] jacMatrix matrix for preconditioner
/// \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(std::shared_ptr<BlockedMatrix> matrix, double *b, 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;
SolverStatus solve_system(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
/// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array

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

@ -360,17 +360,12 @@ void amgclSolverBackend<block_size>::get_result(double *x_) {
template <unsigned int block_size>
SolverStatus amgclSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
SolverStatus amgclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
[[maybe_unused]] WellContributions& wellContribs,
BdaResult &res)
{
return solve_system(matrix, b, wellContribs, res);
}
template <unsigned int block_size>
SolverStatus amgclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions&, BdaResult &res) {
if (initialized == false) {
initialize(matrix->Nb, matrix->nnzbs);
convert_sparsity_pattern(matrix->rowPointers, matrix->colIndices);

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

@ -130,14 +130,12 @@ public:
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] jacMatrix matrix for preconditioner
/// \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(std::shared_ptr<BlockedMatrix> matrix, double *b, 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;
SolverStatus solve_system(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
/// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array

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

@ -478,7 +478,12 @@ void cusparseSolverBackend<block_size>::get_result(double *x) {
template <unsigned int block_size>
SolverStatus cusparseSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions& wellContribs, BdaResult &res) {
SolverStatus cusparseSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
if (initialized == false) {
initialize(matrix->Nb, matrix->nnzbs);
copy_system_to_gpu(matrix->nnzValues, matrix->rowPointers, matrix->colIndices, b);
@ -499,15 +504,6 @@ SolverStatus cusparseSolverBackend<block_size>::solve_system(std::shared_ptr<Blo
return SolverStatus::BDA_SOLVER_SUCCESS;
}
template <unsigned int block_size>
SolverStatus cusparseSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix,
double *b,
[[maybe_unused]] std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
return solve_system(matrix, b, wellContribs, res);
}
#define INSTANTIATE_BDA_FUNCTIONS(n) \
template cusparseSolverBackend<n>::cusparseSolverBackend(int, int, double, unsigned int); \

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

@ -127,14 +127,12 @@ public:
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] jacMatrix matrix for preconditioner
/// \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(std::shared_ptr<BlockedMatrix> matrix, double *b, 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;
SolverStatus solve_system(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
/// \param[inout] x resulting x vector, caller must guarantee that x points to a valid array

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

@ -383,78 +383,26 @@ void openclSolverBackend<block_size>::gpu_pbicgstab(WellContributions& wellContr
template <unsigned int block_size>
void openclSolverBackend<block_size>::initialize(std::shared_ptr<BlockedMatrix> matrix) {
void openclSolverBackend<block_size>::initialize(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;
std::ostringstream out;
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());
out.str("");
out.clear();
try {
prec->setOpencl(context, queue);
#if COPY_ROW_BY_ROW
vals_contiguous.resize(nnz);
#endif
// 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);
d_rb = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_r = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_rw = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_p = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_pw = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_s = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_t = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_v = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_tmp = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * N);
d_Avals = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(double) * nnz);
d_Acols = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(int) * nnzb);
d_Arows = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(int) * (Nb + 1));
bool reorder = (opencl_ilu_reorder != ILUReorder::NONE);
if (reorder) {
rb = new double[N];
d_toOrder = cl::Buffer(*context, CL_MEM_READ_WRITE, sizeof(int) * Nb);
}
} catch (const cl::Error& error) {
std::ostringstream oss;
oss << "OpenCL Error: " << error.what() << "(" << error.err() << ")\n";
oss << getErrorString(error.err());
// rethrow exception
OPM_THROW(std::logic_error, oss.str());
} catch (const std::logic_error& error) {
// rethrow exception by OPM_THROW in the try{}, without this, a segfault occurs
throw error;
if (jacMatrix) {
useJacMatrix = true;
}
initialized = true;
} // end initialize()
template <unsigned int 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;
if (useJacMatrix) {
this->jac_nnz = jacMatrix->nnzbs;
this->jac_nnzb = jac_nnz * block_size * block_size;
}
std::ostringstream out;
out << "Initializing GPU, matrix size: " << N << " blocks, nnzb: " << nnzb << "\n";
out << "Blocks in ILU matrix: " << jac_nnzb << "\n";
if (useJacMatrix) {
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());
@ -594,7 +542,7 @@ bool openclSolverBackend<block_size>::analyze_matrix() {
Timer t;
bool success;
if (blockJacVersion)
if (useJacMatrix)
success = prec->analyze_matrix(mat.get(), jacMat.get());
else
success = prec->analyze_matrix(mat.get());
@ -649,7 +597,7 @@ bool openclSolverBackend<block_size>::create_preconditioner() {
Timer t;
bool result;
if (blockJacVersion)
if (useJacMatrix)
result = prec->create_preconditioner(mat.get(), jacMat.get());
else
result = prec->create_preconditioner(mat.get());
@ -712,37 +660,14 @@ void openclSolverBackend<block_size>::get_result(double *x) {
template <unsigned int block_size>
SolverStatus openclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix, double *b, WellContributions& wellContribs, BdaResult &res) {
if (initialized == false) {
initialize(matrix);
if (analysis_done == false) {
if (!analyze_matrix()) {
return SolverStatus::BDA_SOLVER_ANALYSIS_FAILED;
}
}
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
copy_system_to_gpu();
} else {
update_system(matrix->nnzValues, b, wellContribs);
if (!create_preconditioner()) {
return SolverStatus::BDA_SOLVER_CREATE_PRECONDITIONER_FAILED;
}
update_system_on_gpu();
}
solve_system(wellContribs, res);
return SolverStatus::BDA_SOLVER_SUCCESS;
}
template <unsigned int block_size>
SolverStatus openclSolverBackend<block_size>::solve_system2(std::shared_ptr<BlockedMatrix> matrix, double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs, BdaResult &res)
SolverStatus openclSolverBackend<block_size>::solve_system(std::shared_ptr<BlockedMatrix> matrix,
double *b,
std::shared_ptr<BlockedMatrix> jacMatrix,
WellContributions& wellContribs,
BdaResult &res)
{
if (initialized == false) {
blockJacVersion = true;
initialize2(matrix, jacMatrix);
initialize(matrix, jacMatrix);
if (analysis_done == false) {
if (!analyze_matrix()) {
return SolverStatus::BDA_SOLVER_ANALYSIS_FAILED;
@ -764,6 +689,7 @@ SolverStatus openclSolverBackend<block_size>::solve_system2(std::shared_ptr<Bloc
return SolverStatus::BDA_SOLVER_SUCCESS;
}
#define INSTANTIATE_BDA_FUNCTIONS(n) \
template openclSolverBackend<n>::openclSolverBackend( \
int, int, double, unsigned int, unsigned int, ILUReorder, std::string); \

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

@ -65,7 +65,7 @@ private:
int jac_nnz;
int jac_nnzb;
bool blockJacVersion = false;
bool useJacMatrix = false;
std::unique_ptr<Preconditioner<block_size> > prec;
// can perform blocked ILU0 and AMG on pressure component
@ -137,8 +137,7 @@ private:
/// Initialize GPU and allocate memory
/// \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);
void initialize(std::shared_ptr<BlockedMatrix> matrix, std::shared_ptr<BlockedMatrix> jacMatrix);
/// Clean memory
void finalize();
@ -193,14 +192,12 @@ public:
/// Solve linear system, A*x = b, matrix A must be in blocked-CSR format
/// \param[in] matrix matrix A
/// \param[in] b input vector, contains N values
/// \param[in] jacMatrix matrix for preconditioner
/// \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(std::shared_ptr<BlockedMatrix> matrix, double *b, 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;
SolverStatus solve_system(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