OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 18:21:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

Bugfix: ensure the matrix passed to parallel preconditioner is ok.

Browse Source 
Done with the makeOverlapRowsInvalid() approach.
This commit is contained in:
Atgeirr Flø Rasmussen 2019-06-05 14:04:41 +02:00
parent 866db1c66d
commit d44b974dc8
1 changed files with 32 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
opm/simulators/linalg/ISTLSolverEbosFlexible.hpp
View File

@ -21,6 +21,7 @@
#define OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED #define OPM_ISTLSOLVEREBOSFLEXIBLE_HEADER_INCLUDED
#include <ewoms/linear/matrixblock.hh> #include <ewoms/linear/matrixblock.hh>
#include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp> #include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/simulators/linalg/setupPropertyTree.hpp> #include <opm/simulators/linalg/setupPropertyTree.hpp>
@ -57,6 +58,7 @@ class ISTLSolverEbosFlexible
using SparseMatrixAdapter = typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter); using SparseMatrixAdapter = typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter);
using VectorType = typename GET_PROP_TYPE(TypeTag, GlobalEqVector); using VectorType = typename GET_PROP_TYPE(TypeTag, GlobalEqVector);
using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator); using Simulator = typename GET_PROP_TYPE(TypeTag, Simulator);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using MatrixType = typename SparseMatrixAdapter::IstlMatrix; using MatrixType = typename SparseMatrixAdapter::IstlMatrix;
#if HAVE_MPI #if HAVE_MPI
using Communication = Dune::OwnerOverlapCopyCommunication<int, int>; using Communication = Dune::OwnerOverlapCopyCommunication<int, int>;
@ -77,6 +79,7 @@ public:
parameters_.template init<TypeTag>(); parameters_.template init<TypeTag>();
prm_ = setupPropertyTree(parameters_); prm_ = setupPropertyTree(parameters_);
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
detail::findOverlapRowsAndColumns(simulator_.vanguard().grid(), overlapRowAndColumns_);
#if HAVE_MPI #if HAVE_MPI
if (parallelInformation_.type() == typeid(ParallelISTLInformation)) { if (parallelInformation_.type() == typeid(ParallelISTLInformation)) {
// Parallel case. // Parallel case.
@ -91,7 +94,7 @@ public:
{ {
} }
void prepare(const SparseMatrixAdapter& mat, VectorType& b) void prepare(SparseMatrixAdapter& mat, VectorType& b)
{ {
#if HAVE_MPI #if HAVE_MPI
static bool firstcall = true; static bool firstcall = true;
@ -103,6 +106,7 @@ public:
parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1); parinfo->copyValuesTo(comm_->indexSet(), comm_->remoteIndices(), size, 1);
firstcall = false; firstcall = false;
} }
makeOverlapRowsInvalid(mat.istlMatrix());
#endif #endif
// Decide if we should recreate the solver or just do // Decide if we should recreate the solver or just do
// a minimal preconditioner update. // a minimal preconditioner update.
@ -167,6 +171,32 @@ public:
} }
protected: protected:
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1e100).
void makeOverlapRowsInvalid(MatrixType& matrix) const
{
// Value to set on diagonal
const int numEq = MatrixType::block_type::rows;
typename MatrixType::block_type diag_block(0.0);
for (int eq = 0; eq < numEq; ++eq)
diag_block[eq][eq] = 1.0e100;
// loop over precalculated overlap rows and columns
for (auto row = overlapRowAndColumns_.begin(); row != overlapRowAndColumns_.end(); row++) {
int lcell = row->first;
// diagonal block set to large value diagonal
matrix[lcell][lcell] = diag_block;
// loop over off diagonal blocks in overlap row
for (auto col = row->second.begin(); col != row->second.end(); ++col) {
int ncell = *col;
// zero out block
matrix[lcell][ncell] = 0.0;
}
}
}
const Simulator& simulator_; const Simulator& simulator_;
std::unique_ptr<SolverType> solver_; std::unique_ptr<SolverType> solver_;
@ -176,6 +206,7 @@ protected:
Dune::InverseOperatorResult res_; Dune::InverseOperatorResult res_;
boost::any parallelInformation_; boost::any parallelInformation_;
std::unique_ptr<Communication> comm_; std::unique_ptr<Communication> comm_;
std::vector<std::pair<int, std::vector<int>>> overlapRowAndColumns_;
}; // end ISTLSolverEbosFlexible }; // end ISTLSolverEbosFlexible
} // namespace Opm } // namespace Opm