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

Added CuSeqILU0 with the preconditioner adapter.

Browse Source
This commit is contained in:
Kjetil Olsen Lye 2023-03-31 13:36:27 +02:00
parent ed91f7528b
commit b30e6d79d5
6 changed files with 883 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
CMakeLists_files.cmake
View File

@ -150,7 +150,7 @@ if(CUDA_FOUND)
list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/cuistl/CuVector.cpp)
list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/cuistl/detail/vector_operations.cu)
list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/cuistl/CuSparseMatrix.cpp)
list (APPEND MAIN_SOURCE_FILES opm/simulators/linalg/cuistl/CuSeqILU0.cpp)
# CUISTL HEADERS
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/cuda_safe_call.hpp)
@ -164,12 +164,14 @@ if(CUDA_FOUND)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/CuMatrixDescription.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/CuSparseResource.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/CuSparseResource_impl.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/safe_conversion.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/safe_conversion.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/cublas_wrapper.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/cusparse_wrapper.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/cusparse_constants.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/cuistl/detail/vector_operations.hpp)
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)
endif()
@ -258,6 +260,7 @@ if(CUDA_FOUND)
list(APPEND TEST_SOURCE_FILES tests/cuistl/test_cuvector.cpp)
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)
endif()
if(OPENCL_FOUND)

364
opm/simulators/linalg/cuistl/CuSeqILU0.cpp Normal file
View File

