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