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 branch 'PR/performance-avoid-zero-matrix-entries' of ssh://github.com/dr-robertk/opm-autodiff into PR/performance-avoid-zero-matrix-entries

Browse Source 
Conflicts:
	opm/autodiff/ConservativeSparseSparseProduct.h
This commit is contained in:
Robert K 2014-12-03 15:50:41 +01:00
commit 5583585063
20 changed files with 1096 additions and 114 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -29,3 +29,6 @@ test_vec
# Mac OS X debug info
*.dSYM
# emacs directory setting:
.dir-locals.el

2
CMakeLists_files.cmake
View File

@ -77,6 +77,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/sim_2p_comp_ad.cpp
examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp
examples/opm_init_check.cpp
)
# programs listed here will not only be compiled, but also marked for
@ -84,6 +85,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
list (APPEND PROGRAM_SOURCE_FILES
examples/sim_2p_incomp_ad.cpp
examples/sim_fibo_ad.cpp
examples/opm_init_check.cpp
)
# originally generated with the command:

363
examples/opm_init_check.cpp Normal file
View File

@ -0,0 +1,363 @@
/*
Copyright 2014 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 <config.h>
#include <iostream>
#include <fstream>
#include <memory>
#include <vector>
#include <string>
#include <ert/ecl/ecl_file.h>
#include <ert/ecl/ecl_kw.h>
#include <ert/ecl/ecl_grid.h>
#include <ert/ecl/ecl_nnc_export.h>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
using namespace Opm;
class CellTrans {
public:
CellTrans(const std::tuple<int,int,int>& ijk) :
m_ijk(ijk)
{
}
void update(const std::tuple<int,int,int>& ijk , double trans) {
auto iter = m_trans.find( ijk );
if (iter == m_trans.end())
m_trans.insert( std::pair<std::tuple<int,int,int> , double>(ijk , trans));
else
iter->second *= trans;
}
int numConnections() const {
return m_trans.size();
}
static void jsonDumpIJK(std::ostream& os , const std::tuple<int,int,int>& ijk ) {
os << "[" << std::get<0>(ijk) << "," << std::get<1>(ijk) << "," << std::get<2>(ijk) << "]";
}
void jsonDump(std::ostream& os , bool first) {
if (!first)
os <<",";
os << "[";
jsonDumpIJK(os , m_ijk );
os << ",[";
{
size_t count = 0;
for (auto transIter = m_trans.begin(); transIter != m_trans.end(); ++transIter) {
std::tuple<int,int,int> ijk = transIter->first;
double t = transIter->second;
os << "[";
jsonDumpIJK( os , ijk );
os << "," << t << "]";
count++;
if (count < m_trans.size())
os << ",";
else
os << "]";
}
}
os << "]" << std::endl;
}
std::tuple<int,int,int> m_ijk;
std::map<std::tuple<int,int,int> , double> m_trans;
};
class TransGraph {
public:
TransGraph(int nx , int ny , int nz) :
m_nx(nx),
m_ny(ny),
m_nz(nz)
{
m_transVector.resize( nx*ny*nz );
}
void update(int global_index1 , int global_index2 , double trans) {
int size = m_nx * m_ny * m_nz;
if ((global_index1 >= 0) &&
(global_index2 >= 0) &&
(global_index1 < size) &&
(global_index2 < size)) {
size_t g1 = std::min( global_index1 , global_index2 );
size_t g2 = std::max( global_index1 , global_index2 );
std::shared_ptr<CellTrans> cellTrans = m_transVector[g1];
if (!cellTrans) {
cellTrans = std::make_shared<CellTrans>( getIJK(g1) );
m_transVector[g1] = cellTrans;
}
cellTrans->update( getIJK(g2) , trans );
}
}
int activeCells() const {
int count = 0;
for (size_t g= 0; g < m_transVector.size(); g++) {
std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
if (cellTrans)
count++;
}
return count;
}
int activeConnections() const {
int count = 0;
for (size_t g= 0; g < m_transVector.size(); g++) {
std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
if (cellTrans)
count += cellTrans->numConnections();
}
return count;
}
std::tuple<int,int,int> getIJK(int g) const {
int k = g / (m_nx * m_ny);
int j = (g - k*m_nx*m_ny) / m_nx;
int i = g - k*m_nx*m_ny - j*m_nx;
return std::tuple<int,int,int>(i,j,k);
}
void jsonDump(const std::string& outputFile) {
std::ofstream os;
bool first = true;
os.open(outputFile.c_str());
os << "{ \"dims\" : [" << m_nx << "," << m_ny << "," << m_nz << "]" << " , \"graph\":[";
for (size_t g= 0; g < m_transVector.size(); g++) {
std::shared_ptr<CellTrans> cellTrans = m_transVector[g];
if (cellTrans) {
cellTrans->jsonDump(os , first);
first = false;
}
}
os << "]}" << std::endl;
os.close();
}
int m_nx;
int m_ny;
int m_nz;
std::vector<std::shared_ptr<CellTrans> > m_transVector;
};
/*****************************************************************/
void initOPMTrans(TransGraph& opmTrans , DeckConstPtr deck , std::shared_ptr<const EclipseState> eclipseState) {
std::shared_ptr<GridManager> grid = std::make_shared<GridManager>( eclipseState->getEclipseGrid(), eclipseState->getDoubleGridProperty("PORV")->getData());
const struct UnstructuredGrid * cGrid = grid->c_grid();
std::shared_ptr<BlackoilPropsAdInterface> props;
props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
DerivedGeology geology(*grid->c_grid() , *props, eclipseState);
const double * opm_trans_data = geology.transmissibility().data();
double SIconversion = Opm::unit::cubic(Opm::unit::meter) * Opm::unit::day * Opm::unit::barsa / (Opm::prefix::centi * Opm::unit::Poise);
{
for (int face_index = 0; face_index < cGrid->number_of_faces; face_index++ ) {
int global_index1 = cGrid->global_cell[ cGrid->face_cells[2*face_index] ];
int global_index2 = cGrid->global_cell[ cGrid->face_cells[2*face_index + 1] ];
opmTrans.update( global_index1 , global_index2 , opm_trans_data[ face_index ] * SIconversion );
}
}
}
void initEclipseTrans(TransGraph& eclipseTrans , const ecl_grid_type * ecl_grid , const ecl_file_type * ecl_init) {
int nx = ecl_grid_get_nx( ecl_grid );
int ny = ecl_grid_get_ny( ecl_grid );
int nz = ecl_grid_get_nz( ecl_grid );
if (ecl_file_has_kw( ecl_init , "TRANX")) {
ecl_kw_type * tranx_kw = ecl_file_iget_named_kw( ecl_init , "TRANX" , 0 );
ecl_kw_type * trany_kw = ecl_file_iget_named_kw( ecl_init , "TRANY" , 0 );
ecl_kw_type * tranz_kw = ecl_file_iget_named_kw( ecl_init , "TRANZ" , 0 );
for (int k=0; k < nz; k++) {
for (int j= 0; j < ny; j++) {
for (int i=0; i < nx; i++) {
if (ecl_grid_cell_active3( ecl_grid , i , j , k )) {
size_t g1 = ecl_grid_get_global_index3( ecl_grid , i , j , k );
int a = ecl_grid_get_active_index1( ecl_grid , g1 );
if (a >= 0) {
if (i < (nx - 1) && ecl_grid_cell_active3( ecl_grid , i + 1 , j , k)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i + 1, j , k );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( tranx_kw , a ));
}
if (j < (ny - 1) && ecl_grid_cell_active3( ecl_grid , i , j + 1, k)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j + 1, k );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( trany_kw , a ));
}
if (k < (nz - 1) && ecl_grid_cell_active3( ecl_grid , i , j , k + 1)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j , k + 1 );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( tranz_kw , a ));
}
}
}
}
}
}
} else
std::cerr << "Init file does not have TRAN[XYZ] keywords" << std::endl;
if (ecl_file_has_kw( ecl_init , "TRANX-")) {
ecl_kw_type * tranxm_kw = ecl_file_iget_named_kw( ecl_init , "TRANX-" , 0 );
ecl_kw_type * tranym_kw = ecl_file_iget_named_kw( ecl_init , "TRANY-" , 0 );
ecl_kw_type * tranzm_kw = ecl_file_iget_named_kw( ecl_init , "TRANZ-" , 0 );
for (int k=0; k < nz; k++) {
for (int j= 0; j < ny; j++) {
for (int i=0; i < nx; i++) {
if (ecl_grid_cell_active3( ecl_grid , i , j , k )) {
size_t g1 = ecl_grid_get_global_index3( ecl_grid , i , j , k );
int a = ecl_grid_get_active_index1( ecl_grid , g1 );
if (a >= 0) {
if (i > 0 && ecl_grid_cell_active3( ecl_grid , i - 1 , j , k)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i - 1, j , k );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( tranxm_kw , a ));
}
if (j > 0 && ecl_grid_cell_active3( ecl_grid , i , j - 1, k)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j - 1, k );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( tranym_kw , a ));
}
if (k > 0 && ecl_grid_cell_active3( ecl_grid , i , j , k - 1)) {
size_t g2 = ecl_grid_get_global_index3( ecl_grid , i , j , k - 1 );
eclipseTrans.update( g1 , g2 , ecl_kw_iget_float( tranzm_kw , a ));
}
}
}
}
}
}
}
// NNC
{
size_t num_nnc = static_cast<size_t>( ecl_nnc_export_get_size( ecl_grid ));
std::vector<ecl_nnc_type> nnc(num_nnc);
ecl_nnc_export( ecl_grid , ecl_init , nnc.data());
for (auto nnc_iter = nnc.begin(); nnc_iter != nnc.end(); ++nnc_iter)
eclipseTrans.update( nnc_iter->global_index1 , nnc_iter->global_index2 , nnc_iter->trans );
}
}
void dump_transGraph( DeckConstPtr deck , std::shared_ptr<const EclipseState> eclipseState , const ecl_grid_type * ecl_grid , const ecl_file_type * ecl_init , size_t verbosity) {
int nx = ecl_grid_get_nx( ecl_grid );
int ny = ecl_grid_get_ny( ecl_grid );
int nz = ecl_grid_get_nz( ecl_grid );
TransGraph opmTrans(nx , ny , nz );
TransGraph eclipseTrans( nx , ny , nz);
initOPMTrans( opmTrans , deck , eclipseState );
initEclipseTrans( eclipseTrans , ecl_grid , ecl_init );
opmTrans.jsonDump("opm_trans.json");
eclipseTrans.jsonDump("eclipse_trans.json");
}
int main(int argc, char** argv) {
if (argc < 4) {
std::cerr << "The opm_init_check program needs three arguments:" << std::endl << std::endl;;
std::cerr << " ECLIPSE.DATA ECLIPSE.INIT ECLIPSE.EGRID" << std::endl << std::endl;
std::cerr << "Where the ECLIPSE.INIT and ECLIPSE.EGRID are existing binary files";
exit(1);
}
std::string input_file = argv[1];
std::string init_file = argv[2];
std::string grid_file = argv[3];
ParserPtr parser(new Parser());
std::cout << "Parsing input file ............: " << input_file << std::endl;
DeckConstPtr deck = parser->parseFile(input_file, false);
std::shared_ptr<EclipseState> state = std::make_shared<EclipseState>( deck );
std::cout << "Loading eclipse INIT file .....: " << init_file << std::endl;
ecl_file_type * ecl_init = ecl_file_open( init_file.c_str() , 0 );
std::cout << "Loading eclipse EGRID file ....: " << grid_file << std::endl;
ecl_grid_type * ecl_grid = ecl_grid_alloc( grid_file.c_str() );
dump_transGraph( deck , state , ecl_grid , ecl_init , 3);
ecl_file_close( ecl_init );
ecl_grid_free( ecl_grid );
return 0;
}

