Merge remote-tracking branch 'upstream/master'

Conflicts:
	opm/core/simulator/SimulatorIncompTwophase.cpp

Makefile.am

@@ -92,9 +92,9 @@ opm/core/pressure/tpfa/compr_source.c			\
 opm/core/pressure/tpfa/ifs_tpfa.c			\
 opm/core/pressure/tpfa/trans_tpfa.c			\
 opm/core/pressure/well.c				\
+opm/core/simulator/SimulatorIncompTwophase.cpp		\
 opm/core/simulator/SimulatorReport.cpp			\
 opm/core/simulator/SimulatorTimer.cpp			\
-opm/core/simulator/SimulatorTwophase.cpp		\
 opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp	\
 opm/core/transport/reorder/TransportModelInterface.cpp	\
 opm/core/transport/reorder/TransportModelTwophase.cpp	\
@@ -196,8 +196,8 @@ opm/core/pressure/tpfa/ifs_tpfa.h			\
 opm/core/pressure/tpfa/trans_tpfa.h			\
 opm/core/simulator/BlackoilState.hpp			\
 opm/core/simulator/SimulatorReport.hpp			\
+opm/core/simulator/SimulatorIncompTwophase.hpp		\
 opm/core/simulator/SimulatorTimer.hpp			\
-opm/core/simulator/SimulatorTwophase.hpp		\
 opm/core/simulator/TwophaseState.hpp			\
 opm/core/simulator/WellState.hpp			\
 opm/core/transport/CSRMatrixBlockAssembler.hpp		\

README

@@ -0,0 +1,41 @@
+These are the release notes for opm-core
+
+Open Porous Media Core Library
+
+CONTENT
+opm-core is the core library within OPM and contains the following 
+
+* Eclipse preprosessing
+* Fluid properties (basic PVT models and rock properties)
+* Grid (generates different grids)
+* Linear Algerbra (interface to different linear solvers)
+* Pressure solvers (different discretization schemes, different flow models)
+* Simulator (some basic examples of simulators based on sequential schemes)
+* Transport solvers (different discretization schemes, different flow models)
+* Utility (input and output processing, unit conversion)
+* Wells (basic well handling)
+
+ON WHAT PLATFORMS DOES IT RUN?
+The opm-core module is designed to run on linux platforms. No efforts have
+been made to ensure that the code will compile and run on windows platforms.
+
+DOCUMENTATION
+Efforts have been made to document the code with Doxygen.
+
+DOWNLOAD opm-core
+git clone git://github.com/OPM/opm-core.git
+
+BUILDING opm-core
+cd ../opm-core
+  autoreconf -i
+ ./configure
+ make
+ sudo make install
+
+
+DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS
+
+
+DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
+blas libblas3 lapack liblapack3 libboost libxml2 gcc automake autoconf git
+doxygen umfpack

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_incomp_reorder.cpp

@@ -43,7 +43,7 @@
 
 #include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/simulator/WellState.hpp>
-#include <opm/core/simulator/SimulatorTwophase.hpp>
+#include <opm/core/simulator/SimulatorIncompTwophase.hpp>
 
 #include <boost/scoped_ptr.hpp>
 #include <boost/filesystem.hpp>
@@ -188,7 +188,6 @@ main(int argc, char** argv)
         filename = output_dir + "/step_timing.param";
         std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
         step_os.close();
-      
         param.writeParam(output_dir + "/simulation.param");
     }
 
@@ -200,11 +199,12 @@ main(int argc, char** argv)
     SimulatorReport rep;
     if (!use_deck) {
         // Simple simulation without a deck.
-        SimulatorTwophase simulator(param,
+        WellsManager wells; // no wells.
+        SimulatorIncompTwophase simulator(param,
                                           *grid->c_grid(),
                                           *props,
                                           rock_comp->isActive() ? rock_comp.get() : 0,
-                                    0, // wells
+                                          wells,
                                           src,
                                           bcs.c_bcs(),
                                           linsolver,
@@ -255,11 +255,11 @@ main(int argc, char** argv)
             }
 
             // Create and run simulator.
-            SimulatorTwophase simulator(param,
+            SimulatorIncompTwophase simulator(param,
                                               *grid->c_grid(),
                                               *props,
                                               rock_comp->isActive() ? rock_comp.get() : 0,
-                                        wells.c_wells(),
+                                              wells,
                                               src,
                                               bcs.c_bcs(),
                                               linsolver,
@@ -277,7 +277,6 @@ main(int argc, char** argv)
         }
     }
 
-    epoch_os.close();
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     rep.report(std::cout);
 
@@ -285,7 +284,6 @@ main(int argc, char** argv)
       std::string filename = output_dir + "/walltime.param";
       std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
       rep.reportParam(tot_os);
-      tot_os.close();
     }
 
 }

opm/core/grid/cpgpreprocess/preprocess.c

@@ -587,8 +587,8 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
                     z1 = sign*zcorn[i+2*nx*(j+2*ny*(k))];
                     z2 = sign*zcorn[i+2*nx*(j+2*ny*(k+1))];
 
-                    c1 = i/2 + nx*(j/2 + ny*k/2);
+                    c1 = i/2 + nx*(j/2 + ny*(k/2));
-                    c2 = i/2 + nx*(j/2 + ny*(k+1)/2);
+                    c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
 
                     if (((actnum == NULL) ||
                          (actnum[c1] && actnum[c2]))

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -21,7 +21,7 @@
 #include "config.h"
 #endif // HAVE_CONFIG_H
 
-#include <opm/core/simulator/SimulatorTwophase.hpp>
+#include <opm/core/simulator/SimulatorIncompTwophase.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 
@@ -37,6 +37,8 @@
 #include <opm/core/utility/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 
+#include <opm/core/wells/WellsManager.hpp>
+
 #include <opm/core/fluid/IncompPropertiesInterface.hpp>
 #include <opm/core/fluid/RockCompressibility.hpp>
 
@@ -56,14 +58,14 @@
 namespace Opm
 {
 
-    class SimulatorTwophase::Impl
+    class SimulatorIncompTwophase::Impl
     {
     public:
         Impl(const parameter::ParameterGroup& param,
              const UnstructuredGrid& grid,
              const IncompPropertiesInterface& props,
              const RockCompressibility* rock_comp,
-             const Wells* wells,
+             WellsManager& wells_manager,
              const std::vector<double>& src,
              const FlowBoundaryConditions* bcs,
              LinearSolverInterface& linsolver,
@@ -80,6 +82,9 @@ namespace Opm
         bool output_vtk_;
         std::string output_dir_;
         int output_interval_;
+        // Parameters for well control
+        bool check_well_controls_;
+        int max_well_control_iterations_;
         // Parameters for transport solver.
         int num_transport_substeps_;
         bool use_segregation_split_;
@@ -87,11 +92,10 @@ namespace Opm
         const UnstructuredGrid& grid_;
         const IncompPropertiesInterface& props_;
         const RockCompressibility* rock_comp_;
+        WellsManager& wells_manager_;
         const Wells* wells_;
         const std::vector<double>& src_;
-        //const FlowBoundaryConditions* bcs_;
+        const FlowBoundaryConditions* bcs_;
-        //const LinearSolverInterface& linsolver_;
-        //const double* gravity_;
         // Solvers
         IncompTpfa psolver_;
         TransportModelTwophase tsolver_;
@@ -104,24 +108,24 @@ namespace Opm
 
 
 
-    SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param,
+    SimulatorIncompTwophase::SimulatorIncompTwophase(const parameter::ParameterGroup& param,
                                                      const UnstructuredGrid& grid,
                                                      const IncompPropertiesInterface& props,
                                                      const RockCompressibility* rock_comp,
-                                         const Wells* wells,
+                                                     WellsManager& wells_manager,
                                                      const std::vector<double>& src,
                                                      const FlowBoundaryConditions* bcs,
                                                      LinearSolverInterface& linsolver,
                                                      const double* gravity)
     {
-        pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity));
+        pimpl_.reset(new Impl(param, grid, props, rock_comp, wells_manager, src, bcs, linsolver, gravity));
     }
 
 
 
 
 
-    SimulatorReport SimulatorTwophase::run(SimulatorTimer& timer,
+    SimulatorReport SimulatorIncompTwophase::run(SimulatorTimer& timer,
                                                  TwophaseState& state,
                                                  WellState& well_state)
     {
@@ -293,11 +297,43 @@ namespace Opm
     }
 
 
-    SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param,
+    static bool allNeumannBCs(const FlowBoundaryConditions* bcs)
+    {
+        if (bcs == NULL) {
+            return true;
+        } else {
+            return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE)
+                == bcs->type + bcs->nbc;
+        }
+    }
+
+
+    static bool allRateWells(const Wells* wells)
+    {
+        if (wells == NULL) {
+            return true;
+        }
+        const int nw = wells->number_of_wells;
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells->ctrls[w];
+            if (wc->current >= 0) {
+                if (wc->type[wc->current] == BHP) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+
+
+
+
+    SimulatorIncompTwophase::Impl::Impl(const parameter::ParameterGroup& param,
                                         const UnstructuredGrid& grid,
                                         const IncompPropertiesInterface& props,
                                         const RockCompressibility* rock_comp,
-                                  const Wells* wells,
+                                        WellsManager& wells_manager,
                                         const std::vector<double>& src,
                                         const FlowBoundaryConditions* bcs,
                                         LinearSolverInterface& linsolver,
@@ -305,16 +341,15 @@ namespace Opm
         : grid_(grid),
           props_(props),
           rock_comp_(rock_comp),
-          wells_(wells),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
           src_(src),
-          //bcs_(bcs),
+          bcs_(bcs),
-          //linsolver_(linsolver),
-          //gravity_(gravity),
           psolver_(grid, props, rock_comp, linsolver,
                    param.getDefault("nl_pressure_residual_tolerance", 0.0),
                    param.getDefault("nl_pressure_change_tolerance", 1.0),
                    param.getDefault("nl_pressure_maxiter", 10),
-                   gravity, wells, src, bcs),
+                   gravity, wells_manager.c_wells(), src, bcs),
           tsolver_(grid, props,
                    param.getDefault("nl_tolerance", 1e-9),
                    param.getDefault("nl_maxiter", 30))
@@ -335,6 +370,10 @@ namespace Opm
             output_interval_ = param.getDefault("output_interval", 1);
         }
 
+        // Well control related init.
+        check_well_controls_ = param.getDefault("check_well_controls", false);
+        max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
+
         // Transport related init.
         num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
         use_segregation_split_ = param.getDefault("use_segregation_split", false);
@@ -354,7 +393,7 @@ namespace Opm
 
 
 
-    SimulatorReport SimulatorTwophase::Impl::run(SimulatorTimer& timer,
+    SimulatorReport SimulatorIncompTwophase::Impl::run(SimulatorTimer& timer,
                                                        TwophaseState& state,
                                                        WellState& well_state)
     {
@@ -414,18 +453,77 @@ namespace Opm
                                    tsolver_.getReorderIterations(),
                                    timer.currentStepNum(), output_dir_);
             }
+
             SimulatorReport sreport;
 
-            // Solve pressure.
+            // Solve pressure equation.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, allcells_, state.saturation(), fractional_flows);
+                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
+            }
+            bool well_control_passed = !check_well_controls_;
+            int well_control_iteration = 0;
             do {
+                // Run solver.
                 pressure_timer.start();
+                std::vector<double> initial_pressure = state.pressure();
                 psolver_.solve(timer.currentStepLength(), state, well_state);
+
+                // Renormalize pressure if rock is incompressible, and
+                // there are no pressure conditions (bcs or wells).
+                // It is deemed sufficient for now to renormalize
+                // using geometric volume instead of pore volume.
+                if ((rock_comp_ == NULL || !rock_comp_->isActive())
+                    && allNeumannBCs(bcs_) && allRateWells(wells_)) {
+                    // Compute average pressures of previous and last
+                    // step, and total volume.
+                    double av_prev_press = 0.0;
+                    double av_press = 0.0;
+                    double tot_vol = 0.0;
+                    const int num_cells = grid_.number_of_cells;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                        tot_vol       += grid_.cell_volumes[cell];
+                    }
+                    // Renormalization constant
+                    const double ren_const = (av_prev_press - av_press)/tot_vol;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        state.pressure()[cell] += ren_const;
+                    }
+                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                    for (int well = 0; well < num_wells; ++well) {
+                        well_state.bhp()[well] += ren_const;
+                    }
+                }
+
+                // Stop timer and report.
                 pressure_timer.stop();
                 double pt = pressure_timer.secsSinceStart();
                 std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
                 ptime += pt;
                 sreport.pressure_time = pt;
-            } while (false);
+
+                // Optionally, check if well controls are satisfied.
+                if (check_well_controls_) {
+                    Opm::computePhaseFlowRatesPerWell(*wells_,
+                                                      fractional_flows,
+                                                      well_state.perfRates(),
+                                                      well_resflows_phase);
+                    std::cout << "Checking well conditions." << std::endl;
+                    // For testing we set surface := reservoir
+                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
+                    ++well_control_iteration;
+                    if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
+                        THROW("Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
+                    }
+                    if (!well_control_passed) {
+                        std::cout << "Well controls not passed, solving again." << std::endl;
+                    } else {
+                        std::cout << "Well conditions met." << std::endl;
+                    }
+                }
+            } while (!well_control_passed);
 
             // Update pore volumes if rock is compressible.
             if (rock_comp_ && rock_comp_->isActive()) {
@@ -479,8 +577,6 @@ namespace Opm
             if (output_) {
                 sreport.reportParam(tstep_os);
             }
-       
-            
         }
 
         if (output_) {

opm/core/simulator/SimulatorIncompTwophase.hpp

@@ -17,8 +17,8 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED
+#ifndef OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED
-#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED
+#define OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED
 
 #include <boost/shared_ptr.hpp>
 #include <vector>
@@ -32,6 +32,7 @@ namespace Opm
     namespace parameter { class ParameterGroup; }
     class IncompPropertiesInterface;
     class RockCompressibility;
+    class WellsManager;
     class LinearSolverInterface;
     class SimulatorTimer;
     class TwophaseState;
@@ -39,7 +40,7 @@ namespace Opm
     struct SimulatorReport;
 
     /// Class collecting all necessary components for a two-phase simulation.
-    class SimulatorTwophase
+    class SimulatorIncompTwophase
     {
     public:
         /// Initialise from parameters and objects to observe.
@@ -61,16 +62,16 @@ namespace Opm
         /// \param[in] grid          grid data structure
         /// \param[in] props         fluid and rock properties
         /// \param[in] rock_comp     if non-null, rock compressibility properties
-        /// \param[in] wells       if non-null, wells data structure
+        /// \param[in] well_manager  well manager, may manage no (null) wells
         /// \param[in] src           source terms
         /// \param[in] bcs           boundary conditions, treat as all noflow if null
         /// \param[in] linsolver     linear solver
         /// \param[in] gravity       if non-null, gravity vector
-       SimulatorTwophase(const parameter::ParameterGroup& param,
+       SimulatorIncompTwophase(const parameter::ParameterGroup& param,
                                const UnstructuredGrid& grid,
                                const IncompPropertiesInterface& props,
                                const RockCompressibility* rock_comp,
-                         const Wells* wells,
+                               WellsManager& wells_manager,
                                const std::vector<double>& src,
                                const FlowBoundaryConditions* bcs,
                                LinearSolverInterface& linsolver,