Loading...
Searching...
No Matches
IncompFlowSolverHybrid.hpp
// -*- 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
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.
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< GridInterface, typename RelpermPolicy::template ResProp< GridInterface::Dimension >::Type, BoundaryConditions, RelpermPolicy::template InnerProduct >::SolutionType
const FlowSolution & SolutionType
Definition IncompFlowSolverHybrid.hpp:811
Opm::IncompFlowSolverHybrid::printSystem
void printSystem(const std::string &prefix)
Output current system of linear equations to permanent storage in files.
Definition IncompFlowSolverHybrid.hpp:888
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 ...
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.
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.
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::GridInterface
Opm
Inverting small matrices.
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.
Definition Matrix.hpp:829
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).
Definition Matrix.hpp:913
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.
Definition Matrix.hpp:1252