IncompFlowSolverHybrid.hpp
Go to the documentation of this file.
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 <boost/unordered_map.hpp>
42 
43 #include <boost/bind.hpp>
44 #include <boost/scoped_ptr.hpp>
45 #include <boost/lexical_cast.hpp>
46 
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
49 #include <opm/common/ErrorMacros.hpp>
50 #include <opm/core/utility/SparseTable.hpp>
51 
52 #include <dune/istl/bvector.hh>
53 #include <dune/istl/bcrsmatrix.hh>
54 #include <dune/istl/operators.hh>
55 #include <dune/istl/io.hh>
56 
57 #include <dune/istl/overlappingschwarz.hh>
58 #include <dune/istl/schwarz.hh>
59 #include <dune/istl/preconditioners.hh>
60 #include <dune/istl/solvers.hh>
61 #include <dune/istl/owneroverlapcopy.hh>
62 #include <dune/istl/paamg/amg.hh>
63 #include <dune/common/version.hh>
64 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
65 #include <dune/istl/paamg/fastamg.hh>
66 #endif
67 #include <dune/istl/paamg/kamg.hh>
68 #include <dune/istl/paamg/pinfo.hh>
69 
70 #include <opm/porsol/common/BoundaryConditions.hpp>
71 #include <opm/porsol/common/Matrix.hpp>
72 
73 #include <algorithm>
74 #include <functional>
75 #include <map>
76 #include <numeric>
77 #include <ostream>
78 #include <stdexcept>
79 #include <utility>
80 #include <vector>
81 #include <iostream>
82 
83 namespace Opm {
84  namespace {
108  template<class GI>
109  bool topologyIsSane(const GI& g)
110  {
111  typedef typename GI::CellIterator CI;
112  typedef typename CI::FaceIterator FI;
113 
114  bool sane = g.numberOfCells() >= 0;
115 
116  for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
117  std::vector<int> n;
118 
119  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
120  if (!f->boundary()) {
121  n.push_back(f->neighbourCellIndex());
122  }
123  }
124  std::sort(n.begin(), n.end());
125 
126  sane = std::unique(n.begin(), n.end()) == n.end();
127  }
128 
129  return sane;
130  }
131 
132 
158  template<typename T>
159  class axpby : public std::binary_function<T,T,T> {
160  public:
166  axpby(const T& a, const T& b) : a_(a), b_(b) {}
167 
180  T operator()(const T& x, const T& y)
181  {
182  return a_*x + b_*y;
183  }
184  private:
185  T a_, b_;
186  };
187  }
188 
189 
362  template <class GridInterface,
363  class RockInterface,
364  class BCInterface,
365  template <class GridIF, class RockIF> class InnerProduct>
366  class IncompFlowSolverHybrid {
372  typedef typename GridInterface::Scalar Scalar;
373 
378  enum FaceType { Internal, Dirichlet, Neumann, Periodic };
379 
382  class FlowSolution {
383  public:
389  typedef typename GridInterface::Scalar Scalar;
390 
394  typedef typename GridInterface::CellIterator CI;
395 
398  typedef typename CI ::FaceIterator FI;
399 
400  friend class IncompFlowSolverHybrid;
401 
411  Scalar pressure(const CI& c) const
412  {
413  return pressure_[cellno_[c->index()]];
414  }
415 
426  Scalar outflux (const FI& f) const
427  {
428  return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
429  }
430  Scalar outflux (int hf) const
431  {
432  return outflux_.data(hf);
433  }
434  private:
435  std::vector< int > cellno_;
436  Opm::SparseTable< int > cellFaces_;
437  std::vector<Scalar> pressure_;
438  Opm::SparseTable<Scalar> outflux_;
439 
440  void clear() {
441  std::vector<int>().swap(cellno_);
442  cellFaces_.clear();
443 
444  std::vector<Scalar>().swap(pressure_);
445  outflux_.clear();
446  }
447  };
448 
449  public:
476  template<class Point>
477  void init(const GridInterface& g,
478  const RockInterface& r,
479  const Point& grav,
480  const BCInterface& bc)
481  {
482  clear();
483 
484  if (g.numberOfCells() > 0) {
485  initSystemStructure(g, bc);
486  computeInnerProducts(r, grav);
487  }
488  }
489 
490 
498  void clear()
499  {
500  pgrid_ = 0;
501  max_ncf_ = -1;
502  num_internal_faces_ = 0;
503  total_num_faces_ = 0;
504  matrix_structure_valid_ = false;
505  do_regularization_ = true; // Assume pure Neumann by default.
506 
507  bdry_id_map_.clear();
508 
509  std::vector<Scalar>().swap(L_);
510  std::vector<Scalar>().swap(g_);
511  F_.clear();
512 
513  flowSolution_.clear();
514 
515  cleared_state_ = true;
516  }
517 
518 
534  void initSystemStructure(const GridInterface& g,
535  const BCInterface& bc)
536  {
537  // You must call clear() prior to initSystemStructure()
538  assert (cleared_state_);
539 
540  assert (topologyIsSane(g));
541 
542  enumerateDof(g, bc);
543  allocateConnections(bc);
544  setConnections(bc);
545  }
546 
547 
564  template<class Point>
565  void computeInnerProducts(const RockInterface& r,
566  const Point& grav)
567  {
568  // You must call connectionsCompleted() prior to computeInnerProducts()
569  assert (matrix_structure_valid_);
570 
571  typedef typename GridInterface::CellIterator CI;
572  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
573 
574  int i = 0;
575  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
576  ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
577  }
578  }
579 
580 
657  template<class FluidInterface>
658  void solve(const FluidInterface& r ,
659  const std::vector<double>& sat,
660  const BCInterface& bc ,
661  const std::vector<double>& src,
662  double residual_tolerance = 1e-8,
663  int linsolver_verbosity = 1,
664  int linsolver_type = 1,
665  bool same_matrix = false,
666  int linsolver_maxit = 0,
667  double prolongate_factor = 1.6,
668  int smooth_steps = 1)
669  {
670  assembleDynamic(r, sat, bc, src);
671 // static int count = 0;
672 // ++count;
673 // printSystem(std::string("linsys_mimetic-") + boost::lexical_cast<std::string>(count));
674  switch (linsolver_type) {
675  case 0: // ILU0 preconditioned CG
676  solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
677  break;
678  case 1: // AMG preconditioned CG
679  solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
680  linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
681  break;
682 
683  case 2: // KAMG preconditioned CG
684  solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
685  linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
686  break;
687  case 3: // CG preconditioned with AMG that uses less memory badwidth
688 #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
689  solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
690  linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
691 #else
692  if(linsolver_verbosity)
693  std::cerr<<"Fast AMG is not available; falling back to CG preconditioned with the normal one."<<std::endl;
694  solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
695  linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
696 #endif
697  break;
698  default:
699  std::cerr << "Unknown linsolver_type: " << linsolver_type << '\n';
700  throw std::runtime_error("Unknown linsolver_type");
701  }
702  computePressureAndFluxes(r, sat);
703  }
704 
705  private:
707  class FaceFluxes
708  {
709  public:
710  FaceFluxes(int sz)
711  : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
712  {
713  }
714  void put(double flux, int f_ix) {
715  assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
716  double sign = visited_[f_ix] ? -1.0 : 1.0;
717  fluxes_[f_ix] += sign*flux;
718  ++visited_[f_ix];
719  }
720  void get(double& flux, int f_ix) {
721  assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
722  double sign = visited_[f_ix] ? -1.0 : 1.0;
723  double new_flux = 0.5*sign*fluxes_[f_ix];
724  double diff = std::fabs(flux - new_flux);
725  max_modification_ = std::max(max_modification_, diff);
726  flux = new_flux;
727  ++visited_[f_ix];
728  }
729  void resetVisited()
730  {
731  std::fill(visited_.begin(), visited_.end(), 0);
732  }
733 
734  double maxMod() const
735  {
736  return max_modification_;
737  }
738  private:
739  std::vector<double> fluxes_;
740  std::vector<int> visited_;
741  double max_modification_;
742 
743  };
744 
745  public:
754  double postProcessFluxes()
755  {
756  typedef typename GridInterface::CellIterator CI;
757  typedef typename CI ::FaceIterator FI;
758  const std::vector<int>& cell = flowSolution_.cellno_;
759  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
760  Opm::SparseTable<double>& cflux = flowSolution_.outflux_;
761 
762  FaceFluxes face_fluxes(pgrid_->numberOfFaces());
763  // First pass: compute projected fluxes.
764  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
765  const int cell_index = cell[c->index()];
766  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
767  int f_ix = cf[cell_index][f->localIndex()];
768  double flux = cflux[cell_index][f->localIndex()];
769  if (f->boundary()) {
770  if (ppartner_dof_.empty()) {
771  continue;
772  }
773  int partner_f_ix = ppartner_dof_[f_ix];
774  if (partner_f_ix != -1) {
775  face_fluxes.put(flux, f_ix);
776  face_fluxes.put(flux, partner_f_ix);
777  }
778  } else {
779  face_fluxes.put(flux, f_ix);
780  }
781  }
782  }
783  face_fluxes.resetVisited();
784  // Second pass: set all fluxes to the projected ones.
785  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
786  const int cell_index = cell[c->index()];
787  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
788  int f_ix = cf[cell_index][f->localIndex()];
789  double& flux = cflux[cell_index][f->localIndex()];
790  if (f->boundary()) {
791  if (ppartner_dof_.empty()) {
792  continue;
793  }
794  int partner_f_ix = ppartner_dof_[f_ix];
795  if (partner_f_ix != -1) {
796  face_fluxes.get(flux, f_ix);
797  double dummy = flux;
798  face_fluxes.get(dummy, partner_f_ix);
799  assert(dummy == flux);
800  }
801  } else {
802  face_fluxes.get(flux, f_ix);
803  }
804  }
805  }
806  return face_fluxes.maxMod();
807  }
808 
809 
818  typedef const FlowSolution& SolutionType;
819 
828  SolutionType getSolution()
829  {
830  return flowSolution_;
831  }
832 
833 
848  template<typename charT, class traits>
849  void printStats(std::basic_ostream<charT,traits>& os)
850  {
851  os << "IncompFlowSolverHybrid<>:\n"
852  << "\tMaximum number of cell faces = " << max_ncf_ << '\n'
853  << "\tNumber of internal faces = " << num_internal_faces_ << '\n'
854  << "\tTotal number of faces = " << total_num_faces_ << '\n';
855 
856  const std::vector<int>& cell = flowSolution_.cellno_;
857  os << "cell index map = [";
858  std::copy(cell.begin(), cell.end(),
859  std::ostream_iterator<int>(os, " "));
860  os << "\b]\n";
861 
862  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
863  os << "cell faces =\n";
864  for (int i = 0; i < cf.size(); ++i)
865  {
866  os << "\t[" << i << "] -> [";
867  std::copy(cf[i].begin(), cf[i].end(),
868  std::ostream_iterator<int>(os, ","));
869  os << "\b]\n";
870  }
871  }
872 
873 
895  void printSystem(const std::string& prefix)
896  {
897  writeMatrixToMatlab(S_, prefix + "-mat.dat");
898 
899  std::string rhsfile(prefix + "-rhs.dat");
900  std::ofstream rhs(rhsfile.c_str());
901  rhs.precision(15);
902  rhs.setf(std::ios::scientific | std::ios::showpos);
903  std::copy(rhs_.begin(), rhs_.end(),
904  std::ostream_iterator<VectorBlockType>(rhs, "\n"));
905  }
906 
907  private:
908  typedef std::pair<int,int> DofID;
909  typedef boost::unordered_map<int,DofID> BdryIdMapType;
910  typedef BdryIdMapType::const_iterator BdryIdMapIterator;
911 
912  const GridInterface* pgrid_;
913  BdryIdMapType bdry_id_map_;
914  std::vector<int> ppartner_dof_;
915 
916  InnerProduct<GridInterface, RockInterface> ip_;
917 
918  // ----------------------------------------------------------------
919  bool cleared_state_;
920  int max_ncf_;
921  int num_internal_faces_;
922  int total_num_faces_;
923 
924  // ----------------------------------------------------------------
925  std::vector<Scalar> L_, g_;
926  Opm::SparseTable<Scalar> F_ ;
927 
928  // ----------------------------------------------------------------
929  // Actual, assembled system of linear equations
930  typedef Dune::FieldVector<Scalar, 1 > VectorBlockType;
931  typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
932 
933  Dune::BCRSMatrix <MatrixBlockType> S_; // System matrix
934  Dune::BlockVector<VectorBlockType> rhs_; // System RHS
935  Dune::BlockVector<VectorBlockType> soln_; // System solution (contact pressure)
936  bool matrix_structure_valid_;
937  bool do_regularization_;
938 
939  // ----------------------------------------------------------------
940  // Physical quantities (derived)
941  FlowSolution flowSolution_;
942 
943 
944  // ----------------------------------------------------------------
945  void enumerateDof(const GridInterface& g, const BCInterface& bc)
946  // ----------------------------------------------------------------
947  {
948  enumerateGridDof(g);
949  enumerateBCDof(g, bc);
950 
951  pgrid_ = &g;
952  cleared_state_ = false;
953  }
954 
955  // ----------------------------------------------------------------
956  void enumerateGridDof(const GridInterface& g)
957  // ----------------------------------------------------------------
958  {
959  typedef typename GridInterface::CellIterator CI;
960  typedef typename CI ::FaceIterator FI;
961 
962  const int nc = g.numberOfCells();
963  std::vector<int> fpos ; fpos.reserve(nc + 1);
964  std::vector<int> num_cf(nc) ;
965  std::vector<int> faces ;
966 
967  // Allocate cell structures.
968  std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
969 
970  std::vector<int>& cell = flowSolution_.cellno_;
971 
972  // First pass: enumerate internal faces.
973  int cellno = 0; fpos.push_back(0);
974  int tot_ncf = 0, tot_ncf2 = 0;
975  for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++cellno) {
976  const int c0 = c->index();
977  assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
978 
979  cell[c0] = cellno;
980 
981  int& ncf = num_cf[c0];
982 
983  for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
984  if (!f->boundary()) {
985  const int c1 = f->neighbourCellIndex();
986  assert((0 <= c1) && (c1 < nc) && (c1 != c0));
987 
988  if (cell[c1] == -1) {
989  // Previously undiscovered internal face.
990  faces.push_back(c1);
991  }
992  }
993  ++ncf;
994  }
995 
996  fpos.push_back(int(faces.size()));
997  max_ncf_ = std::max(max_ncf_, ncf);
998  tot_ncf += ncf;
999  tot_ncf2 += ncf * ncf;
1000  }
1001  assert (cellno == nc);
1002 
1003  total_num_faces_ = num_internal_faces_ = int(faces.size());
1004 
1005  ip_.init(max_ncf_); ip_.reserveMatrices(num_cf);
1006  F_ .reserve(nc, tot_ncf);
1007 
1008  flowSolution_.cellFaces_.reserve(nc, tot_ncf);
1009  flowSolution_.outflux_ .reserve(nc, tot_ncf);
1010 
1011  Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1012 
1013  // Avoid (most) allocation(s) inside 'c' loop.
1014  std::vector<int> l2g; l2g .reserve(max_ncf_);
1015  std::vector<Scalar> F_alloc; F_alloc .reserve(max_ncf_);
1016 
1017  // Second pass: build cell-to-face mapping, including boundary.
1018  typedef std::vector<int>::iterator VII;
1019  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1020  const int c0 = c->index();
1021 
1022  assert ((0 <= c0 ) && ( c0 < nc) &&
1023  (0 <= cell[c0]) && (cell[c0] < nc));
1024 
1025  const int ncf = num_cf[cell[c0]];
1026  l2g .resize(ncf , 0 );
1027  F_alloc .resize(ncf , Scalar(0.0));
1028 
1029  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1030  if (f->boundary()) {
1031  // External, not counted before. Add new face...
1032  l2g[f->localIndex()] = total_num_faces_++;
1033  } else {
1034  // Internal face. Need to determine during
1035  // traversal of which cell we discovered this
1036  // face first, and extract the face number
1037  // from the 'faces' table range of that cell.
1038 
1039  // Note: std::find() below is potentially
1040  // *VERY* expensive (e.g., large number of
1041  // seeks in moderately sized data in case of
1042  // faulted cells).
1043 
1044  const int c1 = f->neighbourCellIndex();
1045  assert ((0 <= c1 ) && ( c1 < nc) &&
1046  (0 <= cell[c1]) && (cell[c1] < nc));
1047 
1048  int t = c0, seek = c1;
1049  if (cell[seek] < cell[t])
1050  std::swap(t, seek);
1051 
1052  int s = fpos[cell[t]], e = fpos[cell[t] + 1];
1053 
1054  VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
1055  assert(p != faces.begin() + e);
1056 
1057  l2g[f->localIndex()] = s + (p - (faces.begin() + s));
1058  }
1059  }
1060 
1061  cf.appendRow (l2g .begin(), l2g .end());
1062  F_.appendRow (F_alloc.begin(), F_alloc.end());
1063 
1064  flowSolution_.outflux_
1065  .appendRow(F_alloc.begin(), F_alloc.end());
1066  }
1067  }
1068 
1069 
1070  // ----------------------------------------------------------------
1071  void enumerateBCDof(const GridInterface& g, const BCInterface& bc)
1072  // ----------------------------------------------------------------
1073  {
1074  typedef typename GridInterface::CellIterator CI;
1075  typedef typename CI ::FaceIterator FI;
1076 
1077  const std::vector<int>& cell = flowSolution_.cellno_;
1078  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1079 
1080  bdry_id_map_.clear();
1081  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1082  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1083  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1084  const int bid = f->boundaryId();
1085  DofID dof(cell[c->index()], f->localIndex());
1086  bdry_id_map_.insert(std::make_pair(bid, dof));
1087  }
1088  }
1089  }
1090 
1091  ppartner_dof_.clear();
1092  if (!bdry_id_map_.empty()) {
1093  ppartner_dof_.assign(total_num_faces_, -1);
1094  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
1095  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1096  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1097  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1098 
1099  BdryIdMapIterator j =
1100  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1101  assert (j != bdry_id_map_.end());
1102  const int dof2 = cf[j->second.first][j->second.second];
1103 
1104  ppartner_dof_[dof1] = dof2;
1105  ppartner_dof_[dof2] = dof1;
1106  }
1107  }
1108  }
1109  }
1110  }
1111 
1112 
1113 
1114  // ----------------------------------------------------------------
1115  void allocateConnections(const BCInterface& bc)
1116  // ----------------------------------------------------------------
1117  {
1118  // You must call enumerateDof() prior to allocateConnections()
1119  assert(!cleared_state_);
1120 
1121  assert (!matrix_structure_valid_);
1122 
1123  // Clear any residual data, prepare for assembling structure.
1124  S_.setSize(total_num_faces_, total_num_faces_);
1125  S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
1126 
1127  // Compute row sizes
1128  for (int f = 0; f < total_num_faces_; ++f) {
1129  S_.setrowsize(f, 1);
1130  }
1131 
1132  allocateGridConnections();
1133  allocateBCConnections(bc);
1134 
1135  S_.endrowsizes();
1136 
1137  rhs_ .resize(total_num_faces_);
1138  soln_.resize(total_num_faces_);
1139  }
1140 
1141 
1142  // ----------------------------------------------------------------
1143  void allocateGridConnections()
1144  // ----------------------------------------------------------------
1145  {
1146  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1147  typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1148 
1149  for (int c = 0; c < cf.size(); ++c) {
1150  const int nf = cf[c].size();
1151  fi fb = cf[c].begin(), fe = cf[c].end();
1152 
1153  for (fi f = fb; f != fe; ++f) {
1154  S_.incrementrowsize(*f, nf - 1);
1155  }
1156  }
1157  }
1158 
1159 
1160  // ----------------------------------------------------------------
1161  void allocateBCConnections(const BCInterface& bc)
1162  // ----------------------------------------------------------------
1163  {
1164  // Include system connections for periodic boundary
1165  // conditions, if any. We make an arbitrary choice in
1166  // that the face/degree-of-freedom with the minimum
1167  // numerical id of the two periodic partners represents
1168  // the coupling. Suppose <i_p> is this minimum dof-id.
1169  // We then need to introduce a *symmetric* coupling to
1170  // <i_p> to each of the dof's of the cell *NOT* containing
1171  // <i_p>. This choice is implemented in the following
1172  // loop by introducing couplings only when dof1 (self) is
1173  // less than dof2 (other).
1174  //
1175  // See also: setBCConnections() and addCellContrib().
1176  //
1177  typedef typename GridInterface::CellIterator CI;
1178  typedef typename CI ::FaceIterator FI;
1179 
1180  const std::vector<int>& cell = flowSolution_.cellno_;
1181  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1182 
1183  if (!bdry_id_map_.empty()) {
1184  // At least one periodic BC. Allocate corresponding
1185  // connections.
1186  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1187  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1188  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1189  // dof-id of self
1190  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1191 
1192  // dof-id of other
1193  BdryIdMapIterator j =
1194  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1195  assert (j != bdry_id_map_.end());
1196  const int c2 = j->second.first;
1197  const int dof2 = cf[c2][j->second.second];
1198 
1199  if (dof1 < dof2) {
1200  // Allocate storage for the actual
1201  // couplings.
1202  //
1203  const int ndof = cf.rowSize(c2);
1204  S_.incrementrowsize(dof1, ndof); // self->other
1205  for (int dof = 0; dof < ndof; ++dof) {
1206  int ii = cf[c2][dof];
1207  int pp = ppartner_dof_[ii];
1208  if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1209  S_.incrementrowsize(pp, 1);
1210  }
1211  S_.incrementrowsize(ii, 1); // other->self
1212  }
1213  }
1214  }
1215  }
1216  }
1217  }
1218  }
1219 
1220 
1221 
1222  // ----------------------------------------------------------------
1223  void setConnections(const BCInterface& bc)
1224  // ----------------------------------------------------------------
1225  {
1226  setGridConnections();
1227  setBCConnections(bc);
1228 
1229  S_.endindices();
1230 
1231  const int nc = pgrid_->numberOfCells();
1232  std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
1233  std::vector<Scalar>(nc).swap(g_);
1234  std::vector<Scalar>(nc).swap(L_);
1235 
1236  matrix_structure_valid_ = true;
1237  }
1238 
1239 
1240  // ----------------------------------------------------------------
1241  void setGridConnections()
1242  // ----------------------------------------------------------------
1243  {
1244  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1245  typedef Opm::SparseTable<int>::row_type::const_iterator fi;
1246 
1247  // Compute actual connections (the non-zero structure).
1248  for (int c = 0; c < cf.size(); ++c) {
1249  fi fb = cf[c].begin(), fe = cf[c].end();
1250 
1251  for (fi i = fb; i != fe; ++i) {
1252  for (fi j = fb; j != fe; ++j) {
1253  S_.addindex(*i, *j);
1254  }
1255  }
1256  }
1257  }
1258 
1259 
1260  // ----------------------------------------------------------------
1261  void setBCConnections(const BCInterface& bc)
1262  // ----------------------------------------------------------------
1263  {
1264  // Include system connections for periodic boundary
1265  // conditions, if any. We make an arbitrary choice in
1266  // that the face/degree-of-freedom with the minimum
1267  // numerical id of the two periodic partners represents
1268  // the coupling. Suppose <i_p> is this minimum dof-id.
1269  // We then need to introduce a *symmetric* coupling to
1270  // <i_p> to each of the dof's of the cell *NOT* containing
1271  // <i_p>. This choice is implemented in the following
1272  // loop by introducing couplings only when dof1 (self) is
1273  // less than dof2 (other).
1274  //
1275  // See also: allocateBCConnections() and addCellContrib().
1276  //
1277  typedef typename GridInterface::CellIterator CI;
1278  typedef typename CI ::FaceIterator FI;
1279 
1280  const std::vector<int>& cell = flowSolution_.cellno_;
1281  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1282 
1283  if (!bdry_id_map_.empty()) {
1284  // At least one periodic BC. Assign periodic
1285  // connections.
1286  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1287  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
1288  if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
1289  // dof-id of self
1290  const int dof1 = cf[cell[c->index()]][f->localIndex()];
1291 
1292  // dof-id of other
1293  BdryIdMapIterator j =
1294  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1295  assert (j != bdry_id_map_.end());
1296  const int c2 = j->second.first;
1297  const int dof2 = cf[c2][j->second.second];
1298 
1299  if (dof1 < dof2) {
1300  // Assign actual couplings.
1301  //
1302  const int ndof = cf.rowSize(c2);
1303  for (int dof = 0; dof < ndof; ++dof) {
1304  int ii = cf[c2][dof];
1305  int pp = ppartner_dof_[ii];
1306  if ((pp != -1) && (pp != dof1) && (pp < ii)) {
1307  ii = pp;
1308  }
1309  S_.addindex(dof1, ii); // self->other
1310  S_.addindex(ii, dof1); // other->self
1311  S_.addindex(dof2, ii);
1312  S_.addindex(ii, dof2);
1313  }
1314  }
1315  }
1316  }
1317  }
1318  }
1319  }
1320 
1321 
1322 
1323  // ----------------------------------------------------------------
1324  template<class FluidInterface>
1325  void assembleDynamic(const FluidInterface& fl ,
1326  const std::vector<double>& sat,
1327  const BCInterface& bc ,
1328  const std::vector<double>& src)
1329  // ----------------------------------------------------------------
1330  {
1331  typedef typename GridInterface::CellIterator CI;
1332 
1333  const std::vector<int>& cell = flowSolution_.cellno_;
1334  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1335 
1336  std::vector<Scalar> data_store(max_ncf_ * max_ncf_);
1337  std::vector<Scalar> e (max_ncf_);
1338  std::vector<Scalar> rhs (max_ncf_);
1339  std::vector<Scalar> gflux (max_ncf_);
1340 
1341  std::vector<FaceType> facetype (max_ncf_);
1342  std::vector<Scalar> condval (max_ncf_);
1343  std::vector<int> ppartner (max_ncf_);
1344 
1345  // Clear residual data
1346  S_ = 0.0;
1347  rhs_ = 0.0;
1348 
1349  std::fill(g_.begin(), g_.end(), Scalar(0.0));
1350  std::fill(L_.begin(), L_.end(), Scalar(0.0));
1351  std::fill(e .begin(), e .end(), Scalar(1.0));
1352 
1353  // We will have to regularize resulting system if there
1354  // are no prescribed pressures (i.e., Dirichlet BC's).
1355  do_regularization_ = true;
1356 
1357  // Assemble dynamic contributions for each cell
1358  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1359  const int ci = c->index();
1360  const int c0 = cell[ci]; assert (c0 < cf.size());
1361  const int nf = cf[c0].size();
1362 
1363  rhs .resize(nf);
1364  gflux.resize(nf);
1365 
1366  setExternalContrib(c, c0, bc, src[ci], rhs,
1367  facetype, condval, ppartner);
1368 
1369  ip_.computeDynamicParams(c, fl, sat);
1370 
1371  SharedFortranMatrix S(nf, nf, &data_store[0]);
1372  ip_.getInverseMatrix(c, S);
1373 
1374  std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
1375  ip_.gravityFlux(c, gflux);
1376 
1377  ImmutableFortranMatrix one(nf, 1, &e[0]);
1378  buildCellContrib(c0, one, gflux, S, rhs);
1379 
1380  addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
1381  }
1382  }
1383 
1384 
1385 
1386  // ----------------------------------------------------------------
1387  void solveLinearSystem(double residual_tolerance, int verbosity_level, int maxit)
1388  // ----------------------------------------------------------------
1389  {
1390  // Adapted from DuMux...
1391  Scalar residTol = residual_tolerance;
1392 
1393  typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1394  typedef Dune::BlockVector<VectorBlockType> Vector;
1395  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
1396 
1397  // Regularize the matrix (only for pure Neumann problems...)
1398  if (do_regularization_) {
1399  S_[0][0] *= 2;
1400  }
1401  Adapter opS(S_);
1402 
1403  // Construct preconditioner.
1404  Dune::SeqILU0<Matrix,Vector,Vector> precond(S_, 1.0);
1405 
1406  // Construct solver for system of linear equations.
1407  Dune::CGSolver<Vector> linsolve(opS, precond, residTol,
1408  (maxit>0)?maxit:S_.N(), verbosity_level);
1409 
1410  Dune::InverseOperatorResult result;
1411  soln_ = 0.0;
1412 
1413  // Solve system of linear equations to recover
1414  // face/contact pressure values (soln_).
1415  linsolve.apply(soln_, rhs_, result);
1416  if (!result.converged) {
1417  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1418  << "Residual reduction achieved is " << result.reduction << '\n');
1419  }
1420  }
1421 
1422 
1423 
1424  // ------------------ AMG typedefs --------------------
1425 
1426  // Representation types for linear system.
1427  typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1428  typedef Dune::BlockVector<VectorBlockType> Vector;
1429  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1430 
1431  // AMG specific types.
1432  // Old: FIRST_DIAGONAL 1, SYMMETRIC 1, SMOOTHER_ILU 1, ANISOTROPIC_3D 0
1433  // SPE10: FIRST_DIAGONAL 0, SYMMETRIC 1, SMOOTHER_ILU 0, ANISOTROPIC_3D 1
1434 #ifndef FIRST_DIAGONAL
1435 #define FIRST_DIAGONAL 1
1436 #endif
1437 #ifndef SYMMETRIC
1438 #define SYMMETRIC 1
1439 #endif
1440 #ifndef SMOOTHER_ILU
1441 #define SMOOTHER_ILU 1
1442 #endif
1443 #ifndef SMOOTHER_BGS
1444 #define SMOOTHER_BGS 0
1445 #endif
1446 #ifndef ANISOTROPIC_3D
1447 #define ANISOTROPIC_3D 0
1448 #endif
1449 
1450 #if FIRST_DIAGONAL
1451  typedef Dune::Amg::FirstDiagonal CouplingMetric;
1452 #else
1453  typedef Dune::Amg::RowSum CouplingMetric;
1454 #endif
1455 
1456 #if SYMMETRIC
1457  typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1458 #else
1459  typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1460 #endif
1461 
1462 #if SMOOTHER_BGS
1463  typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1464 #else
1465 #if SMOOTHER_ILU
1466  typedef Dune::SeqILU0<Matrix,Vector,Vector> Smoother;
1467 #else
1468  typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1469 #endif
1470 #endif
1471  typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1472 
1473 
1474  // --------- storing the AMG operator and preconditioner --------
1475  boost::scoped_ptr<Operator> opS_;
1476  typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1477  boost::scoped_ptr<PrecondBase> precond_;
1478 
1479 
1480  // ----------------------------------------------------------------
1481  void solveLinearSystemAMG(double residual_tolerance, int verbosity_level,
1482  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1483  // ----------------------------------------------------------------
1484  {
1485  typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1486  Precond;
1487 
1488  // Adapted from upscaling.cc by Arne Rekdal, 2009
1489  Scalar residTol = residual_tolerance;
1490 
1491  if (!same_matrix) {
1492  // Regularize the matrix (only for pure Neumann problems...)
1493  if (do_regularization_) {
1494  S_[0][0] *= 2;
1495  }
1496  opS_.reset(new Operator(S_));
1497 
1498  // Construct preconditioner.
1499  double relax = 1;
1500  typename Precond::SmootherArgs smootherArgs;
1501  smootherArgs.relaxationFactor = relax;
1502 #if SMOOTHER_BGS
1503  smootherArgs.overlap = Precond::SmootherArgs::none;
1504  smootherArgs.onthefly = false;
1505 #endif
1506  Criterion criterion;
1507  criterion.setDebugLevel(verbosity_level);
1508 #if ANISOTROPIC_3D
1509  criterion.setDefaultValuesAnisotropic(3, 2);
1510 #endif
1511  criterion.setProlongationDampingFactor(prolong_factor);
1512  criterion.setBeta(1e-10);
1513  precond_.reset(new Precond(*opS_, criterion, smootherArgs,
1514  1, smooth_steps, smooth_steps));
1515  }
1516  // Construct solver for system of linear equations.
1517  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1518 
1519  Dune::InverseOperatorResult result;
1520  soln_ = 0.0;
1521  // Adapt initial guess such Dirichlet boundary conditions are
1522  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1523  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1524  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1525  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1526  bool isDirichlet=true;
1527  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1528  if(ci.index()!=ri.index() && *ci!=0.0)
1529  isDirichlet=false;
1530  if(isDirichlet)
1531  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1532  }
1533  // Solve system of linear equations to recover
1534  // face/contact pressure values (soln_).
1535  linsolve.apply(soln_, rhs_, result);
1536  if (!result.converged) {
1537  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1538  << "Residual reduction achieved is " << result.reduction << '\n');
1539  }
1540 
1541  }
1542 
1543 #if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
1544 
1545  // ----------------------------------------------------------------
1546  void solveLinearSystemFastAMG(double residual_tolerance, int verbosity_level,
1547  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1548  // ----------------------------------------------------------------
1549  {
1550  typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1551 
1552  // Adapted from upscaling.cc by Arne Rekdal, 2009
1553  Scalar residTol = residual_tolerance;
1554 
1555  if (!same_matrix) {
1556  // Regularize the matrix (only for pure Neumann problems...)
1557  if (do_regularization_) {
1558  S_[0][0] *= 2;
1559  }
1560  opS_.reset(new Operator(S_));
1561 
1562  // Construct preconditioner.
1563  typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CriterionBase;
1564 
1565  typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1566  Criterion criterion;
1567  criterion.setDebugLevel(verbosity_level);
1568 #if ANISOTROPIC_3D
1569  criterion.setDefaultValuesAnisotropic(3, 2);
1570 #endif
1571  criterion.setProlongationDampingFactor(prolong_factor);
1572  criterion.setBeta(1e-10);
1573  Dune::Amg::Parameters parms;
1574  parms.setDebugLevel(verbosity_level);
1575  parms.setNoPreSmoothSteps(smooth_steps);
1576  parms.setNoPostSmoothSteps(smooth_steps);
1577  precond_.reset(new Precond(*opS_, criterion, parms));
1578  }
1579  // Construct solver for system of linear equations.
1580  Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1581 
1582  Dune::InverseOperatorResult result;
1583  soln_ = 0.0;
1584 
1585  // Adapt initial guess such Dirichlet boundary conditions are
1586  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1587  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1588  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1589  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1590  bool isDirichlet=true;
1591  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1592  if(ci.index()!=ri.index() && *ci!=0.0)
1593  isDirichlet=false;
1594  if(isDirichlet)
1595  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1596  }
1597  // Solve system of linear equations to recover
1598  // face/contact pressure values (soln_).
1599  linsolve.apply(soln_, rhs_, result);
1600  if (!result.converged) {
1601  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1602  << "Residual reduction achieved is " << result.reduction << '\n');
1603  }
1604 
1605  }
1606 #endif
1607 
1608  // ----------------------------------------------------------------
1609  void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level,
1610  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1611  // ----------------------------------------------------------------
1612  {
1613 
1614  typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation,
1615  Dune::CGSolver<Vector> > Precond;
1616  // Adapted from upscaling.cc by Arne Rekdal, 2009
1617  Scalar residTol = residual_tolerance;
1618  if (!same_matrix) {
1619  // Regularize the matrix (only for pure Neumann problems...)
1620  if (do_regularization_) {
1621  S_[0][0] *= 2;
1622  }
1623  opS_.reset(new Operator(S_));
1624 
1625  // Construct preconditioner.
1626  double relax = 1;
1627  typename Precond::SmootherArgs smootherArgs;
1628  smootherArgs.relaxationFactor = relax;
1629 #if SMOOTHER_BGS
1630  smootherArgs.overlap = Precond::SmootherArgs::none;
1631  smootherArgs.onthefly = false;
1632 #endif
1633  Criterion criterion;
1634  criterion.setDebugLevel(verbosity_level);
1635 #if ANISOTROPIC_3D
1636  criterion.setDefaultValuesAnisotropic(3, 2);
1637 #endif
1638  criterion.setProlongationDampingFactor(prolong_factor);
1639  criterion.setBeta(1e-10);
1640  precond_.reset(new Precond(*opS_, criterion, smootherArgs, 2, smooth_steps, smooth_steps));
1641  }
1642  // Construct solver for system of linear equations.
1643  Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1644 
1645  Dune::InverseOperatorResult result;
1646  soln_ = 0.0;
1647  // Adapt initial guess such Dirichlet boundary conditions are
1648  // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1649  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1650  typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1651  for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1652  bool isDirichlet=true;
1653  for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1654  if(ci.index()!=ri.index() && *ci!=0.0)
1655  isDirichlet=false;
1656  if(isDirichlet)
1657  soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1658  }
1659  // Solve system of linear equations to recover
1660  // face/contact pressure values (soln_).
1661  linsolve.apply(soln_, rhs_, result);
1662  if (!result.converged) {
1663  OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
1664  << "Residual reduction achieved is " << result.reduction << '\n');
1665  }
1666 
1667  }
1668 
1669 
1670 
1671  // ----------------------------------------------------------------
1672  template<class FluidInterface>
1673  void computePressureAndFluxes(const FluidInterface& r ,
1674  const std::vector<double>& sat)
1675  // ----------------------------------------------------------------
1676  {
1677  typedef typename GridInterface::CellIterator CI;
1678 
1679  const std::vector<int>& cell = flowSolution_.cellno_;
1680  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1681 
1682  std::vector<Scalar>& p = flowSolution_.pressure_;
1683  Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1684 
1685  //std::vector<double> mob(FluidInterface::NumberOfPhases);
1686  std::vector<double> pi (max_ncf_);
1687  std::vector<double> gflux(max_ncf_);
1688  std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1689 
1690  // Assemble dynamic contributions for each cell
1691  for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1692  const int c0 = cell[c->index()];
1693  const int nf = cf.rowSize(c0);
1694 
1695  pi .resize(nf);
1696  gflux.resize(nf);
1697 
1698  // Extract contact pressures for cell 'c'.
1699  for (int i = 0; i < nf; ++i) {
1700  pi[i] = soln_[cf[c0][i]];
1701  }
1702 
1703  // Compute cell pressure in cell 'c'.
1704  p[c0] = (g_[c0] +
1705  std::inner_product(F_[c0].begin(), F_[c0].end(),
1706  pi.begin(), 0.0)) / L_[c0];
1707 
1708  std::transform(pi.begin(), pi.end(),
1709  pi.begin(),
1710  boost::bind(std::minus<Scalar>(), p[c0], _1));
1711 
1712  // Recover fluxes from local system
1713  // Bv = Bv_g + Cp - D\pi
1714  //
1715  // 1) Solve system Bv = Cp - D\pi
1716  //
1717  ip_.computeDynamicParams(c, r, sat);
1718 
1719  SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1720  ip_.getInverseMatrix(c, Binv);
1721  vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1722 
1723  // 2) Add gravity flux contributions (v <- v + v_g)
1724  //
1725  ip_.gravityFlux(c, gflux);
1726  std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1727  v[c0].begin(),
1728  std::plus<Scalar>());
1729  }
1730  }
1731 
1732 
1733 
1734 
1735  // ----------------------------------------------------------------
1736  void setExternalContrib(const typename GridInterface::CellIterator c,
1737  const int c0, const BCInterface& bc,
1738  const double src,
1739  std::vector<Scalar>& rhs,
1740  std::vector<FaceType>& facetype,
1741  std::vector<double>& condval,
1742  std::vector<int>& ppartner)
1743  // ----------------------------------------------------------------
1744  {
1745  typedef typename GridInterface::CellIterator::FaceIterator FI;
1746 
1747  const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1748 
1749  std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1750  std::fill(facetype.begin(), facetype.end(), Internal );
1751  std::fill(condval .begin(), condval .end(), Scalar(0.0));
1752  std::fill(ppartner.begin(), ppartner.end(), -1 );
1753 
1754  g_[c0] = src;
1755 
1756  int k = 0;
1757  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1758  if (f->boundary()) {
1759  const FlowBC& bcond = bc.flowCond(*f);
1760  if (bcond.isDirichlet()) {
1761  facetype[k] = Dirichlet;
1762  condval[k] = bcond.pressure();
1763  do_regularization_ = false;
1764  } else if (bcond.isPeriodic()) {
1765  BdryIdMapIterator j =
1766  bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1767  assert (j != bdry_id_map_.end());
1768 
1769  facetype[k] = Periodic;
1770  condval[k] = bcond.pressureDifference();
1771  ppartner[k] = cf[j->second.first][j->second.second];
1772  } else {
1773  assert (bcond.isNeumann());
1774  facetype[k] = Neumann;
1775  rhs[k] = bcond.outflux();
1776  }
1777  }
1778  }
1779  }
1780 
1781 
1782 
1783 
1784  // ----------------------------------------------------------------
1785  void buildCellContrib(const int c ,
1786  const ImmutableFortranMatrix& one ,
1787  const std::vector<Scalar>& gflux,
1788  SharedFortranMatrix& S ,
1789  std::vector<Scalar>& rhs)
1790  // ----------------------------------------------------------------
1791  {
1792  // Ft <- B^{-t} * ones([size(S,2),1])
1793  SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1794  matMulAdd_TN(Scalar(1.0), S, one, Scalar(0.0), Ft);
1795 
1796  L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1797  g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1798 
1799  // rhs <- v_g - rhs (== v_g - h)
1800  std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1801  rhs.begin(),
1802  std::minus<Scalar>());
1803 
1804  // rhs <- rhs + g_[c]/L_[c]*F
1805  std::transform(rhs.begin(), rhs.end(), Ft.data(),
1806  rhs.begin(),
1807  axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1808 
1809  // S <- S - F'*F/L_c
1810  symmetricUpdate(-Scalar(1.0)/L_[c], Ft, Scalar(1.0), S);
1811  }
1812 
1813 
1814 
1815  // ----------------------------------------------------------------
1817  template<class L2G>
1818  void addCellContrib(const SharedFortranMatrix& S ,
1819  const std::vector<Scalar>& rhs ,
1820  const std::vector<FaceType>& facetype,
1821  const std::vector<Scalar>& condval ,
1822  const std::vector<int>& ppartner,
1823  const L2G& l2g)
1824  // ----------------------------------------------------------------
1825  {
1826  typedef typename L2G::const_iterator it;
1827 
1828  int r = 0;
1829  for (it i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1830  // Indirection for periodic BC handling.
1831  int ii = *i;
1832 
1833  switch (facetype[r]) {
1834  case Dirichlet:
1835  // Pressure is already known. Assemble trivial
1836  // equation of the form: a*x = a*p where 'p' is
1837  // the known pressure value (i.e., condval[r]).
1838  //
1839  S_ [ii][ii] = S(r,r);
1840  rhs_[ii] = S(r,r) * condval[r];
1841  continue;
1842  case Periodic:
1843  // Periodic boundary condition. Contact pressures
1844  // linked by relations of the form
1845  //
1846  // a*(x0 - x1) = a*pd
1847  //
1848  // where 'pd' is the known pressure difference
1849  // x0-x1 (condval[r]). Preserve matrix symmetry
1850  // by assembling both of the equations
1851  //
1852  // a*(x0-x1) = a* pd, and (1)
1853  // a*(x1-x0) = a*(-pd) (2)
1854  //
1855  assert ((0 <= ppartner[r]) && (ppartner[r] < int(rhs_.size())));
1856  assert (ii != ppartner[r]);
1857 
1858  {
1859  const double a = S(r,r), b = a * condval[r];
1860 
1861  // Equation (1)
1862  S_ [ ii][ ii] += a;
1863  S_ [ ii][ppartner[r]] -= a;
1864  rhs_[ ii] += b;
1865 
1866  // Equation (2)
1867  S_ [ppartner[r]][ ii] -= a;
1868  S_ [ppartner[r]][ppartner[r]] += a;
1869  rhs_[ppartner[r]] -= b;
1870  }
1871 
1872  ii = std::min(ii, ppartner[r]);
1873 
1874  // INTENTIONAL FALL-THROUGH!
1875  // IOW: Don't insert <break;> here!
1876  //
1877  default:
1878  int c = 0;
1879  for (it j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1880  // Indirection for periodic BC handling.
1881  int jj = *j;
1882 
1883  if (facetype[c] == Dirichlet) {
1884  rhs_[ii] -= S(r,c) * condval[c];
1885  continue;
1886  }
1887  if (facetype[c] == Periodic) {
1888  assert ((0 <= ppartner[c]) && (ppartner[c] < int(rhs_.size())));
1889  assert (jj != ppartner[c]);
1890  if (ppartner[c] < jj) {
1891  rhs_[ii] -= S(r,c) * condval[c];
1892  jj = ppartner[c];
1893  }
1894  }
1895  S_[ii][jj] += S(r,c);
1896  }
1897  break;
1898  }
1899  rhs_[ii] += rhs[r];
1900  }
1901  }
1902  };
1903 } // namespace Opm
1904 
1905 #endif // OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
Opm::symmetricUpdate
void symmetricUpdate(const T &a1, const FullMatrix< T, StoragePolicy, FortranOrdering > &A, const T &a2, FullMatrix< T, StoragePolicy, FortranOrdering > &C)
Symmetric, rank update of symmetric matrix. Specifically, .
Definition: Matrix.hpp:830
Opm::IncompFlowSolverHybrid::printSystem
void printSystem(const std::string &prefix)
Output current system of linear equations to permanent storage in files. One file for the coefficient...
Definition: IncompFlowSolverHybrid.hpp:895
Opm::IncompFlowSolverHybrid::solve
void solve(const FluidInterface &r, const std::vector< double > &sat, const BCInterface &bc, const std::vector< double > &src, double residual_tolerance=1e-8, int linsolver_verbosity=1, int linsolver_type=1, bool same_matrix=false, int linsolver_maxit=0, double prolongate_factor=1.6, int smooth_steps=1)
Construct and solve system of linear equations for the pressure values on each interface/contact betw...
Definition: IncompFlowSolverHybrid.hpp:658
Opm
Definition: BlackoilFluid.hpp:31
Opm::ImmutableFortranMatrix
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix
Definition: Matrix.hpp:591
Opm::IncompFlowSolverHybrid::init
void init(const GridInterface &g, const RockInterface &r, const Point &grav, const BCInterface &bc)
All-in-one initialization routine. Enumerates all grid connections, allocates sufficient space...
Definition: IncompFlowSolverHybrid.hpp:477
a_
T a_
Definition: IncompFlowSolverHybrid.hpp:185
Opm::IncompFlowSolverHybrid::SolutionType
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:818
Opm::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
Opm::IncompFlowSolverHybrid::initSystemStructure
void initSystemStructure(const GridInterface &g, const BCInterface &bc)
Compute structure of coefficient matrix in final system of linear equations for this flow problem...
Definition: IncompFlowSolverHybrid.hpp:534
Opm::IncompFlowSolverHybrid::computeInnerProducts
void computeInnerProducts(const RockInterface &r, const Point &grav)
Compute static (i.e., independent of saturation) parts of the spatially varying inner products for e...
Definition: IncompFlowSolverHybrid.hpp:565
Opm::IncompFlowSolverHybrid::clear
void clear()
Clear all topologic, geometric and rock-dependent information currently held in internal data structu...
Definition: IncompFlowSolverHybrid.hpp:498
Matrix.hpp
Opm::IncompFlowSolverHybrid::postProcessFluxes
double postProcessFluxes()
Postprocess the solution fluxes. This method modifies the solution object so that out-fluxes of twin ...
Definition: IncompFlowSolverHybrid.hpp:754
BoundaryConditions.hpp
b_
T b_
Definition: IncompFlowSolverHybrid.hpp:185
Opm::IncompFlowSolverHybrid
Solve mixed formulation of incompressible flow modelled by Darcy's law ] The solver is based on a hyb...
Definition: IncompFlowSolverHybrid.hpp:366
Opm::matMulAdd_TN
void matMulAdd_TN(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix. Specificlly, .
Definition: Matrix.hpp:1253
Opm::vecMulAdd_N
void vecMulAdd_N(const T &a1, const FullMatrix< T, SP, FortranOrdering > &A, const std::vector< T > &x, const T &a2, std::vector< T > &y)
GEneral Matrix-Vector product (GAXPY operation). Specifically, .
Definition: Matrix.hpp:914
Opm::IncompFlowSolverHybrid::getSolution
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:828
Opm::IncompFlowSolverHybrid::printStats
void printStats(std::basic_ostream< charT, traits > &os)
Print statistics about the connections in the current model. This is mostly for debugging purposes an...
Definition: IncompFlowSolverHybrid.hpp:849