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