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);
966 int tot_ncf = 0, tot_ncf2 = 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);
991 tot_ncf2 += ncf * ncf;
993 assert (cellno == nc);
995 total_num_faces_ = num_internal_faces_ = int(faces.size());
997 ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
998 F_ .reserve(nc, tot_ncf);
1000 flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1001 flowSolution_.outflux_ .reserve(nc, tot_ncf);
1003 Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1006 std::vector<int> l2g; l2g .reserve(max_ncf_);
1007 std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1010 typedef std::vector<int>::iterator VII;
1011 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1012 const int c0 = c->index();
1014 assert ((0 <= c0 ) && ( c0 < nc) &&
1015 (0 <= cell[c0]) && (cell[c0] < nc));
1017 const int ncf = num_cf[cell[c0]];
1018 l2g .resize(ncf , 0 );
1019 F_alloc .resize(ncf , Scalar(0.0));
1021 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1022 if (f->boundary()) {
1024 l2g[f->localIndex()] = total_num_faces_++;
1036 const int c1 = f->neighbourCellIndex();
1037 assert ((0 <= c1 ) && ( c1 < nc) &&
1038 (0 <= cell[c1]) && (cell[c1] < nc));
1040 int t = c0, seek = c1;
1041 if (cell[seek] < cell[t])
1044 int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1046 VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1047 assert(p != faces.begin() + e);
1049 l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1053 cf.appendRow (l2g .begin(), l2g .end());
1054 F_.appendRow (F_alloc.begin(), F_alloc.end());
1056 flowSolution_.outflux_
1057 .appendRow(F_alloc.begin(), F_alloc.end());
1063 void enumerateBCDof(
const GridInterface& g,
const BCInterface& bc)
1066 typedef typename GridInterface::CellIterator CI;
1067 typedef typename CI ::FaceIterator FI;
1069 const std::vector<int>& cell = flowSolution_.cellno_;
1070 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1072 bdry_id_map_.clear();
1073 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1074 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1075 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1076 const int bid = f->boundaryId();
1077 DofID dof(cell[c->index()], f->localIndex());
1078 bdry_id_map_.insert(std::make_pair(bid, dof));
1083 ppartner_dof_.clear();
1084 if (!bdry_id_map_.empty()) {
1085 ppartner_dof_.assign(total_num_faces_, -1);
1086 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1087 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1088 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1089 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1091 BdryIdMapIterator j =
1092 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1093 assert (j != bdry_id_map_.end());
1094 const int dof2 = cf[j->second.first][j->second.second];
1096 ppartner_dof_[dof1] = dof2;
1097 ppartner_dof_[dof2] = dof1;
1107 void allocateConnections(
const BCInterface& bc)
1111 assert(!cleared_state_);
1113 assert (!matrix_structure_valid_);
1116 S_.setSize(total_num_faces_, total_num_faces_);
1117 S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1120 for (
int f = 0; f < total_num_faces_; ++f) {
1121 S_.setrowsize(f, 1);
1124 allocateGridConnections();
1125 allocateBCConnections(bc);
1129 rhs_ .resize(total_num_faces_);
1130 soln_.resize(total_num_faces_);
1135 void allocateGridConnections()
1138 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1140 for (
int c = 0; c < cf.size(); ++c) {
1141 const int nf = cf[c].size();
1142 for (
auto& f : cf[c]) {
1143 S_.incrementrowsize(f, nf - 1);
1150 void allocateBCConnections(
const BCInterface& bc)
1166 typedef typename GridInterface::CellIterator CI;
1167 typedef typename CI ::FaceIterator FI;
1169 const std::vector<int>& cell = flowSolution_.cellno_;
1170 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1172 if (!bdry_id_map_.empty()) {
1175 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1176 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1177 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1179 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1182 BdryIdMapIterator j =
1183 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1184 assert (j != bdry_id_map_.end());
1185 const int c2 = j->second.first;
1186 const int dof2 = cf[c2][j->second.second];
1192 const int ndof = cf.rowSize(c2);
1193 S_.incrementrowsize(dof1, ndof);
1194 for (
int dof = 0; dof < ndof; ++dof) {
1195 int ii = cf[c2][dof];
1196 int pp = ppartner_dof_[ii];
1197 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1198 S_.incrementrowsize(pp, 1);
1200 S_.incrementrowsize(ii, 1);
1212 void setConnections(
const BCInterface& bc)
1215 setGridConnections();
1216 setBCConnections(bc);
1220 const int nc = pgrid_->numberOfCells();
1221 std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1222 std::vector<Scalar>(nc).swap(g_);
1223 std::vector<Scalar>(nc).swap(L_);
1225 matrix_structure_valid_ =
true;
1230 void setGridConnections()
1233 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1236 for (
int c = 0; c < cf.size(); ++c) {
1237 auto fb = cf[c].begin(), fe = cf[c].end();
1239 for (
auto i = fb; i != fe; ++i) {
1240 for (
auto j = fb; j != fe; ++j) {
1241 S_.addindex(*i, *j);
1249 void setBCConnections(
const BCInterface& bc)
1265 typedef typename GridInterface::CellIterator CI;
1266 typedef typename CI ::FaceIterator FI;
1268 const std::vector<int>& cell = flowSolution_.cellno_;
1269 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1271 if (!bdry_id_map_.empty()) {
1274 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1275 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1276 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1278 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1281 BdryIdMapIterator j =
1282 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1283 assert (j != bdry_id_map_.end());
1284 const int c2 = j->second.first;
1285 const int dof2 = cf[c2][j->second.second];
1290 const int ndof = cf.rowSize(c2);
1291 for (
int dof = 0; dof < ndof; ++dof) {
1292 int ii = cf[c2][dof];
1293 int pp = ppartner_dof_[ii];
1294 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1297 S_.addindex(dof1, ii);
1298 S_.addindex(ii, dof1);
1299 S_.addindex(dof2, ii);
1300 S_.addindex(ii, dof2);
1312 template<
class Flu
idInterface>
1313 void assembleDynamic(
const FluidInterface& fl ,
1314 const std::vector<double>& sat,
1315 const BCInterface& bc ,
1316 const std::vector<double>& src)
1319 typedef typename GridInterface::CellIterator CI;
1321 const std::vector<int>& cell = flowSolution_.cellno_;
1322 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1324 std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1325 std::vector<Scalar> e (max_ncf_);
1326 std::vector<Scalar> rhs (max_ncf_);
1327 std::vector<Scalar> gflux (max_ncf_);
1329 std::vector<FaceType> facetype (max_ncf_);
1330 std::vector<Scalar> condval (max_ncf_);
1331 std::vector<int> ppartner (max_ncf_);
1337 std::fill(g_.begin(), g_.end(), Scalar(0.0));
1338 std::fill(L_.begin(), L_.end(), Scalar(0.0));
1339 std::fill(e .begin(), e .end(), Scalar(1.0));
1343 do_regularization_ =
true;
1346 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1347 const int ci = c->index();
1348 const int c0 = cell[ci]; assert (c0 < cf.size());
1349 const int nf = cf[c0].size();
1354 setExternalContrib(c, c0, bc, src[ci], rhs,
1355 facetype, condval, ppartner);
1357 ip_.computeDynamicParams(c, fl, sat);
1359 SharedFortranMatrix S(nf, nf, &data_store[0]);
1360 ip_.getInverseMatrix(c, S);
1362 std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1363 ip_.gravityFlux(c, gflux);
1365 ImmutableFortranMatrix one(nf, 1, &e[0]);
1366 buildCellContrib(c0, one, gflux, S, rhs);
1368 addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1375 void solveLinearSystem(
double residual_tolerance,
int verbosity_level,
int maxit)
1379 Scalar residTol = residual_tolerance;
1381 typedef Dune::BCRSMatrix <MatrixBlockType> MatrixT;
1382 typedef Dune::BlockVector<VectorBlockType> VectorT;
1383 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1386 if (do_regularization_) {
1392 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1393 Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1395 Dune::SeqILU0<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1399 Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
1400 (maxit>0)?maxit:S_.N(), verbosity_level);
1402 Dune::InverseOperatorResult result;
1407 linsolve.apply(soln_, rhs_, result);
1408 if (!result.converged) {
1409 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1410 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1419 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1420 typedef Dune::BlockVector<VectorBlockType> Vector;
1421 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1426 #ifndef FIRST_DIAGONAL
1427 #define FIRST_DIAGONAL 1
1432 #ifndef SMOOTHER_ILU
1433 #define SMOOTHER_ILU 1
1435 #ifndef SMOOTHER_BGS
1436 #define SMOOTHER_BGS 0
1438 #ifndef ANISOTROPIC_3D
1439 #define ANISOTROPIC_3D 0
1443 typedef Dune::Amg::FirstDiagonal CouplingMetric;
1445 typedef Dune::Amg::RowSum CouplingMetric;
1449 typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1451 typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1455 typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1458 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1459 typedef Dune::SeqILU<Matrix,Vector,Vector>
Smoother;
1461 typedef Dune::SeqILU0<Matrix,Vector,Vector>
Smoother;
1464 typedef Dune::SeqSSOR<Matrix,Vector,Vector>
Smoother;
1467 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1471 std::unique_ptr<Operator> opS_;
1472 typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1473 std::unique_ptr<PrecondBase> precond_;
1477 void solveLinearSystemAMG(
double residual_tolerance,
int verbosity_level,
1478 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1481 typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1485 Scalar residTol = residual_tolerance;
1489 if (do_regularization_) {
1496 typename Precond::SmootherArgs smootherArgs;
1497 smootherArgs.relaxationFactor = relax;
1499 smootherArgs.overlap = Precond::SmootherArgs::none;
1500 smootherArgs.onthefly =
false;
1502 Criterion criterion;
1503 criterion.setDebugLevel(verbosity_level);
1505 criterion.setDefaultValuesAnisotropic(3, 2);
1507 criterion.setProlongationDampingFactor(prolong_factor);
1508 criterion.setBeta(1e-10);
1509 criterion.setNoPreSmoothSteps(smooth_steps);
1510 criterion.setNoPostSmoothSteps(smooth_steps);
1511 criterion.setGamma(1);
1512 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1515 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1517 Dune::InverseOperatorResult result;
1521 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1522 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1523 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1524 bool isDirichlet=
true;
1525 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1526 if(ci.index()!=ri.index() && *ci!=0.0)
1529 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1533 linsolve.apply(soln_, rhs_, result);
1534 if (!result.converged) {
1535 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1536 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1543 void solveLinearSystemFastAMG(
double residual_tolerance,
int verbosity_level,
1544 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1547 typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1550 Scalar residTol = residual_tolerance;
1554 if (do_regularization_) {
1560 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
1562 typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
1564 criterion.setDebugLevel(verbosity_level);
1566 criterion.setDefaultValuesAnisotropic(3, 2);
1568 criterion.setProlongationDampingFactor(prolong_factor);
1569 criterion.setBeta(1e-10);
1570 Dune::Amg::Parameters parms;
1571 parms.setDebugLevel(verbosity_level);
1572 parms.setNoPreSmoothSteps(smooth_steps);
1573 parms.setNoPostSmoothSteps(smooth_steps);
1574 precond_.reset(
new Precond(*opS_, criterion, parms));
1577 Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1579 Dune::InverseOperatorResult result;
1584 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1585 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1586 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1587 bool isDirichlet=
true;
1588 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1589 if(ci.index()!=ri.index() && *ci!=0.0)
1592 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1596 linsolve.apply(soln_, rhs_, result);
1597 if (!result.converged) {
1598 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1599 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1605 void solveLinearSystemKAMG(
double residual_tolerance,
int verbosity_level,
1606 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1610 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1612 Scalar residTol = residual_tolerance;
1615 if (do_regularization_) {
1622 typename Precond::SmootherArgs smootherArgs;
1623 smootherArgs.relaxationFactor = relax;
1625 smootherArgs.overlap = Precond::SmootherArgs::none;
1626 smootherArgs.onthefly =
false;
1628 Criterion criterion;
1629 criterion.setDebugLevel(verbosity_level);
1631 criterion.setDefaultValuesAnisotropic(3, 2);
1633 criterion.setProlongationDampingFactor(prolong_factor);
1634 criterion.setBeta(1e-10);
1635 criterion.setNoPreSmoothSteps(smooth_steps);
1636 criterion.setNoPostSmoothSteps(smooth_steps);
1637 criterion.setGamma(2);
1638 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1641 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1643 Dune::InverseOperatorResult result;
1647 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1648 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1649 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1650 bool isDirichlet=
true;
1651 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1652 if(ci.index()!=ri.index() && *ci!=0.0)
1655 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1659 linsolve.apply(soln_, rhs_, result);
1660 if (!result.converged) {
1661 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1662 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1670 template<
class Flu
idInterface>
1671 void computePressureAndFluxes(
const FluidInterface& r ,
1672 const std::vector<double>& sat)
1675 typedef typename GridInterface::CellIterator CI;
1677 const std::vector<int>& cell = flowSolution_.cellno_;
1678 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1680 std::vector<Scalar>& p = flowSolution_.pressure_;
1681 Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1684 std::vector<double> pi (max_ncf_);
1685 std::vector<double> gflux(max_ncf_);
1686 std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1689 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1690 const int c0 = cell[c->index()];
1691 const int nf = cf.rowSize(c0);
1697 for (
int i = 0; i < nf; ++i) {
1698 pi[i] = soln_[cf[c0][i]];
1703 std::inner_product(F_[c0].begin(), F_[c0].end(),
1704 pi.begin(), 0.0)) / L_[c0];
1706 std::transform(pi.begin(), pi.end(),
1708 [&p, c0](
const double& input) { return p[c0] - input; });
1715 ip_.computeDynamicParams(c, r, sat);
1717 SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1718 ip_.getInverseMatrix(c, Binv);
1719 vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1723 ip_.gravityFlux(c, gflux);
1724 std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1726 std::plus<Scalar>());
1734 void setExternalContrib(
const typename GridInterface::CellIterator c,
1735 const int c0,
const BCInterface& bc,
1737 std::vector<Scalar>& rhs,
1738 std::vector<FaceType>& facetype,
1739 std::vector<double>& condval,
1740 std::vector<int>& ppartner)
1743 typedef typename GridInterface::CellIterator::FaceIterator FI;
1745 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1747 std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1748 std::fill(facetype.begin(), facetype.end(), Internal );
1749 std::fill(condval .begin(), condval .end(), Scalar(0.0));
1750 std::fill(ppartner.begin(), ppartner.end(), -1 );
1755 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1756 if (f->boundary()) {
1757 const FlowBC& bcond = bc.flowCond(*f);
1758 if (bcond.isDirichlet()) {
1759 facetype[k] = Dirichlet;
1760 condval[k] = bcond.pressure();
1761 do_regularization_ =
false;
1762 }
else if (bcond.isPeriodic()) {
1763 BdryIdMapIterator j =
1764 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1765 assert (j != bdry_id_map_.end());
1767 facetype[k] = Periodic;
1768 condval[k] = bcond.pressureDifference();
1769 ppartner[k] = cf[j->second.first][j->second.second];
1771 assert (bcond.isNeumann());
1772 facetype[k] = Neumann;
1773 rhs[k] = bcond.outflux();
1783 void buildCellContrib(
const int c ,
1784 const ImmutableFortranMatrix& one ,
1785 const std::vector<Scalar>& gflux,
1786 SharedFortranMatrix& S ,
1787 std::vector<Scalar>& rhs)
1791 SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1794 L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1795 g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1798 std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1800 std::minus<Scalar>());
1803 std::transform(rhs.begin(), rhs.end(), Ft.data(),
1805 axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1816 void addCellContrib(
const SharedFortranMatrix& S ,
1817 const std::vector<Scalar>& rhs ,
1818 const std::vector<FaceType>& facetype,
1819 const std::vector<Scalar>& condval ,
1820 const std::vector<int>& ppartner,
1825 for (
auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1829 switch (facetype[r]) {
1835 S_ [ii][ii] = S(r,r);
1836 rhs_[ii] = S(r,r) * condval[r];
1851 assert ((0 <= ppartner[r]) && (ppartner[r] <
int(rhs_.size())));
1852 assert (ii != ppartner[r]);
1855 const double a = S(r,r), b = a * condval[r];
1859 S_ [ ii][ppartner[r]] -= a;
1863 S_ [ppartner[r]][ ii] -= a;
1864 S_ [ppartner[r]][ppartner[r]] += a;
1865 rhs_[ppartner[r]] -= b;
1868 ii = std::min(ii, ppartner[r]);
1875 for (
auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1879 if (facetype[c] == Dirichlet) {
1880 rhs_[ii] -= S(r,c) * condval[c];
1883 if (facetype[c] == Periodic) {
1884 assert ((0 <= ppartner[c]) && (ppartner[c] <
int(rhs_.size())));
1885 assert (jj != ppartner[c]);
1886 if (ppartner[c] < jj) {
1887 rhs_[ii] -= S(r,c) * condval[c];
1891 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:45
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