136
examples/python/example.py Executable file
View File

@ -0,0 +1,136 @@
#!/usr/bin/env python
import sys
from trans_graph import TransGraph
# This file just contains some example functions for how one can poke
# around in the TransGraph datastructure.
def direction_count(tg):
dir_count = {"X" : 0 ,
"Y" : 0 ,
"Z" : 0 ,
"NNC" : 0 }
for cell in tg:
if cell:
for conn in cell:
dir_count[ conn.dir ] += 1
dir_count["Total"] = dir_count["X"] + dir_count["Y"] + dir_count["Z"] + dir_count["NNC"]
return dir_count
def print_cell( prefix , cell ):
print "%s: Cell: (%d,%d,%d) " % (prefix , cell.i , cell.j , cell.k)
for conn in cell:
print " Connection => (%3d,%3d,%3d) Transmissibility: %g Direction: %s" % (conn.i , conn.j , conn.k , conn.T , conn.dir)
print " cell[\"X\"] => %s" % cell["X"]
print " cell[\"Y\"] => %s" % cell["Y"]
print " cell[\"Z\"] => %s" % cell["Z"]
def connection_to_count(tg , i,j,k):
count = 0
for cell in tg:
if cell.connectsWith(i,j,k):
count += 1
return count
#-----------------------------------------------------------------
def direction_example( opm_tg , ecl_tg ):
opm_count = direction_count( opm_tg )
ecl_count = direction_count( ecl_tg )
print "OPM: %s" % opm_count
print "ECL: %s" % ecl_count
print
def cell_example(opm_tg , ecl_tg ):
opm_cell = opm_tg[21,27,10]
ecl_cell = ecl_tg[21,27,10]
print_cell( "OPM: " , opm_cell )
print_cell( "ECL: " , ecl_cell )
def count_example(opm_tg , ecl_tg ):
i = 5
j = 10
k = 7
print "Opm connections to (%d,%d,%d): %d" % (i,j,k , connection_to_count(opm_tg , i,j,k))
print "Ecl connections to (%d,%d,%d): %d" % (i,j,k , connection_to_count(ecl_tg , i,j,k))
print
def xtrace_example( opm_tg , ecl_tg , j , k):
opm_trace = opm_tg.getXTrace( j , k)
ecl_trace = ecl_tg.getXTrace( j , k)
print "OPM: %s" % opm_trace
print "ECL: %s" % ecl_trace
print
def ytrace_example( opm_tg , ecl_tg , i , k):
opm_trace = opm_tg.getYTrace( i , k )
ecl_trace = ecl_tg.getYTrace( i , k )
print "OPM: %s" % opm_trace
print "ECL: %s" % ecl_trace
print
def ztrace_example( opm_tg , ecl_tg , i , j):
opm_trace = opm_tg.getZTrace( i , j )
ecl_trace = ecl_tg.getZTrace( i , j )
print "OPM: %s" % opm_trace
print "ECL: %s" % ecl_trace
print
#-----------------------------------------------------------------
if len(sys.argv) < 3:
sys.exit("example.py opm_trans.json eclipse_trans.json")
opm_tg = TransGraph.load(sys.argv[1])
ecl_tg = TransGraph.load(sys.argv[2])
direction_example( opm_tg , ecl_tg )
cell_example( opm_tg , ecl_tg )
count_example( opm_tg , ecl_tg )
xtrace_example( opm_tg , ecl_tg , 20 ,20)
ytrace_example( opm_tg , ecl_tg , 10 ,10)
ztrace_example( opm_tg , ecl_tg , 10 ,10)
ztrace_example( opm_tg , ecl_tg , 10 ,20)
ztrace_example( opm_tg , ecl_tg , 20 ,10)
ztrace_example( opm_tg , ecl_tg , 20 ,20)
ztrace_example( opm_tg , ecl_tg , 30 ,70)
ztrace_example( opm_tg , ecl_tg , 30 ,60)

