IncompFlowSolverHybrid.hpp

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
//===========================================================================
//
// File: IncompFlowSolverHybrid.hpp
//
// Created: Tue Jun 30 10:25:40 2009
//
// Author(s): B�rd Skaflestad     <bard.skaflestad@sintef.no>
//            Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
 
/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.
 
  This file is part of The Open Reservoir Simulator Project (OpenRS).
 
  OpenRS is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OpenRS is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
#define OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER
 
#include <opm/common/ErrorMacros.hpp>
#include <opm/grid/utility/SparseTable.hpp>
#include <opm/porsol/common/BoundaryConditions.hpp>
#include <opm/porsol/common/Matrix.hpp>
 
#include <opm/common/utility/platform_dependent/disable_warnings.h>
 
#include <unordered_map>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
 
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/overlappingschwarz.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/common/version.hh>
#include <dune/istl/paamg/fastamg.hh>
#include <dune/istl/paamg/kamg.hh>
#include <dune/istl/paamg/pinfo.hh>
 
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 
#include <algorithm>
#include <functional>
#include <map>
#include <memory>
#include <numeric>
#include <ostream>
#include <stdexcept>
#include <utility>
#include <vector>
#include <iostream>
 
namespace Opm {
    namespace {
        template<class GI>
        bool topologyIsSane(const GI& g)
        {
            typedef typename GI::CellIterator CI;
            typedef typename CI::FaceIterator FI;
 
            bool sane = g.numberOfCells() >= 0;
 
            for (CI c = g.cellbegin(); sane && c != g.cellend(); ++c) {
                std::vector<int> n;
 
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    if (!f->boundary()) {
                        n.push_back(f->neighbourCellIndex());
                    }
                }
                std::sort(n.begin(), n.end());
 
                sane = std::unique(n.begin(), n.end()) == n.end();
            }
 
            return sane;
        }
 
 
        template<typename T>
        class axpby {
        public:
            axpby(const T& a, const T& b) : a_(a), b_(b) {}
 
