/*
  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/>.
*/


#include <config.h>

#include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/grid.h>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/common/ErrorMacros.hpp>
#include <cmath>
#include <list>
#include <iostream>
// Choose error policy for scalar solves here.
typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> RootFinder;


class Opm::TransportSolverTwophasePolymer::ResidualEquation
{
public:
    int cell;
    double s0;
    double c0;
    double cmax0;
    double influx;  // sum_j min(v_ij, 0)*f(s_j)
    double influx_polymer;  // sum_j min(v_ij, 0)*f(s_j)*mc(c_j)
    double outflux; // sum_j max(v_ij, 0)
    double comp_term; //q - sum_j v_ij
    double porosity;
    double dtpv;    // dt/pv(i)
    double dps;
    double rhor;
    double ads0;
    GradientMethod gradient_method;

    TransportSolverTwophasePolymer& tm;

    ResidualEquation(TransportSolverTwophasePolymer& tmodel, int cell_index);
    void computeResidual(const double* x, double* res) const;
    void computeResidual(const double* x, double* res, double& mc, double& ff) const;
    double computeResidualS(const double* x) const;
    double computeResidualC(const double* x) const;
    void computeGradientResS(const double* x, double* res, double* gradient) const;
    void computeGradientResC(const double* x, double* res, double* gradient) const;
    void computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const;

private:
    void computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c,
                             const bool if_dres_s_dsdc, const bool if_dres_c_dsdc,
                             double* res, double* dres_s_dsdc,
                             double* dres_c_dsdc, double& mc, double& ff) const;
};

class Opm::TransportSolverTwophasePolymer::ResidualCGrav {
public:
    const TransportSolverTwophasePolymer& tm;
    const int cell;
    const double s0;
    const double c0;
    const double cmax0;
    const double porosity;
    const double dtpv;    // dt/pv(i)
    const double dps;
    const double rhor;
    double c_ads0;
    double gf[2];
    int nbcell[2];
    mutable double last_s;

    ResidualCGrav(const TransportSolverTwophasePolymer& tmodel,
                  const std::vector<int>& cells,
                  const int pos,
                  const double* gravflux);

    double operator()(double c) const;
    double computeGravResidualS(double s, double c) const;
    double computeGravResidualC(double s, double c) const;
    double lastSaturation() const;
};

class Opm::TransportSolverTwophasePolymer::ResidualSGrav {
public:
    const ResidualCGrav& res_c_eq_;
    double c;

    ResidualSGrav(const ResidualCGrav& res_c_eq, const double c_init = 0.0);
    double operator()(double s) const;
};


namespace
{
    bool check_interval(const double* xmin, const double* xmax, double* x);

    double norm(double* res)
    {
	return std::max(std::abs(res[0]), std::abs(res[1]));
    }

    // Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual.
    // The curve starts at "x", goes along the direction "direction" until it hits the boundary of the box of
    // admissible values for "s" and "x" (which is given by "[x_min[0], x_max[0]]x[x_min[1], x_max[1]]").
    // Then it joins in a straight line the point "end_point".
    class CurveInSCPlane{
    public:
	CurveInSCPlane();
	void setup(const double* x, const double* direction,
		   const double* end_point, const double* x_min,
		   const double* x_max, const double tol,
                   double& t_max_out, double& t_out_out);
	void computeXOfT(double*, const double) const;

    private:
	double direction_[2];
	double end_point_[2];
	double x_max_[2];
	double x_min_[2];
	double t_out_;
	double t_max_; // t_max = t_out + 1
	double x_out_[2];
	double x_[2];
    };


    // Compute the "s" residual along the curve "curve" for a given residual equation "res_eq".
    // The operator() is sent to a root solver.
    class ResSOnCurve
    {
    public:
	ResSOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq);
	double operator()(const double t) const;
	CurveInSCPlane curve;
    private:
	const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
    };

    // Compute the "c" residual along the curve "curve" for a given residual equation "res_eq".
    // The operator() is sent to a root solver.
    class ResCOnCurve
    {
    public:
	ResCOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq);
	double operator()(const double t) const;
	CurveInSCPlane curve;
    private:
	const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
    };

}


