Merge remote-tracking branch 'atgeirr/master'

README

@@ -0,0 +1,73 @@
+Open Porous Media Core Library
+==============================
+
+These are release notes for opm-core.
+
+
+CONTENT
+-------
+
+opm-core is the core library within OPM and contains the following 
+
+* Eclipse deck input and preprosessing
+* Fluid properties (basic PVT models and rock properties)
+* Grid handling (cornerpoint grids, unstructured grid interface)
+* Linear Algebra (interface to different linear solvers)
+* Pressure solvers (various discretization schemes, flow models)
+* Simulators (some basic examples of simulators based on sequential splitting schemes)
+* Transport solvers (various discretization schemes, flow models)
+* Utilities (input and output processing, unit conversion)
+* Wells (basic well handling)
+
+
+LICENSE
+-------
+
+The library is distributed under the GNU General Public License,
+version 3 or later (GPLv3+).
+
+
+PLATFORMS
+---------
+
+The opm-core module is designed to run on Linux platforms. It is also
+regularly run on Mac OS X. No efforts have been made to ensure that
+the code will compile and run on windows platforms.
+
+
+DOWNLOADING
+-----------
+
+git clone git://github.com/OPM/opm-core.git
+
+
+BUILDING
+--------
+
+ cd ../opm-core
+ autoreconf -i
+ ./configure
+ make
+ sudo make install
+
+
+DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS
+-------------------------------------------
+
+(to be updated)
+
+
+DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
+-----------------------------------------
+
+blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf
+git doxygen umfpack
+
+
+DOCUMENTATION
+-------------
+
+Efforts have been made to document the code with Doxygen.
+In order to build the documentation, enter the command
+$ doxygen
+in the topmost directory.

configure.ac

@@ -8,6 +8,10 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects])
 
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 
+# Needed for automake since version 1.12 because extra-portability
+# warnings were then added to -Wall. Ifdef makes it backwards compatible.
+m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
+
 AC_CONFIG_MACRO_DIR([m4])
 AC_CONFIG_SRCDIR([opm/core/grid.h])
 AC_CONFIG_HEADERS([config.h])

examples/sim_2p_comp_reorder.cpp

@@ -277,9 +277,9 @@ main(int argc, char** argv)
     rep.report(std::cout);
 
     if (output) {
-      std::string filename = output_dir + "/walltime.param";
+        std::string filename = output_dir + "/walltime.param";
-      std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+        std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
-      rep.reportParam(tot_os);
+        rep.reportParam(tot_os);
     }
 
 }

opm/core/GridManager.cpp

@@ -166,6 +166,16 @@ namespace Opm
     // Construct tensor grid from deck.
     void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck)
     {
+        // Extract logical cartesian size.
+        std::vector<int> dims;
+        if (deck.hasField("DIMENS")) {
+            dims = deck.getIntegerValue("DIMENS");
+        } else if (deck.hasField("SPECGRID")) {
+            dims = deck.getSPECGRID().dimensions;
+        } else {
+            THROW("Deck must have either DIMENS or SPECGRID.");
+        }
+
         // Extract coordinates (or offsets from top, in case of z).
         const std::vector<double>& dxv = deck.getFloatingPointValue("DXV");
         const std::vector<double>& dyv = deck.getFloatingPointValue("DYV");
@@ -174,6 +184,17 @@ namespace Opm
         std::vector<double> y = coordsFromDeltas(dyv);
         std::vector<double> z = coordsFromDeltas(dzv);
 
+        // Check that number of cells given are consistent with DIMENS/SPECGRID.
+        if (dims[0] != int(dxv.size())) {
+            THROW("Number of DXV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[1] != int(dyv.size())) {
+            THROW("Number of DYV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[2] != int(dzv.size())) {
+            THROW("Number of DZV data points do not match DIMENS or SPECGRID.");
+        }
+
         // Extract top corner depths, if available.
         const double* top_depths = 0;
         std::vector<double> top_depths_vec;

opm/core/fluid/RockFromDeck.cpp

@@ -36,8 +36,6 @@ namespace Opm
         PermeabilityKind fillTensor(const EclipseGridParser&                 parser,
                                     std::vector<const std::vector<double>*>& tensor,
                                     std::tr1::array<int,9>&                     kmap);
-
-        int numGlobalCells(const EclipseGridParser& parser);
     } // anonymous namespace
 
 
@@ -334,26 +332,6 @@ namespace Opm
             return kind;
         }
 
-        int numGlobalCells(const EclipseGridParser& parser)
-        {
-            int ngc = -1;
-
-            if (parser.hasField("DIMENS")) {
-                const std::vector<int>&
-                    dims = parser.getIntegerValue("DIMENS");
-
-                ngc = dims[0] * dims[1] * dims[2];
-            }
-            else if (parser.hasField("SPECGRID")) {
-                const SPECGRID& sgr = parser.getSPECGRID();
-
-                ngc  = sgr.dimensions[ 0 ];
-                ngc *= sgr.dimensions[ 1 ];
-                ngc *= sgr.dimensions[ 2 ];
-            }
-
-            return ngc;
-        }
     } // anonymous namespace
 
 } // namespace Opm

opm/core/pressure/CompressibleTpfa.cpp

@@ -504,7 +504,7 @@ namespace Opm
         cq.phasemobf = &face_phasemob_[0];
         cq.voldiscr = &cell_voldisc_[0];
         int was_adjusted = 0;
-        if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
+        if (! (rock_comp_props_ && rock_comp_props_->isActive())) {
             was_adjusted =
                 cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                       &face_gravcap_[0], cell_press, well_bhp,

opm/core/pressure/CompressibleTpfa.hpp

@@ -61,7 +61,7 @@ namespace Opm
         ///                                   and completions does not change during the
         ///                                   run. However, controls (only) are allowed
         ///                                   to change.
-	CompressibleTpfa(const UnstructuredGrid& grid,
+        CompressibleTpfa(const UnstructuredGrid& grid,
                          const BlackoilPropertiesInterface& props,
                          const RockCompressibility* rock_comp_props,
                          const LinearSolverInterface& linsolver,
@@ -71,8 +71,8 @@ namespace Opm
                          const double* gravity,
                          const Wells* wells);
 
-	/// Destructor.
+        /// Destructor.
-	~CompressibleTpfa();
+        ~CompressibleTpfa();
 
         /// Solve the pressure equation by Newton-Raphson scheme.
         /// May throw an exception if the number of iterations
@@ -111,13 +111,13 @@ namespace Opm
         void solveIncrement();
         double residualNorm() const;
         double incrementNorm() const;
-	void computeResults(BlackoilState& state,
+        void computeResults(BlackoilState& state,
                             WellState& well_state) const;
     protected:
         void computeWellPotentials(const BlackoilState& state);
 
         // ------ Data that will remain unmodified after construction. ------
-	const UnstructuredGrid& grid_;
+        const UnstructuredGrid& grid_;
         const BlackoilPropertiesInterface& props_;
         const RockCompressibility* rock_comp_props_;
         const LinearSolverInterface& linsolver_;
@@ -126,12 +126,12 @@ namespace Opm
         const int maxiter_;
         const double* gravity_; // May be NULL
         const Wells* wells_;    // May be NULL, outside may modify controls (only) between calls to solve().
-	std::vector<double> htrans_;
+        std::vector<double> htrans_;
-	std::vector<double> trans_ ;
+        std::vector<double> trans_ ;
         std::vector<int> allcells_;
 
         // ------ Internal data for the cfs_tpfa_res solver. ------
-	struct cfs_tpfa_res_data* h_;
+        struct cfs_tpfa_res_data* h_;
 
         // ------ Data that will be modified for every solve. ------
         std::vector<double> wellperf_gpot_;

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -77,7 +77,6 @@ namespace Opm
 
     private:
         // Data.
-
         // Parameters for output.
         bool output_;
         bool output_vtk_;
@@ -133,7 +132,38 @@ namespace Opm
         return pimpl_->run(timer, state, well_state);
     }
 
-
+    static void reportVolumes(std::ostream &os, double satvol[2], double tot_porevol_init,
+                              double tot_injected[2], double tot_produced[2],
+                              double injected[2], double produced[2],
+                              double init_satvol[2])
+    {
+        std::cout.precision(5);
+        const int width = 18;
+        os << "\nVolume balance report (all numbers relative to total pore volume).\n";
+        os << "    Saturated volumes:     "
+           << std::setw(width) << satvol[0]/tot_porevol_init
+           << std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
+        os << "    Injected volumes:      "
+           << std::setw(width) << injected[0]/tot_porevol_init
+           << std::setw(width) << injected[1]/tot_porevol_init << std::endl;
+        os << "    Produced volumes:      "
+           << std::setw(width) << produced[0]/tot_porevol_init
+           << std::setw(width) << produced[1]/tot_porevol_init << std::endl;
+        os << "    Total inj volumes:     "
+           << std::setw(width) << tot_injected[0]/tot_porevol_init
+           << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
+        os << "    Total prod volumes:    "
+           << std::setw(width) << tot_produced[0]/tot_porevol_init
+           << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
+        os << "    In-place + prod - inj: "
+           << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
+           << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
+        os << "    Init - now - pr + inj: "
+           << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
+           << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
+           << std::endl;
+        os.precision(8);
+    }
 
     static void outputStateVtk(const UnstructuredGrid& grid,
                                const Opm::TwophaseState& state,
@@ -145,10 +175,10 @@ namespace Opm
         vtkfilename << output_dir << "/vtk_files";
         boost::filesystem::path fpath(vtkfilename.str());
         try {
-          create_directories(fpath);
+            create_directories(fpath);
         }
         catch (...) {
-          THROW("Creating directories failed: " << fpath);
+            THROW("Creating directories failed: " << fpath);
         }
         vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
         std::ofstream vtkfile(vtkfilename.str().c_str());
@@ -164,6 +194,27 @@ namespace Opm
         Opm::writeVtkData(grid, dm, vtkfile);
     }
 
+    static void outputVectorMatlab(const std::string& name,
+                                   const std::vector<int>& vec,
+                                   const int step,
+                                   const std::string& output_dir)
+    {
+        std::ostringstream fname;
+        fname << output_dir << "/" << name;
+        boost::filesystem::path fpath = fname.str();
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            THROW("Creating directories failed: " << fpath);
+        }
+        fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+        std::ofstream file(fname.str().c_str());
+        if (!file) {
+            THROW("Failed to open " << fname.str());
+        }
+        std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(file, "\n"));
+    }
 
     static void outputStateMatlab(const UnstructuredGrid& grid,
                                   const Opm::TwophaseState& state,
@@ -183,10 +234,10 @@ namespace Opm
             fname << output_dir << "/" << it->first;
             boost::filesystem::path fpath = fname.str();
             try {
-              create_directories(fpath);
+                create_directories(fpath);
             }
             catch (...) {
-              THROW("Creating directories failed: " << fpath);
+                THROW("Creating directories failed: " << fpath);
             }
             fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
             std::ofstream file(fname.str().c_str());
@@ -377,6 +428,9 @@ namespace Opm
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+                outputVectorMatlab(std::string("reorder_it"),
+                                   tsolver_.getReorderIterations(),
+                                   timer.currentStepNum(), output_dir_);
             }
 
             SimulatorReport sreport;
@@ -474,6 +528,14 @@ namespace Opm
                 if (use_segregation_split_) {
                     tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
                 }
+                watercut.push(timer.currentTime() + timer.currentStepLength(),
+                              produced[0]/(produced[0] + produced[1]),
+                              tot_produced[0]/tot_porevol_init);
+                if (wells_) {
+                    wellreport.push(props_, *wells_, state.saturation(),
+                                    timer.currentTime() + timer.currentStepLength(),
+                                    well_state.bhp(), well_state.perfRates());
+                }
             }
             transport_timer.stop();
             double tt = transport_timer.secsSinceStart();
@@ -486,41 +548,10 @@ namespace Opm
             tot_injected[1] += injected[1];
             tot_produced[0] += produced[0];
             tot_produced[1] += produced[1];
-            std::cout.precision(5);
+            reportVolumes(std::cout,satvol, tot_porevol_init,
-            const int width = 18;
+                          tot_injected, tot_produced,
-            std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n";
+                          injected, produced,
-            std::cout << "    Saturated volumes:     "
+                          init_satvol);
-                      << std::setw(width) << satvol[0]/tot_porevol_init
-                      << std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
-            std::cout << "    Injected volumes:      "
-                      << std::setw(width) << injected[0]/tot_porevol_init
-                      << std::setw(width) << injected[1]/tot_porevol_init << std::endl;
-            std::cout << "    Produced volumes:      "
-                      << std::setw(width) << produced[0]/tot_porevol_init
-                      << std::setw(width) << produced[1]/tot_porevol_init << std::endl;
-            std::cout << "    Total inj volumes:     "
-                      << std::setw(width) << tot_injected[0]/tot_porevol_init
-                      << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
-            std::cout << "    Total prod volumes:    "
-                      << std::setw(width) << tot_produced[0]/tot_porevol_init
-                      << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
-            std::cout << "    In-place + prod - inj: "
-                      << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
-                      << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
-            std::cout << "    Init - now - pr + inj: "
-                      << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
-                      << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
-                      << std::endl;
-            std::cout.precision(8);
-
-            watercut.push(timer.currentTime() + timer.currentStepLength(),
-                          produced[0]/(produced[0] + produced[1]),
-                          tot_produced[0]/tot_porevol_init);
-            if (wells_) {
-                wellreport.push(props_, *wells_, state.saturation(),
-                                timer.currentTime() + timer.currentStepLength(),
-                                well_state.bhp(), well_state.perfRates());
-            }
             sreport.total_time =  step_timer.secsSinceStart();
             if (output_) {
                 sreport.reportParam(tstep_os);
@@ -532,6 +563,9 @@ namespace Opm
                 outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
             outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            outputVectorMatlab(std::string("reorder_it"),
+                               tsolver_.getReorderIterations(),
+                               timer.currentStepNum(), output_dir_);
             outputWaterCut(watercut, output_dir_);
             if (wells_) {
                 outputWellReport(wellreport, output_dir_);

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -96,7 +96,7 @@ namespace Opm
         props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
         props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
 
-        // Check non-miscibility requirement (only done for first cell).
+        // Check immiscibility requirement (only done for first cell).
         if (A_[1] != 0.0 || A_[2] != 0.0) {
             THROW("TransportModelCompressibleTwophase requires a property object without miscibility.");
         }

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -53,6 +53,7 @@ namespace Opm
           dt_(0.0),
           saturation_(grid.number_of_cells, -1.0),
           fractionalflow_(grid.number_of_cells, -1.0),
+          reorder_iterations_(grid.number_of_cells, 0),
           mob_(2*grid.number_of_cells, -1.0)
 #ifdef EXPERIMENT_GAUSS_SEIDEL
         , ia_upw_(grid.number_of_cells + 1, -1),
@@ -101,11 +102,18 @@ namespace Opm
                                &seq[0], &comp[0], &ncomp,
                                &ia_downw_[0], &ja_downw_[0]);
 #endif
-
+        std::fill(reorder_iterations_.begin(),reorder_iterations_.end(),0);
         reorderAndTransport(grid_, darcyflux);
         toBothSat(saturation_, saturation);
     }
 
