Preconditioner: template Scalar type

This commit is contained in:
Arne Morten Kvarving 2024-04-15 22:38:04 +02:00
parent 05a89d1e96
commit b75ea188ee
10 changed files with 73 additions and 105 deletions

View File

@ -31,17 +31,14 @@
#include <sstream>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
using Opm::OpmLog;
using Dune::Timer;
template <unsigned int block_size>
BILU0<block_size>::BILU0(bool opencl_ilu_parallel_, int verbosity_) :
Preconditioner<block_size>(verbosity_), opencl_ilu_parallel(opencl_ilu_parallel_)
BILU0<block_size>::BILU0(bool opencl_ilu_parallel_, int verbosity_)
: Base(verbosity_)
, opencl_ilu_parallel(opencl_ilu_parallel_)
{
#if CHOW_PATEL
chowPatelIlu.setVerbosity(verbosity);
@ -319,5 +316,4 @@ INSTANTIATE_BDA_FUNCTIONS(6);
#undef INSTANTIATE_BDA_FUNCTIONS
} // namespace Accelerator
} // namespace Opm
} // namespace Opm::Accelerator

View File

@ -29,18 +29,15 @@
#include <opm/simulators/linalg/bda/opencl/ChowPatelIlu.hpp>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
/// This class implements a Blocked ILU0 preconditioner
/// The decomposition is done on GPU, using exact decomposition, or ChowPatel decomposition
/// The preconditioner is applied via two exact triangular solves
template <unsigned int block_size>
class BILU0 : public Preconditioner<block_size>
class BILU0 : public Preconditioner<double,block_size>
{
typedef Preconditioner<block_size> Base;
using Base = Preconditioner<double,block_size>;
using Base::N;
using Base::Nb;
@ -53,9 +50,9 @@ class BILU0 : public Preconditioner<block_size>
using Base::err;
private:
std::unique_ptr<BlockedMatrix<double>> LUmat = nullptr;
std::unique_ptr<BlockedMatrix<double>> LUmat{};
#if CHOW_PATEL
std::unique_ptr<BlockedMatrix<double>> Lmat = nullptr, Umat = nullptr;
std::unique_ptr<BlockedMatrix<double>> Lmat{}, Umat{};
#endif
std::vector<double> invDiagVals;
std::vector<int> diagIndex;
@ -122,8 +119,7 @@ public:
}
};
} // namespace Accelerator
} // namespace Opm
} // namespace Opm::Accelerator
#endif

View File

@ -34,17 +34,14 @@
#include <sstream>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
using Opm::OpmLog;
using Dune::Timer;
template <unsigned int block_size>
BISAI<block_size>::BISAI(bool opencl_ilu_parallel_, int verbosity_) :
Preconditioner<block_size>(verbosity_)
BISAI<block_size>::BISAI(bool opencl_ilu_parallel_, int verbosity_)
: Base(verbosity_)
{
#if CHOW_PATEL
OPM_THROW(std::logic_error, "Error --linear-solver=isai cannot be used if ChowPatelIlu is used, probably defined by CMake\n");
@ -312,5 +309,4 @@ INSTANTIATE_BDA_FUNCTIONS(6);
#undef INSTANTIATE_BDA_FUNCTIONS
}
}
} // namespace Opm::Accelerator

View File

@ -26,19 +26,16 @@
#include <opm/simulators/linalg/bda/opencl/BILU0.hpp>
#include <opm/simulators/linalg/bda/opencl/Preconditioner.hpp>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
template<class Scalar> class BlockedMatrix;
/// This class implements a Blocked version of the Incomplete Sparse Approximate Inverse (ISAI) preconditioner.
/// Inspired by the paper "Incomplete Sparse Approximate Inverses for Parallel Preconditioning" by Anzt et. al.
template <unsigned int block_size>
class BISAI : public Preconditioner<block_size>
class BISAI : public Preconditioner<double,block_size>
{
typedef Preconditioner<block_size> Base;
using Base = Preconditioner<double,block_size>;
using Base::N;
using Base::Nb;
@ -134,7 +131,6 @@ public:
/// in the csrToCscOffsetMap[i]-th position in the CSC representation.
std::vector<int> buildCsrToCscOffsetMap(std::vector<int> colPointers, std::vector<int> rowIndices);
} // namespace Accelerator
} // namespace Opm
} // namespace Opm::Accelerator
#endif

View File

