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
83namespace 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::IncompFlowSolverHybrid
Solve mixed formulation of incompressible flow modelled by Darcy's law.
Definition: IncompFlowSolverHybrid.hpp:366
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
Opm::IncompFlowSolverHybrid::getSolution
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:828
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::IncompFlowSolverHybrid::SolutionType
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:818
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::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
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
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
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
Definition: BlackoilFluid.hpp:32
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::ImmutableFortranMatrix
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix
Definition: Matrix.hpp:591
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::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
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