OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into ert

Browse Source 
Conflicts:
	Makefile.am
This commit is contained in:
Atgeirr Flø Rasmussen 2012-09-26 09:07:57 +02:00
commit 0ceff04194
17 changed files with 745 additions and 151 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -29,6 +29,7 @@ ltmain.sh
m4/libtool.m4
m4/lt*.m4
missing
ar-lib
tutorials/tutorial[1-4]
# Ignoring executables
@ -38,6 +39,7 @@ examples/spu_2p
examples/reorder-qfs
examples/refine_wells
examples/sim_2p_incomp_reorder
examples/sim_2p_comp_reorder
examples/sim_wateroil
examples/wells_example
tests/test_cfs_tpfa

79
Makefile.am
View File

@ -38,22 +38,23 @@ $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
# support below for additional details.
nodist_lib_libopmcore_la_SOURCES =
lib_libopmcore_la_SOURCES = \
opm/core/GridManager.cpp \
opm/core/eclipse/EclipseGridInspector.cpp \
opm/core/eclipse/EclipseGridParser.cpp \
opm/core/fluid/BlackoilPropertiesBasic.cpp \
opm/core/fluid/BlackoilPropertiesFromDeck.cpp \
opm/core/fluid/IncompPropertiesBasic.cpp \
opm/core/fluid/IncompPropertiesFromDeck.cpp \
opm/core/fluid/PvtPropertiesBasic.cpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp \
opm/core/fluid/RockBasic.cpp \
opm/core/fluid/RockCompressibility.cpp \
opm/core/fluid/RockFromDeck.cpp \
opm/core/fluid/SaturationPropsBasic.cpp \
opm/core/fluid/SaturationPropsFromDeck.cpp \
opm/core/fluid/SatFuncStone2.cpp \
lib_libopmcore_la_SOURCES = \
opm/core/GridManager.cpp \
opm/core/eclipse/EclipseGridInspector.cpp \
opm/core/eclipse/EclipseGridParser.cpp \
opm/core/fluid/BlackoilPropertiesBasic.cpp \
opm/core/fluid/BlackoilPropertiesFromDeck.cpp \
opm/core/fluid/IncompPropertiesBasic.cpp \
opm/core/fluid/IncompPropertiesFromDeck.cpp \
opm/core/fluid/PvtPropertiesBasic.cpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp \
opm/core/fluid/RockBasic.cpp \
opm/core/fluid/RockCompressibility.cpp \
opm/core/fluid/RockFromDeck.cpp \
opm/core/fluid/SaturationPropsBasic.cpp \
opm/core/fluid/SaturationPropsFromDeck.cpp \
opm/core/fluid/SatFuncGwseg.cpp \
opm/core/fluid/SatFuncStone2.cpp \
opm/core/fluid/SatFuncSimple.cpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \
opm/core/fluid/blackoil/SinglePvtDead.cpp \
@ -131,27 +132,28 @@ opm/core/wells/WellCollection.cpp \
opm/core/wells/WellsGroup.cpp \
opm/core/wells/WellsManager.cpp
nobase_include_HEADERS = \
opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/eclipse/CornerpointChopper.hpp \
opm/core/eclipse/EclipseGridInspector.hpp \
opm/core/eclipse/EclipseGridParser.hpp \
opm/core/eclipse/EclipseGridParserHelpers.hpp \
opm/core/eclipse/EclipseUnits.hpp \
opm/core/eclipse/SpecialEclipseFields.hpp \
opm/core/fluid/BlackoilPropertiesBasic.hpp \
opm/core/fluid/BlackoilPropertiesFromDeck.hpp \
opm/core/fluid/BlackoilPropertiesInterface.hpp \
opm/core/fluid/IncompPropertiesBasic.hpp \
opm/core/fluid/IncompPropertiesFromDeck.hpp \
opm/core/fluid/IncompPropertiesInterface.hpp \
opm/core/fluid/PvtPropertiesBasic.hpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \
opm/core/fluid/RockBasic.hpp \
opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SatFuncStone2.hpp \
nobase_include_HEADERS = \
opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/eclipse/CornerpointChopper.hpp \
opm/core/eclipse/EclipseGridInspector.hpp \
opm/core/eclipse/EclipseGridParser.hpp \
opm/core/eclipse/EclipseGridParserHelpers.hpp \
opm/core/eclipse/EclipseUnits.hpp \
opm/core/eclipse/SpecialEclipseFields.hpp \
opm/core/fluid/BlackoilPropertiesBasic.hpp \
opm/core/fluid/BlackoilPropertiesFromDeck.hpp \
opm/core/fluid/BlackoilPropertiesInterface.hpp \
opm/core/fluid/IncompPropertiesBasic.hpp \
opm/core/fluid/IncompPropertiesFromDeck.hpp \
opm/core/fluid/IncompPropertiesInterface.hpp \
opm/core/fluid/PvtPropertiesBasic.hpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \
opm/core/fluid/RockBasic.hpp \
opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SatFuncGwseg.hpp \
opm/core/fluid/SatFuncStone2.hpp \
opm/core/fluid/SatFuncSimple.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
@ -316,3 +318,6 @@ opm/core/linalg/LinearSolverAGMG.hpp
lib_libopmcore_la_LDFLAGS += \
$(FCLIBS)
endif
pkgconfigdir = $(libdir)/pkgconfig
pkgconfig_DATA = lib/pkgconfig/opm-core.pc