+
+    const std::vector<int>& TransportModelTwophase::getReorderIterations() const
+    {
+        return reorder_iterations_;
+    }
+
+
     // Residual function r(s) for a single-cell implicit Euler transport
     //
     //     r(s) = s - s0 + dt/pv*( influx + outflux*f(s) )
@@ -170,9 +178,11 @@ namespace Opm
         // if (std::fabs(r0) < tol_) {
         //     return;
         // }
-        int iters_used;
+        int iters_used = 0;
         // saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
         saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
+        // add if it is iteration on an out loop
+        reorder_iterations_[cell] = reorder_iterations_[cell] + iters_used;
         fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
     }
 
@@ -544,8 +554,9 @@ namespace Opm
         const int cell = cells[pos];
         GravityResidual res(*this, cells, pos, gravflux);
         if (std::fabs(res(saturation_[cell])) > tol_) {
-            int iters_used;
+            int iters_used = 0;
             saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
+            reorder_iterations_[cell] = reorder_iterations_[cell] + iters_used;
         }
         saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]);
         mobility(saturation_[cell], cell, &mob_[2*cell]);

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -23,7 +23,7 @@
 #include <opm/core/transport/reorder/TransportModelInterface.hpp>
 #include <vector>
 #include <map>
-
+#include <ostream>
 struct UnstructuredGrid;
 
 namespace Opm
@@ -75,6 +75,10 @@ namespace Opm
                           const double dt,
                           std::vector<double>& saturation);
 
+        //// Return the number of iterations used by the reordering solver.
+        //// \return vector of iteration per cell
+        const std::vector<int>& getReorderIterations() const;
+
     private:
         virtual void solveSingleCell(const int cell);
         virtual void solveMultiCell(const int num_cells, const int* cells);
@@ -83,7 +87,6 @@ namespace Opm
                                     const int pos,
                                     const double* gravflux);
         int solveGravityColumn(const std::vector<int>& cells);
-
     private:
         const UnstructuredGrid& grid_;
         const IncompPropertiesInterface& props_;
@@ -99,6 +102,8 @@ namespace Opm
         double dt_;
         std::vector<double> saturation_;        // one per cell, only water saturation!
         std::vector<double> fractionalflow_;  // = m[0]/(m[0] + m[1]) per cell
+        std::vector<int> reorder_iterations_;
+        //std::vector<double> reorder_fval_;
         // For gravity segregation.
         std::vector<double> gravflux_;
         std::vector<double> mob_;

opm/core/utility/miscUtilitiesBlackoil.cpp

@@ -243,9 +243,9 @@ namespace Opm
         // matrix data.
         std::fill(surfacevol, surfacevol + n*np, 0.0);
         for (int i = 0; i < n; ++i) {
-            for (int row = 0; row < np; ++row) {
+            for (int col = 0; col < np; ++col) {
-                for (int col = 0; col < np; ++col) {
+                for (int row = 0; row < np; ++row) {
-                    surfacevol[i*np + row] += A[i*np*np + row*np + col] * saturation[i*np + col];
+                    surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col];
                 }
             }
         }