diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index 226ca51e0..22d66e159 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -72,7 +72,50 @@ namespace mswellhelpers #else // this is not thread safe OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. " - "Reconfigure opm-simulator with SuiteSparse/UMFPACK support and recompile."); + "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); +#endif // HAVE_UMFPACK + } + + + + /// Applies umfpack and checks for singularity + template + MatrixType + invertWithUMFPack(const MatrixType& D, std::shared_ptr >& linsolver) + { +#if HAVE_UMFPACK + const int sz = D.M(); + const int bsz = D[0][0].M(); + VectorType e(sz); + e = 0.0; + + // Make a block matrix with a full pattern. + MatrixType inv(sz, sz, sz*sz, MatrixType::row_wise); + for (auto row = inv.createbegin(); row != inv.createend(); ++row) { + for (int c = 0; c < sz; ++c) { + row.insert(c); + } + } + + // Create inverse by passing basis vectors to the solver. + for (int ii = 0; ii < sz; ++ii) { + for (int jj = 0; jj < bsz; ++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) { + inv[cc][ii][dd][jj] = col[cc][dd]; + } + } + e[ii][jj] = 0.0; + } + } + + return inv; +#else + // this is not thread safe + OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. " + "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); #endif // HAVE_UMFPACK }