198
examples/python/trans_graph.py Normal file
View File

@ -0,0 +1,198 @@
# The opm_init_check cpp program will produce two json files which
# should describe the transmissibility graph. The structure of the
# main chunk of data is a list:
#
#
# ((i1,j1,k1) , (((i2,j2,k2), T12) , ((i3,j3,k3) , T13))),
# ((iA,jA,kA) , (((iB,jB,kB), TAB) , ((iC,jC,kC) , TAC))),
# ....
#
#
# Cell(i1,j1,k1) is connected to cells (i2,j2,k2) and (i3,j3,k3)
# respectively, with transmissibility T12 and T13
# respectively. Furthermore cell (iA,iB,iC) is connected to cells
# (iB,jB,kB) and (iC,jC,kC) with transmissibilty TAB and TAC
# respectively.
import json
class Connection(object):
"""
The connection class holds connection information for one cell;
including the i,j,k of the cell and the Transmissibility of connection.
"""
def __init__(self , i1, j1 , k1 , i2 , j2 , k2 , T):
self.i = i2
self.j = j2
self.k = k2
self.T = T
dx = abs(i1 - i2)
dy = abs(j1 - j2)
dz = abs(k1 - k2)
if dx == 1 and dy == 0 and dz == 0:
self.dir = "X"
elif dx == 0 and dy == 1 and dz == 0:
self.dir = "Y"
elif dx == 0 and dy == 0 and dz == 1:
self.dir = "Z"
else:
self.dir = "NNC"
def __str__(self):
return "<%d,%d,%d>(T:%g)" % (self.i , self.j , self.k , self.T)
class CellConnections(object):
def __init__(self , i,j,k):
self.i = i
self.j = j
self.k = k
self.connection_list = []
self.connection_map = {}
def __getitem__(self , index):
if isinstance(index,int):
return self.connection_list[index]
else:
return self.connection_map[index]
def __len__(self):
return len(self.connection_list)
def has_key(self , dir_key):
return self.connection_map.has_key(dir_key)
def addConnections(self , connections):
for ijk,T in connections:
new_conn = Connection( self.i , self.j , self.k , ijk[0] , ijk[1] , ijk[2] , T)
self.connection_list.append( new_conn )
if new_conn.dir == "NNC":
if not self.connection_map.has_key("NNC"):
self.connection_map["NNC"] = []
self.connection_map["NNC"].append( new_conn )
else:
self.connection_map[new_conn.dir] = new_conn
def __nonzero__(self):
if len(self.connection_list) > 0:
return True
else:
return False
def connectsWith(self, i , j , k):
for conn in self.connection_list:
if conn.i == i and conn.j == j and conn.k == k:
return True
else:
return False
def __str__(self):
return "<%d,%d,%d> : %s" % (self.i , self.j , self.k , [ conn.__str__() for conn in self.connection_list ])
class TransGraph(object):
def __init__(self , nx , ny , nz):
self.nx = nx
self.ny = ny
self.nz = nz
self.cell_connections = [ None ] * nx * ny * nz
def __getitem__(self, index):
if isinstance(index , tuple):
g = index[0] + index[1] * self.nx + index[2]*self.nx*self.ny
else:
g = index
connections = self.cell_connections[g]
if connections is None:
k = g / (self.nx * self.ny)
j = (g - k * self.nx * self.ny) / self.nx
i = g - k * self.nx * self.ny - j * self.nx
self.cell_connections[g] = CellConnections( i,j,k )
return self.cell_connections[g]
def addCell(self , ijk , new_connections):
g = ijk[0] + ijk[1] * self.nx + ijk[2]*self.nx*self.ny
connections = self[g]
connections.addConnections( new_connections )
def getZTrace(self , i , j):
trace = []
for k in range(self.nz):
cell = self[i,j,k]
if cell.has_key("Z"):
trace.append( cell["Z"].T )
else:
trace.append( 0 )
return trace
def getXTrace(self , j , k):
trace = []
for i in range(self.nx):
cell = self[i,j,k]
if cell.has_key("X"):
trace.append( cell["X"].T )
else:
trace.append( 0 )
return trace
def getYTrace(self , i , k):
trace = []
for j in range(self.ny):
cell = self[i,j,k]
if cell.has_key("Y"):
trace.append( cell["Y"].T )
else:
trace.append( 0 )
return trace
@staticmethod
def load(json_file):
with open(json_file) as fileH:
data = json.load( fileH )
dims = data["dims"]
graph = data["graph"]
trans_graph = TransGraph( dims[0] , dims[1] , dims[2])
for cell,connections in graph:
trans_graph.addCell( cell , connections )
return trans_graph