51
README
View File

@ -46,11 +46,17 @@ sudo apt-get install -y build-essential gfortran pkg-config libtool \
sudo apt-get install -y doxygen ghostscript texlive-latex-recommended pgf
# packages necessary for version control
sudo apt-get install -y git-core git-svn subversion
sudo apt-get install -y git-core
# libraries necessary for DUNE
# basic libraries necessary for both DUNE and OPM
sudo apt-get install -y libboost-all-dev libsuperlu3-dev libsuitesparse-dev
# add this repository for necessary backports (required for Ubuntu Precise)
sudo add-apt-repository -y ppa:opm/ppa
# parts of DUNE needed
sudo apt-get install -y libdune-common-dev libdune-istl-dev
# libraries necessary for OPM
sudo apt-get install -y libxml0-dev
@ -59,33 +65,14 @@ DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
-----------------------------------------
# libraries
sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack
sudo zypper in blas libblas3 lapack liblapack3 libboost libxml2 umfpack
# tools
sudo zypper install gcc automake autoconf git doxygen
sudo zypper in gcc automake autoconf git doxygen
RETRIEVING AND BUILDING DUNE PREREQUISITES
------------------------------------------
(only necessary if you want to use opm-core as a dune module)
# trust DUNE certificate (sic)
echo p | svn list https://svn.dune-project.org/svn/dune-common
# checkout DUNE libraries
for module in common istl geometry grid localfunctions; do
git svn clone -s \
https://svn.dune-project.org/svn/dune-$module/branches/release-2.2/ \
dune-$module
done
# building DUNE libraries
for module in common istl geometry grid localfunctions; do
env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=dune-$module \
--configure-opts="--enable-fieldvector-size-is-method" \
--make-opts="-j -l 0.8" autogen : configure : make
done
# DUNE libraries
sudo zypper ar http://download.opensuse.org/repositories/science/openSUSE_12.2/science.repo
sudo zypper in dune-common dune-istl
DOWNLOADING
@ -100,20 +87,10 @@ If you want to contribute, fork OPM/opm-core on github.
BUILDING
--------
(standalone opm-core:)
cd ../opm-core
cd opm-core
autoreconf -i
./configure
make
sudo make install
(using opm-core as a dune module:)
# note: this is done from the parent directory of opm-core
env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \
--configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make
DOCUMENTATION

9
configure.ac
View File

