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

ISTLSolverBda: use Scalar type from vector

Browse Source
This commit is contained in:
Arne Morten Kvarving 2024-04-16 12:36:10 +02:00
parent 0e537a8ae3
commit 41885f5911
2 changed files with 44 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
opm/simulators/linalg/ISTLSolverBda.cpp
View File

@ -50,15 +50,14 @@
std::shared_ptr<std::thread> copyThread; std::shared_ptr<std::thread> copyThread;
#endif // HAVE_OPENMP #endif // HAVE_OPENMP
namespace Opm { namespace Opm::detail {
namespace detail {
template<class Matrix, class Vector> template<class Matrix, class Vector>
BdaSolverInfo<Matrix,Vector>:: BdaSolverInfo<Matrix,Vector>::
BdaSolverInfo(const std::string& accelerator_mode, BdaSolverInfo(const std::string& accelerator_mode,
const int linear_solver_verbosity, const int linear_solver_verbosity,
const int maxit, const int maxit,
const double tolerance, const Scalar tolerance,
const int platformID, const int platformID,
const int deviceID, const int deviceID,
const bool opencl_ilu_parallel, const bool opencl_ilu_parallel,
@ -104,7 +103,7 @@ apply(Vector& rhs,
{ {
bool use_gpu = bridge_->getUseGpu(); bool use_gpu = bridge_->getUseGpu();
if (use_gpu) { if (use_gpu) {
auto wellContribs = WellContributions<double>::create(accelerator_mode_, useWellConn); auto wellContribs = WellContributions<Scalar>::create(accelerator_mode_, useWellConn);
bridge_->initWellContributions(*wellContribs, x.N() * x[0].N()); bridge_->initWellContributions(*wellContribs, x.N() * x[0].N());
// the WellContributions can only be applied separately with CUDA, OpenCL or rocsparse, not with amgcl or rocalution // the WellContributions can only be applied separately with CUDA, OpenCL or rocsparse, not with amgcl or rocalution
@ -179,8 +178,9 @@ blockJacobiAdjacency(const Grid& grid,
const auto& gridView = grid.leafGridView(); const auto& gridView = grid.leafGridView();
auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows
//Loop over cells // Loop over cells
for (Iter row = blockJacobiForGPUILU0_->createbegin(); row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row) for (Iter row = blockJacobiForGPUILU0_->createbegin();
row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row)
{ {
const auto& elem = *elemIt; const auto& elem = *elemIt;
size_type idx = lid.id(elem); size_type idx = lid.id(elem);
@ -221,25 +221,26 @@ copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
auto outerCol = (*outerRow).begin(); auto outerCol = (*outerRow).begin();
for (auto col = (*row).begin(); col != (*row).end(); ++col) { for (auto col = (*row).begin(); col != (*row).end(); ++col) {
// outerRow is guaranteed to have all column entries that row has! // outerRow is guaranteed to have all column entries that row has!
while(outerCol.index() < col.index()) ++outerCol; while (outerCol.index() < col.index()) {
++outerCol;
}
assert(outerCol.index() == col.index()); assert(outerCol.index() == col.index());
*col = *outerCol; // copy nonzero block *col = *outerCol; // copy nonzero block
} }
} }
} }
template<int Dim> template<class Scalar, int Dim>
using BM = Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>; using BM = Dune::BCRSMatrix<MatrixBlock<Scalar,Dim,Dim>>;
template<int Dim> template<class Scalar, int Dim>
using BV = Dune::BlockVector<Dune::FieldVector<double,Dim>>; using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Dim>>;
#define INSTANTIATE_GRID(T, Dim, Grid) \
#define INSTANCE_GRID(Dim, Grid) \ template void BdaSolverInfo<BM<T,Dim>,BV<T,Dim>>:: \
template void BdaSolverInfo<BM<Dim>,BV<Dim>>:: \ prepare(const Grid&, \
prepare(const Grid&, \ const Dune::CartesianIndexMapper<Grid>&, \
const Dune::CartesianIndexMapper<Grid>&, \ const std::vector<Well>&, \
const std::vector<Well>&, \ const std::vector<int>&, \
const std::vector<int>&, \
const std::size_t, const bool); const std::size_t, const bool);
using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>; using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>;
#if HAVE_DUNE_ALUGRID #if HAVE_DUNE_ALUGRID
@ -248,23 +249,26 @@ using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>;
#else #else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI #endif //HAVE_MPI
#define INSTANCE(Dim) \ #define INSTANTIATE(T,Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \ template struct BdaSolverInfo<BM<T,Dim>,BV<T,Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \ INSTANTIATE_GRID(T,Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,ALUGrid3CN) \ INSTANTIATE_GRID(T,Dim,ALUGrid3CN) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) INSTANTIATE_GRID(T,Dim,PolyHedralGrid3D)
#else #else
#define INSTANCE(Dim) \ #define INSTANTIATE(T,Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \ template struct BdaSolverInfo<BM<T,Dim>,BV<T,Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \ INSTANTIATE_GRID(T,Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) INSTANTIATE_GRID(T,Dim,PolyHedralGrid3D)
#endif #endif
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)
INSTANCE(4)
INSTANCE(5)
INSTANCE(6)
} // namespace detail #define INSTANTIATE_TYPE(T) \
} // namespace Opm INSTANTIATE(T,1) \
INSTANTIATE(T,2) \
INSTANTIATE(T,3) \
INSTANTIATE(T,4) \
INSTANTIATE(T,5) \
INSTANTIATE(T,6)
INSTANTIATE_TYPE(double)
} // namespace Opm::detail

9
opm/simulators/linalg/ISTLSolverBda.hpp
View File

@ -41,13 +41,14 @@ namespace detail {
template<class Matrix, class Vector> template<class Matrix, class Vector>
struct BdaSolverInfo struct BdaSolverInfo
{ {
using WellContribFunc = std::function<void(WellContributions<double>&)>; using Scalar = typename Vector::field_type;
using WellContribFunc = std::function<void(WellContributions<Scalar>&)>;
using Bridge = BdaBridge<Matrix,Vector,Matrix::block_type::rows>; using Bridge = BdaBridge<Matrix,Vector,Matrix::block_type::rows>;
BdaSolverInfo(const std::string& accelerator_mode, BdaSolverInfo(const std::string& accelerator_mode,
const int linear_solver_verbosity, const int linear_solver_verbosity,
const int maxit, const int maxit,
const double tolerance, const Scalar tolerance,
const int platformID, const int platformID,
const int deviceID, const int deviceID,
const bool opencl_ilu_parallel, const bool opencl_ilu_parallel,
@ -249,8 +250,8 @@ public:
// Solve system. // Solve system.
Dune::InverseOperatorResult result; Dune::InverseOperatorResult result;
std::function<void(WellContributions<double>&)> getContribs = std::function<void(WellContributions<Scalar>&)> getContribs =
[this](WellContributions<double>& w) [this](WellContributions<Scalar>& w)
{ {
this->simulator_.problem().wellModel().getWellContributions(w); this->simulator_.problem().wellModel().getWellContributions(w);
}; };