@ -35,17 +35,15 @@
#include <opm/simulators/linalg/bda/opencl/openclKernels.hpp>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
using Opm::OpmLog;
using Dune::Timer;
template <unsigned int block_size>
CPR<block_size>::CPR(bool opencl_ilu_parallel_, int verbosity_) :
Preconditioner<block_size>(verbosity_), opencl_ilu_parallel(opencl_ilu_parallel_)
CPR<block_size>::CPR(bool opencl_ilu_parallel_, int verbosity_)
: Base(verbosity_)
, opencl_ilu_parallel(opencl_ilu_parallel_)
{
bilu0 = std::make_unique<BILU0<block_size> >(opencl_ilu_parallel, verbosity_);
diagIndices.resize(1);
@ -570,7 +568,6 @@ INSTANTIATE_BDA_FUNCTIONS(6);
#undef INSTANTIATE_BDA_FUNCTIONS
} // namespace Accelerator
} // namespace Opm
} // namespace Opm::Accelerator

View File

@ -33,18 +33,15 @@
#include <opm/simulators/linalg/bda/opencl/openclSolverBackend.hpp>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
template<class Scalar> class BlockedMatrix;
/// This class implements a Constrained Pressure Residual (CPR) preconditioner
template <unsigned int block_size>
class CPR : public Preconditioner<block_size>
class CPR : public Preconditioner<double,block_size>
{
typedef Preconditioner<block_size> Base;
using Base = Preconditioner<double,block_size>;
using Base::N;
using Base::Nb;
@ -138,8 +135,7 @@ public:
// x and b are vectors with 3 elements
void solve_transposed_3x3(const double *A, const double *b, double *x);
} // namespace Accelerator
} // namespace Opm
} // namespace Opm::Accelerator
#endif

View File

@ -30,21 +30,20 @@
#include <memory>
#include <string>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
template <unsigned int block_size>
void Preconditioner<block_size>::setOpencl(std::shared_ptr<cl::Context>& context_, std::shared_ptr<cl::CommandQueue>& queue_) {
template<class Scalar, unsigned int block_size>
void Preconditioner<Scalar,block_size>::
setOpencl(std::shared_ptr<cl::Context>& context_,
std::shared_ptr<cl::CommandQueue>& queue_)
{
context = context_;
queue = queue_;
}
template <unsigned int block_size>
std::unique_ptr<Preconditioner<block_size>>
Preconditioner<block_size>::create(Type type, bool opencl_ilu_parallel, int verbosity)
template<class Scalar, unsigned int block_size>
std::unique_ptr<Preconditioner<Scalar,block_size>>
Preconditioner<Scalar,block_size>::create(Type type, bool opencl_ilu_parallel, int verbosity)
{
switch (type ) {
case Type::BILU0:
@ -59,37 +58,30 @@ Preconditioner<block_size>::create(Type type, bool opencl_ilu_parallel, int verb
"Invalid preconditioner type " + std::to_string(static_cast<int>(type)));
}
template <unsigned int block_size>
bool Preconditioner<block_size>::
analyze_matrix(BlockedMatrix<double>* mat,
[[maybe_unused]] BlockedMatrix<double>* jacMat)
template<class Scalar, unsigned int block_size>
bool Preconditioner<Scalar,block_size>::
analyze_matrix(BlockedMatrix<Scalar>* mat,
[[maybe_unused]] BlockedMatrix<Scalar>* jacMat)
{
return analyze_matrix(mat);
}
template <unsigned int block_size>
bool Preconditioner<block_size>::
create_preconditioner(BlockedMatrix<double>* mat,
[[maybe_unused]] BlockedMatrix<double>* jacMat)
template<class Scalar, unsigned int block_size>
bool Preconditioner<Scalar,block_size>::
create_preconditioner(BlockedMatrix<Scalar>* mat,
[[maybe_unused]] BlockedMatrix<Scalar>* jacMat)
{
return create_preconditioner(mat);
}
#define INSTANTIATE_BDA_FUNCTIONS(n) \
template std::unique_ptr<Preconditioner<n> > Preconditioner<n>::create(Type, bool, int); \
template void Preconditioner<n>::setOpencl(std::shared_ptr<cl::Context>&, std::shared_ptr<cl::CommandQueue>&); \
template bool Preconditioner<n>::analyze_matrix(BlockedMatrix<double> *, BlockedMatrix<double> *); \
template bool Preconditioner<n>::create_preconditioner(BlockedMatrix<double> *, BlockedMatrix<double> *);
#define INSTANCE_TYPE(T) \
template class Preconditioner<T,1>; \
template class Preconditioner<T,2>; \
template class Preconditioner<T,3>; \
template class Preconditioner<T,4>; \
template class Preconditioner<T,5>; \
template class Preconditioner<T,6>;
INSTANTIATE_BDA_FUNCTIONS(1);
INSTANTIATE_BDA_FUNCTIONS(2);
INSTANTIATE_BDA_FUNCTIONS(3);
INSTANTIATE_BDA_FUNCTIONS(4);
INSTANTIATE_BDA_FUNCTIONS(5);
INSTANTIATE_BDA_FUNCTIONS(6);
#undef INSTANTIATE_BDA_FUNCTIONS
} //namespace Accelerator
} //namespace Opm
INSTANCE_TYPE(double)
} // namespace Opm::Accelerator

