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;
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  }
992  assert (cellno == nc);
993 
994  total_num_faces_ = num_internal_faces_ = int(faces.size());
995 
996  ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
997  F_ .reserve(nc, tot_ncf);
998 
999  flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1000  flowSolution_.outflux_ .reserve(nc, tot_ncf);
1001 
1002  Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1003 
1004  // Avoid (most) allocation(s) inside 'c' loop.
1005  std::vector<int> l2g; l2g .reserve(max_ncf_);
1006  std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1007 
1008  // Second pass: build cell-to-face mapping, including boundary.
1009  typedef std::vector<int>::iterator VII;
1010  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1011  const int c0 = c->index();
1012 
1013  assert ((0 <= c0 ) && ( c0 < nc) &&
1014  (0 <= cell[c0]) && (cell[c0] < nc));
1015 
1016  const int ncf = num_cf[cell[c0]];
1017  l2g .resize(ncf , 0 );
1018  F_alloc .resize(ncf , Scalar(0.0));
1019 
1020  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1021  if (f->boundary()) {
1022  // External, not counted before. Add new face...
1023  l2g[f->localIndex()] = total_num_faces_++;
1024  } else {
1025  // Internal face. Need to determine during
1026  // traversal of which cell we discovered this
1027  // face first, and extract the face number
1028  // from the 'faces' table range of that cell.
1029 
1030  // Note: std::find() below is potentially
1031  // *VERY* expensive (e.g., large number of
1032  // seeks in moderately sized data in case of
1033  // faulted cells).
1034 
1035  const int c1 = f->neighbourCellIndex();
1036  assert ((0 <= c1 ) && ( c1 < nc) &&
1037  (0 <= cell[c1]) && (cell[c1] < nc));
1038 
1039  int t = c0, seek = c1;
1040  if (cell[seek] < cell[t])
1041  std::swap(t, seek);
1042 
1043  int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1044 
1045  VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1046  assert(p != faces.begin() + e);
1047 
1048  l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1049  }
1050  }
1051 
1052  cf.appendRow (l2g .begin(), l2g .end());
1053  F_.appendRow (F_alloc.begin(), F_alloc.end());
1054 
1055  flowSolution_.outflux_
1056  .appendRow(F_alloc.begin(), F_alloc.end());
1057  }
1058  }
1059 
1060 
1061  // ----------------------------------------------------------------
1062  void enumerateBCDof(const GridInterface& g, const BCInterface& bc)
1063  // ----------------------------------------------------------------
1064  {
1065  typedef typename GridInterface::CellIterator CI;
1066  typedef typename CI ::FaceIterator FI;
1067 
1068  const std::vector<int>& cell = flowSolution_.cellno_;
1069  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1070 
1071  bdry_id_map_.clear();
1072  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1073  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1074  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1075  const int bid = f->boundaryId();
1076  DofID dof(cell[c->index()], f->localIndex());
1077  bdry_id_map_.insert(std::make_pair(bid, dof));
1078  }
1079  }
1080  }
1081 
1082  ppartner_dof_.clear();
1083  if (!bdry_id_map_.empty()) {
1084  ppartner_dof_.assign(total_num_faces_, -1);
1085  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1086  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1087  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1088  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1089 
1090  BdryIdMapIterator j =
1091  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1092  assert (j != bdry_id_map_.end());
1093  const int dof2 = cf[j->second.first][j->second.second];
1094 
1095  ppartner_dof_[dof1] = dof2;
1096  ppartner_dof_[dof2] = dof1;
1097  }
1098  }
1099  }
1100  }
1101  }
1102 
1103 
1104 
1105  // ----------------------------------------------------------------
1106  void allocateConnections(const BCInterface& bc)
1107  // ----------------------------------------------------------------
1108  {
1109  // You must call enumerateDof() prior to allocateConnections()
1110  assert(!cleared_state_);
1111 
1112  assert (!matrix_structure_valid_);
1113 
1114  // Clear any residual data, prepare for assembling structure.
1115  S_.setSize(total_num_faces_, total_num_faces_);
1116  S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1117 
1118  // Compute row sizes
1119  for (int f = 0; f < total_num_faces_; ++f) {
1120  S_.setrowsize(f, 1);
1121  }
1122 
1123  allocateGridConnections();
1124  allocateBCConnections(bc);
1125 
1126  S_.endrowsizes();
1127 
1128  rhs_ .resize(total_num_faces_);
1129  soln_.resize(total_num_faces_);
1130  }
1131 
1132 
1133  // ----------------------------------------------------------------
1134  void allocateGridConnections()
1135  // ----------------------------------------------------------------
1136  {
1137  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1138 
1139  for (int c = 0; c < cf.size(); ++c) {
1140  const int nf = cf[c].size();
1141  for (auto& f : cf[c]) {
1142  S_.incrementrowsize(f, nf - 1);
1143  }
1144  }
1145  }
1146 
1147 
1148  // ----------------------------------------------------------------
1149  void allocateBCConnections(const BCInterface& bc)
1150  // ----------------------------------------------------------------
1151  {
1152  // Include system connections for periodic boundary
1153  // conditions, if any. We make an arbitrary choice in
1154  // that the face/degree-of-freedom with the minimum
1155  // numerical id of the two periodic partners represents
1156  // the coupling. Suppose <i_p> is this minimum dof-id.
1157  // We then need to introduce a *symmetric* coupling to
1158  // <i_p> to each of the dof's of the cell *NOT* containing
1159  // <i_p>. This choice is implemented in the following
1160  // loop by introducing couplings only when dof1 (self) is
1161  // less than dof2 (other).
1162  //
1163  // See also: setBCConnections() and addCellContrib().
1164  //
1165  typedef typename GridInterface::CellIterator CI;
1166  typedef typename CI ::FaceIterator FI;
1167 
1168  const std::vector<int>& cell = flowSolution_.cellno_;
1169  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1170 
1171  if (!bdry_id_map_.empty()) {
1172  // At least one periodic BC. Allocate corresponding
1173  // connections.
1174  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1175  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1176  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1177  // dof-id of self
1178  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1179 
1180  // dof-id of other
1181  BdryIdMapIterator j =
1182  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1183  assert (j != bdry_id_map_.end());
1184  const int c2 = j->second.first;
1185  const int dof2 = cf[c2][j->second.second];
1186 
1187  if (dof1 < dof2) {
1188  // Allocate storage for the actual
1189  // couplings.
1190  //
1191  const int ndof = cf.rowSize(c2);
1192  S_.incrementrowsize(dof1, ndof); // self->other
1193  for (int dof = 0; dof < ndof; ++dof) {
1194  int ii = cf[c2][dof];
1195  int pp = ppartner_dof_[ii];
1196  if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1197  S_.incrementrowsize(pp, 1);
1198  }
1199  S_.incrementrowsize(ii, 1); // other->self
1200  }
1201  }
1202  }
1203  }
1204  }
1205  }
1206  }
1207 
1208 
1209 
1210  // ----------------------------------------------------------------
1211  void setConnections(const BCInterface& bc)
1212  // ----------------------------------------------------------------
1213  {
1214  setGridConnections();
1215  setBCConnections(bc);
1216 
1217  S_.endindices();
1218 
1219  const int nc = pgrid_->numberOfCells();
1220  std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1221  std::vector<Scalar>(nc).swap(g_);
1222  std::vector<Scalar>(nc).swap(L_);
1223 
1224  matrix_structure_valid_ = true;
1225  }
1226 
1227 
1228  // ----------------------------------------------------------------
1229  void setGridConnections()
1230  // ----------------------------------------------------------------
1231  {
1232  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1233 
1234  // Compute actual connections (the non-zero structure).
1235  for (int c = 0; c < cf.size(); ++c) {
1236  auto fb = cf[c].begin(), fe = cf[c].end();
1237 
1238  for (auto i = fb; i != fe; ++i) {
1239  for (auto j = fb; j != fe; ++j) {
1240  S_.addindex(*i, *j);
1241  }
1242  }
1243  }
1244  }
1245 
1246 
1247  // ----------------------------------------------------------------
1248  void setBCConnections(const BCInterface& bc)
1249  // ----------------------------------------------------------------
1250  {
1251  // Include system connections for periodic boundary
1252  // conditions, if any. We make an arbitrary choice in
1253  // that the face/degree-of-freedom with the minimum
1254  // numerical id of the two periodic partners represents
1255  // the coupling. Suppose <i_p> is this minimum dof-id.
1256  // We then need to introduce a *symmetric* coupling to
1257  // <i_p> to each of the dof's of the cell *NOT* containing
1258  // <i_p>. This choice is implemented in the following
1259  // loop by introducing couplings only when dof1 (self) is
1260  // less than dof2 (other).
1261  //
1262  // See also: allocateBCConnections() and addCellContrib().
1263  //
1264  typedef typename GridInterface::CellIterator CI;
1265  typedef typename CI ::FaceIterator FI;
1266 
1267  const std::vector<int>& cell = flowSolution_.cellno_;
1268  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1269 
1270  if (!bdry_id_map_.empty()) {
1271  // At least one periodic BC. Assign periodic
1272  // connections.
1273  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1274  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1275  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1276  // dof-id of self
1277  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1278 
1279  // dof-id of other
1280  BdryIdMapIterator j =
1281  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1282  assert (j != bdry_id_map_.end());
1283  const int c2 = j->second.first;
1284  const int dof2 = cf[c2][j->second.second];
1285 
1286  if (dof1 < dof2) {
1287  // Assign actual couplings.
1288  //
1289  const int ndof = cf.rowSize(c2);
1290  for (int dof = 0; dof < ndof; ++dof) {
1291  int ii = cf[c2][dof];
1292  int pp = ppartner_dof_[ii];
1293  if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1294  ii = pp;
1295  }
1296  S_.addindex(dof1, ii); // self->other
1297  S_.addindex(ii, dof1); // other->self
1298  S_.addindex(dof2, ii);
1299  S_.addindex(ii, dof2);
1300  }
1301  }
1302  }
1303  }
1304  }
1305  }
1306  }
1307 
1308 
1309 
1310  // ----------------------------------------------------------------
1311  template<class FluidInterface>
1312  void assembleDynamic(const FluidInterface& fl ,
1313  const std::vector<double>& sat,
1314  const BCInterface& bc ,
1315  const std::vector<double>& src)
1316  // ----------------------------------------------------------------
1317  {
1318  typedef typename GridInterface::CellIterator CI;
1319 
1320  const std::vector<int>& cell = flowSolution_.cellno_;
1321  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1322 
1323  std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1324  std::vector<Scalar> e (max_ncf_);
1325  std::vector<Scalar> rhs (max_ncf_);
1326  std::vector<Scalar> gflux (max_ncf_);
1327 
1328  std::vector<FaceType> facetype (max_ncf_);
1329  std::vector<Scalar> condval (max_ncf_);
1330  std::vector<int> ppartner (max_ncf_);
1331 
1332  // Clear residual data
1333  S_ = 0.0;
1334  rhs_ = 0.0;
1335 
1336  std::fill(g_.begin(), g_.end(), Scalar(0.0));
1337  std::fill(L_.begin(), L_.end(), Scalar(0.0));
1338  std::fill(e .begin(), e .end(), Scalar(1.0));
1339 
1340  // We will have to regularize resulting system if there
1341  // are no prescribed pressures (i.e., Dirichlet BC's).
1342  do_regularization_ = true;
1343 
1344  // Assemble dynamic contributions for each cell
1345  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1346  const int ci = c->index();
1347  const int c0 = cell[ci]; assert (c0 < cf.size());
1348  const int nf = cf[c0].size();
1349 
1350  rhs .resize(nf);
1351  gflux.resize(nf);
1352 
1353  setExternalContrib(c, c0, bc, src[ci], rhs,
1354  facetype, condval, ppartner);
1355 
1356  ip_.computeDynamicParams(c, fl, sat);
1357 
1358  SharedFortranMatrix S(nf, nf, &data_store[0]);
1359  ip_.getInverseMatrix(c, S);
1360 
1361  std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1362  ip_.gravityFlux(c, gflux);
1363 
1364  ImmutableFortranMatrix one(nf, 1, &e[0]);
1365  buildCellContrib(c0, one, gflux, S, rhs);
1366 
1367  addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1368  }
1369  }
1370 
1371 
1372 
1373  // ----------------------------------------------------------------
1374  void solveLinearSystem(double residual_tolerance, int verbosity_level, int maxit)
1375  // ----------------------------------------------------------------
1376  {
1377  // Adapted from DuMux...
1378  Scalar residTol = residual_tolerance;
1379 
1380  typedef Dune::BCRSMatrix <MatrixBlockType> MatrixT;
1381  typedef Dune::BlockVector<VectorBlockType> VectorT;
1382  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1383 
1384  // Regularize the matrix (only for pure Neumann problems...)
1385  if (do_regularization_) {
1386  S_[0][0] *= 2;
1387  }
1388  Adapter opS(S_);
1389 
1390  // Construct preconditioner.
1391 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1392  Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1393 #else
1394  Dune::SeqILU0<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1395 #endif
1396 
1397  // Construct solver for system of linear equations.
1398  Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
1399  (maxit>0)?maxit:S_.N(), verbosity_level);
1400 
1401  Dune::InverseOperatorResult result;
1402  soln_ = 0.0;
1403 
1404  // Solve system of linear equations to recover
1405  // face/contact pressure values (soln_).
1406  linsolve.apply(soln_, rhs_, result);
1407  if (!result.converged) {
1408  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1409  << "Residual reduction achieved is " << result.reduction << '\n');
1410  }
1411  }
1412 
1413 
1414 
1415  // ------------------ AMG typedefs --------------------
1416 
1417  // Representation types for linear system.
1418  typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1419  typedef Dune::BlockVector<VectorBlockType> Vector;
1420  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1421 
1422  // AMG specific types.
1423  // Old: FIRST_DIAGONAL 1, SYMMETRIC 1, SMOOTHER_ILU 1, ANISOTROPIC_3D 0
1424  // SPE10: FIRST_DIAGONAL 0, SYMMETRIC 1, SMOOTHER_ILU 0, ANISOTROPIC_3D 1
1425 #ifndef FIRST_DIAGONAL
1426 #define FIRST_DIAGONAL 1
1427 #endif
1428 #ifndef SYMMETRIC
1429 #define SYMMETRIC 1
1430 #endif
1431 #ifndef SMOOTHER_ILU
1432 #define SMOOTHER_ILU 1
1433 #endif
1434 #ifndef SMOOTHER_BGS
1435 #define SMOOTHER_BGS 0
1436 #endif
1437 #ifndef ANISOTROPIC_3D
1438 #define ANISOTROPIC_3D 0
1439 #endif
1440 
1441 #if FIRST_DIAGONAL
1442  typedef Dune::Amg::FirstDiagonal CouplingMetric;
1443 #else
1444  typedef Dune::Amg::RowSum CouplingMetric;
1445 #endif
1446 
1447 #if SYMMETRIC
1448  typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1449 #else
1450  typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1451 #endif
1452 
1453 #if SMOOTHER_BGS
1454  typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1455 #else
1456 #if SMOOTHER_ILU
1457 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
1458  typedef Dune::SeqILU<Matrix,Vector,Vector> Smoother;
1459 #else
1460  typedef Dune::SeqILU0<Matrix,Vector,Vector> Smoother;
1461 #endif
1462 #else
1463  typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1464 #endif
1465 #endif
1466  typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1467 
1468 
1469  // --------- storing the AMG operator and preconditioner --------
1470  std::unique_ptr<Operator> opS_;
1471  typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1472  std::unique_ptr<PrecondBase> precond_;
1473 
1474 
1475  // ----------------------------------------------------------------
1476  void solveLinearSystemAMG(double residual_tolerance, int verbosity_level,
1477  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1478  // ----------------------------------------------------------------
1479  {
1480  typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1481  Precond;
1482 
1483  // Adapted from upscaling.cc by Arne Rekdal, 2009
1484  Scalar residTol = residual_tolerance;
1485 
1486  if (!same_matrix) {
1487  // Regularize the matrix (only for pure Neumann problems...)
1488  if (do_regularization_) {
1489  S_[0][0] *= 2;
1490  }
1491  opS_.reset(new Operator(S_));
1492 
1493  // Construct preconditioner.
1494  double relax = 1;
1495  typename Precond::SmootherArgs smootherArgs;
1496  smootherArgs.relaxationFactor = relax;
1497 #if SMOOTHER_BGS
1498  smootherArgs.overlap = Precond::SmootherArgs::none;
1499  smootherArgs.onthefly = false;
1500 #endif
1501  Criterion criterion;
1502  criterion.setDebugLevel(verbosity_level);
1503 #if ANISOTROPIC_3D
1504  criterion.setDefaultValuesAnisotropic(3, 2);
1505 #endif
1506  criterion.setProlongationDampingFactor(prolong_factor);
1507  criterion.setBeta(1e-10);
1508  criterion.setNoPreSmoothSteps(smooth_steps);
1509  criterion.setNoPostSmoothSteps(smooth_steps);
1510  criterion.setGamma(1); // V-cycle; this is the default
1511  precond_.reset(new Precond(*opS_, criterion, smootherArgs));
1512  }
1513  // Construct solver for system of linear equations.
1514  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1515 
1516  Dune::InverseOperatorResult result;
1517  soln_ = 0.0;
1518  // Adapt initial guess such Dirichlet boundary conditions are
1519  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1520  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1521  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1522  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1523  bool isDirichlet=true;
1524  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1525  if(ci.index()!=ri.index() && *ci!=0.0)
1526  isDirichlet=false;
1527  if(isDirichlet)
1528  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1529  }
1530  // Solve system of linear equations to recover
1531  // face/contact pressure values (soln_).
1532  linsolve.apply(soln_, rhs_, result);
1533  if (!result.converged) {
1534  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1535  << "Residual reduction achieved is " << result.reduction << '\n');
1536  }
1537 
1538  }
1539 
1540 
1541  // ----------------------------------------------------------------
1542  void solveLinearSystemFastAMG(double residual_tolerance, int verbosity_level,
1543  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1544  // ----------------------------------------------------------------
1545  {
1546  typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1547 
1548  // Adapted from upscaling.cc by Arne Rekdal, 2009
1549  Scalar residTol = residual_tolerance;
1550 
1551  if (!same_matrix) {
1552  // Regularize the matrix (only for pure Neumann problems...)
1553  if (do_regularization_) {
1554  S_[0][0] *= 2;
1555  }
1556  opS_.reset(new Operator(S_));
1557 
1558  // Construct preconditioner.
1559  typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
1560 
1561  typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
1562  Crit criterion;
1563  criterion.setDebugLevel(verbosity_level);
1564 #if ANISOTROPIC_3D
1565  criterion.setDefaultValuesAnisotropic(3, 2);
1566 #endif
1567  criterion.setProlongationDampingFactor(prolong_factor);
1568  criterion.setBeta(1e-10);
1569  Dune::Amg::Parameters parms;
1570  parms.setDebugLevel(verbosity_level);
1571  parms.setNoPreSmoothSteps(smooth_steps);
1572  parms.setNoPostSmoothSteps(smooth_steps);
1573  precond_.reset(new Precond(*opS_, criterion, parms));
1574  }
1575  // Construct solver for system of linear equations.
1576  Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1577 
1578  Dune::InverseOperatorResult result;
1579  soln_ = 0.0;
1580 
1581  // Adapt initial guess such Dirichlet boundary conditions are
1582  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1583  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1584  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1585  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1586  bool isDirichlet=true;
1587  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1588  if(ci.index()!=ri.index() && *ci!=0.0)
1589  isDirichlet=false;
1590  if(isDirichlet)
1591  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1592  }
1593  // Solve system of linear equations to recover
1594  // face/contact pressure values (soln_).
1595  linsolve.apply(soln_, rhs_, result);
1596  if (!result.converged) {
1597  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1598  << "Residual reduction achieved is " << result.reduction << '\n');
1599  }
1600 
1601  }
1602 
1603  // ----------------------------------------------------------------
1604  void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level,
1605  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1606  // ----------------------------------------------------------------
1607  {
1608 
1609  typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1610  // Adapted from upscaling.cc by Arne Rekdal, 2009
1611  Scalar residTol = residual_tolerance;
1612  if (!same_matrix) {
1613  // Regularize the matrix (only for pure Neumann problems...)
1614  if (do_regularization_) {
1615  S_[0][0] *= 2;
1616  }
1617  opS_.reset(new Operator(S_));
1618 
1619  // Construct preconditioner.
1620  double relax = 1;
1621  typename Precond::SmootherArgs smootherArgs;
1622  smootherArgs.relaxationFactor = relax;
1623 #if SMOOTHER_BGS
1624  smootherArgs.overlap = Precond::SmootherArgs::none;
1625  smootherArgs.onthefly = false;
1626 #endif
1627  Criterion criterion;
1628  criterion.setDebugLevel(verbosity_level);
1629 #if ANISOTROPIC_3D
1630  criterion.setDefaultValuesAnisotropic(3, 2);
1631 #endif
1632  criterion.setProlongationDampingFactor(prolong_factor);
1633  criterion.setBeta(1e-10);
1634  criterion.setNoPreSmoothSteps(smooth_steps);
1635  criterion.setNoPostSmoothSteps(smooth_steps);
1636  criterion.setGamma(2);
1637  precond_.reset(new Precond(*opS_, criterion, smootherArgs));
1638  }
1639  // Construct solver for system of linear equations.
1640  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1641 
1642  Dune::InverseOperatorResult result;
1643  soln_ = 0.0;
1644  // Adapt initial guess such Dirichlet boundary conditions are
1645  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1646  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1647  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1648  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1649  bool isDirichlet=true;
1650  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1651  if(ci.index()!=ri.index() && *ci!=0.0)
1652  isDirichlet=false;
1653  if(isDirichlet)
1654  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1655  }
1656  // Solve system of linear equations to recover
1657  // face/contact pressure values (soln_).
1658  linsolve.apply(soln_, rhs_, result);
1659  if (!result.converged) {
1660  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1661  << "Residual reduction achieved is " << result.reduction << '\n');
1662  }
1663 
1664  }
1665 
1666 
1667 
1668  // ----------------------------------------------------------------
1669  template<class FluidInterface>
1670  void computePressureAndFluxes(const FluidInterface& r ,
1671  const std::vector<double>& sat)
1672  // ----------------------------------------------------------------
1673  {
1674  typedef typename GridInterface::CellIterator CI;
1675 
1676  const std::vector<int>& cell = flowSolution_.cellno_;
1677  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1678 
1679  std::vector<Scalar>& p = flowSolution_.pressure_;
1680  Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1681 
1682  //std::vector<double> mob(FluidInterface::NumberOfPhases);
1683  std::vector<double> pi (max_ncf_);
1684  std::vector<double> gflux(max_ncf_);
1685  std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1686 
1687  // Assemble dynamic contributions for each cell
1688  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1689  const int c0 = cell[c->index()];
1690  const int nf = cf.rowSize(c0);
1691 
1692  pi .resize(nf);
1693  gflux.resize(nf);
1694 
1695  // Extract contact pressures for cell 'c'.
1696  for (int i = 0; i < nf; ++i) {
1697  pi[i] = soln_[cf[c0][i]];
1698  }
1699 
1700  // Compute cell pressure in cell 'c'.
1701  p[c0] = (g_[c0] +
1702  std::inner_product(F_[c0].begin(), F_[c0].end(),
1703  pi.begin(), 0.0)) / L_[c0];
1704 
1705  std::transform(pi.begin(), pi.end(),
1706  pi.begin(),
1707  [&p, c0](const double& input) { return p[c0] - input; });
1708 
1709  // Recover fluxes from local system
1710  // Bv = Bv_g + Cp - D\pi
1711  //
1712  // 1) Solve system Bv = Cp - D\pi
1713  //
1714  ip_.computeDynamicParams(c, r, sat);
1715 
1716  SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1717  ip_.getInverseMatrix(c, Binv);
1718  vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1719 
1720  // 2) Add gravity flux contributions (v <- v + v_g)
1721  //
1722  ip_.gravityFlux(c, gflux);
1723  std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1724  v[c0].begin(),
1725  std::plus<Scalar>());
1726  }
1727  }
1728 
1729 
1730 
1731 
1732  // ----------------------------------------------------------------
1733  void setExternalContrib(const typename GridInterface::CellIterator c,
1734  const int c0, const BCInterface& bc,
1735  const double src,
1736  std::vector<Scalar>& rhs,
1737  std::vector<FaceType>& facetype,
1738  std::vector<double>& condval,
1739  std::vector<int>& ppartner)
1740  // ----------------------------------------------------------------
1741  {
1742  typedef typename GridInterface::CellIterator::FaceIterator FI;
1743 
1744  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1745 
1746  std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1747  std::fill(facetype.begin(), facetype.end(), Internal );
1748  std::fill(condval .begin(), condval .end(), Scalar(0.0));
1749  std::fill(ppartner.begin(), ppartner.end(), -1 );
1750 
1751  g_[c0] = src;
1752 
1753  int k = 0;
1754  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1755  if (f->boundary()) {
1756  const FlowBC& bcond = bc.flowCond(*f);
1757  if (bcond.isDirichlet()) {
1758  facetype[k] = Dirichlet;
1759  condval[k] = bcond.pressure();
1760  do_regularization_ = false;
1761  } else if (bcond.isPeriodic()) {
1762  BdryIdMapIterator j =
1763  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1764  assert (j != bdry_id_map_.end());
1765 
1766  facetype[k] = Periodic;
1767  condval[k] = bcond.pressureDifference();
1768  ppartner[k] = cf[j->second.first][j->second.second];
1769  } else {
1770  assert (bcond.isNeumann());
1771  facetype[k] = Neumann;
1772  rhs[k] = bcond.outflux();
1773  }
1774  }
1775  }
1776  }
1777 
1778 
1779 
1780 
1781  // ----------------------------------------------------------------
1782  void buildCellContrib(const int c ,
1783  const ImmutableFortranMatrix& one ,
1784  const std::vector<Scalar>& gflux,
1785  SharedFortranMatrix& S ,
1786  std::vector<Scalar>& rhs)
1787  // ----------------------------------------------------------------
1788  {
1789  // Ft <- B^{-t} * ones([size(S,2),1])
1790  SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1791  matMulAdd_TN(Scalar(1.0), S, one, Scalar(0.0), Ft);
1792 
1793  L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1794  g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1795 
1796  // rhs <- v_g - rhs (== v_g - h)
1797  std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1798  rhs.begin(),
1799  std::minus<Scalar>());
1800 
1801  // rhs <- rhs + g_[c]/L_[c]*F
1802  std::transform(rhs.begin(), rhs.end(), Ft.data(),
1803  rhs.begin(),
1804  axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1805 
1806  // S <- S - F'*F/L_c
1807  symmetricUpdate(-Scalar(1.0)/L_[c], Ft, Scalar(1.0), S);
1808  }
1809 
1810 
1811 
1812  // ----------------------------------------------------------------
1814  template<class L2G>
1815  void addCellContrib(const SharedFortranMatrix& S ,
1816  const std::vector<Scalar>& rhs ,
1817  const std::vector<FaceType>& facetype,
1818  const std::vector<Scalar>& condval ,
1819  const std::vector<int>& ppartner,
1820  const L2G& l2g)
1821  // ----------------------------------------------------------------
1822  {
1823  int r = 0;
1824  for (auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1825  // Indirection for periodic BC handling.
1826  int ii = *i;
1827 
1828  switch (facetype[r]) {
1829  case Dirichlet:
1830  // Pressure is already known. Assemble trivial
1831  // equation of the form: a*x = a*p where 'p' is
1832  // the known pressure value (i.e., condval[r]).
1833  //
1834  S_ [ii][ii] = S(r,r);
1835  rhs_[ii] = S(r,r) * condval[r];
1836  continue;
1837  case Periodic:
1838  // Periodic boundary condition. Contact pressures
1839  // linked by relations of the form
1840  //
1841  // a*(x0 - x1) = a*pd
1842  //
1843  // where 'pd' is the known pressure difference
1844  // x0-x1 (condval[r]). Preserve matrix symmetry
1845  // by assembling both of the equations
1846  //
1847  // a*(x0-x1) = a* pd, and (1)
1848  // a*(x1-x0) = a*(-pd) (2)
1849  //
1850  assert ((0 <= ppartner[r]) && (ppartner[r] < int(rhs_.size())));
1851  assert (ii != ppartner[r]);
1852 
1853  {
1854  const double a = S(r,r), b = a * condval[r];
1855 
1856  // Equation (1)
1857  S_ [ ii][ ii] += a;
1858  S_ [ ii][ppartner[r]] -= a;
1859  rhs_[ ii] += b;
1860 
1861  // Equation (2)
1862  S_ [ppartner[r]][ ii] -= a;
1863  S_ [ppartner[r]][ppartner[r]] += a;
1864  rhs_[ppartner[r]] -= b;
1865  }
1866 
1867  ii = std::min(ii, ppartner[r]);
1868 
1869  // INTENTIONAL FALL-THROUGH!
1870  // IOW: Don't insert <break;> here!
1871  //
1872  default:
1873  int c = 0;
1874  for (auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1875  // Indirection for periodic BC handling.
1876  int jj = *j;
1877 
1878  if (facetype[c] == Dirichlet) {
1879  rhs_[ii] -= S(r,c) * condval[c];
1880  continue;
1881  }
1882  if (facetype[c] == Periodic) {
1883  assert ((0 <= ppartner[c]) && (ppartner[c] < int(rhs_.size())));
1884  assert (jj != ppartner[c]);
1885  if (ppartner[c] < jj) {
1886  rhs_[ii] -= S(r,c) * condval[c];
1887  jj = ppartner[c];
1888  }
1889  }
1890  S_[ii][jj] += S(r,c);
1891  }
1892  break;
1893  }
1894  rhs_[ii] += rhs[r];
1895  }
1896  }
1897  };
1898 } // namespace Opm
1899 
1900 #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