From 11ea882ed8067301d0c52d717b9103add7e4e355 Mon Sep 17 00:00:00 2001
From: Liu Ming <miliu@statoil.com>
Date: Wed, 11 Dec 2013 18:27:05 +0800
Subject: [PATCH] remove some output commits. get values from command line for
 exapmes.

---
 examples/sim_2p_fincomp_ad.cpp                | 33 +++++++++++-------
 examples/sim_poly2p_fincomp_ad.cpp            | 20 ++++++++---
 .../FullyImplicitTwoPhaseSolver.cpp           |  3 +-
 .../FullyImplicitTwophasePolymerSolver.cpp    | 34 +++++++++----------
 4 files changed, 54 insertions(+), 36 deletions(-)

diff --git a/examples/sim_2p_fincomp_ad.cpp b/examples/sim_2p_fincomp_ad.cpp
index 1b7f1679f..17f75425b 100644
--- a/examples/sim_2p_fincomp_ad.cpp
+++ b/examples/sim_2p_fincomp_ad.cpp
@@ -24,14 +24,14 @@
 int main (int argc, char** argv)
 try
 {
-    int nx = 30;
-    int ny = 30;
-    int nz = 1;
-    double dx = 2.0;
-    double dy = 2.0;
-    double dz = 0.5;
     using namespace Opm;
     parameter::ParameterGroup param(argc, argv, false);
+    int nx = param.getDefault("nx", 30);
+    int ny = param.getDefault("ny", 1);
+    int nz = 1;
+    double dx = 10./nx;
+    double dy = 1.0;
+    double dz = 1.0;
     GridManager grid_manager(nx, ny, nz, dx, dy, dz);
     const UnstructuredGrid& grid = *grid_manager.c_grid();
     int num_cells = grid.number_of_cells;
@@ -49,15 +49,15 @@ try
                                 porosity, permeability, grid.dimensions, num_cells);
     std::vector<double> omega;
     std::vector<double> src(num_cells, 0.0);
-    src[0] = 10. / day;
-    src[num_cells-1] = -10. / day;
+    src[0] = 1. / day;
+    src[num_cells-1] = -1. / day;
 
     FlowBCManager bcs;
     LinearSolverUmfpack linsolver;
     FullyImplicitTwoPhaseSolver solver(grid, props, linsolver);
     std::vector<double> porevol;
     Opm::computePorevolume(grid, props.porosity(), porevol);
-    const double dt = param.getDefault("dt", 10) * day;
+    const double dt = param.getDefault("dt", 10.) * day;
     const int num_time_steps = param.getDefault("nsteps", 10);
     std::vector<int> allcells(num_cells);
     for (int cell = 0; cell < num_cells; ++cell) {
@@ -68,18 +68,25 @@ try
 
     //initial sat
     for (int c = 0; c < num_cells; ++c) {
-        state.saturation()[2*c] = 0;
-        state.saturation()[2*c+1] = 1;
+        state.saturation()[2*c] = 0.2;
+        state.saturation()[2*c+1] = 0.8;
     }
     std::vector<double> p(num_cells, 100*Opm::unit::barsa);
     state.pressure() = p;
-//    state.setFirstSat(allcells, props, TwophaseState::MinSat);
     std::ostringstream vtkfilename;
 
+    // Write the initial state.
+    vtkfilename.str("");
+    vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << 0 << ".vtu";
+    std::ofstream vtkfile(vtkfilename.str().c_str());
+    Opm::DataMap dm;
+    dm["saturation"] = &state.saturation();
+    dm["pressure"] = &state.pressure();
+    Opm::writeVtkData(grid, dm, vtkfile);
     for (int i = 0; i < num_time_steps; ++i) {
         solver.step(dt, state, src);
         vtkfilename.str("");
-        vtkfilename << "sim_2p_fincomp-" << std::setw(3) << std::setfill('0') << i << ".vtu";
+        vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << i + 1 << ".vtu";
         std::ofstream vtkfile(vtkfilename.str().c_str());
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
diff --git a/examples/sim_poly2p_fincomp_ad.cpp b/examples/sim_poly2p_fincomp_ad.cpp
index f1262e1f9..68efdb751 100644
--- a/examples/sim_poly2p_fincomp_ad.cpp
+++ b/examples/sim_poly2p_fincomp_ad.cpp
@@ -39,7 +39,7 @@ try
     int nx = param.getDefault("nx", 30);
     int ny = param.getDefault("ny", 30);
     int nz = 1;
-    double dx = 10.0;
+    double dx = 10./nx;
     double dy = 1.0;
     double dz = 1.0;
     GridManager grid_manager(nx, ny, nz, dx, dy, dz);
@@ -60,7 +60,7 @@ try
 
     // Init polymer properties.
     // Setting defaults to provide a simple example case.
-        PolymerProperties polymer_props(deck);
+    PolymerProperties polymer_props(deck);
   #if 0
     if (use_poly_deck) {
     } else {
@@ -96,8 +96,8 @@ try
     std::vector<double> omega;
     std::vector<double> src(num_cells, 0.0);
     std::vector<double> src_polymer(num_cells);
-    src[0] = 10. / day;
-    src[num_cells-1] = -10. / day;
+    src[0] = param.getDefault("insrc", 1.) / day;
+    src[num_cells-1] = -param.getDefault("insrc", 1.) / day;
 
     PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
                                       param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
@@ -127,12 +127,22 @@ try
     state.concentration() = c;
     std::ostringstream vtkfilename;
     double currentime = 0;
+
+    // Write the initial state.
+    vtkfilename.str("");
+    vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << 0<< ".vtu";
+    std::ofstream vtkfile(vtkfilename.str().c_str());
+    Opm::DataMap dm;
+    dm["saturation"] = &state.saturation();
+    dm["pressure"] = &state.pressure();
+    dm["concentration"] = &state.concentration();
+    Opm::writeVtkData(grid, dm, vtkfile);
     for (int i = 0; i < num_time_steps; ++i) {
         currentime += dt;
         polymer_inflow.getInflowValues(currentime, currentime+dt, src_polymer);
         solver.step(dt, state, src, src_polymer);
         vtkfilename.str("");
-        vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i << ".vtu";
+        vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i + 1<< ".vtu";
         std::ofstream vtkfile(vtkfilename.str().c_str());
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp
index 60e7a11be..cf03d56f3 100644
--- a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp
+++ b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp
@@ -396,7 +396,8 @@ typedef Eigen::Array<double,
 //        const ADB tr_mult = transMult(state.pressure);
         const double* mus = fluid_.viscosity();
         ADB  mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
-
+        if (phase ==0)
+            std::cout << "watetr mob\n" << mob.value() << std::endl;
         const ADB dp = ops_.ngrad * state.pressure;
         const ADB head = trans * dp;
 
diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp
index f6ba27bcb..5e16ce992 100644
--- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp
+++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp
@@ -165,8 +165,8 @@ typedef Eigen::Array<double,
         // Concentration
         assert(not x.concentration().empty());
         const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
-
         state.concentration = ADB::constant(c);
+
         return state;
     }
 
@@ -246,25 +246,16 @@ typedef Eigen::Array<double,
         for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
             const ADB mflux = computeMassFlux(phase, trans, kr, state);
             ADB source = accumSource(phase, kr, src);
+ //           std::cout << "phase-"<<phase<<"-source\n"<< source.value() << std::endl;
             residual_[phase] =
                 pvdt*(state.saturation[phase] - old_state.saturation[phase])
                 + ops_.div*mflux - source;
         }
       // Mass balance equation for polymer
         const ADB src_polymer = polymerSource(kr, src, polymer_inflow, state);
-  //      std::cout << "polymer src\n";
-   //     std::cout << src_polymer << std::endl;
- //       const ADB src_polymer =ADB::function(V::Zero(grid_.number_of_cells), state.concentration.derivative());
         ADB mc = computeMc(state);
-        
         ADB poly_mflux  = computePolymerMassFlux(trans, mc, kr, state);
-#if 0
-        std::cout << "polymer mass\n";
-        std::cout << pvdt * (state.saturation[0] * state.concentration 
-                               - old_state.saturation[0] * old_state.concentration) <<std::endl;
-        std::cout << "polymer mass flux\n";
-        std::cout << poly_mflux << std::endl;
-#endif     
+//      std::cout << "polymer source \n" << src_polymer.value() << std::endl;
         residual_[2] = pvdt * (state.saturation[0] * state.concentration 
                       - old_state.saturation[0] * old_state.concentration) 
                       + ops_.div * poly_mflux - src_polymer;
@@ -418,7 +409,6 @@ typedef Eigen::Array<double,
         const V p = p_old - dp;
         std::copy(&p[0], &p[0] + nc, state.pressure().begin());
 
-
         // Saturation updates.
         const double dsmax = 0.3;
         const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
@@ -472,8 +462,18 @@ typedef Eigen::Array<double,
     {
 //        const ADB tr_mult = transMult(state.pressure);
         const double* mus = fluid_.viscosity();
-        ADB  mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
-
+        ADB mob = ADB::null();
+        if (phase == 0) {
+            ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
+//            for (int i = 0; i < grid_.number_of_cells; ++i) {
+//             std::cout << state.concentration.value()(i) << "  " << 1./inv_wat_eff_vis.value()(i) << std::endl;
+//            std::cout << "water effective vis\n"<<1./inv_wat_eff_v1./inv_wat_eff_vis.value()is.value()<<std::endl;
+//            }
+            mob = kr[0] * inv_wat_eff_vis;
+//            std::cout << "watetr mob\n" << mob.value() << std::endl;
+        } else if (phase == 1) {
+            mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
+        }
         const ADB dp = ops_.ngrad * state.pressure;
         const ADB head = trans * dp;
 
@@ -490,8 +490,8 @@ typedef Eigen::Array<double,
 
    {
         const double* mus = fluid_.viscosity();
-        ADB water_mob = kr[0] / V::Constant(kr[0].size(), 1, mus[0]);
-        ADB poly_mob = state.concentration * mc * water_mob;
+        ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
+        ADB poly_mob =  mc * kr[0] * inv_wat_eff_vis;
         
         const ADB dp = ops_.ngrad * state.pressure;
         const ADB head = trans * dp;