Merge pull request #738 from joakim-hove/write-tran

Write TRAN? to the INIT file.

CMakeLists.txt

@@ -149,5 +149,11 @@ if (HAVE_OPM_DATA)
     set_tests_properties(compare_restart_files PROPERTIES DEPENDS flow_SPE1CASE2_restart) # Compares the restart files from tests flow_SPE1CASE2_restart and flow_SPE1CASE2
 
     add_test( NAME flow_sequential_SPE1 COMMAND flow_sequential ${OPM_DATA_ROOT}/spe1/SPE1CASE1.DATA )
+    if (ERT_PYTHON_PATH)
+        include(OpmPythonTest)
+
+        opm_add_python_test( check_INIT_SPE1 ${PROJECT_SOURCE_DIR}/tests/compare_INIT.py $<TARGET_FILE:flow> ${OPM_DATA_ROOT}/spe1/SPE1CASE1.DATA ${OPM_DATA_ROOT}/spe1/eclipse-simulation/SPE1CASE1.INIT TRANX TRANY TRANZ PORO PORV PERMX )
+        opm_add_python_test( check_INIT_NORNE ${PROJECT_SOURCE_DIR}/tests/compare_INIT.py $<TARGET_FILE:flow> ${OPM_DATA_ROOT}/norne/NORNE_ATW2013.DATA ${OPM_DATA_ROOT}/norne/ECL.2014.2/NORNE_ATW2013.INIT PORO PORV PERMX )
+    endif()
 
 endif()

opm/autodiff/FlowMain.hpp

@@ -746,7 +746,7 @@ namespace Opm
                     fullReport.reportParam(tot_os);
                 }
             } else {
-                output_writer_->writeInit();
+                output_writer_->writeInit(geoprops_->simProps(grid_init_->grid()) , geoprops_->nonCartesianConnections( ));
                 if (output_cout_) {
                     std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
                 }

opm/autodiff/GeoProps.hpp

@@ -33,6 +33,7 @@
 #include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
 #include <opm/core/grid/PinchProcessor.hpp>
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <opm/output/Cells.hpp>
 
 #include <Eigen/Eigen>
 
@@ -218,6 +219,71 @@ namespace Opm
         const NNC& nnc() const { return nnc_;}
         const NNC& nonCartesianConnections() const { return noncartesian_;}
 
+
+        /// Most properties are loaded by the parser, and managed by
+        /// the EclipseState class in the opm-parser. However - some
+        /// properties must be calculated by the simulator, the
+        /// purpose of this method is to calculate these properties in
+        /// a form suitable for output. Currently the transmissibility
+        /// is the only property calculated this way:
+        ///
+        /// The grid properties TRANX, TRANY and TRANZ are initialized
+        /// in a form suitable for writing to the INIT file. These
+        /// properties should be interpreted with a
+        /// 'the-grid-is-nearly-cartesian' mindset:
+        ///
+        ///   TRANX[i,j,k] = T on face between cells (i,j,k) and (i+1,j  ,k  )
+        ///   TRANY[i,j,k] = T on face between cells (i,j,k) and (i  ,j+1,k  )
+        ///   TRANZ[i,j,k] = T on face between cells (i,j,k) and (i  ,j  ,k+1)
+        ///
+        /// If the grid structure has no resemblance to a cartesian
+        /// grid the whole TRAN keyword is quite meaningless.
+
+        template <class Grid>
+        const std::vector<data::CellData> simProps( const Grid& grid ) const {
+            using namespace UgGridHelpers;
+            const int* dims = cartDims( grid );
+            const int globalSize = dims[0] * dims[1] * dims[2];
+            const auto& trans = this->transmissibility( );
+
+            data::CellData tranx = {"TRANX" , UnitSystem::measure::transmissibility, std::vector<double>( globalSize )};
+            data::CellData trany = {"TRANY" , UnitSystem::measure::transmissibility, std::vector<double>( globalSize )};
+            data::CellData tranz = {"TRANZ" , UnitSystem::measure::transmissibility, std::vector<double>( globalSize )};
+
+            size_t num_faces = numFaces(grid);
+            auto fc = faceCells(grid);
+            for (size_t i = 0; i < num_faces; ++i) {
+                auto c1 = std::min( fc(i,0) , fc(i,1));
+                auto c2 = std::max( fc(i,0) , fc(i,1));
+
+                if (c1 == -1 || c2 == -1)
+                    continue;
+
+                c1 = globalCell(grid) ? globalCell(grid)[c1] : c1;
+                c2 = globalCell(grid) ? globalCell(grid)[c2] : c2;
+
+                if ((c2 - c1) == 1) {
+                    tranx.data[c1] = trans[i];
+                }
+
+                if ((c2 - c1) == dims[0]) {
+                    trany.data[c1] = trans[i];
+                }
+
+                if ((c2 - c1) == dims[0]*dims[1]) {
+                    tranz.data[c1] = trans[i];
+                }
+            }
+
+            std::vector<data::CellData> tran;
+            tran.push_back( std::move( tranx ));
+            tran.push_back( std::move( trany ));
+            tran.push_back( std::move( tranz ));
+
+            return tran;
+        }
+
+
     private:
         template <class Grid>
         void multiplyHalfIntersections_(const Grid &grid,

opm/autodiff/SimulatorBase_impl.hpp

@@ -19,6 +19,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <utility>
 #include <algorithm>
 #include <locale>
 #include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
@@ -109,7 +110,9 @@ namespace Opm
             adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
         }
 
-        output_writer_.writeInit();
+
+
+        output_writer_.writeInit( geo_.simProps(grid_) , geo_.nonCartesianConnections( ) );
 
         std::string restorefilename = param_.getDefault("restorefile", std::string("") );
         if( ! restorefilename.empty() )

opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp

@@ -73,7 +73,7 @@ namespace Opm
         }
 
         // init output writer