            T operator()(const T& x, const T& y)
            {
                return a_*x + b_*y;
            }
        private:
            T a_, b_;
        };
    }
 
 
    template <class GridInterface,
              class RockInterface,
              class BCInterface,
              template <class GridIF, class RockIF> class InnerProduct>
    class IncompFlowSolverHybrid {
        typedef typename GridInterface::Scalar Scalar;
 
        enum FaceType { Internal, Dirichlet, Neumann, Periodic };
 
        class FlowSolution {
        public:
            typedef typename GridInterface::Scalar       Scalar;
 
            typedef typename GridInterface::CellIterator CI;
 
            typedef typename CI           ::FaceIterator FI;
 
            friend class IncompFlowSolverHybrid;
 
            Scalar pressure(const CI& c) const
            {
                return pressure_[cellno_[c->index()]];
            }
 
            Scalar outflux (const FI& f) const
            {
                return outflux_[cellno_[f->cellIndex()]][f->localIndex()];
            }
            Scalar outflux (int hf) const
            {
                            return outflux_.data(hf);
            }
        private:
            std::vector< int  > cellno_;
        Opm::SparseTable< int  > cellFaces_;
            std::vector<Scalar> pressure_;
        Opm::SparseTable<Scalar> outflux_;
 
            void clear() {
                std::vector<int>().swap(cellno_);
                cellFaces_.clear();
 
                std::vector<Scalar>().swap(pressure_);
                outflux_.clear();
            }
        };
 
    public:
        template<class Point>
        void init(const GridInterface&      g,
                  const RockInterface&      r,
                  const Point&              grav,
                  const BCInterface&        bc)
        {
            clear();
 
            if (g.numberOfCells() > 0) {
                initSystemStructure(g, bc);
                computeInnerProducts(r, grav);
            }
        }
 
 
        void clear()
        {
            pgrid_                  =  0;
            max_ncf_                = -1;
            num_internal_faces_     =  0;
            total_num_faces_        =  0;
            matrix_structure_valid_ = false;
            do_regularization_      = true; // Assume pure Neumann by default.
 
            bdry_id_map_.clear();
 
            std::vector<Scalar>().swap(L_);
            std::vector<Scalar>().swap(g_);
            F_.clear();
 
            flowSolution_.clear();
 
            cleared_state_ = true;
        }
 
 
        void initSystemStructure(const GridInterface& g,
                                 const BCInterface&   bc)
        {
            // You must call clear() prior to initSystemStructure()
            assert (cleared_state_);
 
            assert  (topologyIsSane(g));
 
            enumerateDof(g, bc);
            allocateConnections(bc);
            setConnections(bc);
        }
 
 
        template<class Point>
        void computeInnerProducts(const RockInterface& r,
                                  const Point& grav)
        {
            // You must call connectionsCompleted() prior to computeInnerProducts()
            assert (matrix_structure_valid_);
 
            typedef typename GridInterface::CellIterator CI;
            const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
 
            int i = 0;
            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++i) {
                ip_.buildStaticContrib(c, r, grav, cf.rowSize(i));
            }
        }
 
 
        template<class FluidInterface>
        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)
        {
            assembleDynamic(r, sat, bc, src);
//             static int count = 0;
//             ++count;
//             printSystem(std::string("linsys_mimetic-") + std::to_string(count));
            switch (linsolver_type) {
            case 0: // ILU0 preconditioned CG
              solveLinearSystem(residual_tolerance, linsolver_verbosity, linsolver_maxit);
                break;
            case 1: // AMG preconditioned CG
                solveLinearSystemAMG(residual_tolerance, linsolver_verbosity,
                     linsolver_maxit, prolongate_factor, same_matrix, smooth_steps);
                break;
 
            case 2: // KAMG preconditioned CG
                solveLinearSystemKAMG(residual_tolerance, linsolver_verbosity,
                      linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
                break;
            case 3: // CG preconditioned with AMG that uses less memory badwidth
                solveLinearSystemFastAMG(residual_tolerance, linsolver_verbosity,
                             linsolver_maxit, prolongate_factor, same_matrix,smooth_steps);
                break;
            default:
                std::cerr << "Unknown linsolver_type: " << linsolver_type << '\n';
                throw std::runtime_error("Unknown linsolver_type");
            }
            computePressureAndFluxes(r, sat);
        }
 
    private:
        class FaceFluxes
        {
        public:
            explicit FaceFluxes(int sz)
                : fluxes_(sz, 0.0), visited_(sz, 0), max_modification_(0.0)
            {
            }
            void put(double flux, int f_ix) {
                assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
                double sign = visited_[f_ix] ? -1.0 : 1.0;
                fluxes_[f_ix] += sign*flux;
                ++visited_[f_ix];
            }
            void get(double& flux, int f_ix) {
                assert(visited_[f_ix] == 0 || visited_[f_ix] == 1);
                double sign = visited_[f_ix] ? -1.0 : 1.0;
                double new_flux = 0.5*sign*fluxes_[f_ix];
                double diff = std::fabs(flux - new_flux);
                max_modification_ = std::max(max_modification_, diff);
                flux = new_flux;
                ++visited_[f_ix];
            }
            void resetVisited()
            {
                std::fill(visited_.begin(), visited_.end(), 0);
            }
 
            double maxMod() const
            {
                return max_modification_;
            }
        private:
            std::vector<double> fluxes_;
            std::vector<int> visited_;
            double max_modification_;
 
        };
 
    public:
        double postProcessFluxes()
        {
            typedef typename GridInterface::CellIterator CI;
            typedef typename CI           ::FaceIterator FI;
            const std::vector<int>& cell     = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf       = flowSolution_.cellFaces_;
        Opm::SparseTable<double>& cflux = flowSolution_.outflux_;
 
            FaceFluxes face_fluxes(pgrid_->numberOfFaces());
            // First pass: compute projected fluxes.
            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                const int cell_index = cell[c->index()];
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    int f_ix = cf[cell_index][f->localIndex()];
                    double flux = cflux[cell_index][f->localIndex()];
                    if (f->boundary()) {
                        if (ppartner_dof_.empty()) {
                            continue;
                        }
                        int partner_f_ix = ppartner_dof_[f_ix];
                        if (partner_f_ix != -1) {
                            face_fluxes.put(flux, f_ix);
                            face_fluxes.put(flux, partner_f_ix);
                        }
                    } else {
                        face_fluxes.put(flux, f_ix);
                    }
                }
            }
            face_fluxes.resetVisited();
            // Second pass: set all fluxes to the projected ones.
            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                const int cell_index = cell[c->index()];
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    int f_ix = cf[cell_index][f->localIndex()];
                    double& flux = cflux[cell_index][f->localIndex()];
                    if (f->boundary()) {
                        if (ppartner_dof_.empty()) {
                            continue;
                        }
                        int partner_f_ix = ppartner_dof_[f_ix];
                        if (partner_f_ix != -1) {
                            face_fluxes.get(flux, f_ix);
                            double dummy = flux;
                            face_fluxes.get(dummy, partner_f_ix);
                            assert(dummy == flux);
                        }
                    } else {
                        face_fluxes.get(flux, f_ix);
                    }
                }
            }
            return face_fluxes.maxMod();
        }
 
 
        typedef const FlowSolution& SolutionType;
 
        SolutionType getSolution()
        {
            return flowSolution_;
        }
 
 
        template<typename charT, class traits>
        void printStats(std::basic_ostream<charT,traits>& os)
        {
            os << "IncompFlowSolverHybrid<>:\n"
               << "\tMaximum number of cell faces = " << max_ncf_ << '\n'
               << "\tNumber of internal faces     = " << num_internal_faces_ << '\n'
               << "\tTotal number of faces        = " << total_num_faces_ << '\n';
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            os << "cell index map = [";
            std::copy(cell.begin(), cell.end(),
                      std::ostream_iterator<int>(os, " "));
            os << "\b]\n";
 
            const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
            os << "cell faces     =\n";
            for (int i = 0; i < cf.size(); ++i)
            {
                os << "\t[" << i << "] -> [";
                std::copy(cf[i].begin(), cf[i].end(),
                          std::ostream_iterator<int>(os, ","));
                os << "\b]\n";
            }
        }
 
 
        void printSystem(const std::string& prefix)
        {
            writeMatrixToMatlab(S_, prefix + "-mat.dat");
 
            std::string rhsfile(prefix + "-rhs.dat");
            std::ofstream rhs(rhsfile.c_str());
            rhs.precision(15);
            rhs.setf(std::ios::scientific | std::ios::showpos);
            std::copy(rhs_.begin(), rhs_.end(),
                      std::ostream_iterator<VectorBlockType>(rhs, "\n"));
        }
 
    private:
        typedef std::pair<int,int>                 DofID;
        typedef std::unordered_map<int,DofID> BdryIdMapType;
        typedef BdryIdMapType::const_iterator      BdryIdMapIterator;
 
        const GridInterface* pgrid_;
        BdryIdMapType        bdry_id_map_;
        std::vector<int>     ppartner_dof_;
 
        InnerProduct<GridInterface, RockInterface> ip_;
 
        // ----------------------------------------------------------------
        bool cleared_state_;
        int  max_ncf_;
        int  num_internal_faces_;
        int  total_num_faces_;
 
        // ----------------------------------------------------------------
        std::vector<Scalar> L_, g_;
    Opm::SparseTable<Scalar> F_    ;
 
        // ----------------------------------------------------------------
        // Actual, assembled system of linear equations
        typedef Dune::FieldVector<Scalar, 1   > VectorBlockType;
        typedef Dune::FieldMatrix<Scalar, 1, 1> MatrixBlockType;
 
        Dune::BCRSMatrix <MatrixBlockType>      S_;    // System matrix
        Dune::BlockVector<VectorBlockType>      rhs_;  // System RHS
        Dune::BlockVector<VectorBlockType>      soln_; // System solution (contact pressure)
        bool                              matrix_structure_valid_;
        bool                              do_regularization_;
 
        // ----------------------------------------------------------------
        // Physical quantities (derived)
        FlowSolution flowSolution_;
 
 
        // ----------------------------------------------------------------
        void enumerateDof(const GridInterface& g, const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            enumerateGridDof(g);
            enumerateBCDof(g, bc);
 
            pgrid_ = &g;
            cleared_state_ = false;
        }
 
        // ----------------------------------------------------------------
        void enumerateGridDof(const GridInterface& g)
        // ----------------------------------------------------------------
        {
            typedef typename GridInterface::CellIterator CI;
            typedef typename CI           ::FaceIterator FI;
 
            const int nc = g.numberOfCells();
            std::vector<int> fpos           ;   fpos.reserve(nc + 1);
            std::vector<int> num_cf(nc)     ;
            std::vector<int> faces          ;
 
            // Allocate cell structures.
            std::vector<int>(nc, -1).swap(flowSolution_.cellno_);
 
            std::vector<int>& cell = flowSolution_.cellno_;
 
            // First pass: enumerate internal faces.
            int cellno = 0; fpos.push_back(0);
            int tot_ncf = 0;
            for (CI c = g.cellbegin(); c != g.cellend(); ++c, ++cellno) {
                const int c0 = c->index();
                assert((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
 
                cell[c0] = cellno;
 
                int& ncf = num_cf[c0];
 
                for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
                    if (!f->boundary()) {
                        const int c1 = f->neighbourCellIndex();
                        assert((0 <= c1) && (c1 < nc) && (c1 != c0));
 
                        if (cell[c1] == -1) {
                            // Previously undiscovered internal face.
                            faces.push_back(c1);
                        }
                    }
                    ++ncf;
                }
 
                fpos.push_back(int(faces.size()));
                max_ncf_  = std::max(max_ncf_, ncf);
                tot_ncf  += ncf;
            }
            assert (cellno == nc);
 
            total_num_faces_ = num_internal_faces_ = int(faces.size());
 
            ip_.init(max_ncf_);        ip_.reserveMatrices(num_cf);
            F_ .reserve(nc, tot_ncf);
 
            flowSolution_.cellFaces_.reserve(nc, tot_ncf);
            flowSolution_.outflux_  .reserve(nc, tot_ncf);
 
        Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
 
            // Avoid (most) allocation(s) inside 'c' loop.
            std::vector<int>    l2g;        l2g       .reserve(max_ncf_);
            std::vector<Scalar> F_alloc;    F_alloc   .reserve(max_ncf_);
 
            // Second pass: build cell-to-face mapping, including boundary.
            typedef std::vector<int>::iterator VII;
            for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
                const int c0 = c->index();
 
                assert ((0 <=      c0 ) && (     c0  < nc) &&
                        (0 <= cell[c0]) && (cell[c0] < nc));
 
                const int ncf = num_cf[cell[c0]];
                l2g       .resize(ncf      ,        0   );
                F_alloc   .resize(ncf      , Scalar(0.0));
 
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    if (f->boundary()) {
                        // External, not counted before.  Add new face...
                        l2g[f->localIndex()] = total_num_faces_++;
                    } else {
                        // Internal face.  Need to determine during
                        // traversal of which cell we discovered this
                        // face first, and extract the face number
                        // from the 'faces' table range of that cell.
 
                        // Note: std::find() below is potentially
                        // *VERY* expensive (e.g., large number of
                        // seeks in moderately sized data in case of
                        // faulted cells).
 
                        const int c1 = f->neighbourCellIndex();
                        assert ((0 <=      c1 ) && (     c1  < nc) &&
                                (0 <= cell[c1]) && (cell[c1] < nc));
 
                        int t = c0, seek = c1;
                        if (cell[seek] < cell[t])
                            std::swap(t, seek);
 
                        int s = fpos[cell[t]], e = fpos[cell[t] + 1];
 
                        VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
                        assert(p != faces.begin() + e);
 
                        l2g[f->localIndex()] = s + (p - (faces.begin() + s));
                    }
                }
 
                cf.appendRow  (l2g    .begin(), l2g    .end());
                F_.appendRow  (F_alloc.begin(), F_alloc.end());
 
                flowSolution_.outflux_
                    .appendRow(F_alloc.begin(), F_alloc.end());
            }
        }
 
 
        // ----------------------------------------------------------------
        void enumerateBCDof(const GridInterface& g, const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            typedef typename GridInterface::CellIterator CI;
            typedef typename CI           ::FaceIterator FI;
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf   = flowSolution_.cellFaces_;
 
            bdry_id_map_.clear();
            for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
                        const int bid = f->boundaryId();
                        DofID dof(cell[c->index()], f->localIndex());
                        bdry_id_map_.insert(std::make_pair(bid, dof));
                    }
                }
            }
 
            ppartner_dof_.clear();
            if (!bdry_id_map_.empty()) {
                ppartner_dof_.assign(total_num_faces_, -1);
                for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
                    for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                        if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
                            const int dof1 = cf[cell[c->index()]][f->localIndex()];
 
                            BdryIdMapIterator j =
                                bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
                            assert (j != bdry_id_map_.end());
                            const int dof2 = cf[j->second.first][j->second.second];
 
                            ppartner_dof_[dof1] = dof2;
                            ppartner_dof_[dof2] = dof1;
                        }
                    }
                }
            }
        }
 
 
 
        // ----------------------------------------------------------------
        void allocateConnections(const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            // You must call enumerateDof() prior to allocateConnections()
            assert(!cleared_state_);
 
            assert  (!matrix_structure_valid_);
 
            // Clear any residual data, prepare for assembling structure.
            S_.setSize(total_num_faces_, total_num_faces_);
            S_.setBuildMode(Dune::BCRSMatrix<MatrixBlockType>::random);
 
            // Compute row sizes
            for (int f = 0; f < total_num_faces_; ++f) {
                S_.setrowsize(f, 1);
            }
 
            allocateGridConnections();
            allocateBCConnections(bc);
 
            S_.endrowsizes();
 
            rhs_ .resize(total_num_faces_);
            soln_.resize(total_num_faces_);
        }
 
 
        // ----------------------------------------------------------------
        void allocateGridConnections()
        // ----------------------------------------------------------------
        {
            const   Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
 
            for (int c = 0; c < cf.size(); ++c) {
                const int nf = cf[c].size();
                for (auto& f : cf[c]) {
                    S_.incrementrowsize(f, nf - 1);
                }
            }
        }
 
 
        // ----------------------------------------------------------------
        void allocateBCConnections(const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            // Include system connections for periodic boundary
            // conditions, if any.  We make an arbitrary choice in
            // that the face/degree-of-freedom with the minimum
            // numerical id of the two periodic partners represents
            // the coupling.  Suppose <i_p> is this minimum dof-id.
            // We then need to introduce a *symmetric* coupling to
            // <i_p> to each of the dof's of the cell *NOT* containing
            // <i_p>.  This choice is implemented in the following
            // loop by introducing couplings only when dof1 (self) is
            // less than dof2 (other).
            //
            // See also: setBCConnections() and addCellContrib().
            //
            typedef typename GridInterface::CellIterator CI;
            typedef typename CI           ::FaceIterator FI;
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf   = flowSolution_.cellFaces_;
 
            if (!bdry_id_map_.empty()) {
                // At least one periodic BC.  Allocate corresponding
                // connections.
                for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                    for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                        if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
                            // dof-id of self
                            const int dof1 = cf[cell[c->index()]][f->localIndex()];
 
                            // dof-id of other
                            BdryIdMapIterator j =
                                bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
                            assert (j != bdry_id_map_.end());
                            const int c2   = j->second.first;
                            const int dof2 = cf[c2][j->second.second];
 
                            if (dof1 < dof2) {
                                // Allocate storage for the actual
                                // couplings.
                                //
                                const int ndof = cf.rowSize(c2);
                                S_.incrementrowsize(dof1, ndof); // self->other
                                for (int dof = 0; dof < ndof; ++dof) {
                                    int ii = cf[c2][dof];
                                    int pp = ppartner_dof_[ii];
                                    if ((pp != -1) && (pp != dof1) && (pp < ii)) {
                                        S_.incrementrowsize(pp, 1);
                                    }
                                    S_.incrementrowsize(ii, 1);  // other->self
                                }
                            }
                        }
                    }
                }
            }
        }
 
 
 
        // ----------------------------------------------------------------
        void setConnections(const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            setGridConnections();
            setBCConnections(bc);
 
            S_.endindices();
 
            const int nc = pgrid_->numberOfCells();
            std::vector<Scalar>(nc).swap(flowSolution_.pressure_);
            std::vector<Scalar>(nc).swap(g_);
            std::vector<Scalar>(nc).swap(L_);
 
            matrix_structure_valid_ = true;
        }
 
 
        // ----------------------------------------------------------------
        void setGridConnections()
        // ----------------------------------------------------------------
        {
            const   Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
 
            // Compute actual connections (the non-zero structure).
            for (int c = 0; c < cf.size(); ++c) {
                auto fb = cf[c].begin(), fe = cf[c].end();
 
                for (auto i = fb; i != fe; ++i) {
                    for (auto j = fb; j != fe; ++j) {
                        S_.addindex(*i, *j);
                    }
                }
            }
        }
 
 
        // ----------------------------------------------------------------
        void setBCConnections(const BCInterface& bc)
        // ----------------------------------------------------------------
        {
            // Include system connections for periodic boundary
            // conditions, if any.  We make an arbitrary choice in
            // that the face/degree-of-freedom with the minimum
            // numerical id of the two periodic partners represents
            // the coupling.  Suppose <i_p> is this minimum dof-id.
            // We then need to introduce a *symmetric* coupling to
            // <i_p> to each of the dof's of the cell *NOT* containing
            // <i_p>.  This choice is implemented in the following
            // loop by introducing couplings only when dof1 (self) is
            // less than dof2 (other).
            //
            // See also: allocateBCConnections() and addCellContrib().
            //
            typedef typename GridInterface::CellIterator CI;
            typedef typename CI           ::FaceIterator FI;
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf   = flowSolution_.cellFaces_;
 
            if (!bdry_id_map_.empty()) {
                // At least one periodic BC.  Assign periodic
                // connections.
                for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                    for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                        if (f->boundary() && bc.flowCond(*f).isPeriodic()) {
                            // dof-id of self
                            const int dof1 = cf[cell[c->index()]][f->localIndex()];
 
                            // dof-id of other
                            BdryIdMapIterator j =
                                bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
                            assert (j != bdry_id_map_.end());
                            const int c2   = j->second.first;
                            const int dof2 = cf[c2][j->second.second];
 
                            if (dof1 < dof2) {
                                // Assign actual couplings.
                                //
                                const int ndof = cf.rowSize(c2);
                                for (int dof = 0; dof < ndof; ++dof) {
                                    int ii = cf[c2][dof];
                                    int pp = ppartner_dof_[ii];
                                    if ((pp != -1) && (pp != dof1) && (pp < ii)) {
                                        ii = pp;
                                    }
                                    S_.addindex(dof1, ii);  // self->other
                                    S_.addindex(ii, dof1);  // other->self
                                    S_.addindex(dof2, ii);
                                    S_.addindex(ii, dof2);
                                }
                            }
                        }
                    }
                }
            }
        }
 
 
 
        // ----------------------------------------------------------------
        template<class FluidInterface>
        void assembleDynamic(const FluidInterface&      fl ,
                             const std::vector<double>& sat,
                             const BCInterface&         bc ,
                             const std::vector<double>& src)
        // ----------------------------------------------------------------
        {
            typedef typename GridInterface::CellIterator CI;
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf   = flowSolution_.cellFaces_;
 
            std::vector<Scalar>   data_store(max_ncf_ * max_ncf_);
            std::vector<Scalar>   e         (max_ncf_);
            std::vector<Scalar>   rhs       (max_ncf_);
            std::vector<Scalar>   gflux     (max_ncf_);
 
            std::vector<FaceType> facetype  (max_ncf_);
            std::vector<Scalar>   condval   (max_ncf_);
            std::vector<int>      ppartner  (max_ncf_);
 
            // Clear residual data
            S_   = 0.0;
            rhs_ = 0.0;
 
            std::fill(g_.begin(), g_.end(), Scalar(0.0));
            std::fill(L_.begin(), L_.end(), Scalar(0.0));
            std::fill(e .begin(), e .end(), Scalar(1.0));
 
            // We will have to regularize resulting system if there
            // are no prescribed pressures (i.e., Dirichlet BC's).
            do_regularization_ = true;
 
            // Assemble dynamic contributions for each cell
            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                const int ci = c->index();
                const int c0 = cell[ci];            assert (c0 < cf.size());
                const int nf = cf[c0].size();
 
                rhs  .resize(nf);
                gflux.resize(nf);
 
                setExternalContrib(c, c0, bc, src[ci], rhs,
                                   facetype, condval, ppartner);
 
                ip_.computeDynamicParams(c, fl, sat);
 
                SharedFortranMatrix    S(nf, nf, &data_store[0]);
                ip_.getInverseMatrix(c, S);
 
                std::fill(gflux.begin(), gflux.end(), Scalar(0.0));
                ip_.gravityFlux(c, gflux);
 
                ImmutableFortranMatrix one(nf, 1, &e[0]);
                buildCellContrib(c0, one, gflux, S, rhs);
 
                addCellContrib(S, rhs, facetype, condval, ppartner, cf[c0]);
            }
        }
 
 
 
        // ----------------------------------------------------------------
        void solveLinearSystem(double residual_tolerance, int verbosity_level, int maxit)
        // ----------------------------------------------------------------
        {
            // Adapted from DuMux...
            Scalar residTol = residual_tolerance;
 
            typedef Dune::BCRSMatrix <MatrixBlockType>        MatrixT;
            typedef Dune::BlockVector<VectorBlockType>        VectorT;
            typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Adapter;
 
            // Regularize the matrix (only for pure Neumann problems...)
            if (do_regularization_) {
                S_[0][0] *= 2;
            }
            Adapter opS(S_);
 
            // Construct preconditioner.
            Dune::SeqILU<MatrixT,VectorT,VectorT> precond(S_, 1.0);
 
            // Construct solver for system of linear equations.
            Dune::CGSolver<VectorT> linsolve(opS, precond, residTol,
                                             (maxit>0)?maxit:S_.N(), verbosity_level);
 
            Dune::InverseOperatorResult result;
            soln_ = 0.0;
 
            // Solve system of linear equations to recover
            // face/contact pressure values (soln_).
            linsolve.apply(soln_, rhs_, result);
            if (!result.converged) {
                OPM_THROW(std::runtime_error,
                          "Linear solver failed to converge in " +
                          std::to_string(result.iterations) + " iterations.\n"
                          "Residual reduction achieved is " +
                          std::to_string(result.reduction) + '\n');
            }
        }
 
 
 
        // ------------------ AMG typedefs --------------------
 
        // Representation types for linear system.
        typedef Dune::BCRSMatrix <MatrixBlockType>        Matrix;
        typedef Dune::BlockVector<VectorBlockType>        Vector;
        typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
 
        // AMG specific types.
        // Old:   FIRST_DIAGONAL 1, SYMMETRIC 1, SMOOTHER_ILU 1, ANISOTROPIC_3D 0
        // SPE10: FIRST_DIAGONAL 0, SYMMETRIC 1, SMOOTHER_ILU 0, ANISOTROPIC_3D 1