@ -0,0 +1,364 @@
/*
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 <cuda.h>
#include <cuda_runtime.h>
#include <cusparse.h>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <fmt/core.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/cuistl/CuSeqILU0.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/cusparse_wrapper.hpp>
#include <opm/simulators/linalg/cuistl/detail/fix_zero_diagonal.hpp>
#include <opm/simulators/linalg/cuistl/detail/safe_conversion.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
// This file is based on the guide at https://docs.nvidia.com/cuda/cusparse/index.html#csrilu02_solve ,
// it highly recommended to read that before proceeding.
namespace Opm::cuistl
{
template <class M, class X, class Y, int l>
CuSeqILU0<M, X, Y, l>::CuSeqILU0(const M& A, field_type w)
: m_underlyingMatrix(A)
, m_w(w)
, m_LU(CuSparseMatrix<field_type>::fromMatrix(detail::makeMatrixWithNonzeroDiagonal(A)))
, m_temporaryStorage(m_LU.N() * m_LU.blockSize())
, m_descriptionL(detail::createLowerDiagonalDescription())
, m_descriptionU(detail::createUpperDiagonalDescription())
, m_cuSparseHandle(detail::CuSparseHandle::getInstance())
{
// Some sanity check
OPM_ERROR_IF(A.N() != m_LU.N(),
fmt::format("CuSparse matrix not same size as DUNE matrix. {} vs {}.", m_LU.N(), A.N()));
OPM_ERROR_IF(
A[0][0].N() != m_LU.blockSize(),
fmt::format("CuSparse matrix not same blocksize as DUNE matrix. {} vs {}.", m_LU.blockSize(), A[0][0].N()));
OPM_ERROR_IF(
A.N() * A[0][0].N() != m_LU.dim(),
fmt::format("CuSparse matrix not same dimension as DUNE matrix. {} vs {}.", m_LU.dim(), A.N() * A[0][0].N()));
OPM_ERROR_IF(A.nonzeroes() != m_LU.nonzeroes(),
fmt::format("CuSparse matrix not same number of non zeroes as DUNE matrix. {} vs {}. ",
m_LU.nonzeroes(),
A.nonzeroes()));
updateILUConfiguration();
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::pre([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
{
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::apply(X& v, const Y& d)
{
// We need to pass the solve routine a scalar to multiply.
// In our case this scalar is 1.0
const field_type one = 1.0;
const auto numberOfRows = detail::to_int(m_LU.N());
const auto numberOfNonzeroBlocks = detail::to_int(m_LU.nonzeroes());
const auto blockSize = detail::to_int(m_LU.blockSize());
auto nonZeroValues = m_LU.getNonZeroValues().data();
auto rowIndices = m_LU.getRowIndices().data();
auto columnIndices = m_LU.getColumnIndices().data();
// Solve L m_temporaryStorage = d
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_solve(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
&one,
m_descriptionL->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoL.get(),
d.data(),
m_temporaryStorage.data(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
// Solve U v = m_temporaryStorage
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_solve(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
&one,
m_descriptionU->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoU.get(),
m_temporaryStorage.data(),
v.data(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
v *= m_w;
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::post([[maybe_unused]] X& x)
{
}
template <class M, class X, class Y, int l>
Dune::SolverCategory::Category
CuSeqILU0<M, X, Y, l>::category() const
{
return Dune::SolverCategory::sequential;
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::update()
{
m_LU.updateNonzeroValues(detail::makeMatrixWithNonzeroDiagonal(m_underlyingMatrix));
// updateILUConfiguration();
createILU();
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::analyzeMatrix()
{
if (!m_buffer) {
OPM_THROW(std::runtime_error,
"Buffer not initialized. Call findBufferSize() then initialize with the appropiate size.");
}
const auto numberOfRows = detail::to_int(m_LU.N());
const auto numberOfNonzeroBlocks = detail::to_int(m_LU.nonzeroes());
const auto blockSize = detail::to_int(m_LU.blockSize());
auto nonZeroValues = m_LU.getNonZeroValues().data();
auto rowIndices = m_LU.getRowIndices().data();
auto columnIndices = m_LU.getColumnIndices().data();
// analysis of ilu LU decomposition
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrilu02_analysis(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
numberOfRows,
numberOfNonzeroBlocks,
m_LU.getDescription().get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoM.get(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
// Make sure we can decompose the matrix.
int structuralZero;
auto statusPivot = cusparseXbsrilu02_zeroPivot(m_cuSparseHandle.get(), m_infoM.get(), &structuralZero);
OPM_ERROR_IF(statusPivot != CUSPARSE_STATUS_SUCCESS,
fmt::format("Found a structucal zero at A({}, {}). Could not decompuse LU \\approx A.\n\n A has "
"dimension {}, and has {} nonzeroes.",
structuralZero,
structuralZero,
m_LU.N(),
m_LU.nonzeroes()));
// analysis of ilu apply
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_analysis(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
m_descriptionL->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoL.get(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_analysis(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
m_descriptionU->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoU.get(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
m_analyzisDone = true;
}
template <class M, class X, class Y, int l>
size_t
CuSeqILU0<M, X, Y, l>::findBufferSize()
{
// We have three calls that need buffers:
// 1) LU decomposition
// 2) solve Lv = y
// 3) solve Ux = z
// we combine these buffers into one since it is not used across calls,
const auto numberOfRows = detail::to_int(m_LU.N());
const auto numberOfNonzeroBlocks = detail::to_int(m_LU.nonzeroes());
const auto blockSize = detail::to_int(m_LU.blockSize());
auto nonZeroValues = m_LU.getNonZeroValues().data();
auto rowIndices = m_LU.getRowIndices().data();
auto columnIndices = m_LU.getColumnIndices().data();
int bufferSizeM = 0;
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrilu02_bufferSize(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
numberOfRows,
numberOfNonzeroBlocks,
m_LU.getDescription().get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoM.get(),
&bufferSizeM));
int bufferSizeL = 0;
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_bufferSize(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
m_descriptionL->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoL.get(),
&bufferSizeL));
int bufferSizeU = 0;
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrsv2_bufferSize(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
CUSPARSE_OPERATION_NON_TRANSPOSE,
numberOfRows,
numberOfNonzeroBlocks,
m_descriptionL->get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoU.get(),
&bufferSizeU));
OPM_ERROR_IF(bufferSizeL <= 0, fmt::format("bufferSizeL is non-positive. Given value is {}.", bufferSizeL));
OPM_ERROR_IF(bufferSizeU <= 0, fmt::format("bufferSizeU is non-positive. Given value is {}.", bufferSizeU));
OPM_ERROR_IF(bufferSizeM <= 0, fmt::format("bufferSizeM is non-positive. Given value is {}.", bufferSizeM));
return size_t(std::max(bufferSizeL, std::max(bufferSizeU, bufferSizeM)));
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::createILU()
{
OPM_ERROR_IF(!m_buffer, "Buffer not initialized. Call findBufferSize() then initialize with the appropiate size.");
OPM_ERROR_IF(!m_analyzisDone, "Analyzis of matrix not done. Call analyzeMatrix() first.");
const auto numberOfRows = detail::to_int(m_LU.N());
const auto numberOfNonzeroBlocks = detail::to_int(m_LU.nonzeroes());
const auto blockSize = detail::to_int(m_LU.blockSize());
auto nonZeroValues = m_LU.getNonZeroValues().data();
auto rowIndices = m_LU.getRowIndices().data();
auto columnIndices = m_LU.getColumnIndices().data();
OPM_CUSPARSE_SAFE_CALL(detail::cusparseBsrilu02(m_cuSparseHandle.get(),
detail::CUSPARSE_MATRIX_ORDER,
numberOfRows,
numberOfNonzeroBlocks,
m_LU.getDescription().get(),
nonZeroValues,
rowIndices,
columnIndices,
blockSize,
m_infoM.get(),
CUSPARSE_SOLVE_POLICY_USE_LEVEL,
m_buffer->data()));
// We need to do this here as well. The first call was to check that we could decompose the system A=LU
// the second call here is to make sure we can solve LUx=y
int structuralZero;
// cusparseXbsrilu02_zeroPivot() calls cudaDeviceSynchronize()
auto statusPivot = cusparseXbsrilu02_zeroPivot(m_cuSparseHandle.get(), m_infoM.get(), &structuralZero);
OPM_ERROR_IF(
statusPivot != CUSPARSE_STATUS_SUCCESS,
fmt::format("Found a structucal zero at LU({}, {}). Could not solve LUx = y.", structuralZero, structuralZero));
}
template <class M, class X, class Y, int l>
void
CuSeqILU0<M, X, Y, l>::updateILUConfiguration()
{
auto bufferSize = findBufferSize();
if (!m_buffer || m_buffer->dim() < bufferSize) {
m_buffer.reset(new CuVector<field_type>((bufferSize + sizeof(field_type) - 1) / sizeof(field_type)));
}
analyzeMatrix();
createILU();
}
} // namespace Opm::cuistl
#define INSTANTIATE_CUSEQILU0_DUNE(realtype, blockdim) \
template class ::Opm::cuistl::CuSeqILU0<Dune::BCRSMatrix<Dune::FieldMatrix<realtype, blockdim, blockdim>>, \
::Opm::cuistl::CuVector<realtype>, \
::Opm::cuistl::CuVector<realtype>>; \
template class ::Opm::cuistl::CuSeqILU0<Dune::BCRSMatrix<Opm::MatrixBlock<realtype, blockdim, blockdim>>, \
::Opm::cuistl::CuVector<realtype>, \
::Opm::cuistl::CuVector<realtype>>
INSTANTIATE_CUSEQILU0_DUNE(double, 1);
INSTANTIATE_CUSEQILU0_DUNE(double, 2);
INSTANTIATE_CUSEQILU0_DUNE(double, 3);
INSTANTIATE_CUSEQILU0_DUNE(double, 4);
INSTANTIATE_CUSEQILU0_DUNE(double, 5);
INSTANTIATE_CUSEQILU0_DUNE(double, 6);
INSTANTIATE_CUSEQILU0_DUNE(float, 1);
INSTANTIATE_CUSEQILU0_DUNE(float, 2);
INSTANTIATE_CUSEQILU0_DUNE(float, 3);
INSTANTIATE_CUSEQILU0_DUNE(float, 4);
INSTANTIATE_CUSEQILU0_DUNE(float, 5);
INSTANTIATE_CUSEQILU0_DUNE(float, 6);

126
opm/simulators/linalg/cuistl/PreconditionerAdapter.hpp Normal file
View File

@ -0,0 +1,126 @@
/*
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_PRECONDITIONERADAPTER_HPP
#define OPM_PRECONDITIONERADAPTER_HPP
#include <cusparse.h>
#include <dune/istl/preconditioner.hh>
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
#include <opm/simulators/linalg/cuistl/CuVector.hpp>
#include <opm/simulators/linalg/cuistl/detail/preconditioner_should_call_post_pre.hpp>
namespace Opm::cuistl
{
//!\brief Makes a CUDA preconditioner available to a CPU simulator.
//!
//! The use case for this adapter is to use a CUDA preconditioner during a linear
//! solver that works on the CPU. The motivation for this is benchmarking new preconditioners on the GPU.
//!
//! \tparam X the domain type (should be on the CPU). Typicall a Dune::BlockVector
//! \tparam Y the range type (should be on the CPU). Typicall a Dune::BlockVector
//! \tparam CudaPreconditionerType the preconditioner taking CuVector<real_type> as arguments to apply
template <class X, class Y, class CudaPreconditionerType>
class PreconditionerAdapter : public Dune::PreconditionerWithUpdate<X, Y>
{
public:
//! \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;
//! \brief Constructor.
//!
//! Constructor gets all parameters to operate the prec.
//! \param A The matrix to operate on.
//! \param w The relaxation factor.
//!
explicit PreconditionerAdapter(std::shared_ptr<CudaPreconditionerType> preconditioner)
: m_underlyingPreconditioner(preconditioner)
{
}
//! \brief Prepare the preconditioner.
//!
//! Currently not supported.
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.
//!
//! \copydoc Preconditioner::apply(X&,const Y&)
virtual void apply(X& v, const Y& d) override
{
if (!m_inputBuffer) {
m_inputBuffer.reset(new CuVector<field_type>(v.dim()));
m_outputBuffer.reset(new CuVector<field_type>(v.dim()));
}
m_inputBuffer->copyFromHost(d);
m_underlyingPreconditioner->apply(*m_outputBuffer, *m_inputBuffer);
m_outputBuffer->copyToHost(v);
}
//! \brief Clean up.
//!
//! Currently not supported.
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
{
return m_underlyingPreconditioner->category();
}
//! Calls update on the underlying CUDA preconditioner
virtual void update() override
{
m_underlyingPreconditioner->update();
}
static constexpr bool shouldCallPre()
{
return detail::shouldCallPreconditionerPost<CudaPreconditionerType>();
}
static constexpr bool shouldCallPost()
{
return detail::shouldCallPreconditionerPre<CudaPreconditionerType>();
}
private:
//! \brief the underlying preconditioner to use
std::shared_ptr<CudaPreconditionerType> m_underlyingPreconditioner;
std::unique_ptr<CuVector<field_type>> m_inputBuffer;
std::unique_ptr<CuVector<field_type>> m_outputBuffer;
};
} // end namespace Opm::cuistl
#endif

88
opm/simulators/linalg/cuistl/detail/has_function.hpp Normal file
View File

@ -0,0 +1,88 @@
/*
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_CUISTL_DETAIL_HAS_FUNCTION_HPP
#define OPM_CUISTL_DETAIL_HAS_FUNCTION_HPP
#include <type_traits>
/**
* Simple utility structs to test for the existence of functions in types.
*
* Note that there are alternatives to this, see for instance
* https://stackoverflow.com/questions/257288/templated-check-for-the-existence-of-a-class-member-function , however,
* this is by far the cleanest approach for where this is going to be used for now.
*
* TODO: Use the requires-keyword once C++20 becomes availble (https://en.cppreference.com/w/cpp/language/constraints ).
* With C++20 this file can be removed.
*/
namespace Opm::cuistl::detail
{
/**
* @brief The has_should_call_pre class detects the presence of the method shouldCallPre
*
* Usage:
*
* @code{.cpp}
* if constexpr (has_should_call_pre<decltype(preconditioner)>::value) {
* // We know that the function shouldCallPre is present:
* auto shouldCallPre = preconditioner.shouldCallPre();
* }
* @endcode
*
* @note This is mainly done in the GPU preconditioner to avoid having to copy data in the pre step.
*/
template <typename T>
class has_should_call_pre
{
template <typename U>
static std::true_type test(decltype(&U::shouldCallPre));
template <typename U>
static std::false_type test(...);
public:
static constexpr bool value = std::is_same_v<decltype(test<T>(0)), std::true_type>;
};
/**
* @brief The has_should_call_post class detects the presence of the method shouldCallPost
*
* Usage:
*
* @code{.cpp}
* if constexpr (has_should_call_post<decltype(preconditioner)>::value) {
* // We know that the function shouldCallPost is present:
* auto shouldCallPost = preconditioner.shouldCallPost();
* }
* @endcode
*
* @note This is mainly done in the GPU preconditioner to avoid having to copy data in the post step.
*/
template <typename T>
class has_should_call_post
{
template <typename U>
static std::true_type test(decltype(&U::shouldCallPost));
template <typename U>
static std::false_type test(...);
public:
static constexpr bool value = std::is_same_v<decltype(test<T>(0)), std::true_type>;
};
} // namespace Opm::cuistl::detail
#endif