namespace Opm
{
    TransportSolverTwophasePolymer::TransportSolverTwophasePolymer(const UnstructuredGrid& grid,
						 const IncompPropertiesInterface& props,
						 const PolymerProperties& polyprops,
						 const SingleCellMethod method,
						 const double tol,
						 const int maxit)
	: grid_(grid),
	  porosity_(props.porosity()),
	  porevolume_(NULL),
	  props_(props),
	  polyprops_(polyprops),
	  tol_(tol),
	  maxit_(maxit),
	  darcyflux_(0),
	  source_(0),
	  polymer_inflow_c_(0),
	  dt_(0.0),
	  concentration_(0),
	  cmax_(0),
	  fractionalflow_(grid.number_of_cells, -1.0),
	  mc_(grid.number_of_cells, -1.0),
	  method_(method),
	  adhoc_safety_(1.1) 
    {
	if (props.numPhases() != 2) {
	    OPM_THROW(std::runtime_error, "Property object must have 2 phases");
	}
	visc_ = props.viscosity();

#ifdef PROFILING
        res_counts.clear();
#endif

	// Set up smin_ and smax_
	int num_cells = props.numCells();
	smin_.resize(props.numPhases()*num_cells);
	smax_.resize(props.numPhases()*num_cells);
	std::vector<int> cells(num_cells);
	for (int i = 0; i < num_cells; ++i) {
	    cells[i] = i;
	}
	props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]);
    }




    void TransportSolverTwophasePolymer::setPreferredMethod(SingleCellMethod method)
    {
        method_ = method;
    }




    void TransportSolverTwophasePolymer::solve(const double* darcyflux,
                                      const double* porevolume,
				      const double* source,
				      const double* polymer_inflow_c,
				      const double dt,
				      std::vector<double>& saturation,
				      std::vector<double>& concentration,
				      std::vector<double>& cmax)
    {
	darcyflux_ = darcyflux;
        porevolume_ = porevolume;
	source_ = source;
	polymer_inflow_c_ = polymer_inflow_c;
	dt_ = dt;
        toWaterSat(saturation, saturation_);
	concentration_ = &concentration[0];
	cmax_ = &cmax[0];
#if PROFILING
        res_counts.clear();
#endif
        reorderAndTransport(grid_, darcyflux);
        toBothSat(saturation_, saturation);
    }




    // Residual for saturation equation, single-cell implicit Euler transport
    //
    //     r(s) = s - s0 + dt/pv*( influx + outflux*f(s) )
    //
    // where influx is water influx, outflux is total outflux.
    // Influxes are negative, outfluxes positive.
    struct TransportSolverTwophasePolymer::ResidualS
    {
        TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
	const double c_;
	explicit ResidualS(TransportSolverTwophasePolymer::ResidualEquation& res_eq,
			   const double c)
	    : res_eq_(res_eq),
	      c_(c)
	{
	}

	double operator()(double s) const
	{
            double x[2];
            x[0] = s;
            x[1] = c_;
            return res_eq_.computeResidualS(x);
	}
    };

    // Residual for concentration equation, single-cell implicit Euler transport
    //
    //  \TODO doc me
    // where ...
    // Influxes are negative, outfluxes positive.
    struct TransportSolverTwophasePolymer::ResidualC
    {
	mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value.
        TransportSolverTwophasePolymer::ResidualEquation& res_eq_;
	explicit ResidualC(TransportSolverTwophasePolymer::ResidualEquation& res_eq)
	    : res_eq_(res_eq)
	{}

	void computeBothResiduals(const double s_arg, const double c_arg, double& res_s, double& res_c, double& mc, double& ff) const
	{
            double x[2];
            double res[2];
            x[0] = s_arg;
            x[1] = c_arg;
            res_eq_.computeResidual(x, res, mc, ff);
            res_s = res[0];
            res_c = res[1];
	}

	double operator()(double c) const
	{
	    ResidualS res_s(res_eq_, c);
	    int iters_used;
	    // Solve for s first.
	    // s = modifiedRegulaFalsi(res_s, std::max(tm.smin_[2*cell], dps), tm.smax_[2*cell],
	    //     		    tm.maxit_, tm.tol_, iters_used);
	    s = RootFinder::solve(res_s, res_eq_.s0, 0.0, 1.0,
                                  res_eq_.tm.maxit_, res_eq_.tm.tol_, iters_used);
            double x[2];
            x[0] = s;
            x[1] = c;
            double res = res_eq_.computeResidualC(x);
#ifdef EXTRA_DEBUG_OUTPUT
	    std::cout << "c = " << c << "    s = " << s << "    c-residual = " << res << std::endl;
#endif
 	    return res;
	}

	double lastSaturation() const
	{
	    return s;
	}
    };


    // ResidualEquation gathers parameters to construct the residual, computes its
    // value and the values of its derivatives.

    TransportSolverTwophasePolymer::ResidualEquation::ResidualEquation(TransportSolverTwophasePolymer& tmodel, int cell_index)
	: tm(tmodel)
    {
	gradient_method = Analytic;
	cell    = cell_index;
	s0      = tm.saturation_[cell];
	c0      = tm.concentration_[cell];
	cmax0   = tm.cmax_[cell];
	dps = tm.polyprops_.deadPoreVol();
	rhor = tm.polyprops_.rockDensity();
        tm.polyprops_.adsorption(c0, cmax0, ads0);
	double dflux       = -tm.source_[cell];
	bool src_is_inflow = dflux < 0.0;
	influx  =  src_is_inflow ? dflux : 0.0;
        double mc;
        tm.computeMc(tm.polymer_inflow_c_[cell_index], mc);
	influx_polymer = src_is_inflow ? dflux*mc : 0.0;
	outflux = !src_is_inflow ? dflux : 0.0;
	comp_term = tm.source_[cell];   // Note: this assumes that all source flux is water.
	dtpv    = tm.dt_/tm.porevolume_[cell];
	porosity = tm.porosity_[cell];
	for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
	    int f = tm.grid_.cell_faces[i];
	    double flux;
	    int other;
	    // Compute cell flux
	    if (cell == tm.grid_.face_cells[2*f]) {
		flux  = tm.darcyflux_[f];
		other = tm.grid_.face_cells[2*f+1];
	    } else {
		flux  =-tm.darcyflux_[f];
		other = tm.grid_.face_cells[2*f];
	    }
	    // Add flux to influx or outflux, if interior.
	    if (other != -1) {
		if (flux < 0.0) {
		    influx  += flux*tm.fractionalflow_[other];
		    influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other];
		} else {
		    outflux += flux;
		}
		comp_term -= flux;
	    }
	}
    }


    void TransportSolverTwophasePolymer::ResidualEquation::computeResidual(const double* x, double* res) const
    {
        double dres_s_dsdc[2];
        double dres_c_dsdc[2];
        double mc;
        double ff;
        computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
    }

    void TransportSolverTwophasePolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const
    {
        double dres_s_dsdc[2];
        double dres_c_dsdc[2];
        computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
    }


    double TransportSolverTwophasePolymer::ResidualEquation::computeResidualS(const double* x) const
    {
        double res[2];
        double dres_s_dsdc[2];
        double dres_c_dsdc[2];
        double mc;
        double ff;
        computeResAndJacobi(x, true, false, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
        return res[0];
    }

    double TransportSolverTwophasePolymer::ResidualEquation::computeResidualC(const double* x) const
    {
        double res[2];
        double dres_s_dsdc[2];
        double dres_c_dsdc[2];
        double mc;
        double ff;
        computeResAndJacobi(x, false, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
        return res[1];
    }

    void TransportSolverTwophasePolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const
    // If gradient_method == FinDif, use finite difference
    // If gradient_method == Analytic, use analytic expresions
    {
        double dres_c_dsdc[2];
        double mc;
        double ff;
        computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff);
    }

    void TransportSolverTwophasePolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const
    // If gradient_method == FinDif, use finite difference
    // If gradient_method == Analytic, use analytic expresions
    {
        double dres_s_dsdc[2];
        double mc;
        double ff;
        computeResAndJacobi(x, true, true, false, true, res, dres_s_dsdc, gradient, mc, ff);
    }

    // Compute the Jacobian of the residual equations.
    void TransportSolverTwophasePolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const
    {
        double res[2];
        double mc;
        double ff;
        computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff);
    }

    void TransportSolverTwophasePolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c,
                                                                      const bool if_dres_s_dsdc, const bool if_dres_c_dsdc,
                                                                      double* res, double* dres_s_dsdc,
                                                                      double* dres_c_dsdc, double& mc, double& ff) const
    {
        if ((if_dres_s_dsdc || if_dres_c_dsdc) && gradient_method == Analytic) {
            double s = x[0];
            double c = x[1];
            double  dff_dsdc[2];
            double mc_dc;
            double ads_dc;
            double ads;
            tm.fracFlowWithDer(s, c, cmax0, cell, ff, dff_dsdc);
            if (if_dres_c_dsdc) {
                tm.polyprops_.adsorptionWithDer(c, cmax0, ads, ads_dc);
                tm.computeMcWithDer(c, mc, mc_dc);
            } else {
                tm.polyprops_.adsorption(c, cmax0, ads);
                tm.computeMc(c, mc);
            }
            if (if_res_s) {
                res[0] = s - s0 +  dtpv*(outflux*ff + influx + s*comp_term);
#if PROFILING
                tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1]));
#endif
            }
            if (if_res_c) {
                res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
                    + rhor*((1.0 - porosity)/porosity)*(ads - ads0)
                    + dtpv*(outflux*ff*mc + influx_polymer)
		    + dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term;
#if PROFILING
                tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1]));
#endif
            }
            if (if_dres_s_dsdc) {
                dres_s_dsdc[0] = 1 + dtpv*(outflux*dff_dsdc[0] + comp_term);
                dres_s_dsdc[1] = dtpv*outflux*dff_dsdc[1];
            }
            if (if_dres_c_dsdc) {
                dres_c_dsdc[0] = (1.0 - dps)*c + dtpv*outflux*(dff_dsdc[0])*mc
		    + dtpv*c*(1.0 - dps)*comp_term;
                dres_c_dsdc[1] = (1 - dps)*s + rhor*((1.0 - porosity)/porosity)*ads_dc
                    + dtpv*outflux*(dff_dsdc[1]*mc + ff*mc_dc)
		    + dtpv*(s*(1.0 - dps) - rhor*ads_dc)*comp_term;
            }

        } else if (if_res_c || if_res_s) {
            double s = x[0];
            double c = x[1];
            tm.fracFlow(s, c, cmax0, cell, ff);
            if (if_res_s) {
                res[0] = s - s0 +  dtpv*(outflux*ff + influx + s*comp_term);
#if PROFILING
                tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1]));
#endif
            }
            if (if_res_c) {
                tm.computeMc(c, mc);
                double ads;
                tm.polyprops_.adsorption(c, cmax0, ads);
                res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0
                    + rhor*((1.0 - porosity)/porosity)*(ads - ads0)
                    + dtpv*(outflux*ff*mc + influx_polymer)
		    + dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term;
#if PROFILING
                tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1]));
