diff --git a/.gitignore b/.gitignore index 2cf3eaafe..585d691bc 100644 --- a/.gitignore +++ b/.gitignore @@ -29,3 +29,6 @@ test_vec # Mac OS X debug info *.dSYM + +# emacs directory setting: +.dir-locals.el diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 45f86373f..3f8c45a06 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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: diff --git a/examples/opm_init_check.cpp b/examples/opm_init_check.cpp new file mode 100644 index 000000000..7b04c1f31 --- /dev/null +++ b/examples/opm_init_check.cpp @@ -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 . +*/ + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include +#include + + +using namespace Opm; + + + + +class CellTrans { +public: + CellTrans(const std::tuple& ijk) : + m_ijk(ijk) + { + + } + + void update(const std::tuple& ijk , double trans) { + auto iter = m_trans.find( ijk ); + if (iter == m_trans.end()) + m_trans.insert( std::pair , double>(ijk , trans)); + else + iter->second *= trans; + } + + int numConnections() const { + return m_trans.size(); + } + + + static void jsonDumpIJK(std::ostream& os , const std::tuple& 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 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 m_ijk; + std::map , 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 = m_transVector[g1]; + if (!cellTrans) { + cellTrans = std::make_shared( 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 = 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 = m_transVector[g]; + if (cellTrans) + count += cellTrans->numConnections(); + } + return count; + } + + + + std::tuple 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(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 = 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 > m_transVector; +}; + + +/*****************************************************************/ + +void initOPMTrans(TransGraph& opmTrans , DeckConstPtr deck , std::shared_ptr eclipseState) { + std::shared_ptr grid = std::make_shared( eclipseState->getEclipseGrid(), eclipseState->getDoubleGridProperty("PORV")->getData()); + const struct UnstructuredGrid * cGrid = grid->c_grid(); + std::shared_ptr 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( ecl_nnc_export_get_size( ecl_grid )); + std::vector 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 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 state = std::make_shared( 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; +} + diff --git a/examples/python/example.py b/examples/python/example.py new file mode 100755 index 000000000..7b181eb38 --- /dev/null +++ b/examples/python/example.py @@ -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) diff --git a/examples/python/trans_graph.py b/examples/python/trans_graph.py new file mode 100644 index 000000000..7a3112d29 --- /dev/null +++ b/examples/python/trans_graph.py @@ -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 + diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 068996267..8036582d9 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -70,7 +70,7 @@ struct HelperOps TwoColInt nbi; extractInternalFaces(grid, internal_faces, nbi); int num_internal=internal_faces.size(); - + // std::cout << "nbi = \n" << nbi << std::endl; // Create matrices. ngrad.resize(num_internal, nc); @@ -189,11 +189,11 @@ namespace { template - Eigen::SparseMatrix + typename AutoDiffBlock::M constructSupersetSparseMatrix(const int full_size, const IntVec& indices) { const int subset_size = indices.size(); - Eigen::SparseMatrix mat(full_size, subset_size); + typename AutoDiffBlock::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 -AutoDiffBlock -subset(const AutoDiffBlock& x, - const IntVec& indices) -{ - Eigen::SparseMatrix sub - = constructSupersetSparseMatrix(x.value().size(), indices).transpose(); - return sub * x; -} - - /// Returns x(indices). template @@ -223,9 +211,25 @@ Eigen::Array subset(const Eigen::Array& x, const IntVec& indices) { - return (constructSupersetSparseMatrix(x.size(), indices).transpose() * x.matrix()).array(); + typedef typename Eigen::Array::Index Index; + const Index size = indices.size(); + Eigen::Array ret( size ); + for( Index i=0; i +AutoDiffBlock +subset(const AutoDiffBlock& x, + const IntVec& indices) +{ + const typename AutoDiffBlock::M sub + = constructSupersetSparseMatrix(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::V& d) /// Returns the input expression, but with all Jacobians collapsed to one. +template inline -AutoDiffBlock -collapseJacs(const AutoDiffBlock& x) +void +collapseJacs(const AutoDiffBlock& x, Matrix& jacobian) { typedef AutoDiffBlock ADB; const int nb = x.numBlocks(); @@ -383,9 +388,21 @@ collapseJacs(const AutoDiffBlock& 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 +collapseJacs(const AutoDiffBlock& x) +{ + typedef AutoDiffBlock ADB; + // Build final jacobian. std::vector 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); } diff --git a/opm/autodiff/BlackoilPropsAd.cpp b/opm/autodiff/BlackoilPropsAd.cpp index 6d323531d..e7da7df3f 100644 --- a/opm/autodiff/BlackoilPropsAd.cpp +++ b/opm/autodiff/BlackoilPropsAd.cpp @@ -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& /*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& /*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 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& 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& 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& /*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& /*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& /*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& /*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)); diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index e7656f6ad..dea27c36f 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -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& 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& 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& 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& 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& 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& 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& 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& cond, const Cells& cells) const; diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index c9ac43dbd..d41a9be6e 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -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& 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& 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& 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& 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& 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& 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& 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& 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); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 16c3af45e..88f2db987 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -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& 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& 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& 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& 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& 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& 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& 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& cond, const Cells& cells) const; diff --git a/opm/autodiff/BlackoilPropsAdInterface.hpp b/opm/autodiff/BlackoilPropsAdInterface.hpp index f15b33ca1..b59720e35 100644 --- a/opm/autodiff/BlackoilPropsAdInterface.hpp +++ b/opm/autodiff/BlackoilPropsAdInterface.hpp @@ -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& 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& 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& 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& 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& 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& 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& cond, const Cells& cells) const = 0; diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h index 03828f10b..1917dddf8 100644 --- a/opm/autodiff/ConservativeSparseSparseProduct.h +++ b/opm/autodiff/ConservativeSparseSparseProduct.h @@ -15,7 +15,10 @@ #include #include #include +<<<<<<< HEAD #include +======= +>>>>>>> 54feee5987a5b84c4f3d3b283d03f574c325d466 #include #include @@ -64,11 +67,14 @@ struct QuickSort< 0 > template 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::type::Scalar Scalar; typedef typename remove_all::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>>>>>> 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 { diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 329980dd9..a419adbd0 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -142,6 +142,7 @@ namespace Opm { struct SolutionState { SolutionState(const int np); ADB pressure; + ADB temperature; std::vector 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& 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& 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& cond, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 560855e8d..445a507af 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -357,6 +357,7 @@ namespace { template FullyImplicitBlackoilSolver::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(& x.temperature()[0], x.temperature().size()); + state.temperature = ADB::constant(temp); + // Saturations const std::vector& 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& 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 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 perf_cond(nperf); const std::vector& pc = phaseCondition(); for (int perf = 0; perf < nperf; ++perf) { @@ -630,21 +638,21 @@ namespace { std::vector rssat_perf(nperf, 0.0); std::vector 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 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; @@ -1910,7 +1918,8 @@ namespace { template ADB FullyImplicitBlackoilSolver::fluidViscosity(const int phase, - const ADB& p , + const ADB& p , + const ADB& temp , const ADB& rs , const ADB& rv , const std::vector& 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::fluidReciprocFVF(const int phase, const ADB& p , + const ADB& temp , const ADB& rs , const ADB& rv , const std::vector& 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); } @@ -1962,14 +1972,15 @@ namespace { template ADB FullyImplicitBlackoilSolver::fluidDensity(const int phase, - const ADB& p , + const ADB& p , + const ADB& temp , const ADB& rs , const ADB& rv , const std::vector& cond, const std::vector& 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. diff --git a/opm/autodiff/ImpesTPFAAD.cpp b/opm/autodiff/ImpesTPFAAD.cpp index e25e78fed..697f80a39 100644 --- a/opm/autodiff/ImpesTPFAAD.cpp +++ b/opm/autodiff/ImpesTPFAAD.cpp @@ -280,8 +280,9 @@ namespace { } V cell_rho_total = V::Zero(nc,1); const Eigen::Map p(state.pressure().data(), nc, 1); + const Eigen::Map 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(&state.pressure()[0], nc, 1); + const V T0 = Eigen::Map(&state.temperature()[0], nc, 1); const V bhp0 = Eigen::Map(&well_state.bhp()[0], nw, 1); std::vector vars0 = { p0, bhp0 }; std::vector vars = ADB::variables(vars0); const ADB& p = vars[0]; + const ADB T = ADB::constant(T0); const ADB& bhp = vars[1]; std::vector 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 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 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(wells_.WI, nperf, 1); const V p = Eigen::Map(&state.pressure()[0], nc, 1); + const V T = Eigen::Map(&state.temperature()[0], nc, 1); const V bhp = Eigen::Map(&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 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& cells) const + V ImpesTPFAAD::fluidMu(const int phase, const V& p, const V& T, const std::vector& 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 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& cells) const + ADB ImpesTPFAAD::fluidMu(const int phase, const ADB& p, const ADB& T, const std::vector& 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 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& cells) const + V ImpesTPFAAD::fluidFvf(const int phase, const V& p, const V& T, const std::vector& 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 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& cells) const + ADB ImpesTPFAAD::fluidFvf(const int phase, const ADB& p, const ADB& T, const std::vector& 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 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& cells) const + V ImpesTPFAAD::fluidRho(const int phase, const V& p, const V& T, const std::vector& 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& cells) const + ADB ImpesTPFAAD::fluidRho(const int phase, const ADB& p, const ADB& T, const std::vector& 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; } diff --git a/opm/autodiff/ImpesTPFAAD.hpp b/opm/autodiff/ImpesTPFAAD.hpp index edb2480e3..aa137f4ea 100644 --- a/opm/autodiff/ImpesTPFAAD.hpp +++ b/opm/autodiff/ImpesTPFAAD.hpp @@ -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& cells) const; - ADB fluidMu(const int phase, const ADB& p, const std::vector& cells) const; - V fluidFvf(const int phase, const V& p, const std::vector& cells) const; - ADB fluidFvf(const int phase, const ADB& p, const std::vector& cells) const; - V fluidRho(const int phase, const V& p, const std::vector& cells) const; - ADB fluidRho(const int phase, const ADB& p, const std::vector& cells) const; + V fluidMu(const int phase, const V& p, const V& T, const std::vector& cells) const; + ADB fluidMu(const int phase, const ADB& p, const ADB& T, const std::vector& cells) const; + V fluidFvf(const int phase, const V& p, const V& T, const std::vector& cells) const; + ADB fluidFvf(const int phase, const ADB& p, const ADB& T, const std::vector& cells) const; + V fluidRho(const int phase, const V& p, const V& T, const std::vector& cells) const; + ADB fluidRho(const int phase, const ADB& p, const ADB& T, const std::vector& cells) const; V fluidKr(const int phase) const; V fluidKrWell(const int phase) const; }; diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 6d6058a59..7a4d3500a 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -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 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()); diff --git a/opm/autodiff/RateConverter.hpp b/opm/autodiff/RateConverter.hpp index 23d69ae51..d349bae06 100644 --- a/opm/autodiff/RateConverter.hpp +++ b/opm/autodiff/RateConverter.hpp @@ -267,6 +267,7 @@ namespace Opm { , repcells_(Details::representative(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& T = state.temperature(); + for (std::vector::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 * diff --git a/opm/autodiff/SimulatorCompressibleAd.cpp b/opm/autodiff/SimulatorCompressibleAd.cpp index 80f778fa1..a97d2118a 100644 --- a/opm/autodiff/SimulatorCompressibleAd.cpp +++ b/opm/autodiff/SimulatorCompressibleAd.cpp @@ -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 }; diff --git a/tests/test_boprops_ad.cpp b/tests/test_boprops_ad.cpp index d92fe0992..9c06b5c32 100644 --- a/tests/test_boprops_ad.cpp +++ b/tests/test_boprops_ad.cpp @@ -124,7 +124,10 @@ BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture) 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) 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 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]); }