64
opm/simulators/linalg/cuistl/detail/preconditioner_should_call_post_pre.hpp Normal file
View File

@ -0,0 +1,64 @@
/*
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_CUISTL_PRECONDIDTIONER_SHOULD_CALL_POST_PRE_HPP
#define OPM_CUISTL_PRECONDIDTIONER_SHOULD_CALL_POST_PRE_HPP
#include <opm/simulators/linalg/cuistl/detail/has_function.hpp>
namespace Opm::cuistl::detail
{
//! @brief Tests (compile time) if the preconditioner type needs to call pre
//!
//! @note This is mostly used to avoid unneeded copying back and front to the GPU, as well
//! as avoiding communication.
template <class PreconditionerType>
constexpr bool
shouldCallPreconditionerPre()
{
if constexpr (has_should_call_pre<PreconditionerType>::value) {
return PreconditionerType::shouldCallPre();
} else {
// If the preconditioner type does not have a way of signalling the need
// to call pre(), we should always call it. This is because it could be an
// old design.
return true;
}
}
//! @brief Tests (compile time) if the preconditioner type needs to call pre
//!
//! @note This is mostly used to avoid unneeded copying back and front to the GPU, as well
//! as avoiding communication.
template <class PreconditionerType>
constexpr bool
shouldCallPreconditionerPost()
{
if constexpr (has_should_call_post<PreconditionerType>::value) {
return PreconditionerType::shouldCallPost();
} else {
// If the preconditioner type does not have a way of signalling the need
// to call post(), we should always call it. This is because it could be an
// old design.
return true;
}
}
} // namespace Opm::cuistl::detail
#endif

235
tests/cuistl/test_cuseqilu0.cpp Normal file
View File

@ -0,0 +1,235 @@
/*
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 TestCuSeqILU0
#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 <opm/simulators/linalg/cuistl/CuSeqILU0.hpp>
#include <opm/simulators/linalg/cuistl/CuVector.hpp>
#include <opm/simulators/linalg/cuistl/PreconditionerAdapter.hpp>
#include <opm/simulators/linalg/cuistl/detail/cuda_safe_call.hpp>
#include <limits>
#include <memory>
using NumericTypes = boost::mpl::list<double, float>;
BOOST_AUTO_TEST_CASE_TEMPLATE(TestFiniteDifference1D, T, NumericTypes)
{
// Here we will test a simple 1D finite difference scheme for
// the Laplace equation:
//
// -\Delta u = f on [0,1]
//
// Using a central difference approximation of \Delta u, this can
// be approximated by
//
// -(u_{i+1}-2u_i+u_{i-1})/Dx^2 = f(x_i)
//
// giving rise to the matrix
//
// -2 1 0 0 ... 0 0
// 1 -2 1 0 0 ... 0
// ....
// 0 0 0 ...1 -2 1
// 0 0 0 ... 1 -2
const int N = 5;
const int nonZeroes = N * 3 - 2;
using M = Dune::FieldMatrix<T, 1, 1>;
using SpMatrix = Dune::BCRSMatrix<M>;
using Vector = Dune::BlockVector<Dune::FieldVector<T, 1>>;
using CuILU0 = Opm::cuistl::CuSeqILU0<SpMatrix, Opm::cuistl::CuVector<T>, Opm::cuistl::CuVector<T>>;
SpMatrix B(N, N, nonZeroes, SpMatrix::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;
}
}
auto duneILU = Dune::SeqILU<SpMatrix, Vector, Vector>(B, 1.0);
auto cuILU = Opm::cuistl::PreconditionerAdapter<Vector, Vector, CuILU0>(std::make_shared<CuILU0>(B, 1.0));
// check for the standard basis {e_i}
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
for (int i = 0; i < N; ++i) {
Vector inputVector(N);
inputVector[i][0] = 1.0;
Vector outputVectorDune(N);
Vector outputVectorCuistl(N);
duneILU.apply(outputVectorDune, inputVector);
cuILU.apply(outputVectorCuistl, inputVector);
for (int component = 0; component < N; ++component) {
BOOST_CHECK_CLOSE(outputVectorDune[component][0],
outputVectorCuistl[component][0],
std::numeric_limits<T>::epsilon() * 1000);
}
}
// Now we check that we can update the matrix. We basically just negate B
B *= -1.0;
auto duneILUNew = Dune::SeqILU<SpMatrix, Vector, Vector>(B, 1.0);
cuILU.update();
// check for the standard basis {e_i}
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
for (int i = 0; i < N; ++i) {
Vector inputVector(N);
inputVector[i][0] = 1.0;
Vector outputVectorDune(N);
Vector outputVectorCuistl(N);
duneILUNew.apply(outputVectorDune, inputVector);
cuILU.apply(outputVectorCuistl, inputVector);
for (int component = 0; component < N; ++component) {
BOOST_CHECK_CLOSE(outputVectorDune[component][0],
outputVectorCuistl[component][0],
std::numeric_limits<T>::epsilon() * 1000);
}
}
}
BOOST_AUTO_TEST_CASE_TEMPLATE(TestFiniteDifferenceBlock2, T, NumericTypes)
{
// Here we will test a simple 1D finite difference scheme for
// the Laplace equation:
//
// -\Delta u = f on [0,1]
//
// Using a central difference approximation of \Delta u, this can
// be approximated by
//
// -(u_{i+1}-2u_i+u_{i-1})/Dx^2 = f(x_i)
//
// giving rise to the matrix
//
// -2 1 0 0 ... 0 0
// 1 -2 1 0 0 ... 0
// ....
// 0 0 0 ...1 -2 1
// 0 0 0 ... 1 -2
const int N = 5;
const int nonZeroes = N * 3 - 2;
using M = Dune::FieldMatrix<T, 2, 2>;
using SpMatrix = Dune::BCRSMatrix<M>;
using Vector = Dune::BlockVector<Dune::FieldVector<T, 2>>;
using CuILU0 = Opm::cuistl::CuSeqILU0<SpMatrix, Opm::cuistl::CuVector<T>, Opm::cuistl::CuVector<T>>;
SpMatrix B(N, N, nonZeroes, SpMatrix::row_wise);
for (auto row = B.createbegin(); row != B.createend(); ++row) {
row.insert(row.index());
if (row.index() < N - 1) {
row.insert(row.index() + 1);
}
if (row.index() > 0) {
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][0][0] = -2;
B[i][i][1][1] = -2;
B[i][i][0][1] = 1;
B[i][i][1][0] = 1;
}
auto duneILU = Dune::SeqILU<SpMatrix, Vector, Vector>(B, 1.0);
auto cuILU = Opm::cuistl::PreconditionerAdapter<Vector, Vector, CuILU0>(std::make_shared<CuILU0>(B, 1.0));
// check for the standard basis {e_i}
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
for (int i = 0; i < N; ++i) {
Vector inputVector(N);
inputVector[i][0] = 1.0;
Vector outputVectorDune(N);
Vector outputVectorCuistl(N);
duneILU.apply(outputVectorDune, inputVector);
cuILU.apply(outputVectorCuistl, inputVector);
for (int component = 0; component < N; ++component) {
BOOST_CHECK_CLOSE(outputVectorDune[component][0],
outputVectorCuistl[component][0],
std::numeric_limits<T>::epsilon() * 1000);
}
}
// Now we check that we can update the matrix. We basically just negate B
B *= -1.0;
auto duneILUNew = Dune::SeqILU<SpMatrix, Vector, Vector>(B, 1.0);
cuILU.update();
// check for the standard basis {e_i}
// (e_i=(0,...,0, 1 (i-th place), 0, ..., 0))
for (int i = 0; i < N; ++i) {
Vector inputVector(N);
inputVector[i][0] = 1.0;
Vector outputVectorDune(N);
Vector outputVectorCuistl(N);
duneILUNew.apply(outputVectorDune, inputVector);
cuILU.apply(outputVectorCuistl, inputVector);
for (int component = 0; component < N; ++component) {
BOOST_CHECK_CLOSE(outputVectorDune[component][0],
outputVectorCuistl[component][0],
std::numeric_limits<T>::epsilon() * 1000);
}
}
}
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);
}