View File

@ -24,17 +24,13 @@
#include <memory>
namespace Opm
{
namespace Accelerator
{
namespace Opm::Accelerator {
template<class Scalar> class BlockedMatrix;
template <unsigned int block_size>
template<class Scalar, unsigned int block_size>
class Preconditioner
{
protected:
int N = 0; // number of rows of the matrix
int Nb = 0; // number of blockrows of the matrix
@ -65,7 +61,8 @@ public:
virtual ~Preconditioner() = default;
// nested Preconditioners might need to override this
virtual void setOpencl(std::shared_ptr<cl::Context>& context, std::shared_ptr<cl::CommandQueue>& queue);
virtual void setOpencl(std::shared_ptr<cl::Context>& context,
std::shared_ptr<cl::CommandQueue>& queue);
// apply preconditioner, x = prec(y)
virtual void apply(const cl::Buffer& y, cl::Buffer& x) = 0;
@ -73,18 +70,17 @@ public:
// analyze matrix, e.g. the sparsity pattern
// probably only called once
// the version with two params can be overloaded, if not, it will default to using the one param version
virtual bool analyze_matrix(BlockedMatrix<double>* mat) = 0;
virtual bool analyze_matrix(BlockedMatrix<double>* mat,
BlockedMatrix<double>* jacMat);
virtual bool analyze_matrix(BlockedMatrix<Scalar>* mat) = 0;
virtual bool analyze_matrix(BlockedMatrix<Scalar>* mat,
BlockedMatrix<Scalar>* jacMat);
// create/update preconditioner, probably used every linear solve
// the version with two params can be overloaded, if not, it will default to using the one param version
virtual bool create_preconditioner(BlockedMatrix<double>* mat) = 0;
virtual bool create_preconditioner(BlockedMatrix<double>* mat,
BlockedMatrix<double>* jacMat);
virtual bool create_preconditioner(BlockedMatrix<Scalar>* mat) = 0;
virtual bool create_preconditioner(BlockedMatrix<Scalar>* mat,
BlockedMatrix<Scalar>* jacMat);
};
} //namespace Accelerator
} //namespace Opm
} // namespace Opm::Accelerator
#endif

View File

@ -65,13 +65,16 @@ openclSolverBackend<block_size>::openclSolverBackend(int verbosity_, int maxit_,
OPM_THROW(std::logic_error, "Error unknown value for argument --linear-solver, " + linsolver);
}
using PreconditionerType = typename Preconditioner<block_size>::Type;
using PreconditionerType = Preconditioner<double,block_size>;
if (use_cpr) {
prec = Preconditioner<block_size>::create(PreconditionerType::CPR, opencl_ilu_parallel, verbosity);
prec = PreconditionerType::create(PreconditionerType::Type::CPR,
opencl_ilu_parallel, verbosity);
} else if (use_isai) {
prec = Preconditioner<block_size>::create(PreconditionerType::BISAI, opencl_ilu_parallel, verbosity);
prec = PreconditionerType::create(PreconditionerType::Type::BISAI,
opencl_ilu_parallel, verbosity);
} else {
prec = Preconditioner<block_size>::create(PreconditionerType::BILU0, opencl_ilu_parallel, verbosity);
prec = PreconditionerType::create(PreconditionerType::Type::BILU0,
opencl_ilu_parallel, verbosity);
}
std::ostringstream out;

View File

@ -63,7 +63,7 @@ private:
bool useJacMatrix = false;
std::unique_ptr<Preconditioner<block_size> > prec;
std::unique_ptr<Preconditioner<double,block_size> > prec;
// can perform blocked ILU0 and AMG on pressure component
bool is_root; // allow for nested solvers, the root solver is called by BdaBridge
bool analysis_done = false;