#ifndef FIRST_DIAGONAL
#define FIRST_DIAGONAL 1
#endif
#ifndef SYMMETRIC
#define SYMMETRIC 1
#endif
#ifndef SMOOTHER_ILU
#define SMOOTHER_ILU 1
#endif
#ifndef SMOOTHER_BGS
#define SMOOTHER_BGS 0
#endif
#ifndef ANISOTROPIC_3D
#define ANISOTROPIC_3D 0
#endif
 
#if FIRST_DIAGONAL
        typedef Dune::Amg::FirstDiagonal CouplingMetric;
#else
        typedef Dune::Amg::RowSum        CouplingMetric;
#endif
 
#if SYMMETRIC
        typedef Dune::Amg::SymmetricCriterion<Matrix,CouplingMetric>   CriterionBase;
#else
        typedef Dune::Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
#endif
 
#if SMOOTHER_BGS
      typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
#else
#if SMOOTHER_ILU
        typedef Dune::SeqILU<Matrix,Vector,Vector>         Smoother;
#else
        typedef Dune::SeqSSOR<Matrix,Vector,Vector>        Smoother;
#endif
#endif
        typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
 
 
        // --------- storing the AMG operator and preconditioner --------
        std::unique_ptr<Operator> opS_;
        typedef Dune::Preconditioner<Vector,Vector>   PrecondBase;
        std::unique_ptr<PrecondBase> precond_;
 
 
        // ----------------------------------------------------------------
        void solveLinearSystemAMG(double residual_tolerance, int verbosity_level,
                                  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
        // ----------------------------------------------------------------
        {
            typedef Dune::Amg::AMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation>
                Precond;
 
            // Adapted from upscaling.cc by Arne Rekdal, 2009
            Scalar residTol = residual_tolerance;
 
            if (!same_matrix) {
                // Regularize the matrix (only for pure Neumann problems...)
                if (do_regularization_) {
                    S_[0][0] *= 2;
                }
                opS_.reset(new Operator(S_));
 
                // Construct preconditioner.
                double relax = 1;
                typename Precond::SmootherArgs smootherArgs;
                smootherArgs.relaxationFactor = relax;
#if SMOOTHER_BGS
                smootherArgs.overlap =  Precond::SmootherArgs::none;
                smootherArgs.onthefly = false;
#endif
                Criterion criterion;
                criterion.setDebugLevel(verbosity_level);
#if ANISOTROPIC_3D
                criterion.setDefaultValuesAnisotropic(3, 2);
#endif
                criterion.setProlongationDampingFactor(prolong_factor);
                criterion.setBeta(1e-10);
                criterion.setNoPreSmoothSteps(smooth_steps);
                criterion.setNoPostSmoothSteps(smooth_steps);
                criterion.setGamma(1); // V-cycle; this is the default
                precond_.reset(new Precond(*opS_, criterion, smootherArgs));
            }
            // Construct solver for system of linear equations.
            Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
 
            Dune::InverseOperatorResult result;
            soln_ = 0.0;
            // Adapt initial guess such Dirichlet boundary conditions are
            // represented, i.e. soln_i=A_{ii}^-1 rhs_i
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
            for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
                bool isDirichlet=true;
                for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
                    if(ci.index()!=ri.index() && *ci!=0.0)
                        isDirichlet=false;
                if(isDirichlet)
                    soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
            }
            // Solve system of linear equations to recover
            // face/contact pressure values (soln_).
            linsolve.apply(soln_, rhs_, result);
            if (!result.converged) {
                OPM_THROW(std::runtime_error,
                          "Linear solver failed to converge in " +
                          std::to_string(result.iterations) + " iterations.\n"
                          "Residual reduction achieved is " +
                          std::to_string(result.reduction) + '\n');
            }
 
        }
 
 
        // ----------------------------------------------------------------
        void solveLinearSystemFastAMG(double residual_tolerance, int verbosity_level,
                                  int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
        // ----------------------------------------------------------------
        {
            typedef Dune::Amg::FastAMG<Operator,Vector> Precond;
 
            // Adapted from upscaling.cc by Arne Rekdal, 2009
            Scalar residTol = residual_tolerance;
 
            if (!same_matrix) {
                // Regularize the matrix (only for pure Neumann problems...)
                if (do_regularization_) {
                    S_[0][0] *= 2;
                }
                opS_.reset(new Operator(S_));
 
                // Construct preconditioner.
                typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<Matrix,CouplingMetric> > CritBase;
 
                typedef Dune::Amg::CoarsenCriterion<CritBase> Crit;
                Crit criterion;
                criterion.setDebugLevel(verbosity_level);
#if ANISOTROPIC_3D
                criterion.setDefaultValuesAnisotropic(3, 2);
#endif
                criterion.setProlongationDampingFactor(prolong_factor);
                criterion.setBeta(1e-10);
                Dune::Amg::Parameters parms;
                parms.setDebugLevel(verbosity_level);
                parms.setNoPreSmoothSteps(smooth_steps);
                parms.setNoPostSmoothSteps(smooth_steps);
                precond_.reset(new Precond(*opS_, criterion, parms));
            }
            // Construct solver for system of linear equations.
            Dune::GeneralizedPCGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
 
            Dune::InverseOperatorResult result;
            soln_ = 0.0;
 
            // Adapt initial guess such Dirichlet boundary conditions are
            // represented, i.e. soln_i=A_{ii}^-1 rhs_i
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
            for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
                bool isDirichlet=true;
                for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
                    if(ci.index()!=ri.index() && *ci!=0.0)
                        isDirichlet=false;
                if(isDirichlet)
                    soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
            }
            // Solve system of linear equations to recover
            // face/contact pressure values (soln_).
            linsolve.apply(soln_, rhs_, result);
            if (!result.converged) {
                OPM_THROW(std::runtime_error,
                          "Linear solver failed to converge in " +
                          std::to_string(result.iterations) + " iterations.\n"
                          "Residual reduction achieved is " +
                          std::to_string(result.reduction) + '\n');
            }
 
        }
 
        // ----------------------------------------------------------------
        void solveLinearSystemKAMG(double residual_tolerance, int verbosity_level,
                                   int maxit, double prolong_factor, bool same_matrix, int smooth_steps)
        // ----------------------------------------------------------------
        {
 
            typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> Precond;
            // Adapted from upscaling.cc by Arne Rekdal, 2009
            Scalar residTol = residual_tolerance;
            if (!same_matrix) {
                // Regularize the matrix (only for pure Neumann problems...)
                if (do_regularization_) {
                    S_[0][0] *= 2;
                }
                opS_.reset(new Operator(S_));
 
                // Construct preconditioner.
                double relax = 1;
                typename Precond::SmootherArgs smootherArgs;
                smootherArgs.relaxationFactor = relax;
#if SMOOTHER_BGS
                smootherArgs.overlap =  Precond::SmootherArgs::none;
                smootherArgs.onthefly = false;
#endif
                Criterion criterion;
                criterion.setDebugLevel(verbosity_level);
#if ANISOTROPIC_3D
                criterion.setDefaultValuesAnisotropic(3, 2);
#endif
                criterion.setProlongationDampingFactor(prolong_factor);
                criterion.setBeta(1e-10);
                criterion.setNoPreSmoothSteps(smooth_steps);
                criterion.setNoPostSmoothSteps(smooth_steps);
                criterion.setGamma(2);
                precond_.reset(new Precond(*opS_, criterion, smootherArgs));
            }
            // Construct solver for system of linear equations.
            Dune::CGSolver<Vector> linsolve(*opS_, dynamic_cast<Precond&>(*precond_), residTol, (maxit>0)?maxit:S_.N(), verbosity_level);
 
            Dune::InverseOperatorResult result;
            soln_ = 0.0;
            // Adapt initial guess such Dirichlet boundary conditions are
            // represented, i.e. soln_i=A_{ii}^-1 rhs_i
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstRowIterator RowIter;
            typedef typename Dune::BCRSMatrix <MatrixBlockType>::ConstColIterator ColIter;
            for(RowIter ri=S_.begin(); ri!=S_.end(); ++ri){
                bool isDirichlet=true;
                for(ColIter ci=ri->begin(); ci!=ri->end(); ++ci)
                    if(ci.index()!=ri.index() && *ci!=0.0)
                        isDirichlet=false;
                if(isDirichlet)
                    soln_[ri.index()]=rhs_[ri.index()]/S_[ri.index()][ri.index()];
            }
            // Solve system of linear equations to recover
            // face/contact pressure values (soln_).
            linsolve.apply(soln_, rhs_, result);
            if (!result.converged) {
                OPM_THROW(std::runtime_error,
                          "Linear solver failed to converge in " +
                          std::to_string(result.iterations) + " iterations.\n"
                          "Residual reduction achieved is " +
                          std::to_string(result.reduction) + '\n');
            }
 
        }
 
 
 
        // ----------------------------------------------------------------
        template<class FluidInterface>
        void computePressureAndFluxes(const FluidInterface&      r  ,
                                      const std::vector<double>& sat)
        // ----------------------------------------------------------------
        {
            typedef typename GridInterface::CellIterator CI;
 
            const std::vector<int>& cell = flowSolution_.cellno_;
            const Opm::SparseTable<int>& cf   = flowSolution_.cellFaces_;
 
            std::vector<Scalar>& p = flowSolution_.pressure_;
        Opm::SparseTable<Scalar>& v = flowSolution_.outflux_;
 
            //std::vector<double> mob(FluidInterface::NumberOfPhases);
            std::vector<double> pi   (max_ncf_);
            std::vector<double> gflux(max_ncf_);
            std::vector<double> Binv_storage(max_ncf_ * max_ncf_);
 
            // Assemble dynamic contributions for each cell
            for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
                const int c0 = cell[c->index()];
                const int nf = cf.rowSize(c0);
 
                pi   .resize(nf);
                gflux.resize(nf);
 
                // Extract contact pressures for cell 'c'.
                for (int i = 0; i < nf; ++i) {
                    pi[i] = soln_[cf[c0][i]];
                }
 
                // Compute cell pressure in cell 'c'.
                p[c0] = (g_[c0] +
                         std::inner_product(F_[c0].begin(), F_[c0].end(),
                                            pi.begin(), 0.0)) / L_[c0];
 
                std::transform(pi.begin(), pi.end(),
                               pi.begin(),
                               [&p, c0](const double& input) { return p[c0] - input; });
 
                // Recover fluxes from local system
                //    Bv = Bv_g + Cp - D\pi
                //
                // 1) Solve system Bv = Cp - D\pi
                //
                ip_.computeDynamicParams(c, r, sat);
 
                SharedFortranMatrix Binv(nf, nf, &Binv_storage[0]);
                ip_.getInverseMatrix(c, Binv);
                vecMulAdd_N(Scalar(1.0), Binv, &pi[0], Scalar(0.0), &v[c0][0]);
 
                // 2) Add gravity flux contributions (v <- v + v_g)
                //
                ip_.gravityFlux(c, gflux);
                std::transform(gflux.begin(), gflux.end(), v[c0].begin(),
                               v[c0].begin(),
                               std::plus<Scalar>());
            }
        }
 
 
 
 
        // ----------------------------------------------------------------
        void setExternalContrib(const typename GridInterface::CellIterator c,
                                const int c0, const BCInterface& bc,
                                const double           src,
                                std::vector<Scalar>&   rhs,
                                std::vector<FaceType>& facetype,
                                std::vector<double>&   condval,
                                std::vector<int>&      ppartner)
        // ----------------------------------------------------------------
        {
            typedef typename GridInterface::CellIterator::FaceIterator FI;
 
            const Opm::SparseTable<int>& cf = flowSolution_.cellFaces_;
 
            std::fill(rhs     .begin(), rhs     .end(), Scalar(0.0));
            std::fill(facetype.begin(), facetype.end(), Internal   );
            std::fill(condval .begin(), condval .end(), Scalar(0.0));
            std::fill(ppartner.begin(), ppartner.end(), -1         );
 
            g_[c0] = src;
 
            int k = 0;
            for (FI f = c->facebegin(); f != c->faceend(); ++f, ++k) {
                if (f->boundary()) {
                    const FlowBC& bcond = bc.flowCond(*f);
                    if (bcond.isDirichlet()) {
                        facetype[k]        = Dirichlet;
                        condval[k]         = bcond.pressure();
                        do_regularization_ = false;
                    } else if (bcond.isPeriodic()) {
                        BdryIdMapIterator j =
                            bdry_id_map_.find(bc.getPeriodicPartner(f->boundaryId()));
                        assert (j != bdry_id_map_.end());
 
                        facetype[k] = Periodic;
                        condval[k]  = bcond.pressureDifference();
                        ppartner[k] = cf[j->second.first][j->second.second];
                    } else {
                        assert (bcond.isNeumann());
                        facetype[k] = Neumann;
                        rhs[k]      = bcond.outflux();
                    }
                }
            }
        }
 
 
 
 
        // ----------------------------------------------------------------
        void buildCellContrib(const int                     c    ,
                              const ImmutableFortranMatrix& one  ,
                              const std::vector<Scalar>&    gflux,
                              SharedFortranMatrix&          S    ,
                              std::vector<Scalar>&          rhs)
        // ----------------------------------------------------------------
        {
            // Ft <- B^{-t} * ones([size(S,2),1])
            SharedFortranMatrix Ft(S.numRows(), 1, &F_[c][0]);
            matMulAdd_TN(Scalar(1.0), S, one, Scalar(0.0), Ft);
 
            L_[c]  = std::accumulate(Ft.data(), Ft.data() + Ft.numRows(), 0.0);
            g_[c] -= std::accumulate(gflux.begin(), gflux.end(), Scalar(0.0));
 
            // rhs <- v_g - rhs (== v_g - h)
            std::transform(gflux.begin(), gflux.end(), rhs.begin(),
                           rhs.begin(),
                           std::minus<Scalar>());
 
            // rhs <- rhs + g_[c]/L_[c]*F
            std::transform(rhs.begin(), rhs.end(), Ft.data(),
                           rhs.begin(),
                           axpby<Scalar>(Scalar(1.0), Scalar(g_[c] / L_[c])));
 
            // S <- S - F'*F/L_c
            symmetricUpdate(-Scalar(1.0)/L_[c], Ft, Scalar(1.0), S);
        }
 
 
 
        // ----------------------------------------------------------------
        template<class L2G>
        void addCellContrib(const SharedFortranMatrix&   S       ,
                            const std::vector<Scalar>&   rhs     ,
                            const std::vector<FaceType>& facetype,
                            const std::vector<Scalar>&   condval ,
                            const std::vector<int>&      ppartner,
                            const L2G&                   l2g)
        // ----------------------------------------------------------------
        {
            int r = 0;
            for (auto i = l2g.begin(); i != l2g.end(); ++i, ++r) {
                // Indirection for periodic BC handling.
                int ii = *i;
 
                switch (facetype[r]) {
                case Dirichlet:
                    // Pressure is already known.  Assemble trivial
                    // equation of the form: a*x = a*p where 'p' is
                    // the known pressure value (i.e., condval[r]).
                    //
                    S_  [ii][ii] = S(r,r);
                    rhs_[ii]     = S(r,r) * condval[r];
                    continue;
                case Periodic:
                    // Periodic boundary condition.  Contact pressures
                    // linked by relations of the form
                    //
                    //      a*(x0 - x1) = a*pd
                    //
                    // where 'pd' is the known pressure difference
                    // x0-x1 (condval[r]).  Preserve matrix symmetry
                    // by assembling both of the equations
                    //
                    //      a*(x0-x1) = a*  pd, and  (1)
                    //      a*(x1-x0) = a*(-pd)      (2)
                    //
                    assert ((0 <= ppartner[r]) && (ppartner[r] < int(rhs_.size())));
                    assert (ii != ppartner[r]);
 
                    {
                        const double a = S(r,r), b = a * condval[r];
 
                        // Equation (1)
                        S_  [         ii][         ii] += a;
                        S_  [         ii][ppartner[r]] -= a;
                        rhs_[         ii]              += b;
 
                        // Equation (2)
                        S_  [ppartner[r]][         ii] -= a;
                        S_  [ppartner[r]][ppartner[r]] += a;
                        rhs_[ppartner[r]]              -= b;
                    }
 
                    ii = std::min(ii, ppartner[r]);
 
                    // INTENTIONAL FALL-THROUGH!
                    // IOW: Don't insert <break;> here!
                    //
                default:
                    int c = 0;
                    for (auto j = l2g.begin(); j != l2g.end(); ++j, ++c) {
                        // Indirection for periodic BC handling.
                        int jj = *j;
 
                        if (facetype[c] == Dirichlet) {
                            rhs_[ii] -= S(r,c) * condval[c];
                            continue;
                        }
                        if (facetype[c] == Periodic) {
                            assert ((0 <= ppartner[c]) && (ppartner[c] < int(rhs_.size())));
                            assert (jj != ppartner[c]);
                            if (ppartner[c] < jj) {
                                rhs_[ii] -= S(r,c) * condval[c];
                                jj = ppartner[c];
                            }
                        }
                        S_[ii][jj] += S(r,c);
                    }
                    break;
                }
                rhs_[ii] += rhs[r];
            }
        }
    };
} // namespace Opm
 
#endif // OPENRS_INCOMPFLOWSOLVERHYBRID_HEADER