My Project
IncompFlowSolverHybrid.hpp
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set ts=8 sw=4 et sts=4:
3 //===========================================================================
4 //
5 // File: IncompFlowSolverHybrid.hpp
6 //
7 // Created: Tue Jun 30 10:25:40 2009
8 //
9 // Author(s): B�rd Skaflestad <bard.skaflestad@sintef.no>
10 // Atgeirr F Rasmussen <atgeirr@sintef.no>
11 //
12 // $Date$
13 //
14 // $Revision$
15 //
16 //===========================================================================
17 
18 /*
19  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
20  Copyright 2009, 2010 Statoil ASA.
21 
22  This file is part of The Open Reservoir Simulator Project (OpenRS).
23 
24  OpenRS is free software: you can redistribute it and/or modify
25  it under the terms of the GNU General Public License as published by
26  the Free Software Foundation, either version 3 of the License, or
27  (at your option) any later version.
28 
29  OpenRS is distributed in the hope that it will be useful,
30  but WITHOUT ANY WARRANTY; without even the implied warranty of
31  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32  GNU General Public License for more details.
33 
34  You should have received a copy of the GNU General Public License
35  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
36 */
37 
38 #ifndef OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
39 #define OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
40 
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>
45 
46 #include <opm/common/utility/platform_dependent/disable_warnings.h>
47 
48 #include <unordered_map>
49 
50 #include <dune/common/fvector.hh>
51 #include <dune/common/fmatrix.hh>
52 
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>
67 
68 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
69 
70 
71 #include <algorithm>
72 #include <functional>
73 #include <map>
74 #include <memory>
75 #include <numeric>
76 #include <ostream>
77 #include <stdexcept>
78 #include <utility>
79 #include <vector>
80 #include <iostream>
81 
82 namespace Opm {
83  namespace {
107  template<class GI>
108  bool topologyIsSane(const GI& g)
109  {
110  typedef typename GI::CellIterator CI;
111  typedef typename CI::FaceIterator FI;
112 
113  bool sane = g.numberOfCells() >= 0;
114 
115  for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
116  std::vector<int> n;
117 
118  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
119  if (!f->boundary()) {
120  n.push_back(f->neighbourCellIndex());
121  }
122  }
123  std::sort(n.begin(), n.end());
124 
125  sane = std::unique(n.begin(), n.end()) == n.end();
126  }
127 
128  return sane;
129  }
130 
131 
157  template<typename T>
158  class axpby : public std::binary_function<T,T,T> {
159  public:
165  axpby(const T& a, const T& b) : a_(a), b_(b) {}
166 
179  T operator()(const T& x, const T& y)
180  {
181  return a_*x + b_*y;
182  }
183  private:
184  T a_, b_;
185  };
186  }
187 
188 
361  template <class GridInterface,
362  class RockInterface,
363  class BCInterface,
364  template <class GridIF, class RockIF> class InnerProduct>
371  typedef typename GridInterface::Scalar Scalar;
372 
377  enum FaceType { Internal, Dirichlet, Neumann, Periodic };
378 
381  class FlowSolution {
382  public:
388  typedef typename GridInterface::Scalar Scalar;
389 
393  typedef typename GridInterface::CellIterator CI;
394 
397  typedef typename CI ::FaceIterator FI;
398 
399  friend class IncompFlowSolverHybrid;
400 
410  Scalar pressure(const CI& c) const
411  {
412  return pressure_[cellno_[c->index()]];
413  }
414 
425  Scalar outflux (const FI& f) const
426  {
427  return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
428  }
429  Scalar outflux (int hf) const
430  {
431  return outflux_.data(hf);
432  }
433  private:
434  std::vector< int > cellno_;
435  Opm::SparseTable< int > cellFaces_;
436  std::vector<Scalar> pressure_;
437  Opm::SparseTable<Scalar> outflux_;
438 
439  void clear() {
440  std::vector<int>().swap(cellno_);
441  cellFaces_.clear();
442 
443  std::vector<Scalar>().swap(pressure_);
444  outflux_.clear();
445  }
446  };
447 
448  public:
475  template<class Point>
476  void init(const GridInterface& g,
477  const RockInterface& r,
478  const Point& grav,
479  const BCInterface& bc)
480  {
481  clear();
482 
483  if (g.numberOfCells() > 0) {
484  initSystemStructure(g, bc);
485  computeInnerProducts(r, grav);
486  }
487  }
488 
489 
497  void clear()
498  {
499  pgrid_ = 0;
500  max_ncf_ = -1;
501  num_internal_faces_ = 0;
502  total_num_faces_ = 0;
503  matrix_structure_valid_ = false;
504  do_regularization_ = true; // Assume pure Neumann by default.
505 
506  bdry_id_map_.clear();
507 
508  std::vector<Scalar>().swap(L_);
509  std::vector<Scalar>().swap(g_);
510  F_.clear();
511 
512  flowSolution_.clear();
513 
514  cleared_state_ = true;
515  }
516 
517 
533  void initSystemStructure(const GridInterface& g,
534  const BCInterface& bc)
535  {
536  // You must call clear() prior to initSystemStructure()
537  assert (cleared_state_);
538 
539  assert (topologyIsSane(g));
540 
541  enumerateDof(g, bc);
542  allocateConnections(bc);
543  setConnections(bc);
544  }
545 
546 
563  template<class Point>
564  void computeInnerProducts(const RockInterface& r,
565  const Point& grav)
566  {
567  // You must call connectionsCompleted() prior to computeInnerProducts()
568  assert (matrix_structure_valid_);
569 
570  typedef typename GridInterface::CellIterator CI;
571  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
572 
573  int i = 0;
574  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
575  ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
576  }
577  }
578 
579 
656  template<class FluidInterface>
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)
668  {
669  assembleDynamic(r, sat, bc, src);
670 // static int count = 0;
671 // ++count;
672 // printSystem(std::string("linsys_mimetic-") + std::to_string(count));
673  switch (linsolver_type) {
674  case 0: // ILU0 preconditioned CG
675  solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
676  break;
677  case 1: // AMG preconditioned CG
678  solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
679  linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
680  break;
681 
682  case 2: // KAMG preconditioned CG
683  solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
684  linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
685  break;
686  case 3: // CG preconditioned with AMG that uses less memory badwidth
687  solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
688  linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
689  break;
690  default:
691  std::cerr << "Unknown linsolver_type: " << linsolver_type << '\n';
692  throw std::runtime_error("Unknown linsolver_type");
693  }
694  computePressureAndFluxes(r, sat);
695  }
696 
697  private:
699  class FaceFluxes
700  {
701  public:
702  FaceFluxes(int sz)
703  : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
704  {
705  }
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;
710  ++visited_[f_ix];
711  }
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);
718  flux = new_flux;
719  ++visited_[f_ix];
720  }
721  void resetVisited()
722  {
723  std::fill(visited_.begin(), visited_.end(), 0);
724  }
725 
726  double maxMod() const
727  {
728  return max_modification_;
729  }
730  private:
731  std::vector<double> fluxes_;
732  std::vector<int> visited_;
733  double max_modification_;
734 
735  };
736 
737  public:
747  {
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_;
753 
754  FaceFluxes face_fluxes(pgrid_->numberOfFaces());
755  // First pass: compute projected fluxes.
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()];
761  if (f->boundary()) {
762  if (ppartner_dof_.empty()) {
763  continue;
764  }
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);
769  }
770  } else {
771  face_fluxes.put(flux, f_ix);
772  }
773  }
774  }
775  face_fluxes.resetVisited();
776  // Second pass: set all fluxes to the projected ones.
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()];
782  if (f->boundary()) {
783  if (ppartner_dof_.empty()) {
784  continue;
785  }
786  int partner_f_ix = ppartner_dof_[f_ix];
787  if (partner_f_ix != -1) {
788  face_fluxes.get(flux, f_ix);
789  double dummy = flux;
790  face_fluxes.get(dummy, partner_f_ix);
791  assert(dummy == flux);
792  }
793  } else {
794  face_fluxes.get(flux, f_ix);
795  }
796  }
797  }
798  return face_fluxes.maxMod();
799  }
800 
801 
810  typedef const FlowSolution& SolutionType;
811 
821  {
822  return flowSolution_;
823  }
824 
825 
840  template<typename charT, class traits>
841  void printStats(std::basic_ostream<charT,traits>& os)
842  {
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';
847 
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, " "));
852  os << "\b]\n";
853 
854  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
855  os << "cell faces =\n";
856  for (int i = 0; i < cf.size(); ++i)
857  {
858  os << "\t[" << i << "] -> [";
859  std::copy(cf[i].begin(), cf[i].end(),
860  std::ostream_iterator<int>(os, ","));
861  os << "\b]\n";
862  }
863  }
864 
865 
887  void printSystem(const std::string& prefix)
888  {
889  writeMatrixToMatlab(S_, prefix + "-mat.dat");
890 
891  std::string rhsfile(prefix + "-rhs.dat");
892  std::ofstream rhs(rhsfile.c_str());
893  rhs.precision(15);
894  rhs.setf(std::ios::scientific | std::ios::showpos);
895  std::copy(rhs_.begin(), rhs_.end(),
896  std::ostream_iterator<VectorBlockType>(rhs, "\n"));
897  }
898 
899  private:
900  typedef std::pair<int,int> DofID;
901  typedef std::unordered_map<int,DofID> BdryIdMapType;
902  typedef BdryIdMapType::const_iterator BdryIdMapIterator;
903 
904  const GridInterface* pgrid_;
905  BdryIdMapType bdry_id_map_;
906  std::vector<int> ppartner_dof_;
907 
908  InnerProduct<GridInterface, RockInterface> ip_;
909 
910  // ----------------------------------------------------------------
911  bool cleared_state_;
912  int max_ncf_;
913  int num_internal_faces_;
914  int total_num_faces_;
915 
916  // ----------------------------------------------------------------
917  std::vector<Scalar> L_, g_;
918  Opm::SparseTable<Scalar> F_ ;
919 
920  // ----------------------------------------------------------------
921  // Actual, assembled system of linear equations
922  typedef Dune::FieldVector<Scalar, 1 > VectorBlockType;
923  typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
924 
925  Dune::BCRSMatrix <MatrixBlockType> S_; // System matrix
926  Dune::BlockVector<VectorBlockType> rhs_; // System RHS
927  Dune::BlockVector<VectorBlockType> soln_; // System solution (contact pressure)
928  bool matrix_structure_valid_;
929  bool do_regularization_;
930 
931  // ----------------------------------------------------------------
932  // Physical quantities (derived)
933  FlowSolution flowSolution_;
934 
935 
936  // ----------------------------------------------------------------
937  void enumerateDof(const GridInterface& g, const BCInterface& bc)
938  // ----------------------------------------------------------------
939  {
940  enumerateGridDof(g);
941  enumerateBCDof(g, bc);
942 
943  pgrid_ = &g;
944  cleared_state_ = false;
945  }
946 
947  // ----------------------------------------------------------------
948  void enumerateGridDof(const GridInterface& g)
949  // ----------------------------------------------------------------
950  {
951  typedef typename GridInterface::CellIterator CI;
952  typedef typename CI ::FaceIterator FI;
953 
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 ;
958 
959  // Allocate cell structures.
960  std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
961 
962  std::vector<int>& cell = flowSolution_.cellno_;
963 
964  // First pass: enumerate internal faces.
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));
970 
971  cell[c0] = cellno;
972 
973  int& ncf = num_cf[c0];
974 
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));
979 
980  if (cell[c1] == -1) {
981  // Previously undiscovered internal face.
982  faces.push_back(c1);
983  }
984  }
985  ++ncf;
986  }
987 
988  fpos.push_back(int(faces.size()));
989  max_ncf_ = std::max(max_ncf_, ncf);
990  tot_ncf += ncf;
991  tot_ncf2 += ncf * ncf;
992  }
993  assert (cellno == nc);
994 
995  total_num_faces_ = num_internal_faces_ = int(faces.size());
996 
997  ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
998  F_ .reserve(nc, tot_ncf);
999 
1000  flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1001  flowSolution_.outflux_ .reserve(nc, tot_ncf);
1002 
1003  Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1004 
1005  // Avoid (most) allocation(s) inside 'c' loop.
1006  std::vector<int> l2g; l2g .reserve(max_ncf_);
1007  std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1008 
1009  // Second pass: build cell-to-face mapping, including boundary.
1010  typedef std::vector<int>::iterator VII;
1011  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1012  const int c0 = c->index();
1013 
1014  assert ((0 <= c0 ) && ( c0 < nc) &&
1015  (0 <= cell[c0]) && (cell[c0] < nc));
1016 
1017  const int ncf = num_cf[cell[c0]];
1018  l2g .resize(ncf , 0 );
1019  F_alloc .resize(ncf , Scalar(0.0));
1020 
1021  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1022  if (f->boundary()) {
1023  // External, not counted before. Add new face...
1024  l2g[f->localIndex()] = total_num_faces_++;
1025  } else {
1026  // Internal face. Need to determine during
1027  // traversal of which cell we discovered this
1028  // face first, and extract the face number
1029  // from the 'faces' table range of that cell.
1030 
1031  // Note: std::find() below is potentially
1032  // *VERY* expensive (e.g., large number of
1033  // seeks in moderately sized data in case of
1034  // faulted cells).
1035 
1036  const int c1 = f->neighbourCellIndex();
1037  assert ((0 <= c1 ) && ( c1 < nc) &&
1038  (0 <= cell[c1]) && (cell[c1] < nc));
1039 
1040  int t = c0, seek = c1;
1041  if (cell[seek] < cell[t])
1042  std::swap(t, seek);
1043 
1044  int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1045 
1046  VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1047  assert(p != faces.begin() + e);
1048 
1049  l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1050  }
1051  }
1052 
1053  cf.appendRow (l2g .begin(), l2g .end());
1054  F_.appendRow (F_alloc.begin(), F_alloc.end());
1055 
1056  flowSolution_.outflux_
1057  .appendRow(F_alloc.begin(), F_alloc.end());
1058  }
1059  }
1060 
1061 
1062  // ----------------------------------------------------------------
1063  void enumerateBCDof(const GridInterface& g, const BCInterface& bc)
1064  // ----------------------------------------------------------------
1065  {
1066  typedef typename GridInterface::CellIterator CI;
1067  typedef typename CI ::FaceIterator FI;
1068 
1069  const std::vector<int>& cell = flowSolution_.cellno_;
1070  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1071 
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));
1079  }
1080  }
1081  }
1082 
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()];
1090 
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];
1095 
1096  ppartner_dof_[dof1] = dof2;
1097  ppartner_dof_[dof2] = dof1;
1098  }
1099  }
1100  }
1101  }
1102  }
1103 
1104 
1105 
1106  // ----------------------------------------------------------------
1107  void allocateConnections(const BCInterface& bc)
1108  // ----------------------------------------------------------------
1109  {
1110  // You must call enumerateDof() prior to allocateConnections()
1111  assert(!cleared_state_);
1112 
1113  assert (!matrix_structure_valid_);
1114 
1115  // Clear any residual data, prepare for assembling structure.
1116  S_.setSize(total_num_faces_, total_num_faces_);
1117  S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1118 
1119  // Compute row sizes
1120  for (int f = 0; f < total_num_faces_; ++f) {
1121  S_.setrowsize(f, 1);
1122  }
1123 
1124  allocateGridConnections();
1125  allocateBCConnections(bc);
1126 
1127  S_.endrowsizes();
1128 
1129  rhs_ .resize(total_num_faces_);
1130  soln_.resize(total_num_faces_);
1131  }
1132 
1133 
1134  // ----------------------------------------------------------------
1135  void allocateGridConnections()
1136  // ----------------------------------------------------------------
1137  {
1138  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1139 
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);
1144  }
1145  }
1146  }
1147 
1148 
1149  // ----------------------------------------------------------------
1150  void allocateBCConnections(const BCInterface& bc)
1151  // ----------------------------------------------------------------
1152  {
1153  // Include system connections for periodic boundary
1154  // conditions, if any. We make an arbitrary choice in
1155  // that the face/degree-of-freedom with the minimum
1156  // numerical id of the two periodic partners represents
1157  // the coupling. Suppose <i_p> is this minimum dof-id.
1158  // We then need to introduce a *symmetric* coupling to
1159  // <i_p> to each of the dof's of the cell *NOT* containing
1160  // <i_p>. This choice is implemented in the following
1161  // loop by introducing couplings only when dof1 (self) is
1162  // less than dof2 (other).
1163  //
1164  // See also: setBCConnections() and addCellContrib().
1165  //
1166  typedef typename GridInterface::CellIterator CI;
1167  typedef typename CI ::FaceIterator FI;
1168 
1169  const std::vector<int>& cell = flowSolution_.cellno_;
1170  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1171 
1172  if (!bdry_id_map_.empty()) {
1173  // At least one periodic BC. Allocate corresponding
1174  // connections.
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()) {
1178  // dof-id of self
1179  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1180 
1181  // dof-id of other
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];
1187 
1188  if (dof1 < dof2) {
1189  // Allocate storage for the actual
1190  // couplings.
1191  //
1192  const int ndof = cf.rowSize(c2);
1193  S_.incrementrowsize(dof1, ndof); // self->other
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);
1199  }
1200  S_.incrementrowsize(ii, 1); // other->self
1201  }
1202  }
1203  }
1204  }
1205  }
1206  }
1207  }
1208 
1209 
1210 
1211  // ----------------------------------------------------------------
1212  void setConnections(const BCInterface& bc)
1213  // ----------------------------------------------------------------
1214  {
1215  setGridConnections();
1216  setBCConnections(bc);
1217 
1218  S_.endindices();
1219 
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_);
1224 
1225  matrix_structure_valid_ = true;
1226  }
1227 
1228 
1229  // ----------------------------------------------------------------
1230  void setGridConnections()
1231  // ----------------------------------------------------------------
1232  {
1233  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1234 
1235  // Compute actual connections (the non-zero structure).
1236  for (int c = 0; c < cf.size(); ++c) {
1237  auto fb = cf[c].begin(), fe = cf[c].end();
1238 
1239  for (auto i = fb; i != fe; ++i) {
1240  for (auto j = fb; j != fe; ++j) {
1241  S_.addindex(*i, *j);
1242  }
1243  }
1244  }
1245  }
1246 
1247 
1248  // ----------------------------------------------------------------
1249  void setBCConnections(const BCInterface& bc)
1250  // ----------------------------------------------------------------
1251  {
1252  // Include system connections for periodic boundary
1253  // conditions, if any. We make an arbitrary choice in
1254  // that the face/degree-of-freedom with the minimum
1255  // numerical id of the two periodic partners represents
1256  // the coupling. Suppose <i_p> is this minimum dof-id.
1257  // We then need to introduce a *symmetric* coupling to
1258  // <i_p> to each of the dof's of the cell *NOT* containing
1259  // <i_p>. This choice is implemented in the following
1260  // loop by introducing couplings only when dof1 (self) is
1261  // less than dof2 (other).
1262  //
1263  // See also: allocateBCConnections() and addCellContrib().
1264  //
1265  typedef typename GridInterface::CellIterator CI;
1266  typedef typename CI ::FaceIterator FI;
1267 
1268  const std::vector<int>& cell = flowSolution_.cellno_;
1269  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1270 
1271  if (!bdry_id_map_.empty()) {
1272  // At least one periodic BC. Assign periodic
1273  // connections.
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()) {
1277  // dof-id of self
1278  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1279 
1280  // dof-id of other
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];
1286 
1287  if (dof1 < dof2) {
1288  // Assign actual couplings.
1289  //
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)) {
1295  ii = pp;
1296  }
1297  S_.addindex(dof1, ii); // self->other
1298  S_.addindex(ii, dof1); // other->self
1299  S_.addindex(dof2, ii);
1300  S_.addindex(ii, dof2);
1301  }
1302  }
1303  }
1304  }
1305  }
1306  }
1307  }
1308 
1309 
1310 
1311  // ----------------------------------------------------------------
1312  template<class FluidInterface>
1313  void assembleDynamic(const FluidInterface& fl ,
1314  const std::vector<double>& sat,
1315  const BCInterface& bc ,
1316  const std::vector<double>& src)
1317  // ----------------------------------------------------------------
1318  {
1319  typedef typename GridInterface::CellIterator CI;
1320 
1321  const std::vector<int>& cell = flowSolution_.cellno_;
1322  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1323 
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_);
1328 
1329  std::vector<FaceType> facetype (max_ncf_);
1330  std::vector<Scalar> condval (max_ncf_);
1331  std::vector<int> ppartner (max_ncf_);
1332 
1333  // Clear residual data
1334  S_ = 0.0;
1335  rhs_ = 0.0;
1336 
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));
1340 
1341  // We will have to regularize resulting system if there
1342  // are no prescribed pressures (i.e., Dirichlet BC's).
1343  do_regularization_ = true;
1344 
1345  // Assemble dynamic contributions for each cell
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();
1350 
1351  rhs .resize(nf);
1352  gflux.resize(nf);
1353 
1354  setExternalContrib(c, c0, bc, src[ci], rhs,
1355  facetype, condval, ppartner);
1356 
1357  ip_.computeDynamicParams(c, fl, sat);
1358 
1359  SharedFortranMatrix S(nf, nf, &data_store[0]);
1360  ip_.getInverseMatrix(c, S);
1361 
1362  std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1363  ip_.gravityFlux(c, gflux);
1364 
1365  ImmutableFortranMatrix one(nf, 1, &e[0]);
1366  buildCellContrib(c0, one, gflux, S, rhs);
1367 
1368  addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1369  }
1370  }
1371 
1372 
1373 
1374  // ----------------------------------------------------------------
1375  void solveLinearSystem(double residual_tolerance, int verbosity_level, int maxit)
1376  // ----------------------------------------------------------------
1377  {
1378  // Adapted from DuMux...
1379  Scalar residTol = residual_tolerance;
1380 
1381  typedef Dune::BCRSMatrix <MatrixBlockType> MatrixT;
1382  typedef Dune::BlockVector<VectorBlockType> VectorT;
1383  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1384 
1385  // Regularize the matrix (only for pure Neumann problems...)
1386  if (do_regularization_) {
1387  S_[0][0] *= 2;
1388  }
1389  Adapter opS(S_);
1390 
1391  // Construct preconditioner.
1392 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1393  Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1394 #else
1395  Dune::SeqILU0<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1396 #endif
1397 
1398  // Construct solver for system of linear equations.
1399  Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
1400  (maxit>0)?maxit:S_.N(), verbosity_level);
1401 
1402  Dune::InverseOperatorResult result;
1403  soln_ = 0.0;
1404 
1405  // Solve system of linear equations to recover
1406  // face/contact pressure values (soln_).
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');
1411  }
1412  }
1413 
1414 
1415 
1416  // ------------------ AMG typedefs --------------------
1417 
1418  // Representation types for linear system.
1419  typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1420  typedef Dune::BlockVector<VectorBlockType> Vector;
1421  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1422 
1423  // AMG specific types.
1424  // Old: FIRST_DIAGONAL 1, SYMMETRIC 1, SMOOTHER_ILU 1, ANISOTROPIC_3D 0
1425  // SPE10: FIRST_DIAGONAL 0, SYMMETRIC 1, SMOOTHER_ILU 0, ANISOTROPIC_3D 1
1426 #ifndef FIRST_DIAGONAL
1427 #define FIRST_DIAGONAL 1
1428 #endif
1429 #ifndef SYMMETRIC
1430 #define SYMMETRIC 1
1431 #endif
1432 #ifndef SMOOTHER_ILU
1433 #define SMOOTHER_ILU 1
1434 #endif
1435 #ifndef SMOOTHER_BGS
1436 #define SMOOTHER_BGS 0
1437 #endif
1438 #ifndef ANISOTROPIC_3D
1439 #define ANISOTROPIC_3D 0
1440 #endif
1441 
1442 #if FIRST_DIAGONAL
1443  typedef Dune::Amg::FirstDiagonal CouplingMetric;
1444 #else
1445  typedef Dune::Amg::RowSum CouplingMetric;
1446 #endif
1447 
1448 #if SYMMETRIC
1449  typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1450 #else
1451  typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1452 #endif
1453 
1454 #if SMOOTHER_BGS
1455  typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1456 #else
1457 #if SMOOTHER_ILU
1458 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1459  typedef Dune::SeqILU<Matrix,Vector,Vector> Smoother;
1460 #else
1461  typedef Dune::SeqILU0<Matrix,Vector,Vector> Smoother;
1462 #endif
1463 #else
1464  typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1465 #endif
1466 #endif
1467  typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1468 
1469 
1470  // --------- storing the AMG operator and preconditioner --------
1471  std::unique_ptr<Operator> opS_;
1472  typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1473  std::unique_ptr<PrecondBase> precond_;
1474 
1475 
1476  // ----------------------------------------------------------------
1477  void solveLinearSystemAMG(double residual_tolerance, int verbosity_level,
1478  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1479  // ----------------------------------------------------------------
1480  {
1481  typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1482  Precond;
1483 
1484  // Adapted from upscaling.cc by Arne Rekdal, 2009
1485  Scalar residTol = residual_tolerance;
1486 
1487  if (!same_matrix) {
1488  // Regularize the matrix (only for pure Neumann problems...)
1489  if (do_regularization_) {
1490  S_[0][0] *= 2;
1491  }
1492  opS_.reset(new Operator(S_));
1493 
1494  // Construct preconditioner.
1495  double relax = 1;
1496  typename Precond::SmootherArgs smootherArgs;
1497  smootherArgs.relaxationFactor = relax;
1498 #if SMOOTHER_BGS
1499  smootherArgs.overlap = Precond::SmootherArgs::none;
1500  smootherArgs.onthefly = false;
1501 #endif
1502  Criterion criterion;
1503  criterion.setDebugLevel(verbosity_level);
1504 #if ANISOTROPIC_3D
1505  criterion.setDefaultValuesAnisotropic(3, 2);
1506 #endif
1507  criterion.setProlongationDampingFactor(prolong_factor);
1508  criterion.setBeta(1e-10);
1509  criterion.setNoPreSmoothSteps(smooth_steps);
1510  criterion.setNoPostSmoothSteps(smooth_steps);
1511  criterion.setGamma(1); // V-cycle; this is the default
1512  precond_.reset(new Precond(*opS_, criterion, smootherArgs));
1513  }
1514  // Construct solver for system of linear equations.
1515  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1516 
1517  Dune::InverseOperatorResult result;
1518  soln_ = 0.0;
1519  // Adapt initial guess such Dirichlet boundary conditions are
1520  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
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)
1527  isDirichlet=false;
1528  if(isDirichlet)
1529  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1530  }
1531  // Solve system of linear equations to recover
1532  // face/contact pressure values (soln_).
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');
1537  }
1538 
1539  }
1540 
1541 
1542  // ----------------------------------------------------------------
1543  void solveLinearSystemFastAMG(double residual_tolerance, int verbosity_level,
1544  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1545  // ----------------------------------------------------------------
1546  {
1547  typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1548 
1549  // Adapted from upscaling.cc by Arne Rekdal, 2009
1550  Scalar residTol = residual_tolerance;
1551 
1552  if (!same_matrix) {
1553  // Regularize the matrix (only for pure Neumann problems...)
1554  if (do_regularization_) {
1555  S_[0][0] *= 2;
1556  }
1557  opS_.reset(new Operator(S_));
1558 
1559  // Construct preconditioner.
1560  typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
1561 
1562  typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
1563  Crit criterion;
1564  criterion.setDebugLevel(verbosity_level);
1565 #if ANISOTROPIC_3D
1566  criterion.setDefaultValuesAnisotropic(3, 2);
1567 #endif
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));
1575  }
1576  // Construct solver for system of linear equations.
1577  Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1578 
1579  Dune::InverseOperatorResult result;
1580  soln_ = 0.0;
1581 
1582  // Adapt initial guess such Dirichlet boundary conditions are
1583  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
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)
1590  isDirichlet=false;
1591  if(isDirichlet)
1592  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1593  }
1594  // Solve system of linear equations to recover
1595  // face/contact pressure values (soln_).
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');
1600  }
1601 
1602  }
1603 
1604  // ----------------------------------------------------------------
1605  void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level,
1606  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1607  // ----------------------------------------------------------------
1608  {
1609 
1610  typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1611  // Adapted from upscaling.cc by Arne Rekdal, 2009
1612  Scalar residTol = residual_tolerance;
1613  if (!same_matrix) {
1614  // Regularize the matrix (only for pure Neumann problems...)
1615  if (do_regularization_) {
1616  S_[0][0] *= 2;
1617  }
1618  opS_.reset(new Operator(S_));
1619 
1620  // Construct preconditioner.
1621  double relax = 1;
1622  typename Precond::SmootherArgs smootherArgs;
1623  smootherArgs.relaxationFactor = relax;
1624 #if SMOOTHER_BGS
1625  smootherArgs.overlap = Precond::SmootherArgs::none;
1626  smootherArgs.onthefly = false;
1627 #endif
1628  Criterion criterion;
1629  criterion.setDebugLevel(verbosity_level);
1630 #if ANISOTROPIC_3D
1631  criterion.setDefaultValuesAnisotropic(3, 2);
1632 #endif
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));
1639  }
1640  // Construct solver for system of linear equations.
1641  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1642 
1643  Dune::InverseOperatorResult result;
1644  soln_ = 0.0;
1645  // Adapt initial guess such Dirichlet boundary conditions are
1646  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
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)
1653  isDirichlet=false;
1654  if(isDirichlet)
1655  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1656  }
1657  // Solve system of linear equations to recover
1658  // face/contact pressure values (soln_).
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');
1663  }
1664 
1665  }
1666 
1667 
1668 
1669  // ----------------------------------------------------------------
1670  template<class FluidInterface>
1671  void computePressureAndFluxes(const FluidInterface& r ,
1672  const std::vector<double>& sat)
1673  // ----------------------------------------------------------------
1674  {
1675  typedef typename GridInterface::CellIterator CI;
1676 
1677  const std::vector<int>& cell = flowSolution_.cellno_;
1678  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1679 
1680  std::vector<Scalar>& p = flowSolution_.pressure_;
1681  Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1682 
1683  //std::vector<double> mob(FluidInterface::NumberOfPhases);
1684  std::vector<double> pi (max_ncf_);
1685  std::vector<double> gflux(max_ncf_);
1686  std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1687 
1688  // Assemble dynamic contributions for each cell
1689  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1690  const int c0 = cell[c->index()];
1691  const int nf = cf.rowSize(c0);
1692 
1693  pi .resize(nf);
1694  gflux.resize(nf);
1695 
1696  // Extract contact pressures for cell 'c'.
1697  for (int i = 0; i < nf; ++i) {
1698  pi[i] = soln_[cf[c0][i]];
1699  }
1700 
1701  // Compute cell pressure in cell 'c'.
1702  p[c0] = (g_[c0] +
1703  std::inner_product(F_[c0].begin(), F_[c0].end(),
1704  pi.begin(), 0.0)) / L_[c0];
1705 
1706  std::transform(pi.begin(), pi.end(),
1707  pi.begin(),
1708  [&p, c0](const double& input) { return p[c0] - input; });
1709 
1710  // Recover fluxes from local system
1711  // Bv = Bv_g + Cp - D\pi
1712  //
1713  // 1) Solve system Bv = Cp - D\pi
1714  //
1715  ip_.computeDynamicParams(c, r, sat);
1716 
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]);
1720 
1721  // 2) Add gravity flux contributions (v <- v + v_g)
1722  //
1723  ip_.gravityFlux(c, gflux);
1724  std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1725  v[c0].begin(),
1726  std::plus<Scalar>());
1727  }
1728  }
1729 
1730 
1731 
1732 
1733  // ----------------------------------------------------------------
1734  void setExternalContrib(const typename GridInterface::CellIterator c,
1735  const int c0, const BCInterface& bc,
1736  const double src,
1737  std::vector<Scalar>& rhs,
1738  std::vector<FaceType>& facetype,
1739  std::vector<double>& condval,
1740  std::vector<int>& ppartner)
1741  // ----------------------------------------------------------------
1742  {
1743  typedef typename GridInterface::CellIterator::FaceIterator FI;
1744 
1745  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1746 
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 );
1751 
1752  g_[c0] = src;
1753 
1754  int k = 0;
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());
1766 
1767  facetype[k] = Periodic;
1768  condval[k] = bcond.pressureDifference();
1769  ppartner[k] = cf[j->second.first][j->second.second];
1770  } else {
1771  assert (bcond.isNeumann());
1772  facetype[k] = Neumann;
1773  rhs[k] = bcond.outflux();
1774  }
1775  }
1776  }
1777  }
1778 
1779 
1780 
1781 
1782  // ----------------------------------------------------------------
1783  void buildCellContrib(const int c ,
1784  const ImmutableFortranMatrix& one ,
1785  const std::vector<Scalar>& gflux,
1786  SharedFortranMatrix& S ,
1787  std::vector<Scalar>& rhs)
1788  // ----------------------------------------------------------------
1789  {
1790  // Ft <- B^{-t} * ones([size(S,2),1])
1791  SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1792  matMulAdd_TN(Scalar(1.0), S, one, Scalar(0.0), Ft);
1793 
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));
1796 
1797  // rhs <- v_g - rhs (== v_g - h)
1798  std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1799  rhs.begin(),
1800  std::minus<Scalar>());
1801 
1802  // rhs <- rhs + g_[c]/L_[c]*F
1803  std::transform(rhs.begin(), rhs.end(), Ft.data(),
1804  rhs.begin(),
1805  axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1806 
1807  // S <- S - F'*F/L_c
1808  symmetricUpdate(-Scalar(1.0)/L_[c], Ft, Scalar(1.0), S);
1809  }
1810 
1811 
1812 
1813  // ----------------------------------------------------------------
1815  template<class L2G>
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,
1821  const L2G& l2g)
1822  // ----------------------------------------------------------------
1823  {
1824  int r = 0;
1825  for (auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1826  // Indirection for periodic BC handling.
1827  int ii = *i;
1828 
1829  switch (facetype[r]) {
1830  case Dirichlet:
1831  // Pressure is already known. Assemble trivial
1832  // equation of the form: a*x = a*p where 'p' is
1833  // the known pressure value (i.e., condval[r]).
1834  //
1835  S_ [ii][ii] = S(r,r);
1836  rhs_[ii] = S(r,r) * condval[r];
1837  continue;
1838  case Periodic:
1839  // Periodic boundary condition. Contact pressures
1840  // linked by relations of the form
1841  //
1842  // a*(x0 - x1) = a*pd
1843  //
1844  // where 'pd' is the known pressure difference
1845  // x0-x1 (condval[r]). Preserve matrix symmetry
1846  // by assembling both of the equations
1847  //
1848  // a*(x0-x1) = a* pd, and (1)
1849  // a*(x1-x0) = a*(-pd) (2)
1850  //
1851  assert ((0 <= ppartner[r]) && (ppartner[r] < int(rhs_.size())));
1852  assert (ii != ppartner[r]);
1853 
1854  {
1855  const double a = S(r,r), b = a * condval[r];
1856 
1857  // Equation (1)
1858  S_ [ ii][ ii] += a;
1859  S_ [ ii][ppartner[r]] -= a;
1860  rhs_[ ii] += b;
1861 
1862  // Equation (2)
1863  S_ [ppartner[r]][ ii] -= a;
1864  S_ [ppartner[r]][ppartner[r]] += a;
1865  rhs_[ppartner[r]] -= b;
1866  }
1867 
1868  ii = std::min(ii, ppartner[r]);
1869 
1870  // INTENTIONAL FALL-THROUGH!
1871  // IOW: Don't insert <break;> here!
1872  //
1873  default:
1874  int c = 0;
1875  for (auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1876  // Indirection for periodic BC handling.
1877  int jj = *j;
1878 
1879  if (facetype[c] == Dirichlet) {
1880  rhs_[ii] -= S(r,c) * condval[c];
1881  continue;
1882  }
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];
1888  jj = ppartner[c];
1889  }
1890  }
1891  S_[ii][jj] += S(r,c);
1892  }
1893  break;
1894  }
1895  rhs_[ii] += rhs[r];
1896  }
1897  }
1898  };
1899 } // namespace Opm
1900 
1901 #endif // OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
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