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

- moded all bda spesific tings to separete class

Browse Source
This commit is contained in:
hnil 2023-08-09 10:22:33 +02:00
parent 659d5efa3e
commit d623695d2a
2 changed files with 3 additions and 324 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

197
opm/simulators/linalg/ISTLSolverEbos.cpp
View File

@ -202,169 +202,6 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
}
}
#if COMPILE_BDA_BRIDGE
template<class Matrix, class Vector>
BdaSolverInfo<Matrix,Vector>::
BdaSolverInfo(const std::string& accelerator_mode,
const int linear_solver_verbosity,
const int maxit,
const double tolerance,
const int platformID,
const int deviceID,
const bool opencl_ilu_parallel,
const std::string& linsolver)
: bridge_(std::make_unique<Bridge>(accelerator_mode,
linear_solver_verbosity, maxit,
tolerance, platformID, deviceID,
opencl_ilu_parallel, linsolver))
, accelerator_mode_(accelerator_mode)
{}
template<class Matrix, class Vector>
BdaSolverInfo<Matrix,Vector>::~BdaSolverInfo() = default;
template<class Matrix, class Vector>
template<class Grid>
void BdaSolverInfo<Matrix,Vector>::
prepare(const Grid& grid,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const std::vector<Well>& wellsForConn,
const std::vector<int>& cellPartition,
const size_t nonzeroes,
const bool useWellConn)
{
if (numJacobiBlocks_ > 1) {
detail::setWellConnections(grid, cartMapper, wellsForConn,
useWellConn,
wellConnectionsGraph_,
numJacobiBlocks_);
this->blockJacobiAdjacency(grid, cellPartition, nonzeroes);
}
}
template<class Matrix, class Vector>
bool BdaSolverInfo<Matrix,Vector>::
apply(Vector& rhs,
const bool useWellConn,
WellContribFunc getContribs,
const int rank,
Matrix& matrix,
Vector& x,
Dune::InverseOperatorResult& result)
{
bool use_gpu = bridge_->getUseGpu();
if (use_gpu) {
auto wellContribs = WellContributions::create(accelerator_mode_, useWellConn);
bridge_->initWellContributions(*wellContribs, x.N() * x[0].N());
// the WellContributions can only be applied separately with CUDA or OpenCL, not with amgcl or rocalution
#if HAVE_CUDA || HAVE_OPENCL
if (!useWellConn) {
getContribs(*wellContribs);
}
#endif
if (numJacobiBlocks_ > 1) {
this->copyMatToBlockJac(matrix, *blockJacobiForGPUILU0_);
// Const_cast needed since the CUDA stuff overwrites values for better matrix condition..
bridge_->solve_system(&matrix, blockJacobiForGPUILU0_.get(),
numJacobiBlocks_, rhs, *wellContribs, result);
}
else
bridge_->solve_system(&matrix, &matrix,
numJacobiBlocks_, rhs, *wellContribs, result);
if (result.converged) {
// get result vector x from non-Dune backend, iff solve was successful
bridge_->get_result(x);
return true;
} else {
// warn about CPU fallback
// BdaBridge might have disabled its BdaSolver for this simulation due to some error
// in that case the BdaBridge is disabled and flexibleSolver is always used
// or maybe the BdaSolver did not converge in time, then it will be used next linear solve
if (rank == 0) {
OpmLog::warning(bridge_->getAccleratorName() + " did not converge, now trying Dune to solve current linear system...");
}
}
}
return false;
}
template<class Matrix, class Vector>
bool BdaSolverInfo<Matrix,Vector>::
gpuActive()
{
return bridge_->getUseGpu();
}
template<class Matrix, class Vector>
template<class Grid>
void BdaSolverInfo<Matrix,Vector>::
blockJacobiAdjacency(const Grid& grid,
const std::vector<int>& cell_part,
size_t nonzeroes)
{
using size_type = typename Matrix::size_type;
using Iter = typename Matrix::CreateIterator;
size_type numCells = grid.size(0);
blockJacobiForGPUILU0_ = std::make_unique<Matrix>(numCells, numCells,
nonzeroes, Matrix::row_wise);
const auto& lid = grid.localIdSet();
const auto& gridView = grid.leafGridView();
auto elemIt = gridView.template begin<0>(); // should never overrun, since blockJacobiForGPUILU0_ is initialized with numCells rows
//Loop over cells
for (Iter row = blockJacobiForGPUILU0_->createbegin(); row != blockJacobiForGPUILU0_->createend(); ++elemIt, ++row)
{
const auto& elem = *elemIt;
size_type idx = lid.id(elem);
row.insert(idx);
// Add well non-zero connections
for (const auto wc : wellConnectionsGraph_[idx]) {
row.insert(wc);
}
int locPart = cell_part[idx];
//Add neighbor if it is on the same part
auto isend = gridView.iend(elem);
for (auto is = gridView.ibegin(elem); is!=isend; ++is)
{
//check if face has neighbor
if (is->neighbor())
{
size_type nid = lid.id(is->outside());
int nabPart = cell_part[nid];
if (locPart == nabPart) {
row.insert(nid);
}
}
}
}
}
template<class Matrix, class Vector>
void BdaSolverInfo<Matrix,Vector>::
copyMatToBlockJac(const Matrix& mat, Matrix& blockJac)
{
auto rbegin = blockJac.begin();
auto rend = blockJac.end();
auto outerRow = mat.begin();
for (auto row = rbegin; row != rend; ++row, ++outerRow) {
auto outerCol = (*outerRow).begin();
for (auto col = (*row).begin(); col != (*row).end(); ++col) {
// outerRow is guaranteed to have all column entries that row has!
while(outerCol.index() < col.index()) ++outerCol;
assert(outerCol.index() == col.index());
*col = *outerCol; // copy nonzero block
}
}
}
#endif // COMPILE_BDA_BRIDGE
template<int Dim>
using BM = Dune::BCRSMatrix<MatrixBlock<double,Dim,Dim>>;
template<int Dim>
@ -380,40 +217,6 @@ using CommunicationType = Dune::CollectiveCommunication<int>;
template void makeOverlapRowsInvalid<BM<Dim>>(BM<Dim>&, const std::vector<int>&); \
template struct FlexibleSolverInfo<BM<Dim>,BV<Dim>,CommunicationType>;
#if COMPILE_BDA_BRIDGE
#define INSTANCE_GRID(Dim, Grid) \
template void BdaSolverInfo<BM<Dim>,BV<Dim>>:: \
prepare(const Grid&, \
const Dune::CartesianIndexMapper<Grid>&, \
const std::vector<Well>&, \
const std::vector<int>&, \
const size_t, const bool);
using PolyHedralGrid3D = Dune::PolyhedralGrid<3, 3>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
#define INSTANCE(Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,ALUGrid3CN) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) \
INSTANCE_FLEX(Dim)
#else
#define INSTANCE(Dim) \
template struct BdaSolverInfo<BM<Dim>,BV<Dim>>; \
INSTANCE_GRID(Dim,Dune::CpGrid) \
INSTANCE_GRID(Dim,PolyHedralGrid3D) \
INSTANCE_FLEX(Dim)
#endif
#else
#define INSTANCE(Dim) \
INSTANCE_FLEX(Dim)
#endif // COMPILE_BDA_BRIDGE
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)

130
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -87,10 +87,6 @@ public:
namespace Opm
{
#if COMPILE_BDA_BRIDGE
template<class Matrix, class Vector, int block_size> class BdaBridge;
class WellContributions;
#endif
namespace detail
{
@ -116,60 +112,6 @@ struct FlexibleSolverInfo
size_t interiorCellNum_ = 0;
};
#if COMPILE_BDA_BRIDGE
template<class Matrix, class Vector>
struct BdaSolverInfo
{
using WellContribFunc = std::function<void(WellContributions&)>;
using Bridge = BdaBridge<Matrix,Vector,Matrix::block_type::rows>;
BdaSolverInfo(const std::string& accelerator_mode,
const int linear_solver_verbosity,
const int maxit,
const double tolerance,
const int platformID,
const int deviceID,
const bool opencl_ilu_parallel,
const std::string& linsolver);
~BdaSolverInfo();
template<class Grid>
void prepare(const Grid& grid,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const std::vector<Well>& wellsForConn,
const std::vector<int>& cellPartition,
const size_t nonzeroes,
const bool useWellConn);
bool apply(Vector& rhs,
const bool useWellConn,
WellContribFunc getContribs,
const int rank,
Matrix& matrix,
Vector& x,
Dune::InverseOperatorResult& result);
bool gpuActive();
int numJacobiBlocks_ = 0;
private:
/// Create sparsity pattern for block-Jacobi matrix based on partitioning of grid.
/// Do not initialize the values, that is done in copyMatToBlockJac()
template<class Grid>
void blockJacobiAdjacency(const Grid& grid,
const std::vector<int>& cell_part,
size_t nonzeroes);
void copyMatToBlockJac(const Matrix& mat, Matrix& blockJac);
std::unique_ptr<Bridge> bridge_;
std::string accelerator_mode_;
std::unique_ptr<Matrix> blockJacobiForGPUILU0_;
std::vector<std::set<int>> wellConnectionsGraph_;
};
#endif
#ifdef HAVE_MPI
/// Copy values in parallel.
@ -270,37 +212,6 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
#if HAVE_MPI
comm_.reset( new CommunicationType( simulator_.vanguard().grid().comm() ) );
#endif
#if COMPILE_BDA_BRIDGE
{
std::string accelerator_mode = EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode);
if ((simulator_.vanguard().grid().comm().size() > 1) && (accelerator_mode != "none")) {
if (on_io_rank) {
OpmLog::warning("Cannot use GPU with MPI, GPU are disabled");
}
accelerator_mode = "none";
}
const int platformID = EWOMS_GET_PARAM(TypeTag, int, OpenclPlatformId);
const int deviceID = EWOMS_GET_PARAM(TypeTag, int, BdaDeviceId);
const int maxit = EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIter);
const double tolerance = EWOMS_GET_PARAM(TypeTag, double, LinearSolverReduction);
const bool opencl_ilu_parallel = EWOMS_GET_PARAM(TypeTag, bool, OpenclIluParallel);
const int linear_solver_verbosity = parameters_.linear_solver_verbosity_;
std::string linsolver = EWOMS_GET_PARAM(TypeTag, std::string, LinearSolver);
bdaBridge = std::make_unique<detail::BdaSolverInfo<Matrix,Vector>>(accelerator_mode,
linear_solver_verbosity,
maxit,
tolerance,
platformID,
deviceID,
opencl_ilu_parallel,
linsolver);
}
#else
if (EWOMS_GET_PARAM(TypeTag, std::string, AcceleratorMode) != "none") {
OPM_THROW(std::logic_error,"Cannot use accelerated solver since CUDA, OpenCL and amgcl were not found by cmake");
}
#endif
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
// For some reason simulator_.model().elementMapper() is not initialized at this stage
@ -341,7 +252,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
void prepare(const Matrix& M, Vector& b)
{
OPM_TIMEBLOCK(istlSolverEbosPrepare);
OPM_TIMEBLOCK(istlSolverEbosPrepare);
const bool firstcall = (matrix_ == nullptr);
#if HAVE_MPI
if (firstcall) {
@ -359,14 +270,6 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
// setup sparsity pattern for jacobi matrix for preconditioner (only used for openclSolver)
#if HAVE_OPENCL
bdaBridge->numJacobiBlocks_ = EWOMS_GET_PARAM(TypeTag, int, NumJacobiBlocks);
bdaBridge->prepare(simulator_.vanguard().grid(),
simulator_.vanguard().cartesianIndexMapper(),
simulator_.vanguard().schedule().getWellsatEnd(),
simulator_.vanguard().cellPartition(),
getMatrix().nonzeroes(), useWellConn_);
#endif
} else {
// Pointers should not change
if ( &M != matrix_ ) {
@ -379,13 +282,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
if (isParallel() && prm_.get<std::string>("preconditioner.type") != "ParOverILU0") {
detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
}
#if COMPILE_BDA_BRIDGE
if(!bdaBridge->gpuActive()){
prepareFlexibleSolver();
}
#else
prepareFlexibleSolver();
#endif
}
@ -406,7 +303,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
bool solve(Vector& x)
{
OPM_TIMEBLOCK(istlSolverEbosSolve);
OPM_TIMEBLOCK(istlSolverEbosSolve);
calls_ += 1;
// Write linear system if asked for.
const int verbosity = prm_.get<int>("verbosity", 0);
@ -420,25 +317,8 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
// Solve system.
Dune::InverseOperatorResult result;
#if COMPILE_BDA_BRIDGE
std::function<void(WellContributions&)> getContribs =
[this](WellContributions& w)
{
this->simulator_.problem().wellModel().getWellContributions(w);
};
if (!bdaBridge->apply(*rhs_, useWellConn_, getContribs,
simulator_.gridView().comm().rank(),
const_cast<Matrix&>(this->getMatrix()),
x, result))
#endif
{
OPM_TIMEBLOCK(flexibleSolverApply);
#if COMPILE_BDA_BRIDGE
if(bdaBridge->gpuActive()){
prepareFlexibleSolver();
}
#endif
assert(flexibleSolver_.solver_);
flexibleSolver_.solver_->apply(x, *rhs_, result);
}
@ -521,7 +401,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
}
else
{
OPM_TIMEBLOCK(flexibleSolverUpdate);
OPM_TIMEBLOCK(flexibleSolverUpdate);
flexibleSolver_.pre_->update();
}
}
@ -609,10 +489,6 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
Matrix* matrix_;
Vector *rhs_;
#if COMPILE_BDA_BRIDGE
std::unique_ptr<detail::BdaSolverInfo<Matrix, Vector>> bdaBridge;
#endif
detail::FlexibleSolverInfo<Matrix,Vector,CommunicationType> flexibleSolver_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;