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

MSWellHelpers: simplify interfaces

Browse Source 
- avoid passing the matrix, only pass the solver. possible
  since setting up the solver is not done in here any more.
- avoid passing shared_ptr's
This commit is contained in:
Arne Morten Kvarving
2022-12-01 11:14:55 +01:00
parent e67e58d0c8
commit a7cb444328
5 changed files with 26 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
opm/simulators/wells/MSWellHelpers.cpp
View File

@@ -100,16 +100,10 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename MatrixType, typename VectorType>
VectorType
applyUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType>>& linsolver,
applyUMFPack(Dune::UMFPack<MatrixType>& linsolver,
VectorType x)
{
#if HAVE_UMFPACK
if (!linsolver)
{
linsolver = std::make_shared<Dune::UMFPack<MatrixType>>(D, 0);
}
// The copy of x seems mandatory for calling UMFPack!
VectorType y(x.size());
y = 0.;
@@ -118,7 +112,7 @@ applyUMFPack(const MatrixType& D,
Dune::InverseOperatorResult res;
// Solve
linsolver->apply(y, x, res);
linsolver.apply(y, x, res);
// Checking if there is any inf or nan in y
// it will be the solution before we find a way to catch the singularity of the matrix
@@ -141,24 +135,24 @@ applyUMFPack(const MatrixType& D,
template <typename VectorType, typename MatrixType>
Dune::Matrix<typename MatrixType::block_type>
invertWithUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver)
invertWithUMFPack(const int size,
const int bsize,
Dune::UMFPack<MatrixType>& linsolver)
{
#if HAVE_UMFPACK
const int sz = D.M();
const int bsz = D[0][0].M();
VectorType e(sz);
VectorType e(size);
e = 0.0;
// Make a full block matrix.
Dune::Matrix<typename MatrixType::block_type> inv(sz, sz);
Dune::Matrix<typename MatrixType::block_type> inv(size, size);
// Create inverse by passing basis vectors to the solver.
for (int ii = 0; ii < sz; ++ii) {
for (int jj = 0; jj < bsz; ++jj) {
for (int ii = 0; ii < size; ++ii) {
for (int jj = 0; jj < bsize; ++jj) {
e[ii][jj] = 1.0;
auto col = applyUMFPack(D, linsolver, e);
for (int cc = 0; cc < sz; ++cc) {
for (int dd = 0; dd < bsz; ++dd) {
auto col = applyUMFPack(linsolver, e);
for (int cc = 0; cc < size; ++cc) {
for (int dd = 0; dd < bsize; ++dd) {
inv[cc][ii][dd][jj] = col[cc][dd];
}
}
@@ -328,12 +322,10 @@ template<int Dim>
using Mat = Dune::BCRSMatrix<Dune::FieldMatrix<double,Dim,Dim>>;
#define INSTANCE_UMF(Dim) \
template Vec<Dim> applyUMFPack<Mat<Dim>,Vec<Dim>>(const Mat<Dim>&, \
std::shared_ptr<Dune::UMFPack<Mat<Dim>>>&, \
template Vec<Dim> applyUMFPack<Mat<Dim>,Vec<Dim>>(Dune::UMFPack<Mat<Dim>>&, \
Vec<Dim>); \
template Dune::Matrix<typename Mat<Dim>::block_type> \
invertWithUMFPack<Vec<Dim>,Mat<Dim>>(const Mat<Dim>& D, \
std::shared_ptr<Dune::UMFPack<Mat<Dim>>>&);
invertWithUMFPack<Vec<Dim>,Mat<Dim>>(const int, const int, Dune::UMFPack<Mat<Dim>>&);
INSTANCE_UMF(2)
INSTANCE_UMF(3)

8
opm/simulators/wells/MSWellHelpers.hpp
View File

@@ -42,8 +42,7 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename MatrixType, typename VectorType>
VectorType
applyUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType>>& linsolver,
applyUMFPack(Dune::UMFPack<MatrixType>& linsolver,
VectorType x);
@@ -51,8 +50,9 @@ namespace mswellhelpers
/// Applies umfpack and checks for singularity
template <typename VectorType, typename MatrixType>
Dune::Matrix<typename MatrixType::block_type>
invertWithUMFPack(const MatrixType& D,
std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver);
invertWithUMFPack(const int size,
const int bsize,
Dune::UMFPack<MatrixType>& linsolver);

4
opm/simulators/wells/MultisegmentWellEquations.cpp
View File

@@ -130,7 +130,7 @@ apply(const BVector& x, BVector& Ax) const
duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx);
const BVectorWell invDBx = mswellhelpers::applyUMFPack(*duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
duneC_.mmtv(invDBx,Ax);
@@ -141,7 +141,7 @@ void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
apply(BVector& r) const
{
// invDrw_ = duneD^-1 * resWell_
const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);
const BVectorWell invDrw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell_);
// r = r - duneC_^T * invDrw
duneC_.mmtv(invDrw, r);
}

2
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@@ -400,7 +400,7 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const
// resWell = resWell - B * x
linSys_.duneB_.mmv(x, resWell);
// xw = D^-1 * resWell
xw = mswellhelpers::applyUMFPack(linSys_.duneD_, linSys_.duneDSolver_, resWell);
xw = mswellhelpers::applyUMFPack(*linSys_.duneDSolver_, resWell);
}
template<typename FluidSystem, typename Indices, typename Scalar>

10
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@@ -519,8 +519,7 @@ namespace Opm
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_,
const BVectorWell dx_well = mswellhelpers::applyUMFPack(*this->linSys_.duneDSolver_,
this->linSys_.resWell_);
updateWellState(dx_well, well_state, deferred_logger);
@@ -722,7 +721,9 @@ namespace Opm
MultisegmentWell<TypeTag>::
addWellContributions(SparseMatrixAdapter& jacobian) const
{
const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(this->linSys_.duneD_, this->linSys_.duneDSolver_);
const auto invDuneD = mswellhelpers::invertWithUMFPack<BVectorWell>(numWellEq,
Indices::numEq,
*this->linSys_.duneDSolver_);
// We need to change matrix A as follows
// A -= C^T D^-1 B
@@ -1494,8 +1495,7 @@ namespace Opm
assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_,
const BVectorWell dx_well = mswellhelpers::applyUMFPack(*this->linSys_.duneDSolver_,
this->linSys_.resWell_);
if (it > this->param_.strict_inner_iter_wells_) {