38 #ifndef OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
39 #define OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
41 #include <opm/common/ErrorMacros.hpp>
42 #include <opm/core/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 <boost/unordered_map.hpp>
49 #include <boost/bind.hpp>
50 #include <boost/scoped_ptr.hpp>
51 #include <boost/lexical_cast.hpp>
53 #include <dune/common/fvector.hh>
54 #include <dune/common/fmatrix.hh>
56 #include <dune/istl/bvector.hh>
57 #include <dune/istl/bcrsmatrix.hh>
58 #include <dune/istl/operators.hh>
59 #include <dune/istl/io.hh>
60 #include <dune/istl/overlappingschwarz.hh>
61 #include <dune/istl/schwarz.hh>
62 #include <dune/istl/preconditioners.hh>
63 #include <dune/istl/solvers.hh>
64 #include <dune/istl/owneroverlapcopy.hh>
65 #include <dune/istl/paamg/amg.hh>
66 #include <dune/common/version.hh>
67 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
68 #include <dune/istl/paamg/fastamg.hh>
70 #include <dune/istl/paamg/kamg.hh>
71 #include <dune/istl/paamg/pinfo.hh>
73 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
112 bool topologyIsSane(
const GI& g)
114 typedef typename GI::CellIterator CI;
115 typedef typename CI::FaceIterator FI;
117 bool sane = g.numberOfCells() >= 0;
119 for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
122 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
123 if (!f->boundary()) {
124 n.push_back(f->neighbourCellIndex());
127 std::sort(n.begin(), n.end());
129 sane = std::unique(n.begin(), n.end()) == n.end();
162 class axpby :
public std::binary_function<T,T,T> {
169 axpby(
const T& a,
const T& b) : a_(a), b_(b) {}
183 T operator()(
const T& x,
const T& y)
365 template <
class GridInterface,
368 template <
class Gr
idIF,
class RockIF>
class InnerProduct>
375 typedef typename GridInterface::Scalar Scalar;
381 enum FaceType { Internal, Dirichlet, Neumann, Periodic };
392 typedef typename GridInterface::Scalar Scalar;
397 typedef typename GridInterface::CellIterator CI;
401 typedef typename CI ::FaceIterator FI;
414 Scalar pressure(
const CI& c)
const
416 return pressure_[cellno_[c->index()]];
429 Scalar outflux (
const FI& f)
const
431 return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
433 Scalar outflux (
int hf)
const
435 return outflux_.data(hf);
438 std::vector< int > cellno_;
439 Opm::SparseTable< int > cellFaces_;
440 std::vector<Scalar> pressure_;
441 Opm::SparseTable<Scalar> outflux_;
444 std::vector<int>().swap(cellno_);
447 std::vector<Scalar>().swap(pressure_);
479 template<
class Po
int>
480 void init(
const GridInterface& g,
481 const RockInterface& r,
483 const BCInterface& bc)
487 if (g.numberOfCells() > 0) {
505 num_internal_faces_ = 0;
506 total_num_faces_ = 0;
507 matrix_structure_valid_ =
false;
508 do_regularization_ =
true;
510 bdry_id_map_.clear();
512 std::vector<Scalar>().swap(L_);
513 std::vector<Scalar>().swap(g_);
516 flowSolution_.clear();
518 cleared_state_ =
true;
538 const BCInterface& bc)
541 assert (cleared_state_);
543 assert (topologyIsSane(g));
546 allocateConnections(bc);
567 template<
class Po
int>
572 assert (matrix_structure_valid_);
574 typedef typename GridInterface::CellIterator CI;
575 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
578 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
579 ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
660 template<
class Flu
idInterface>
661 void solve(
const FluidInterface& r ,
662 const std::vector<double>& sat,
663 const BCInterface& bc ,
664 const std::vector<double>& src,
665 double residual_tolerance = 1e-8,
666 int linsolver_verbosity = 1,
667 int linsolver_type = 1,
668 bool same_matrix =
false,
669 int linsolver_maxit = 0,
670 double prolongate_factor = 1.6,
671 int smooth_steps = 1)
673 assembleDynamic(r, sat, bc, src);
677 switch (linsolver_type) {
679 solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
682 solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
683 linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
687 solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
688 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
691 #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
692 solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
693 linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
695 if(linsolver_verbosity)
696 std::cerr<<
"Fast AMG is not available; falling back to CG preconditioned with the normal one."<<std::endl;
697 solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
698 linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
702 std::cerr <<
"Unknown linsolver_type: " << linsolver_type <<
'\n';
703 throw std::runtime_error(
"Unknown linsolver_type");
705 computePressureAndFluxes(r, sat);
714 : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
717 void put(
double flux,
int f_ix) {
718 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
719 double sign = visited_[f_ix] ? -1.0 : 1.0;
720 fluxes_[f_ix] += sign*flux;
723 void get(
double& flux,
int f_ix) {
724 assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
725 double sign = visited_[f_ix] ? -1.0 : 1.0;
726 double new_flux = 0.5*sign*fluxes_[f_ix];
727 double diff = std::fabs(flux - new_flux);
728 max_modification_ = std::max(max_modification_, diff);
734 std::fill(visited_.begin(), visited_.end(), 0);
737 double maxMod()
const
739 return max_modification_;
742 std::vector<double> fluxes_;
743 std::vector<int> visited_;
744 double max_modification_;
759 typedef typename GridInterface::CellIterator CI;
760 typedef typename CI ::FaceIterator FI;
761 const std::vector<int>& cell = flowSolution_.cellno_;
762 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
763 Opm::SparseTable<double>& cflux = flowSolution_.outflux_;
765 FaceFluxes face_fluxes(pgrid_->numberOfFaces());
767 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
768 const int cell_index = cell[c->index()];
769 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
770 int f_ix = cf[cell_index][f->localIndex()];
771 double flux = cflux[cell_index][f->localIndex()];
773 if (ppartner_dof_.empty()) {
776 int partner_f_ix = ppartner_dof_[f_ix];
777 if (partner_f_ix != -1) {
778 face_fluxes.put(flux, f_ix);
779 face_fluxes.put(flux, partner_f_ix);
782 face_fluxes.put(flux, f_ix);
786 face_fluxes.resetVisited();
788 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
789 const int cell_index = cell[c->index()];
790 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
791 int f_ix = cf[cell_index][f->localIndex()];
792 double& flux = cflux[cell_index][f->localIndex()];
794 if (ppartner_dof_.empty()) {
797 int partner_f_ix = ppartner_dof_[f_ix];
798 if (partner_f_ix != -1) {
799 face_fluxes.get(flux, f_ix);
801 face_fluxes.get(dummy, partner_f_ix);
802 assert(dummy == flux);
805 face_fluxes.get(flux, f_ix);
809 return face_fluxes.maxMod();
833 return flowSolution_;
851 template<
typename charT,
class traits>
854 os <<
"IncompFlowSolverHybrid<>:\n"
855 <<
"\tMaximum number of cell faces = " << max_ncf_ <<
'\n'
856 <<
"\tNumber of internal faces = " << num_internal_faces_ <<
'\n'
857 <<
"\tTotal number of faces = " << total_num_faces_ <<
'\n';
859 const std::vector<int>& cell = flowSolution_.cellno_;
860 os <<
"cell index map = [";
861 std::copy(cell.begin(), cell.end(),
862 std::ostream_iterator<int>(os,
" "));
865 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
866 os <<
"cell faces =\n";
867 for (
int i = 0; i < cf.size(); ++i)
869 os <<
"\t[" << i <<
"] -> [";
870 std::copy(cf[i].begin(), cf[i].end(),
871 std::ostream_iterator<int>(os,
","));
900 writeMatrixToMatlab(S_, prefix +
"-mat.dat");
902 std::string rhsfile(prefix +
"-rhs.dat");
903 std::ofstream rhs(rhsfile.c_str());
905 rhs.setf(std::ios::scientific | std::ios::showpos);
906 std::copy(rhs_.begin(), rhs_.end(),
907 std::ostream_iterator<VectorBlockType>(rhs,
"\n"));
911 typedef std::pair<int,int> DofID;
912 typedef boost::unordered_map<int,DofID> BdryIdMapType;
913 typedef BdryIdMapType::const_iterator BdryIdMapIterator;
915 const GridInterface* pgrid_;
916 BdryIdMapType bdry_id_map_;
917 std::vector<int> ppartner_dof_;
919 InnerProduct<GridInterface, RockInterface> ip_;
924 int num_internal_faces_;
925 int total_num_faces_;
928 std::vector<Scalar> L_, g_;
929 Opm::SparseTable<Scalar> F_ ;
933 typedef Dune::FieldVector<Scalar, 1 > VectorBlockType;
934 typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
936 Dune::BCRSMatrix <MatrixBlockType> S_;
937 Dune::BlockVector<VectorBlockType> rhs_;
938 Dune::BlockVector<VectorBlockType> soln_;
939 bool matrix_structure_valid_;
940 bool do_regularization_;
944 FlowSolution flowSolution_;
948 void enumerateDof(
const GridInterface& g,
const BCInterface& bc)
952 enumerateBCDof(g, bc);
955 cleared_state_ =
false;
959 void enumerateGridDof(
const GridInterface& g)
962 typedef typename GridInterface::CellIterator CI;
963 typedef typename CI ::FaceIterator FI;
965 const int nc = g.numberOfCells();
966 std::vector<int> fpos ; fpos.reserve(nc + 1);
967 std::vector<int> num_cf(nc) ;
968 std::vector<int> faces ;
971 std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
973 std::vector<int>& cell = flowSolution_.cellno_;
976 int cellno = 0; fpos.push_back(0);
977 int tot_ncf = 0, tot_ncf2 = 0;
978 for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++cellno) {
979 const int c0 = c->index();
980 assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
984 int& ncf = num_cf[c0];
986 for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
987 if (!f->boundary()) {
988 const int c1 = f->neighbourCellIndex();
989 assert((0 <= c1) && (c1 < nc) && (c1 != c0));
991 if (cell[c1] == -1) {
999 fpos.push_back(
int(faces.size()));
1000 max_ncf_ = std::max(max_ncf_, ncf);
1002 tot_ncf2 += ncf * ncf;
1004 assert (cellno == nc);
1006 total_num_faces_ = num_internal_faces_ = int(faces.size());
1008 ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
1009 F_ .reserve(nc, tot_ncf);
1011 flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1012 flowSolution_.outflux_ .reserve(nc, tot_ncf);
1014 Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1017 std::vector<int> l2g; l2g .reserve(max_ncf_);
1018 std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1021 typedef std::vector<int>::iterator VII;
1022 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1023 const int c0 = c->index();
1025 assert ((0 <= c0 ) && ( c0 < nc) &&
1026 (0 <= cell[c0]) && (cell[c0] < nc));
1028 const int ncf = num_cf[cell[c0]];
1029 l2g .resize(ncf , 0 );
1030 F_alloc .resize(ncf , Scalar(0.0));
1032 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1033 if (f->boundary()) {
1035 l2g[f->localIndex()] = total_num_faces_++;
1047 const int c1 = f->neighbourCellIndex();
1048 assert ((0 <= c1 ) && ( c1 < nc) &&
1049 (0 <= cell[c1]) && (cell[c1] < nc));
1051 int t = c0, seek = c1;
1052 if (cell[seek] < cell[t])
1055 int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1057 VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1058 assert(p != faces.begin() + e);
1060 l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1064 cf.appendRow (l2g .begin(), l2g .end());
1065 F_.appendRow (F_alloc.begin(), F_alloc.end());
1067 flowSolution_.outflux_
1068 .appendRow(F_alloc.begin(), F_alloc.end());
1074 void enumerateBCDof(
const GridInterface& g,
const BCInterface& bc)
1077 typedef typename GridInterface::CellIterator CI;
1078 typedef typename CI ::FaceIterator FI;
1080 const std::vector<int>& cell = flowSolution_.cellno_;
1081 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1083 bdry_id_map_.clear();
1084 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1085 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1086 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1087 const int bid = f->boundaryId();
1088 DofID dof(cell[c->index()], f->localIndex());
1089 bdry_id_map_.insert(std::make_pair(bid, dof));
1094 ppartner_dof_.clear();
1095 if (!bdry_id_map_.empty()) {
1096 ppartner_dof_.assign(total_num_faces_, -1);
1097 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1098 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1099 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1100 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1102 BdryIdMapIterator j =
1103 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1104 assert (j != bdry_id_map_.end());
1105 const int dof2 = cf[j->second.first][j->second.second];
1107 ppartner_dof_[dof1] = dof2;
1108 ppartner_dof_[dof2] = dof1;
1118 void allocateConnections(
const BCInterface& bc)
1122 assert(!cleared_state_);
1124 assert (!matrix_structure_valid_);
1127 S_.setSize(total_num_faces_, total_num_faces_);
1128 S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1131 for (
int f = 0; f < total_num_faces_; ++f) {
1132 S_.setrowsize(f, 1);
1135 allocateGridConnections();
1136 allocateBCConnections(bc);
1140 rhs_ .resize(total_num_faces_);
1141 soln_.resize(total_num_faces_);
1146 void allocateGridConnections()
1149 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1150 typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1152 for (
int c = 0; c < cf.size(); ++c) {
1153 const int nf = cf[c].size();
1154 fi fb = cf[c].begin(), fe = cf[c].end();
1156 for (fi f = fb; f != fe; ++f) {
1157 S_.incrementrowsize(*f, nf - 1);
1164 void allocateBCConnections(
const BCInterface& bc)
1180 typedef typename GridInterface::CellIterator CI;
1181 typedef typename CI ::FaceIterator FI;
1183 const std::vector<int>& cell = flowSolution_.cellno_;
1184 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1186 if (!bdry_id_map_.empty()) {
1189 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1190 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1191 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1193 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1196 BdryIdMapIterator j =
1197 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1198 assert (j != bdry_id_map_.end());
1199 const int c2 = j->second.first;
1200 const int dof2 = cf[c2][j->second.second];
1206 const int ndof = cf.rowSize(c2);
1207 S_.incrementrowsize(dof1, ndof);
1208 for (
int dof = 0; dof < ndof; ++dof) {
1209 int ii = cf[c2][dof];
1210 int pp = ppartner_dof_[ii];
1211 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1212 S_.incrementrowsize(pp, 1);
1214 S_.incrementrowsize(ii, 1);
1226 void setConnections(
const BCInterface& bc)
1229 setGridConnections();
1230 setBCConnections(bc);
1234 const int nc = pgrid_->numberOfCells();
1235 std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1236 std::vector<Scalar>(nc).swap(g_);
1237 std::vector<Scalar>(nc).swap(L_);
1239 matrix_structure_valid_ =
true;
1244 void setGridConnections()
1247 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1248 typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1251 for (
int c = 0; c < cf.size(); ++c) {
1252 fi fb = cf[c].begin(), fe = cf[c].end();
1254 for (fi i = fb; i != fe; ++i) {
1255 for (fi j = fb; j != fe; ++j) {
1256 S_.addindex(*i, *j);
1264 void setBCConnections(
const BCInterface& bc)
1280 typedef typename GridInterface::CellIterator CI;
1281 typedef typename CI ::FaceIterator FI;
1283 const std::vector<int>& cell = flowSolution_.cellno_;
1284 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1286 if (!bdry_id_map_.empty()) {
1289 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1290 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1291 if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1293 const int dof1 = cf[cell[c->index()]][f->localIndex()];
1296 BdryIdMapIterator j =
1297 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1298 assert (j != bdry_id_map_.end());
1299 const int c2 = j->second.first;
1300 const int dof2 = cf[c2][j->second.second];
1305 const int ndof = cf.rowSize(c2);
1306 for (
int dof = 0; dof < ndof; ++dof) {
1307 int ii = cf[c2][dof];
1308 int pp = ppartner_dof_[ii];
1309 if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1312 S_.addindex(dof1, ii);
1313 S_.addindex(ii, dof1);
1314 S_.addindex(dof2, ii);
1315 S_.addindex(ii, dof2);
1327 template<
class Flu
idInterface>
1328 void assembleDynamic(
const FluidInterface& fl ,
1329 const std::vector<double>& sat,
1330 const BCInterface& bc ,
1331 const std::vector<double>& src)
1334 typedef typename GridInterface::CellIterator CI;
1336 const std::vector<int>& cell = flowSolution_.cellno_;
1337 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1339 std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1340 std::vector<Scalar> e (max_ncf_);
1341 std::vector<Scalar> rhs (max_ncf_);
1342 std::vector<Scalar> gflux (max_ncf_);
1344 std::vector<FaceType> facetype (max_ncf_);
1345 std::vector<Scalar> condval (max_ncf_);
1346 std::vector<int> ppartner (max_ncf_);
1352 std::fill(g_.begin(), g_.end(), Scalar(0.0));
1353 std::fill(L_.begin(), L_.end(), Scalar(0.0));
1354 std::fill(e .begin(), e .end(), Scalar(1.0));
1358 do_regularization_ =
true;
1361 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1362 const int ci = c->index();
1363 const int c0 = cell[ci]; assert (c0 < cf.size());
1364 const int nf = cf[c0].size();
1369 setExternalContrib(c, c0, bc, src[ci], rhs,
1370 facetype, condval, ppartner);
1372 ip_.computeDynamicParams(c, fl, sat);
1374 SharedFortranMatrix S(nf, nf, &data_store[0]);
1375 ip_.getInverseMatrix(c, S);
1377 std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1378 ip_.gravityFlux(c, gflux);
1380 ImmutableFortranMatrix one(nf, 1, &e[0]);
1381 buildCellContrib(c0, one, gflux, S, rhs);
1383 addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1390 void solveLinearSystem(
double residual_tolerance,
int verbosity_level,
int maxit)
1394 Scalar residTol = residual_tolerance;
1396 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1397 typedef Dune::BlockVector<VectorBlockType> Vector;
1398 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1401 if (do_regularization_) {
1407 Dune::SeqILU0<Matrix,Vector,Vector> precond(S_, 1.0);
1410 Dune::CGSolver<Vector> linsolve(opS, precond, residTol,
1411 (maxit>0)?maxit:S_.N(), verbosity_level);
1413 Dune::InverseOperatorResult result;
1418 linsolve.apply(soln_, rhs_, result);
1419 if (!result.converged) {
1420 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1421 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1430 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1431 typedef Dune::BlockVector<VectorBlockType> Vector;
1432 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1437 #ifndef FIRST_DIAGONAL
1438 #define FIRST_DIAGONAL 1
1443 #ifndef SMOOTHER_ILU
1444 #define SMOOTHER_ILU 1
1446 #ifndef SMOOTHER_BGS
1447 #define SMOOTHER_BGS 0
1449 #ifndef ANISOTROPIC_3D
1450 #define ANISOTROPIC_3D 0
1454 typedef Dune::Amg::FirstDiagonal CouplingMetric;
1456 typedef Dune::Amg::RowSum CouplingMetric;
1460 typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1462 typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1466 typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1469 typedef Dune::SeqILU0<Matrix,Vector,Vector> Smoother;
1471 typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1474 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1478 boost::scoped_ptr<Operator> opS_;
1479 typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1480 boost::scoped_ptr<PrecondBase> precond_;
1484 void solveLinearSystemAMG(
double residual_tolerance,
int verbosity_level,
1485 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1488 typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1492 Scalar residTol = residual_tolerance;
1496 if (do_regularization_) {
1499 opS_.reset(
new Operator(S_));
1503 typename Precond::SmootherArgs smootherArgs;
1504 smootherArgs.relaxationFactor = relax;
1506 smootherArgs.overlap = Precond::SmootherArgs::none;
1507 smootherArgs.onthefly =
false;
1509 Criterion criterion;
1510 criterion.setDebugLevel(verbosity_level);
1512 criterion.setDefaultValuesAnisotropic(3, 2);
1514 criterion.setProlongationDampingFactor(prolong_factor);
1515 criterion.setBeta(1e-10);
1516 criterion.setNoPreSmoothSteps(smooth_steps);
1517 criterion.setNoPostSmoothSteps(smooth_steps);
1518 criterion.setGamma(1);
1519 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1522 Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1524 Dune::InverseOperatorResult result;
1528 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1529 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1530 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1531 bool isDirichlet=
true;
1532 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1533 if(ci.index()!=ri.index() && *ci!=0.0)
1536 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1540 linsolve.apply(soln_, rhs_, result);
1541 if (!result.converged) {
1542 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1543 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1548 #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
1551 void solveLinearSystemFastAMG(
double residual_tolerance,
int verbosity_level,
1552 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1555 typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1558 Scalar residTol = residual_tolerance;
1562 if (do_regularization_) {
1565 opS_.reset(
new Operator(S_));
1568 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CriterionBase;
1570 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1571 Criterion criterion;
1572 criterion.setDebugLevel(verbosity_level);
1574 criterion.setDefaultValuesAnisotropic(3, 2);
1576 criterion.setProlongationDampingFactor(prolong_factor);
1577 criterion.setBeta(1e-10);
1578 Dune::Amg::Parameters parms;
1579 parms.setDebugLevel(verbosity_level);
1580 parms.setNoPreSmoothSteps(smooth_steps);
1581 parms.setNoPostSmoothSteps(smooth_steps);
1582 precond_.reset(
new Precond(*opS_, criterion, parms));
1585 Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1587 Dune::InverseOperatorResult result;
1592 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1593 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1594 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1595 bool isDirichlet=
true;
1596 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1597 if(ci.index()!=ri.index() && *ci!=0.0)
1600 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1604 linsolve.apply(soln_, rhs_, result);
1605 if (!result.converged) {
1606 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1607 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1614 void solveLinearSystemKAMG(
double residual_tolerance,
int verbosity_level,
1615 int maxit,
double prolong_factor,
bool same_matrix,
int smooth_steps)
1619 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1621 Scalar residTol = residual_tolerance;
1624 if (do_regularization_) {
1627 opS_.reset(
new Operator(S_));
1631 typename Precond::SmootherArgs smootherArgs;
1632 smootherArgs.relaxationFactor = relax;
1634 smootherArgs.overlap = Precond::SmootherArgs::none;
1635 smootherArgs.onthefly =
false;
1637 Criterion criterion;
1638 criterion.setDebugLevel(verbosity_level);
1640 criterion.setDefaultValuesAnisotropic(3, 2);
1642 criterion.setProlongationDampingFactor(prolong_factor);
1643 criterion.setBeta(1e-10);
1644 criterion.setNoPreSmoothSteps(smooth_steps);
1645 criterion.setNoPostSmoothSteps(smooth_steps);
1646 criterion.setGamma(2);
1647 precond_.reset(
new Precond(*opS_, criterion, smootherArgs));
1650 Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1652 Dune::InverseOperatorResult result;
1656 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1657 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1658 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1659 bool isDirichlet=
true;
1660 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1661 if(ci.index()!=ri.index() && *ci!=0.0)
1664 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1668 linsolve.apply(soln_, rhs_, result);
1669 if (!result.converged) {
1670 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
1671 <<
"Residual reduction achieved is " << result.reduction <<
'\n');
1679 template<
class Flu
idInterface>
1680 void computePressureAndFluxes(
const FluidInterface& r ,
1681 const std::vector<double>& sat)
1684 typedef typename GridInterface::CellIterator CI;
1686 const std::vector<int>& cell = flowSolution_.cellno_;
1687 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1689 std::vector<Scalar>& p = flowSolution_.pressure_;
1690 Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1693 std::vector<double> pi (max_ncf_);
1694 std::vector<double> gflux(max_ncf_);
1695 std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1698 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1699 const int c0 = cell[c->index()];
1700 const int nf = cf.rowSize(c0);
1706 for (
int i = 0; i < nf; ++i) {
1707 pi[i] = soln_[cf[c0][i]];
1712 std::inner_product(F_[c0].begin(), F_[c0].end(),
1713 pi.begin(), 0.0)) / L_[c0];
1715 std::transform(pi.begin(), pi.end(),
1717 boost::bind(std::minus<Scalar>(), p[c0], _1));
1724 ip_.computeDynamicParams(c, r, sat);
1726 SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1727 ip_.getInverseMatrix(c, Binv);
1728 vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1732 ip_.gravityFlux(c, gflux);
1733 std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1735 std::plus<Scalar>());
1743 void setExternalContrib(
const typename GridInterface::CellIterator c,
1744 const int c0,
const BCInterface& bc,
1746 std::vector<Scalar>& rhs,
1747 std::vector<FaceType>& facetype,
1748 std::vector<double>& condval,
1749 std::vector<int>& ppartner)
1752 typedef typename GridInterface::CellIterator::FaceIterator FI;
1754 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1756 std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1757 std::fill(facetype.begin(), facetype.end(), Internal );
1758 std::fill(condval .begin(), condval .end(), Scalar(0.0));
1759 std::fill(ppartner.begin(), ppartner.end(), -1 );
1764 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1765 if (f->boundary()) {
1766 const FlowBC& bcond = bc.flowCond(*f);
1767 if (bcond.isDirichlet()) {
1768 facetype[k] = Dirichlet;
1769 condval[k] = bcond.pressure();
1770 do_regularization_ =
false;
1771 }
else if (bcond.isPeriodic()) {
1772 BdryIdMapIterator j =
1773 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1774 assert (j != bdry_id_map_.end());
1776 facetype[k] = Periodic;
1777 condval[k] = bcond.pressureDifference();
1778 ppartner[k] = cf[j->second.first][j->second.second];
1780 assert (bcond.isNeumann());
1781 facetype[k] = Neumann;
1782 rhs[k] = bcond.outflux();
1792 void buildCellContrib(
const int c ,
1793 const ImmutableFortranMatrix& one ,
1794 const std::vector<Scalar>& gflux,
1795 SharedFortranMatrix& S ,
1796 std::vector<Scalar>& rhs)
1800 SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1803 L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1804 g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1807 std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1809 std::minus<Scalar>());
1812 std::transform(rhs.begin(), rhs.end(), Ft.data(),
1814 axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1825 void addCellContrib(
const SharedFortranMatrix& S ,
1826 const std::vector<Scalar>& rhs ,
1827 const std::vector<FaceType>& facetype,
1828 const std::vector<Scalar>& condval ,
1829 const std::vector<int>& ppartner,
1833 typedef typename L2G::const_iterator it;
1836 for (it i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1840 switch (facetype[r]) {
1846 S_ [ii][ii] = S(r,r);
1847 rhs_[ii] = S(r,r) * condval[r];
1862 assert ((0 <= ppartner[r]) && (ppartner[r] <
int(rhs_.size())));
1863 assert (ii != ppartner[r]);
1866 const double a = S(r,r), b = a * condval[r];
1870 S_ [ ii][ppartner[r]] -= a;
1874 S_ [ppartner[r]][ ii] -= a;
1875 S_ [ppartner[r]][ppartner[r]] += a;
1876 rhs_[ppartner[r]] -= b;
1879 ii = std::min(ii, ppartner[r]);
1886 for (it j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1890 if (facetype[c] == Dirichlet) {
1891 rhs_[ii] -= S(r,c) * condval[c];
1894 if (facetype[c] == Periodic) {
1895 assert ((0 <= ppartner[c]) && (ppartner[c] <
int(rhs_.size())));
1896 assert (jj != ppartner[c]);
1897 if (ppartner[c] < jj) {
1898 rhs_[ii] -= S(r,c) * condval[c];
1902 S_[ii][jj] += S(r,c);
1912 #endif // OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
void printSystem(const std::string &prefix)
Output current system of linear equations to permanent storage in files.
Definition: IncompFlowSolverHybrid.hpp:898
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:537
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:821
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:568
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:661
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:1253
void printStats(std::basic_ostream< charT, traits > &os)
Print statistics about the connections in the current model.
Definition: IncompFlowSolverHybrid.hpp:852
double postProcessFluxes()
Postprocess the solution fluxes.
Definition: IncompFlowSolverHybrid.hpp:757
void clear()
Clear all topologic, geometric and rock-dependent information currently held in internal data structu...
Definition: IncompFlowSolverHybrid.hpp:501
Solve mixed formulation of incompressible flow modelled by Darcy's law ] The solver is based on a ...
Definition: IncompFlowSolverHybrid.hpp:369
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:831
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:830
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:914
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine.
Definition: IncompFlowSolverHybrid.hpp:480