From ffd33be55d204d413769d2c7d0e838a0f704411f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Wed, 13 Jun 2012 14:18:32 +0200
Subject: [PATCH] Changed interface of TransportModelPolymer to expect full
 saturation, not just water.

---
 examples/polymer_reorder.cpp          | 13 +++----------
 opm/polymer/TransportModelPolymer.cpp | 20 ++++++++++++--------
 opm/polymer/TransportModelPolymer.hpp |  8 ++++----
 3 files changed, 19 insertions(+), 22 deletions(-)

diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp
index 03793ccf9..8406b6a1a 100644
--- a/examples/polymer_reorder.cpp
+++ b/examples/polymer_reorder.cpp
@@ -509,12 +509,8 @@ main(int argc, char** argv)
     }
     double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
 
-    // We need a separate reorder_sat, because the reorder
-    // code expects a scalar sw, not both sw and so.
-    std::vector<double> reorder_sat(num_cells);
-    std::vector<double> src(num_cells, 0.0);
-
     // Initialising src
+    std::vector<double> src(num_cells, 0.0);
     if (wells->c_wells()) {
         // Do nothing, wells will be the driving force, not source terms.
         // Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
@@ -746,19 +742,16 @@ main(int argc, char** argv)
         }
         for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
             if (use_reorder) {
-                Opm::toWaterSat(state.saturation(), reorder_sat);
                 reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0], stepsize, inflow_c,
-                                    &reorder_sat[0], &state.concentration()[0], &state.maxconcentration()[0]);
-                Opm::toBothSat(reorder_sat, state.saturation());
+                                    state.saturation(), state.concentration(), state.maxconcentration());
                 Opm::computeInjectedProduced(*props, polyprop, state.saturation(), state.concentration(), state.maxconcentration(),
                                              reorder_src, simtimer.currentStepLength(), inflow_c,
                                              injected, produced, polyinj, polyprod);
                 if (use_segregation_split) {
                     if (use_column_solver) {
                         if (use_gauss_seidel_gravity) {
-                            reorder_model.solveGravity(columns, &porevol[0], stepsize, reorder_sat,
+                            reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation(),
                                                        state.concentration(), state.maxconcentration());
-                            Opm::toBothSat(reorder_sat, state.saturation());
                         } else {
                             colsolver.solve(columns, stepsize, state.saturation(), state.concentration(),
                                             state.maxconcentration());
diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp
index 37da58de8..932e3b23a 100644
--- a/opm/polymer/TransportModelPolymer.cpp
+++ b/opm/polymer/TransportModelPolymer.cpp
@@ -22,6 +22,7 @@
 #include <opm/core/fluid/IncompPropertiesInterface.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/utility/RootFinders.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 #include <cmath>
 
@@ -188,7 +189,6 @@ namespace Opm
 	  source_(0),
 	  dt_(0.0),
 	  inflow_c_(0.0),
-	  saturation_(0),
 	  concentration_(0),
 	  cmax_(0),
 	  fractionalflow_(grid.number_of_cells, -1.0),
@@ -219,19 +219,20 @@ namespace Opm
 				      const double* source,
 				      const double dt,
 				      const double inflow_c,
-				      double* saturation,
-				      double* concentration,
-				      double* cmax)
+				      std::vector<double>& saturation,
+				      std::vector<double>& concentration,
+				      std::vector<double>& cmax)
     {
 	darcyflux_ = darcyflux;
         porevolume_ = porevolume;
 	source_ = source;
 	dt_ = dt;
 	inflow_c_ = inflow_c;
-	saturation_ = saturation;
-	concentration_ = concentration;
-	cmax_ = cmax;
+        toWaterSat(saturation, saturation_);
+	concentration_ = &concentration[0];
+	cmax_ = &cmax[0];
 	reorderAndTransport(grid_, darcyflux);
+        toBothSat(saturation_, saturation);
     }
 
 
@@ -1198,7 +1199,7 @@ namespace Opm
         // initialize variables.
         porevolume_ = porevolume;
         dt_ = dt;
-        saturation_ = &saturation[0];
+        toWaterSat(saturation, saturation_);
         concentration_ = &concentration[0];
         cmax_ = &cmax[0];
         const int nc = grid_.number_of_cells;
@@ -1215,12 +1216,15 @@ namespace Opm
 
         // Solve on all columns.
         int num_iters = 0;
+        // std::cout << "Gauss-Seidel column solver # columns: " << columns.size() << std::endl;
         for (std::vector<std::vector<int> >::size_type i = 0; i < columns.size(); i++) {
             // std::cout << "==== new column" << std::endl;
             num_iters += solveGravityColumn(columns[i]);
         }
         std::cout << "Gauss-Seidel column solver average iterations: "
                   << double(num_iters)/double(columns.size()) << std::endl;
+
+        toBothSat(saturation_, saturation);
     }
 
 } // namespace Opm
diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp
index 0f1017787..3d2c12324 100644
--- a/opm/polymer/TransportModelPolymer.hpp
+++ b/opm/polymer/TransportModelPolymer.hpp
@@ -58,9 +58,9 @@ namespace Opm
 		   const double* source,
 		   const double dt,
 		   const double inflow_c,
-		   double* saturation,
-		   double* concentration,
-		   double* cmax);
+		   std::vector<double>& saturation,
+                   std::vector<double>& concentration,
+                   std::vector<double>& cmax);
 
 	virtual void solveSingleCell(const int cell);
 	virtual void solveMultiCell(const int num_cells, const int* cells);
@@ -96,7 +96,7 @@ namespace Opm
 	const double* source_;      // one source per cell
 	double dt_;
 	double inflow_c_;
-	double* saturation_;        // one per cell
+        std::vector<double> saturation_; // one per cell, only water saturation!
 	double* concentration_;
 	double* cmax_;
 	std::vector<double> fractionalflow_;  // one per cell