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 <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 {
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>
365 class IncompFlowSolverHybrid {
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:
746 double postProcessFluxes()
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
820 SolutionType getSolution()
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 Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
1392
1393 // Construct solver for system of linear equations.
1394 Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
1395 (maxit>0)?maxit:S_.N(), verbosity_level);
1396
1397 Dune::InverseOperatorResult result;
1398 soln_ = 0.0;
1399
1400 // Solve system of linear equations to recover
1401 // face/contact pressure values (soln_).
1402 linsolve.apply(soln_, rhs_, result);
1403 if (!result.converged) {
1404 OPM_THROW(std::runtime_error,
1405 "Linear solver failed to converge in " +
1406 std::to_string(result.iterations) + " iterations.\n"
1407 "Residual reduction achieved is " +
1408 std::to_string(result.reduction) + '\n');
1409 }
1410 }
1411
1412
1413
1414 // ------------------ AMG typedefs --------------------
1415
1416 // Representation types for linear system.
1417 typedef Dune::BCRSMatrix <MatrixBlockType> Matrix;
1418 typedef Dune::BlockVector<VectorBlockType> Vector;
1419 typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
1420
1421 // AMG specific types.
1422 // Old: FIRST_DIAGONAL 1, SYMMETRIC 1, SMOOTHER_ILU 1, ANISOTROPIC_3D 0
1423 // SPE10: FIRST_DIAGONAL 0, SYMMETRIC 1, SMOOTHER_ILU 0, ANISOTROPIC_3D 1
1424#ifndef FIRST_DIAGONAL
1425#define FIRST_DIAGONAL 1
1426#endif
1427#ifndef SYMMETRIC
1428#define SYMMETRIC 1
1429#endif
1430#ifndef SMOOTHER_ILU
1431#define SMOOTHER_ILU 1
1432#endif
1433#ifndef SMOOTHER_BGS
1434#define SMOOTHER_BGS 0
1435#endif
1436#ifndef ANISOTROPIC_3D
1437#define ANISOTROPIC_3D 0
1438#endif
1439
1440#if FIRST_DIAGONAL
1441 typedef Dune::Amg::FirstDiagonal CouplingMetric;
1442#else
1443 typedef Dune::Amg::RowSum CouplingMetric;
1444#endif
1445
1446#if SYMMETRIC
1447 typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1448#else
1449 typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
1450#endif
1451
1452#if SMOOTHER_BGS
1453 typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
1454#else
1455#if SMOOTHER_ILU
1456 typedef Dune::SeqILU<Matrix,Vector,Vector> Smoother;
1457#else
1458 typedef Dune::SeqSSOR<Matrix,Vector,Vector> Smoother;
1459#endif
1460#endif
1461 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
1462
1463
1464 // --------- storing the AMG operator and preconditioner --------
1465 std::unique_ptr<Operator> opS_;
1466 typedef Dune::Preconditioner<Vector,Vector> PrecondBase;
1467 std::unique_ptr<PrecondBase> precond_;
1468
1469
1470 // ----------------------------------------------------------------
1471 void solveLinearSystemAMG(double residual_tolerance, int verbosity_level,
1472 int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1473 // ----------------------------------------------------------------
1474 {
1475 typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
1476 Precond;
1477
1478 // Adapted from upscaling.cc by Arne Rekdal, 2009
1479 Scalar residTol = residual_tolerance;
1480
1481 if (!same_matrix) {
1482 // Regularize the matrix (only for pure Neumann problems...)
1483 if (do_regularization_) {
1484 S_[0][0] *= 2;
1485 }
1486 opS_.reset(new Operator(S_));
1487
1488 // Construct preconditioner.
1489 double relax = 1;
1490 typename Precond::SmootherArgs smootherArgs;
1491 smootherArgs.relaxationFactor = relax;
1492#if SMOOTHER_BGS
1493 smootherArgs.overlap = Precond::SmootherArgs::none;
1494 smootherArgs.onthefly = false;
1495#endif
1496 Criterion criterion;
1497 criterion.setDebugLevel(verbosity_level);
1498#if ANISOTROPIC_3D
1499 criterion.setDefaultValuesAnisotropic(3, 2);
1500#endif
1501 criterion.setProlongationDampingFactor(prolong_factor);
1502 criterion.setBeta(1e-10);
1503 criterion.setNoPreSmoothSteps(smooth_steps);
1504 criterion.setNoPostSmoothSteps(smooth_steps);
1505 criterion.setGamma(1); // V-cycle; this is the default
1506 precond_.reset(new Precond(*opS_, criterion, smootherArgs));
1507 }
1508 // Construct solver for system of linear equations.
1509 Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1510
1511 Dune::InverseOperatorResult result;
1512 soln_ = 0.0;
1513 // Adapt initial guess such Dirichlet boundary conditions are
1514 // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1515 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1516 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1517 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1518 bool isDirichlet=true;
1519 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1520 if(ci.index()!=ri.index() && *ci!=0.0)
1521 isDirichlet=false;
1522 if(isDirichlet)
1523 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1524 }
1525 // Solve system of linear equations to recover
1526 // face/contact pressure values (soln_).
1527 linsolve.apply(soln_, rhs_, result);
1528 if (!result.converged) {
1529 OPM_THROW(std::runtime_error,
1530 "Linear solver failed to converge in " +
1531 std::to_string(result.iterations) + " iterations.\n"
1532 "Residual reduction achieved is " +
1533 std::to_string(result.reduction) + '\n');
1534 }
1535
1536 }
1537
1538
1539 // ----------------------------------------------------------------
1540 void solveLinearSystemFastAMG(double residual_tolerance, int verbosity_level,
1541 int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1542 // ----------------------------------------------------------------
1543 {
1544 typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
1545
1546 // Adapted from upscaling.cc by Arne Rekdal, 2009
1547 Scalar residTol = residual_tolerance;
1548
1549 if (!same_matrix) {
1550 // Regularize the matrix (only for pure Neumann problems...)
1551 if (do_regularization_) {
1552 S_[0][0] *= 2;
1553 }
1554 opS_.reset(new Operator(S_));
1555
1556 // Construct preconditioner.
1557 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
1558
1559 typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
1560 Crit criterion;
1561 criterion.setDebugLevel(verbosity_level);
1562#if ANISOTROPIC_3D
1563 criterion.setDefaultValuesAnisotropic(3, 2);
1564#endif
1565 criterion.setProlongationDampingFactor(prolong_factor);
1566 criterion.setBeta(1e-10);
1567 Dune::Amg::Parameters parms;
1568 parms.setDebugLevel(verbosity_level);
1569 parms.setNoPreSmoothSteps(smooth_steps);
1570 parms.setNoPostSmoothSteps(smooth_steps);
1571 precond_.reset(new Precond(*opS_, criterion, parms));
1572 }
1573 // Construct solver for system of linear equations.
1574 Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1575
1576 Dune::InverseOperatorResult result;
1577 soln_ = 0.0;
1578
1579 // Adapt initial guess such Dirichlet boundary conditions are
1580 // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1581 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1582 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1583 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1584 bool isDirichlet=true;
1585 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1586 if(ci.index()!=ri.index() && *ci!=0.0)
1587 isDirichlet=false;
1588 if(isDirichlet)
1589 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1590 }
1591 // Solve system of linear equations to recover
1592 // face/contact pressure values (soln_).
1593 linsolve.apply(soln_, rhs_, result);
1594 if (!result.converged) {
1595 OPM_THROW(std::runtime_error,
1596 "Linear solver failed to converge in " +
1597 std::to_string(result.iterations) + " iterations.\n"
1598 "Residual reduction achieved is " +
1599 std::to_string(result.reduction) + '\n');
1600 }
1601
1602 }
1603
1604 // ----------------------------------------------------------------
1605 void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level,
1606 int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
1607 // ----------------------------------------------------------------
1608 {
1609
1610 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
1611 // Adapted from upscaling.cc by Arne Rekdal, 2009
1612 Scalar residTol = residual_tolerance;
1613 if (!same_matrix) {
1614 // Regularize the matrix (only for pure Neumann problems...)
1615 if (do_regularization_) {
1616 S_[0][0] *= 2;
1617 }
1618 opS_.reset(new Operator(S_));
1619
1620 // Construct preconditioner.
1621 double relax = 1;
1622 typename Precond::SmootherArgs smootherArgs;
1623 smootherArgs.relaxationFactor = relax;
1624#if SMOOTHER_BGS
1625 smootherArgs.overlap = Precond::SmootherArgs::none;
1626 smootherArgs.onthefly = false;
1627#endif
1628 Criterion criterion;
1629 criterion.setDebugLevel(verbosity_level);
1630#if ANISOTROPIC_3D
1631 criterion.setDefaultValuesAnisotropic(3, 2);
1632#endif
1633 criterion.setProlongationDampingFactor(prolong_factor);
1634 criterion.setBeta(1e-10);
1635 criterion.setNoPreSmoothSteps(smooth_steps);
1636 criterion.setNoPostSmoothSteps(smooth_steps);
1637 criterion.setGamma(2);
1638 precond_.reset(new Precond(*opS_, criterion, smootherArgs));
1639 }
1640 // Construct solver for system of linear equations.
1641 Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
1642
1643 Dune::InverseOperatorResult result;
1644 soln_ = 0.0;
1645 // Adapt initial guess such Dirichlet boundary conditions are
1646 // represented, i.e. soln_i=A_{ii}^-1 rhs_i
1647 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
1648 typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
1649 for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
1650 bool isDirichlet=true;
1651 for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
1652 if(ci.index()!=ri.index() && *ci!=0.0)
1653 isDirichlet=false;
1654 if(isDirichlet)
1655 soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
1656 }
1657 // Solve system of linear equations to recover
1658 // face/contact pressure values (soln_).
1659 linsolve.apply(soln_, rhs_, result);
1660 if (!result.converged) {
1661 OPM_THROW(std::runtime_error,
1662 "Linear solver failed to converge in " +
1663 std::to_string(result.iterations) + " iterations.\n"
1664 "Residual reduction achieved is " +
1665 std::to_string(result.reduction) + '\n');
1666 }
1667
1668 }
1669
1670
1671
1672 // ----------------------------------------------------------------
1673 template<class FluidInterface>
1674 void computePressureAndFluxes(const FluidInterface& r ,
1675 const std::vector<double>& sat)
1676 // ----------------------------------------------------------------
1677 {
1678 typedef typename GridInterface::CellIterator CI;
1679
1680 const std::vector<int>& cell = flowSolution_.cellno_;
1681 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1682
1683 std::vector<Scalar>& p = flowSolution_.pressure_;
1684 Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
1685
1686 //std::vector<double> mob(FluidInterface::NumberOfPhases);
1687 std::vector<double> pi (max_ncf_);
1688 std::vector<double> gflux(max_ncf_);
1689 std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
1690
1691 // Assemble dynamic contributions for each cell
1692 for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
1693 const int c0 = cell[c->index()];
1694 const int nf = cf.rowSize(c0);
1695
1696 pi .resize(nf);
1697 gflux.resize(nf);
1698
1699 // Extract contact pressures for cell 'c'.
1700 for (int i = 0; i < nf; ++i) {
1701 pi[i] = soln_[cf[c0][i]];
1702 }
1703
1704 // Compute cell pressure in cell 'c'.
1705 p[c0] = (g_[c0] +
1706 std::inner_product(F_[c0].begin(), F_[c0].end(),
1707 pi.begin(), 0.0)) / L_[c0];
1708
1709 std::transform(pi.begin(), pi.end(),
1710 pi.begin(),
1711 [&p, c0](const double& input) { return p[c0] - input; });
1712
1713 // Recover fluxes from local system
1714 // Bv = Bv_g + Cp - D\pi
1715 //
1716 // 1) Solve system Bv = Cp - D\pi
1717 //
1718 ip_.computeDynamicParams(c, r, sat);
1719
1720 SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
1721 ip_.getInverseMatrix(c, Binv);
1722 vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
1723
1724 // 2) Add gravity flux contributions (v <- v + v_g)
1725 //
1726 ip_.gravityFlux(c, gflux);
1727 std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
1728 v[c0].begin(),
1729 std::plus<Scalar>());
1730 }
1731 }
1732
1733
1734
1735
1736 // ----------------------------------------------------------------
1737 void setExternalContrib(const typename GridInterface::CellIterator c,
1738 const int c0, const BCInterface& bc,
1739 const double src,
1740 std::vector<Scalar>& rhs,
1741 std::vector<FaceType>& facetype,
1742 std::vector<double>& condval,
1743 std::vector<int>& ppartner)
1744 // ----------------------------------------------------------------
1745 {
1746 typedef typename GridInterface::CellIterator::FaceIterator FI;
1747
1748 const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
1749
1750 std::fill(rhs .begin(), rhs .end(), Scalar(0.0));
1751 std::fill(facetype.begin(), facetype.end(), Internal );
1752 std::fill(condval .begin(), condval .end(), Scalar(0.0));
1753 std::fill(ppartner.begin(), ppartner.end(), -1 );
1754
1755 g_[c0] = src;
1756
1757 int k = 0;
1758 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
1759 if (f->boundary()) {
1760 const FlowBC& bcond = bc.flowCond(*f);
1761 if (bcond.isDirichlet()) {
1762 facetype[k] = Dirichlet;
1763 condval[k] = bcond.pressure();
1764 do_regularization_ = false;
1765 } else if (bcond.isPeriodic()) {
1766 BdryIdMapIterator j =
1767 bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
1768 assert (j != bdry_id_map_.end());
1769
1770 facetype[k] = Periodic;
1771 condval[k] = bcond.pressureDifference();
1772 ppartner[k] = cf[j->second.first][j->second.second];
1773 } else {
1774 assert (bcond.isNeumann());
1775 facetype[k] = Neumann;
1776 rhs[k] = bcond.outflux();
1777 }
1778 }
1779 }
1780 }
1781
1782
1783
1784
1785 // ----------------------------------------------------------------
1786 void buildCellContrib(const int c ,
1787 const ImmutableFortranMatrix& one ,
1788 const std::vector<Scalar>& gflux,
1789 SharedFortranMatrix& S ,
1790 std::vector<Scalar>& rhs)
1791 // ----------------------------------------------------------------
1792 {
1793 // Ft <- B^{-t} * ones([size(S,2),1])
1794 SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
1795 matMulAdd_TN(Scalar(1.0), S, one, Scalar(0.0), Ft);
1796
1797 L_[c] = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
1798 g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
1799
1800 // rhs <- v_g - rhs (== v_g - h)
1801 std::transform(gflux.begin(), gflux.end(), rhs.begin(),
1802 rhs.begin(),
1803 std::minus<Scalar>());
1804
1805 // rhs <- rhs + g_[c]/L_[c]*F
1806 std::transform(rhs.begin(), rhs.end(), Ft.data(),
1807 rhs.begin(),
1808 axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
1809
1810 // S <- S - F'*F/L_c
1811 symmetricUpdate(-Scalar(1.0)/L_[c], Ft, Scalar(1.0), S);
1812 }
1813
1814
1815
1816 // ----------------------------------------------------------------
1818 template<class L2G>
1819 void addCellContrib(const SharedFortranMatrix& S ,
1820 const std::vector<Scalar>& rhs ,
1821 const std::vector<FaceType>& facetype,
1822 const std::vector<Scalar>& condval ,
1823 const std::vector<int>& ppartner,
1824 const L2G& l2g)
1825 // ----------------------------------------------------------------
1826 {
1827 int r = 0;
1828 for (auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
1829 // Indirection for periodic BC handling.
1830 int ii = *i;
1831
1832 switch (facetype[r]) {
1833 case Dirichlet:
1834 // Pressure is already known. Assemble trivial
1835 // equation of the form: a*x = a*p where 'p' is
1836 // the known pressure value (i.e., condval[r]).
1837 //
1838 S_ [ii][ii] = S(r,r);
1839 rhs_[ii] = S(r,r) * condval[r];
1840 continue;
1841 case Periodic:
1842 // Periodic boundary condition. Contact pressures
1843 // linked by relations of the form
1844 //
1845 // a*(x0 - x1) = a*pd
1846 //
1847 // where 'pd' is the known pressure difference
1848 // x0-x1 (condval[r]). Preserve matrix symmetry
1849 // by assembling both of the equations
1850 //
1851 // a*(x0-x1) = a* pd, and (1)
1852 // a*(x1-x0) = a*(-pd) (2)
1853 //
1854 assert ((0 <= ppartner[r]) && (ppartner[r] < int(rhs_.size())));
1855 assert (ii != ppartner[r]);
1856
1857 {
1858 const double a = S(r,r), b = a * condval[r];
1859
1860 // Equation (1)
1861 S_ [ ii][ ii] += a;
1862 S_ [ ii][ppartner[r]] -= a;
1863 rhs_[ ii] += b;
1864
1865 // Equation (2)
1866 S_ [ppartner[r]][ ii] -= a;
1867 S_ [ppartner[r]][ppartner[r]] += a;
1868 rhs_[ppartner[r]] -= b;
1869 }
1870
1871 ii = std::min(ii, ppartner[r]);
1872
1873 // INTENTIONAL FALL-THROUGH!
1874 // IOW: Don't insert <break;> here!
1875 //
1876 default:
1877 int c = 0;
1878 for (auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
1879 // Indirection for periodic BC handling.
1880 int jj = *j;
1881
1882 if (facetype[c] == Dirichlet) {
1883 rhs_[ii] -= S(r,c) * condval[c];
1884 continue;
1885 }
1886 if (facetype[c] == Periodic) {
1887 assert ((0 <= ppartner[c]) && (ppartner[c] < int(rhs_.size())));
1888 assert (jj != ppartner[c]);
1889 if (ppartner[c] < jj) {
1890 rhs_[ii] -= S(r,c) * condval[c];
1891 jj = ppartner[c];
1892 }
1893 }
1894 S_[ii][jj] += S(r,c);
1895 }
1896 break;
1897 }
1898 rhs_[ii] += rhs[r];
1899 }
1900 }
1901 };
1902} // namespace Opm
1903
1904#endif // OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
Opm::IncompFlowSolverHybrid
Solve mixed formulation of incompressible flow modelled by Darcy's law.
Definition: IncompFlowSolverHybrid.hpp:365
Opm::IncompFlowSolverHybrid::postProcessFluxes
double postProcessFluxes()
Postprocess the solution fluxes. This method modifies the solution object so that out-fluxes of twin ...
Definition: IncompFlowSolverHybrid.hpp:746
Opm::IncompFlowSolverHybrid::getSolution
SolutionType getSolution()
Recover the solution to the problem defined by the parameters to method.
Definition: IncompFlowSolverHybrid.hpp:820
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:657
Opm::IncompFlowSolverHybrid::SolutionType
const FlowSolution & SolutionType
Type representing the solution to the problem defined by the parameters to.
Definition: IncompFlowSolverHybrid.hpp:810
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:887
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:564
Opm::IncompFlowSolverHybrid::clear
void clear()
Clear all topologic, geometric and rock-dependent information currently held in internal data structu...
Definition: IncompFlowSolverHybrid.hpp:497
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:841
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:476
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:533
Opm::Elasticity::Operator
Dune::MatrixAdapter< Matrix, Vector, Vector > Operator
A linear operator.
Definition: elasticity_preconditioners.hpp:49
Opm::Elasticity::min
min[0]
Definition: elasticity_upscale_impl.hpp:146
Opm::Elasticity::j
int j
Definition: elasticity_upscale_impl.hpp:658
Opm
Definition: ImplicitAssembly.hpp:43
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:829
Opm::ImmutableFortranMatrix
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix
Definition: Matrix.hpp:590
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:913
Opm::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:589
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:1252