Merge pull request #102 from atgeirr/separate-wellstate

Specialized well state class for the fully implicit blackoil solver.

CMakeLists_files.cmake

@@ -106,4 +106,5 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp
+	opm/autodiff/WellStateFullyImplicitBlackoil.hpp
 	)

examples/sim_fibo_ad.cpp

@@ -39,7 +39,7 @@
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
@@ -205,7 +205,7 @@ try
     std::cout << "\n\n================    Starting main simulation loop     ===============\n"
               << std::flush;
 
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
     SimulatorTimer simtimer;
 
@@ -248,7 +248,7 @@ try
 
     SimulatorReport rep;
     // With a deck, we may have more epochs etc.
-    WellState well_state;
+    WellStateFullyImplicitBlackoil well_state;
     int step = 0;
     SimulatorTimer simtimer;
     // Use timer for last epoch to obtain total time.

examples/test_implicit_ad.cpp

@@ -23,6 +23,7 @@
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAd.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -37,7 +38,6 @@
 #include <opm/core/utility/Units.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/initState.hpp>
 
 #include <boost/shared_ptr.hpp>
@@ -110,7 +110,7 @@ try
     initStateBasic(*g, props0, param, 0.0, state);
     initBlackoilSurfvol(*g, props0, state);
 
-    Opm::WellState well_state;
+    Opm::WellStateFullyImplicitBlackoil well_state;
     well_state.init(wells.get(), state);
 
     solver.step(1.0, state, well_state);

opm/autodiff/FullyImplicitBlackoilSolver.cpp

@@ -24,12 +24,12 @@
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/well_controls.h>
 
