38#ifndef OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
39#define OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
41#include <opm/common/ErrorMacros.hpp>
42#include <opm/grid/utility/SparseTable.hpp>
43#include <opm/porsol/common/BoundaryConditions.hpp>
44#include <opm/porsol/common/Matrix.hpp>
46#include <opm/common/utility/platform_dependent/disable_warnings.h>
48#include <unordered_map>
50#include <dune/common/fvector.hh>
51#include <dune/common/fmatrix.hh>
53#include <dune/istl/bvector.hh>
54#include <dune/istl/bcrsmatrix.hh>
55#include <dune/istl/operators.hh>
56#include <dune/istl/io.hh>
57#include <dune/istl/overlappingschwarz.hh>
58#include <dune/istl/schwarz.hh>
59#include <dune/istl/preconditioners.hh>
60#include <dune/istl/solvers.hh>
61#include <dune/istl/owneroverlapcopy.hh>
62#include <dune/istl/paamg/amg.hh>
63#include <dune/common/version.hh>
64#include <dune/istl/paamg/fastamg.hh>
65#include <dune/istl/paamg/kamg.hh>
66#include <dune/istl/paamg/pinfo.hh>
68#include <opm/common/utility/platform_dependent/reenable_warnings.h>
108 bool topologyIsSane(
const GI& g)
110 typedef typename GI::CellIterator CI;
111 typedef typename CI::FaceIterator FI;
113 bool sane = g.numberOfCells() >= 0;
115 for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
118 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
119 if (!f->boundary()) {
120 n.push_back(f->neighbourCellIndex());
123 std::sort(n.begin(), n.end());
125 sane = std::unique(n.begin(), n.end()) == n.end();
158 class axpby :
public std::binary_function<T,T,T> {
165 axpby(
const T& a,
const T& b) : a_(a), b_(b) {}
179 T operator()(
const T& x,
const T& y)
361 template <
class GridInterface,
364 template <
class Gr
idIF,
class RockIF>
class InnerProduct>
371 typedef typename GridInterface::Scalar Scalar;
377 enum FaceType { Internal, Dirichlet, Neumann, Periodic };
388 typedef typename GridInterface::Scalar Scalar;
393 typedef typename GridInterface::CellIterator CI;
397 typedef typename CI ::FaceIterator FI;
410 Scalar pressure(
const CI& c)
const
412 return pressure_[cellno_[c->index()]];
425 Scalar outflux (
const FI& f)
const
427 return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
429 Scalar outflux (
int hf)
const
431 return outflux_.data(hf);
434 std::vector< int > cellno_;
435 Opm::SparseTable< int > cellFaces_;
436 std::vector<Scalar> pressure_;
437 Opm::SparseTable<Scalar> outflux_;
440 std::vector<int>().swap(cellno_);
443 std::vector<Scalar>().swap(pressure_);
475 template<
class Po
int>
476 void init(
const GridInterface& g,
477 const RockInterface& r,
479 const BCInterface& bc)
483 if (g.numberOfCells() > 0) {
501 num_internal_faces_ = 0;
502 total_num_faces_ = 0;
503 matrix_structure_valid_ =
false;
504 do_regularization_ =
true;
506 bdry_id_map_.clear();
508 std::vector<Scalar>().swap(L_);
509 std::vector<Scalar>().swap(g_);
512 flowSolution_.clear();
514 cleared_state_ =
true;
534 const BCInterface& bc)
537 assert (cleared_state_);
539 assert (topologyIsSane(g));
542 allocateConnections(bc);
563 template<
class Po
int>
568 assert (matrix_structure_valid_);
570 typedef typename GridInterface::CellIterator CI;
571 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
574 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
575 ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
656 template<
class Flu
idInterface>
657 void solve(
const FluidInterface& r ,
658 const std::vector<double>& sat,
659 const BCInterface& bc ,
660 const std::vector<double>& src,
661 double residual_tolerance = 1e-8,
662 int linsolver_verbosity = 1,
663 int linsolver_type = 1,
664 bool same_matrix =
false,
665 int linsolver_maxit = 0,
666 double prolongate_factor = 1.6,
667 int smooth_steps = 1)
669 assembleDynamic(r, sat, bc, src);
673 switch (linsolver_type) {
675 solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
678 solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
679 linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
683 solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
684 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
687 solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
688 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
691 std::cerr <<
"Unknown linsolver_type: " << linsolver_type <<
'\n';
692 throw std::runtime_error(
"Unknown linsolver_type");
694 computePressureAndFluxes(r, sat);
703 : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
706 void put(
double flux,
int f_ix) {
707 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
708 double sign = visited_[f_ix] ? -1.0 : 1.0;
709 fluxes_[f_ix] += sign*flux;
712 void get(
double& flux,
int f_ix) {
713 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
714 double sign = visited_[f_ix] ? -1.0 : 1.0;
715 double new_flux = 0.5*sign*fluxes_[f_ix];
716 double diff = std::fabs(flux - new_flux);
717 max_modification_ = std::max(max_modification_, diff);
723 std::fill(visited_.begin(), visited_.end(), 0);
726 double maxMod()
const
728 return max_modification_;
731 std::vector<double> fluxes_;
732 std::vector<int> visited_;
733 double max_modification_;
748 typedef typename GridInterface::CellIterator CI;
749 typedef typename CI ::FaceIterator FI;
750 const std::vector<int>& cell = flowSolution_.cellno_;
751 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
752 Opm::SparseTable<double>& cflux = flowSolution_.outflux_;
754 FaceFluxes face_fluxes(pgrid_->numberOfFaces());
756 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
757 const int cell_index = cell[c->index()];
758 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
759 int f_ix = cf[cell_index][f->localIndex()];
760 double flux = cflux[cell_index][f->localIndex()];
762 if (ppartner_dof_.empty()) {
765 int partner_f_ix = ppartner_dof_[f_ix];
766 if (partner_f_ix != -1) {
767 face_fluxes.put(flux, f_ix);
768 face_fluxes.put(flux, partner_f_ix);
771 face_fluxes.put(flux, f_ix);
775 face_fluxes.resetVisited();
777 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
778 const int cell_index = cell[c->index()];
779 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
780 int f_ix = cf[cell_index][f->localIndex()];
781 double& flux = cflux[cell_index][f->localIndex()];
783 if (ppartner_dof_.empty()) {
786 int partner_f_ix = ppartner_dof_[f_ix];
787 if (partner_f_ix != -1) {
788 face_fluxes.get(flux, f_ix);
790 face_fluxes.get(dummy, partner_f_ix);
791 assert(dummy == flux);
794 face_fluxes.get(flux, f_ix);
798 return face_fluxes.maxMod();
822 return flowSolution_;
840 template<
typename charT,
class traits>
843 os <<
"IncompFlowSolverHybrid<>:\n"
844 <<
"\tMaximum number of cell faces = " << max_ncf_ <<
'\n'
845 <<
"\tNumber of internal faces = " << num_internal_faces_ <<
'\n'
846 <<
"\tTotal number of faces = " << total_num_faces_ <<
'\n';
848 const std::vector<int>& cell = flowSolution_.cellno_;
849 os <<
"cell index map = [";
850 std::copy(cell.begin(), cell.end(),
851 std::ostream_iterator<int>(os,
" "));
854 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
855 os <<
"cell faces =\n";
856 for (
int i = 0; i < cf.size(); ++i)
858 os <<
"\t[" << i <<
"] -> [";
859 std::copy(cf[i].begin(), cf[i].end(),
860 std::ostream_iterator<int>(os,
","));
889 writeMatrixToMatlab(S_, prefix +
"-mat.dat");
891 std::string rhsfile(prefix +
"-rhs.dat");
892 std::ofstream rhs(rhsfile.c_str());
894 rhs.setf(std::ios::scientific | std::ios::showpos);
895 std::copy(rhs_.begin(), rhs_.end(),
896 std::ostream_iterator<VectorBlockType>(rhs,
"\n"));
900 typedef std::pair<int,int> DofID;
901 typedef std::unordered_map<int,DofID> BdryIdMapType;
902 typedef BdryIdMapType::const_iterator BdryIdMapIterator;
904 const GridInterface* pgrid_;
905 BdryIdMapType bdry_id_map_;
906 std::vector<int> ppartner_dof_;
908 InnerProduct<GridInterface, RockInterface> ip_;
913 int num_internal_faces_;
914 int total_num_faces_;
917 std::vector<Scalar> L_, g_;
918 Opm::SparseTable<Scalar> F_ ;
922 typedef Dune::FieldVector<Scalar, 1 > VectorBlockType;
923 typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
925 Dune::BCRSMatrix <MatrixBlockType> S_;
926 Dune::BlockVector<VectorBlockType> rhs_;
927 Dune::BlockVector<VectorBlockType> soln_;
928 bool matrix_structure_valid_;
929 bool do_regularization_;
933 FlowSolution flowSolution_;
937 void enumerateDof(
const GridInterface& g,
const BCInterface& bc)
941 enumerateBCDof(g, bc);
944 cleared_state_ =
false;
948 void enumerateGridDof(
const GridInterface& g)
951 typedef typename GridInterface::CellIterator CI;
952 typedef typename CI ::FaceIterator FI;
954 const int nc = g.numberOfCells();
955 std::vector<int> fpos ; fpos.reserve(nc + 1);
956 std::vector<int> num_cf(nc) ;
957 std::vector<int> faces ;
960 std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
962 std::vector<int>& cell = flowSolution_.cellno_;
965 int cellno = 0; fpos.push_back(0);
967 for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++cellno) {
968 const int c0 = c->index();
969 assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
973 int& ncf = num_cf[c0];
975 for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
976 if (!f->boundary()) {
977 const int c1 = f->neighbourCellIndex();
978 assert((0 <= c1) && (c1 < nc) && (c1 != c0));
980 if (cell[c1] == -1) {
988 fpos.push_back(
int(faces.size()));
989 max_ncf_ = std::max(max_ncf_, ncf);
992 assert (cellno == nc);
994 total_num_faces_ = num_internal_faces_ = int(faces.size());
996 ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
997 F_ .reserve(nc, tot_ncf);
999 flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1000 flowSolution_.outflux_ .reserve(nc, tot_ncf);
1002 Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1005 std::vector<int> l2g; l2g .reserve(max_ncf_);
1006 std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1009 typedef std::vector<int>::iterator VII;
1010 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1011 const int c0 = c->index();
1013 assert ((0 <= c0 ) && ( c0 < nc) &&
1014 (0 <= cell[c0]) && (cell[c0] < nc));
1016 const int ncf = num_cf[cell[c0]];
1017 l2g .resize(ncf , 0 );
1018 F_alloc .resize(ncf , Scalar(0.0));
1020 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1021 if (f->boundary()) {
1023 l2g[f->localIndex()] = total_num_faces_++;
1035 const int c1 = f->neighbourCellIndex();
1036 assert ((0 <= c1 ) && ( c1 < nc) &&
1037 (0 <= cell[c1]) && (cell[c1] < nc));
1039 int t = c0, seek = c1;
1040 if (cell[seek] < cell[t])
1043 int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1045 VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1046 assert(p != faces.begin() + e);
1048 l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1052 cf.appendRow (l2g .begin(), l2g .end());
1053 F_.appendRow (F_alloc.begin(), F_alloc.end());
1055 flowSolution_.outflux_
1056 .appendRow(F_alloc.begin(), F_alloc.end());
1062 void enumerateBCDof(
const GridInterface& g,
const BCInterface& bc)
1065 typedef typename GridInterface::CellIterator CI;
1066 typedef typename CI ::FaceIterator FI;
1068 const std::vector<int>& cell = flowSolution_.cellno_;
1069 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1071 bdry_id_map_.clear();
1072 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1073 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1074 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1075 const int bid = f->boundaryId();
1076 DofID dof(cell[c->index()], f->localIndex());
1077 bdry_id_map_.insert(std::make_pair(bid, dof));
1082 ppartner_dof_.clear();
1083 if (!bdry_id_map_.empty()) {
1084 ppartner_dof_.assign(total_num_faces_, -1);
1085 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1086 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1087 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1088 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1090 BdryIdMapIterator j =
1091 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1092 assert (j != bdry_id_map_.end());
1093 const int dof2 = cf[j->second.first][j->second.second];
1095 ppartner_dof_[dof1] = dof2;
1096 ppartner_dof_[dof2] = dof1;
1106 void allocateConnections(
const BCInterface& bc)
1110 assert(!cleared_state_);
1112 assert (!matrix_structure_valid_);
1115 S_.setSize(total_num_faces_, total_num_faces_);
1116 S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1119 for (
int f = 0; f < total_num_faces_; ++f) {
1120 S_.setrowsize(f, 1);
1123 allocateGridConnections();
1124 allocateBCConnections(bc);
1128 rhs_ .resize(total_num_faces_);
1129 soln_.resize(total_num_faces_);
1134 void allocateGridConnections()
1137 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1139 for (
int c = 0; c < cf.size(); ++c) {
1140 const int nf = cf[c].size();
1141 for (
auto& f : cf[c]) {
1142 S_.incrementrowsize(f, nf - 1);
1149 void allocateBCConnections(
const BCInterface& bc)
1165 typedef typename GridInterface::CellIterator CI;
1166 typedef typename CI ::FaceIterator FI;
1168 const std::vector<int>& cell = flowSolution_.cellno_;
1169 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1171 if (!bdry_id_map_.empty()) {
1174 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1175 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1176 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1178 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1181 BdryIdMapIterator j =
1182 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1183 assert (j != bdry_id_map_.end());
1184 const int c2 = j->second.first;
1185 const int dof2 = cf[c2][j->second.second];
1191 const int ndof = cf.rowSize(c2);
1192 S_.incrementrowsize(dof1, ndof);
1193 for (
int dof = 0; dof < ndof; ++dof) {
1194 int ii = cf[c2][dof];
1195 int pp = ppartner_dof_[ii];
1196 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1197 S_.incrementrowsize(pp, 1);
1199 S_.incrementrowsize(ii, 1);
1211 void setConnections(
const BCInterface& bc)
1214 setGridConnections();
1215 setBCConnections(bc);
1219 const int nc = pgrid_->numberOfCells();
1220 std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1221 std::vector<Scalar>(nc).swap(g_);
1222 std::vector<Scalar>(nc).swap(L_);
1224 matrix_structure_valid_ =
true;
1229 void setGridConnections()
1232 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1235 for (
int c = 0; c < cf.size(); ++c) {
1236 auto fb = cf[c].begin(), fe = cf[c].end();
1238 for (
auto i = fb; i != fe; ++i) {
1239 for (
auto j = fb; j != fe; ++j) {
1240 S_.addindex(*i, *j);
1248 void setBCConnections(
const BCInterface& bc)
1264 typedef typename GridInterface::CellIterator CI;
1265 typedef typename CI ::FaceIterator FI;
1267 const std::vector<int>& cell = flowSolution_.cellno_;
1268 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1270 if (!bdry_id_map_.empty()) {
1273 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1274 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1275 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1277 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1280 BdryIdMapIterator j =
1281 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1282 assert (j != bdry_id_map_.end());
1283 const int c2 = j->second.first;
1284 const int dof2 = cf[c2][j->second.second];
1289 const int ndof = cf.rowSize(c2);
1290 for (
int dof = 0; dof < ndof; ++dof) {
1291 int ii = cf[c2][dof];
1292 int pp = ppartner_dof_[ii];
1293 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1296 S_.addindex(dof1, ii);
1297 S_.addindex(ii, dof1);
1298 S_.addindex(dof2, ii);
1299 S_.addindex(ii, dof2);
1311 template<
class Flu
idInterface>
1312 void assembleDynamic(
const FluidInterface& fl ,
1313 const std::vector<double>& sat,
1314 const BCInterface& bc ,
1315 const std::vector<double>& src)
1318 typedef typename GridInterface::CellIterator CI;
1320 const std::vector<int>& cell = flowSolution_.cellno_;
1321 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1323 std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1324 std::vector<Scalar> e (max_ncf_);
1325 std::vector<Scalar> rhs (max_ncf_);
1326 std::vector<Scalar> gflux (max_ncf_);
1328 std::vector<FaceType> facetype (max_ncf_);
1329 std::vector<Scalar> condval (max_ncf_);
1330 std::vector<int> ppartner (max_ncf_);
1336 std::fill(g_.begin(), g_.end(), Scalar(0.0));
1337 std::fill(L_.begin(), L_.end(), Scalar(0.0));
1338 std::fill(e .begin(), e .end(), Scalar(1.0));
1342 do_regularization_ =
true;
1345 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1346 const int ci = c->index();
1347 const int c0 = cell[ci]; assert (c0 < cf.size());
1348 const int nf = cf[c0].size();
1353 setExternalContrib(c, c0, bc, src[ci], rhs,
1354 facetype, condval, ppartner);
1356 ip_.computeDynamicParams(c, fl, sat);
1358 SharedFortranMatrix S(nf, nf, &data_store[0]);
1359 ip_.getInverseMatrix(c, S);
1361 std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1362 ip_.gravityFlux(c, gflux);
1364 ImmutableFortranMatrix one(nf, 1, &e[0]);
1365 buildCellContrib(c0, one, gflux, S, rhs);
1367 addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1374 void solveLinearSystem(
double residual_tolerance,
int verbosity_level,
int maxit)
1378 Scalar residTol = residual_tolerance;
1380 typedef Dune::BCRSMatrix <MatrixBlockType> MatrixT;
1381 typedef Dune::BlockVector<VectorBlockType> VectorT;
1382 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1385 if (do_regularization_) {
1391#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1392 Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1394 Dune::SeqILU0<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1398 Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
1399 (maxit>0)?maxit:S_.N(), verbosity_level);
1401 Dune::InverseOperatorResult result;
1406 linsolve.apply(soln_, rhs_, result);
1407 if (!result.converged) {
1408 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1409 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1418 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1419 typedef Dune::BlockVector<VectorBlockType> Vector;
1420 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1425#ifndef FIRST_DIAGONAL
1426#define FIRST_DIAGONAL 1
1432#define SMOOTHER_ILU 1
1435#define SMOOTHER_BGS 0
1437#ifndef ANISOTROPIC_3D
1438#define ANISOTROPIC_3D 0
1442 typedef Dune::Amg::FirstDiagonal CouplingMetric;
1444 typedef Dune::Amg::RowSum CouplingMetric;
1448 typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1450 typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1454 typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1457#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1458 typedef Dune::SeqILU<Matrix,Vector,Vector>
Smoother;
1460 typedef Dune::SeqILU0<Matrix,Vector,Vector>
Smoother;
1463 typedef Dune::SeqSSOR<Matrix,Vector,Vector>
Smoother;
1466 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1470 std::unique_ptr<Operator> opS_;
1471 typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1472 std::unique_ptr<PrecondBase> precond_;
1476 void solveLinearSystemAMG(
double residual_tolerance,
int verbosity_level,
1477 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1480 typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1484 Scalar residTol = residual_tolerance;
1488 if (do_regularization_) {
1495 typename Precond::SmootherArgs smootherArgs;
1496 smootherArgs.relaxationFactor = relax;
1498 smootherArgs.overlap = Precond::SmootherArgs::none;
1499 smootherArgs.onthefly =
false;
1501 Criterion criterion;
1502 criterion.setDebugLevel(verbosity_level);
1504 criterion.setDefaultValuesAnisotropic(3, 2);
1506 criterion.setProlongationDampingFactor(prolong_factor);
1507 criterion.setBeta(1e-10);
1508 criterion.setNoPreSmoothSteps(smooth_steps);
1509 criterion.setNoPostSmoothSteps(smooth_steps);
1510 criterion.setGamma(1);
1511 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1514 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1516 Dune::InverseOperatorResult result;
1520 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1521 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1522 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1523 bool isDirichlet=
true;
1524 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1525 if(ci.index()!=ri.index() && *ci!=0.0)
1528 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1532 linsolve.apply(soln_, rhs_, result);
1533 if (!result.converged) {
1534 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1535 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1542 void solveLinearSystemFastAMG(
double residual_tolerance,
int verbosity_level,
1543 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1546 typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1549 Scalar residTol = residual_tolerance;
1553 if (do_regularization_) {
1559 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
1561 typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
1563 criterion.setDebugLevel(verbosity_level);
1565 criterion.setDefaultValuesAnisotropic(3, 2);
1567 criterion.setProlongationDampingFactor(prolong_factor);
1568 criterion.setBeta(1e-10);
1569 Dune::Amg::Parameters parms;
1570 parms.setDebugLevel(verbosity_level);
1571 parms.setNoPreSmoothSteps(smooth_steps);
1572 parms.setNoPostSmoothSteps(smooth_steps);
1573 precond_.reset(
new Precond(*opS_, criterion, parms));
1576 Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1578 Dune::InverseOperatorResult result;
1583 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1584 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1585 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1586 bool isDirichlet=
true;
1587 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1588 if(ci.index()!=ri.index() && *ci!=0.0)
1591 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1595 linsolve.apply(soln_, rhs_, result);
1596 if (!result.converged) {
1597 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1598 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1604 void solveLinearSystemKAMG(
double residual_tolerance,
int verbosity_level,
1605 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1609 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1611 Scalar residTol = residual_tolerance;
1614 if (do_regularization_) {
1621 typename Precond::SmootherArgs smootherArgs;
1622 smootherArgs.relaxationFactor = relax;
1624 smootherArgs.overlap = Precond::SmootherArgs::none;
1625 smootherArgs.onthefly =
false;
1627 Criterion criterion;
1628 criterion.setDebugLevel(verbosity_level);
1630 criterion.setDefaultValuesAnisotropic(3, 2);
1632 criterion.setProlongationDampingFactor(prolong_factor);
1633 criterion.setBeta(1e-10);
1634 criterion.setNoPreSmoothSteps(smooth_steps);
1635 criterion.setNoPostSmoothSteps(smooth_steps);
1636 criterion.setGamma(2);
1637 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1640 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1642 Dune::InverseOperatorResult result;
1646 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1647 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1648 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1649 bool isDirichlet=
true;
1650 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1651 if(ci.index()!=ri.index() && *ci!=0.0)
1654 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1658 linsolve.apply(soln_, rhs_, result);
1659 if (!result.converged) {
1660 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1661 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1669 template<
class Flu
idInterface>
1670 void computePressureAndFluxes(
const FluidInterface& r ,
1671 const std::vector<double>& sat)
1674 typedef typename GridInterface::CellIterator CI;
1676 const std::vector<int>& cell = flowSolution_.cellno_;
1677 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1679 std::vector<Scalar>& p = flowSolution_.pressure_;
1680 Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1683 std::vector<double> pi (max_ncf_);
1684 std::vector<double> gflux(max_ncf_);
1685 std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1688 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1689 const int c0 = cell[c->index()];
1690 const int nf = cf.rowSize(c0);
1696 for (
int i = 0; i < nf; ++i) {
1697 pi[i] = soln_[cf[c0][i]];
1702 std::inner_product(F_[c0].begin(), F_[c0].end(),
1703 pi.begin(), 0.0)) / L_[c0];
1705 std::transform(pi.begin(), pi.end(),
1707 [&p, c0](
const double& input) { return p[c0] - input; });
1714 ip_.computeDynamicParams(c, r, sat);
1716 SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1717 ip_.getInverseMatrix(c, Binv);
1718 vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1722 ip_.gravityFlux(c, gflux);
1723 std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1725 std::plus<Scalar>());
1733 void setExternalContrib(
const typename GridInterface::CellIterator c,
1734 const int c0,
const BCInterface& bc,
1736 std::vector<Scalar>& rhs,
1737 std::vector<FaceType>& facetype,
1738 std::vector<double>& condval,
1739 std::vector<int>& ppartner)
1742 typedef typename GridInterface::CellIterator::FaceIterator FI;
1744 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1746 std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1747 std::fill(facetype.begin(), facetype.end(), Internal );
1748 std::fill(condval .begin(), condval .end(), Scalar(0.0));
1749 std::fill(ppartner.begin(), ppartner.end(), -1 );
1754 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1755 if (f->boundary()) {
1756 const FlowBC& bcond = bc.flowCond(*f);
1757 if (bcond.isDirichlet()) {
1758 facetype[k] = Dirichlet;
1759 condval[k] = bcond.pressure();
1760 do_regularization_ =
false;
1761 }
else if (bcond.isPeriodic()) {
1762 BdryIdMapIterator j =
1763 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1764 assert (j != bdry_id_map_.end());
1766 facetype[k] = Periodic;
1767 condval[k] = bcond.pressureDifference();
1768 ppartner[k] = cf[j->second.first][j->second.second];
1770 assert (bcond.isNeumann());
1771 facetype[k] = Neumann;
1772 rhs[k] = bcond.outflux();
1782 void buildCellContrib(
const int c ,
1783 const ImmutableFortranMatrix& one ,
1784 const std::vector<Scalar>& gflux,
1785 SharedFortranMatrix& S ,
1786 std::vector<Scalar>& rhs)
1790 SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1793 L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1794 g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1797 std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1799 std::minus<Scalar>());
1802 std::transform(rhs.begin(), rhs.end(), Ft.data(),
1804 axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1815 void addCellContrib(
const SharedFortranMatrix& S ,
1816 const std::vector<Scalar>& rhs ,
1817 const std::vector<FaceType>& facetype,
1818 const std::vector<Scalar>& condval ,
1819 const std::vector<int>& ppartner,
1824 for (
auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1828 switch (facetype[r]) {
1834 S_ [ii][ii] = S(r,r);
1835 rhs_[ii] = S(r,r) * condval[r];
1850 assert ((0 <= ppartner[r]) && (ppartner[r] <
int(rhs_.size())));
1851 assert (ii != ppartner[r]);
1854 const double a = S(r,r), b = a * condval[r];
1858 S_ [ ii][ppartner[r]] -= a;
1862 S_ [ppartner[r]][ ii] -= a;
1863 S_ [ppartner[r]][ppartner[r]] += a;
1864 rhs_[ppartner[r]] -= b;
1867 ii = std::min(ii, ppartner[r]);
1874 for (
auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1878 if (facetype[c] == Dirichlet) {
1879 rhs_[ii] -= S(r,c) * condval[c];
1882 if (facetype[c] == Periodic) {
1883 assert ((0 <= ppartner[c]) && (ppartner[c] <
int(rhs_.size())));
1884 assert (jj != ppartner[c]);
1885 if (ppartner[c] < jj) {
1886 rhs_[ii] -= S(r,c) * condval[c];
1890 S_[ii][jj] += S(r,c);
Solve mixed formulation of incompressible flow modelled by Darcy's law.
Definition: IncompFlowSolverHybrid.hpp:365
double postProcessFluxes()
Postprocess the solution fluxes.
Definition: IncompFlowSolverHybrid.hpp:746
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:820
void solve(const FluidInterface &r, const std::vector< double > &sat, const BCInterface &bc, const std::vector< double > &src, double residual_tolerance=1e-8, int linsolver_verbosity=1, int linsolver_type=1, bool same_matrix=false, int linsolver_maxit=0, double prolongate_factor=1.6, int smooth_steps=1)
Construct and solve system of linear equations for the pressure values on each interface/contact betw...
Definition: IncompFlowSolverHybrid.hpp:657
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:810
void printSystem(const std::string &prefix)
Output current system of linear equations to permanent storage in files.
Definition: IncompFlowSolverHybrid.hpp:887
void computeInnerProducts(const RockInterface &r, const Point &grav)
Compute static (i.e., independent of saturation) parts of the spatially varying inner products for e...
Definition: IncompFlowSolverHybrid.hpp:564
void clear()
Clear all topologic, geometric and rock-dependent information currently held in internal data structu...
Definition: IncompFlowSolverHybrid.hpp:497
void printStats(std::basic_ostream< charT, traits > &os)
Print statistics about the connections in the current model.
Definition: IncompFlowSolverHybrid.hpp:841
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine.
Definition: IncompFlowSolverHybrid.hpp:476
void initSystemStructure(const GridInterface &g, const BCInterface &bc)
Compute structure of coefficient matrix in final system of linear equations for this flow problem.
Definition: IncompFlowSolverHybrid.hpp:533
Dune::MatrixAdapter< Matrix, Vector, Vector > Operator
A linear operator.
Definition: elasticity_preconditioners.hpp:49
Smoother
Smoother used in the AMG.
Definition: elasticity_upscale.hpp:75
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
void symmetricUpdate(const T &a1, const FullMatrix< T, StoragePolicy, FortranOrdering > &A, const T &a2, FullMatrix< T, StoragePolicy, FortranOrdering > &C)
Symmetric, rank update of symmetric matrix.
Definition: Matrix.hpp:829
void vecMulAdd_N(const T &a1, const FullMatrix< T, SP, FortranOrdering > &A, const std::vector< T > &x, const T &a2, std::vector< T > &y)
GEneral Matrix-Vector product (GAXPY operation).
Definition: Matrix.hpp:913
void matMulAdd_TN(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition: Matrix.hpp:1252