38#ifndef OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
39#define OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
41#include <boost/unordered_map.hpp>
43#include <boost/bind.hpp>
44#include <boost/scoped_ptr.hpp>
45#include <boost/lexical_cast.hpp>
47#include <dune/common/fvector.hh>
48#include <dune/common/fmatrix.hh>
49#include <opm/common/ErrorMacros.hpp>
50#include <opm/core/utility/SparseTable.hpp>
52#include <dune/istl/bvector.hh>
53#include <dune/istl/bcrsmatrix.hh>
54#include <dune/istl/operators.hh>
55#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#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
65#include <dune/istl/paamg/fastamg.hh>
67#include <dune/istl/paamg/kamg.hh>
68#include <dune/istl/paamg/pinfo.hh>
109 bool topologyIsSane(
const GI& g)
111 typedef typename GI::CellIterator CI;
112 typedef typename CI::FaceIterator FI;
114 bool sane = g.numberOfCells() >= 0;
116 for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
119 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
120 if (!f->boundary()) {
121 n.push_back(f->neighbourCellIndex());
124 std::sort(n.begin(), n.end());
126 sane = std::unique(n.begin(), n.end()) == n.end();
159 class axpby :
public std::binary_function<T,T,T> {
166 axpby(
const T& a,
const T& b) : a_(a), b_(b) {}
180 T operator()(
const T& x,
const T& y)
362 template <
class GridInterface,
365 template <
class Gr
idIF,
class RockIF>
class InnerProduct>
372 typedef typename GridInterface::Scalar Scalar;
378 enum FaceType { Internal, Dirichlet, Neumann, Periodic };
389 typedef typename GridInterface::Scalar Scalar;
394 typedef typename GridInterface::CellIterator CI;
398 typedef typename CI ::FaceIterator FI;
411 Scalar pressure(
const CI& c)
const
413 return pressure_[cellno_[c->index()]];
426 Scalar outflux (
const FI& f)
const
428 return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
430 Scalar outflux (
int hf)
const
432 return outflux_.data(hf);
435 std::vector< int > cellno_;
436 Opm::SparseTable< int > cellFaces_;
437 std::vector<Scalar> pressure_;
438 Opm::SparseTable<Scalar> outflux_;
441 std::vector<int>().swap(cellno_);
444 std::vector<Scalar>().swap(pressure_);
476 template<
class Po
int>
477 void init(
const GridInterface& g,
478 const RockInterface& r,
480 const BCInterface& bc)
484 if (g.numberOfCells() > 0) {
502 num_internal_faces_ = 0;
503 total_num_faces_ = 0;
504 matrix_structure_valid_ =
false;
505 do_regularization_ =
true;
507 bdry_id_map_.clear();
509 std::vector<Scalar>().swap(L_);
510 std::vector<Scalar>().swap(g_);
513 flowSolution_.clear();
515 cleared_state_ =
true;
535 const BCInterface& bc)
538 assert (cleared_state_);
540 assert (topologyIsSane(g));
543 allocateConnections(bc);
564 template<
class Po
int>
569 assert (matrix_structure_valid_);
571 typedef typename GridInterface::CellIterator CI;
572 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
575 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
576 ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
657 template<
class Flu
idInterface>
658 void solve(
const FluidInterface& r ,
659 const std::vector<double>& sat,
660 const BCInterface& bc ,
661 const std::vector<double>& src,
662 double residual_tolerance = 1e-8,
663 int linsolver_verbosity = 1,
664 int linsolver_type = 1,
665 bool same_matrix =
false,
666 int linsolver_maxit = 0,
667 double prolongate_factor = 1.6,
668 int smooth_steps = 1)
670 assembleDynamic(r, sat, bc, src);
674 switch (linsolver_type) {
676 solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
679 solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
680 linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
684 solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
685 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
688#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
689 solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
690 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
692 if(linsolver_verbosity)
693 std::cerr<<
"Fast AMG is not available; falling back to CG preconditioned with the normal one."<<std::endl;
694 solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
695 linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
699 std::cerr <<
"Unknown linsolver_type: " << linsolver_type <<
'\n';
700 throw std::runtime_error(
"Unknown linsolver_type");
702 computePressureAndFluxes(r, sat);
711 : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
714 void put(
double flux,
int f_ix) {
715 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
716 double sign = visited_[f_ix] ? -1.0 : 1.0;
717 fluxes_[f_ix] += sign*flux;
720 void get(
double& flux,
int f_ix) {
721 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
722 double sign = visited_[f_ix] ? -1.0 : 1.0;
723 double new_flux = 0.5*sign*fluxes_[f_ix];
724 double diff = std::fabs(flux - new_flux);
725 max_modification_ = std::max(max_modification_, diff);
731 std::fill(visited_.begin(), visited_.end(), 0);
734 double maxMod()
const
736 return max_modification_;
739 std::vector<double> fluxes_;
740 std::vector<int> visited_;
741 double max_modification_;
756 typedef typename GridInterface::CellIterator CI;
757 typedef typename CI ::FaceIterator FI;
758 const std::vector<int>& cell = flowSolution_.cellno_;
759 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
760 Opm::SparseTable<double>& cflux = flowSolution_.outflux_;
762 FaceFluxes face_fluxes(pgrid_->numberOfFaces());
764 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
765 const int cell_index = cell[c->index()];
766 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
767 int f_ix = cf[cell_index][f->localIndex()];
768 double flux = cflux[cell_index][f->localIndex()];
770 if (ppartner_dof_.empty()) {
773 int partner_f_ix = ppartner_dof_[f_ix];
774 if (partner_f_ix != -1) {
775 face_fluxes.put(flux, f_ix);
776 face_fluxes.put(flux, partner_f_ix);
779 face_fluxes.put(flux, f_ix);
783 face_fluxes.resetVisited();
785 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
786 const int cell_index = cell[c->index()];
787 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
788 int f_ix = cf[cell_index][f->localIndex()];
789 double& flux = cflux[cell_index][f->localIndex()];
791 if (ppartner_dof_.empty()) {
794 int partner_f_ix = ppartner_dof_[f_ix];
795 if (partner_f_ix != -1) {
796 face_fluxes.get(flux, f_ix);
798 face_fluxes.get(dummy, partner_f_ix);
799 assert(dummy == flux);
802 face_fluxes.get(flux, f_ix);
806 return face_fluxes.maxMod();
830 return flowSolution_;
848 template<
typename charT,
class traits>
851 os <<
"IncompFlowSolverHybrid<>:\n"
852 <<
"\tMaximum number of cell faces = " << max_ncf_ <<
'\n'
853 <<
"\tNumber of internal faces = " << num_internal_faces_ <<
'\n'
854 <<
"\tTotal number of faces = " << total_num_faces_ <<
'\n';
856 const std::vector<int>& cell = flowSolution_.cellno_;
857 os <<
"cell index map = [";
858 std::copy(cell.begin(), cell.end(),
859 std::ostream_iterator<int>(os,
" "));
862 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
863 os <<
"cell faces =\n";
864 for (
int i = 0; i < cf.size(); ++i)
866 os <<
"\t[" << i <<
"] -> [";
867 std::copy(cf[i].begin(), cf[i].end(),
868 std::ostream_iterator<int>(os,
","));
897 writeMatrixToMatlab(S_, prefix +
"-mat.dat");
899 std::string rhsfile(prefix +
"-rhs.dat");
900 std::ofstream rhs(rhsfile.c_str());
902 rhs.setf(std::ios::scientific | std::ios::showpos);
903 std::copy(rhs_.begin(), rhs_.end(),
904 std::ostream_iterator<VectorBlockType>(rhs,
"\n"));
908 typedef std::pair<int,int> DofID;
909 typedef boost::unordered_map<int,DofID> BdryIdMapType;
910 typedef BdryIdMapType::const_iterator BdryIdMapIterator;
912 const GridInterface* pgrid_;
913 BdryIdMapType bdry_id_map_;
914 std::vector<int> ppartner_dof_;
916 InnerProduct<GridInterface, RockInterface> ip_;
921 int num_internal_faces_;
922 int total_num_faces_;
925 std::vector<Scalar> L_, g_;
926 Opm::SparseTable<Scalar> F_ ;
930 typedef Dune::FieldVector<Scalar, 1 > VectorBlockType;
931 typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
933 Dune::BCRSMatrix <MatrixBlockType> S_;
934 Dune::BlockVector<VectorBlockType> rhs_;
935 Dune::BlockVector<VectorBlockType> soln_;
936 bool matrix_structure_valid_;
937 bool do_regularization_;
941 FlowSolution flowSolution_;
945 void enumerateDof(
const GridInterface& g,
const BCInterface& bc)
949 enumerateBCDof(g, bc);
952 cleared_state_ =
false;
956 void enumerateGridDof(
const GridInterface& g)
959 typedef typename GridInterface::CellIterator CI;
960 typedef typename CI ::FaceIterator FI;
962 const int nc = g.numberOfCells();
963 std::vector<int> fpos ; fpos.reserve(nc + 1);
964 std::vector<int> num_cf(nc) ;
965 std::vector<int> faces ;
968 std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
970 std::vector<int>& cell = flowSolution_.cellno_;
973 int cellno = 0; fpos.push_back(0);
974 int tot_ncf = 0, tot_ncf2 = 0;
975 for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++cellno) {
976 const int c0 = c->index();
977 assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
981 int& ncf = num_cf[c0];
983 for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
984 if (!f->boundary()) {
985 const int c1 = f->neighbourCellIndex();
986 assert((0 <= c1) && (c1 < nc) && (c1 != c0));
988 if (cell[c1] == -1) {
996 fpos.push_back(
int(faces.size()));
997 max_ncf_ = std::max(max_ncf_, ncf);
999 tot_ncf2 += ncf * ncf;
1001 assert (cellno == nc);
1003 total_num_faces_ = num_internal_faces_ = int(faces.size());
1005 ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
1006 F_ .reserve(nc, tot_ncf);
1008 flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1009 flowSolution_.outflux_ .reserve(nc, tot_ncf);
1011 Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1014 std::vector<int> l2g; l2g .reserve(max_ncf_);
1015 std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1018 typedef std::vector<int>::iterator VII;
1019 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1020 const int c0 = c->index();
1022 assert ((0 <= c0 ) && ( c0 < nc) &&
1023 (0 <= cell[c0]) && (cell[c0] < nc));
1025 const int ncf = num_cf[cell[c0]];
1026 l2g .resize(ncf , 0 );
1027 F_alloc .resize(ncf , Scalar(0.0));
1029 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1030 if (f->boundary()) {
1032 l2g[f->localIndex()] = total_num_faces_++;
1044 const int c1 = f->neighbourCellIndex();
1045 assert ((0 <= c1 ) && ( c1 < nc) &&
1046 (0 <= cell[c1]) && (cell[c1] < nc));
1048 int t = c0, seek = c1;
1049 if (cell[seek] < cell[t])
1052 int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1054 VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1055 assert(p != faces.begin() + e);
1057 l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1061 cf.appendRow (l2g .begin(), l2g .end());
1062 F_.appendRow (F_alloc.begin(), F_alloc.end());
1064 flowSolution_.outflux_
1065 .appendRow(F_alloc.begin(), F_alloc.end());
1071 void enumerateBCDof(
const GridInterface& g,
const BCInterface& bc)
1074 typedef typename GridInterface::CellIterator CI;
1075 typedef typename CI ::FaceIterator FI;
1077 const std::vector<int>& cell = flowSolution_.cellno_;
1078 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1080 bdry_id_map_.clear();
1081 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1082 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1083 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1084 const int bid = f->boundaryId();
1085 DofID dof(cell[c->index()], f->localIndex());
1086 bdry_id_map_.insert(std::make_pair(bid, dof));
1091 ppartner_dof_.clear();
1092 if (!bdry_id_map_.empty()) {
1093 ppartner_dof_.assign(total_num_faces_, -1);
1094 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1095 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1096 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1097 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1099 BdryIdMapIterator j =
1100 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1101 assert (j != bdry_id_map_.end());
1102 const int dof2 = cf[j->second.first][j->second.second];
1104 ppartner_dof_[dof1] = dof2;
1105 ppartner_dof_[dof2] = dof1;
1115 void allocateConnections(
const BCInterface& bc)
1119 assert(!cleared_state_);
1121 assert (!matrix_structure_valid_);
1124 S_.setSize(total_num_faces_, total_num_faces_);
1125 S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1128 for (
int f = 0; f < total_num_faces_; ++f) {
1129 S_.setrowsize(f, 1);
1132 allocateGridConnections();
1133 allocateBCConnections(bc);
1137 rhs_ .resize(total_num_faces_);
1138 soln_.resize(total_num_faces_);
1143 void allocateGridConnections()
1146 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1147 typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1149 for (
int c = 0; c < cf.size(); ++c) {
1150 const int nf = cf[c].size();
1151 fi fb = cf[c].begin(), fe = cf[c].end();
1153 for (fi f = fb; f != fe; ++f) {
1154 S_.incrementrowsize(*f, nf - 1);
1161 void allocateBCConnections(
const BCInterface& bc)
1177 typedef typename GridInterface::CellIterator CI;
1178 typedef typename CI ::FaceIterator FI;
1180 const std::vector<int>& cell = flowSolution_.cellno_;
1181 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1183 if (!bdry_id_map_.empty()) {
1186 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1187 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1188 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1190 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1193 BdryIdMapIterator j =
1194 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1195 assert (j != bdry_id_map_.end());
1196 const int c2 = j->second.first;
1197 const int dof2 = cf[c2][j->second.second];
1203 const int ndof = cf.rowSize(c2);
1204 S_.incrementrowsize(dof1, ndof);
1205 for (
int dof = 0; dof < ndof; ++dof) {
1206 int ii = cf[c2][dof];
1207 int pp = ppartner_dof_[ii];
1208 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1209 S_.incrementrowsize(pp, 1);
1211 S_.incrementrowsize(ii, 1);
1223 void setConnections(
const BCInterface& bc)
1226 setGridConnections();
1227 setBCConnections(bc);
1231 const int nc = pgrid_->numberOfCells();
1232 std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1233 std::vector<Scalar>(nc).swap(g_);
1234 std::vector<Scalar>(nc).swap(L_);
1236 matrix_structure_valid_ =
true;
1241 void setGridConnections()
1244 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1245 typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1248 for (
int c = 0; c < cf.size(); ++c) {
1249 fi fb = cf[c].begin(), fe = cf[c].end();
1251 for (fi i = fb; i != fe; ++i) {
1252 for (fi j = fb; j != fe; ++j) {
1253 S_.addindex(*i, *j);
1261 void setBCConnections(
const BCInterface& bc)
1277 typedef typename GridInterface::CellIterator CI;
1278 typedef typename CI ::FaceIterator FI;
1280 const std::vector<int>& cell = flowSolution_.cellno_;
1281 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1283 if (!bdry_id_map_.empty()) {
1286 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1287 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1288 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1290 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1293 BdryIdMapIterator j =
1294 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1295 assert (j != bdry_id_map_.end());
1296 const int c2 = j->second.first;
1297 const int dof2 = cf[c2][j->second.second];
1302 const int ndof = cf.rowSize(c2);
1303 for (
int dof = 0; dof < ndof; ++dof) {
1304 int ii = cf[c2][dof];
1305 int pp = ppartner_dof_[ii];
1306 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1309 S_.addindex(dof1, ii);
1310 S_.addindex(ii, dof1);
1311 S_.addindex(dof2, ii);
1312 S_.addindex(ii, dof2);
1324 template<
class Flu
idInterface>
1325 void assembleDynamic(
const FluidInterface& fl ,
1326 const std::vector<double>& sat,
1327 const BCInterface& bc ,
1328 const std::vector<double>& src)
1331 typedef typename GridInterface::CellIterator CI;
1333 const std::vector<int>& cell = flowSolution_.cellno_;
1334 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1336 std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1337 std::vector<Scalar> e (max_ncf_);
1338 std::vector<Scalar> rhs (max_ncf_);
1339 std::vector<Scalar> gflux (max_ncf_);
1341 std::vector<FaceType> facetype (max_ncf_);
1342 std::vector<Scalar> condval (max_ncf_);
1343 std::vector<int> ppartner (max_ncf_);
1349 std::fill(g_.begin(), g_.end(), Scalar(0.0));
1350 std::fill(L_.begin(), L_.end(), Scalar(0.0));
1351 std::fill(e .begin(), e .end(), Scalar(1.0));
1355 do_regularization_ =
true;
1358 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1359 const int ci = c->index();
1360 const int c0 = cell[ci]; assert (c0 < cf.size());
1361 const int nf = cf[c0].size();
1366 setExternalContrib(c, c0, bc, src[ci], rhs,
1367 facetype, condval, ppartner);
1369 ip_.computeDynamicParams(c, fl, sat);
1372 ip_.getInverseMatrix(c, S);
1374 std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1375 ip_.gravityFlux(c, gflux);
1378 buildCellContrib(c0, one, gflux, S, rhs);
1380 addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1387 void solveLinearSystem(
double residual_tolerance,
int verbosity_level,
int maxit)
1391 Scalar residTol = residual_tolerance;
1393 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1394 typedef Dune::BlockVector<VectorBlockType> Vector;
1395 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1398 if (do_regularization_) {
1404 Dune::SeqILU0<Matrix,Vector,Vector> precond(S_, 1.0);
1407 Dune::CGSolver<Vector> linsolve(opS, precond, residTol,
1408 (maxit>0)?maxit:S_.N(), verbosity_level);
1410 Dune::InverseOperatorResult result;
1415 linsolve.apply(soln_, rhs_, result);
1416 if (!result.converged) {
1417 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1418 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1427 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1428 typedef Dune::BlockVector<VectorBlockType> Vector;
1429 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1434#ifndef FIRST_DIAGONAL
1435#define FIRST_DIAGONAL 1
1441#define SMOOTHER_ILU 1
1444#define SMOOTHER_BGS 0
1446#ifndef ANISOTROPIC_3D
1447#define ANISOTROPIC_3D 0
1451 typedef Dune::Amg::FirstDiagonal CouplingMetric;
1453 typedef Dune::Amg::RowSum CouplingMetric;
1457 typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1459 typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1463 typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1466 typedef Dune::SeqILU0<Matrix,Vector,Vector> Smoother;
1468 typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1471 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1475 boost::scoped_ptr<Operator> opS_;
1476 typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1477 boost::scoped_ptr<PrecondBase> precond_;
1481 void solveLinearSystemAMG(
double residual_tolerance,
int verbosity_level,
1482 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1485 typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1489 Scalar residTol = residual_tolerance;
1493 if (do_regularization_) {
1496 opS_.reset(
new Operator(S_));
1500 typename Precond::SmootherArgs smootherArgs;
1501 smootherArgs.relaxationFactor = relax;
1503 smootherArgs.overlap = Precond::SmootherArgs::none;
1504 smootherArgs.onthefly =
false;
1506 Criterion criterion;
1507 criterion.setDebugLevel(verbosity_level);
1509 criterion.setDefaultValuesAnisotropic(3, 2);
1511 criterion.setProlongationDampingFactor(prolong_factor);
1512 criterion.setBeta(1e-10);
1513 precond_.reset(
new Precond(*opS_, criterion, smootherArgs,
1514 1, smooth_steps, smooth_steps));
1517 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1519 Dune::InverseOperatorResult result;
1523 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1524 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1525 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1526 bool isDirichlet=
true;
1527 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1528 if(ci.index()!=ri.index() && *ci!=0.0)
1531 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1535 linsolve.apply(soln_, rhs_, result);
1536 if (!result.converged) {
1537 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1538 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1543#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
1546 void solveLinearSystemFastAMG(
double residual_tolerance,
int verbosity_level,
1547 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1550 typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1553 Scalar residTol = residual_tolerance;
1557 if (do_regularization_) {
1560 opS_.reset(
new Operator(S_));
1563 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CriterionBase;
1565 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1566 Criterion criterion;
1567 criterion.setDebugLevel(verbosity_level);
1569 criterion.setDefaultValuesAnisotropic(3, 2);
1571 criterion.setProlongationDampingFactor(prolong_factor);
1572 criterion.setBeta(1e-10);
1573 Dune::Amg::Parameters parms;
1574 parms.setDebugLevel(verbosity_level);
1575 parms.setNoPreSmoothSteps(smooth_steps);
1576 parms.setNoPostSmoothSteps(smooth_steps);
1577 precond_.reset(
new Precond(*opS_, criterion, parms));
1580 Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1582 Dune::InverseOperatorResult result;
1587 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1588 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1589 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1590 bool isDirichlet=
true;
1591 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1592 if(ci.index()!=ri.index() && *ci!=0.0)
1595 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1599 linsolve.apply(soln_, rhs_, result);
1600 if (!result.converged) {
1601 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1602 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1609 void solveLinearSystemKAMG(
double residual_tolerance,
int verbosity_level,
1610 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1614 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation,
1615 Dune::CGSolver<Vector> > Precond;
1617 Scalar residTol = residual_tolerance;
1620 if (do_regularization_) {
1623 opS_.reset(
new Operator(S_));
1627 typename Precond::SmootherArgs smootherArgs;
1628 smootherArgs.relaxationFactor = relax;
1630 smootherArgs.overlap = Precond::SmootherArgs::none;
1631 smootherArgs.onthefly =
false;
1633 Criterion criterion;
1634 criterion.setDebugLevel(verbosity_level);
1636 criterion.setDefaultValuesAnisotropic(3, 2);
1638 criterion.setProlongationDampingFactor(prolong_factor);
1639 criterion.setBeta(1e-10);
1640 precond_.reset(
new Precond(*opS_, criterion, smootherArgs, 2, smooth_steps, smooth_steps));
1643 Dune::CGSolver<Vector> linsolve(*opS_,
dynamic_cast<Precond&
>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1645 Dune::InverseOperatorResult result;
1649 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1650 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1651 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1652 bool isDirichlet=
true;
1653 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1654 if(ci.index()!=ri.index() && *ci!=0.0)
1657 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1661 linsolve.apply(soln_, rhs_, result);
1662 if (!result.converged) {
1663 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1664 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1672 template<
class Flu
idInterface>
1673 void computePressureAndFluxes(
const FluidInterface& r ,
1674 const std::vector<double>& sat)
1677 typedef typename GridInterface::CellIterator CI;
1679 const std::vector<int>& cell = flowSolution_.cellno_;
1680 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1682 std::vector<Scalar>& p = flowSolution_.pressure_;
1683 Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1686 std::vector<double> pi (max_ncf_);
1687 std::vector<double> gflux(max_ncf_);
1688 std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1691 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1692 const int c0 = cell[c->index()];
1693 const int nf = cf.rowSize(c0);
1699 for (
int i = 0; i < nf; ++i) {
1700 pi[i] = soln_[cf[c0][i]];
1705 std::inner_product(F_[c0].begin(), F_[c0].end(),
1706 pi.begin(), 0.0)) / L_[c0];
1708 std::transform(pi.begin(), pi.end(),
1710 boost::bind(std::minus<Scalar>(), p[c0], _1));
1717 ip_.computeDynamicParams(c, r, sat);
1720 ip_.getInverseMatrix(c, Binv);
1721 vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1725 ip_.gravityFlux(c, gflux);
1726 std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1728 std::plus<Scalar>());
1736 void setExternalContrib(
const typename GridInterface::CellIterator c,
1737 const int c0,
const BCInterface& bc,
1739 std::vector<Scalar>& rhs,
1740 std::vector<FaceType>& facetype,
1741 std::vector<double>& condval,
1742 std::vector<int>& ppartner)
1745 typedef typename GridInterface::CellIterator::FaceIterator FI;
1747 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1749 std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1750 std::fill(facetype.begin(), facetype.end(), Internal );
1751 std::fill(condval .begin(), condval .end(), Scalar(0.0));
1752 std::fill(ppartner.begin(), ppartner.end(), -1 );
1757 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1758 if (f->boundary()) {
1759 const FlowBC& bcond = bc.flowCond(*f);
1760 if (bcond.isDirichlet()) {
1761 facetype[k] = Dirichlet;
1762 condval[k] = bcond.pressure();
1763 do_regularization_ =
false;
1764 }
else if (bcond.isPeriodic()) {
1765 BdryIdMapIterator j =
1766 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1767 assert (j != bdry_id_map_.end());
1769 facetype[k] = Periodic;
1770 condval[k] = bcond.pressureDifference();
1771 ppartner[k] = cf[j->second.first][j->second.second];
1773 assert (bcond.isNeumann());
1774 facetype[k] = Neumann;
1775 rhs[k] = bcond.outflux();
1785 void buildCellContrib(
const int c ,
1787 const std::vector<Scalar>& gflux,
1789 std::vector<Scalar>& rhs)
1796 L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1797 g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1800 std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1802 std::minus<Scalar>());
1805 std::transform(rhs.begin(), rhs.end(), Ft.data(),
1807 axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1819 const std::vector<Scalar>& rhs ,
1820 const std::vector<FaceType>& facetype,
1821 const std::vector<Scalar>& condval ,
1822 const std::vector<int>& ppartner,
1826 typedef typename L2G::const_iterator it;
1829 for (it i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1833 switch (facetype[r]) {
1839 S_ [ii][ii] = S(r,r);
1840 rhs_[ii] = S(r,r) * condval[r];
1855 assert ((0 <= ppartner[r]) && (ppartner[r] <
int(rhs_.size())));
1856 assert (ii != ppartner[r]);
1859 const double a = S(r,r), b = a * condval[r];
1863 S_ [ ii][ppartner[r]] -= a;
1867 S_ [ppartner[r]][ ii] -= a;
1868 S_ [ppartner[r]][ppartner[r]] += a;
1869 rhs_[ppartner[r]] -= b;
1872 ii = std::min(ii, ppartner[r]);
1879 for (it j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1883 if (facetype[c] == Dirichlet) {
1884 rhs_[ii] -= S(r,c) * condval[c];
1887 if (facetype[c] == Periodic) {
1888 assert ((0 <= ppartner[c]) && (ppartner[c] <
int(rhs_.size())));
1889 assert (jj != ppartner[c]);
1890 if (ppartner[c] < jj) {
1891 rhs_[ii] -= S(r,c) * condval[c];
1895 S_[ii][jj] += S(r,c);
Solve mixed formulation of incompressible flow modelled by Darcy's law.
Definition: IncompFlowSolverHybrid.hpp:366
double postProcessFluxes()
Postprocess the solution fluxes. This method modifies the solution object so that out-fluxes of twin ...
Definition: IncompFlowSolverHybrid.hpp:754
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:828
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:658
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:818
void printSystem(const std::string &prefix)
Output current system of linear equations to permanent storage in files. One file for the coefficient...
Definition: IncompFlowSolverHybrid.hpp:895
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:565
void clear()
Clear all topologic, geometric and rock-dependent information currently held in internal data structu...
Definition: IncompFlowSolverHybrid.hpp:498
void printStats(std::basic_ostream< charT, traits > &os)
Print statistics about the connections in the current model. This is mostly for debugging purposes an...
Definition: IncompFlowSolverHybrid.hpp:849
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine. Enumerates all grid connections, allocates sufficient space,...
Definition: IncompFlowSolverHybrid.hpp:477
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:534
Definition: BlackoilFluid.hpp:32
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. Specifically, .
Definition: Matrix.hpp:830
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix
Definition: Matrix.hpp:591
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). Specifically, .
Definition: Matrix.hpp:914
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
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. Specificlly, .
Definition: Matrix.hpp:1253