Simulators handle WPOLYMER keyword.

They do this by using the class PolymerInflowFromDeck
instead of PolymerInflowBasic if there is a WPOLYMER
keyword somewhere in the deck epochs. If there is no
WPOLYMER, the parameters 'poly_start_days' etc. will
be used to construct an instance of PolymerInflowBasic
instead.

examples/sim_poly2p_comp_reorder.cpp

@@ -44,6 +44,7 @@
 #include <opm/polymer/PolymerBlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/polymer/SimulatorCompressiblePolymer.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
 #include <opm/polymer/PolymerProperties.hpp>
 
 #include <boost/scoped_ptr.hpp>
@@ -237,12 +238,16 @@ main(int argc, char** argv)
     SimulatorReport rep;
     if (!use_deck) {
         // Simple simulation without a deck.
+        PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                          param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                          param.getDefault("poly_amount", poly_props.cMax()));
         SimulatorCompressiblePolymer simulator(param,
                                                *grid->c_grid(),
                                                *props,
                                                poly_props,
                                                rock_comp->isActive() ? rock_comp.get() : 0,
                                                0, // wells
+                                               polymer_inflow,
                                                src,
                                                bcs.c_bcs(),
                                                linsolver,
@@ -262,6 +267,15 @@ main(int argc, char** argv)
         deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
         simtimer.init(*deck);
         const double total_time = simtimer.totalTime();
+        // Check for WPOLYMER presence in last epoch to decide
+        // polymer injection control type.
+        const bool use_wpolymer = deck->hasField("WPOLYMER");
+        if (use_wpolymer) {
+            if (param.has("poly_start_days")) {
+                MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
+                        "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
+            }
+        }
         for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
             // Set epoch index.
             deck->setCurrentEpoch(epoch);
@@ -283,8 +297,20 @@ main(int argc, char** argv)
                       << "\n                  (number of steps: "
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
-            // Create new wells, well_state
+            // Create new wells, polymer inflow controls.
             WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
+            if (use_wpolymer) {
+                if (wells.c_wells() == 0) {
+                    THROW("Cannot control polymer injection via WPOLYMER without wells.");
+                }
+                polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells()));
+            } else {
+                polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                                            param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                                            param.getDefault("poly_amount", poly_props.cMax())));
+            }
+
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.
@@ -299,6 +325,7 @@ main(int argc, char** argv)
                                                    poly_props,
                                                    rock_comp->isActive() ? rock_comp.get() : 0,
                                                    wells.c_wells(),
+                                                   *polymer_inflow,
                                                    src,
                                                    bcs.c_bcs(),
                                                    linsolver,

examples/sim_poly2p_incomp_reorder.cpp

@@ -44,6 +44,7 @@
 #include <opm/polymer/PolymerState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/polymer/SimulatorPolymer.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
 #include <opm/polymer/PolymerProperties.hpp>
 
 #include <boost/scoped_ptr.hpp>
@@ -241,12 +242,16 @@ main(int argc, char** argv)
     SimulatorReport rep;
     if (!use_deck) {
         // Simple simulation without a deck.
+        PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                          param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                          param.getDefault("poly_amount", poly_props.cMax()));
         SimulatorPolymer simulator(param,
                                    *grid->c_grid(),
                                    *props,
                                    poly_props,
                                    rock_comp->isActive() ? rock_comp.get() : 0,
                                    0, // wells
+                                   polymer_inflow,
                                    src,
                                    bcs.c_bcs(),
                                    linsolver,
@@ -266,6 +271,15 @@ main(int argc, char** argv)
         deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
         simtimer.init(*deck);
         const double total_time = simtimer.totalTime();
+        // Check for WPOLYMER presence in last epoch to decide
+        // polymer injection control type.
+        const bool use_wpolymer = deck->hasField("WPOLYMER");
+        if (use_wpolymer) {
+            if (param.has("poly_start_days")) {
+                MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
+                        "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
+            }
+        }
         for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
             // Set epoch index.
             deck->setCurrentEpoch(epoch);
@@ -287,8 +301,20 @@ main(int argc, char** argv)
                       << "\n                  (number of steps: "
                       << simtimer.numSteps() - step << ")\n\n" << std::flush;
 
-            // Create new wells, well_state
+            // Create new wells, polymer inflow controls.
             WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
+            if (use_wpolymer) {
+                if (wells.c_wells() == 0) {
+                    THROW("Cannot control polymer injection via WPOLYMER without wells.");
+                }
+                polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells()));
+            } else {
+                polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
+                                                            param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
+                                                            param.getDefault("poly_amount", poly_props.cMax())));
+            }
+
             // @@@ HACK: we should really make a new well state and
             // properly transfer old well state to it every epoch,
             // since number of wells may change etc.
@@ -303,6 +329,7 @@ main(int argc, char** argv)
                                        poly_props,
                                        rock_comp->isActive() ? rock_comp.get() : 0,
                                        wells.c_wells(),
+                                       *polymer_inflow,
                                        src,
                                        bcs.c_bcs(),
                                        linsolver,

opm/polymer/SimulatorCompressiblePolymer.cpp

@@ -93,6 +93,7 @@ namespace Opm
              const PolymerProperties& poly_props,
              const RockCompressibility* rock_comp_props,
              const Wells* wells,
+             const PolymerInflowInterface& polymer_inflow,
              const std::vector<double>& src,
              const FlowBoundaryConditions* bcs,
              LinearSolverInterface& linsolver,
@@ -119,6 +120,7 @@ namespace Opm
         const PolymerProperties& poly_props_;
         const RockCompressibility* rock_comp_props_;
         const Wells* wells_;
+        const PolymerInflowInterface& polymer_inflow_;
         const std::vector<double>& src_;
         const FlowBoundaryConditions* bcs_;
         const LinearSolverInterface& linsolver_;
@@ -130,7 +132,6 @@ namespace Opm
         std::vector< std::vector<int> > columns_;
         // Misc. data
         std::vector<int> allcells_;
-        boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
     };
 
 
@@ -142,12 +143,14 @@ namespace Opm
                                                                const PolymerProperties& poly_props,
                                                                const RockCompressibility* rock_comp_props,
                                                                const Wells* wells,
+                                                               const PolymerInflowInterface& polymer_inflow,
                                                                const std::vector<double>& src,
                                                                const FlowBoundaryConditions* bcs,
                                                                LinearSolverInterface& linsolver,
                                                                const double* gravity)
     {
-        pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props, wells, src, bcs, linsolver, gravity));
+        pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props,
+                              wells, polymer_inflow, src, bcs, linsolver, gravity));
     }
 
 
@@ -171,6 +174,7 @@ namespace Opm
                                              const PolymerProperties& poly_props,
                                              const RockCompressibility* rock_comp_props,
                                              const Wells* wells,
+                                             const PolymerInflowInterface& polymer_inflow,
                                              const std::vector<double>& src,
                                              const FlowBoundaryConditions* bcs,
                                              LinearSolverInterface& linsolver,
@@ -180,6 +184,7 @@ namespace Opm
           poly_props_(poly_props),
           rock_comp_props_(rock_comp_props),
           wells_(wells),
+          polymer_inflow_(polymer_inflow),
           src_(src),
           bcs_(bcs),
           linsolver_(linsolver),
@@ -228,11 +233,6 @@ namespace Opm
             extractColumn(grid_, columns_);
         }
 
-        // Polymer inflow control.
-        poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
-                                                  param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
-                                                  param.getDefault("poly_amount", poly_props.cMax())));
-
         // Misc init.
         const int num_cells = grid.number_of_cells;
         allcells_.resize(num_cells);
@@ -324,7 +324,7 @@ namespace Opm
             // Find inflow rate.
             const double current_time = timer.currentTime();
             double stepsize = timer.currentStepLength();
-            poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
+            polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
 
             // Solve transport.
             transport_timer.start();

opm/polymer/SimulatorCompressiblePolymer.hpp

@@ -33,6 +33,7 @@ namespace Opm
     class BlackoilPropertiesInterface;
     class PolymerProperties;
     class RockCompressibility;
+    class PolymerInflowInterface;
     class LinearSolverInterface;
     class SimulatorTimer;
     class PolymerBlackoilState;
@@ -64,6 +65,7 @@ namespace Opm
         /// \param[in] poly_props  polymer properties
         /// \param[in] rock_comp   if non-null, rock compressibility properties
         /// \param[in] wells       if non-null, wells data structure
+        /// \param[in] polymer_inflow  polymer inflow controls
         /// \param[in] src         source terms
         /// \param[in] bcs         boundary conditions, treat as all noflow if null
         /// \param[in] linsolver   linear solver
@@ -74,6 +76,7 @@ namespace Opm
                                      const PolymerProperties& poly_props,
                                      const RockCompressibility* rock_comp_props,
                                      const Wells* wells,
+                                     const PolymerInflowInterface& polymer_inflow,
                                      const std::vector<double>& src,
                                      const FlowBoundaryConditions* bcs,
                                      LinearSolverInterface& linsolver,

opm/polymer/SimulatorPolymer.cpp

@@ -90,6 +90,7 @@ namespace Opm
              const PolymerProperties& poly_props,
              const RockCompressibility* rock_comp_props,
              const Wells* wells,
+             const PolymerInflowInterface& polymer_inflow,
              const std::vector<double>& src,
              const FlowBoundaryConditions* bcs,
              LinearSolverInterface& linsolver,
@@ -116,6 +117,7 @@ namespace Opm
         const PolymerProperties& poly_props_;
         const RockCompressibility* rock_comp_props_;
         const Wells* wells_;
+        const PolymerInflowInterface& polymer_inflow_;
         const std::vector<double>& src_;
         const FlowBoundaryConditions* bcs_;
         const LinearSolverInterface& linsolver_;
@@ -127,7 +129,6 @@ namespace Opm
         std::vector< std::vector<int> > columns_;
         // Misc. data
         std::vector<int> allcells_;
-        boost::scoped_ptr<PolymerInflowInterface> poly_inflow_;
     };
 
 
@@ -139,12 +140,14 @@ namespace Opm
                                        const PolymerProperties& poly_props,
                                        const RockCompressibility* rock_comp_props,
                                        const Wells* wells,
+                                       const PolymerInflowInterface& polymer_inflow,
                                        const std::vector<double>& src,
                                        const FlowBoundaryConditions* bcs,
                                        LinearSolverInterface& linsolver,
                                        const double* gravity)
     {
-        pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props, wells, src, bcs, linsolver, gravity));
+        pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props,
+                              wells, polymer_inflow, src, bcs, linsolver, gravity));
     }
 
 
@@ -168,6 +171,7 @@ namespace Opm
                                  const PolymerProperties& poly_props,
                                  const RockCompressibility* rock_comp_props,
                                  const Wells* wells,
+                                 const PolymerInflowInterface& polymer_inflow,
                                  const std::vector<double>& src,
                                  const FlowBoundaryConditions* bcs,
                                  LinearSolverInterface& linsolver,
@@ -177,6 +181,7 @@ namespace Opm
           poly_props_(poly_props),
           rock_comp_props_(rock_comp_props),
           wells_(wells),
+          polymer_inflow_(polymer_inflow),
           src_(src),
           bcs_(bcs),
           linsolver_(linsolver),
@@ -224,11 +229,6 @@ namespace Opm
             extractColumn(grid_, columns_);
         }
 
-        // Polymer inflow control.
-        poly_inflow_.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
-                                                  param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
-                                                  param.getDefault("poly_amount", poly_props.cMax())));
-
         // Misc init.
         const int num_cells = grid.number_of_cells;
         allcells_.resize(num_cells);
@@ -320,7 +320,7 @@ namespace Opm
             // Find inflow rate.
             const double current_time = timer.currentTime();
             double stepsize = timer.currentStepLength();
-            poly_inflow_->getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
+            polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
 
             // Solve transport.
             transport_timer.start();

opm/polymer/SimulatorPolymer.hpp

@@ -33,6 +33,7 @@ namespace Opm
     class IncompPropertiesInterface;
     class PolymerProperties;
     class RockCompressibility;
+    class PolymerInflowInterface;
     class LinearSolverInterface;
     class SimulatorTimer;
     class PolymerState;
@@ -64,6 +65,7 @@ namespace Opm
         /// \param[in] poly_props  polymer properties
         /// \param[in] rock_comp   if non-null, rock compressibility properties
         /// \param[in] wells       if non-null, wells data structure
+        /// \param[in] polymer_inflow  polymer inflow controls
         /// \param[in] src         source terms
         /// \param[in] bcs         boundary conditions, treat as all noflow if null
         /// \param[in] linsolver   linear solver
@@ -74,6 +76,7 @@ namespace Opm
                         const PolymerProperties& poly_props,
                         const RockCompressibility* rock_comp_props,
                         const Wells* wells,
+                        const PolymerInflowInterface& polymer_inflow,
                         const std::vector<double>& src,
                         const FlowBoundaryConditions* bcs,
                         LinearSolverInterface& linsolver,