Merge pull request #6 from atgeirr/master

Improved output facilities

Makefile.am

@@ -7,28 +7,29 @@ SUBDIRS = . examples
 lib_LTLIBRARIES = libopmpolymer.la
 
 libopmpolymer_la_SOURCES =				\
-opm/polymer/IncompTpfaPolymer.cpp			\
 opm/polymer/CompressibleTpfaPolymer.cpp			\
-opm/polymer/SimulatorPolymer.cpp			\
-opm/polymer/SimulatorCompressiblePolymer.cpp	        \
-opm/polymer/TransportModelPolymer.cpp			\
-opm/polymer/TransportModelCompressiblePolymer.cpp       \
+opm/polymer/IncompTpfaPolymer.cpp			\
 opm/polymer/PolymerInflow.cpp				\
 opm/polymer/PolymerProperties.cpp			\
+opm/polymer/SimulatorCompressiblePolymer.cpp		\
+opm/polymer/SimulatorPolymer.cpp			\
+opm/polymer/TransportModelCompressiblePolymer.cpp	\
+opm/polymer/TransportModelPolymer.cpp			\
 opm/polymer/polymerUtilities.cpp
 
 nobase_include_HEADERS =				\
+opm/polymer/CompressibleTpfaPolymer.hpp			\
 opm/polymer/GravityColumnSolverPolymer.hpp		\
 opm/polymer/GravityColumnSolverPolymer_impl.hpp		\
 opm/polymer/IncompPropertiesDefaultPolymer.hpp		\
 opm/polymer/IncompTpfaPolymer.hpp			\
-opm/polymer/CompressibleTpfaPolymer.hpp			\
+opm/polymer/PolymerBlackoilState.hpp			\
 opm/polymer/PolymerInflow.hpp				\
 opm/polymer/PolymerProperties.hpp			\
 opm/polymer/PolymerState.hpp				\
-opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp	\
-opm/polymer/SimulatorPolymer.hpp			\
 opm/polymer/SimulatorCompressiblePolymer.hpp		\
+opm/polymer/SimulatorPolymer.hpp			\
+opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp	\
+opm/polymer/TransportModelCompressiblePolymer.hpp	\
 opm/polymer/TransportModelPolymer.hpp			\
-opm/polymer/TransportModelCompressiblePolymer.hpp       \
 opm/polymer/polymerUtilities.hpp

configure.ac

@@ -39,6 +39,8 @@ AX_BOOST_DATE_TIME
 AX_BOOST_FILESYSTEM
 AX_BOOST_UNIT_TEST_FRAMEWORK
 
+ERT
+
 # Checks for header files.
 AC_CHECK_HEADERS([opm/core/grid.h], [],
 [AC_MSG_ERROR(

m4/ert.m4

@@ -0,0 +1,67 @@
+AC_DEFUN([_ERT_SOURCE_TEXT],
+[
+AC_LANG_PROGRAM(
+[[
+#include <stddef.h>
+#include <ecl_util.h>
+]],dnl
+[[
+int sz;
+sz = ecl_util_get_sizeof_ctype(ECL_INT_TYPE);
+]])[]dnl
+])[]dnl
+
+# ----------------------------------------------------------------------
+
+AC_DEFUN([ERT],
+[
+AC_ARG_WITH([ert],
+            [AS_HELP_STRING([--with-ert=<root>], [Use ERT libraries])],
+            [], [with_ert=no])
+
+use_ert=no
+
+AS_IF([test x"${with_ert}" != x"no"],
+[
+  _ert_LDFLAGS_SAVE="${LDFLAGS}"
+  _ert_LIBS_SAVE="${LIBS}"
+  _ert_CPPFLAGS_SAVE="${CPPFLAGS}"
+  _ert_CFLAGS_SAVE="${CFLAGS}"
+
+  ERT_CPPFLAGS=
+  ERT_LDFLAGS=
+  ERT_LIBS="-lecl -lgeometry -lert_util -lpthread -lz -lgomp"
+  AS_IF([test x"${with_ert}" != x"yes"],
+    [ERT_LDFLAGS="-L${with_ert}/lib"
+     ERT_CPPFLAGS="-I${with_ert}/include"], [:])[]dnl
+
+  CFLAGS="-std=gnu99"
+  CPPFLAGS="${ERT_CPPFLAGS} ${CPPFLAGS}"
+  LDFLAGS="${ERT_LDFLAGS} ${LDFLAGS}"
+  LIBS="${ERT_LIBS} ${LIBS}"
+
+  AC_LINK_IFELSE([_ERT_SOURCE_TEXT],
+    [use_ert=yes], [use_ert=no])
+
+  LIBS="${_ert_LIBS_SAVE}"
+  CPPFLAGS="${_ert_CPPFLAGS_SAVE}"
+  LDFLAGS="${_ert_LDFLAGS_SAVE}"
+  CFLAGS="${_ert_CFLAGS_SAVE}"
+
+  AS_IF([test x"${use_ert}" = x"yes"],
+   [AC_SUBST([ERT_CPPFLAGS])
+    AC_SUBST([ERT_LDFLAGS])
+    AC_SUBST([ERT_LIBS])
+    AC_DEFINE([HAVE_ERT], [1],
+              [Are the ERT libraries available for reading and writing ECLIPSE files.])],dnl
+   [:])[]dnl
+], [:])[]dnl
+
+AM_CONDITIONAL([HAVE_ERT], [test x"${use_ert}" != x"no"])
+
+# AC_MSG_ERROR(
+# [**** ERT_CPPFLAGS = ${ERT_CPPFLAGS} ****
+# **** ERT_LDFLAGS = ${ERT_LDFLAGS} ****
+# **** ERT_LIBS = ${ERT_LIBS} ****
+# ])
+])

opm/polymer/SimulatorCompressiblePolymer.cpp

@@ -541,7 +541,7 @@ namespace Opm
             catch (...) {
                 THROW("Creating directories failed: " << fpath);
             }
-            vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+            vtkfilename << "/output-" << std::setw(5) << std::setfill('0') << step << ".vtu";
             std::ofstream vtkfile(vtkfilename.str().c_str());
             if (!vtkfile) {
                 THROW("Failed to open " << vtkfilename.str());
@@ -584,7 +584,7 @@ namespace Opm
                 catch (...) {
                     THROW("Creating directories failed: " << fpath);
                 }
-                fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+                fname << "/" << std::setw(5) << std::setfill('0') << step << ".txt";
                 std::ofstream file(fname.str().c_str());
                 if (!file) {
                     THROW("Failed to open " << fname.str());

opm/polymer/SimulatorPolymer.cpp

@@ -59,6 +59,10 @@
 #include <numeric>
 #include <fstream>
 
+#ifdef HAVE_ERT
+#include <opm/core/utility/writeECLData.hpp>
+#endif
+
 
 namespace Opm
 {
@@ -71,6 +75,10 @@ namespace Opm
                             const Opm::PolymerState& state,
                             const int step,
                             const std::string& output_dir);
+        void outputStateBinary(const UnstructuredGrid& grid,
+                               const Opm::PolymerState& state,
+                               const SimulatorTimer& simtimer,
+                               const std::string& output_dir);
         void outputStateMatlab(const UnstructuredGrid& grid,
                                const Opm::PolymerState& state,
                                const int step,
@@ -111,6 +119,7 @@ namespace Opm
         // Parameters for output.
         bool output_;
         bool output_vtk_;
+        bool output_binary_;
         std::string output_dir_;
         int output_interval_;
         // Parameters for well control
@@ -207,6 +216,12 @@ namespace Opm
         output_ = param.getDefault("output", true);
         if (output_) {
             output_vtk_ = param.getDefault("output_vtk", true);
+            output_binary_ = param.getDefault("output_binary", false);
+#ifndef HAVE_ERT
+            if (output_binary_) {
+                THROW("Cannot make binary output without ert library support. Reconfigure opm-core and opm-polymer with --with-ert and recompile.");
+            }
+#endif
             output_dir_ = param.getDefault("output_dir", std::string("output"));
             // Ensure that output dir exists
             boost::filesystem::path fpath(output_dir_);
@@ -307,6 +322,9 @@ namespace Opm
                 if (output_vtk_) {
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
+                if (output_binary_) {
+                    outputStateBinary(grid_, state, timer, output_dir_);
+                }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
             }
 
@@ -492,6 +510,9 @@ namespace Opm
             if (output_vtk_) {
                 outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
+            if (output_binary_) {
+                outputStateBinary(grid_, state, timer, output_dir_);
+            }
             outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
             outputWaterCut(watercut, output_dir_);
             if (wells_) {
@@ -532,7 +553,7 @@ namespace Opm
             catch (...) {
                 THROW("Creating directories failed: " << fpath);
             }
-            vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+            vtkfilename << "/output-" << std::setw(5) << std::setfill('0') << step << ".vtu";
             std::ofstream vtkfile(vtkfilename.str().c_str());
             if (!vtkfile) {
                 THROW("Failed to open " << vtkfilename.str());
@@ -573,7 +594,7 @@ namespace Opm
                 catch (...) {
                     THROW("Creating directories failed: " << fpath);
                 }
-                fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+                fname << "/" << std::setw(5) << std::setfill('0') << step << ".txt";
                 std::ofstream file(fname.str().c_str());
                 if (!file) {
                     THROW("Failed to open " << fname.str());
@@ -583,6 +604,28 @@ namespace Opm
             }
         }
 
+        void outputStateBinary(const UnstructuredGrid& grid,
+                               const Opm::PolymerState& state,
+                               const SimulatorTimer& simtimer,
+                               const std::string& output_dir)
+        {
+#ifdef HAVE_ERT
+            Opm::DataMap dm;
+            dm["saturation"] = &state.saturation();
+            dm["pressure"] = &state.pressure();
+            dm["concentration"] = &state.concentration();
+            dm["cmax"] = &state.maxconcentration();
+            std::vector<double> cell_velocity;
+            Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+            dm["velocity"] = &cell_velocity;
+
+            writeECLData(grid, dm, simtimer.currentStepNum(), simtimer.currentTime(), simtimer.currentDateTime(),
+                         output_dir, "polymer_ecl");
+#else
+        THROW("Cannot call outputStateBinary() without ert library support. Reconfigure with --with-ert and recompile.");
+#endif
+        }
+
         void outputWaterCut(const Opm::Watercut& watercut,
                             const std::string& output_dir)
         {

opm/polymer/TransportModelCompressiblePolymer.cpp

@@ -1028,8 +1028,8 @@ namespace Opm
             THROW("In solveMultiCell(), we did not converge after "
                   << num_iters << " iterations. Delta c = " << max_c_change);
         }
-        std::cout << "Solved " << num_cells << " cell multicell problem in "
-                  << num_iters << " iterations." << std::endl;
+        // std::cout << "Solved " << num_cells << " cell multicell problem in "
+        //           << num_iters << " iterations." << std::endl;
     }
 
     void TransportModelCompressiblePolymer::fracFlow(double s, double c, double cmax,