-        output_writer_.writeInit();
+        output_writer_.writeInit( geo_.simProps(grid_) , geo_.nonCartesianConnections( ) );
 
         std::string restorefilename = param_.getDefault("restorefile", std::string("") );
         if( ! restorefilename.empty() )

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -247,14 +247,15 @@ namespace Opm
 
     void
     BlackoilOutputWriter::
-    writeInit()
+    writeInit(const std::vector<data::CellData>& simProps,  const NNC& nnc)
     {
         if( eclWriter_ ) {
-            eclWriter_->writeInit();
+            eclWriter_->writeInitAndEgrid(simProps, nnc);
         }
     }
 
 
+
     namespace detail {
 
         struct WriterCall : public ThreadHandle :: ObjectInterface

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -29,6 +29,7 @@
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
+#include <opm/output/Cells.hpp>
 #include <opm/output/OutputWriter.hpp>
 #include <opm/output/eclipse/EclipseWriter.hpp>
 
@@ -215,7 +216,7 @@ namespace Opm
                              const double* permeability );
 
         /** \copydoc Opm::OutputWriter::writeInit */
-        void writeInit();
+        void writeInit(const std::vector<data::CellData>& simProps, const NNC& nnc);
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
@@ -301,8 +302,7 @@ namespace Opm
                     param.getDefault("output_ecl", true) ?
                     new EclipseWriter(eclipseState,
                                       parallelOutput_->numCells(),
-                                      parallelOutput_->globalCell(),
+                                      parallelOutput_->globalCell())
-                                      nnc )
                    : 0 ),
         eclipseState_(eclipseState),
         asyncOutput_()

tests/compare_INIT.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+import subprocess
+import sys
+import os.path
+
+from ert.ecl import EclFile
+from ert.test import TestAreaContext
+
+
+def compare_files( flow_file , ref_file , kw_list):
+    flow = EclFile( flow_file )
+    ref = EclFile( ref_file )
+
+    for kw in kw_list:
+        flow_kw = flow[kw][0]
+        ref_kw = ref[kw][0]
+
+        if not flow_kw.equal_numeric( ref_kw , epsilon = 1e-3 ):
+            sys.exit("Keyword:%s was different in flow simulation and reference")
+            
+
+        
+
+
+
+
+flow = sys.argv[1]
+data_file = sys.argv[2]
+ref_init_file = sys.argv[3]
+kw_list = sys.argv[4:]
+
+with TestAreaContext("flow_init") as ta:
+    subprocess.check_call( [flow , "nosim=true" , data_file ] )
+
+    compare_files( os.path.splitext( os.path.basename( data_file ) )[0] + ".INIT" , ref_init_file , kw_list)
+
+