#endif
            }
        }

	if ((if_dres_c_dsdc || if_dres_s_dsdc) && gradient_method == FinDif) {
	    double epsi = 1e-8;
	    double res_epsi[2];
            double res_0[2];
	    double x_epsi[2];
	    computeResidual(x, res_0);
            if (if_dres_s_dsdc) {
                x_epsi[0] = x[0] + epsi;
                x_epsi[1] = x[1];
                computeResidual(x_epsi, res_epsi);
                dres_s_dsdc[0] =  (res_epsi[0] - res_0[0])/epsi;
                x_epsi[0] = x[0];
                x_epsi[1] = x[1] + epsi;
                computeResidual(x_epsi, res_epsi);
                dres_s_dsdc[1] = (res_epsi[0] - res_0[0])/epsi;
            }
            if (if_dres_c_dsdc) {
                x_epsi[0] = x[0] + epsi;
                x_epsi[1] = x[1];
                computeResidual(x_epsi, res_epsi);
                dres_c_dsdc[0] =  (res_epsi[1] - res_0[1])/epsi;
                x_epsi[0] = x[0];
                x_epsi[1] = x[1] + epsi;
                computeResidual(x_epsi, res_epsi);
                dres_c_dsdc[1] = (res_epsi[1] - res_0[1])/epsi;
            }
        }
    }

    void TransportSolverTwophasePolymer::solveSingleCell(const int cell)
    {
	switch (method_) {
	case Bracketing:
	    solveSingleCellBracketing(cell);
	    break;
	case Newton:
	    solveSingleCellNewton(cell);
	    break;
	case Gradient:
	    solveSingleCellGradient(cell);
	    break;
	case NewtonSimpleSC:
	    solveSingleCellNewtonSimple(cell,true);
	    break;
	case NewtonSimpleC:
	    solveSingleCellNewtonSimple(cell,false);
	    break;	    
	default:
	    OPM_THROW(std::runtime_error, "Unknown method " << method_);
	}
    }


    void TransportSolverTwophasePolymer::solveSingleCellBracketing(int cell)
    {
        
	ResidualEquation res_eq(*this, cell);
	ResidualC res(res_eq);
	const double a = 0.0;
	const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system.
	int iters_used;

	// Check if current state is an acceptable solution.
	double res_sc[2];
	double mc, ff;
	res.computeBothResiduals(saturation_[cell], concentration_[cell], res_sc[0], res_sc[1], mc, ff);
	if (norm(res_sc) < tol_) {
	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	    return;
	}

	concentration_[cell] = RootFinder::solve(res, a, b, maxit_, tol_, iters_used);
	cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
	saturation_[cell] = res.lastSaturation();
	fracFlow(saturation_[cell], concentration_[cell], cmax_[cell], cell,
                 fractionalflow_[cell]);
	computeMc(concentration_[cell], mc_[cell]);
    }



    // Newton method, where we first try a Newton step. Then, if it does not work well, we look for
    // the zero of either the residual in s or the residual in c along a specified piecewise linear
    // curve. In these cases, we can use a robust 1d solver.
    void TransportSolverTwophasePolymer::solveSingleCellGradient(int cell)
    {
	int iters_used_falsi = 0;
	const int max_iters_split = maxit_;
	int iters_used_split = 0;

	// Check if current state is an acceptable solution.
	ResidualEquation res_eq(*this, cell);
	double x[2] = {saturation_[cell], saturation_[cell]*concentration_[cell]};
	double res[2];
	double mc;
	double ff;
        double x_c[2];
        scToc(x, x_c);
	res_eq.computeResidual(x_c, res, mc, ff);
	if (norm(res) <= tol_) {
	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
 	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	    return;
	} 

        double x_min[2] = { 0.0, 0.0 };
	double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ };
        double x_min_res_s[2] = { x_min[0], x_min[1] };
        double x_max_res_s[2] = { x_max[0], x_min[0] };
        double x_min_res_sc[2] = { x_min[0], x_min[1] };
        double x_max_res_sc[2] = { x_max[0], x_max[1] };
	double t;
	double t_max;
	double t_out;
	double direction[2];
	double end_point[2];
        double gradient[2];
	ResSOnCurve res_s_on_curve(res_eq);
	ResCOnCurve res_c_on_curve(res_eq);
	bool if_res_s;

 	while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) {
            if (std::abs(res[0]) < std::abs(res[1])) {
                if (res[0] < -tol_) {
                    direction[0] = x_max_res_s[0] - x[0];
                    direction[1] = x_max_res_s[1] - x[1];
                    if_res_s = true;
                } else if (res[0] > tol_) {
                    direction[0] = x_min_res_s[0] - x[0];
                    direction[1] = x_min_res_s[1] - x[1];
                    if_res_s = true;
                } else {
                    scToc(x, x_c);
                    res_eq.computeGradientResS(x_c, res, gradient);
                    // dResS/d(s_) = dResS/ds - c/s*dResS/ds
                    // dResS/d(sc_) = -1/s*dResS/dc
                    if (x[0] > 1e-2*tol_) {
                        // With s,c variables, we should have
                        // direction[0] = -gradient[1];
                        // direction[1] = gradient[0];
                        // With s, sc variables, we get:
                        scToc(x, x_c);
                        direction[0] = 1.0/x[0]*gradient[1];
                        direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1];
                    } else {
                        // acceptable approximation for nonlinear relative permeability.
                        direction[0] = 0.0;
                        direction[1] = gradient[0];
                    }
                    if_res_s = false;
                }
            } else {
                if (res[1] < -tol_) {
                    direction[0] = x_max_res_sc[0] - x[0];
                    direction[1] = x_max_res_sc[1] - x[1];
                    if_res_s = false;
                } else if (res[1] > tol_) {
                    direction[0] = x_min_res_sc[0] - x[0];
                    direction[1] = x_min_res_sc[1] - x[1];
                    if_res_s = false;
                } else {
                    res_eq.computeGradientResC(x, res, gradient);
                    // dResC/d(s_) = dResC/ds - c/s*dResC/ds
                    // dResC/d(sc_) = -1/s*dResC/dc
                    if (x[0] > 1e-2*tol_) {
                        // With s,c variables, we should have
                        // direction[0] = -gradient[1];
                        // direction[1] = gradient[0];
                        // With s, sc variables, we get:
                        scToc(x, x_c);
                        direction[0] = 1.0/x[0]*gradient[1];
                        direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1];
                    } else {
                        // We take 1.0/s*gradient[1]: wrong for linear permeability,
                        // acceptable for nonlinear relative permeability.
                        direction[0] = 1.0 - res_eq.dps;
                        direction[1] = gradient[0];
                    }
                    if_res_s = true;
                }
            }
            if (if_res_s) {
                if (res[0] < 0) {
                    end_point[0] = x_max_res_s[0];
                    end_point[1] = x_max_res_s[1];
                    res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
                    if (res_s_on_curve(t_out) >= 0) {
                        t_max = t_out;
                    }
                } else {
                    end_point[0] = x_min_res_s[0];
                    end_point[1] = x_min_res_s[1];
                    res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
                    if (res_s_on_curve(t_out) <= 0) {
                        t_max = t_out;
                    }
                }
                // Note: In some experiments modifiedRegularFalsi does not yield a result under the given tolerance.
                t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
                res_s_on_curve.curve.computeXOfT(x, t);
            } else {
                if (res[1] < 0) {
                    end_point[0] = x_max_res_sc[0];
                    end_point[1] = x_max_res_sc[1];
                    res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
                    if (res_c_on_curve(t_out) >= 0) {
                        t_max = t_out;
                    }
                } else {
                    end_point[0] = x_min_res_sc[0];
                    end_point[1] = x_min_res_sc[1];
                    res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out);
                    if (res_c_on_curve(t_out) <= 0) {
                        t_max = t_out;
                    }
                }
                t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi);
                res_c_on_curve.curve.computeXOfT(x, t);

            }
            scToc(x, x_c);
            res_eq.computeResidual(x_c, res, mc, ff);
            iters_used_split += 1;
        }
	    


        if ((iters_used_split >=  max_iters_split) && (norm(res) > tol_)) {
            OPM_MESSAGE("Newton for single cell did not work in cell number " << cell);
            solveSingleCellBracketing(cell);
        } else {
            scToc(x, x_c);
            concentration_[cell] = x_c[1];
            cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
            saturation_[cell] = x[0];
            fractionalflow_[cell] = ff;
            mc_[cell] = mc;
        }
    }
    
    void TransportSolverTwophasePolymer::solveSingleCellNewton(int cell)
    {
        const int max_iters_split = maxit_;
	int iters_used_split = 0;

	// Check if current state is an acceptable solution.
	ResidualEquation res_eq(*this, cell);
	double x[2] = {saturation_[cell], concentration_[cell]};
	double res[2];
	double mc;
	double ff;
	res_eq.computeResidual(x, res, mc, ff);
	if (norm(res) <= tol_) {
	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
 	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	    return;
	} else  if (0.99 < x[0]) {
            // x[0] = 0.5; 
            // x[1] = polyprops_.cMax()/2.0;
            // res_eq.computeResidual(x, res, mc, ff);
	}

        const double x_min[2] = { 0.0, 0.0 };
	const double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ };
	bool successfull_newton_step = true;
        // initialize x_new to avoid warning
	double x_new[2] = {0.0, 0.0};
	double res_new[2];
	ResSOnCurve res_s_on_curve(res_eq);
	ResCOnCurve res_c_on_curve(res_eq);

        // We switch to s-sc variable
        x[1] = x[0]*x[1];
        // x_c will contain the s-c variable 
        double x_c[2];
	
 	while ((norm(res) > tol_) &&
	       (iters_used_split < max_iters_split)  &&
	       successfull_newton_step) {
            double dres_s_dsdc[2];
            double dres_c_dsdc[2];
            scToc(x, x_c);
            res_eq.computeJacobiRes(x_c, dres_s_dsdc, dres_c_dsdc);
            double dFx_dx = (dres_s_dsdc[0]-x_c[1]*dres_s_dsdc[1]);
            double dFx_dy;
            if (x[0] < 1e-2*tol_) {
                dFx_dy = 0.0;
            } else {
                dFx_dy = (dres_s_dsdc[1]/x[0]);
            }
            double dFy_dx = (dres_c_dsdc[0]-x_c[1]*dres_c_dsdc[1]);
            double dFy_dy;
            if (x[0] < 1e-2*tol_) {
                dFy_dy = 1.0 - res_eq.dps;
            } else {
                dFy_dy = (dres_c_dsdc[1]/x[0]);
            }
            double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
            double alpha = 1.0;
            int max_lin_it = 100;
            int lin_it = 0;
            res_new[0] = res[0]*2;
            res_new[1] = res[1]*2;
            while((norm(res_new)>norm(res)) && (lin_it<max_lin_it)) {
                x_new[0] = x[0] - alpha*(res[0]*dFy_dy - res[1]*dFx_dy)/det;
                x_new[1] = x[1] - alpha*(res[1]*dFx_dx - res[0]*dFy_dx)/det;
                check_interval(x_min, x_max, x_new);
                scToc(x_new, x_c);
                res_eq.computeResidual(x_c, res_new, mc, ff);
                alpha = alpha/2.0;
                lin_it = lin_it + 1;
            }
            if (lin_it>=max_lin_it) {
                successfull_newton_step = false;
            } else  {
                scToc(x_new, x_c);
                x[0] = x_c[0];
                x[1] = x_c[1];
                res[0] = res_new[0];
                res[1] = res_new[1];
                iters_used_split += 1;
                successfull_newton_step = true;;		    
            }
	}
		
	if ((iters_used_split >=  max_iters_split) && (norm(res) > tol_)) {
	    OPM_MESSAGE("Newton for single cell did not work in cell number " << cell);
	    solveSingleCellBracketing(cell);
	} else {
	    concentration_[cell] = x[1];
	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
	    saturation_[cell] = x[0];
	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	}
    }

    void TransportSolverTwophasePolymer::solveSingleCellNewtonSimple(int cell,bool use_sc)
    {
	const int max_iters_split = maxit_;
	int iters_used_split = 0;

	// Check if current state is an acceptable solution.
	ResidualEquation res_eq(*this, cell);
	double x[2] = {saturation_[cell], concentration_[cell]};
	double res[2];
	double mc;
	double ff;
	res_eq.computeResidual(x, res, mc, ff);
	if (norm(res) <= tol_) {
	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
 	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	    return;
	}else{
	    //*
	    x[0] = saturation_[cell]-res[0];
	    if((x[0]>1) || (x[0]<0)){
		x[0] = 0.5;
	        x[1] = x[1];
	    }
	    if(x[0]>0){
                x[1] =  concentration_[cell]*saturation_[cell]-res[1];
                x[1] = x[1]/x[0];
                if(x[1]> polyprops_.cMax()){
                    x[1]= polyprops_.cMax()/2.0;		
                }
                if(x[1]<0){
                    x[1]=0;
                }
	    }else{
		x[1]=0;
	    }	    	    
	    //x[0]=0.5;x[1]=polyprops_.cMax()/2.0;
	    res_eq.computeResidual(x, res, mc, ff);
	    //*/
	}


	// double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 };
        double x_min[2] = { 0.0, 0.0 };
	double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ };
	bool successfull_newton_step = true;
	double x_new[2];
	double res_new[2];
	ResSOnCurve res_s_on_curve(res_eq);
	ResCOnCurve res_c_on_curve(res_eq);
	
 	while ((norm(res) > tol_) &&
	       (iters_used_split < max_iters_split)  &&
	       successfull_newton_step) {
	    // We first try a Newton step
	    double dres_s_dsdc[2];
	    double dres_c_dsdc[2];
	    double dx=tol_;
	    double tmp_x[2];
	    if(!(x[0]>0)){
		tmp_x[0]=dx;
		tmp_x[1]=0;
	    }else{
		tmp_x[0]=x[0];
		tmp_x[1]=x[1];
	    }
	    res_eq.computeJacobiRes(tmp_x, dres_s_dsdc, dres_c_dsdc);
	    double dFx_dx(0),dFx_dy(0),dFy_dx(0),dFy_dy(0);
	    double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
	    if(use_sc){
		dFx_dx=(dres_s_dsdc[0]-tmp_x[1]*dres_s_dsdc[1]/tmp_x[0]);
		dFx_dy= (dres_s_dsdc[1]/tmp_x[0]);
		dFy_dx=(dres_c_dsdc[0]-tmp_x[1]*dres_c_dsdc[1]/tmp_x[0]);
		dFy_dy= (dres_c_dsdc[1]/tmp_x[0]);
	    }else{
		dFx_dx= dres_s_dsdc[0];
		dFx_dy= dres_s_dsdc[1];
		dFy_dx= dres_c_dsdc[0];
		dFy_dy= dres_c_dsdc[1];
	    }
	    det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
	    if(det==0){
		OPM_THROW(std::runtime_error, "det is zero");
	    }

	    double alpha=1;
	    int max_lin_it=10;
	    int lin_it=0;
	    /*
	      std::cout << "Nonlinear " << iters_used_split
	      << "  " << norm(res_new)
	      << "  " << norm(res) << std::endl;*/
	    /*
	      std::cout  << "x" << std::endl;
	      std::cout  << " " <<  x[0] << " " <<  x[1] << std::endl;
	      std::cout  << "dF" << std::endl;
	      std::cout  << " " <<  dFx_dx << " " <<  dFx_dy << std::endl;
	      std::cout  << " " <<  dFy_dx << " " <<  dFy_dy << std::endl;
	      std::cout  << "dres" << std::endl;
	      std::cout  << " " <<  dres_s_dsdc[0] << " " <<  dres_s_dsdc[1] << std::endl;
	      std::cout  << " " <<  dres_c_dsdc[0] << " " <<  dres_c_dsdc[1] << std::endl;
	      std::cout  << std::endl;
	    */
	    res_new[0]=res[0]*2;
	    res_new[1]=res[1]*2;
	    double update[2]={(res[0]*dFy_dy - res[1]*dFx_dy)/det,
			      (res[1]*dFx_dx - res[0]*dFy_dx)/det};
	    while((norm(res_new)>norm(res)) && (lin_it<max_lin_it)) {
		x_new[0] = x[0] - alpha*update[0];
		if(use_sc){
		    x_new[1] = x[1]*x[0] - alpha*update[1];
		    x_new[1] = (x_new[0]>0) ?  x_new[1]/x_new[0] : 0.0;
		}else{
		    x_new[1] = x[1] - alpha*update[1];
		}
		check_interval(x_min, x_max, x_new);
		res_eq.computeResidual(x_new, res_new, mc, ff);
		//  std::cout << " " << res_new[0] << "  " << res_new[1] << std::endl;
		alpha=alpha/2.0;
		lin_it=lin_it+1;
		//  std::cout << "Linear iterations" << lin_it << "  " << norm(res_new) << std::endl;
	    }
	    if (lin_it>=max_lin_it) {
		successfull_newton_step = false;
	    } else  {
		x[0] = x_new[0];
		x[1] = x_new[1];
		res[0] = res_new[0];
		res[1] = res_new[1];
		iters_used_split += 1;
		successfull_newton_step = true;;		    
	    }
	    //	    std::cout << "Nonlinear " << iters_used_split << "  " << norm(res) << std::endl;
	}
		
	if ((iters_used_split >=  max_iters_split) || (norm(res) > tol_)) {
	    OPM_MESSAGE("NewtonSimple for single cell did not work in cell number " << cell);
	    solveSingleCellBracketing(cell);
	} else {
	    concentration_[cell] = x[1];
	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
	    saturation_[cell] = x[0];
	    fractionalflow_[cell] = ff;
	    mc_[cell] = mc;
	}
    }



    void TransportSolverTwophasePolymer::solveMultiCell(const int num_cells, const int* cells)
    {
	double max_s_change = 0.0;
	double max_c_change = 0.0;
	int num_iters = 0;
	// Must store state variables before we start.
	std::vector<double> s0(num_cells);
	std::vector<double> c0(num_cells);
	std::vector<double> cmax0(num_cells);
	// Must set initial fractional flows etc. before we start.
	for (int i = 0; i < num_cells; ++i) {
	    const int cell = cells[i];
	    fracFlow(saturation_[cell], concentration_[cell], cmax_[cell],
                     cell, fractionalflow_[cell]);
            computeMc(concentration_[cell], mc_[cell]);
	    s0[i] = saturation_[cell];
	    c0[i] = concentration_[cell];
	    cmax0[i] = cmax_[i];
	}
	do {
	    // int max_s_change_cell = -1;
	    // int max_c_change_cell = -1;
	    max_s_change = 0.0;
	    max_c_change = 0.0;
	    for (int i = 0; i < num_cells; ++i) {
		const int cell = cells[i];
		const double old_s = saturation_[cell];
		const double old_c = concentration_[cell];
		saturation_[cell] = s0[i];
		concentration_[cell] = c0[i];
		cmax_[cell] = cmax0[i];
		solveSingleCell(cell);
		// std::cout << "cell = " << cell << "    delta s = " << saturation_[cell] - old_s << std::endl;
		// if (max_s_change < std::fabs(saturation_[cell] - old_s)) {
		//     max_s_change_cell = cell;
		// }
		// if (max_c_change < std::fabs(concentration_[cell] - old_c)) {
		//     max_c_change_cell = cell;
		// }
		max_s_change = std::max(max_s_change, std::fabs(saturation_[cell] - old_s));
		max_c_change = std::max(max_c_change, std::fabs(concentration_[cell] - old_c));
	    }
	    // std::cout << "Iter = " << num_iters << "    max_s_change = " << max_s_change
	    // 	      << "    in cell " << max_change_cell << std::endl;
	} while (((max_s_change > tol_) || (max_c_change > tol_)) && ++num_iters < maxit_);
	if (max_s_change > tol_) {
	    OPM_THROW(std::runtime_error, "In solveMultiCell(), we did not converge after "
		  << num_iters << " iterations. Delta s = " << max_s_change);
	}
	if (max_c_change > tol_) {
	    OPM_THROW(std::runtime_error, "In solveMultiCell(), we did not converge after "
		  << num_iters << " iterations. Delta c = " << max_c_change);
	}
	std::cout << "Solved " << num_cells << " cell multicell problem in "
		  << num_iters << " iterations." << std::endl;
    }

    void TransportSolverTwophasePolymer::fracFlow(double s, double c, double cmax,
                                         int cell, double& ff) const
    {
        double dummy[2];
        fracFlowBoth(s, c, cmax, cell, ff,  dummy, false);
    }

    void TransportSolverTwophasePolymer::fracFlowWithDer(double s, double c, double cmax,
                                                int cell, double& ff,
                                                double* dff_dsdc) const
    {
        fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true);
    }

    void TransportSolverTwophasePolymer::fracFlowBoth(double s, double c, double cmax, int cell,
                                             double& ff, double* dff_dsdc,
                                             bool if_with_der) const
    {
	double relperm[2];
	double drelperm_ds[4];
        double sat[2] = {s, 1 - s};
        if (if_with_der) {
            props_.relperm(1, sat, &cell, relperm, drelperm_ds);
        } else {
            props_.relperm(1, sat, &cell, relperm, 0);
        }
        double mob[2];
        double dmob_ds[4];
        double dmob_dc[2];
	double dmobwat_dc;
        polyprops_.effectiveMobilitiesBoth(c, cmax, visc_, relperm, drelperm_ds,
                                           mob, dmob_ds, dmobwat_dc, if_with_der);
	
 	ff = mob[0]/(mob[0] + mob[1]);
        if (if_with_der) {
            dmob_dc[0] = dmobwat_dc;
            dmob_dc[1] = 0.;
            //dff_dsdc[0] = (dmob_ds[0]*mob[1] + dmob_ds[3]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s
	    // at the moment the dmob_ds only have diagonal elements since the saturation is derivated out in effectiveMobilitiesBouth
	    dff_dsdc[0] = ((dmob_ds[0]-dmob_ds[2])*mob[1] - (dmob_ds[1]-dmob_ds[3])*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s 
            dff_dsdc[1] = (dmob_dc[0]*mob[1] - dmob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to c
        }
    }

    void TransportSolverTwophasePolymer::computeMc(double c, double& mc) const
    {
        polyprops_.computeMc(c, mc);
    }

    void TransportSolverTwophasePolymer::computeMcWithDer(double c, double& mc,
                                                 double &dmc_dc) const
    {
        polyprops_.computeMcWithDer(c, mc, dmc_dc);
    }



    TransportSolverTwophasePolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq,
                                                        const double c_init)
        : res_c_eq_(res_c_eq),
          c(c_init)
    {
    }

    double TransportSolverTwophasePolymer::ResidualSGrav::operator()(double s) const
    {
        return res_c_eq_.computeGravResidualS(s, c);
    }


    // Residual for concentration equation for gravity segregation
    //
    // res_c = s*(1 - dps)*c - s0*( - dps)*c0
    //     +  dtpv*sum_{j adj i}( mc * gravmod_ij * gf_ij ).
    //  \TODO doc me
    // where ...
    // Influxes are negative, outfluxes positive.

    TransportSolverTwophasePolymer::ResidualCGrav::ResidualCGrav(const TransportSolverTwophasePolymer& tmodel,
                                                        const std::vector<int>& cells,
                                                        const int pos,
                                                        const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
        : tm(tmodel),
          cell(cells[pos]),
          s0(tm.saturation_[cell]),
          c0(tm.concentration_[cell]),
          cmax0(tm.cmax0_[cell]),
          porosity(tm.porosity_[cell]),
          dtpv(tm.dt_/tm.porevolume_[cell]),
          dps(tm.polyprops_.deadPoreVol()),
          rhor(tm.polyprops_.rockDensity())

    {
        last_s = s0;
        nbcell[0] = -1;
        gf[0] = 0.0;
        if (pos > 0) {
            nbcell[0] = cells[pos - 1];
            gf[0] = -gravflux[pos - 1];
        }
        nbcell[1] = -1;
        gf[1] = 0.0;
        if (pos < int(cells.size() - 1)) {
            nbcell[1] = cells[pos + 1];
            gf[1] = gravflux[pos];
        }

        tm.polyprops_.adsorption(c0, cmax0, c_ads0);
    }

    double TransportSolverTwophasePolymer::ResidualCGrav::operator()(double c) const
    {

        ResidualSGrav res_s(*this);
        res_s.c = c;
        int iters_used;
        last_s =  RootFinder::solve(res_s, last_s, 0.0, 1.0,
                                    tm.maxit_, tm.tol_,
                                    iters_used);
        return computeGravResidualC(last_s, c);

    }

    double TransportSolverTwophasePolymer::ResidualCGrav::computeGravResidualS(double s, double c) const
    {

        double mobcell[2];
        tm.mobility(s, c, cell, mobcell);

        double res = s - s0;

        for (int nb = 0; nb < 2; ++nb) {
            if (nbcell[nb] != -1) {
                double m[2];
                if (gf[nb] < 0.0) {
                    m[0] = mobcell[0];
                    m[1] = tm.mob_[2*nbcell[nb] + 1];
                } else {
                    m[0] = tm.mob_[2*nbcell[nb]];
                    m[1] = mobcell[1];
                }
                if (m[0] + m[1] > 0.0) {
                    res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]);
                }
            }
        }
        return res;
    }

    double TransportSolverTwophasePolymer::ResidualCGrav::computeGravResidualC(double s, double c) const
    {

        double mobcell[2];
        tm.mobility(s, c, cell,  mobcell);
        double c_ads;
        tm.polyprops_.adsorption(c, cmax0, c_ads);

        double res = (1 - dps)*s*c - (1 - dps)*s0*c0
            + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0);
        for (int nb = 0; nb < 2; ++nb) {
            if (nbcell[nb] != -1) {
                double m[2];
                double mc;
                if (gf[nb] < 0.0) {
                    m[0] = mobcell[0];
                    tm.computeMc(c, mc);
                    m[1] = tm.mob_[2*nbcell[nb] + 1];
                } else {
                    m[0] = tm.mob_[2*nbcell[nb]];
                    mc = tm.mc_[nbcell[nb]];
                    m[1] = mobcell[1];
                }
                if (m[0] + m[1] > 0.0) {
                    res += -dtpv*gf[nb]*mc*m[0]*m[1]/(m[0] + m[1]);
                }
            }
        }
        return res;
    }

    double TransportSolverTwophasePolymer::ResidualCGrav::lastSaturation() const
    {
        return last_s;
    }


    void TransportSolverTwophasePolymer::mobility(double s, double c, int cell, double* mob) const
    {
	double sat[2] = { s, 1.0 - s };
        double relperm[2];
	props_.relperm(1, sat, &cell, relperm, 0);
        polyprops_.effectiveMobilities(c, cmax0_[cell], visc_, relperm, mob);
    }


    void TransportSolverTwophasePolymer::initGravity(const double* grav)
    {
        // Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
        std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
        const int nf = grid_.number_of_faces;
        const int dim = grid_.dimensions;
        gravflux_.resize(nf);
        tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
        tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &gravflux_[0]);
        const double delta_rho = props_.density()[0] - props_.density()[1];
        for (int f = 0; f < nf; ++f) {
            const int* c = &grid_.face_cells[2*f];
            double gdz = 0.0;
            if (c[0] != -1 && c[1] != -1) {
                for (int d = 0; d < dim; ++d) {
                    gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]);
                }
            }
            gravflux_[f] *= delta_rho*gdz;
        }
    }


    void TransportSolverTwophasePolymer::solveSingleCellGravity(const std::vector<int>& cells,
                                                       const int pos,
                                                       const double* gravflux)
    {
        const int cell = cells[pos];
        ResidualCGrav res_c(*this, cells, pos, gravflux);

        // Check if current state is an acceptable solution.
	double res_sc[2];
	res_sc[0]=res_c.computeGravResidualS(saturation_[cell], concentration_[cell]);
	res_sc[1]=res_c.computeGravResidualC(saturation_[cell], concentration_[cell]);

        if (norm(res_sc) < tol_) {
            return;
        }

	const double a = 0.0;
	const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system.
        int iters_used;
        concentration_[cell] = RootFinder::solve(res_c, a, b, maxit_, tol_, iters_used);
        saturation_[cell] = res_c.lastSaturation();
        cmax_[cell] = std::max(cmax0_[cell], concentration_[cell]);
        saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]);
        computeMc(concentration_[cell], mc_[cell]);
        mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]);
    }

    int TransportSolverTwophasePolymer::solveGravityColumn(const std::vector<int>& cells)
    {
        // Set up column gravflux.
        const int nc = cells.size();
        std::vector<double> col_gravflux(nc - 1);
        for (int ci = 0; ci < nc - 1; ++ci) {
	    const int cell = cells[ci];
	    const int next_cell = cells[ci + 1];
	    for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
		const int face = grid_.cell_faces[j];
		const int c1 = grid_.face_cells[2*face + 0];
                const int c2 = grid_.face_cells[2*face + 1];
		if (c1 == next_cell || c2 == next_cell) {
                    const double gf = gravflux_[face];
                    col_gravflux[ci] = (c1 == cell) ? gf : -gf;
		}
	    }
        }

        // Store initial saturation s0
        s0_.resize(nc);
        c0_.resize(nc);
        for (int ci = 0; ci < nc; ++ci) {
            s0_[ci] = saturation_[cells[ci]];
            c0_[ci] = concentration_[cells[ci]];
        }

        // Solve single cell problems, repeating if necessary.
	double max_sc_change = 0.0;
        int num_iters = 0;
        do {
            max_sc_change = 0.0;
            for (int ci = 0; ci < nc; ++ci) {
                const int ci2 = nc - ci - 1;
                double old_s[2] = { saturation_[cells[ci]],
                                    saturation_[cells[ci2]] };
                double old_c[2] = { concentration_[cells[ci]],
                                    concentration_[cells[ci2]] };
                saturation_[cells[ci]] = s0_[ci];
                concentration_[cells[ci]] = c0_[ci];
                solveSingleCellGravity(cells, ci, &col_gravflux[0]);
                saturation_[cells[ci2]] = s0_[ci2];
                concentration_[cells[ci2]] = c0_[ci2];
                solveSingleCellGravity(cells, ci2, &col_gravflux[0]);
                max_sc_change = std::max(max_sc_change, 0.25*(std::fabs(saturation_[cells[ci]] - old_s[0]) + 
                                                              std::fabs(concentration_[cells[ci]] - old_c[0]) +
                                                              std::fabs(saturation_[cells[ci2]] - old_s[1]) +
                                                              std::fabs(concentration_[cells[ci2]] - old_c[1])));
            }
            // std::cout << "Iter = " << num_iters << "    max_s_change = " << max_s_change << std::endl;
	} while (max_sc_change > tol_ && ++num_iters < maxit_);

	if (max_sc_change > tol_) {
	    OPM_THROW(std::runtime_error, "In solveGravityColumn(), we did not converge after "
	    	  << num_iters << " iterations. Delta s = " << max_sc_change);
	}
        return num_iters + 1;
    }


    void TransportSolverTwophasePolymer::solveGravity(const std::vector<std::vector<int> >& columns,
                                             const double* porevolume,
                                             const double dt,
                                             std::vector<double>& saturation,
                                             std::vector<double>& concentration,
                                             std::vector<double>& cmax)
    {
        // initialize variables.
        porevolume_ = porevolume;
        dt_ = dt;
        toWaterSat(saturation, saturation_);
        concentration_ = &concentration[0];
        cmax_ = &cmax[0];
        const int nc = grid_.number_of_cells;
        cmax0_.resize(nc);
        std::copy(cmax.begin(), cmax.end(), &cmax0_[0]);

        // Initialize mobilities.
        mob_.resize(2*nc);

        for (int cell = 0; cell < nc; ++cell) {
            mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]);
        }


        // 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);
    }

    void TransportSolverTwophasePolymer::scToc(const double* x, double* x_c) const {
        x_c[0] = x[0];
        if (x[0] < 1e-2*tol_) {
            x_c[1] = 0.5*polyprops_.cMax();
        } else {
            x_c[1] = x[1]/x[0];
        }
    }

} // namespace Opm


namespace
{
    bool check_interval(const double* xmin, const double* xmax, double* x) {
        bool test = false;
        if (x[0] < xmin[0]) {
            test = true;
            x[0] = xmin[0];
        } else if (x[0] > xmax[0]) {
            test = true;
            x[0] = xmax[0];
        }
        if (x[1] < xmin[1]) {
            test = true;
            x[1] = xmin[1];
        } else if (x[1] > xmax[1]) {
            test = true;
            x[1] = xmax[1];
        }
        return test;
    }


    CurveInSCPlane::CurveInSCPlane()
    {
    }

    // Setup the curve (see comment above).
    // The curve is parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding
    // rectangle. x_out=(s_out, c_out) denotes the values of s and c at that point.
    void CurveInSCPlane::setup(const double* x, const double* direction,
                               const double* end_point, const double* x_min,
                               const double* x_max, const double tol,
                               double& t_max_out, double& t_out_out)
    {
	x_[0] = x[0];
	x_[1] = x[1];
	x_max_[0] = x_max[0];
	x_max_[1] = x_max[1];
	x_min_[0] = x_min[0];
	x_min_[1] = x_min[1];
	direction_[0] = direction[0];
	direction_[1] = direction[1];
	end_point_[0] = end_point[0];
	end_point_[1] = end_point[1];
        const double size_direction = std::abs(direction_[0]) + std::abs(direction_[1]);
	if (size_direction < tol) {
	    direction_[0] = end_point_[0]-x_[0];
	    direction_[1] = end_point_[1]-x_[1];
	} else if ((end_point_[0]-x_[0])*direction_[0] + (end_point_[1]-x_[1])*direction_[1] < 0) {
            direction_[0] *= -1.0;
            direction_[1] *= -1.0;
        }
	bool t0_exists = true;
	double t0 = 0; // dummy default value (so that compiler does not complain).
	if (direction_[0] > 0) {
	    t0 = (x_max_[0] - x_[0])/direction_[0];
	} else if (direction_[0] < 0) {
	    t0 = (x_min_[0] - x_[0])/direction_[0];
	} else {
	    t0_exists = false;
	}
	bool t1_exists = true;
	double t1 = 0; // dummy default value.
	if (direction_[1] > 0) {
	    t1 = (x_max_[1] - x_[1])/direction_[1];
	} else if (direction_[1] < 0) {
	    t1 = (x_min_[1] - x_[1])/direction_[1];
	} else {
	    t1_exists = false;
	}
	if (t0_exists) {
	    if (t1_exists) {
		t_out_ = std::min(t0, t1);
	    } else {
		t_out_ = t0;
	    }
	} else if (t1_exists) {
	    t_out_ = t1;
	} else {
	    OPM_THROW(std::runtime_error, "Direction illegal: is a zero vector.");
	}
	x_out_[0] = x_[0] + t_out_*direction_[0];
	x_out_[1] = x_[1] + t_out_*direction_[1];
	t_max_ = t_out_ + 1;
	t_max_out = t_max_;
	t_out_out = t_out_;
    }


    // Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve)
    void CurveInSCPlane::computeXOfT(double* x_of_t, const double t) const {
	if (t <= t_out_) {
	    x_of_t[0] = x_[0] + t*direction_[0];
            x_of_t[1] = x_[1] + t*direction_[1];
	} else {
	    x_of_t[0] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[0] + end_point_[0]*(t - t_out_));
	    x_of_t[1] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[1] + end_point_[1]*(t - t_out_));
	}
    }


    ResSOnCurve::ResSOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq)
	: res_eq_(res_eq)
    {
    }

    double ResSOnCurve::operator()(const double t) const
    {
	double x_of_t[2];
        double x_c[2];
	curve.computeXOfT(x_of_t, t);
        res_eq_.tm.scToc(x_of_t, x_c);
	return res_eq_.computeResidualS(x_c);
    }

    ResCOnCurve::ResCOnCurve(const Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq)
	: res_eq_(res_eq)
    {
    }

    double ResCOnCurve::operator()(const double t) const
    {
	double x_of_t[2];
        double x_c[2];
	curve.computeXOfT(x_of_t, t);
        res_eq_.tm.scToc(x_of_t, x_c);
	return res_eq_.computeResidualC(x_c);
    }

} // Anonymous namespace


/* Local Variables:    */
/* c-basic-offset:4    */
/* End:                */