Merge pull request #5 from atgeirr/master

Add well management to polymer simulators.

examples/sim_poly2p_comp_reorder.cpp

@@ -76,7 +76,7 @@ main(int argc, char** argv)
 {
     using namespace Opm;
 
-    std::cout << "\n================    Test program for incompressible two-phase flow with polymer    ===============\n\n";
+    std::cout << "\n================    Test program for weakly compressible two-phase flow with polymer    ===============\n\n";
     parameter::ParameterGroup param(argc, argv, false);
     std::cout << "---------------    Reading parameters     ---------------" << std::endl;
 
@@ -241,12 +241,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;
         SimulatorCompressiblePolymer simulator(param,
                                                *grid->c_grid(),
                                                *props,
                                                poly_props,
                                                rock_comp->isActive() ? rock_comp.get() : 0,
-                                               0, // wells
+                                               wells,
                                                polymer_inflow,
                                                src,
                                                bcs.c_bcs(),
@@ -324,7 +325,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(),

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/SimulatorCompressiblePolymer.cpp

@@ -17,7 +17,6 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-
 #if HAVE_CONFIG_H
 #include "config.h"
 #endif // HAVE_CONFIG_H
@@ -39,6 +38,7 @@
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
+#include <opm/core/wells/WellsManager.hpp>
 
 #include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
 #include <opm/core/fluid/RockCompressibility.hpp>
@@ -52,7 +52,7 @@
 #include <opm/polymer/PolymerProperties.hpp>
 #include <opm/polymer/polymerUtilities.hpp>
 
-#include <boost/filesystem/convenience.hpp>
+#include <boost/filesystem.hpp>
 #include <boost/scoped_ptr.hpp>
 #include <boost/lexical_cast.hpp>
 
@@ -92,7 +92,7 @@ namespace Opm
              const BlackoilPropertiesInterface& 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,
@@ -111,6 +111,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_;
@@ -119,11 +122,11 @@ namespace Opm
         const BlackoilPropertiesInterface& 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_;
         const FlowBoundaryConditions* bcs_;
-        const LinearSolverInterface& linsolver_;
         const double* gravity_;
         // Solvers
         CompressibleTpfaPolymer psolver_;
@@ -142,7 +145,7 @@ namespace Opm
                                                                const BlackoilPropertiesInterface& 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,
@@ -150,7 +153,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));
     }
 
 
@@ -173,7 +176,7 @@ namespace Opm
                                              const BlackoilPropertiesInterface& 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,
@@ -183,18 +186,18 @@ 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),
-          linsolver_(linsolver),
           gravity_(gravity),
           psolver_(grid, props, rock_comp_props, poly_props, 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),
+                   gravity, wells_manager.c_wells()),
-          tsolver_(grid, props, poly_props, *rock_comp_props,
+          tsolver_(grid, props, poly_props,
                    TransportModelCompressiblePolymer::Bracketing,
                    param.getDefault("nl_tolerance", 1e-9),
                    param.getDefault("nl_maxiter", 30))
@@ -215,6 +218,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.
         TransportModelCompressiblePolymer::SingleCellMethod method;
         std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
@@ -228,7 +235,7 @@ namespace Opm
         tsolver_.setPreferredMethod(method);
         num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
         use_segregation_split_ = param.getDefault("use_segregation_split", false);
-        if (gravity_ != 0 && use_segregation_split_){
+        if (gravity != 0 && use_segregation_split_) {
             tsolver_.initGravity(gravity);
             extractColumn(grid_, columns_);
         }
@@ -261,7 +268,6 @@ namespace Opm
         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;
         double ptime = 0.0;
@@ -301,15 +307,75 @@ namespace Opm
 
             initial_pressure = state.pressure();
 
-            // Solve pressure.
+            // Solve pressure equation.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, poly_props_, allcells_,
+                                      state.pressure(), state.surfacevol(), 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();
                 psolver_.solve(timer.currentStepLength(), state, well_state);
+
+                // Renormalize pressure if both fluids and rock are
+                // 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 (psolver_.singularPressure()) {
+                    // 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()) {

opm/polymer/SimulatorCompressiblePolymer.hpp

@@ -33,6 +33,7 @@ namespace Opm
     class BlackoilPropertiesInterface;
     class PolymerProperties;
     class RockCompressibility;
+    class WellsManager;
     class PolymerInflowInterface;
     class LinearSolverInterface;
     class SimulatorTimer;
@@ -63,8 +64,8 @@ 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] wells_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 +76,7 @@ namespace Opm
                                      const BlackoilPropertiesInterface& 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/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>
 
@@ -50,7 +52,7 @@
 #include <opm/polymer/PolymerProperties.hpp>
 #include <opm/polymer/polymerUtilities.hpp>
 
-#include <boost/filesystem/convenience.hpp>
+#include <boost/filesystem.hpp>
 #include <boost/scoped_ptr.hpp>
 #include <boost/lexical_cast.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,11 +124,11 @@ 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_;
         const FlowBoundaryConditions* bcs_;
-        const LinearSolverInterface& linsolver_;
         const double* gravity_;
         // Solvers
         IncompTpfaPolymer psolver_;
@@ -139,7 +147,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 +155,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 +178,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,17 +188,17 @@ 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),
-          linsolver_(linsolver),
           gravity_(gravity),
           psolver_(grid, props, rock_comp_props, poly_props, 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, poly_props, TransportModelPolymer::Bracketing,
                    param.getDefault("nl_tolerance", 1e-9),
                    param.getDefault("nl_maxiter", 30))
