Added well management to SimulatorPolymer.

A polymer-specific overload of computeFractionalFlow()
has been added in support of this.

Also added a minor feature: pressure normalization for
situations with arbitrary absolute pressure.

examples/sim_poly2p_incomp_reorder.cpp

@@ -245,12 +245,13 @@ main(int argc, char** argv)
         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()));
+        WellsManager wells;
         SimulatorPolymer simulator(param,
                                    *grid->c_grid(),
                                    *props,
                                    poly_props,
                                    rock_comp->isActive() ? rock_comp.get() : 0,
-                                   0, // wells
+                                   wells,
                                    polymer_inflow,
                                    src,
                                    bcs.c_bcs(),
@@ -328,7 +329,7 @@ main(int argc, char** argv)
                                        *props,
                                        poly_props,
                                        rock_comp->isActive() ? rock_comp.get() : 0,
-                                       wells.c_wells(),
+                                       wells,
                                        *polymer_inflow,
                                        src,
                                        bcs.c_bcs(),

opm/polymer/SimulatorPolymer.cpp

@@ -38,6 +38,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>
 
@@ -77,6 +79,9 @@ namespace Opm
                             const std::string& output_dir);
         void outputWellReport(const Opm::WellReport& wellreport,
                               const std::string& output_dir);
+        bool allNeumannBCs(const FlowBoundaryConditions* bcs);
+        bool allRateWells(const Wells* wells);
+
     } // anonymous namespace
 
 
@@ -89,7 +94,7 @@ namespace Opm
              const IncompPropertiesInterface& props,
              const PolymerProperties& poly_props,
              const RockCompressibility* rock_comp_props,
-             const Wells* wells,
+             WellsManager& wells_manager,
              const PolymerInflowInterface& polymer_inflow,
              const std::vector<double>& src,
              const FlowBoundaryConditions* bcs,
@@ -108,6 +113,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_;
@@ -116,6 +124,7 @@ namespace Opm
         const IncompPropertiesInterface& props_;
         const PolymerProperties& poly_props_;
         const RockCompressibility* rock_comp_props_;
+        WellsManager& wells_manager_;
         const Wells* wells_;
         const PolymerInflowInterface& polymer_inflow_;
         const std::vector<double>& src_;
@@ -139,7 +148,7 @@ namespace Opm
                                        const IncompPropertiesInterface& props,
                                        const PolymerProperties& poly_props,
                                        const RockCompressibility* rock_comp_props,
-                                       const Wells* wells,
+                                       WellsManager& wells_manager,
                                        const PolymerInflowInterface& polymer_inflow,
                                        const std::vector<double>& src,
                                        const FlowBoundaryConditions* bcs,
@@ -147,7 +156,7 @@ namespace Opm
                                        const double* gravity)
     {
         pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props,
-                              wells, polymer_inflow, src, bcs, linsolver, gravity));
+                              wells_manager, polymer_inflow, src, bcs, linsolver, gravity));
     }
 
 
@@ -170,7 +179,7 @@ namespace Opm
                                  const IncompPropertiesInterface& props,
                                  const PolymerProperties& poly_props,
                                  const RockCompressibility* rock_comp_props,
-                                 const Wells* wells,
+                                 WellsManager& wells_manager,
                                  const PolymerInflowInterface& polymer_inflow,
                                  const std::vector<double>& src,
                                  const FlowBoundaryConditions* bcs,
@@ -180,7 +189,8 @@ namespace Opm
           props_(props),
           poly_props_(poly_props),
           rock_comp_props_(rock_comp_props),
-          wells_(wells),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
           polymer_inflow_(polymer_inflow),
           src_(src),
           bcs_(bcs),
@@ -190,7 +200,7 @@ namespace Opm
                    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, poly_props, TransportModelPolymer::Bracketing,
                    param.getDefault("nl_tolerance", 1e-9),
                    param.getDefault("nl_maxiter", 30))
@@ -211,6 +221,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.
         TransportModelPolymer::SingleCellMethod method;
         std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
@@ -254,7 +268,7 @@ namespace Opm
             computePorevolume(grid_, props_.porosity(), porevol);
         }
         const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
-
+        std::vector<double> initial_porevol = porevol;
 
         // Main simulation loop.
         Opm::time::StopWatch pressure_timer;
@@ -299,17 +313,78 @@ namespace Opm
             }
 
             // Solve pressure.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, poly_props_, allcells_,
+                                      state.saturation(), state.concentration(), state.maxconcentration(),
+                                      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_props_ == NULL || !rock_comp_props_->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;
-            } while (false);
+
+                // 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_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_props_ && rock_comp_props_->isActive()) {
+                initial_porevol = porevol;
                 computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
             }
 
@@ -334,7 +409,7 @@ namespace Opm
             double substep_polyprod = 0.0;
             injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
             for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
+                tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
                                state.saturation(), state.concentration(), state.maxconcentration());
                 Opm::computeInjectedProduced(props_, poly_props_,
                                              state,
@@ -536,6 +611,35 @@ namespace Opm
         }
 
 
+        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;
+            }
+        }
+
+
+        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;
+        }
+
+
     } // anonymous namespace
 
 

opm/polymer/SimulatorPolymer.hpp

@@ -33,6 +33,7 @@ namespace Opm
     class IncompPropertiesInterface;
     class PolymerProperties;
     class RockCompressibility;
+    class WellsManager;
     class PolymerInflowInterface;
     class LinearSolverInterface;
     class SimulatorTimer;
@@ -63,8 +64,9 @@ namespace Opm
         /// \param[in] grid        grid data structure
         /// \param[in] props       fluid and rock properties
         /// \param[in] poly_props  polymer properties
-        /// \param[in] rock_comp   if non-null, rock compressibility properties
+        /// \param[in] rock_comp_props   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] polymer_inflow  polymer inflow controls
         /// \param[in] src         source terms
         /// \param[in] bcs         boundary conditions, treat as all noflow if null
@@ -75,7 +77,7 @@ namespace Opm
                         const IncompPropertiesInterface& props,
                         const PolymerProperties& poly_props,
                         const RockCompressibility* rock_comp_props,
-                        const Wells* wells,
+                        WellsManager& wells_manager,
                         const PolymerInflowInterface& polymer_inflow,
                         const std::vector<double>& src,
                         const FlowBoundaryConditions* bcs,

opm/polymer/polymerUtilities.cpp

@@ -91,6 +91,42 @@ namespace Opm
     }
 
 
+    /// Computes the fractional flow for each cell in the cells argument
+    /// @param[in]  props            rock and fluid properties
+    /// @param[in]  polyprops        polymer properties
+    /// @param[in]  cells            cells with which the saturation values are associated
+    /// @param[in]  s                saturation values (for all phases)
+    /// @param[in]  c                concentration values
+    /// @param[in]  cmax             max polymer concentration experienced by cell
+    /// @param[out] fractional_flow  the fractional flow for each phase for each cell.
+    void computeFractionalFlow(const Opm::IncompPropertiesInterface& props,
+                               const Opm::PolymerProperties& polyprops,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& s,
+                               const std::vector<double>& c,
+                               const std::vector<double>& cmax,
+                               std::vector<double>& fractional_flows)
+    {
+	int num_cells = cells.size();
+	int num_phases = props.numPhases();
+        if (num_phases != 2) {
+            THROW("computeFractionalFlow() assumes 2 phases.");
+        }
+	fractional_flows.resize(num_cells*num_phases);
+	ASSERT(int(s.size()) == num_cells*num_phases);
+	std::vector<double> kr(num_cells*num_phases);
+	props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
+	const double* visc = props.viscosity();
+        double mob[2]; // here we assume num_phases=2
+	for (int cell = 0; cell < num_cells; ++cell) {
+            double* kr_cell = &kr[2*cell];
+            polyprops.effectiveMobilities(c[cell], cmax[cell], visc, kr_cell, mob);
+            fractional_flows[2*cell]     = mob[0] / (mob[0] + mob[1]);
+            fractional_flows[2*cell + 1] = mob[1] / (mob[0] + mob[1]);
+        }
+    }
+
+
     /// @brief Computes injected and produced volumes of all phases,
     ///        and injected and produced polymer mass.
     /// Note 1: assumes that only the first phase is injected.

opm/polymer/polymerUtilities.hpp

@@ -56,6 +56,7 @@ namespace Opm
     /// @param[in]  cells     cells with which the saturation values are associated
     /// @param[in]  s         saturation values (for all phases)
     /// @param[in]  c         polymer concentration
+    /// @param[in]  cmax      max polymer concentration experienced by cell
     /// @param[out] totmob    total mobility
     /// @param[out] omega     mobility-weighted (or fractional-flow weighted)
     ///                       fluid densities.
@@ -68,6 +69,22 @@ namespace Opm
 				   std::vector<double>& totmob,
 				   std::vector<double>& omega);
 
+    /// Computes the fractional flow for each cell in the cells argument
+    /// @param[in]  props            rock and fluid properties
+    /// @param[in]  polyprops        polymer properties
+    /// @param[in]  cells            cells with which the saturation values are associated
+    /// @param[in]  s                saturation values (for all phases)
+    /// @param[in]  c                concentration values
+    /// @param[in]  cmax             max polymer concentration experienced by cell
+    /// @param[out] fractional_flow  the fractional flow for each phase for each cell.
+    void computeFractionalFlow(const Opm::IncompPropertiesInterface& props,
+                               const Opm::PolymerProperties& polyprops,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& s,
+                               const std::vector<double>& c,
+                               const std::vector<double>& cmax,
+                               std::vector<double>& fractional_flows);
+
     /// @brief Computes injected and produced volumes of all phases,
     ///        and injected and produced polymer mass.
     /// Note 1: assumes that only the first phase is injected.