Merge remote-tracking branch 'bska/master'

This commit is contained in:
Atgeirr Flø Rasmussen
2013-05-16 09:55:11 +02:00
83 changed files with 7373 additions and 23 deletions

303
opm/autodiff/AutoDiff.hpp Normal file
View File

@@ -0,0 +1,303 @@
/*===========================================================================
//
// File: AutoDiff.hpp
//
// Created: 2013-04-29 10:51:23+0200
//
// Authors: Knut-Andreas Lie <Knut-Andreas.Lie@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Xavier Raynaud <Xavier.Raynaud@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cmath>
#ifndef OPM_AUTODIFF_HPP_HEADER
#define OPM_AUTODIFF_HPP_HEADER
namespace AutoDiff {
template <typename Scalar>
class Forward {
public:
static Forward
constant(const Scalar x)
{
// Constant is function with zero derivative.
return function(x, Scalar(0));
}
static Forward
variable(const Scalar x)
{
// Variable is function with unit derivative (wrt. itself).
return function(x, Scalar(1));
}
static Forward
function(const Scalar x, const Scalar dx)
{
return Forward(x, dx);
}
void
operator +=(const Scalar& rhs)
{
val_ += rhs;
}
void
operator +=(const Forward& rhs)
{
val_ += rhs.val_;
der_ += rhs.der_;
}
void
operator -=(const Scalar& rhs)
{
val_ -= rhs;
}
void
operator -=(const Forward& rhs)
{
val_ -= rhs.val_;
der_ -= rhs.der_;
}
void
operator *=(const Scalar& rhs)
{
val_ *= rhs;
der_ *= rhs;
}
void
operator *=(const Forward& rhs)
{
der_ = der_*rhs.val_ + val_*rhs.der_;
val_ *= rhs.val_;
}
void
operator /=(const Scalar& rhs)
{
val_ /= rhs;
der_ /= rhs;
}
void
operator /=(const Forward& rhs)
{
der_ = (der_*rhs.val_ - val_*rhs.der_) / (rhs.val_ * rhs.val_);
val_ /= rhs.val_;
}
template <class Ostream>
Ostream&
print(Ostream& os) const
{
os << "(x,dx) = (" << val_ << ',' << der_ << ")";
return os;
}
const Scalar val() const { return val_; }
const Scalar der() const { return der_; }
private:
Forward(const Scalar x, const Scalar dx)
: val_(x), der_(dx)
{}
Scalar val_;
Scalar der_;
};
template <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const Forward<Scalar>& fw)
{
return fw.print(os);
}
template <typename Scalar>
Forward<Scalar>
operator +(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
ret += rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator +(const T lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = rhs;
ret += Scalar(lhs);
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator +(const Forward<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
ret += Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator -(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
ret -= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator -(const T lhs,
const Forward<Scalar>& rhs)
{
return Forward<Scalar>::function(Scalar(lhs) - rhs.val(), -rhs.der());
}
template <typename Scalar, typename T>
Forward<Scalar>
operator -(const Forward<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
ret -= Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator *(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
ret *= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator *(const T lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = rhs;
ret *= Scalar(lhs);
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator *(const Forward<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
ret *= Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator /(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
ret /= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator /(const T lhs,
const Forward<Scalar>& rhs)
{
Scalar a = Scalar(lhs) / rhs.val();
Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val());
return Forward<Scalar>::function(a, b);
}
template <typename Scalar, typename T>
Forward<Scalar>
operator /(const Forward<Scalar>& lhs,
const T rhs)
{
Scalar a = lhs.val() / Scalar(rhs);
Scalar b = lhs.der() / Scalar(rhs);
return Forward<Scalar>::function(a, b);
}
template <typename Scalar>
Forward<Scalar>
cos(const Forward<Scalar>& x)
{
Scalar a = std::cos(x.val());
Scalar b = -std::sin(x.val()) * x.der();
return Forward<Scalar>::function(a, b);
}
template <typename Scalar>
Forward<Scalar>
sqrt(const Forward<Scalar>& x)
{
Scalar a = std::sqrt(x.val());
Scalar b = (Scalar(1.0) / (Scalar(2.0) * a)) * x.der();
return Forward<Scalar>::function(a, b);
}
} // namespace AutoDiff
namespace std {
using AutoDiff::cos;
using AutoDiff::sqrt;
}
#endif /* OPM_AUTODIFF_HPP_HEADER */

View File

@@ -0,0 +1,322 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#include <opm/autodiff/AutoDiff.hpp>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <vector>
#include <cassert>
namespace AutoDiff
{
template <typename Scalar>
class ForwardBlock
{
public:
/// Underlying types for scalar vectors and jacobians.
typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
typedef Eigen::SparseMatrix<Scalar> M;
/// Named constructor pattern used here.
static ForwardBlock null()
{
V val;
std::vector<M> jac;
return ForwardBlock(val, jac);
}
static ForwardBlock constant(const V& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
const int num_blocks = blocksizes.size();
// For constants, all jacobian blocks are zero.
jac.resize(num_blocks);
for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]);
}
return ForwardBlock(val, jac);
}
static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
const int num_blocks = blocksizes.size();
// First, set all jacobian blocks to zero...
jac.resize(num_blocks);
for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]);
}
// ... then set the one corrresponding to this variable to identity.
assert(blocksizes[index] == num_elem);
jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1));
for (typename M::Index row = 0; row < val.size(); ++row) {
jac[index].insert(row, row) = Scalar(1.0);
}
return ForwardBlock(val, jac);
}
static ForwardBlock function(const V& val, const std::vector<M>& jac)
{
return ForwardBlock(val, jac);
}
/// Construct a set of primary variables,
/// each initialized to a given vector.
static std::vector<ForwardBlock> variables(const std::vector<V>& initial_values)
{
const int num_vars = initial_values.size();
std::vector<int> bpat;
for (int v = 0; v < num_vars; ++v) {
bpat.push_back(initial_values[v].size());
}
std::vector<ForwardBlock> vars;
for (int v = 0; v < num_vars; ++v) {
vars.emplace_back(variable(v, initial_values[v], bpat));
}
return vars;
}
/// Operator +
ForwardBlock operator+(const ForwardBlock& rhs) const
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks();
for (int block = 0; block < num_blocks; ++block) {
assert(jac[block].rows() == rhs.jac_[block].rows());
assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] += rhs.jac_[block];
}
return function(val_ + rhs.val_, jac);
}
/// Operator -
ForwardBlock operator-(const ForwardBlock& rhs) const
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks();
for (int block = 0; block < num_blocks; ++block) {
assert(jac[block].rows() == rhs.jac_[block].rows());
assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] -= rhs.jac_[block];
}
return function(val_ - rhs.val_, jac);
}
/// Operator *
ForwardBlock operator*(const ForwardBlock& rhs) const
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols());
jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
}
return function(val_ * rhs.val_, jac);
}
/// Operator /
ForwardBlock operator/(const ForwardBlock& rhs) const
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols());
jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
}
return function(val_ / rhs.val_, jac);
}
/// I/O.
template <class Ostream>
Ostream&
print(Ostream& os) const
{
int num_blocks = jac_.size();
os << "Value =\n" << val_ << "\n\nJacobian =\n";
for (int i = 0; i < num_blocks; ++i) {
os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
}
return os;
}
/// Number of elements
int size() const
{
return val_.size();
}
/// Number of Jacobian blocks.
int numBlocks() const
{
return jac_.size();
}
/// Sizes (number of columns) of Jacobian blocks.
std::vector<int> blockPattern() const
{
const int nb = numBlocks();
std::vector<int> bp(nb);
for (int block = 0; block < nb; ++block) {
bp[block] = jac_[block].cols();
}
return bp;
}
/// Function value
const V& value() const
{
return val_;
}
/// Function derivatives
const std::vector<M>& derivative() const
{
return jac_;
}
private:
ForwardBlock(const V& val,
const std::vector<M>& jac)
: val_(val), jac_(jac)
{
#ifndef NDEBUG
const int num_elem = val_.size();
const int num_blocks = jac_.size();
for (int block = 0; block < num_blocks; ++block) {
assert(num_elem == jac_[block].rows());
}
#endif
}
V val_;
std::vector<M> jac_;
};
template <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const ForwardBlock<Scalar>& fw)
{
return fw.print(os);
}
/// Multiply with sparse matrix from the left.
template <typename Scalar>
ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::M& lhs,
const ForwardBlock<Scalar>& rhs)
{
int num_blocks = rhs.numBlocks();
std::vector<typename ForwardBlock<Scalar>::M> jac(num_blocks);
assert(lhs.cols() == rhs.value().rows());
for (int block = 0; block < num_blocks; ++block) {
jac[block] = lhs*rhs.derivative()[block];
}
typename ForwardBlock<Scalar>::V val = lhs*rhs.value().matrix();
return ForwardBlock<Scalar>::function(val, jac);
}
template <typename Scalar>
ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
}
template <typename Scalar>
ForwardBlock<Scalar> operator*(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
{
return rhs * lhs; // Commutative operation.
}
template <typename Scalar>
ForwardBlock<Scalar> operator+(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
}
template <typename Scalar>
ForwardBlock<Scalar> operator+(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
{
return rhs + lhs; // Commutative operation.
}
template <typename Scalar>
ForwardBlock<Scalar> operator-(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
}
template <typename Scalar>
ForwardBlock<Scalar> operator-(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
{
return lhs - ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
}
template <typename Scalar>
ForwardBlock<Scalar> operator/(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
}
template <typename Scalar>
ForwardBlock<Scalar> operator/(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
{
return lhs / ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
}
} // namespace Autodiff
#endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED

View File

@@ -0,0 +1,216 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED
#define OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED
// #include "AutoDiff.hpp"
// #include <Eigen/Eigen>
// #include <Eigen/Sparse>
#include <armadillo>
#include <vector>
#include <cassert>
namespace AutoDiff
{
template <typename Scalar>
class ForwardBlock
{
public:
/// Underlying types for scalar vectors and jacobians.
typedef arma::Col<Scalar> V;
typedef arma::SpMat<Scalar> M;
/// Named constructor pattern used here.
static ForwardBlock constant(const V& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
const int num_blocks = blocksizes.size();
// For constants, all jacobian blocks are zero.
jac.resize(num_blocks);
for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]);
}
return ForwardBlock(val, jac);
}
static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
const int num_blocks = blocksizes.size();
// First, set all jacobian blocks to zero...
jac.resize(num_blocks);
for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]);
}
// ... then set the one corrresponding to this variable to identity.
assert(blocksizes[index] == num_elem);
jac[index].eye(num_elem, num_elem);
return ForwardBlock(val, jac);
}
static ForwardBlock function(const V& val, const std::vector<M>& jac)
{
return ForwardBlock(val, jac);
}
/// Operator +
ForwardBlock operator+(const ForwardBlock& rhs)
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks();
for (int block = 0; block < num_blocks; ++block) {
assert(jac[block].n_rows == rhs.jac_[block].n_rows);
assert(jac[block].n_cols == rhs.jac_[block].n_cols);
jac[block] += rhs.jac_[block];
}
return function(val_ + rhs.val_, jac);
}
/// Operator -
ForwardBlock operator-(const ForwardBlock& rhs)
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks();
for (int block = 0; block < num_blocks; ++block) {
assert(jac[block].n_rows == rhs.jac_[block].n_rows);
assert(jac[block].n_cols == rhs.jac_[block].n_cols);
jac[block] -= rhs.jac_[block];
}
return function(val_ - rhs.val_, jac);
}
/// Operator *
ForwardBlock operator*(const ForwardBlock& rhs)
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].n_rows == rhs.jac_[block].n_rows);
assert(jac_[block].n_cols == rhs.jac_[block].n_cols);
jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
}
return function(val_ * rhs.val_, jac);
}
/// Operator /
ForwardBlock operator/(const ForwardBlock& rhs)
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].n_rows == rhs.jac_[block].n_rows);
assert(jac_[block].n_cols == rhs.jac_[block].n_cols);
jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
}
return function(val_ / rhs.val_, jac);
}
/// I/O.
template <class Ostream>
Ostream&
print(Ostream& os) const
{
int num_blocks = jac_.size();
os << "Value =\n" << val_ << "\n\nJacobian =\n";
for (int i = 0; i < num_blocks; ++i) {
os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
}
return os;
}
/// Number of variables or Jacobian blocks.
int numBlocks() const
{
return jac_.size();
}
/// Function value
const V& value() const
{
return val_;
}
/// Function derivatives
const std::vector<M>& derivative() const
{
return jac_;
}
private:
ForwardBlock(const V& val,
const std::vector<M>& jac)
: val_(val), jac_(jac)
{
#ifndef NDEBUG
const int num_elem = val_.size();
const int num_blocks = jac_.size();
for (int block = 0; block < num_blocks; ++block) {
assert(num_elem == jac_[block].n_rows);
}
#endif
}
V val_;
std::vector<M> jac_;
};
template <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const ForwardBlock<Scalar>& fw)
{
return fw.print(os);
}
/// Multiply with sparse matrix from the left.
template <typename Scalar>
ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::M& lhs,
const ForwardBlock<Scalar>& rhs)
{
int num_blocks = rhs.numBlocks();
std::vector<typename ForwardBlock<Scalar>::M> jac(num_blocks);
assert(lhs.n_cols == rhs.value().n_rows);
for (int block = 0; block < num_blocks; ++block) {
jac[block] = lhs*rhs.derivative()[block];
}
typename ForwardBlock<Scalar>::V val = lhs*rhs.value().matrix();
return ForwardBlock<Scalar>::function(val, jac);
}
} // namespace Autodiff
#endif // OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED

View File

