From e0d39027d009cacb4b732928bb611240e5db89fc Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Wed, 26 Nov 2014 14:53:35 +0100 Subject: [PATCH 01/12] Added ignore to emacs per directory configuration file. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) 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 From ad92edab89b7e63c3375db255e8fc859de926302 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Wed, 26 Nov 2014 14:54:20 +0100 Subject: [PATCH 02/12] Added opm_init_check to compare OPM and Eclipse. --- CMakeLists_files.cmake | 2 + examples/opm_init_check.cpp | 362 ++++++++++++++++++++++++++++++++++++ 2 files changed, 364 insertions(+) create mode 100644 examples/opm_init_check.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index d3ee985be..cb5456350 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..f1995e7fe --- /dev/null +++ b/examples/opm_init_check.cpp @@ -0,0 +1,362 @@ +/* + 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 + + +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; +} + From 29a0c723ef5cd0402b31211943973d903e5e295e Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Wed, 26 Nov 2014 14:55:45 +0100 Subject: [PATCH 03/12] Python class TransGraph to load json from opm_init_check --- examples/python/example.py | 136 ++++++++++++++++++++++ examples/python/trans_graph.py | 198 +++++++++++++++++++++++++++++++++ 2 files changed, 334 insertions(+) create mode 100755 examples/python/example.py create mode 100644 examples/python/trans_graph.py 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 + From c69b2fd8edc75c027b6e08839350bbe0ac943f78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 27 Nov 2014 15:12:56 +0100 Subject: [PATCH 04/12] Include to fix build on GCC 4.8 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit GCC 4.8 provides 'nullptr' in C++11 mode. That capability is detected at configuration time and stored in and we need it when targeting Dune 2.2.1 to avoid diagnostics of the form [...]/dune/common/nullptr.hh:27:1: error: expected ‘;’ after class definition } nullptr = {}; // and whose name is nullptr ^ [...]/dune/common/nullptr.hh:27:1: error: qualifiers can only be specified for objects and functions [...]/dune/common/nullptr.hh:27:3: error: expected unqualified-id before ‘nullptr’ } nullptr = {}; // and whose name is nullptr --- examples/opm_init_check.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/opm_init_check.cpp b/examples/opm_init_check.cpp index f1995e7fe..7b04c1f31 100644 --- a/examples/opm_init_check.cpp +++ b/examples/opm_init_check.cpp @@ -17,6 +17,7 @@ along with OPM. If not, see . */ +#include #include #include From 4e3a69cc907499b4fac93f94f8887b58df9ecdb9 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 20 Nov 2014 12:31:50 +0100 Subject: [PATCH 05/12] PVT properties: allow them to be temperature dependent Note that this patch does not introduce any real temperature dependence but only changes the APIs for the viscosity and for the density related methods. Note that I also don't like the fact that this requires so many changes to so many files, but with the current design of the property classes I cannot see a way to avoid this... --- opm/autodiff/BlackoilPropsAd.cpp | 66 ++++++++++++++----- opm/autodiff/BlackoilPropsAd.hpp | 32 +++++++++ opm/autodiff/BlackoilPropsAdFromDeck.cpp | 64 +++++++++++++----- opm/autodiff/BlackoilPropsAdFromDeck.hpp | 32 +++++++++ opm/autodiff/BlackoilPropsAdInterface.hpp | 30 +++++++++ opm/autodiff/FullyImplicitBlackoilSolver.hpp | 4 ++ .../FullyImplicitBlackoilSolver_impl.hpp | 41 +++++++----- opm/autodiff/ImpesTPFAAD.cpp | 64 ++++++++++-------- opm/autodiff/ImpesTPFAAD.hpp | 12 ++-- opm/autodiff/RateConverter.hpp | 53 +++++++++++++-- opm/autodiff/SimulatorCompressibleAd.cpp | 4 +- tests/test_boprops_ad.cpp | 15 ++++- 12 files changed, 324 insertions(+), 93 deletions(-) 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/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/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..0dd598bb9 100644 --- a/tests/test_boprops_ad.cpp +++ b/tests/test_boprops_ad.cpp @@ -124,7 +124,11 @@ 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; + T.resize(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 +153,21 @@ BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture) Vpw[3] = 8*Opm::unit::barsa; Vpw[4] = 16*Opm::unit::barsa; + // standard temperature + V T; + T.resize(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]); } From 7b2334b75844c12f17ad9f3cdb1cd08908f1605d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 2 Dec 2014 09:24:03 +0100 Subject: [PATCH 06/12] Fix wrong resize() usage for Eigen::Array. --- tests/test_boprops_ad.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_boprops_ad.cpp b/tests/test_boprops_ad.cpp index 0dd598bb9..9c06b5c32 100644 --- a/tests/test_boprops_ad.cpp +++ b/tests/test_boprops_ad.cpp @@ -125,8 +125,7 @@ BOOST_FIXTURE_TEST_CASE(ViscosityValue, TestFixture) Vpw[4] = 16*Opm::unit::barsa; // standard temperature - V T; - T.resize(cells.size(), 273.15+20); + V T = V::Constant(cells.size(), 273.15+20); const Opm::BlackoilPropsAd::V VmuWat = boprops_ad.muWat(Vpw, T, cells); @@ -154,8 +153,7 @@ BOOST_FIXTURE_TEST_CASE(ViscosityAD, TestFixture) Vpw[4] = 16*Opm::unit::barsa; // standard temperature - V T; - T.resize(cells.size(), 273.15+20); + V T = V::Constant(cells.size(), 273.15+20); typedef Opm::BlackoilPropsAd::ADB ADB; From c51a794cac3028222b8bb44df1afc5a247cc3c46 Mon Sep 17 00:00:00 2001 From: Robert K Date: Tue, 2 Dec 2014 10:38:21 +0100 Subject: [PATCH 07/12] overloaded ConservativeSparseSparseProduct to speed up matrix-matrix multiplication. --- CMakeLists_files.cmake | 1 + opm/autodiff/AutoDiffBlock.hpp | 9 +- .../ConservativeSparseSparseProduct.h | 318 ++++++++++++++++++ 3 files changed, 324 insertions(+), 4 deletions(-) create mode 100644 opm/autodiff/ConservativeSparseSparseProduct.h diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index cb5456350..3f8c45a06 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -98,6 +98,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilPropsAdFromDeck.hpp opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/CPRPreconditioner.hpp + opm/autodiff/ConservativeSparseSparseProduct.h opm/autodiff/DuneMatrix.hpp opm/autodiff/GeoProps.hpp opm/autodiff/GridHelpers.hpp diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index c053c6536..19dfd5061 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -22,6 +22,7 @@ #include +#include #include #include @@ -102,7 +103,7 @@ namespace Opm } /// Create an AutoDiffBlock representing a constant. - /// \param[in] val values + /// \param[in] val values static AutoDiffBlock constant(const V& val) { return AutoDiffBlock(val); @@ -112,7 +113,7 @@ namespace Opm /// This variant requires specifying the block sizes used /// for the Jacobians even though the Jacobian matrices /// themselves will be zero. - /// \param[in] val values + /// \param[in] val values /// \param[in] blocksizes block pattern static AutoDiffBlock constant(const V& val, const std::vector& blocksizes) { @@ -129,7 +130,7 @@ namespace Opm /// Create an AutoDiffBlock representing a single variable block. /// \param[in] index index of the variable you are constructing - /// \param[in] val values + /// \param[in] val values /// \param[in] blocksizes block pattern /// The resulting object will have size() equal to block_pattern[index]. /// Its jacobians will all be zero, except for derivative()[index], which @@ -154,7 +155,7 @@ namespace Opm } /// Create an AutoDiffBlock by directly specifying values and jacobians. - /// \param[in] val values + /// \param[in] val values /// \param[in] jac vector of jacobians static AutoDiffBlock function(const V& val, const std::vector& jac) { diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h new file mode 100644 index 000000000..b1b036913 --- /dev/null +++ b/opm/autodiff/ConservativeSparseSparseProduct.h @@ -0,0 +1,318 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2011 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H +#define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H + +#warning "Using overloaded Eigen::ConservativeSparseSparseProduct.h" + +#include +#include +#include + +#include + +namespace Eigen { + +// forward declaration of SparseMatrix +template +class SparseMatrix; + + +namespace internal { + +template < unsigned int depth > +struct QuickSort +{ + template + static inline void sort(T begin, T end) + { + if (begin != end) + { + T middle = std::partition (begin, end, + std::bind2nd(std::less::value_type>(), *begin) + ); + QuickSort< depth-1 >::sort(begin, middle); + + // std::sort (max(begin + 1, middle), end); + T new_middle = begin; + QuickSort< depth-1 >::sort(++new_middle, end); + } + } +}; + +template <> +struct QuickSort< 0 > +{ + template + static inline void sort(T begin, T end) + { + std::sort( begin, end ); + } +}; + + +template +static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) +{ + typedef typename remove_all::type::Scalar Scalar; + typedef typename remove_all::type::Index Index; + + // make sure to call innerSize/outerSize since we fake the storage order. + Index rows = lhs.innerSize(); + Index cols = rhs.outerSize(); + eigen_assert(lhs.outerSize() == rhs.innerSize()); + + std::vector mask(rows,false); + Matrix values(rows); + Matrix indices(rows); + + // estimate the number of non zero entries + // given a rhs column containing Y non zeros, we assume that the respective Y columns + // of the lhs differs in average of one non zeros, thus the number of non zeros for + // the product of a rhs column with the lhs is X+Y where X is the average number of non zero + // per column of the lhs. + // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs) + Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros(); + + res.setZero(); + res.reserve(Index(estimated_nnz_prod)); + // we compute each column of the result, one after the other + for (Index j=0; j 1 ) + { + // sort indices for sorted insertion to avoid later copying + QuickSort< 1 >::sort( indices.data(), indices.data()+nnz ); + } + + // ordered insertion + // still using insertBackByOuterInnerUnordered since we know what we are doing + for(Index k=0; k use a quick sort + // otherwise => loop through the entire vector + // In order to avoid to perform an expensive log2 when the + // result is clearly very sparse we use a linear bound up to 200. + //if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); + for(Index k=0; k::Flags&RowMajorBit) ? RowMajor : ColMajor, + int RhsStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor, + int ResStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor> +struct conservative_sparse_sparse_product_selector; + +template +struct conservative_sparse_sparse_product_selector +{ + typedef typename remove_all::type LhsCleaned; + typedef typename LhsCleaned::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + //typedef SparseMatrix RowMajorMatrix; + typedef SparseMatrix ColMajorMatrix; + //ColMajorMatrix resCol(lhs.rows(),rhs.cols()); + res = ColMajorMatrix(lhs.rows(),rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhs, res); + //internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); + // sort the non zeros: + //RowMajorMatrix resRow(resCol); + //res = resRow; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + //RowMajorMatrix rhsRow = rhs; + //RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + ColMajorMatrix lhsCol = lhs; + res = ResultType( lhs.rows(), rhs.cols() ); + internal::conservative_sparse_sparse_product_impl( lhsCol, rhs, res ); + //internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, resRow); + //res = resRow; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix rhsCol = rhs; + res = ResultType( lhs.rows(), rhs.cols() ); + internal::conservative_sparse_sparse_product_impl( lhs, rhsCol, res); + /* + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix lhsRow = lhs; + RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, resRow); + res = resRow; + */ + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); + res = resRow; + } +}; + + +template +struct conservative_sparse_sparse_product_selector +{ + typedef typename traits::type>::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); + res = resCol; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix rhsRow = rhs; + res = ResultType( lhs.rows(), rhs.cols() ); + internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, res); + /* + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix lhsCol = lhs; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhsCol, rhs, resCol); + res = resCol; + */ + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix lhsRow = lhs; + res = RowMajorMatrix( lhs.rows(), rhs.cols() ); + internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, res); + /* + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix rhsCol = rhs; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhsCol, resCol); + res = resCol; + */ + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + //typedef SparseMatrix ColMajorMatrix; + res = RowMajorMatrix( lhs.rows(),rhs.cols() ); + //RowMajorMatrix resRow(lhs.rows(),rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhs, res); + // sort the non zeros: + //ColMajorMatrix resCol(resRow); + //res = resCol; + } +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H From a266e98bac1af082209bd7f2b6f6672bab6dca1f Mon Sep 17 00:00:00 2001 From: Robert K Date: Tue, 2 Dec 2014 10:49:42 +0100 Subject: [PATCH 08/12] added some comment about sort. --- opm/autodiff/ConservativeSparseSparseProduct.h | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h index b1b036913..7320c8e90 100644 --- a/opm/autodiff/ConservativeSparseSparseProduct.h +++ b/opm/autodiff/ConservativeSparseSparseProduct.h @@ -53,6 +53,7 @@ struct QuickSort< 0 > template static inline void sort(T begin, T end) { + // fall back to standard insertion sort std::sort( begin, end ); } }; From 2ac6a211b232c0697d53c5396b520ee333ef735d Mon Sep 17 00:00:00 2001 From: Robert K Date: Tue, 2 Dec 2014 10:12:27 +0100 Subject: [PATCH 09/12] use correct types of SparseMatrices. --- opm/autodiff/AutoDiffHelpers.hpp | 34 ++++++++++++--------- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 2 +- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 068996267..7bc3c87ca 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -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. 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()); From efe8ee79f9128ed41b2f0a7d98d6461c95d2dd64 Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 1 Dec 2014 10:27:44 +0100 Subject: [PATCH 10/12] added collapseJacs method. --- opm/autodiff/AutoDiffHelpers.hpp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 7bc3c87ca..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); @@ -361,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(); @@ -387,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); } From 054d4f4dcbcc25f6c8bddbde6ecc3e8384506224 Mon Sep 17 00:00:00 2001 From: Robert K Date: Tue, 2 Dec 2014 14:15:17 +0100 Subject: [PATCH 11/12] include vector to make compile on proprietary systems. --- opm/autodiff/ConservativeSparseSparseProduct.h | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h index 7320c8e90..ca4221622 100644 --- a/opm/autodiff/ConservativeSparseSparseProduct.h +++ b/opm/autodiff/ConservativeSparseSparseProduct.h @@ -15,6 +15,7 @@ #include #include #include +#include #include From 54feee5987a5b84c4f3d3b283d03f574c325d466 Mon Sep 17 00:00:00 2001 From: Robert K Date: Tue, 2 Dec 2014 14:57:00 +0100 Subject: [PATCH 12/12] avoid the multiplication with zero matrix entries. --- .../ConservativeSparseSparseProduct.h | 28 +++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/opm/autodiff/ConservativeSparseSparseProduct.h b/opm/autodiff/ConservativeSparseSparseProduct.h index ca4221622..42c6ef750 100644 --- a/opm/autodiff/ConservativeSparseSparseProduct.h +++ b/opm/autodiff/ConservativeSparseSparseProduct.h @@ -88,24 +88,29 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r // we compute each column of the result, one after the other for (Index j=0; j::sort( indices.data(), indices.data()+nnz ); } + res.startVec(j); // ordered insertion // still using insertBackByOuterInnerUnordered since we know what we are doing for(Index k=0; k