@@ -211,6 +219,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 +266,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 +311,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 +407,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 +609,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,8 @@ 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] wells_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 +76,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/TransportModelCompressiblePolymer.cpp

@@ -154,14 +154,12 @@ namespace Opm
     TransportModelCompressiblePolymer::TransportModelCompressiblePolymer(const UnstructuredGrid& grid,
                                                                          const BlackoilPropertiesInterface& props,
                                                                          const PolymerProperties& polyprops,
-                                                                         const RockCompressibility& rock_comp,
                                                                          const SingleCellMethod method,
                                                                          const double tol,
                                                                          const int maxit)
         : grid_(grid),
           props_(props),
           polyprops_(polyprops),
-          rock_comp_(rock_comp),
           darcyflux_(0),
           porevolume0_(0),
           porevolume_(0),

opm/polymer/TransportModelCompressiblePolymer.hpp

@@ -20,7 +20,6 @@
 #ifndef OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
 #define OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
 
-#include <opm/core/fluid/RockCompressibility.hpp>
 #include <opm/polymer/PolymerProperties.hpp>
 #include <opm/core/transport/reorder/TransportModelInterface.hpp>
 #include <opm/core/utility/linearInterpolation.hpp>
@@ -63,7 +62,6 @@ namespace Opm
 	TransportModelCompressiblePolymer(const UnstructuredGrid& grid,
                                           const BlackoilPropertiesInterface& props,
                                           const PolymerProperties& polyprops,
-                                          const RockCompressibility& rock_comp,
                                           const SingleCellMethod method,
                                           const double tol,
                                           const int maxit);
@@ -134,7 +132,6 @@ namespace Opm
 	const UnstructuredGrid& grid_;
 	const BlackoilPropertiesInterface& props_;
 	const PolymerProperties& polyprops_;
-        const RockCompressibility& rock_comp_;
 	const double* darcyflux_;   // one flux per grid face
         const double* porevolume0_; // one volume per cell
         const double* porevolume_;  // one volume per cell

opm/polymer/polymerUtilities.cpp

@@ -91,6 +91,84 @@ 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]);
+        }
+    }
+
+
+    /// 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]  p                pressure (one value per cell)
+    /// @param[in]  z                surface-volume values (for all P phases)
+    /// @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::BlackoilPropertiesInterface& props,
+                               const Opm::PolymerProperties& polyprops,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& p,
+                               const std::vector<double>& z,
+                               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);
+	std::vector<double> mu(num_cells*num_phases);
+	props.viscosity(num_phases, &p[0], &z[0], &cells[0], &mu[0], 0);
+        double mob[2]; // here we assume num_phases=2
+	for (int cell = 0; cell < num_cells; ++cell) {
+            double* kr_cell = &kr[2*cell];
+            double* mu_cell = &mu[2*cell];
+            polyprops.effectiveMobilities(c[cell], cmax[cell], mu_cell, 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,42 @@ 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);
+
+    /// 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]  p                pressure (one value per cell)
+    /// @param[in]  z                surface-volume values (for all P phases)
+    /// @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::BlackoilPropertiesInterface& props,
+                               const Opm::PolymerProperties& polyprops,
+                               const std::vector<int>& cells,
+                               const std::vector<double>& p,
+                               const std::vector<double>& z,
+                               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.