@@ -0,0 +1,306 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
#define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/grid.h>
// -------------------- class HelperOps --------------------
/// Contains vectors and sparse matrices that represent subsets or
/// operations on (AD or regular) vectors of data.
struct HelperOps
{
typedef AutoDiff::ForwardBlock<double>::M M;
typedef AutoDiff::ForwardBlock<double>::V V;
/// A list of internal faces.
typedef Eigen::Array<int, Eigen::Dynamic, 1> IFaces;
IFaces internal_faces;
/// Extract for each face the difference of its adjacent cells'values.
M ngrad;
/// Extract for each face the average of its adjacent cells' values.
M caver;
/// Extract for each cell the sum of its adjacent faces' (signed) values.
M div;
/// Constructs all helper vectors and matrices.
HelperOps(const UnstructuredGrid& grid)
{
const int nc = grid.number_of_cells;
const int nf = grid.number_of_faces;
// Define some neighbourhood-derived helper arrays.
typedef Eigen::Array<int, Eigen::Dynamic, 1> OneColInt;
typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
TwoColInt nb = Eigen::Map<TwoColInt>(grid.face_cells, nf, 2);
// std::cout << "nb = \n" << nb << std::endl;
TwoColBool nbib = nb >= 0;
OneColBool ifaces = nbib.rowwise().all();
const int num_internal = ifaces.cast<int>().sum();
// std::cout << num_internal << " internal faces." << std::endl;
TwoColInt nbi(num_internal, 2);
internal_faces.resize(num_internal);
int fi = 0;
for (int f = 0; f < nf; ++f) {
if (ifaces[f]) {
internal_faces[fi] = f;
nbi.row(fi) = nb.row(f);
++fi;
}
}
// std::cout << "nbi = \n" << nbi << std::endl;
// Create matrices.
ngrad.resize(num_internal, nc);
caver.resize(num_internal, nc);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> ngrad_tri;
std::vector<Tri> caver_tri;
ngrad_tri.reserve(2*num_internal);
caver_tri.reserve(2*num_internal);
for (int i = 0; i < num_internal; ++i) {
ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
caver_tri.emplace_back(i, nbi(i,0), 0.5);
caver_tri.emplace_back(i, nbi(i,1), 0.5);
}
ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
div = ngrad.transpose();
}
};
// -------------------- debugger output helpers --------------------
#if !defined(NDEBUG)
#include <cstdio>
#endif // !defined(NDEBUG)
namespace {
#if !defined(NDEBUG)
void
printSparseMatrix(const Eigen::SparseMatrix<double>& A,
std::FILE* fp)
{
typedef Eigen::SparseMatrix<double>::Index Index;
const Index osize = A.outerSize();
for (Index k = 0; k < osize; ++k) {
for (Eigen::SparseMatrix<double>::InnerIterator
i(A, k); i ; ++i) {
std::fprintf(fp, "%lu %lu %26.18e\n",
static_cast<unsigned long>(i.row() + 1),
static_cast<unsigned long>(i.col() + 1),
i.value());
}
}
}
void
printSparseMatrix(const Eigen::SparseMatrix<double>& A ,
const char* const fn)
{
std::FILE* fp;
fp = std::fopen(fn, "w");
if (fp != 0) {
printSparseMatrix(A, fp);
}
std::fclose(fp);
}
#endif // !defined(NDEBUG)
// -------------------- upwinding helper class --------------------
/// Upwind selection in absence of counter-current flow (i.e.,
/// without effects of gravity and/or capillary pressure).
template <typename Scalar>
class UpwindSelector {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
UpwindSelector(const UnstructuredGrid& g,
const HelperOps& h,
const typename ADB::V& ifaceflux)
{
typedef HelperOps::IFaces::Index IFIndex;
const IFIndex nif = h.internal_faces.size();
assert(nif == ifaceflux.size());
// Define selector structure.
typedef typename Eigen::Triplet<Scalar> Triplet;
std::vector<Triplet> s; s.reserve(nif);
for (IFIndex iface = 0; iface < nif; ++iface) {
const int f = h.internal_faces[iface];
const int c1 = g.face_cells[2*f + 0];
const int c2 = g.face_cells[2*f + 1];
assert ((c1 >= 0) && (c2 >= 0));
// Select upwind cell.
const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
s.push_back(Triplet(iface, c, Scalar(1)));
}
// Assemble explicit selector operator.
select_.resize(nif, g.number_of_cells);
select_.setFromTriplets(s.begin(), s.end());
}
/// Apply selector to multiple per-cell quantities.
std::vector<ADB>
select(const std::vector<ADB>& xc) const
{
// Absence of counter-current flow means that the same
// selector applies to all quantities, 'x', defined per
// cell.
std::vector<ADB> xf; xf.reserve(xc.size());
for (typename std::vector<ADB>::const_iterator
b = xc.begin(), e = xc.end(); b != e; ++b)
{
xf.push_back(select_ * (*b));
}
return xf;
}
/// Apply selector to single per-cell ADB quantity.
ADB select(const ADB& xc) const
{
return select_*xc;
}
/// Apply selector to single per-cell constant quantity.
typename ADB::V select(const typename ADB::V& xc) const
{
return (select_*xc.matrix()).array();
}
private:
typename ADB::M select_;
};
}
namespace {
template <typename Scalar, class IntVec>
Eigen::SparseMatrix<Scalar>
constructSubsetSparseMatrix(const int full_size, const IntVec& indices)
{
typedef Eigen::Triplet<Scalar> Tri;
const int subset_size = indices.size();
std::vector<Tri> triplets(subset_size);
for (int i = 0; i < subset_size; ++i) {
triplets[i] = Tri(i, indices[i], 1);
}
Eigen::SparseMatrix<Scalar> sub(subset_size, full_size);
sub.setFromTriplets(triplets.begin(), triplets.end());
return sub;
}
template <typename Scalar, class IntVec>
Eigen::SparseMatrix<Scalar>
constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
{
return constructSubsetSparseMatrix<Scalar>(full_size, indices).transpose();
}
} // anon namespace
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiff::ForwardBlock<Scalar>
subset(const AutoDiff::ForwardBlock<Scalar>& x,
const IntVec& indices)
{
return ::constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
}
/// Returns x(indices).
template <typename Scalar, class IntVec>
Eigen::Array<Scalar, Eigen::Dynamic, 1>
subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices)
{
return (::constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
}
/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
template <typename Scalar, class IntVec>
AutoDiff::ForwardBlock<Scalar>
superset(const AutoDiff::ForwardBlock<Scalar>& x,
const IntVec& indices,
const int n)
{
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x;
}
/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
template <typename Scalar, class IntVec>
Eigen::Array<Scalar, Eigen::Dynamic, 1>
superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices,
const int n)
{
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x;
}
/// Construct square sparse matrix with the
/// elements of d on the diagonal.
/// Need to mark this as inline since it is defined in a header and not a template.
inline
AutoDiff::ForwardBlock<double>::M
spdiag(const AutoDiff::ForwardBlock<double>::V& d)
{
typedef AutoDiff::ForwardBlock<double>::M M;
const int n = d.size();
M mat(n, n);
mat.reserve(Eigen::ArrayXi::Ones(n, 1));
for (M::Index i = 0; i < n; ++i) {
mat.insert(i, i) = d[i];
}
return mat;
}
#endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

View File

@@ -0,0 +1,134 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_AUTODIFFVEC_HEADER_INCLUDED
#define OPM_AUTODIFFVEC_HEADER_INCLUDED
#include <opm/autodiff/AutoDiff.hpp>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
namespace AutoDiff
{
template <typename Scalar>
class ForwardVec
{
public:
/// Underlying types for scalar vectors and jacobians.
typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
typedef Eigen::SparseMatrix<Scalar> M;
/// Named constructor pattern used here.
static ForwardVec constant(const V& val)
{
return ForwardVec(val);
}
static ForwardVec variable(const V& val)
{
ForwardVec ret(val);
// typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
// D ones = V::Ones(val.size()).matrix().asDiagonal();
// ret.jac_ = ones;
ret.jac_.reserve(Eigen::VectorXi::Constant(val.size(), 1));
for (typename M::Index row = 0; row < val.size(); ++row) {
ret.jac_.insert(row, row) = Scalar(1.0);
}
ret.jac_.makeCompressed();
return ret;
}
static ForwardVec function(const V& val, const M& jac)
{
return ForwardVec(val, jac);
}
/// Operators.
ForwardVec operator+(const ForwardVec& rhs)
{
return function(val_ + rhs.val_, jac_ + rhs.jac_);
}
/// Operators.
ForwardVec operator-(const ForwardVec& rhs)
{
return function(val_ - rhs.val_, jac_ - rhs.jac_);
}
/// Operators.
ForwardVec operator*(const ForwardVec& rhs)
{
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
return function(val_ * rhs.val_, D2*jac_ + D1*rhs.jac_);
}
/// Operators.
ForwardVec operator/(const ForwardVec& rhs)
{
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal();
D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
return function(val_ / rhs.val_, D3 * (D2*jac_ - D1*rhs.jac_));
}
/// I/O.
template <class Ostream>
Ostream&
print(Ostream& os) const
{
os << "val =\n" << val_ << "\n\njac =\n" << jac_ << "\n";
return os;
}
private:
explicit ForwardVec(const V& val)
: val_(val), jac_(val.size(), val.size())
{
}
ForwardVec(const V& val, const M& jac)
: val_(val), jac_(jac)
{
}
V val_;
M jac_;
};
template <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const ForwardVec<Scalar>& fw)
{
return fw.print(os);
}
} // namespace Autodiff
#endif // OPM_AUTODIFFVEC_HEADER_INCLUDED

View File

@@ -0,0 +1,557 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
namespace Opm
{
// Making these typedef to make the code more readable.
typedef BlackoilPropsAd::ADB ADB;
typedef BlackoilPropsAd::V V;
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAd::BlackoilPropsAd(const BlackoilPropertiesInterface& props)
: props_(props),
pu_(props.phaseUsage())
{
}
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int BlackoilPropsAd::numDimensions() const
{
return props_.numDimensions();
}
/// \return N, the number of cells.
int BlackoilPropsAd::numCells() const
{
return props_.numCells();
}
/// \return Array of N porosity values.
const double* BlackoilPropsAd::porosity() const
{
return props_.porosity();
}
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* BlackoilPropsAd::permeability() const
{
return props_.permeability();
}
////////////////////////////
// Fluid interface //
////////////////////////////
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
const double* BlackoilPropsAd::surfaceDensity() const
{
return props_.surfaceDensity();
}
// ------ Viscosity ------
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAd::muWat(const V& pw,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
THROW("Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
ASSERT(pw.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block mu(n, np);
props_.viscosity(n, pw.data(), z.data(), cells.data(), mu.data(), 0);
return mu.col(pu_.phase_pos[Water]);
}
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAd::muOil(const V& po,
const V& rs,
const Cells& cells) const
{
if (!pu_.phase_used[Oil]) {
THROW("Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
ASSERT(po.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
if (pu_.phase_used[Gas]) {
// Faking a z with the right ratio:
// rs = zg/zo
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
z.col(pu_.phase_pos[Gas]) = rs;
}
Block mu(n, np);
props_.viscosity(n, po.data(), z.data(), cells.data(), mu.data(), 0);
return mu.col(pu_.phase_pos[Oil]);
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAd::muGas(const V& pg,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
THROW("Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
ASSERT(pg.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block mu(n, np);
props_.viscosity(n, pg.data(), z.data(), cells.data(), mu.data(), 0);
return mu.col(pu_.phase_pos[Gas]);
}
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAd::muWat(const ADB& pw,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
THROW("Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
ASSERT(pw.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block mu(n, np);
Block dmu(n, np);
props_.viscosity(n, pw.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Water]));
const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmu_diag * pw.derivative()[block];
}
return ADB::function(mu.col(pu_.phase_pos[Water]), jacs);
}
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAd::muOil(const ADB& po,
const ADB& rs,
const Cells& cells) const
{
if (!pu_.phase_used[Oil]) {
THROW("Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
ASSERT(po.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
if (pu_.phase_used[Gas]) {
// Faking a z with the right ratio:
// rs = zg/zo
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
z.col(pu_.phase_pos[Gas]) = rs.value();
}
Block mu(n, np);
Block dmu(n, np);
props_.viscosity(n, po.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Oil]));
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
// For now, we deliberately ignore the derivative with respect to rs,
// since the BlackoilPropertiesInterface class does not evaluate it.
// We would add to the next line: + dmu_drs_diag * rs.derivative()[block]
jacs[block] = dmu_diag * po.derivative()[block];
}
return ADB::function(mu.col(pu_.phase_pos[Oil]), jacs);
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAd::muGas(const ADB& pg,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
THROW("Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
ASSERT(pg.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block mu(n, np);
Block dmu(n, np);
props_.viscosity(n, pg.value().data(), z.data(), cells.data(), mu.data(), dmu.data());
ADB::M dmu_diag = spdiag(dmu.col(pu_.phase_pos[Gas]));
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmu_diag * pg.derivative()[block];
}
return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
}
// ------ Formation volume factor (b) ------
// These methods all call the matrix() method, after which the variable
// (also) called 'matrix' contains, in each row, the A = RB^{-1} matrix for
// a cell. For three-phase black oil:
// A = [ bw 0 0
// 0 bo 0
// 0 b0*rs bw ]
// Where b = B^{-1}.
// Therefore, we extract the correct diagonal element, and are done.
// When we need the derivatives (w.r.t. p, since we don't do w.r.t. rs),
// we also get the following derivative matrix:
// A = [ dbw 0 0
// 0 dbo 0
// 0 db0*rs dbw ]
// Again, we just extract a diagonal element.
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V BlackoilPropsAd::bWat(const V& pw,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
THROW("Cannot call bWat(): water phase not present.");
}
const int n = cells.size();
ASSERT(pw.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block matrix(n, np*np);
props_.matrix(n, pw.data(), z.data(), cells.data(), matrix.data(), 0);
const int wi = pu_.phase_pos[Water];
return matrix.col(wi*np + wi);
}
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V BlackoilPropsAd::bOil(const V& po,
const V& rs,
const Cells& cells) const
{
if (!pu_.phase_used[Oil]) {
THROW("Cannot call bOil(): oil phase not present.");
}
const int n = cells.size();
ASSERT(po.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
if (pu_.phase_used[Gas]) {
// Faking a z with the right ratio:
// rs = zg/zo
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
z.col(pu_.phase_pos[Gas]) = rs;
}
Block matrix(n, np*np);
props_.matrix(n, po.data(), z.data(), cells.data(), matrix.data(), 0);
const int oi = pu_.phase_pos[Oil];
return matrix.col(oi*np + oi);
}
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V BlackoilPropsAd::bGas(const V& pg,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
THROW("Cannot call bGas(): gas phase not present.");
}
const int n = cells.size();
ASSERT(pg.size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block matrix(n, np*np);
props_.matrix(n, pg.data(), z.data(), cells.data(), matrix.data(), 0);
const int gi = pu_.phase_pos[Gas];
return matrix.col(gi*np + gi);
}
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB BlackoilPropsAd::bWat(const ADB& pw,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
THROW("Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
ASSERT(pw.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, pw.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
const int phase_ind = pu_.phase_pos[Water];
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
ADB::M db_diag = spdiag(dmatrix.col(column));
const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = db_diag * pw.derivative()[block];
}
return ADB::function(matrix.col(column), jacs);
}
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB BlackoilPropsAd::bOil(const ADB& po,
const ADB& rs,
const Cells& cells) const
{
if (!pu_.phase_used[Oil]) {
THROW("Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
ASSERT(po.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
if (pu_.phase_used[Gas]) {
// Faking a z with the right ratio:
// rs = zg/zo
z.col(pu_.phase_pos[Oil]) = V::Ones(n, 1);
z.col(pu_.phase_pos[Gas]) = rs.value();
}
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, po.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
const int phase_ind = pu_.phase_pos[Oil];
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
ADB::M db_diag = spdiag(dmatrix.col(column));
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
// For now, we deliberately ignore the derivative with respect to rs,
// since the BlackoilPropertiesInterface class does not evaluate it.
// We would add to the next line: + db_drs_diag * rs.derivative()[block]
jacs[block] = db_diag * po.derivative()[block];
}
return ADB::function(matrix.col(column), jacs);
}
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB BlackoilPropsAd::bGas(const ADB& pg,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
THROW("Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
ASSERT(pg.value().size() == n);
const int np = props_.numPhases();
Block z = Block::Zero(n, np);
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
const int phase_ind = pu_.phase_pos[Gas];
const int column = phase_ind*np + phase_ind; // Index of our sought diagonal column.
ADB::M db_diag = spdiag(dmatrix.col(column));
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = db_diag * pg.derivative()[block];
}
return ADB::function(matrix.col(column), jacs);
}
#if 0
// ------ Rs bubble point curve ------
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n bubble point values for Rs.
V BlackoilPropsAd::rsMax(const V& po,
const Cells& cells) const
{
}
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n bubble point values for Rs.
ADB BlackoilPropsAd::rsMax(const ADB& po,
const Cells& cells) const
{
}
#endif
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] sg Array of n gas saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n relperm values,
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
std::vector<V> BlackoilPropsAd::relperm(const V& sw,
const V& so,
const V& sg,
const Cells& cells) const
{
const int n = cells.size();
const int np = props_.numPhases();
Block s_all(n, np);
if (pu_.phase_used[Water]) {
ASSERT(sw.size() == n);
s_all.col(pu_.phase_pos[Water]) = sw;
}
if (pu_.phase_used[Oil]) {
ASSERT(so.size() == n);
s_all.col(pu_.phase_pos[Oil]) = so;
}
if (pu_.phase_used[Gas]) {
ASSERT(sg.size() == n);
s_all.col(pu_.phase_pos[Gas]) = sg;
}
Block kr(n, np);
props_.relperm(n, s_all.data(), cells.data(), kr.data(), 0);
std::vector<V> relperms;
relperms.reserve(3);
for (int phase = 0; phase < 3; ++phase) {
if (pu_.phase_used[phase]) {
relperms.emplace_back(kr.col(pu_.phase_pos[phase]));
} else {
relperms.emplace_back();
}
}
return relperms;
}
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] sg Array of n gas saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n relperm values,
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
std::vector<ADB> BlackoilPropsAd::relperm(const ADB& sw,
const ADB& so,
const ADB& sg,
const Cells& cells) const
{
const int n = cells.size();
const int np = props_.numPhases();
Block s_all(n, np);
if (pu_.phase_used[Water]) {
ASSERT(sw.value().size() == n);
s_all.col(pu_.phase_pos[Water]) = sw.value();
}
if (pu_.phase_used[Oil]) {
ASSERT(so.value().size() == n);
s_all.col(pu_.phase_pos[Oil]) = so.value();
} else {
THROW("BlackoilPropsAd::relperm() assumes oil phase is active.");
}
if (pu_.phase_used[Gas]) {
ASSERT(sg.value().size() == n);
s_all.col(pu_.phase_pos[Gas]) = sg.value();
}
Block kr(n, np);
Block dkr(n, np*np);
props_.relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
const int num_blocks = so.numBlocks();
std::vector<ADB> relperms;
relperms.reserve(3);
typedef const ADB* ADBPtr;
ADBPtr s[3] = { &sw, &so, &sg };
for (int phase1 = 0; phase1 < 3; ++phase1) {
if (pu_.phase_used[phase1]) {
const int phase1_pos = pu_.phase_pos[phase1];
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
}
for (int phase2 = 0; phase2 < 3; ++phase2) {
if (!pu_.phase_used[phase2]) {
continue;
}
const int phase2_pos = pu_.phase_pos[phase2];
// Assemble dkr1/ds2.
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
for (int block = 0; block < num_blocks; ++block) {
jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
}
}
relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
} else {
relperms.emplace_back(ADB::null());
}
}
return relperms;
}
} // namespace Opm

View File

@@ -0,0 +1,235 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILPROPSAD_HEADER_INCLUDED
#define OPM_BLACKOILPROPSAD_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm
{
class BlackoilPropertiesInterface;
/// This class is intended to present a fluid interface for
/// three-phase black-oil that is easy to use with the AD-using
/// simulators.
///
/// Most methods are available in two overloaded versions, one
/// taking a constant vector and returning the same, and one
/// taking an AD type and returning the same. Derivatives are not
/// returned separately by any method, only implicitly with the AD
/// version of the methods.
class BlackoilPropsAd
{
public:
/// Constructor wrapping an opm-core black oil interface.
explicit BlackoilPropsAd(const BlackoilPropertiesInterface& props);
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int numDimensions() const;
/// \return N, the number of cells.
int numCells() const;
/// \return Array of N porosity values.
const double* porosity() const;
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* permeability() const;
////////////////////////////
// Fluid interface //
////////////////////////////
typedef AutoDiff::ForwardBlock<double> ADB;
typedef ADB::V V;
typedef std::vector<int> Cells;
// ------ Canonical named indices for each phase ------
/// Canonical named indices for each phase.
enum PhaseIndex { Water = 0, Oil = 1, Gas = 2 };
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
const double* surfaceDensity() const;
// ------ Viscosity ------
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V muWat(const V& pw,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V muOil(const V& po,
const V& rs,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V muGas(const V& pg,
const Cells& cells) const;
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB muWat(const ADB& pw,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB muOil(const ADB& po,
const ADB& rs,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB muGas(const ADB& pg,
const Cells& cells) const;
// ------ Formation volume factor (b) ------
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V bWat(const V& pw,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V bOil(const V& po,
const V& rs,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V bGas(const V& pg,
const Cells& cells) const;
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB bWat(const ADB& pw,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB bOil(const ADB& po,
const ADB& rs,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB bGas(const ADB& pg,
const Cells& cells) const;
// ------ Rs bubble point curve ------
#if 0
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n bubble point values for Rs.
V rsMax(const V& po,
const Cells& cells) const;
/// Bubble point curve for Rs as function of oil pressure.
/// \param[in] po Array of n oil pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n bubble point values for Rs.
ADB rsMax(const ADB& po,
const Cells& cells) const;
#endif
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] sg Array of n gas saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n relperm values,
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
std::vector<V> relperm(const V& sw,
const V& so,
const V& sg,
const Cells& cells) const;
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] sg Array of n gas saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n relperm values,
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
std::vector<ADB> relperm(const ADB& sw,
const ADB& so,
const ADB& sg,
const Cells& cells) const;
private:
const BlackoilPropertiesInterface& props_;
PhaseUsage pu_;
};
} // namespace Opm
#endif // OPM_BLACKOILPROPSAD_HEADER_INCLUDED

View File

@@ -0,0 +1,460 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_IMPESTPFAAD_HEADER_INCLUDED
#define OPM_IMPESTPFAAD_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/wells.h>
#include <algorithm>
#include <cassert>
#include <vector>
#include <boost/shared_ptr.hpp>
struct UnstructuredGrid;
namespace {
std::vector<int>
buildAllCells(const int nc)
{
std::vector<int> all_cells(nc);
for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
return all_cells;
}
template <class GeoProps>
AutoDiff::ForwardBlock<double>::M
gravityOperator(const UnstructuredGrid& grid,
const HelperOps& ops ,
const GeoProps& geo )
{
const int nc = grid.number_of_cells;
std::vector<int> f2hf(2 * grid.number_of_faces, -1);
for (int c = 0, i = 0; c < nc; ++c) {
for (; i < grid.cell_facepos[c + 1]; ++i) {
const int f = grid.cell_faces[ i ];
const int p = 0 + (grid.face_cells[2*f + 0] != c);
f2hf[2*f + p] = i;
}
}
typedef AutoDiff::ForwardBlock<double>::V V;
typedef AutoDiff::ForwardBlock<double>::M M;
const V& gpot = geo.gravityPotential();
const V& trans = geo.transmissibility();
const HelperOps::IFaces::Index ni = ops.internal_faces.size();
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> grav; grav.reserve(2 * ni);
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
const int f = ops.internal_faces[ i ];
const int c1 = grid.face_cells[2*f + 0];
const int c2 = grid.face_cells[2*f + 1];
assert ((c1 >= 0) && (c2 >= 0));
const double dG1 = gpot[ f2hf[2*f + 0] ];
const double dG2 = gpot[ f2hf[2*f + 1] ];
const double t = trans[ f ];
grav.push_back(Tri(i, c1, t * dG1));
grav.push_back(Tri(i, c2, - t * dG2));
}
M G(ni, nc); G.setFromTriplets(grav.begin(), grav.end());
return G;
}
}
namespace Opm {
template <typename Scalar, class BOFluid>
class PressureDependentFluidData {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
typedef typename AutoDiff::ForwardBlock<Scalar>::V V;
typedef typename AutoDiff::ForwardBlock<Scalar>::M M;
PressureDependentFluidData(const int nc,
const BOFluid& fluid)
: nc_ (nc)
, np_ (fluid.numPhases())
, cells_(buildAllCells(nc))
, fluid_(fluid)
, A_ (nc_, np_ * np_)
, dA_ (nc_, np_ * np_)
, mu_ (nc_, np_ )
, dmu_ (nc_, np_ )
, kr_ (nc_, np_ )
, zero_ (ADB::V::Zero(nc, 1))
, one_ (ADB::V::Ones(nc, 1))
{
}
void
computeSatQuant(const BlackoilState& state)
{
const std::vector<double>& s = state.saturation();
assert (s.size() == std::vector<double>::size_type(nc_ * np_));
double* dkrds = 0; // Ignore rel-perm derivatives
fluid_.relperm(nc_, & s[0], & cells_[0],
kr_.data(), dkrds);
}
void
computePressQuant(const BlackoilState& state)
{
const std::vector<double>& p = state.pressure();
const std::vector<double>& z = state.surfacevol();
assert (p.size() == std::vector<double>::size_type(nc_ * 1 ));
assert (z.size() == std::vector<double>::size_type(nc_ * np_));
fluid_.matrix (nc_, & p[0], & z[0], & cells_[0],
A_ .data(), dA_ .data());
fluid_.viscosity(nc_, & p[0], & z[0], & cells_[0],
mu_.data(), /*dmu_.data()*/ 0);
dmu_.fill(0.0);
}
ADB
fvf(const int phase, const ADB& p) const
{
assert (0 <= phase);
assert (phase < np_ );
typedef typename ADB::V V;
const V A = A_ .block(0, phase * (np_ + 1), nc_, 1);
const V dA = dA_.block(0, phase * (np_ + 1), nc_, 1);
std::vector<typename ADB::M> jac(p.numBlocks());
assert(p.numBlocks() == 2);
jac[0] = spdiag(dA);
assert(jac[0].cols() == p.blockPattern()[0]);
jac[1] = M(A.rows(), p.blockPattern()[1]);
return one_ / ADB::function(A, jac);
}
typename ADB::V
phaseRelPerm(const int phase) const
{
const typename ADB::V kr = kr_.block(0, phase, nc_, 1);
return kr;
}
ADB
phaseViscosity(const int phase, const ADB& p) const
{
assert (0 <= phase);
assert (phase < np_ );
typedef typename ADB::V V;
const V mu = mu_ .block(0, phase, nc_, 1);
const V dmu = dmu_.block(0, phase, nc_, 1);
std::vector<typename ADB::M> jac(p.numBlocks());
assert(p.numBlocks() == 2);
jac[0] = spdiag(dmu);
assert(jac[0].cols() == p.blockPattern()[0]);
jac[1] = M(mu.rows(), p.blockPattern()[1]);
return ADB::function(mu, jac);
}
ADB
phaseDensity(const int phase, const ADB& p) const
{
typedef typename ADB::V V;
const double* rho0 = fluid_.surfaceDensity();
V rho = V::Zero(nc_, 1);
V drho = V::Zero(nc_, 1);
for (int i = 0; i < np_; ++i) {
rho += rho0[i] * A_ .block(0, phase*np_ + i, nc_, 1);
drho += rho0[i] * dA_.block(0, phase*np_ + i, nc_, 1);
}
assert (p.numBlocks() == 2);
std::vector<typename ADB::M> jac(p.numBlocks());
jac[0] = spdiag(drho);
jac[1] = M(rho.rows(), p.blockPattern()[1]);
assert (jac[0].cols() == p.blockPattern()[0]);
return ADB::function(rho, jac);
}
private:
const int nc_;
const int np_;
const std::vector<int> cells_;
const BOFluid& fluid_;
typedef Eigen::Array<Scalar,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DerivedQuant;
// Pressure dependent quantities (essentially B and \mu)
DerivedQuant A_ ;
DerivedQuant dA_ ;
DerivedQuant mu_ ;
DerivedQuant dmu_;
// Saturation dependent quantities (rel-perm only)
DerivedQuant kr_;
const typename ADB::V zero_;
const typename ADB::V one_ ;
};
template <class BOFluid, class GeoProps>
class ImpesTPFAAD {
public:
ImpesTPFAAD(const UnstructuredGrid& grid ,
const BOFluid& fluid,
const GeoProps& geo ,
const Wells& wells,
const LinearSolverInterface& linsolver)
: grid_ (grid)
, geo_ (geo)
, wells_ (wells)
, linsolver_(linsolver)
, pdepfdata_(grid.number_of_cells, fluid)
, ops_ (grid)
, grav_ (gravityOperator(grid_, ops_, geo_))
, cell_residual_ (ADB::null())
, well_residual_ (ADB::null())
{
}
void
solve(const double dt,
BlackoilState& state,
WellState& well_state)
{
pdepfdata_.computeSatQuant(state);
const double atol = 1.0e-15;
const double rtol = 5.0e-10;
const int maxit = 15;
assemble(dt, state, well_state);
const double r0 = residualNorm();
int it = 0;
bool resTooLarge = r0 > atol;
while (resTooLarge && (it < maxit)) {
solveJacobianSystem(state);
assemble(dt, state, well_state);
const double r = residualNorm();
resTooLarge = (r > atol) && (r > rtol*r0);
it += 1;
}
if (resTooLarge) {
THROW("Failed to compute converged pressure solution");
}
else {
computeFluxes();
}
}
private:
// Disallow copying and assignment
ImpesTPFAAD(const ImpesTPFAAD& rhs);
ImpesTPFAAD& operator=(const ImpesTPFAAD& rhs);
typedef PressureDependentFluidData<double, BOFluid> PDepFData;
typedef typename PDepFData::ADB ADB;
typedef typename ADB::V V;
typedef typename ADB::M M;
const UnstructuredGrid& grid_;
const GeoProps& geo_ ;
const Wells& wells_;
const LinearSolverInterface& linsolver_;
PDepFData pdepfdata_;
HelperOps ops_;
const M grav_;
ADB cell_residual_;
ADB well_residual_;
void
assemble(const double dt,
const BlackoilState& state,
const WellState& well_state)
{
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
const V& pv = geo_.poreVolume();
const int nc = grid_.number_of_cells;
const int np = state.numPhases();
const int nw = wells_.number_of_wells;
pdepfdata_.computePressQuant(state);
const Eigen::Map<const DataBlock> z0all(&state.surfacevol()[0], nc, np);
const DataBlock qall = DataBlock::Zero(nc, np);
const V delta_t = dt * V::Ones(nc, 1);
const V transi = subset(geo_.transmissibility(),
ops_.internal_faces);
const int num_perf = wells_.well_connpos[nw];
const std::vector<int> well_cells(wells_.well_cells,
wells_.well_cells + num_perf);
const V transw = Eigen::Map<const V>(wells_.WI, num_perf, 1);
// Initialize AD variables: p (cell pressures) and bhp (well bhp).
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V bhp0 = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
std::vector<V> vars0 = { p0, bhp0 };
std::vector<ADB> vars= ADB::variables(vars0);
const ADB& p = vars[0];
const ADB& bhp = vars[1];
std::vector<int> bpat = p.blockPattern();
// Compute T_ij * (p_i - p_j) and use for upwinding.
const ADB nkgradp = transi * (ops_.ngrad * p);
const UpwindSelector<double> upwind(grid_, ops_, nkgradp.value());
// Extract variables for perforation cell pressures
// and corresponding perforation well pressures.
const ADB p_perfcell = subset(p, well_cells);
// Construct matrix to map wells->perforations.
M well_to_perf(well_cells.size(), nw);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> w2p;
for (int w = 0; w < nw; ++w) {
for (int perf = wells_.well_connpos[w]; perf < wells_.well_connpos[w+1]; ++perf) {
w2p.emplace_back(perf, w, 1.0);
}
}
well_to_perf.setFromTriplets(w2p.begin(), w2p.end());
// Construct pressure difference vector for wells.
const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet!
// Finally construct well perforation pressures.
const ADB p_perfwell = well_to_perf*bhp + well_perf_dp;
const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
cell_residual_ = ADB::constant(pv, bpat);
#define COMPENSATE_FLOAT_PRECISION 0
#if COMPENSATE_FLOAT_PRECISION
ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat);
#endif
for (int phase = 0; phase < np; ++phase) {
const ADB cell_B = pdepfdata_.fvf(phase, p);
const ADB cell_rho = pdepfdata_.phaseDensity(phase, p);
const V kr = pdepfdata_.phaseRelPerm(phase);
const ADB mu = pdepfdata_.phaseViscosity(phase, p);
const ADB mf = upwind.select(kr / mu);
const ADB flux = mf * (nkgradp + (grav_ * cell_rho));
const ADB face_B = upwind.select(cell_B);
const V z0 = z0all.block(0, phase, nc, 1);
const V q = qall .block(0, phase, nc, 1);
#if COMPENSATE_FLOAT_PRECISION
const ADB divcontrib = delta_t * (ops_.div * (flux / face_B));
const V qcontrib = delta_t * q;
const ADB pvcontrib = ADB::constant(pv*z0, bpat);
const ADB component_contrib = pvcontrib + qcontrib;
divcontrib_sum = divcontrib_sum - cell_B*divcontrib;
cell_residual_ = cell_residual_ - (cell_B * component_contrib);
#else
const ADB component_contrib = pv*z0 + delta_t*(q - (ops_.div * (flux / face_B)));
cell_residual_ = cell_residual_ - (cell_B * component_contrib);
#endif
}
#if COMPENSATE_FLOAT_PRECISION
cell_residual_ = cell_residual_ + divcontrib_sum;
#endif
}
void
solveJacobianSystem(BlackoilState& state) const
{
const int nc = grid_.number_of_cells;
Eigen::SparseMatrix<double, Eigen::RowMajor> matr = cell_residual_.derivative()[0];
#if HACK_INCOMPRESSIBLE_GRAVITY
matr.coeffRef(0, 0) *= 2;
#endif
V dp(nc);
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(nc, matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
cell_residual_.value().data(), dp.data());
if (!rep.converged) {
THROW("ImpesTPFAAD::solve(): Linear solver convergence failure.");
}
const V p = p0 - dp;
std::copy(&p[0], &p[0] + nc, state.pressure().begin());
}
double
residualNorm() const
{
return cell_residual_.value().matrix().norm();
}
void
computeFluxes() const
{
}
};
} // namespace Opm
#endif /* OPM_IMPESTPFAAD_HEADER_INCLUDED */

View File

@@ -0,0 +1,638 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/autodiff/SimulatorIncompTwophaseAdfi.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
#include <opm/autodiff/TransportSolverTwophaseAd.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <numeric>
#include <fstream>
namespace Opm
{
class SimulatorIncompTwophaseAdfi::Impl
{
public:
Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
TwophaseState& state,
WellState& well_state);
private:
// Data.
// Parameters for output.
bool output_;
bool output_vtk_;
std::string output_dir_;
int output_interval_;
// Parameters for well control
bool check_well_controls_;
int max_well_control_iterations_;
// Parameters for transport solver.
int num_transport_substeps_;
std::string transport_solver_type_;
bool use_segregation_split_;
// Observed objects.
const UnstructuredGrid& grid_;
const IncompPropertiesInterface& props_;
const RockCompressibility* rock_comp_props_;
WellsManager& wells_manager_;
const Wells* wells_;
const std::vector<double>& src_;
const FlowBoundaryConditions* bcs_;
// Solvers
IncompTpfa psolver_;
boost::scoped_ptr<TransportSolverTwophaseInterface> tsolver_;
// Misc. data
std::vector<int> allcells_;
};
SimulatorIncompTwophaseAdfi::SimulatorIncompTwophaseAdfi(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity));
}
SimulatorReport SimulatorIncompTwophaseAdfi::run(SimulatorTimer& timer,
TwophaseState& state,
WellState& well_state)
{
return pimpl_->run(timer, state, well_state);
}
static void reportVolumes(std::ostream &os, double satvol[2], double tot_porevol_init,
double tot_injected[2], double tot_produced[2],
double injected[2], double produced[2],
double init_satvol[2])
{
std::cout.precision(5);
const int width = 18;
os << "\nVolume balance report (all numbers relative to total pore volume).\n";
os << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol_init
<< std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
os << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol_init
<< std::setw(width) << injected[1]/tot_porevol_init << std::endl;
os << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol_init
<< std::setw(width) << produced[1]/tot_porevol_init << std::endl;
os << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol_init
<< std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
os << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol_init
<< std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
os << " In-place + prod - inj: "
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
os << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
<< std::endl;
os.precision(8);
}
static void outputStateVtk(const UnstructuredGrid& grid,
const Opm::TwophaseState& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/vtk_files";
boost::filesystem::path fpath(vtkfilename.str());
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
}
static void outputVectorMatlab(const std::string& name,
const std::vector<int>& vec,
const int step,
const std::string& output_dir)
{
std::ostringstream fname;
fname << output_dir << "/" << name;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(file, "\n"));
}
static void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::TwophaseState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
watercut.write(os);
}
static void outputWellReport(const Opm::WellReport& wellreport,
const std::string& output_dir)
{
// Write well report.
std::string fname = output_dir + "/wellreport.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
wellreport.write(os);
}
static bool allNeumannBCs(const FlowBoundaryConditions* bcs)
{
if (bcs == NULL) {
return true;
} else {
return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE)
== bcs->type + bcs->nbc;
}
}
static bool allRateWells(const Wells* wells)
{
if (wells == NULL) {
return true;
}
const int nw = wells->number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells->ctrls[w];
if (wc->current >= 0) {
if (wc->type[wc->current] == BHP) {
return false;
}
}
}
return true;
}
SimulatorIncompTwophaseAdfi::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
: transport_solver_type_(param.getDefault<std::string>("transport_solver_type", "ad")),
use_segregation_split_(param.getDefault("use_segregation_split", false)),
grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
src_(src),
bcs_(bcs),
psolver_(grid, props, rock_comp_props, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs)
{
// Initialize transport solver.
if (transport_solver_type_ == "reorder") {
tsolver_.reset(new Opm::TransportSolverTwophaseReorder(grid,
props,
use_segregation_split_ ? gravity : NULL,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)));
} else if (transport_solver_type_ == "implicit") {
if (rock_comp_props && rock_comp_props->isActive()) {
THROW("The implicit transport solver cannot handle rock compressibility.");
}
if (use_segregation_split_) {
THROW("The implicit transport solver is not set up to use segregation splitting.");
}
std::vector<double> porevol;
computePorevolume(grid, props.porosity(), porevol);
tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(grid,
props,
porevol,
gravity,
psolver_.getHalfTrans(),
param));
} else if (transport_solver_type_ == "ad") {
if (rock_comp_props && rock_comp_props->isActive()) {
THROW("The implicit ad transport solver cannot handle rock compressibility.");
}
if (use_segregation_split_) {
THROW("The implicit ad transport solver is not set up to use segregation splitting.");
}
std::vector<double> porevol;
computePorevolume(grid, props.porosity(), porevol);
tsolver_.reset(new Opm::TransportSolverTwophaseAd(grid,
props,
linsolver,
gravity,
param));
} else {
THROW("Unknown transport solver type: " << transport_solver_type_);
}
// For output.
output_ = param.getDefault("output", true);
if (output_) {
output_vtk_ = param.getDefault("output_vtk", true);
output_dir_ = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir_);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
output_interval_ = param.getDefault("output_interval", 1);
}
// Well control related init.
check_well_controls_ = param.getDefault("check_well_controls", false);
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
// Misc init.
const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells_[cell] = cell;
}
}
SimulatorReport SimulatorIncompTwophaseAdfi::Impl::run(SimulatorTimer& timer,
TwophaseState& state,
WellState& well_state)
{
std::vector<double> transport_src;
// Initialisation.
std::vector<double> porevol;
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
} else {
computePorevolume(grid_, props_.porosity(), porevol);
}
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
std::vector<double> initial_porevol = porevol;
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
Opm::time::StopWatch transport_timer;
double ttime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init
<< " " << init_satvol[1]/tot_porevol_init << std::endl;
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
if (wells_) {
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
}
std::fstream tstep_os;
if (output_) {
std::string filename = output_dir_ + "/step_timing.param";
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
}
for (; !timer.done(); ++timer) {
// Report timestep and (optionally) write state to disk.
step_timer.start();
timer.report(std::cout);
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
if (transport_solver_type_ == "reorder") {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
}
SimulatorReport sreport;
// Solve pressure equation.
if (check_well_controls_) {
computeFractionalFlow(props_, allcells_, state.saturation(), fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
bool well_control_passed = !check_well_controls_;
int well_control_iteration = 0;
do {
// Run solver.
pressure_timer.start();
std::vector<double> initial_pressure = state.pressure();
psolver_.solve(timer.currentStepLength(), state, well_state);
// Renormalize pressure if rock is incompressible, and
// there are no pressure conditions (bcs or wells).
// It is deemed sufficient for now to renormalize
// using geometric volume instead of pore volume.
if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
&& allNeumannBCs(bcs_) && allRateWells(wells_)) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.0;
double av_press = 0.0;
double tot_vol = 0.0;
const int num_cells = grid_.number_of_cells;
for (int cell = 0; cell < num_cells; ++cell) {
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
tot_vol += grid_.cell_volumes[cell];
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
}
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
for (int well = 0; well < num_wells; ++well) {
well_state.bhp()[well] += ren_const;
}
}
// Stop timer and report.
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
sreport.pressure_time = pt;
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
if (rock_comp_props_ && rock_comp_props_->isActive()) {
initial_porevol = porevol;
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
}
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
wells_, well_state.perfRates(), transport_src);
// Solve transport.
transport_timer.start();
double stepsize = timer.currentStepLength();
if (num_transport_substeps_ != 1) {
stepsize /= double(num_transport_substeps_);
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_->solve(&initial_porevol[0], &transport_src[0], stepsize, state);
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize,
substep_injected, substep_produced);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (transport_solver_type_ == "reorder" && use_segregation_split_) {
// Again, unfortunate but safe use of dynamic_cast.
// Possible solution: refactor gravity solver to its own class.
dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
.solveGravity(&initial_porevol[0], stepsize, state);
}
watercut.push(timer.currentTime() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells_) {
wellreport.push(props_, *wells_, state.saturation(),
timer.currentTime() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
sreport.transport_time = tt;
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
reportVolumes(std::cout,satvol, tot_porevol_init,
tot_injected, tot_produced,
injected, produced,
init_satvol);
sreport.total_time = step_timer.secsSinceStart();
if (output_) {
sreport.reportParam(tstep_os);
}
}
if (output_) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
if (transport_solver_type_ == "reorder") {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
outputWaterCut(watercut, output_dir_);
if (wells_) {
outputWellReport(wellreport, output_dir_);
}
tstep_os.close();
}
total_timer.stop();
SimulatorReport report;
report.pressure_time = ptime;
report.transport_time = ttime;
report.total_time = total_timer.secsSinceStart();
return report;
}
} // namespace Opm

View File

@@ -0,0 +1,99 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED
#define OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED
#include <boost/shared_ptr.hpp>
#include <vector>
struct UnstructuredGrid;
struct Wells;
struct FlowBoundaryConditions;
namespace Opm
{
namespace parameter { class ParameterGroup; }
class IncompPropertiesInterface;
class RockCompressibility;
class WellsManager;
class LinearSolverInterface;
class SimulatorTimer;
class TwophaseState;
class WellState;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation.
class SimulatorIncompTwophaseAdfi
{
public:
/// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
/// parameter (default) effect
/// -----------------------------------------------------------
/// output (true) write output to files?
/// output_dir ("output") output directoty
/// output_interval (1) output every nth step
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal)
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure
/// nl_maxiter (30) max nonlinear iterations in transport
/// nl_tolerance (1e-9) transport solver absolute residual tolerance
/// num_transport_substeps (1) number of transport steps per pressure step
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] src source terms
/// \param[in] bcs boundary conditions, treat as all noflow if null
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorIncompTwophaseAdfi(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will
/// modify the reservoir and well states.
/// \param[in,out] timer governs the requested reporting timesteps
/// \param[in,out] state state of reservoir: pressure, fluxes
/// \param[in,out] well_state state of wells: bhp, perforation rates
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
TwophaseState& state,
WellState& well_state);
private:
class Impl;
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
boost::shared_ptr<Impl> pimpl_;
};
} // namespace Opm
#endif // OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED

View File

@@ -0,0 +1,254 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/autodiff/TransportSolverTwophaseAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver for Newton-Raphson scheme.
/// \param[in] gravity Gravity vector (null for no gravity).
/// \param[in] param Parameters for the solver.
TransportSolverTwophaseAd::TransportSolverTwophaseAd(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const LinearSolverInterface& linsolver,
const double* gravity,
const parameter::ParameterGroup& param)
: grid_(grid),
props_(props),
linsolver_(linsolver),
ops_(grid),
gravity_(0.0),
tol_(param.getDefault("nl_tolerance", 1e-9)),
maxit_(param.getDefault("nl_maxiter", 30))
{
const int nc = grid_.number_of_cells;
allcells_.resize(nc);
for (int i = 0; i < nc; ++i) {
allcells_[i] = i;
}
if (gravity && gravity[grid_.dimensions - 1] != 0.0) {
gravity_ = gravity[grid_.dimensions - 1];
for (int dd = 0; dd < grid_.dimensions - 1; ++dd) {
if (gravity[dd] != 0.0) {
THROW("TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
}
}
V htrans(grid.cell_facepos[grid.number_of_cells]);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
V trans(grid_.number_of_faces);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans.data());
transi_ = subset(trans, ops_.internal_faces);
}
}
// Virtual destructor.
TransportSolverTwophaseAd::~TransportSolverTwophaseAd()
{
}
namespace
{
template <class ADB>
std::vector<ADB>
phaseMobility(const Opm::IncompPropertiesInterface& props,
const std::vector<int>& cells,
const typename ADB::V& sw)
{
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
typedef Eigen::Array<double, Eigen::Dynamic, 4, Eigen::RowMajor> FourCol;
typedef typename ADB::V V;
typedef typename ADB::M M;
const int nc = props.numCells();
TwoCol s(nc, 2);
s.leftCols<1>() = sw;
s.rightCols<1>() = 1.0 - s.leftCols<1>();
TwoCol kr(nc, 2);
FourCol dkr(nc, 4);
props.relperm(nc, s.data(), cells.data(), kr.data(), dkr.data());
V krw = kr.leftCols<1>();
V kro = kr.rightCols<1>();
// In dkr, columns col(0..3) are:
// dkrw/dsw dkro/dsw dkrw/dso dkrw/dso <-- partial derivatives, really.
// If we want the derivatives with respect to some variable x,
// we must apply the chain rule:
// dkrw/dx = dkrw/dsw*dsw/dx + dkrw/dso*dso/dx.
// If x is sw as in our case we are left with.
// dkrw/dsw = col(0) - col(2)
// dkro/dsw = col(1) - col(3)
V dkrw = dkr.leftCols<1>() - dkr.rightCols<2>().leftCols<1>();
V dkro = dkr.leftCols<2>().rightCols<1>() - dkr.rightCols<1>();
M krwjac(nc,nc);
M krojac(nc,nc);
auto sizes = Eigen::ArrayXi::Ones(nc);
krwjac.reserve(sizes);
krojac.reserve(sizes);
for (int c = 0; c < nc; ++c) {
krwjac.insert(c,c) = dkrw(c);
krojac.insert(c,c) = dkro(c);
}
const double* mu = props.viscosity();
std::vector<M> dmw = { krwjac/mu[0] };
std::vector<M> dmo = { krojac/mu[1] };
std::vector<ADB> pmobc = { ADB::function(krw / mu[0], dmw) ,
ADB::function(kro / mu[1], dmo) };
return pmobc;
}
/// Returns fw(sw).
template <class ADB>
ADB
fluxFunc(const std::vector<ADB>& m)
{
assert (m.size() == 2);
ADB f = m[0] / (m[0] + m[1]);
return f;
}
} // anonymous namespace
/// Solve for saturation at next timestep.
/// Note that this only performs advection by total velocity, and
/// no gravity segregation.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term. For interpretation see Opm::computeTransportSource().
/// \param[in] dt Time step.
/// \param[in, out] state Reservoir state. Calling solve() will read state.faceflux() and
/// read and write state.saturation().
void TransportSolverTwophaseAd::solve(const double* porevolume,
const double* source,
const double dt,
TwophaseState& state)
{
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
typedef Eigen::Map<const V> Vec;
const int nc = grid_.number_of_cells;
const TwoCol s0 = Eigen::Map<const TwoCol>(state.saturation().data(), nc, 2);
double res_norm = 1e100;
const V sw0 = s0.leftCols<1>();
// sw1 is the object that will be changed every Newton iteration.
// V sw1 = sw0;
V sw1 = 0.5*V::Ones(nc,1);
const V dflux_all = Vec(state.faceflux().data(), grid_.number_of_faces, 1);
const int num_internal = ops_.internal_faces.size();
V dflux = subset(dflux_all, ops_.internal_faces);
// Upwind selection of mobilities by phase.
// We have that for a phase P
// v_P = lambda_P K (-grad p + rho_P g grad z)
// and we assume that this has the same direction as
// dh_P = -grad p + rho_P g grad z.
// This may not be true for arbitrary anisotropic situations,
// but for scalar lambda and using TPFA it holds.
const V p1 = Vec(state.pressure().data(), nc, 1);
const V ndp = (ops_.ngrad * p1.matrix()).array();
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DynArr;
const V z = Eigen::Map<DynArr>(grid_.cell_centroids, nc, grid_.dimensions).rightCols<1>();
const V ndz = (ops_.ngrad * z.matrix()).array();
ASSERT(num_internal == ndp.size());
const double* density = props_.density();
const V dhw = ndp - ndz*(gravity_*density[0]);
const V dho = ndp - ndz*(gravity_*density[1]);
const UpwindSelector<double> upwind_w(grid_, ops_, dhw);
const UpwindSelector<double> upwind_o(grid_, ops_, dho);
// Compute more explicit and constant terms used in the equations.
const V pv = Vec(porevolume, nc, 1);
const V dtpv = dt/pv;
const V q = Vec(source, nc, 1);
const V qneg = q.min(V::Zero(nc,1));
const V qpos = q.max(V::Zero(nc,1));
const double gfactor = gravity_*(density[0] - density[1]);
const V gravflux = (gravity_ == 0.0) ? V(V::Zero(num_internal, 1))
: ndz*transi_*gfactor;
// Block pattern for variables.
// Primary variables:
// sw : one per cell
std::vector<int> bpat = { nc };
// Newton-Raphson loop.
int it = 0;
do {
// Assemble linear system.
const ADB sw = ADB::variable(0, sw1, bpat);
const std::vector<ADB> pmobc = phaseMobility<ADB>(props_, allcells_, sw.value());
const ADB fw_cell = fluxFunc(pmobc);
const std::vector<ADB> pmobf = { upwind_w.select(pmobc[0]),
upwind_o.select(pmobc[1]) };
const ADB fw_face = fluxFunc(pmobf);
const ADB flux = fw_face * (dflux - pmobf[1]*gravflux);
// const ADB fw_face = upwind_w.select(fw_cell);
// const ADB flux = fw_face * dflux;
const ADB qtr_ad = qpos + fw_cell*qneg;
const ADB transport_residual = sw - sw0 + dtpv*(ops_.div*flux - qtr_ad);
res_norm = transport_residual.value().matrix().norm();
std::cout << "Residual l2-norm = " << res_norm << std::endl;
// Solve linear system.
Eigen::SparseMatrix<double, Eigen::RowMajor> smatr = transport_residual.derivative()[0];
ASSERT(smatr.isCompressed());
V ds(nc);
LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(nc, smatr.nonZeros(),
smatr.outerIndexPtr(), smatr.innerIndexPtr(), smatr.valuePtr(),
transport_residual.value().data(), ds.data());
if (!rep.converged) {
THROW("Linear solver convergence error in TransportSolverTwophaseAd::solve()");
}
// Update (possible clamp) sw1.
sw1 = sw.value() - ds;
sw1 = sw1.min(V::Ones(nc,1)).max(V::Zero(nc,1));
it += 1;
} while (res_norm > tol_ && it < maxit_);
// Write to output data structure.
Eigen::Map<TwoCol> sref(state.saturation().data(), nc, 2);
sref.leftCols<1>() = sw1;
sref.rightCols<1>() = 1.0 - sw1;
}
} // namespace Opm

View File

@@ -0,0 +1,89 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTSOLVERTWOPHASEAD_HEADER_INCLUDED
#define OPM_TRANSPORTSOLVERTWOPHASEAD_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/transport/TransportSolverTwophaseInterface.hpp>
#include <vector>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
class LinearSolverInterface;
namespace parameter { class ParameterGroup; }
/// Implements an implicit transport solver for incompressible two-phase flow,
/// using automatic differentiation.
class TransportSolverTwophaseAd : public TransportSolverTwophaseInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver for Newton-Raphson scheme.
/// \param[in] gravity Gravity vector (null for no gravity).
/// \param[in] param Parameters for the solver.
TransportSolverTwophaseAd(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const LinearSolverInterface& linsolver,
const double* gravity,
const parameter::ParameterGroup& param);
// Virtual destructor.
virtual ~TransportSolverTwophaseAd();
/// Solve for saturation at next timestep.
/// Note that this only performs advection by total velocity, and
/// no gravity segregation.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term. For interpretation see Opm::computeTransportSource().
/// \param[in] dt Time step.
/// \param[in, out] state Reservoir state. Calling solve() will read state.faceflux() and
/// read and write state.saturation().
virtual void solve(const double* porevolume,
const double* source,
const double dt,
TwophaseState& state);
private:
typedef AutoDiff::ForwardBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
private:
const UnstructuredGrid& grid_;
const IncompPropertiesInterface& props_;
const LinearSolverInterface& linsolver_;
const HelperOps ops_;
double gravity_;
double tol_;
int maxit_;
std::vector<int> allcells_;
V transi_;
};
} // namespace Opm
#endif // OPM_TRANSPORTSOLVERTWOPHASEAD_HEADER_INCLUDED