mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Added conversion preconditioner.
This commit is contained in:
Kjetil Olsen Lye
parent
b30e6d79d5
commit
84305a7a8d
@ -172,6 +172,7 @@ if(CUDA_FOUND)
|
||||
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/has_function.hpp)
|
||||
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/preconditioner_should_call_post_pre.hpp)
|
||||
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/PreconditionerAdapter.hpp)
|
||||
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/PreconditionerConvertFieldTypeAdapter.hpp)
|
||||
|
||||
|
||||
endif()
|
||||
@ -261,6 +262,7 @@ if(CUDA_FOUND)
|
||||
list(APPEND TEST_SOURCE_FILES tests/cuistl/test_cusparsematrix.cpp)
|
||||
list(APPEND TEST_SOURCE_FILES tests/cuistl/test_safe_conversion.cpp)
|
||||
list(APPEND TEST_SOURCE_FILES tests/cuistl/test_cuseqilu0.cpp)
|
||||
list(APPEND TEST_SOURCE_FILES tests/cuistl/test_converttofloatadapter.cpp)
|
||||
|
||||
endif()
|
||||
if(OPENCL_FOUND)
|
||||
|
||||
@ -0,0 +1,240 @@
|
||||
/*
|
||||
Copyright 2022-2023 SINTEF AS
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef OPM_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
|
||||
#define OPM_PRECONDITIONERCONVERTOFLOATADAPTER_HPP
|
||||
#include <cusparse.h>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/preconditioner.hh>
|
||||
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/CuSparseMatrix.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/CuMatrixDescription.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/CuSparseHandle.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/CuSparseResource.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/cusparse_constants.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/cusparse_safe_call.hpp>
|
||||
#include <opm/simulators/linalg/cuistl/detail/preconditioner_should_call_post_pre.hpp>
|
||||
|
||||
|
||||
namespace Opm::cuistl
|
||||
{
|
||||
//! \brief Converts the field type (eg. double to float) to benchmark single precision preconditioners
|
||||
//!
|
||||
//! \note This is not a fast conversion, it is simply meant to benchmark the potential of some
|
||||
//! preconditioners on consumer grade GPUs where the double precision performance is often
|
||||
//! artificially limited.
|
||||
//!
|
||||
//! \note In theory this can do any field_type conversion that is meaningful, but it is only tested
|
||||
//! on double to float conversion.
|
||||
//!
|
||||
//! \note Remember to set the underlying preconditioner with setUnderlyingPreconditioner (should use the matrix from
|
||||
//! getConvertedMatrix())
|
||||
//!
|
||||
//! \note One could simply change the constructor design by accepting a creator function for the underlying
|
||||
//! preconditioner. For the current use cases this is however not needed.
|
||||
//!
|
||||
//! To use this, use something like the following code:
|
||||
//! \code{.cpp}
|
||||
//! #include <opm/simulators/linalg/cuistl/PreconditionerConvertFieldTypeAdapter.hpp>
|
||||
//! #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
|
||||
//!
|
||||
//! using XDouble = Dune::BlockVector<Dune::FieldVector<double, 2>>;
|
||||
//! using MDouble = Dune::FieldMatrix<double, 2, 2>;
|
||||
//! using SpMatrixDouble = Dune::BCRSMatrix<MDouble>;
|
||||
//! using XFloat = Dune::BlockVector<Dune::FieldVector<float, 2>>;
|
||||
//! using MFloat = Dune::FieldMatrix<float, 2, 2>;
|
||||
//! using SpMatrixFloat = Dune::BCRSMatrix<MFloat>;
|
||||
//!
|
||||
//! template<class ParallelInfo>
|
||||
//! void applyILU0AsFloat(const MDouble& matrix, const XDouble& x, XDouble& y) {
|
||||
//!
|
||||
//! using FloatILU0 = typename Opm::ParallelOverlappingILU0<MFloat, XFloat, XFloat, ParallelInfo>;
|
||||
//! using DoubleToFloatConverter = typename Opm::cuistl::PreconditionerConvertFieldTypeAdapter<FloatILU0, MDouble,
|
||||
//! XDouble, XDouble>;
|
||||
//!
|
||||
//! // Note that we do not need to make a new instance for every invocation, this
|
||||
//! // is just done for this example
|
||||
//! auto doubleToFloatConverter = DoubleToFloatConverter(matrix);
|
||||
//! const auto& convertedMatrix = doubleToFloatConverter.getConvertedMatrix();
|
||||
//!
|
||||
//! auto floatILU0 = std::make_shared<FloatILU0>(convertedMatrix, 0, 1.0, 0);
|
||||
//!
|
||||
//! doubleToFloatConverter.setUnderlyingPreconditioner(floatILU0);
|
||||
//!
|
||||
//! // This will convert x and y to float, thenn call floatILU0.apply on the converted arguments
|
||||
//! doubleToFloatConverter.apply(x, y);
|
||||
//! }
|
||||
//!
|
||||
//! \endcode
|
||||
template <class CudaPreconditionerType, class M, class X, class Y, int l = 1>
|
||||
class PreconditionerConvertFieldTypeAdapter : public Dune::PreconditionerWithUpdate<X, Y>
|
||||
{
|
||||
public:
|
||||
//! \brief The matrix type the preconditioner is for.
|
||||
using matrix_type = typename std::remove_const<M>::type;
|
||||
|
||||
//! \brief The domain type of the preconditioner.
|
||||
using domain_type = X;
|
||||
//! \brief The range type of the preconditioner.
|
||||
using range_type = Y;
|
||||
//! \brief The field type of the preconditioner.
|
||||
using field_type = typename X::field_type;
|
||||
|
||||
|
||||
using domain_type_to = typename CudaPreconditionerType::domain_type;
|
||||
//! \brief The range type of the preconditioner.
|
||||
using range_type_to = typename CudaPreconditionerType::range_type;
|
||||
//! \brief The field type of the preconditioner.
|
||||
using field_type_to = typename domain_type_to::field_type;
|
||||
|
||||
|
||||
using block_type = typename domain_type::block_type;
|
||||
|
||||
using XTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
|
||||
using YTo = Dune::BlockVector<Dune::FieldVector<field_type_to, block_type::dimension>>;
|
||||
using matrix_type_to =
|
||||
typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type_to, block_type::dimension, block_type::dimension>>;
|
||||
|
||||
//! \brief Constructor.
|
||||
//!
|
||||
//! \param A The matrix to operate on.
|
||||
//!
|
||||
//! \note After the PreconditionerConvertFieldTypeAdapter you can get the converted matrix
|
||||
//! by calling getConvertedMatrix(), which in turn can be used to create the underlying preconditioner.
|
||||
//! Once the underlying precondtioner has been called, this must be supplied to setUnderlyingPreconditioner.
|
||||
PreconditionerConvertFieldTypeAdapter(const M& matrix)
|
||||
: m_matrix(matrix)
|
||||
, m_convertedMatrix(createConvertedMatrix())
|
||||
{
|
||||
}
|
||||
|
||||
//! \brief Not used at the moment
|
||||
virtual void pre([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
|
||||
{
|
||||
static_assert(!detail::shouldCallPreconditionerPre<CudaPreconditionerType>(),
|
||||
"We currently do not support Preconditioner::pre().");
|
||||
}
|
||||
|
||||
//! \brief Apply the preconditoner.
|
||||
virtual void apply(X& v, const Y& d) override
|
||||
{
|
||||
OPM_ERROR_IF(!m_underlyingPreconditioner,
|
||||
"You need to set the underlying preconditioner with setUnderlyingPreconditioner.");
|
||||
XTo convertedV(v.N());
|
||||
for (size_t i = 0; i < v.N(); ++i) {
|
||||
for (size_t j = 0; j < block_type::dimension; ++j) {
|
||||
// This is probably unnecessary, but doing it anyway:
|
||||
convertedV[i][j] = field_type_to(v[i][j]);
|
||||
}
|
||||
}
|
||||
YTo convertedD(d.N());
|
||||
for (size_t i = 0; i < d.N(); ++i) {
|
||||
for (size_t j = 0; j < block_type::dimension; ++j) {
|
||||
convertedD[i][j] = field_type_to(d[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
m_underlyingPreconditioner->apply(convertedV, convertedD);
|
||||
|
||||
for (size_t i = 0; i < v.N(); ++i) {
|
||||
for (size_t j = 0; j < block_type::dimension; ++j) {
|
||||
v[i][j] = field_type(convertedV[i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//! \brief Not used at the moment
|
||||
virtual void post([[maybe_unused]] X& x) override
|
||||
{
|
||||
static_assert(!detail::shouldCallPreconditionerPost<CudaPreconditionerType>(),
|
||||
"We currently do not support Preconditioner::post().");
|
||||
}
|
||||
|
||||
//! Category of the preconditioner (see SolverCategory::Category)
|
||||
virtual Dune::SolverCategory::Category category() const override
|
||||
{
|
||||
return m_underlyingPreconditioner->category();
|
||||
}
|
||||
|
||||
virtual void update() override
|
||||
{
|
||||
OPM_ERROR_IF(!m_underlyingPreconditioner,
|
||||
"You need to set the underlying preconditioner with setUnderlyingPreconditioner.");
|
||||
updateMatrix();
|
||||
m_underlyingPreconditioner->update();
|
||||
}
|
||||
|
||||
const matrix_type_to& getConvertedMatrix() const
|
||||
{
|
||||
return m_convertedMatrix;
|
||||
}
|
||||
|
||||
void setUnderlyingPreconditioner(const std::shared_ptr<CudaPreconditionerType>& conditioner)
|
||||
{
|
||||
m_underlyingPreconditioner = conditioner;
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
void updateMatrix()
|
||||
{
|
||||
const auto nnz = m_matrix.nonzeroes() * m_matrix[0][0].N() * m_matrix[0][0].N();
|
||||
const auto dataPointerIn = static_cast<const field_type*>(&((m_matrix[0][0][0][0])));
|
||||
auto dataPointerOut = static_cast<field_type_to*>(&((m_convertedMatrix[0][0][0][0])));
|
||||
|
||||
std::vector<field_type_to> buffer(nnz, 0);
|
||||
for (size_t i = 0; i < nnz; ++i) {
|
||||
dataPointerOut[i] = field_type_to(dataPointerIn[i]);
|
||||
}
|
||||
}
|
||||
matrix_type_to createConvertedMatrix()
|
||||
{
|
||||
// TODO: Check if this whole conversion can be done more efficiently.
|
||||
const auto N = m_matrix.N();
|
||||
matrix_type_to matrixBuilder(N, N, m_matrix.nonzeroes(), matrix_type_to::row_wise);
|
||||
{
|
||||
auto rowIn = m_matrix.begin();
|
||||
for (auto rowOut = matrixBuilder.createbegin(); rowOut != matrixBuilder.createend(); ++rowOut) {
|
||||
for (auto column = rowIn->begin(); column != rowIn->end(); ++column) {
|
||||
rowOut.insert(column.index());
|
||||
}
|
||||
++rowIn;
|
||||
}
|
||||
}
|
||||
|
||||
for (auto row = m_matrix.begin(); row != m_matrix.end(); ++row) {
|
||||
for (auto column = row->begin(); column != row->end(); ++column) {
|
||||
for (size_t i = 0; i < block_type::dimension; ++i) {
|
||||
for (size_t j = 0; j < block_type::dimension; ++j) {
|
||||
matrixBuilder[row.index()][column.index()][i][j]
|
||||
= field_type_to(m_matrix[row.index()][column.index()][i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return matrixBuilder;
|
||||
}
|
||||
const M& m_matrix;
|
||||
matrix_type_to m_convertedMatrix;
|
||||
//! \brief the underlying preconditioner to use
|
||||
std::shared_ptr<CudaPreconditionerType> m_underlyingPreconditioner;
|
||||
};
|
||||
} // end namespace Opm::cuistl
|
||||
|
||||
#endif
|
||||
192
tests/cuistl/test_converttofloatadapter.cpp
Normal file
192
tests/cuistl/test_converttofloatadapter.cpp
Normal file
@ -0,0 +1,192 @@
|
||||
/*
|
||||
Copyright 2022-2023 SINTEF AS
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#include <config.h>
|
||||
|
||||
#define BOOST_TEST_MODULE TestConvertToFloatAdapter
|
||||
#define BOOST_TEST_NO_MAIN
|
||||
|
||||
|
||||
#include <boost/mpl/list.hpp>
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <dune/common/parallel/mpihelper.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <limits>
|
||||
#include <memory>
|
||||
#include <opm/simulators/linalg/cuistl/PreconditionerConvertFieldTypeAdapter.hpp>
|
||||
|
||||
|
||||
using XDouble = Dune::BlockVector<Dune::FieldVector<double, 2>>;
|
||||
using MDouble = Dune::FieldMatrix<double, 2, 2>;
|
||||
using SpMatrixDouble = Dune::BCRSMatrix<MDouble>;
|
||||
using XFloat = Dune::BlockVector<Dune::FieldVector<float, 2>>;
|
||||
using MFloat = Dune::FieldMatrix<float, 2, 2>;
|
||||
using SpMatrixFloat = Dune::BCRSMatrix<MFloat>;
|
||||
namespace
|
||||
{
|
||||
class TestPreconditioner : Dune::PreconditionerWithUpdate<XFloat, XFloat>
|
||||
{
|
||||
public:
|
||||
using range_type = XFloat;
|
||||
using domain_type = XFloat;
|
||||
TestPreconditioner(const SpMatrixFloat& matrix,
|
||||
const XDouble& expectedInput,
|
||||
const SpMatrixDouble& expectedMatrix,
|
||||
const XDouble& expectedOutputVector)
|
||||
: m_matrix(matrix)
|
||||
, m_expectedInput(expectedInput)
|
||||
, m_expectedMatrix(expectedMatrix)
|
||||
, m_expectedOutputVector(expectedOutputVector)
|
||||
{
|
||||
}
|
||||
|
||||
virtual void pre([[maybe_unused]] XFloat& x, [[maybe_unused]] XFloat& b) override
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
virtual void apply([[maybe_unused]] XFloat& v, const XFloat& d) override
|
||||
{
|
||||
// Make sure the correct input is copied
|
||||
for (size_t i = 0; i < d.N(); ++i) {
|
||||
for (size_t j = 0; j < d[i].N(); ++j) {
|
||||
BOOST_CHECK_EQUAL(d[i][j], float(m_expectedInput[i][j]));
|
||||
v[i][j] = float(m_expectedOutputVector[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
// make sure we get the correct matrix
|
||||
BOOST_CHECK_EQUAL(m_expectedMatrix.N(), m_matrix.N());
|
||||
BOOST_CHECK_EQUAL(m_expectedMatrix.nonzeroes(), m_matrix.nonzeroes());
|
||||
|
||||
for (auto row = m_matrix.begin(); row != m_matrix.end(); ++row) {
|
||||
for (auto column = row->begin(); column != row->end(); ++column) {
|
||||
for (int i = 0; i < MFloat::rows; ++i) {
|
||||
for (int j = 0; j < MFloat::cols; ++j) {
|
||||
BOOST_CHECK_EQUAL(float(m_expectedMatrix[i][j][i][j]), m_matrix[i][j][i][j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
virtual void post([[maybe_unused]] XFloat& x) override
|
||||
{
|
||||
}
|
||||
|
||||
virtual void update() override
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
//! Category of the preconditioner (see SolverCategory::Category)
|
||||
virtual Dune::SolverCategory::Category category() const override
|
||||
{
|
||||
return Dune::SolverCategory::sequential;
|
||||
}
|
||||
|
||||
static constexpr bool shouldCallPre()
|
||||
{
|
||||
return false;
|
||||
}
|
||||
static constexpr bool shouldCallPost()
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
const SpMatrixFloat& m_matrix;
|
||||
const XDouble& m_expectedInput;
|
||||
const SpMatrixDouble& m_expectedMatrix;
|
||||
const XDouble& m_expectedOutputVector;
|
||||
};
|
||||
} // namespace
|
||||
|
||||
using NumericTypes = boost::mpl::list<double, float>;
|
||||
|
||||
BOOST_AUTO_TEST_CASE(TestFiniteDifference1D)
|
||||
{
|
||||
|
||||
|
||||
const int N = 5;
|
||||
const int nonZeroes = N * 3 - 2;
|
||||
|
||||
SpMatrixDouble B(N, N, nonZeroes, SpMatrixDouble::row_wise);
|
||||
for (auto row = B.createbegin(); row != B.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() < B.N() - 1) {
|
||||
row.insert(row.index() + 1);
|
||||
}
|
||||
}
|
||||
// This might not be the most elegant way of filling in a Dune sparse matrix, but it works.
|
||||
for (int i = 0; i < N; ++i) {
|
||||
B[i][i] = -2;
|
||||
if (i < N - 1) {
|
||||
B[i][i + 1] = 1;
|
||||
}
|
||||
|
||||
if (i > 0) {
|
||||
B[i][i - 1] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// check for the standard basis {e_i}
|
||||
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
|
||||
for (int i = 0; i < N; ++i) {
|
||||
XDouble inputVector(N);
|
||||
XDouble outputVector(N);
|
||||
XDouble expectedOutputVector(N);
|
||||
expectedOutputVector[i][0] = 42.0;
|
||||
expectedOutputVector[i][1] = 43.0;
|
||||
inputVector[i][0] = 1.0;
|
||||
auto converter
|
||||
= Opm::cuistl::PreconditionerConvertFieldTypeAdapter<TestPreconditioner, SpMatrixDouble, XDouble, XDouble>(
|
||||
B);
|
||||
auto underlyingPreconditioner = std::make_shared<TestPreconditioner>(
|
||||
converter.getConvertedMatrix(), inputVector, B, expectedOutputVector);
|
||||
converter.setUnderlyingPreconditioner(underlyingPreconditioner);
|
||||
converter.apply(outputVector, inputVector);
|
||||
|
||||
for (size_t j = 0; j < expectedOutputVector.N(); ++j) {
|
||||
for (size_t k = 0; k < expectedOutputVector[i].N(); ++k) {
|
||||
BOOST_CHECK_EQUAL(outputVector[j][k], float(expectedOutputVector[j][k]));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
bool
|
||||
init_unit_test_func()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
int
|
||||
main(int argc, char** argv)
|
||||
{
|
||||
[[maybe_unused]] const auto& helper = Dune::MPIHelper::instance(argc, argv);
|
||||
boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv);
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user