Moved AdHocProps to new Inc.Props.DefaultPolymer class.

Also moved some functions in polymer_reorder.cpp to be more similar to spu_2p.cpp.

Makefile.am

@@ -13,6 +13,7 @@ opm/polymer/polymerUtilities.cpp
 nobase_include_HEADERS =				\
 opm/polymer/GravityColumnSolverPolymer.hpp		\
 opm/polymer/GravityColumnSolverPolymer_impl.hpp		\
+opm/polymer/IncompPropertiesDefaultPolymer.hpp		\
 opm/polymer/PolymerProperties.hpp			\
 opm/polymer/PolymerState.hpp				\
 opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp	\

examples/polymer_reorder.cpp

@@ -38,7 +38,7 @@
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
-#include <opm/core/fluid/IncompPropertiesBasic.hpp>
+#include <opm/polymer/IncompPropertiesDefaultPolymer.hpp>
 #include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
 
 #include <opm/core/linalg/LinearSolverUmfpack.hpp>
@@ -84,90 +84,54 @@
 
 
 
-class AdHocProps : public Opm::IncompPropertiesBasic
+
+static void outputState(const UnstructuredGrid& grid,
+                        const Opm::PolymerState& state,
+                        const int step,
+                        const std::string& output_dir)
 {
-public:
+    // Write data in VTK format.
-    AdHocProps(const Opm::parameter::ParameterGroup& param, int dim, int num_cells)
+    std::ostringstream vtkfilename;
-        : Opm::IncompPropertiesBasic(param, dim, num_cells)
+    vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
-    {
+    std::ofstream vtkfile(vtkfilename.str().c_str());
-        ASSERT(numPhases() == 2);
+    if (!vtkfile) {
-        sw_.resize(3);
+        THROW("Failed to open " << vtkfilename.str());
-        sw_[0] = 0.2;
-        sw_[1] = 0.7;
-        sw_[2] = 1.0;
-        krw_.resize(3);
-        krw_[0] = 0.0;
-        krw_[1] = 0.7;
-        krw_[2] = 1.0;
-        so_.resize(2);
-        so_[0] = 0.3;
-        so_[1] = 0.8;
-        kro_.resize(2);
-        kro_[0] = 0.0;
-        kro_[1] = 1.0;
     }
+    Opm::DataMap dm;
+    dm["saturation"] = &state.saturation();
+    dm["pressure"] = &state.pressure();
+    dm["concentration"] = &state.concentration();
+    dm["cmax"] = &state.maxconcentration();
+    std::vector<double> cell_velocity;
+    Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+    dm["velocity"] = &cell_velocity;
+    Opm::writeVtkData(grid, dm, vtkfile);
 
-    virtual void relperm(const int n,
+    // Write data (not grid) in Matlab format
-                         const double* s,
+    for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
-                         const int* /*cells*/,
+        std::ostringstream fname;
-                         double* kr,
+        fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
-                         double* dkrds) const
+        std::ofstream file(fname.str().c_str());
-    {
+        if (!file) {
-        // ASSERT(dkrds == 0);
+            THROW("Failed to open " << fname.str());
-        // We assume two phases flow
-        for (int i = 0; i < n; ++i) {
-            kr[2*i] = krw(s[2*i]);
-            kr[2*i+1] = kro(s[2*i+1]);
-            if (dkrds != 0) {
-                dkrds[2*i] = krw_dsw(s[2*i]);
-                dkrds[2*i+3] = kro_dso(s[2*i+1]);
-                dkrds[2*i+1] = -dkrds[2*i+3];
-                dkrds[2*i+2] = -dkrds[2*i];
         }
+        const std::vector<double>& d = *(it->second);
+        std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
     }
-    }
+}
-
-    virtual void satRange(const int n,
-                          const int* /*cells*/,
-                          double* smin,
-                          double* smax) const
-    {
-        const int np = 2;
-        for (int i = 0; i < n; ++i) {
-            smin[np*i + 0] = sw_[0];
-            smax[np*i + 0] = sw_.back();
-            smin[np*i + 1] = 1.0 - sw_[0];
-            smax[np*i + 1] = 1.0 - sw_.back();
-        }
-    }
-
-private:
-    double krw(double s) const
-    {
-        return Opm::linearInterpolation(sw_, krw_, s);
-    }
-
-    double krw_dsw(double s) const
-    {
-        return Opm::linearInterpolationDerivative(sw_, krw_, s);
-    }
 
 
-    double kro(double s) const
+static void outputWaterCut(const Opm::Watercut& watercut,
-    {
+                           const std::string& output_dir)
-        return Opm::linearInterpolation(so_, kro_, s);
+{
+    // Write water cut curve.
+    std::string fname = output_dir  + "/watercut.txt";
+    std::ofstream os(fname.c_str());
+    if (!os) {
+        THROW("Failed to open " << fname);
     }
-
+    watercut.write(os);
-    double kro_dso(double s) const
+}
-    {
-        return Opm::linearInterpolationDerivative(so_, kro_, s);
-    }
-
-    std::vector<double> sw_;
-    std::vector<double> krw_;
-    std::vector<double> so_;
-    std::vector<double> kro_;
-};
 
 
 // --------------- Types needed to define transport solver ---------------
@@ -301,42 +265,6 @@ typedef Opm::ImplicitTransport<NewtonPolymerTransportModel,
 
 
 
-static void outputState(const UnstructuredGrid& grid,
-                        const Opm::PolymerState& state,
-                        const int step,
-                        const std::string& output_dir)
-{
-    // Write data in VTK format.
-    std::ostringstream vtkfilename;
-    vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
-    std::ofstream vtkfile(vtkfilename.str().c_str());
-    if (!vtkfile) {
-        THROW("Failed to open " << vtkfilename.str());
-    }
-    Opm::DataMap dm;
-    dm["saturation"] = &state.saturation();
-    dm["pressure"] = &state.pressure();
-    dm["concentration"] = &state.concentration();
-    dm["cmax"] = &state.maxconcentration();
-    std::vector<double> cell_velocity;
-    Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
-    dm["velocity"] = &cell_velocity;
-    Opm::writeVtkData(grid, dm, vtkfile);
-
-    // Write data (not grid) in Matlab format
-    for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
-        std::ostringstream fname;
-        fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
-        std::ofstream file(fname.str().c_str());
-        if (!file) {
-            THROW("Failed to open " << fname.str());
-        }
-        const std::vector<double>& d = *(it->second);
-        std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
-    }
-}
-
-
 
 
 class PolymerInflow
@@ -364,19 +292,6 @@ private:
 
 
 
-static void outputWaterCut(const Opm::Watercut& watercut,
-                           const std::string& output_dir)
-{
-    // Write water cut curve.
-    std::string fname = output_dir  + "/watercut.txt";
-    std::ofstream os(fname.c_str());
-    if (!os) {
-        THROW("Failed to open " << fname);
-    }
-    watercut.write(os);
-}
-
-
 // ----------------- Main program -----------------
 int
 main(int argc, char** argv)
@@ -418,7 +333,6 @@ main(int argc, char** argv)
         const int* gc = grid->c_grid()->global_cell;
         std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
         props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
-        // props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
         // Wells init.
         wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
         // Timer init.
@@ -444,7 +358,7 @@ main(int argc, char** argv)
         grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
         // Rock and fluid init.
         // props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        props.reset(new Opm::IncompPropertiesDefaultPolymer(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
         // Wells init.
         wells.reset(new Opm::WellsManager());
         // Timer init.

opm/polymer/IncompPropertiesDefaultPolymer.hpp

@@ -0,0 +1,119 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED
+#define OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED
+
+#include <opm/core/fluid/IncompPropertiesBasic.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/linearInterpolation.hpp>
+#include <vector>
+
+namespace Opm
+{
+
+    class IncompPropertiesDefaultPolymer : public Opm::IncompPropertiesBasic
+    {
+    public:
+        IncompPropertiesDefaultPolymer(const Opm::parameter::ParameterGroup& param, int dim, int num_cells)
+            : Opm::IncompPropertiesBasic(param, dim, num_cells)
+        {
+            ASSERT(numPhases() == 2);
+            sw_.resize(3);
+            sw_[0] = 0.2;
+            sw_[1] = 0.7;
+            sw_[2] = 1.0;
+            krw_.resize(3);
+            krw_[0] = 0.0;
+            krw_[1] = 0.7;
+            krw_[2] = 1.0;
+            so_.resize(2);
+            so_[0] = 0.3;
+            so_[1] = 0.8;
+            kro_.resize(2);
+            kro_[0] = 0.0;
+            kro_[1] = 1.0;
+        }
+
+        virtual void relperm(const int n,
+                             const double* s,
+                             const int* /*cells*/,
+                             double* kr,
+                             double* dkrds) const
+        {
+            // ASSERT(dkrds == 0);
+            // We assume two phases flow
+            for (int i = 0; i < n; ++i) {
+                kr[2*i] = krw(s[2*i]);
+                kr[2*i+1] = kro(s[2*i+1]);
+                if (dkrds != 0) {
+                    dkrds[2*i] = krw_dsw(s[2*i]);
+                    dkrds[2*i+3] = kro_dso(s[2*i+1]);
+                    dkrds[2*i+1] = -dkrds[2*i+3];
+                    dkrds[2*i+2] = -dkrds[2*i];
+                }
+            }
+        }
+
+        virtual void satRange(const int n,
+                              const int* /*cells*/,
+                              double* smin,
+                              double* smax) const
+        {
+            const int np = 2;
+            for (int i = 0; i < n; ++i) {
+                smin[np*i + 0] = sw_[0];
+                smax[np*i + 0] = sw_.back();
+                smin[np*i + 1] = 1.0 - sw_[0];
+                smax[np*i + 1] = 1.0 - sw_.back();
+            }
+        }
+
+    private:
+        double krw(double s) const
+        {
+            return Opm::linearInterpolation(sw_, krw_, s);
+        }
+
+        double krw_dsw(double s) const
+        {
+            return Opm::linearInterpolationDerivative(sw_, krw_, s);
+        }
+
+
+        double kro(double s) const
+        {
+            return Opm::linearInterpolation(so_, kro_, s);
+        }
+
+        double kro_dso(double s) const
+        {
+            return Opm::linearInterpolationDerivative(so_, kro_, s);
+        }
+
+        std::vector<double> sw_;
+        std::vector<double> krw_;
+        std::vector<double> so_;
+        std::vector<double> kro_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED