OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2 from atgeirr/master

Browse Source 
Multiple changes and advances
This commit is contained in:
Bård Skaflestad
2013-05-22 01:17:51 -07:00
parent f6b46bda8f f2f83fee32
commit 04e51491ee
9 changed files with 444 additions and 653 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -64,7 +64,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/AutoDiffBlock.hpp
opm/autodiff/AutoDiffHelpers.hpp
opm/autodiff/AutoDiff.hpp
opm/autodiff/AutoDiffVec.hpp
opm/autodiff/BlackoilPropsAd.hpp
opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/SimulatorIncompTwophaseAdfi.hpp

22
cmake/Modules/FindEigen3.cmake
View File

@@ -7,7 +7,7 @@
# Once done this will define
#
# EIGEN3_FOUND - system has eigen lib with correct version
# EIGEN3_INCLUDE_DIRS - the eigen include directory
# EIGEN3_INCLUDE_DIR - the eigen include directory
# EIGEN3_VERSION - eigen version
# Copyright (c) 2006, 2007 Montel Laurent, <montel@kde.org>
@@ -30,7 +30,7 @@ if(NOT Eigen3_FIND_VERSION)
endif(NOT Eigen3_FIND_VERSION)
macro(_eigen3_check_version)
file(READ "${EIGEN3_INCLUDE_DIRS}/Eigen/src/Core/util/Macros.h" _eigen3_version_header)
file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header)
string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}")
set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}")
@@ -48,34 +48,34 @@ macro(_eigen3_check_version)
if(NOT EIGEN3_VERSION_OK)
message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in ${EIGEN3_INCLUDE_DIRS}, "
message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in ${EIGEN3_INCLUDE_DIR}, "
"but at least version ${Eigen3_FIND_VERSION} is required")
endif(NOT EIGEN3_VERSION_OK)
endmacro(_eigen3_check_version)
if (EIGEN3_INCLUDE_DIRS)
if (EIGEN3_INCLUDE_DIR)
# in cache already
_eigen3_check_version()
set(EIGEN3_FOUND ${EIGEN3_VERSION_OK})
else (EIGEN3_INCLUDE_DIRS)
else (EIGEN3_INCLUDE_DIR)
find_path(EIGEN3_INCLUDE_DIRS NAMES signature_of_eigen3_matrix_library
find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library
PATHS
${CMAKE_INSTALL_PREFIX}/include
${KDE4_INCLUDE_DIR}
PATH_SUFFIXES eigen3 eigen
)
if(EIGEN3_INCLUDE_DIRS)
if(EIGEN3_INCLUDE_DIR)
_eigen3_check_version()
endif(EIGEN3_INCLUDE_DIRS)
endif(EIGEN3_INCLUDE_DIR)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(Eigen3 DEFAULT_MSG EIGEN3_INCLUDE_DIRS EIGEN3_VERSION_OK)
find_package_handle_standard_args(Eigen3 DEFAULT_MSG EIGEN3_INCLUDE_DIR EIGEN3_VERSION_OK)
mark_as_advanced(EIGEN3_INCLUDE_DIRS)
mark_as_advanced(EIGEN3_INCLUDE_DIR)
endif(EIGEN3_INCLUDE_DIRS)
endif(EIGEN3_INCLUDE_DIR)

6
cmake/Modules/OpmFind.cmake
View File

@@ -158,6 +158,12 @@ macro (find_and_append_package_to prefix name)
elseif (DEFINED ${NAME}_INCLUDE_PATH)
list (APPEND ${prefix}_INCLUDE_DIRS ${${NAME}_INCLUDE_PATH})
endif (DEFINED ${name}_INCLUDE_PATH)
# some packages define only _DIR and not _DIRS (Hi, Eigen3!)
if (DEFINED ${name}_INCLUDE_DIR)
list (APPEND ${prefix}_INCLUDE_DIRS ${${name}_INCLUDE_DIR})
elseif (DEFINED ${NAME}_INCLUDE_DIR)
list (APPEND ${prefix}_INCLUDE_DIRS ${${NAME}_INCLUDE_DIR})
endif (DEFINED ${name}_INCLUDE_DIR)
endif ("${var}" STREQUAL "INCLUDE_DIRS")
# cleanup lists
if ("${var}" STREQUAL "LIBRARIES")

45
examples/test_impestpfa_ad.cpp
View File

@@ -20,7 +20,10 @@
#include <config.h>
#define HACK_INCOMPRESSIBLE_GRAVITY 0
#include <opm/autodiff/ImpesTPFAAD.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
@@ -97,51 +100,65 @@ namespace {
};
}
int
main(int argc, char* argv[])
{
const Opm::parameter::ParameterGroup param(argc, argv, false);
const Opm::GridManager gm(20, 1);
const Opm::GridManager gm(5, 5);
const UnstructuredGrid* g = gm.c_grid();
const int nc = g->number_of_cells;
const Opm::BlackoilPropertiesBasic props(param, 2, nc);
const Opm::BlackoilPropertiesBasic oldprops(param, 2, nc);
const Opm::BlackoilPropsAd props(oldprops);
typedef AutoDiff::ForwardBlock<double> ADB;
typedef Opm::BlackoilPropertiesInterface Geology;
typedef DerivedGeology<Geology, ADB::V> GeoProps;
typedef Opm::BlackoilPropertiesInterface BOFluid;
typedef Opm::BlackoilPropsAd BOFluid;
typedef Opm::ImpesTPFAAD<BOFluid, GeoProps> PSolver;
Wells* wells = create_wells(2, 2, 2);
Wells* wells = create_wells(2, 2, 5);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int inj_cell = 0;
const int prod_cell = g->number_of_cells - 1;
const double WI = 1e-8;
bool ok = add_well(INJECTOR, 0.0, 1, inj_frac, &inj_cell, &WI, "Inj", wells);
ok = ok && add_well(PRODUCER, 0.0, 1, prod_frac, &prod_cell, &WI, "Prod", wells);
const int inj_cells[] = { 0, 1, 2 };
const int prod_cells[] = { 20, 21 };
const double WI[3] = { 1e-14 };
bool ok = add_well(INJECTOR, 0.0, 1, inj_frac, inj_cells, WI, "Inj", wells);
ok = ok && add_well(PRODUCER, 0.0, 1, prod_frac, prod_cells, WI, "Prod", wells);
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells);
ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 8.64297e-05, oildistr, 1, wells);
if (!ok) {
THROW("Something went wrong with well init.");
}
set_current_control(0, 0, wells);
set_current_control(1, 0, wells);
double grav[] = { 1.0, 0.0 };
GeoProps geo(*g, props, grav);
double grav[] = { /*1.0*/ 0.0, 0.0 };
GeoProps geo(*g, oldprops, grav);
Opm::LinearSolverFactory linsolver(param);
PSolver ps (*g, props, geo, *wells, linsolver);
Opm::BlackoilState state;
initStateBasic(*g, props, param, 0.0, state);
initBlackoilSurfvol(*g, props, state);
initStateBasic(*g, oldprops, param, 0.0, state);
initBlackoilSurfvol(*g, oldprops, state);
Opm::WellState well_state;
well_state.init(wells, state);
ps.solve(1.0, state, well_state);
std::cout << "Cell pressure:" << std::endl;
std::copy(state.pressure().begin(), state.pressure().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Face flux:" << std::endl;
std::copy(state.faceflux().begin(), state.faceflux().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Well bhp pressure:" << std::endl;
std::copy(well_state.bhp().begin(), well_state.bhp().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;

216
opm/autodiff/AutoDiffBlockArma.hpp
View File

@@ -1,216 +0,0 @@
/*
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

129
opm/autodiff/AutoDiffHelpers.hpp
View File

@@ -100,8 +100,8 @@ struct HelperOps
#include <cstdio>
#endif // !defined(NDEBUG)
namespace {
#if !defined(NDEBUG)
namespace {
void
printSparseMatrix(const Eigen::SparseMatrix<double>& A,
std::FILE* fp)
@@ -134,6 +134,7 @@ namespace {
std::fclose(fp);
}
} // anonymous namespace
#endif // !defined(NDEBUG)
@@ -209,7 +210,8 @@ namespace {
private:
typename ADB::M select_;
};
}
namespace {
@@ -279,7 +281,7 @@ superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices,
const int n)
{
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x;
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
}
@@ -303,4 +305,125 @@ spdiag(const AutoDiff::ForwardBlock<double>::V& d)
return mat;
}
/// Selection. Choose first of two elements if selection basis element is nonnegative.
template <typename Scalar>
class Selector {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
Selector(const typename ADB::V& selection_basis)
{
// Define selector structure.
const int n = selection_basis.size();
// Over-reserving so we do not have to count.
left_elems_.reserve(n);
right_elems_.reserve(n);
for (int i = 0; i < n; ++i) {
if (selection_basis[i] < 0.0) {
right_elems_.push_back(i);
} else {
left_elems_.push_back(i);
}
}
}
/// Apply selector to ADB quantities.
ADB select(const ADB& x1, const ADB& x2) const
{
if (right_elems_.empty()) {
return x1;
} else if (left_elems_.empty()) {
return x2;
} else {
return superset(subset(x1, left_elems_), left_elems_, x1.size())
+ superset(subset(x2, right_elems_), right_elems_, x2.size());
}
}
/// Apply selector to ADB quantities.
typename ADB::V select(const typename ADB::V& x1, const typename ADB::V& x2) const
{
if (right_elems_.empty()) {
return x1;
} else if (left_elems_.empty()) {
return x2;
} else {
return superset(subset(x1, left_elems_), left_elems_, x1.size())
+ superset(subset(x2, right_elems_), right_elems_, x2.size());
}
}
private:
std::vector<int> left_elems_;
std::vector<int> right_elems_;
};
/// Returns the input expression, but with all Jacobians collapsed to one.
inline
AutoDiff::ForwardBlock<double>
collapseJacs(const AutoDiff::ForwardBlock<double>& x)
{
typedef AutoDiff::ForwardBlock<double> ADB;
const int nb = x.numBlocks();
typedef Eigen::Triplet<double> Tri;
int nnz = 0;
for (int block = 0; block < nb; ++block) {
nnz += x.derivative()[block].nonZeros();
}
std::vector<Tri> t;
t.reserve(nnz);
int block_col_start = 0;
for (int block = 0; block < nb; ++block) {
const ADB::M& jac = x.derivative()[block];
// ADB::M is column major
for (int col = 0; col < jac.cols(); ++col) {
for (int elem = jac.outerIndexPtr()[col];
elem < jac.outerIndexPtr()[col + 1];
++elem) {
const int row = jac.innerIndexPtr()[elem];
t.emplace_back(row, block_col_start + col, jac.valuePtr()[elem]);
}
}
block_col_start += jac.cols();
}
// Build final jacobian.
std::vector<ADB::M> jacs(1);
jacs[0].resize(x.size(), block_col_start);
jacs[0].setFromTriplets(t.begin(), t.end());
return ADB::function(x.value(), jacs);
}
/// Returns the vertical concatenation [ x; y ] of the inputs.
inline
AutoDiff::ForwardBlock<double>
vertcat(const AutoDiff::ForwardBlock<double>& x,
const AutoDiff::ForwardBlock<double>& y)
{
const int nx = x.size();
const int ny = y.size();
const int n = nx + ny;
std::vector<int> xind(nx);
for (int i = 0; i < nx; ++i) {
xind[i] = i;
}
std::vector<int> yind(ny);
for (int i = 0; i < ny; ++i) {
yind[i] = nx + i;
}
return superset(x, xind, n) + superset(y, yind, n);
}
#endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

134
opm/autodiff/AutoDiffVec.hpp
View File

@@ -1,134 +0,0 @@
/*
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

489
opm/autodiff/ImpesTPFAAD.hpp
View File

@@ -33,6 +33,7 @@
#include <algorithm>
#include <cassert>
#include <vector>
#include <iterator>
#include <boost/shared_ptr.hpp>
@@ -101,154 +102,6 @@ namespace {
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:
@@ -258,14 +111,17 @@ namespace Opm {
const Wells& wells,
const LinearSolverInterface& linsolver)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
, wells_ (wells)
, linsolver_(linsolver)
, pdepfdata_(grid.number_of_cells, fluid)
// , pdepfdata_(grid.number_of_cells, fluid)
, ops_ (grid)
, grav_ (gravityOperator(grid_, ops_, geo_))
, cell_residual_ (ADB::null())
, well_flow_residual_ ()
, well_residual_ (ADB::null())
, total_residual_ (ADB::null())
{
}
@@ -274,7 +130,28 @@ namespace Opm {
BlackoilState& state,
WellState& well_state)
{
pdepfdata_.computeSatQuant(state);
const int nc = grid_.number_of_cells;
const int np = state.numPhases();
well_flow_residual_.resize(np, ADB::null());
// Compute relperms once and for all (since saturations are explicit).
DataBlock s = Eigen::Map<const DataBlock>(state.saturation().data(), nc, np);
ASSERT(np == 2);
kr_ = fluid_.relperm(s.col(0), s.col(1), V::Zero(nc,1), buildAllCells(nc));
// Compute relperms for wells. This must be revisited for crossflow.
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
DataBlock well_s(nperf, np);
for (int w = 0; w < nw; ++w) {
const double* comp_frac = &wells_.comp_frac[np*w];
for (int j = wells_.well_connpos[w]; j < wells_.well_connpos[w+1]; ++j) {
well_s.row(j) = Eigen::Map<const DataBlock>(comp_frac, 1, np);
}
}
const std::vector<int> well_cells(wells_.well_cells,
wells_.well_cells + nperf);
well_kr_ = fluid_.relperm(well_s.col(0), well_s.col(1), V::Zero(nperf,1), well_cells);
const double atol = 1.0e-15;
const double rtol = 5.0e-10;
@@ -286,7 +163,7 @@ namespace Opm {
int it = 0;
bool resTooLarge = r0 > atol;
while (resTooLarge && (it < maxit)) {
solveJacobianSystem(state);
solveJacobianSystem(state, well_state);
assemble(dt, state, well_state);
@@ -310,55 +187,67 @@ namespace Opm {
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;
// typedef PressureDependentFluidData<double, BOFluid> PDepFData;
// typedef typename PDepFData::ADB ADB;
typedef AutoDiff::ForwardBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
const UnstructuredGrid& grid_;
const GeoProps& geo_ ;
const Wells& wells_;
const UnstructuredGrid& grid_;
const BOFluid& fluid_;
const GeoProps& geo_ ;
const Wells& wells_;
const LinearSolverInterface& linsolver_;
PDepFData pdepfdata_;
HelperOps ops_;
const M grav_;
ADB cell_residual_;
ADB well_residual_;
HelperOps ops_;
const M grav_;
ADB cell_residual_;
std::vector<ADB> well_flow_residual_;
ADB well_residual_;
ADB total_residual_;
std::vector<V> kr_;
std::vector<V> well_kr_;
enum { Water = BOFluid::Water,
Oil = BOFluid::Oil,
Gas = BOFluid::Gas };
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;
const int nperf = wells_.well_connpos[nw];
pdepfdata_.computePressQuant(state);
const std::vector<int> cells = buildAllCells(nc);
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);
wells_.well_cells + nperf);
const V transw = Eigen::Map<const V>(wells_.WI, nperf, 1);
// Initialize AD variables: p (cell pressures) and bhp (well bhp).
// Initialize AD variables: p (cell pressures), qs (well rates) and bhp (well bhp).
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V qs0 = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np, 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);
std::vector<V> vars0 = { p0, qs0, bhp0 };
std::vector<ADB> vars = ADB::variables(vars0);
const ADB& p = vars[0];
const ADB& bhp = vars[1];
const ADB& qs = vars[1];
const ADB& bhp = vars[2];
std::vector<int> bpat = p.blockPattern();
// Compute T_ij * (p_i - p_j) and use for upwinding.
@@ -378,108 +267,270 @@ namespace Opm {
}
}
well_to_perf.setFromTriplets(w2p.begin(), w2p.end());
const M perf_to_well = well_to_perf.transpose();
// Construct pressure difference vector for wells.
const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet!
// Finally construct well perforation pressures.
// Finally construct well perforation pressures and well flows.
const ADB p_perfwell = well_to_perf*bhp + well_perf_dp;
const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
const Selector<double> cell_to_well_selector(nkgradp_well.value());
cell_residual_ = ADB::constant(pv, bpat);
#define COMPENSATE_FLOAT_PRECISION 0
#if COMPENSATE_FLOAT_PRECISION
well_residual_ = ADB::constant(V::Zero(nw,1), bpat);
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 ADB cell_b = fluidFvf(phase, p, cells);
const ADB cell_rho = fluidRho(phase, p, cells);
const ADB well_b = fluidFvf(phase, p_perfwell, well_cells);
const V kr = fluidKr(phase);
// Explicitly not asking for derivatives of viscosity,
// since they are not available yet.
const V mu = fluidMu(phase, p.value(), cells);
const V cell_mob = kr / mu;
const V face_mob = upwind.select(cell_mob);
const V well_kr = fluidKrWell(phase);
const V well_mu = fluidMu(phase, p_perfwell.value(), well_cells);
const V well_mob = well_kr / well_mu;
const V perf_mob = cell_to_well_selector.select(subset(cell_mob, well_cells), well_mob);
const ADB flux = face_mob * (nkgradp + (grav_ * cell_rho));
const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations.
const ADB face_b = upwind.select(cell_b);
const ADB perf_b = cell_to_well_selector.select(subset(cell_b, well_cells), well_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 ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
const ADB divcontrib = delta_t * (ops_.div * (flux * face_b)
+ well_contrib);
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
divcontrib_sum = divcontrib_sum - divcontrib/cell_b;
cell_residual_ = cell_residual_ - (component_contrib/cell_b);
const ADB well_rates = perf_to_well * (perf_flux*perf_b);
std::vector<int> well_flow_res_phase_idx(nw);
for (int w = 0; w < nw; ++w) {
well_flow_res_phase_idx[w] = w + phase*nw;
}
well_flow_residual_[phase] = well_rates - subset(qs, well_flow_res_phase_idx);
}
#if COMPENSATE_FLOAT_PRECISION
cell_residual_ = cell_residual_ + divcontrib_sum;
#endif
// Handling BHP and SURFACE_RATE wells.
V bhp_targets(nw,1);
V rate_targets(nw,1);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
if (wc->type[wc->current] == BHP) {
bhp_targets[w] = wc->target[wc->current];
rate_targets[w] = -1e100;
} else if (wc->type[wc->current] == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = wc->target[wc->current];
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = wc->distr[phase];
}
} else {
THROW("Can only handle BHP and SURFACE_RATE type controls.");
}
}
const ADB bhp_residual = bhp - bhp_targets;
const ADB rate_residual = rate_distr * qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
well_residual_ = bhp_selector.select(bhp_residual, rate_residual);
ASSERT(np == 2);
const ADB well_flow_res = vertcat(well_flow_residual_[0], well_flow_residual_[1]);
const ADB well_res = vertcat(well_flow_res, well_residual_);
total_residual_ = collapseJacs(vertcat(cell_residual_, well_res));
}
void
solveJacobianSystem(BlackoilState& state) const
solveJacobianSystem(BlackoilState& state,
WellState& well_state) const
{
const int nc = grid_.number_of_cells;
Eigen::SparseMatrix<double, Eigen::RowMajor> matr = cell_residual_.derivative()[0];
const int nw = wells_.number_of_wells;
const int np = state.numPhases();
Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_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);
V dx(total_residual_.size());
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(nc, matr.nonZeros(),
= linsolver_.solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
cell_residual_.value().data(), dp.data());
total_residual_.value().data(), dx.data());
if (!rep.converged) {
THROW("ImpesTPFAAD::solve(): Linear solver convergence failure.");
}
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V dp = subset(dx, buildAllCells(nc));
const V p = p0 - dp;
std::copy(&p[0], &p[0] + nc, state.pressure().begin());
const V qs0 = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np, 1);
std::vector<int> qs_dofs(np*nw);
for (int w = 0; w < nw; ++w) {
for (int phase = 0; phase < np; ++phase) {
qs_dofs[w*np + phase] = nc + w*np + phase;
}
}
const V dqs = subset(dx, qs_dofs);
const V qs = qs0 - dqs;
std::copy(&qs[0], &qs[0] + np*nw, well_state.wellRates().begin());
const V bhp0 = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
std::vector<int> bhp_dofs(nw);
for (int w = 0; w < nw; ++w) {
bhp_dofs[w] = nc + np*nw + w;
}
ASSERT(bhp_dofs.back() + 1 == total_residual_.size());
const V dbhp = subset(dx, bhp_dofs);
const V bhp = bhp0 - dbhp;
std::copy(&bhp[0], &bhp[0] + nw, well_state.bhp().begin());
}
double
residualNorm() const
{
return cell_residual_.value().matrix().norm();
return total_residual_.value().matrix().norm();
}
void
computeFluxes(BlackoilState& state) const
{
const int nc = grid_.number_of_cells;
const int np = state.numPhases();
const std::vector<int> cells = buildAllCells(nc);
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const ADB p = ADB::constant(p0, std::vector<int>(1, nc));
const ADB p = ADB::constant(p0, cell_residual_.blockPattern());
const V transi = subset(geo_.transmissibility(),
ops_.internal_faces);
const V nkgradp = transi * (ops_.ngrad * p0.matrix()).array();
V flux = Eigen::Map<const V>(&state.faceflux()[0],
grid_.number_of_faces, 1);
V flux = V::Zero(ops_.internal_faces.size(), 1);
for (int phase = 0; phase < np; ++phase) {
const ADB cell_rho = pdepfdata_.phaseDensity(phase, p);
const V head = nkgradp +
(grav_ * cell_rho.value().matrix()).array();
const ADB cell_rho = fluidRho(phase, p, cells);
const V head = nkgradp + (grav_ * cell_rho.value().matrix()).array();
const UpwindSelector<double> upwind(grid_, ops_, head);
const V kr = pdepfdata_.phaseRelPerm(phase);
const ADB mu = pdepfdata_.phaseViscosity(phase, p);
const V mf = upwind.select(kr / mu.value());
const V kr = fluidKr(phase);
const V mu = fluidMu(phase, p.value(), cells);
const V mf = upwind.select(kr / mu);
flux += mf * head;
}
V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces);
std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().data());
}
V fluidMu(const int phase, const V& p, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
case Oil: {
V dummy_rs = V::Zero(p.size(), 1) * p;
return fluid_.muOil(p, dummy_rs, cells);
}
case Gas:
return fluid_.muGas(p, cells);
default:
THROW("Unknown phase index " << phase);
}
}
ADB fluidMu(const int phase, const ADB& p, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
case Oil: {
ADB dummy_rs = V::Zero(p.size(), 1) * p;
return fluid_.muOil(p, dummy_rs, cells);
}
case Gas:
return fluid_.muGas(p, cells);
default:
THROW("Unknown phase index " << phase);
}
}
ADB fluidFvf(const int phase, const ADB& p, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
case Oil: {
ADB dummy_rs = V::Zero(p.size(), 1) * p;
return fluid_.bOil(p, dummy_rs, cells);
}
case Gas:
return fluid_.bGas(p, cells);
default:
THROW("Unknown phase index " << phase);
}
}
ADB fluidRho(const int phase, const ADB& p, const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidFvf(phase, p, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
return rho;
}
V fluidKr(const int phase) const
{
return kr_[phase];
}
V fluidKrWell(const int phase) const
{
return well_kr_[phase];
}
};
} // namespace Opm

55
tests/not-unit/test_vec.cpp
View File

@@ -1,55 +0,0 @@
/*
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/AutoDiffVec.hpp>
#include <iostream>
int main()
{
typedef AutoDiff::ForwardVec<double> ADV;
ADV::V v(3);
v << 1.0, 2.2, 3.4;
ADV::V v2(3);
v2 << 0.2, 1.2, 13.4;
// std::cout << v << std::endl;
ADV a = ADV::constant(v2);
ADV x = ADV::variable(v);
ADV::M jac(3,3);
jac.insert(0, 0) = 1.0;
ADV f = ADV::function(v, jac);
ADV xpx = x + x;
std::cout << xpx;
ADV xpxpa = x + x + a;
std::cout << xpxpa;
std::cout << xpxpa - xpx;
ADV sqx = x * x;
std::cout << sqx;
ADV sqxdx = sqx / x;
std::cout << sqxdx;
// std::cout << a << "\n\n" << x << "\n\n" << f << std::endl;
}