fix FullyImplicitTwophasePolymersolver constructor problem.

add some input commits for debugging.

examples/sim_2p_fim.cpp

@@ -0,0 +1,262 @@
+/*
+*/
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/core/props/IncompPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/polymer/SimulatorFullyImplicitTwophase.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdInterface.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdBasic.hpp>
+#include <opm/autodiff/polymer/IncompPropsAdFromDeck.hpp>
+
+#include <boost/scoped_ptr.hpp>
+#include <boost/filesystem.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for incompressible two-phase flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    boost::scoped_ptr<EclipseGridParser> deck;
+    boost::scoped_ptr<GridManager> grid;
+    boost::scoped_ptr<IncompPropsAdInterface> props;
+    TwophaseState state;
+    double gravity[3] = { 0.0 };
+    if (use_deck) {
+        std::string deck_filename = param.get<std::string>("deck_filename");
+        deck.reset(new EclipseGridParser(deck_filename));
+        // Grid init
+        grid.reset(new GridManager(*deck));
+        // Rock and fluid init
+        props.reset(new IncompPropsAdFromDeck(*deck, *grid->c_grid()));
+        // Gravity.
+        gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
+        // Init state variables (saturation and pressure).
+        int num_cells = grid->c_grid()->number_of_cells;
+        if (param.has("init_saturation")) {
+            //initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+            const double init_saturation = param.get<double>("init_saturation");
+            for (int c = 0; c < num_cells; ++c) {
+                state.saturation()[2*c] = init_saturation;
+                state.saturation()[2*c+1] = 1. - init_saturation;
+            }
+        } else {
+            if (deck->hasField("PRESSURE")) {
+                // Set saturations from SWAT/SGAS, pressure from PRESSURE.
+                std::vector<double>& s = state.saturation();
+                std::vector<double>& p = state.pressure();
+                const std::vector<double>& p_deck = deck->getFloatingPointValue("PRESSURE");
+                // water-oil or water-gas: we require SWAT
+                if (!deck->hasField("SWAT")) {
+                    OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
+                }
+                const std::vector<double>& sw_deck = deck->getFloatingPointValue("SWAT");
+                for (int c = 0; c < num_cells; ++c) {
+                    int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c];
+                    s[2*c] = sw_deck[c_deck];
+                    s[2*c + 1] = 1.0 - sw_deck[c_deck];
+                    p[c] = p_deck[c_deck];
+                }
+            }
+        }
+    } else {
+        // Grid init.
+        const int nx = param.getDefault("nx", 100);
+        const int ny = param.getDefault("ny", 100);
+        const int nz = param.getDefault("nz", 1);
+        const double dx = param.getDefault("dx", 1.0);
+        const double dy = param.getDefault("dy", 1.0);
+        const double dz = param.getDefault("dz", 1.0);
+        grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
+        // Rock and fluid init.
+        props.reset(new IncompPropsAdBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        // Rock compressibility.
+        // Gravity.
+        gravity[2] = param.getDefault("gravity", 0.0);
+        int num_cells = grid->c_grid()->number_of_cells;
+    }
+
+    // Warn if gravity but no density difference.
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Initialising src
+    std::vector<double> src(num_cells, 0.0);
+    if (use_deck) {
+        // Do nothing, wells will be the driving force, not source terms.
+        if (deck->hasField("SRC")) {
+            const std::vector<double>& src_deck = deck->getFloatingPointValue("SRC");
+            for (int c = 0; c < num_cells; ++c) {
+                int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c];
+                src[c] = src_deck[c_deck];
+            }
+        }
+    } else {
+        // Compute pore volumes, in order to enable specifying injection rate
+        // terms of total pore volume.
+        std::vector<double> porevol;
+        computePorevolume(*grid->c_grid(), props->porosity(), porevol);
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        const double default_injection = use_gravity ? 0.0 : 0.1;
+        const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
+            *tot_porevol_init/unit::day;
+        src[0] = flow_per_sec;
+        src[num_cells - 1] = -flow_per_sec;
+    }
+
+    in: num_cells = grid->c_grid()->number_of_cells;
+
+    // Linear solver.
+    LinearSolverFactory linsolver(param);
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::ofstream epoch_os;
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        std::string filename = output_dir + "/epoch_timing.param";
+        epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        // open file to clean it. The file is appended to in SimulatorTwophase
+        filename = output_dir + "/step_timing.param";
+        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        step_os.close();
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
+              << "                        (number of epochs: "
+              << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush;
+
+    SimulatorReport rep;
+    if (!use_deck) {
+        // Simple simulation without a deck.
+        SimulatorFullyImplicitTwophase simulator(param,
+                                            *grid->c_grid(),
+                                            *props,
+                                            linsolver,
+                                            src);
+        SimulatorTimer simtimer;
+        simtimer.init(param);
+        warnIfUnusedParams(param);
+        rep = simulator.run(simtimer, state, src);
+    } else {
+        // With a deck, we may have more epochs etc.
+        int step = 0;
+        SimulatorTimer simtimer;
+        // Use timer for last epoch to obtain total time.
+        deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
+        simtimer.init(*deck);
+        const double total_time = simtimer.totalTime();
+        for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
+            // Set epoch index.
+            deck->setCurrentEpoch(epoch);
+
+            // Update the timer.
+            if (deck->hasField("TSTEP")) {
+                simtimer.init(*deck);
+            } else {
+                if (epoch != 0) {
+                    OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch);
+                }
+                simtimer.init(param);
+            }
+            simtimer.setCurrentStepNum(step);
+            simtimer.setTotalTime(total_time);
+
+            // Report on start of epoch.
+            std::cout << "\n\n--------------    Starting epoch " << epoch << "    --------------"
+                      << "\n                  (number of steps: "
+                      << simtimer.numSteps() - step << ")\n\n" << std::flush;
+
+
+            // Create and run simulator.
+            SimulatorFullyImplicitTwophase simulator(param,
+                                                *grid->c_grid(),
+                                                *props,
+                                                linsolver,
+                                                src);
+            if (epoch == 0) {
+                warnIfUnusedParams(param);
+            }
+            SimulatorReport epoch_rep = simulator.run(simtimer, state, src);
+            if (output) {
+                epoch_rep.reportParam(epoch_os);
+            }
+            // Update total timing report and remember step number.
+            rep += epoch_rep;
+            step = simtimer.currentStepNum();
+        }
+    }
+
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    rep.report(std::cout);
+
+    if (output) {
+      std::string filename = output_dir + "/walltime.param";
+      std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+      rep.reportParam(tot_os);
+    }
+
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}
+

examples/sim_2p_fincomp_ad.cpp

@@ -0,0 +1,91 @@
+#include "config.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <vector>
+#include <cassert>
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/linalg/LinearSolverUmfpack.hpp>
+#include <opm/core/pressure/FlowBCManager.hpp>
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp>
+#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
+
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+int main (int argc, char** argv)
+try
+{
+    int nx = 20;
+    int ny = 20;
+    int nz = 1;
+    double dx = 10.0;
+    double dy = 10.0;
+    double dz = 10.0;
+    using namespace Opm;
+    parameter::ParameterGroup param(argc, argv, false);
+    GridManager grid_manager(nx, ny, nz, dx, dy, dz);
+    const UnstructuredGrid& grid = *grid_manager.c_grid();
+    int num_cells = grid.number_of_cells;
+    int num_phases = 2;
+    using namespace Opm::unit;
+    using namespace Opm::prefix;
+    std::vector<double> density(num_phases, 1000.0);
+    std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
+    double porosity = 0.5;
+    double permeability = 10.0*milli*darcy;
+    SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
+    IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity,
+                                porosity, permeability, grid.dimensions, num_cells);
+    std::vector<double> omega;
+    std::vector<double> src(num_cells, 0.0);
+    src[0] = 1.;
+    src[num_cells-1] = -1.;
+
+    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", 0.1) * day;
+    const int num_time_steps = param.getDefault("nsteps", 20);
+    std::vector<int> allcells(num_cells);
+    for (int cell = 0; cell < num_cells; ++cell) {
+        allcells[cell] = cell;
+    }
+    TwophaseState state;
+    state.init(grid, 2);
+
+    //initial sat
+    for (int c = 0; c < num_cells; ++c) {
+        state.saturation()[2*c] = 0.2;
+        state.saturation()[2*c+1] = 0.8;
+    }
+    std::vector<double> p(num_cells, 200*Opm::unit::barsa);
+    state.pressure() = p;
+//    state.setFirstSat(allcells, props, TwophaseState::MinSat);
+    std::ostringstream vtkfilename;
+
+    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";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

examples/sim_poly2p_fincomp_ad.cpp

@@ -0,0 +1,146 @@
+#include "config.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <vector>
+#include <cassert>
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/linalg/LinearSolverUmfpack.hpp>
+#include <opm/core/pressure/FlowBCManager.hpp>
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
+#include <opm/polymer/PolymerState.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/io/eclipse/EclipseGridParser.hpp>
+int main (int argc, char** argv)
+try
+{
+    using namespace Opm;
+    parameter::ParameterGroup param(argc, argv, false);
+    bool use_poly_deck = param.has("deck_filename");
+    if (!use_poly_deck) {
+        OPM_THROW(std::runtime_error, "Polymer Properties must be read from deck_filename\n");
+    }
+    std::string deck_filename = param.get<std::string>("deck_filename");
+    EclipseGridParser deck = EclipseGridParser(deck_filename);
+    int nx = param.getDefault("nx", 20);
+    int ny = param.getDefault("ny", 20);
+    int nz = 1;
+    double dx = 2.0;
+    double dy = 2.0;
+    double dz = 0.5;
+    GridManager grid_manager(nx, ny, nz, dx, dy, dz);
+    const UnstructuredGrid& grid = *grid_manager.c_grid();
+    int num_cells = grid.number_of_cells;
+    int num_phases = 2;
+    using namespace Opm::unit;
+    using namespace Opm::prefix;
+    std::vector<double> density(num_phases, 1000.0);
+    std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
+    viscosity[0] = 0.5 * centi * Poise;
+    viscosity[0] = 5 * centi * Poise;
+    double porosity = 0.35;
+    double permeability = 10.0*milli*darcy;
+    SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
+    IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity,
+                                porosity, permeability, grid.dimensions, num_cells);
+
+    // Init polymer properties.
+    // Setting defaults to provide a simple example case.
+        PolymerProperties polymer_props(deck);
+  #if 0
+    if (use_poly_deck) {
+    } else {
+        double c_max = param.getDefault("c_max_limit", 5.0);
+        double mix_param = param.getDefault("mix_param", 1.0);
+        double rock_density = param.getDefault("rock_density", 1000.0);
+        double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15);
+        double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability
+        double c_max_ads = param.getDefault("c_max_ads", 1.);
+        int ads_index = param.getDefault<int>("ads_index", Opm::PolymerProperties::NoDesorption);
+        std::vector<double> c_vals_visc(2, -1e100);
+        c_vals_visc[0] = 0.0;
+        c_vals_visc[1] = 7.0;
+        std::vector<double> visc_mult_vals(2, -1e100);
+        visc_mult_vals[0] = 1.0;
+        // poly_props.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
+        visc_mult_vals[1] = 20.0;
+        std::vector<double> c_vals_ads(3, -1e100);
+        c_vals_ads[0] = 0.0;
+        c_vals_ads[1] = 2.0;
+        c_vals_ads[2] = 8.0;
+        std::vector<double> ads_vals(3, -1e100);
+        ads_vals[0] = 0.0;
+        ads_vals[1] = 0.0015;
+        ads_vals[2] = 0.0025;
+        PolymerProperties polymer_props;
+        polymer_props.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
+                         static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
+                           c_vals_visc,  visc_mult_vals, c_vals_ads, ads_vals);
+    }
+    #endif
+    PolymerPropsAd polymer_props_ad(polymer_props);
+    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] = -1. / day;
+
+    PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 30.0)*Opm::unit::day,
+                                      param.getDefault("poly_end_days", 80.0)*Opm::unit::day,
+                                      param.getDefault("poly_amount", polymer_props.cMax()));
+    FlowBCManager bcs;
+    LinearSolverUmfpack linsolver;
+    FullyImplicitTwophasePolymerSolver solver(grid, props,polymer_props_ad, linsolver);
+    std::vector<double> porevol;
+    Opm::computePorevolume(grid, props.porosity(), porevol);
+    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) {
+        allcells[cell] = cell;
+    }
+    PolymerState state;
+    state.init(grid, 2);
+    //initial sat
+    for (int c = 0; c < num_cells; ++c) {
+        state.saturation()[2*c] = 0.2;
+        state.saturation()[2*c+1] = 0.8;
+    }
+    std::vector<double> p(num_cells, 200*Opm::unit::barsa);
+    state.pressure() = p;
+
+    std::vector<double> c(num_cells, 0.0);
+    state.concentration() = c;
+    std::ostringstream vtkfilename;
+    double currentime = 0;
+    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";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp

@@ -64,11 +64,11 @@ typedef Eigen::Array<double,
                                 const LinearSolverInterface&    linsolver)
         : grid_ (grid)
         , fluid_(fluid)
-        , polymer_props_ad_ (polymer_props_ad_)
+        , polymer_props_ad_ (polymer_props_ad)
         , linsolver_(linsolver)
         , cells_ (buildAllCells(grid.number_of_cells))
         , ops_(grid)
-        , residual_(std::vector<ADB>(fluid.numPhases() + 1, ADB::null()))
+        , residual_(std::vector<ADB>(3, ADB::null()))
      {
      }
 
@@ -82,7 +82,6 @@ typedef Eigen::Array<double,
          PolymerState& x,
          const std::vector<double>& src,
          const std::vector<double>& polymer_inflow)
-    //     const bool if_polymer_actived)
     {
         
         V pvol(grid_.number_of_cells);
@@ -100,6 +99,10 @@ typedef Eigen::Array<double,
         const double rtol  = 5.0e-8;
         const int    maxit = 15;
 
+        std::cout << "Primary variables:\n";
+        std::cout << "pressure\n"<<old_state.pressure<< std::endl;
+        std::cout << "saturation\n"<<old_state.saturation[0] << std::endl;
+        std::cout << "concentration\n"<<old_state.concentration << std::endl;
         assemble(pvdt, old_state, x, src, polymer_inflow);//, if_polymer_actived);
 
         const double r0  = residualNorm();
@@ -162,11 +165,9 @@ typedef Eigen::Array<double,
         const DataBlock s_all = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
         for (int phase = 0; phase < np; ++phase) {
             state.saturation[phase] = ADB::constant(s_all.col(phase));
-   //     state.saturation[1] = ADB::constant(s_all.col(1));
         }
 
         // Concentration
-
         assert(not x.concentration().empty());
         const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
 
@@ -256,29 +257,28 @@ typedef Eigen::Array<double,
                 + ops_.div*mflux - source;
         }
 
-        double amount = PolymerInjectedAmount(polymer_inflow);
+     //   bool use_polymer = 1;//(amount != 0);
-        bool use_polymer = (amount != 0);
+   //     const int nc = grid_.number_of_cells;
-        std::cout << "\n             Polymer Amount\n"
+//        if (use_polymer) {
-                  << std::setprecision(9)
-                  << std::setw(18) << amount << std::endl;
-        const int nc = grid_.number_of_cells;
-        if (use_polymer) {
             // Mass balance equation for polymer
-            const V src_polymer = Eigen::Map<const V>(&polymer_inflow[0], nc);
+            const ADB src_polymer = polymerSource(kr ,src, polymer_inflow, state);  
             ADB mc = computeMc(state);
             ADB mflux  = computeMassFlux(0, trans, kr, state);
             residual_[2] = pvdt * (state.saturation[0] * state.concentration 
                         - old_state.saturation[0] * old_state.concentration) 
                         + ops_.div * state.concentration * mc * mflux - src_polymer;
-        } else {
+  //      } else {
-            residual_[2] = ADB::constant(V::Zero(nc));
+    //        residual_[2] = ADB::constant(V::Zero(nc));
-        }
+//        }
+        for (int i = 0; i < 3; ++i)
+        std::cout<<"residual_["<<i<<"]\n"<<residual_[i] << std::endl;
+
     }
    
 
     double
     FullyImplicitTwophasePolymerSolver::
-    PolymerInjectedAmount(const std::vector<double>& polymer_inflow) const
+    polymerInjectedAmount(const std::vector<double>& polymer_inflow) const
     {
         double amount = 0;
         for (int i = 0; i < int(polymer_inflow.size()); ++i) {
@@ -292,7 +292,8 @@ typedef Eigen::Array<double,
     ADB 
     FullyImplicitTwophasePolymerSolver::accumSource(const int phase,
                                                  const std::vector<ADB>& kr,
-                                                 const std::vector<double>& src) const
+                                                 const std::vector<double>& src
+                                                 ) const
     {
         //extract the source to out and in source.
         std::vector<double> outsrc;
@@ -313,7 +314,6 @@ typedef Eigen::Array<double,
         const V source = Eigen::Map<const V>(& src[0], grid_.number_of_cells);
         const V outSrc = Eigen::Map<const V>(& outsrc[0], grid_.number_of_cells);
         const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells);
-        
         // compute the out-fracflow.
         ADB f_out = computeFracFlow(phase, kr);
         // compute the in-fracflow.
@@ -323,11 +323,44 @@ typedef Eigen::Array<double,
         } else if (phase == 0) {
             f_in = V::Ones(grid_.number_of_cells);
         }
-        return f_out * outSrc + f_in * inSrc;
+        return f_out * outSrc + f_in * inSrc ;
      }
 
 
-
+    ADB 
+    FullyImplicitTwophasePolymerSolver::
+    polymerSource(
+                  const std::vector<ADB>& kr,
+                  const std::vector<double>& src,
+                  const std::vector<double>& polymer_inflow_c,
+                  const SolutionState& state) const
+    {
+        //extract the source to out and in source.
+        std::vector<double> outsrc;
+        std::vector<double> insrc;
+        std::vector<double>::const_iterator it;
+        for (it = src.begin(); it != src.end(); ++it) {
+            if (*it < 0) {
+                outsrc.push_back(*it);
+                insrc.push_back(0.0);
+            } else if (*it > 0) {
+                insrc.push_back(*it);
+                outsrc.push_back(0.0);
+            } else {
+                outsrc.emplace_back(0);
+                insrc.emplace_back(0);
+            }
+        }
+        const V source = Eigen::Map<const V>(& src[0], grid_.number_of_cells);
+        const V outSrc = Eigen::Map<const V>(& outsrc[0], grid_.number_of_cells);
+        const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells);
+        const V polyin = Eigen::Map<const V>(& polymer_inflow_c[0], grid_.number_of_cells);
+        // compute the out-fracflow.
+        ADB f_out = computeFracFlow(0, kr);
+        // compute the in-fracflow.
+        V f_in = V::Ones(grid_.number_of_cells);
+        return f_out * outSrc * state.concentration + f_in * inSrc * polyin ;
+     }
 
 
     ADB
@@ -355,8 +388,8 @@ typedef Eigen::Array<double,
     	if (np != 2) {
 	        OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwophasePolymerSolver.");
 	    }
-        ADB mass_res = vertcat(residual_[0], residual_[1]);
+        ADB mass_phase_res = vertcat(residual_[0], residual_[1]);
-	    mass_res = collapseJacs(vertcat(mass_res, residual_[2]));
+	    ADB mass_res = collapseJacs(vertcat(mass_phase_res, residual_[2]));
 
         const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = mass_res.derivative()[0];
         V dx(V::Zero(mass_res.size()));
@@ -391,7 +424,6 @@ typedef Eigen::Array<double,
         int varstart = nc;
         const V dsw = subset(dx, Span(nc, 1, varstart));
         varstart += dsw.size();
-        std::cout << dx.size() << " " << varstart << std::endl;
         const V dc = subset(dx, Span(nc, 1, varstart));
         varstart += dc.size();
 

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp

@@ -87,7 +87,11 @@ namespace Opm {
                         const SolutionState&    state) const;
         double
         residualNorm() const;
-        
+        ADB
+        polymerSource(const std::vector<ADB>& kr,
+                      const std::vector<double>& src,
+                      const std::vector<double>& polymer_inflow_c,
+                      const SolutionState& state) const;
 
         ADB 
         computeMc(const SolutionState& state) const;
@@ -100,7 +104,7 @@ namespace Opm {
         ADB
         transMult(const ADB& p) const;
         double
-        PolymerInjectedAmount(const std::vector<double>& polymer_inflow) const;
+        polymerInjectedAmount(const std::vector<double>& polymer_inflow) const;
     };
 } // namespace Opm
 #endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -65,7 +65,6 @@ namespace Opm {
 */
         typedef AutoDiffBlock<double> ADB;
         typedef ADB::V V;
-
         PolymerPropsAd(const PolymerProperties& polymer_props);
 
         ~PolymerPropsAd();