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