OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-11-02 09:16:07 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6241 from kjetilly/adding_gpu_preconditioners

Browse Source 
Instatiate preconditioners directly with GPU matrices
This commit is contained in:
Atgeirr Flø Rasmussen authored and GitHub committed 2025-05-13 11:26:05 +02:00
2 parents 97deff78bb dd7da86e45
commit 46d34a0b52
27 files changed
+745 -241
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

tests/gpuistl/test_GpuDILU.cpp
+6 -3
View File
@@ -30,6 +30,7 @@
#include <opm/simulators/linalg/gpuistl/GpuVector.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpu_safe_call.hpp>
#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
#include <opm/simulators/linalg/gpuistl/PreconditionerCPUMatrixToGPUMatrix.hpp>
#include <random>
#include <vector>
@@ -45,7 +46,9 @@ using CuMatrix = Opm::gpuistl::GpuSparseMatrix<T>;
using CuIntVec = Opm::gpuistl::GpuVector<int>;
using CuFloatingPointVec = Opm::gpuistl::GpuVector<T>;
using GpuDilu1x1 = Opm::gpuistl::GpuDILU<Sp1x1BlockMatrix, CuFloatingPointVec, CuFloatingPointVec>;
using Wrapped1x1 = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<CuFloatingPointVec, CuFloatingPointVec, GpuDilu1x1, Sp1x1BlockMatrix>;
using GpuDilu2x2 = Opm::gpuistl::GpuDILU<Sp2x2BlockMatrix, CuFloatingPointVec, CuFloatingPointVec>;
using Wrapped2x2 = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<CuFloatingPointVec, CuFloatingPointVec, GpuDilu2x2, Sp2x2BlockMatrix>;
Sp1x1BlockMatrix
get1x1BlockTestMatrix()
@@ -211,7 +214,7 @@ BOOST_AUTO_TEST_CASE(TestDiluApply)
// Initialize preconditioner objects
Dune::MultithreadDILU<Sp1x1BlockMatrix, B1x1Vec, B1x1Vec> cpudilu(matA);
auto gpudilu = GpuDilu1x1(matA, true, true, false);
auto gpudilu = Wrapped1x1(matA, matA, true, true, false);
// Use the apply
gpudilu.apply(d_output, d_input);
@@ -235,7 +238,7 @@ BOOST_AUTO_TEST_CASE(TestDiluApplyBlocked)
// init matrix with 2x2 blocks
Sp2x2BlockMatrix matA = get2x2BlockTestMatrix();
auto gpudilu = GpuDilu2x2(matA, true, true, false);
auto gpudilu = Wrapped2x2(matA, matA, true, true, false);
Dune::MultithreadDILU<Sp2x2BlockMatrix, B2x2Vec, B2x2Vec> cpudilu(matA);
// create input/output buffers for the apply
@@ -275,7 +278,7 @@ BOOST_AUTO_TEST_CASE(TestDiluInitAndUpdateLarge)
{
// create gpu dilu preconditioner
Sp1x1BlockMatrix matA = get1x1BlockTestMatrix();
auto gpudilu = GpuDilu1x1(matA, true, true, false);
auto gpudilu = Wrapped1x1(matA, matA, true, true, false);
matA[0][0][0][0] = 11.0;
matA[0][1][0][0] = 12.0;
tests/gpuistl/test_GpuJac.cpp
+11 -4
View File
@@ -30,6 +30,7 @@
#include <opm/simulators/linalg/gpuistl/detail/gpusparse_matrix_operations.hpp>
#include <opm/simulators/linalg/gpuistl/detail/fix_zero_diagonal.hpp>
#include <opm/simulators/linalg/gpuistl/detail/vector_operations.hpp>
#include <opm/simulators/linalg/gpuistl/PreconditionerCPUMatrixToGPUMatrix.hpp>
using NumericTypes = boost::mpl::list<double, float>;
@@ -49,7 +50,10 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize2, T, NumericTypes)
using M = Dune::FieldMatrix<T, blocksize, blocksize>;
using SpMatrix = Dune::BCRSMatrix<M>;
using Vector = Dune::BlockVector<Dune::FieldVector<T, blocksize>>;
using GpuJac = Opm::gpuistl::GpuJac<SpMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
using GpuMatrix = Opm::gpuistl::GpuSparseMatrix<T>;
using GpuJac = Opm::gpuistl::GpuJac<GpuMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
using GpuVector = Opm::gpuistl::GpuVector<T>;
using Wrapped = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<GpuVector, GpuVector, GpuJac, SpMatrix>;
SpMatrix B(N, N, nonZeroes, SpMatrix::row_wise);
for (auto row = B.createbegin(); row != B.createend(); ++row) {
@@ -70,7 +74,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize2, T, NumericTypes)
B[1][1][0][0] = -1.0;
B[1][1][1][1] = -1.0;
auto gpujac = Opm::gpuistl::PreconditionerAdapter<Vector, Vector, GpuJac>(std::make_shared<GpuJac>(B, 0.5));
auto gpujac = Opm::gpuistl::PreconditionerAdapter<Vector, Vector, Wrapped>(std::make_shared<Wrapped>(B, 0.5));
Vector vVector(2);
Vector dVector(2);
@@ -103,7 +107,10 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize1, T, NumericTypes)
using M = Dune::FieldMatrix<T, blocksize, blocksize>;
using SpMatrix = Dune::BCRSMatrix<M>;
using Vector = Dune::BlockVector<Dune::FieldVector<T, blocksize>>;
using GpuJac = Opm::gpuistl::GpuJac<SpMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
using GpuMatrix = Opm::gpuistl::GpuSparseMatrix<T>;
using GpuJac = Opm::gpuistl::GpuJac<GpuMatrix, Opm::gpuistl::GpuVector<T>, Opm::gpuistl::GpuVector<T>>;
using GpuVector = Opm::gpuistl::GpuVector<T>;
using Wrapped = Opm::gpuistl::PreconditionerCPUMatrixToGPUMatrix<GpuVector, GpuVector, GpuJac, SpMatrix>;
SpMatrix B(N, N, nonZeroes, SpMatrix::row_wise);
for (auto row = B.createbegin(); row != B.createend(); ++row) {
@@ -129,7 +136,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(GPUJACApplyBlocksize1, T, NumericTypes)
B[2][2][0][0] = -1.0;
B[3][3][0][0] = -1.0;
auto gpujac = Opm::gpuistl::PreconditionerAdapter<Vector, Vector, GpuJac>(std::make_shared<GpuJac>(B, 0.5));
auto gpujac = Opm::gpuistl::PreconditionerAdapter<Vector, Vector, Wrapped>(std::make_shared<Wrapped>(B, 0.5));
Vector vVector(4);
Vector dVector(4);
tests/gpuistl/test_preconditioner_factory_gpu.cpp
+82 -39
View File
@@ -13,6 +13,7 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "opm/simulators/linalg/PropertyTree.hpp"
#include <boost/test/tools/old/interface.hpp>
#include <config.h>
#include <cstddef>
@@ -27,17 +28,18 @@
#include <opm/simulators/linalg/PreconditionerFactory.hpp>
#include <opm/simulators/linalg/PreconditionerFactory_impl.hpp>
using Scalar = double;
constexpr static int Dim = 3;
using MatrixTypeCPU = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, Dim, Dim>>;
// using VectorType = Dune::BlockVector<Dune::FieldVector<Scalar, Dim>>;
using CommSeq = Dune::Amg::SequentialInformation;
// using GpuOperator = Dune::MatrixAdapter<MatrixType,
// VectorType,
// VectorType>;
using ScalarT = double;
constexpr static int Dim = 1;
using GpuMatrixType = Opm::gpuistl::GpuSparseMatrix<Scalar>;
using GpuVectorType = Opm::gpuistl::GpuVector<Scalar>;
using CommSeq = Dune::Amg::SequentialInformation;
using MatrixTypeCPU = Dune::BCRSMatrix<Dune::FieldMatrix<ScalarT, Dim, Dim>>;
using VectorCPU = Dune::BlockVector<Dune::FieldVector<ScalarT, 1>>;
using CpuOperatorType = Dune::MatrixAdapter<MatrixTypeCPU, VectorCPU, VectorCPU>;
using FactoryTypeCpu = Opm::PreconditionerFactory<CpuOperatorType, CommSeq>;
using GpuMatrixType = Opm::gpuistl::GpuSparseMatrix<ScalarT>;
using GpuVectorType = Opm::gpuistl::GpuVector<ScalarT>;
using GpuOperatorType = Dune::MatrixAdapter<GpuMatrixType, GpuVectorType, GpuVectorType>;
using FactoryTypeGpu = Opm::PreconditionerFactory<GpuOperatorType, CommSeq>;
@@ -45,36 +47,77 @@ BOOST_AUTO_TEST_CASE(TestMatrixAdapter)
{
const std::size_t N = 10;
const std::size_t nonZeroes = N * 3 - 2;
MatrixTypeCPU matrix(N, N, nonZeroes, MatrixTypeCPU::row_wise);
for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) {
// Add nonzeros for left neighbour, diagonal and right neighbour
if (row.index() > 0) {
row.insert(row.index() - 1);
// The name of the preconditioner on the GPU directly (key) and
// the name of the wrapped preconditioner (value).
auto preconditionerPairs = std::map<std::string, std::string> {
{"ILU0", "GPUILU0"},
{"JAC", "GPUJAC"},
{"OPMILU0", "OPMGPUILU0"},
{"DILU", "GPUDILU"},
};
for (const auto& [gpuPreconditionerName, wrappedPreconditionerName] : preconditionerPairs) {
MatrixTypeCPU matrix(N, N, nonZeroes, MatrixTypeCPU::row_wise);
for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) {
// Add nonzeros for left neighbour, diagonal and right neighbour
if (row.index() > 0) {
row.insert(row.index() - 1);
}
row.insert(row.index());
if (row.index() < matrix.N() - 1) {
row.insert(row.index() + 1);
}
}
row.insert(row.index());
if (row.index() < matrix.N() - 1) {
row.insert(row.index() + 1);
for (std::size_t i = 0; i < N; ++i) {
matrix[i][i] = -2;
if (i < int(N) - 1) {
matrix[i][i + 1] = 1;
}
if (i > 0) {
matrix[i][i - 1] = 1;
}
}
GpuMatrixType gpuMatrix = GpuMatrixType::fromMatrix(matrix);
GpuOperatorType gpuOp(gpuMatrix);
CpuOperatorType cpuOp(matrix);
Opm::PropertyTree prm;
prm.put(std::string("type"), gpuPreconditionerName);
Opm::PropertyTree prmCPU;
prmCPU.put(std::string("type"), wrappedPreconditionerName);
auto preconditionerGPU = FactoryTypeGpu::create(gpuOp, prm);
auto preconditionerWrapped = FactoryTypeCpu::create(cpuOp, prmCPU);
// check for the standard basis {e_i}
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
for (std::size_t i = 0; i < N; ++i) {
VectorCPU inputVector(N);
inputVector[i][0] = 1.0;
VectorCPU outputVectorWrapped(N);
auto outputVectorGPU = Opm::gpuistl::GpuVector<ScalarT>(N);
auto inputVectorGPU = Opm::gpuistl::GpuVector<ScalarT>(inputVector);
preconditionerGPU->apply(outputVectorGPU, inputVectorGPU);
preconditionerWrapped->apply(outputVectorWrapped, inputVector);
auto outputVectorgGPUonCPU = outputVectorGPU.asDuneBlockVector<Dim>();
for (std::size_t component = 0; component < N; ++component) {
BOOST_TEST(outputVectorWrapped[component][0] == outputVectorgGPUonCPU[component][0],
"Component " << component << " of the output vector is not equal to the expected value."
<< "\nExpected: " << outputVectorWrapped[component][0]
<< ",\nbut got: " << outputVectorgGPUonCPU[component][0]
<< "\nfor preconditioner: " << gpuPreconditionerName
<< " and wrapped preconditioner: " << wrappedPreconditionerName);
}
}
}
for (std::size_t i = 0; i < N; ++i) {
matrix[i][i] = -2;
if (i < int(N) - 1) {
matrix[i][i + 1] = 1;
}
if (i > 0) {
matrix[i][i - 1] = 1;
}
}
GpuMatrixType gpuMatrix = GpuMatrixType::fromMatrix(matrix);
Dune::MatrixAdapter<GpuMatrixType, GpuVectorType, GpuVectorType> gpuOp(gpuMatrix);
Opm::PropertyTree prm;
prm.put(std::string("preconditioner"), std::string("DILU"));
auto weightsCalculator = []() { return GpuVectorType(N); };
// TODO: Once the GPU preconditioners are implemented, we can remove the check throw
// and actually test the preconditioner creation.
BOOST_CHECK_THROW(FactoryTypeGpu::create(gpuOp, prm, weightsCalculator), std::invalid_argument);
}