diff --git a/.gitignore b/.gitignore
index 2cf3eaafe..585d691bc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,3 +29,6 @@ test_vec
# Mac OS X debug info
*.dSYM
+
+# emacs directory setting:
+.dir-locals.el
diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 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;
+}
+
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
+