@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
AC_INIT([OPM Core Library], [0.1], [atgeirr@sintef.no],dnl
AC_INIT([OPM Core Library], [0.1], [atgeirr@sintef.no],
[opmcore], [https://public.ict.sintef.no/opm/hg/opmcore])
AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects])
@ -20,8 +20,11 @@ AC_CONFIG_HEADERS([config.h])
AC_PROG_CC
AM_PROG_CC_C_O
dnl Initialize libtool; the funny indentation here is to
dnl satisfy libtoolize' check for the presence of this macro
m4_ifdef([LT_INIT],
[LT_INIT[]dnl
[
LT_INIT[]dnl
LT_LANG([C++])dnl
LT_LANG([Fortran 77])dnl
LT_LANG([Fortran])dnl
@ -53,6 +56,8 @@ AC_CONFIG_FILES([
tests/Makefile
examples/Makefile
tutorials/Makefile
opm-core.pc
lib/pkgconfig/opm-core.pc
])
AC_OUTPUT

11
lib/pkgconfig/opm-core.pc.in Normal file
View File

@ -0,0 +1,11 @@
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
includedir=@includedir@
Name: @PACKAGE_NAME@
Description: @PACKAGE_STRING@
Version: @PACKAGE_VERSION@
URL: @PACKAGE_URL@
Libs: -L${libdir} -l@PACKAGE@
Cflags: -I${includedir}

18
opm-core.pc.in Normal file
View File

@ -0,0 +1,18 @@
# This is the configuration for local builds. Use this by putting the
# compilation output path (the directory in which you ran ./configure)
# into the environment variable PKG_CONFIG_PATH. This will enable you
# to use pkg-config in your code while making changes to opm-core.
# This is NOT the file that is installed in the system directories when
# you do `make install`. That is the one in lib/pkgconfig. However, if
# you make changes here, you should consider that one as well.
libdir=@abs_top_builddir@/lib/.libs
includedir=@abs_top_srcdir@
Name: @PACKAGE_NAME@
Description: @PACKAGE_STRING@
Version: @PACKAGE_VERSION@
URL: @PACKAGE_URL@
Libs: -L${libdir} -l@PACKAGE@
Cflags: -I${includedir}

27
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -28,8 +28,8 @@ namespace Opm
{
rock_.init(deck, grid);
pvt_.init(deck, 200);
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, 200);
@ -50,30 +50,43 @@ namespace Opm
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
bool use_stone2 = (threephase_model == "stone2");
if (sat_samples > 1) {
if (use_stone2) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
THROW("Unknown threephase_model: " << threephase_model);
}
} else {
if (use_stone2) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, sat_samples);
} else {
THROW("Unknown threephase_model: " << threephase_model);
}
}

440
opm/core/fluid/SatFuncGwseg.cpp Normal file
View File

@ -0,0 +1,440 @@
/*
Copyright 2012 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 <opm/core/fluid/SatFuncGwseg.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
void SatFuncGwsegUniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegUniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(ssw);
const double dkrgg = krg_.derivative(ssg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(ssg);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww;
dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg;
dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegUniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncGwsegNonuniform ======
void SatFuncGwsegNonuniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int /*samples*/)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
krw_ = NonuniformTableLinear<double>(sw, krw);
krow_ = NonuniformTableLinear<double>(sw, krow);
pcow_ = NonuniformTableLinear<double>(sw, pcow);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
krg_ = NonuniformTableLinear<double>(sg, krg);
krog_ = NonuniformTableLinear<double>(sg, krog);
pcog_ = NonuniformTableLinear<double>(sg, pcog);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegNonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(ssw);
const double dkrgg = krg_.derivative(ssg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(ssg);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww;
dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg;
dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
ASSERT(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegNonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm

80
opm/core/fluid/SatFuncGwseg.hpp Normal file
View File

@ -0,0 +1,80 @@
/*
Copyright 2012 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 SATFUNCGWSEG_HPP
#define SATFUNCGWSEG_HPP
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <vector>
namespace Opm
{
class SatFuncGwsegUniform : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
class SatFuncGwsegNonuniform : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
} // namespace Opm
#endif // SATFUNCGWSEG_HPP

8
opm/core/fluid/SatFuncStone2.cpp
View File

@ -33,9 +33,9 @@ namespace Opm
void SatFuncStone2Uniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
@ -221,7 +221,7 @@ namespace Opm
// ====== Methods for SatFuncSimpleNonuniform ======
// ====== Methods for SatFuncStone2Nonuniform ======

1
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@ -26,6 +26,7 @@
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/SatFuncGwseg.hpp>
#include <vector>
struct UnstructuredGrid;

4
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -26,6 +26,10 @@
#include <opm/core/utility/have_boost_redef.hpp>
// Silence compatibility warning from DUNE headers since we don't use
// the deprecated member anyway (in this compilation unit)
#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1
// TODO: clean up includes.
#include <dune/common/deprecated.hh>
#include <dune/istl/bvector.hh>

25
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -123,7 +123,7 @@ namespace Opm
WellState& well_state)
{
const int nc = grid_.number_of_cells;
const int nw = wells_->number_of_wells;
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
// Set up dynamic data.
computePerSolveDynamicData(dt, state, well_state);
@ -454,8 +454,8 @@ namespace Opm
// std::vector<double> wellperf_A_;
// std::vector<double> wellperf_phasemob_;
const int np = props_.numPhases();
const int nw = wells_->number_of_wells;
const int nperf = wells_->well_connpos[nw];
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
wellperf_A_.resize(nperf*np*np);
wellperf_phasemob_.resize(nperf*np);
// The A matrix is set equal to the perforation grid cells'
@ -512,9 +512,9 @@ namespace Opm
const double* z = &state.surfacevol()[0];
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = &wellperf_gpot_[0];
completion_data.A = &wellperf_A_[0];
completion_data.phasemob = &wellperf_phasemob_[0];
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@ -599,9 +599,9 @@ namespace Opm
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = const_cast<double*>(&wellperf_gpot_[0]);
completion_data.A = const_cast<double*>(&wellperf_A_[0]);
completion_data.phasemob = const_cast<double*>(&wellperf_phasemob_[0]);
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@ -609,6 +609,9 @@ namespace Opm
forces.wells = &wells_tmp;
forces.src = NULL;
double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0;
double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;
cfs_tpfa_res_flux(gg,
&forces,
props_.numPhases(),
@ -617,9 +620,9 @@ namespace Opm
&face_phasemob_[0],
&face_gravcap_[0],
&state.pressure()[0],
&well_state.bhp()[0],
wpress,
&state.faceflux()[0],
&well_state.perfRates()[0]);
wflux);
cfs_tpfa_res_fpress(gg,
props_.numPhases(),
&htrans_[0],

84
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G ,
nnu = G->number_of_cells;
nwperf = 0;
if (wells != NULL) { assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nnu += wells->W->number_of_wells;
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, w, i, nc, nnu;
int wells_present;
size_t nnz;
struct CSRMatrix *A;
nc = nnu = G->number_of_cells;
if (wells != NULL) {
assert (wells->W != NULL);
wells_present = (wells != NULL) && (wells->W != NULL);
if (wells_present) {
nnu += wells->W->number_of_wells;
}
@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Well <-> cell connections */
struct Wells *W = wells->W;
@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Fill well <-> cell connections */
struct Wells *W = wells->W;
@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h,
double *res, double *w2c, double *w2w)
/* ---------------------------------------------------------------------- */
{
int p;
double fwi;
const double *dpflux_w;
/* Sum reservoir phase flux derivatives set by
* compute_darcyflux_and_deriv(). */
dpflux_w = h->pimpl->flux_work + (1 * np);
for (p = 0, fwi = 0.0; p < np; p++) {
fwi += dpflux_w[ p ];
}
*res = 0.0;
*w2c = 0.0;
*w2w = fwi;
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
@ -815,23 +840,25 @@ assemble_completion_to_well(int w, int c, int nc, int np,
ctrl = W->ctrls[ w ];
if (ctrl->current < 0) {
assert (0); /* Shut wells currently not supported */
/* Interpreting a negative current control index to mean a shut well */
welleq_coeff_shut(np, h, &res, &w2c, &w2w);
}
else {
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case SURFACE_RATE:
assert (0); /* Surface rate currently not supported */
break;
case SURFACE_RATE:
assert (0); /* Surface rate currently not supported */
break;
}
}
/* Assemble completion contributions */
@ -854,7 +881,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
struct cfs_tpfa_res_data *h )
{
int w, i, c, np, np2, nc;
int is_neumann;
int is_neumann, is_open;
double pw, dp;
double *WI, *gpot, *pmobp;
const double *Ac, *dAc;
@ -876,6 +903,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[ w ];
is_open = W->ctrls[w]->current >= 0;
for (; i < W->well_connpos[w + 1];
i++, gpot += np, pmobp += np) {
@ -888,7 +916,9 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
assemble_completion_to_cell(c, nc + w, np, dt, h);
if (is_open) {
assemble_completion_to_cell(c, nc + w, np, dt, h);
}
/* Prepare for RESV controls */
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
@ -1127,8 +1157,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
nf = G->number_of_faces;
nwperf = 0;
if (wells != NULL) {
assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@ -1194,7 +1223,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
assemble_cell_contrib(G, c, h);
}
if ((forces != NULL) && (forces->wells != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL)) {
compute_well_compflux_and_deriv(forces->wells, cq->nphases,
cpress, wpress, h->pimpl);
@ -1297,8 +1328,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
{
compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
if ((forces != NULL) && (forces->wells != NULL) &&
(wpress != NULL) && (wflux != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) {
compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
}

1
opm/core/utility/Units.hpp
View File

@ -101,6 +101,7 @@ namespace Opm
/// @{
const double gallon = 231 * cubic(inch);
const double stb = 42 * gallon;
const double liter = 1 * cubic(prefix::deci*meter);
/// @}
/// \name Mass

5
opm/core/utility/miscUtilities.cpp
View File

@ -374,7 +374,10 @@ namespace Opm
if (perf_rate > 0.0) {
// perf_rate is a total inflow rate, we want a water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl;
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);

51
opm/core/wells/WellsManager.cpp
View File

@ -716,32 +716,31 @@ namespace Opm
#endif
if (deck.hasField("WELOPEN")) {
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
int index = it->second;
if (line.openshutflag_ == "SHUT") {
// We currently don't care if the well is open or not.
/// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }?
//ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
//ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
// We revert back to it's original control.
// Note that this is OK as ~~ = id.
w_->ctrls[index]->current = ~w_->ctrls[index]->current;
}
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
const int index = it->second;
if (line.openshutflag_ == "SHUT") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl >= 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl < 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
}
}
// Build the well_collection_ well group hierarchy.