55
opm/autodiff/AutoDiffHelpers.hpp
View File

@ -189,11 +189,11 @@ namespace {
template <typename Scalar, class IntVec>
Eigen::SparseMatrix<Scalar>
typename AutoDiffBlock<Scalar>::M
constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
{
const int subset_size = indices.size();
Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
typename AutoDiffBlock<Scalar>::M mat(full_size, subset_size);
mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
for (int i = 0; i < subset_size; ++i) {
mat.insert(indices[i], i) = 1;
@ -204,18 +204,6 @@ namespace {
} // anon namespace
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices)
{
Eigen::SparseMatrix<Scalar> sub
= constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose();
return sub * x;
}
/// Returns x(indices).
template <typename Scalar, class IntVec>
@ -223,9 +211,25 @@ Eigen::Array<Scalar, Eigen::Dynamic, 1>
subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices)
{
return (constructSupersetSparseMatrix<Scalar>(x.size(), indices).transpose() * x.matrix()).array();
typedef typename Eigen::Array<Scalar, Eigen::Dynamic, 1>::Index Index;
const Index size = indices.size();
Eigen::Array<Scalar, Eigen::Dynamic, 1> ret( size );
for( Index i=0; i<size; ++i )
ret[ i ] = x[ indices[ i ] ];
return std::move(ret);
}
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices)
{
const typename AutoDiffBlock<Scalar>::M sub
= constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose();
return sub * x;
}
/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
@ -357,9 +361,10 @@ spdiag(const AutoDiffBlock<double>::V& d)
/// Returns the input expression, but with all Jacobians collapsed to one.
template <class Matrix>
inline
AutoDiffBlock<double>
collapseJacs(const AutoDiffBlock<double>& x)
void
collapseJacs(const AutoDiffBlock<double>& x, Matrix& jacobian)
{
typedef AutoDiffBlock<double> ADB;
const int nb = x.numBlocks();
@ -383,9 +388,21 @@ collapseJacs(const AutoDiffBlock<double>& x)
block_col_start += jac.cols();
}
// Build final jacobian.
jacobian = Matrix(x.size(), block_col_start);
jacobian.setFromTriplets(t.begin(), t.end());
}
/// Returns the input expression, but with all Jacobians collapsed to one.
inline
AutoDiffBlock<double>
collapseJacs(const AutoDiffBlock<double>& x)
{
typedef AutoDiffBlock<double> ADB;
// Build final jacobian.
std::vector<ADB::M> jacs(1);
jacs[0].resize(x.size(), block_col_start);
jacs[0].setFromTriplets(t.begin(), t.end());
collapseJacs( x, jacs[ 0 ] );
return ADB::function(x.value(), jacs);
}

66
opm/autodiff/BlackoilPropsAd.cpp
View File

@ -103,9 +103,11 @@ namespace Opm
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
@ -116,17 +118,19 @@ namespace Opm
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);
props_.viscosity(n, pw.data(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -145,15 +149,17 @@ namespace Opm
z.col(pu_.phase_pos[Gas]) = rs;
}
Block mu(n, np);
props_.viscosity(n, po.data(), z.data(), cells.data(), mu.data(), 0);
props_.viscosity(n, po.data(), T.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] T Array of n temperature 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 V& T,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
@ -164,17 +170,19 @@ namespace Opm
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);
props_.viscosity(n, pg.data(), T.data(), z.data(), cells.data(), mu.data(), 0);
return mu.col(pu_.phase_pos[Gas]);
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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::muGas(const V& pg,
const V& T,
const V& rv,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -193,19 +201,21 @@ namespace Opm
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
Block mu(n, np);
props_.viscosity(n, pg.data(), z.data(), cells.data(), mu.data(), 0);
props_.viscosity(n, pg.data(), T.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] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
#if 1
return ADB::constant(muWat(pw.value(), cells), pw.blockPattern());
return ADB::constant(muWat(pw.value(), T.value(), cells), pw.blockPattern());
#else
if (!pu_.phase_used[Water]) {
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
@ -216,7 +226,7 @@ namespace Opm
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());
props_.viscosity(n, pw.value().data(), T.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);
@ -229,17 +239,19 @@ namespace Opm
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
{
#if 1
return ADB::constant(muOil(po.value(), rs.value(), cond, cells), po.blockPattern());
return ADB::constant(muOil(po.value(), T.value(), rs.value(), cond, cells), po.blockPattern());
#else
if (!pu_.phase_used[Oil]) {
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
@ -272,13 +284,15 @@ namespace Opm
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
#if 1
return ADB::constant(muGas(pg.value(), cells), pg.blockPattern());
return ADB::constant(muGas(pg.value(), T.value(), cells), pg.blockPattern());
#else
if (!pu_.phase_used[Gas]) {
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
@ -301,17 +315,19 @@ namespace Opm
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
{
#if 1
return ADB::constant(muGas(pg.value(), rv.value(),cond,cells), pg.blockPattern());
return ADB::constant(muGas(pg.value(), T.value(), rv.value(),cond,cells), pg.blockPattern());
#else
if (!pu_.phase_used[Gas]) {
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
@ -359,9 +375,11 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
@ -372,18 +390,20 @@ namespace Opm
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);
props_.matrix(n, pw.data(), T.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] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -402,16 +422,18 @@ namespace Opm
z.col(pu_.phase_pos[Gas]) = rs;
}
Block matrix(n, np*np);
props_.matrix(n, po.data(), z.data(), cells.data(), matrix.data(), 0);
props_.matrix(n, po.data(), T.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] T Array of n temperature 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 V& T,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
@ -422,18 +444,20 @@ namespace Opm
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);
props_.matrix(n, pg.data(), pg.data(), z.data(), cells.data(), matrix.data(), 0);
const int gi = pu_.phase_pos[Gas];
return matrix.col(gi*np + gi);
}
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -452,16 +476,18 @@ namespace Opm
z.col(pu_.phase_pos[Gas]) = V::Ones(n, 1);
}
Block matrix(n, np*np);
props_.matrix(n, pg.data(), z.data(), cells.data(), matrix.data(), 0);
props_.matrix(n, pg.data(), T.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] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
if (!pu_.phase_used[Water]) {
@ -473,7 +499,7 @@ namespace Opm
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());
props_.matrix(n, pw.value().data(), T.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));
@ -487,11 +513,13 @@ namespace Opm
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -511,7 +539,7 @@ namespace Opm
}
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, po.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
props_.matrix(n, po.value().data(), T.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));
@ -528,9 +556,11 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const
{
if (!pu_.phase_used[Gas]) {
@ -542,7 +572,7 @@ namespace Opm
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());
props_.matrix(n, pg.value().data(), T.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));
@ -556,11 +586,13 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& /*cond*/,
const Cells& cells) const
@ -580,7 +612,7 @@ namespace Opm
}
Block matrix(n, np*np);
Block dmatrix(n, np*np);
props_.matrix(n, pg.value().data(), z.data(), cells.data(), matrix.data(), dmatrix.data());
props_.matrix(n, pg.value().data(), T.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));

32
opm/autodiff/BlackoilPropsAd.hpp
View File

@ -101,72 +101,88 @@ namespace Opm
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
@ -175,73 +191,89 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;

64
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -269,9 +269,11 @@ namespace Opm
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAdFromDeck::muWat(const V& pw,
const V& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Water]) {
@ -284,18 +286,20 @@ namespace Opm
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.data(), rs,
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.data(), T.data(), rs,
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAdFromDeck::muOil(const V& po,
const V& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -309,16 +313,18 @@ namespace Opm
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.data(), rs.data(), &cond[0],
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.data(), T.data(), rs.data(), &cond[0],
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAdFromDeck::muGas(const V& pg,
const V& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Gas]) {
@ -331,16 +337,18 @@ namespace Opm
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), rs,
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), T.data(), rs,
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
V BlackoilPropsAdFromDeck::muGas(const V& pg,
const V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -354,16 +362,18 @@ namespace Opm
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), rv.data(),&cond[0],
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), T.data(), rv.data(),&cond[0],
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAdFromDeck::muWat(const ADB& pw,
const ADB& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Water]) {
@ -376,7 +386,7 @@ namespace Opm
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.value().data(), rs,
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.value().data(), T.value().data(), rs,
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
const int num_blocks = pw.numBlocks();
@ -389,11 +399,13 @@ namespace Opm
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAdFromDeck::muOil(const ADB& po,
const ADB& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -407,7 +419,7 @@ namespace Opm
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.value().data(), rs.value().data(),
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.value().data(), T.value().data(), rs.value().data(),
&cond[0], mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -422,9 +434,11 @@ namespace Opm
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAdFromDeck::muGas(const ADB& pg,
const ADB& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Gas]) {
@ -437,7 +451,7 @@ namespace Opm
V dmudr(n);
const double* rv = 0;
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), rv,
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv,
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -451,11 +465,13 @@ namespace Opm
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB BlackoilPropsAdFromDeck::muGas(const ADB& pg,
const ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -469,7 +485,7 @@ namespace Opm
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), rv.value().data(),&cond[0],
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv.value().data(),&cond[0],
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -502,9 +518,11 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 BlackoilPropsAdFromDeck::bWat(const V& pw,
const V& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Water]) {
@ -518,7 +536,7 @@ namespace Opm
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.data(), rs,
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.data(), T.data(), rs,
b.data(), dbdp.data(), dbdr.data());
return b;
@ -526,11 +544,13 @@ namespace Opm
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V BlackoilPropsAdFromDeck::bOil(const V& po,
const V& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -545,7 +565,7 @@ namespace Opm
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.data(), rs.data(), &cond[0],
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.data(), T.data(), rs.data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
return b;
@ -553,9 +573,11 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 BlackoilPropsAdFromDeck::bGas(const V& pg,
const V& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Gas]) {
@ -569,7 +591,7 @@ namespace Opm
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), rs,
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), T.data(), rs,
b.data(), dbdp.data(), dbdr.data());
return b;
@ -577,11 +599,13 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
V BlackoilPropsAdFromDeck::bGas(const V& pg,
const V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -596,7 +620,7 @@ namespace Opm
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), rv.data(), &cond[0],
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), T.data(), rv.data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
return b;
@ -604,9 +628,11 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 BlackoilPropsAdFromDeck::bWat(const ADB& pw,
const ADB& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Water]) {
@ -620,7 +646,7 @@ namespace Opm
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.value().data(), rs,
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.value().data(), T.value().data(), rs,
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -634,11 +660,13 @@ namespace Opm
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n taxonomies classifying fluid condition.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB BlackoilPropsAdFromDeck::bOil(const ADB& po,
const ADB& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -653,7 +681,7 @@ namespace Opm
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.value().data(), rs.value().data(),
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.value().data(), T.value().data(), rs.value().data(),
&cond[0], b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -668,9 +696,11 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 BlackoilPropsAdFromDeck::bGas(const ADB& pg,
const ADB& T,
const Cells& cells) const
{
if (!phase_usage_.phase_used[Gas]) {
@ -684,7 +714,7 @@ namespace Opm
V dbdr(n);
const double* rv = 0;
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), rv,
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv,
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -698,11 +728,13 @@ namespace Opm
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB BlackoilPropsAdFromDeck::bGas(const ADB& pg,
const ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const
@ -717,7 +749,7 @@ namespace Opm
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), rv.value().data(), &cond[0],
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv.value().data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);

32
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@ -126,70 +126,86 @@ namespace Opm
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& po,
const V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
@ -198,72 +214,88 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 V& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature 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 ADB& T,
const Cells& cells) const;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \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 ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const;

30
opm/autodiff/BlackoilPropsAdInterface.hpp
View File

@ -91,66 +91,80 @@ namespace Opm
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
V muWat(const V& pw,
const V& T,
const Cells& cells) const = 0;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
V muOil(const V& po,
const V& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
V muGas(const V& pg,
const V& T,
const Cells& cells) const = 0;
/// Water viscosity.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
ADB muWat(const ADB& pw,
const ADB& T,
const Cells& cells) const = 0;
/// Oil viscosity.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
ADB muOil(const ADB& po,
const ADB& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
ADB muGas(const ADB& pg,
const ADB& T,
const Cells& cells) const = 0;
/// Gas viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
virtual
ADB muGas(const ADB& pg,
const ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
@ -159,80 +173,96 @@ namespace Opm
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
V bWat(const V& pw,
const V& T,
const Cells& cells) const = 0;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
V bOil(const V& po,
const V& T,
const V& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
V bGas(const V& pg,
const V& T,
const Cells& cells) const = 0;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
V bGas(const V& pg,
const V& T,
const V& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
/// Water formation volume factor.
/// \param[in] pw Array of n water pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
ADB bWat(const ADB& pw,
const ADB& T,
const Cells& cells) const = 0;
/// Oil formation volume factor.
/// \param[in] po Array of n oil pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rs Array of n gas solution factor values.
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
ADB bOil(const ADB& po,
const ADB& T,
const ADB& rs,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
ADB bGas(const ADB& pg,
const ADB& T,
const Cells& cells) const = 0;
/// Gas formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] T Array of n temperature values.
/// \param[in] rv Array of n vapor oil/gas ratio
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
virtual
ADB bGas(const ADB& pg,
const ADB& T,
const ADB& rv,
const std::vector<PhasePresence>& cond,
const Cells& cells) const = 0;

34
opm/autodiff/ConservativeSparseSparseProduct.h
View File

@ -15,7 +15,10 @@
#include <algorithm>
#include <iterator>
#include <functional>
<<<<<<< HEAD
#include <limits>
=======
>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466
#include <vector>
#include <Eigen/Core>
@ -64,11 +67,14 @@ struct QuickSort< 0 >
template<typename Lhs, typename Rhs, typename ResultType>
static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
{
<<<<<<< HEAD
// if one of the matrices does not contain non zero elements
// the result will only contain an empty matrix
if( lhs.nonZeros() == 0 || rhs.nonZeros() == 0 )
return ;
=======
>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466
typedef typename remove_all<Lhs>::type::Scalar Scalar;
typedef typename remove_all<Lhs>::type::Index Index;
@ -91,10 +97,13 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
res.setZero();
res.reserve(Index(estimated_nnz_prod));
<<<<<<< HEAD
//const Scalar epsilon = std::numeric_limits< Scalar >::epsilon();
const Scalar epsilon = 1e-15 ;
=======
>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466
// we compute each column of the result, one after the other
for (Index j=0; j<cols; ++j)
{
@ -102,6 +111,7 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
{
const Scalar y = rhsIt.value();
<<<<<<< HEAD
for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
{
const Index i = lhsIt.index();
@ -117,6 +127,26 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
}
else
values[i] += val;
=======
if( y != 0.0 )
{
for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
{
const Index i = lhsIt.index();
const Scalar x = lhsIt.value();
if( x != 0.0 )
{
if(!mask[i])
{
mask[i] = true;
values[i] = x * y;
indices[nnz] = i;
++nnz;
}
else
values[i] += x * y;
}
>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466
}
}
}
@ -179,6 +209,10 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
res.finalize();
}
<<<<<<< HEAD
=======
>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466
} // end namespace internal
namespace internal {

4
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@ -142,6 +142,7 @@ namespace Opm {
struct SolutionState {
SolutionState(const int np);
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
@ -259,6 +260,7 @@ namespace Opm {
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -267,6 +269,7 @@ namespace Opm {
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -275,6 +278,7 @@ namespace Opm {
ADB
fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,

37
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -357,6 +357,7 @@ namespace {
template<class T>
FullyImplicitBlackoilSolver<T>::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, temperature( ADB::null())
, saturation(np, ADB::null())
, rs ( ADB::null())
, rv ( ADB::null())
@ -413,6 +414,7 @@ namespace {
// automatically consistent with variableState() (and doing
// things automatically is all the rage in this module ;)
state.pressure = ADB::constant(state.pressure.value());
state.temperature = ADB::constant(state.temperature.value());
state.rs = ADB::constant(state.rs.value());
state.rv = ADB::constant(state.rv.value());
for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx)
@ -512,6 +514,10 @@ namespace {
int nextvar = 0;
state.pressure = vars[ nextvar++ ];
// temperature
const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
state.temperature = ADB::constant(temp);
// Saturations
const std::vector<int>& bpat = vars[0].blockPattern();
{
@ -573,6 +579,7 @@ namespace {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const ADB& rs = state.rs;
const ADB& rv = state.rv;
@ -586,7 +593,7 @@ namespace {
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
rq_[pos].b = fluidReciprocFVF(phase, pressures[phase], rs, rv, cond, cells_);
rq_[pos].b = fluidReciprocFVF(phase, pressures[phase], temp, rs, rv, cond, cells_);
rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
// DUMP(rq_[pos].b);
// DUMP(rq_[pos].accum[aix]);
@ -620,6 +627,7 @@ namespace {
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
// Compute b, rsmax, rvmax values for perforations.
const ADB perf_press = subset(state.pressure, well_cells);
const ADB perf_temp = subset(state.temperature, well_cells);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
@ -630,21 +638,21 @@ namespace {
std::vector<double> rssat_perf(nperf, 0.0);
std::vector<double> rvsat_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const ADB bw = fluid_.bWat(perf_press, well_cells);
const ADB bw = fluid_.bWat(perf_press, perf_temp, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
}
assert(active_[Oil]);
const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_temp, perf_rs, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
rssat_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_temp, perf_rv, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
@ -1611,11 +1619,11 @@ namespace {
const std::vector<PhasePresence> cond = phaseCondition();
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.rs, state.rv,cond, cells_);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_);
rq_[ actph ].mob = tr_mult * kr / mu;
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure, state.rs, state.rv,cond, cells_);
const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv,cond, cells_);
ADB& head = rq_[ actph ].head;
@ -1911,6 +1919,7 @@ namespace {
ADB
FullyImplicitBlackoilSolver<T>::fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -1918,12 +1927,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, temp, cells);
case Oil: {
return fluid_.muOil(p, rs, cond, cells);
return fluid_.muOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, rv, cond, cells);
return fluid_.muGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -1937,6 +1946,7 @@ namespace {
ADB
FullyImplicitBlackoilSolver<T>::fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
@ -1944,12 +1954,12 @@ namespace {
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, temp, cells);
case Oil: {
return fluid_.bOil(p, rs, cond, cells);
return fluid_.bOil(p, temp, rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, rv, cond, cells);
return fluid_.bGas(p, temp, rv, cond, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -1963,13 +1973,14 @@ namespace {
ADB
FullyImplicitBlackoilSolver<T>::fluidDensity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond,
const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidReciprocFVF(phase, p, rs, rv, cond, cells);
ADB b = fluidReciprocFVF(phase, p, temp, rs, rv, cond, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.

64
opm/autodiff/ImpesTPFAAD.cpp
View File

@ -280,8 +280,9 @@ namespace {
}
V cell_rho_total = V::Zero(nc,1);
const Eigen::Map<const V> p(state.pressure().data(), nc, 1);
const Eigen::Map<const V> T(state.temperature().data(), nc, 1);
for (int phase = 0; phase < np; ++phase) {
const V cell_rho = fluidRho(phase, p, cells);
const V cell_rho = fluidRho(phase, p, T, cells);
const V cell_s = s.col(phase);
cell_rho_total += cell_s * cell_rho;
}
@ -318,10 +319,12 @@ namespace {
// Initialize AD variables: p (cell pressures) and bhp (well bhp).
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V T0 = Eigen::Map<const V>(&state.temperature()[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 T = ADB::constant(T0);
const ADB& bhp = vars[1];
std::vector<int> bpat = p.blockPattern();
@ -331,6 +334,7 @@ namespace {
// Extract variables for perforation cell pressures
// and corresponding perforation well pressures.
const ADB p_perfcell = subset(p, well_cells);
const ADB T_perfcell = subset(T, well_cells);
// Construct matrix to map wells->perforations.
M well_to_perf(well_cells.size(), nw);
typedef Eigen::Triplet<double> Tri;
@ -352,20 +356,20 @@ namespace {
ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat);
qs_ = ADB::constant(V::Zero(nw*np, 1), bpat);
for (int phase = 0; phase < np; ++phase) {
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 ADB cell_b = fluidFvf(phase, p, T, cells);
const ADB cell_rho = fluidRho(phase, p, T, cells);
const ADB well_b = fluidFvf(phase, p_perfwell, T_perfcell, 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 mu = fluidMu(phase, p.value(), T.value(), cells);
const V cell_mob = kr / mu;
const ADB head_diff_grav = (grav_ * cell_rho);
const ADB head = nkgradp + (grav_ * cell_rho);
const UpwindSelector<double> upwind(grid_, ops_, head.value());
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_mu = fluidMu(phase, p_perfwell.value(), T_perfcell.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 * head;
@ -499,9 +503,11 @@ namespace {
const V transw = Eigen::Map<const V>(wells_.WI, nperf, 1);
const V p = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V T = Eigen::Map<const V>(&state.temperature()[0], nc, 1);
const V bhp = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
const V p_perfcell = subset(p, well_cells);
const V T_perfcell = subset(T, well_cells);
const V transi = subset(geo_.transmissibility(),
ops_.internal_faces);
@ -515,15 +521,15 @@ namespace {
V perf_flux = V::Zero(nperf, 1);
for (int phase = 0; phase < np; ++phase) {
const V cell_rho = fluidRho(phase, p, cells);
const V cell_rho = fluidRho(phase, p, T, cells);
const V head = nkgradp + (grav_ * cell_rho.matrix()).array();
const UpwindSelector<double> upwind(grid_, ops_, head);
const V kr = fluidKr(phase);
const V mu = fluidMu(phase, p, cells);
const V mu = fluidMu(phase, p, T, 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, well_cells);
const V well_mu = fluidMu(phase, p_perfwell, T_perfcell, 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);
@ -545,19 +551,19 @@ namespace {
V ImpesTPFAAD::fluidMu(const int phase, const V& p, const std::vector<int>& cells) const
V ImpesTPFAAD::fluidMu(const int phase, const V& p, const V& T, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, T, cells);
case Oil: {
V dummy_rs = V::Zero(p.size(), 1) * p;
std::vector<PhasePresence> cond(dummy_rs.size());
return fluid_.muOil(p, dummy_rs, cond, cells);
return fluid_.muOil(p, T, dummy_rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, cells);
return fluid_.muGas(p, T, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -567,19 +573,19 @@ namespace {
ADB ImpesTPFAAD::fluidMu(const int phase, const ADB& p, const std::vector<int>& cells) const
ADB ImpesTPFAAD::fluidMu(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.muWat(p, cells);
return fluid_.muWat(p, T, cells);
case Oil: {
ADB dummy_rs = V::Zero(p.size(), 1) * p;
std::vector<PhasePresence> cond(dummy_rs.size());
return fluid_.muOil(p, dummy_rs, cond, cells);
return fluid_.muOil(p, T, dummy_rs, cond, cells);
}
case Gas:
return fluid_.muGas(p, cells);
return fluid_.muGas(p, T, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -589,19 +595,19 @@ namespace {
V ImpesTPFAAD::fluidFvf(const int phase, const V& p, const std::vector<int>& cells) const
V ImpesTPFAAD::fluidFvf(const int phase, const V& p, const V& T, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, T, cells);
case Oil: {
V dummy_rs = V::Zero(p.size(), 1) * p;
std::vector<PhasePresence> cond(dummy_rs.size());
return fluid_.bOil(p, dummy_rs, cond, cells);
return fluid_.bOil(p, T, dummy_rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, cells);
return fluid_.bGas(p, T, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -611,19 +617,19 @@ namespace {
ADB ImpesTPFAAD::fluidFvf(const int phase, const ADB& p, const std::vector<int>& cells) const
ADB ImpesTPFAAD::fluidFvf(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const
{
switch (phase) {
case Water:
return fluid_.bWat(p, cells);
return fluid_.bWat(p, T, cells);
case Oil: {
ADB dummy_rs = V::Zero(p.size(), 1) * p;
std::vector<PhasePresence> cond(dummy_rs.size());
return fluid_.bOil(p, dummy_rs, cond, cells);
return fluid_.bOil(p, T, dummy_rs, cond, cells);
}
case Gas:
return fluid_.bGas(p, cells);
return fluid_.bGas(p, T, cells);
default:
OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
}
@ -633,10 +639,10 @@ namespace {
V ImpesTPFAAD::fluidRho(const int phase, const V& p, const std::vector<int>& cells) const
V ImpesTPFAAD::fluidRho(const int phase, const V& p, const V& T, const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
V b = fluidFvf(phase, p, cells);
V b = fluidFvf(phase, p, T, cells);
V rho = V::Constant(p.size(), 1, rhos[phase]) * b;
return rho;
}
@ -645,10 +651,10 @@ namespace {
ADB ImpesTPFAAD::fluidRho(const int phase, const ADB& p, const std::vector<int>& cells) const
ADB ImpesTPFAAD::fluidRho(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
ADB b = fluidFvf(phase, p, cells);
ADB b = fluidFvf(phase, p, T, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
return rho;
}

12
opm/autodiff/ImpesTPFAAD.hpp
View File

@ -104,12 +104,12 @@ namespace Opm {
void computeFluxes(BlackoilState& state, WellState& well_state) const;
// Fluid interface forwarding calls to correct methods of fluid_.
V fluidMu(const int phase, const V& p, const std::vector<int>& cells) const;
ADB fluidMu(const int phase, const ADB& p, const std::vector<int>& cells) const;
V fluidFvf(const int phase, const V& p, const std::vector<int>& cells) const;
ADB fluidFvf(const int phase, const ADB& p, const std::vector<int>& cells) const;
V fluidRho(const int phase, const V& p, const std::vector<int>& cells) const;
ADB fluidRho(const int phase, const ADB& p, const std::vector<int>& cells) const;
V fluidMu(const int phase, const V& p, const V& T, const std::vector<int>& cells) const;
ADB fluidMu(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const;
V fluidFvf(const int phase, const V& p, const V& T, const std::vector<int>& cells) const;
ADB fluidFvf(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const;
V fluidRho(const int phase, const V& p, const V& T, const std::vector<int>& cells) const;
ADB fluidRho(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const;
V fluidKr(const int phase) const;
V fluidKrWell(const int phase) const;
};

2
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -260,7 +260,7 @@ namespace Opm
#endif
M id(Jn[n].rows(), Jn[n].cols());
id.setIdentity();
const M Di = solver.solve(id);
const Eigen::SparseMatrix<typename M::Scalar, Eigen::ColMajor> Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());

53
opm/autodiff/RateConverter.hpp
View File

@ -267,6 +267,7 @@ namespace Opm {
, repcells_(Details::representative<typename Property::Cells>(rmap_))
, ncells_ (Details::countCells(rmap_))
, p_avg_ (rmap_.numRegions())
, T_avg_ (rmap_.numRegions())
, Rmax_ (rmap_.numRegions(), props.numPhases())
{}
@ -327,6 +328,7 @@ namespace Opm {
const PhaseUsage& pu = props_.phaseUsage();
const V& p = getRegPress(r);
const V& T = getRegTemp(r);
const typename Property::Cells& c = getRegCell (r);
const int iw = Details::PhasePos::water(pu);
@ -338,7 +340,7 @@ namespace Opm {
if (Details::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const V& bw = props_.bWat(p, c);
const V& bw = props_.bWat(p, T, c);
coeff[iw] = 1.0 / bw(0);
}
@ -351,7 +353,7 @@ namespace Opm {
if (Details::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const V& bo = props_.bOil(p, m.rs, m.cond, c);
const V& bo = props_.bOil(p, T, m.rs, m.cond, c);
const double den = bo(0) * detR;
coeff[io] += 1.0 / den;
@ -364,7 +366,7 @@ namespace Opm {
if (Details::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const V& bg = props_.bGas(p, m.rv, m.cond, c);
const V& bg = props_.bGas(p, T, m.rv, m.cond, c);
const double den = bg(0) * detR;
coeff[ig] += 1.0 / den;
@ -404,6 +406,11 @@ namespace Opm {
*/
Eigen::ArrayXd p_avg_;
/**
* Average temperature in each FIP region.
*/
Eigen::ArrayXd T_avg_;
/**
* Maximum dissolution and evaporation ratios at average
* hydrocarbon pressure.
@ -474,6 +481,26 @@ namespace Opm {
p_avg_ /= ncells_;
}
/**
* Compute average temperature in all regions.
*
* \param[in] state Dynamic reservoir state.
*/
void
averageTemperature(const BlackoilState& state)
{
T_avg_.setZero();
const std::vector<double>& T = state.temperature();
for (std::vector<double>::size_type
i = 0, n = T.size(); i < n; ++i)
{
T_avg_(rmap_.region(i)) += T[i];
}
T_avg_ /= ncells_;
}
/**
* Compute maximum dissolution and evaporation ratios at
* average hydrocarbon pressure.
@ -499,8 +526,8 @@ namespace Opm {
// pressure into account. This facility uses the
// average *hydrocarbon* pressure rather than
// average phase pressure.
Rmax_.col(io) = props_.rsSat(p_avg_, repcells_);
Rmax_.col(ig) = props_.rvSat(p_avg_, repcells_);
Rmax_.col(io) = props_.rsSat(p_avg_, T_avg_, repcells_);
Rmax_.col(ig) = props_.rvSat(p_avg_, T_avg_, repcells_);
}
}
@ -591,6 +618,22 @@ namespace Opm {
return p;
}
/**
* Retrieve average temperature in region.
*
* \param[in] r Particular region.
*
* \return Average temperature in region \c r.
*/
typename Property::V
getRegTemp(const RegionId r) const
{
typename Property::V T(1);
T << T_avg_(r);
return T;
}
/**
* Retrieve representative cell of region
*

4
opm/autodiff/SimulatorCompressibleAd.cpp
View File

@ -357,7 +357,7 @@ namespace Opm
// Solve pressure equation.
if (check_well_controls_) {
computeFractionalFlow(props_, allcells_,
state.pressure(), state.surfacevol(), state.saturation(),
state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
fractional_flows);
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
@ -446,7 +446,7 @@ namespace Opm
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0],
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.temperature()[0],
&initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
state.saturation(), state.surfacevol());
double substep_injected[2] = { 0.0 };

13
tests/test_boprops_ad.cpp
View File

@ -124,7 +124,10 @@ BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture<SetupSimple>)
Vpw[3] = 8*Opm::unit::barsa;
Vpw[4] = 16*Opm::unit::barsa;
const Opm::BlackoilPropsAd::V VmuWat = boprops_ad.muWat(Vpw, cells);
// standard temperature
V T = V::Constant(cells.size(), 273.15+20);
const Opm::BlackoilPropsAd::V VmuWat = boprops_ad.muWat(Vpw, T, cells);
// Zero pressure dependence in water viscosity
for (V::Index i = 0, n = VmuWat.size(); i < n; ++i) {
@ -149,16 +152,20 @@ BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture<SetupSimple>)
Vpw[3] = 8*Opm::unit::barsa;
Vpw[4] = 16*Opm::unit::barsa;
// standard temperature
V T = V::Constant(cells.size(), 273.15+20);
typedef Opm::BlackoilPropsAd::ADB ADB;
const V VmuWat = boprops_ad.muWat(Vpw, cells);
const V VmuWat = boprops_ad.muWat(Vpw, T, cells);
for (V::Index i = 0, n = Vpw.size(); i < n; ++i) {
const std::vector<int> bp(1, grid.c_grid()->number_of_cells);
const Opm::BlackoilPropsAd::Cells c(1, 0);
const V pw = V(1, 1) * Vpw[i];
const ADB Apw = ADB::variable(0, pw, bp);
const ADB AmuWat = boprops_ad.muWat(Apw, c);
const ADB AT = ADB::constant(T);
const ADB AmuWat = boprops_ad.muWat(Apw, AT, c);
BOOST_CHECK_EQUAL(AmuWat.value()[0], VmuWat[i]);
}