output watercut for incomp polymer solver.

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp

@@ -150,7 +150,8 @@ namespace {
     step(const double   dt,
          PolymerState&  x,
          WellState&     xw,
-         const std::vector<double>& polymer_inflow)
+         const std::vector<double>& polymer_inflow,
+		 std::vector<double>& src)
     {
         
         V pvol(grid_.number_of_cells);
@@ -169,7 +170,7 @@ namespace {
         const double atol  = 1.0e-12;
         const double rtol  = 5.0e-8;
         const int    maxit = 40;
-        assemble(pvdt, old_state, x, xw, polymer_inflow);
+        assemble(pvdt, x, xw, polymer_inflow, src);
 
         const double r0  = residualNorm();
         int          it  = 0;
@@ -181,7 +182,7 @@ namespace {
             const V dx = solveJacobianSystem();
             updateState(dx, x, xw);
 
-            assemble(pvdt, old_state, x, xw, polymer_inflow);
+            assemble(pvdt, x, xw, polymer_inflow, src);
 
             const double r = residualNorm();
 
@@ -393,10 +394,10 @@ namespace {
     void
     FullyImplicitTwophasePolymerSolver::
     assemble(const V&             pvdt,
-             const SolutionState& old_state,
              const PolymerState&  x,
              const WellState&     xw,
-             const std::vector<double>& polymer_inflow)
+             const std::vector<double>& polymer_inflow,
+			 std::vector<double>& src)
     {
         // Create the primary variables.
         const SolutionState state = variableState(x, xw);
@@ -404,27 +405,10 @@ namespace {
         // -------- Mass balance equations for water and oil --------
         const V trans = subset(transmissibility(), ops_.internal_faces);
         const std::vector<ADB> kr = computeRelPerm(state);
-
 		const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
- //       const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
         const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
         const ADB mc = computeMc(state);
         computeMassFlux(trans, mc, kr[1], krw_eff, state);
-        //const std::vector<ADB> source = accumSource(kr[1], krw_eff, state.concentration, src, polymer_inflow);
-       // const std::vector<ADB> source = polymerSource();
- //       const double rho_r = polymer_props_ad_.rockDensity();
-//        const V phi = V::Constant(pvdt.size(), 1, *fluid_.porosity());
-
-//        const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
-//        residual_.mass_balance[0] =  pvdt * (state.saturation[0] - old_state.saturation[0])
-//                        + ops_.div * mflux[0];
-//        residual_.mass_balance[1] =  pvdt * (state.saturation[1] - old_state.saturation[1])
-//                        + ops_.div * mflux[1];
-      // Mass balance equation for polymer
-//        residual_.mass_balance[2] = pvdt * (state.saturation[0] * state.concentration
-//                      - old_state.saturation[0] * old_state.concentration) * (1. - dead_pore_vol)
-//                      + pvdt * rho_r * (1. - phi) / phi * ads
- //                     + ops_.div * mflux[2];
         residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
                                     + ops_.div*rq_[0].mflux;
         residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
@@ -516,14 +500,17 @@ namespace {
             well_contribs[phase] = superset(perf_flux, well_cells, nc);
             // DUMP(well_contribs[phase]);
             residual_.mass_balance[phase] += well_contribs[phase];
+			for (int i = 0; i < nc; ++i) {
+				src[i] += well_contribs[phase].value()[i];
+			}
         }
 
         // well rates contribs to polymer mass balance eqn.
         // for injection wells.
         const V polyin = Eigen::Map<const V>(& polymer_inflow[0], nc);
         const V poly_in_perf = subset(polyin, well_cells);
-        const V poly_c_cell = subset(state.concentration, well_cells).value();
+        const V poly_mc_cell = subset(mc, well_cells).value();
-        const V poly_c = producer.select(poly_c_cell, poly_in_perf);
+        const V poly_c = producer.select(poly_mc_cell, poly_in_perf);
         residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc);
 
         // Set the well flux equation

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp

@@ -29,7 +29,8 @@ namespace Opm {
         void step(const double   dt,
                   PolymerState& state,
                   WellState&    well_state,
-                  const std::vector<double>& polymer_inflow);
+                  const std::vector<double>& polymer_inflow,
+				  std::vector<double>& src);
     private:
         typedef AutoDiffBlock<double> ADB;
         typedef ADB::V V;
@@ -88,10 +89,10 @@ namespace Opm {
                       const WellState&    xw);
         void
         assemble(const V&               pvdt,
-                 const SolutionState&   old_state,
                  const PolymerState&    x,
                  const WellState&       xw,
-                 const std::vector<double>& polymer_inflow);
+                 const std::vector<double>& polymer_inflow,
+				 std::vector<double>& src);
         V solveJacobianSystem() const;
         void updateState(const V&             dx,
                          PolymerState& x,

opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp

@@ -24,6 +24,7 @@
 #include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
 #include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
 #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/polymer/fullyimplicit/utilities.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
 #include <opm/core/wells/WellsManager.hpp>
@@ -146,6 +147,7 @@ namespace Opm
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
+        dm["cmax"] = &state.maxconcentration();
         dm["concentration"] = &state.concentration();
         std::vector<double> cell_velocity;
         Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@@ -162,6 +164,7 @@ namespace Opm
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
+        dm["cmax"] = &state.maxconcentration();
         dm["concentration"] = &state.concentration();
         std::vector<double> cell_velocity;
         Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@@ -191,7 +194,6 @@ namespace Opm
 
 
     
-/*
     static void outputWaterCut(const Opm::Watercut& watercut,
                                const std::string& output_dir)
     {
@@ -203,6 +205,7 @@ namespace Opm
         }
         watercut.write(os);
     }
+/*
     static void outputWellReport(const Opm::WellReport& wellreport,
                                  const std::string& output_dir)
     {
@@ -298,8 +301,9 @@ namespace Opm
         std::vector<double> porevol;
         Opm::computePorevolume(grid_, props_.porosity(), porevol);
 
-        // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
         std::vector<double> polymer_inflow_c(grid_.number_of_cells);
+		std::vector<double> transport_src(grid_.number_of_cells);
 
         // Main simulation loop.
         Opm::time::StopWatch solver_timer;
@@ -307,6 +311,10 @@ namespace Opm
         Opm::time::StopWatch step_timer;
         Opm::time::StopWatch total_timer;
         total_timer.start();
+        double tot_injected[2] = { 0.0 };
+        double tot_produced[2] = { 0.0 };
+        Opm::Watercut watercut;
+        watercut.push(0.0, 0.0, 0.0);
 #if 0
         // These must be changed for three-phase.
         double init_surfvol[2] = { 0.0 };
@@ -333,7 +341,7 @@ namespace Opm
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
-                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+        //        outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
 
             }
 
@@ -348,7 +356,7 @@ namespace Opm
                 polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
                 solver_timer.start();
                 std::vector<double> initial_pressure = state.pressure();
-                solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
+                solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
 
                 // Stop timer and report.
                 solver_timer.stop();
@@ -378,6 +386,38 @@ namespace Opm
                     }
                 }
             } while (!well_control_passed);
+         	double injected[2] = { 0.0 };
+          	double produced[2] = { 0.0 };
+			double polyinj = 0;
+			double polyprod = 0;
+
+            Opm::computeInjectedProduced(props_, polymer_props_,
+                                         state,
+                                         transport_src, polymer_inflow_c, timer.currentStepLength(),
+                                         injected, produced,
+                                         polyinj, polyprod);
+            tot_injected[0] += injected[0];
+            tot_injected[1] += injected[1];
+            tot_produced[0] += produced[0];
+            tot_produced[1] += produced[1];
+         	watercut.push(timer.currentTime() + timer.currentStepLength(),
+                          	  produced[0]/(produced[0] + produced[1]),
+                          	  tot_produced[0]/tot_porevol_init);
+            std::cout.precision(5);
+            const int width = 18;
+            std::cout << "\nMass balance report.\n";
+            std::cout << "    Injected reservoir volumes:      "
+                      << std::setw(width) << injected[0]
+                      << std::setw(width) << injected[1] << std::endl;
+            std::cout << "    Produced reservoir volumes:      "
+                      << std::setw(width) << produced[0]
+                      << std::setw(width) << produced[1] << std::endl;
+            std::cout << "    Total inj reservoir volumes:     "
+                      << std::setw(width) << tot_injected[0]
+                      << std::setw(width) << tot_injected[1] << std::endl;
+            std::cout << "    Total prod reservoir volumes:    "
+                      << std::setw(width) << tot_produced[0]
+                      << std::setw(width) << tot_produced[1] << std::endl;
 
 
             // Update pore volumes if rock is compressible.
@@ -440,9 +480,9 @@ namespace Opm
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
-                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
-#if 0
                 outputWaterCut(watercut, output_dir_);
+#if 0
+                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
                 if (wells_) {
                     outputWellReport(wellreport, output_dir_);
                 }

opm/polymer/fullyimplicit/utilities.cpp

@@ -144,7 +144,7 @@ namespace Opm
 		const V inv_muw_eff = polymer_props.effectiveInvWaterVisc(c, mus);
 		std::vector<V> mob(np);
 		mob[0] = krw_eff * inv_muw_eff;
-		mob[1] = kr[1] / mu[1];
+		mob[1] = kr[1] / mus[1];
 		
 		const V watmob_c = src_selector.select(mob[0], one);
 		const V oilmob_c = src_selector.select(mob[1], zero);