@@ -223,7 +223,7 @@ namespace {
     FullyImplicitBlackoilSolver::
     step(const double   dt,
          BlackoilState& x ,
-         WellState&     xw)
+         WellStateFullyImplicitBlackoil& xw)
     {
         const V pvdt = geo_.poreVolume() / dt;
 
@@ -331,7 +331,7 @@ namespace {
 
     FullyImplicitBlackoilSolver::SolutionState
     FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                               const WellStateFullyImplicitBlackoil& xw)
     {
         const int nc = grid_.number_of_cells;
         const int np = x.numPhases();
@@ -427,7 +427,7 @@ namespace {
 
     FullyImplicitBlackoilSolver::SolutionState
     FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
-                                               const WellState&     xw)
+                                               const WellStateFullyImplicitBlackoil& xw)
     {
         const int nc = grid_.number_of_cells;
         const int np = x.numPhases();
@@ -616,7 +616,7 @@ namespace {
     FullyImplicitBlackoilSolver::
     assemble(const V&             pvdt,
              const BlackoilState& x   ,
-             const WellState&     xw  )
+             const WellStateFullyImplicitBlackoil& xw  )
     {
         // Create the primary variables.
         const SolutionState state = variableState(x, xw);
@@ -878,7 +878,7 @@ namespace {
 
     void FullyImplicitBlackoilSolver::updateState(const V& dx,
                                                   BlackoilState& state,
-                                                  WellState& well_state)
+                                                  WellStateFullyImplicitBlackoil& well_state)
     {
         const int np = fluid_.numPhases();
         const int nc = grid_.number_of_cells;

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -33,7 +33,7 @@ namespace Opm {
     class RockCompressibility;
     class LinearSolverInterface;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
 
 
     /// A fully implicit solver for the black-oil problem.
@@ -76,7 +76,7 @@ namespace Opm {
         void
         step(const double   dt    ,
              BlackoilState& state ,
-             WellState&     wstate);
+             WellStateFullyImplicitBlackoil&     wstate);
 
     private:
         // Types and enums
@@ -148,11 +148,11 @@ namespace Opm {
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         SolutionState
         variableState(const BlackoilState& x,
-                      const WellState&     xw);
+                      const WellStateFullyImplicitBlackoil& xw);
 
         void
         computeAccum(const SolutionState& state,
@@ -161,13 +161,13 @@ namespace Opm {
         void
         assemble(const V&             dtpv,
                  const BlackoilState& x,
-                 const WellState&     xw  );
+                 const WellStateFullyImplicitBlackoil& xw);
 
         V solveJacobianSystem() const;
 
         void updateState(const V& dx,
                          BlackoilState& state,
-                         WellState& well_state);
+                         WellStateFullyImplicitBlackoil& well_state);
 
         std::vector<ADB>
         computePressures(const SolutionState& state) const;

opm/autodiff/SimulatorFullyImplicitBlackoil.cpp

@@ -29,6 +29,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -48,7 +49,6 @@
 
 #include <opm/core/grid/ColumnExtract.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
 #include <boost/filesystem.hpp>
@@ -77,7 +77,7 @@ namespace Opm
 
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         // Data.
@@ -127,7 +127,7 @@ namespace Opm
 
     SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer,
                                                         BlackoilState& state,
-                                                        WellState& well_state)
+                                                        WellStateFullyImplicitBlackoil& well_state)
     {
         return pimpl_->run(timer, state, well_state);
     }
@@ -198,7 +198,7 @@ namespace Opm
             std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
         }
     }
-    static void outputWellStateMatlab(const Opm::WellState& well_state,
+    static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
                                   const int step,
                                   const std::string& output_dir)
     {
@@ -312,10 +312,10 @@ namespace Opm
 
     SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
                                                               BlackoilState& state,
-                                                              WellState& well_state)
+                                                              WellStateFullyImplicitBlackoil& well_state)
     {
-        eclipseWriter_.writeInit(timer, state, well_state);
+        eclipseWriter_.writeInit(timer, state, well_state.basicWellState());
-        eclipseWriter_.writeTimeStep(timer, state, well_state);
+        eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
 
         // Initialisation.
         std::vector<double> porevol;
@@ -384,13 +384,9 @@ namespace Opm
 
                 // Optionally, check if well controls are satisfied.
                 if (check_well_controls_) {
-                    Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                      well_state.perfRates(),
-                                                      fractional_flows,
-                                                      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_passed = wells_manager_.conditionsMet(well_state.bhp(), well_state.wellRates(), well_state.wellRates());
                     ++well_control_iteration;
                     if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
                         OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
@@ -426,7 +422,7 @@ namespace Opm
 
             // write an output file for later inspection
             if (output_) {
-                eclipseWriter_.writeTimeStep(timer, state, well_state);
+                eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
             }
 
             // advance to next timestep before reporting at this location

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -37,7 +37,7 @@ namespace Opm
     class LinearSolverInterface;
     class SimulatorTimer;
     class BlackoilState;
-    class WellState;
+    class WellStateFullyImplicitBlackoil;
     struct SimulatorReport;
 
     /// Class collecting all necessary components for a two-phase simulation.
@@ -84,7 +84,7 @@ namespace Opm
         /// \return                    simulation report, with timing data
         SimulatorReport run(SimulatorTimer& timer,
                             BlackoilState& state,
-                            WellState& well_state);
+                            WellStateFullyImplicitBlackoil& well_state);
 
     private:
         class Impl;

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -0,0 +1,101 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+
+
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
+#include <opm/core/simulator/WellState.hpp>
+#include <vector>
+#include <cassert>
+
+namespace Opm
+{
+
+    /// The state of a set of wells, tailored for use by the fully
+    /// implicit blackoil simulator.
+    class WellStateFullyImplicitBlackoil
+    {
+    public:
+        /// Allocate and initialize if wells is non-null.  Also tries
+        /// to give useful initial values to the bhp(), wellRates()
+        /// and perfPhaseRates() fields, depending on controls
+        template <class State>
+        void init(const Wells* wells, const State& state)
+        {
+            if (wells == 0) {
+                return;
+            }
+
+            // We use the WellState::init() function to do bhp and well rates init.
+            // The alternative would be to copy that function wholesale.
+            basic_well_state_.init(wells, state);
+
+            // Initialize perfphaserates_, which must be done here.
+            const int nw = wells->number_of_wells;
+            const int np = wells->number_of_phases;
+            const int nperf = wells->well_connpos[nw];
+            perfphaserates_.resize(nperf * np, 0.0);
+            for (int w = 0; w < nw; ++w) {
+                assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
+                const WellControls* ctrl = wells->ctrls[w];
+                if (well_controls_well_is_shut(ctrl)) {
+                    // Shut well: perfphaserates_ are all zero.
+                } else {
+                    // Open well: Initialize perfphaserates_ to well
+                    // rates divided by the number of perforations.
+                    const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
+                    for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                        for (int p = 0; p < np; ++p) {
+                            perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
+                        }
+                    }
+                }
+            }
+        }
+
+        /// One bhp pressure per well.
+        std::vector<double>& bhp() { return basic_well_state_.bhp(); }
+        const std::vector<double>& bhp() const { return basic_well_state_.bhp(); }
+
+        /// One rate per well and phase.
+        std::vector<double>& wellRates() { return basic_well_state_.wellRates(); }
+        const std::vector<double>& wellRates() const { return basic_well_state_.wellRates(); }
+
+        /// One rate per phase and well connection.
+        std::vector<double>& perfPhaseRates() { return perfphaserates_; }
+        const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
+
+        /// For interfacing with functions that take a WellState.
+        const WellState& basicWellState() const
+        {
+            return basic_well_state_;
+        }
+
+    private:
+        WellState basic_well_state_;
+        std::vector<double> perfphaserates_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED