TransportSolverTwophasePolymer.hpp

Go to the documentation of this file.

/*
  Copyright 2012 SINTEF ICT, Applied Mathematics.
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED
 
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <vector>
#include <list>
 
struct UnstructuredGrid;
 
namespace Opm
{
 
    class IncompPropertiesInterface;
 
    class TransportSolverTwophasePolymer : public ReorderSolverInterface
    {
    public:
 
    enum SingleCellMethod { Bracketing, Newton, Gradient, NewtonSimpleSC, NewtonSimpleC};
        enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded)
 
    TransportSolverTwophasePolymer(const UnstructuredGrid& grid,
                                       const IncompPropertiesInterface& props,
                                       const PolymerProperties& polyprops,
                                       const SingleCellMethod method,
                                       const double tol,
                                       const int maxit);
 
        void setPreferredMethod(SingleCellMethod method);
 
    void 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);
 
        void 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);
 
    public: // But should be made private...
    virtual void solveSingleCell(const int cell);
    virtual void solveMultiCell(const int num_cells, const int* cells);
    void solveSingleCellBracketing(int cell);
    void solveSingleCellNewton(int cell);
    void solveSingleCellGradient(int cell);
    void solveSingleCellNewtonSimple(int cell,bool use_sc);
    class ResidualEquation;
 
        void initGravity(const double* grav);
        void solveSingleCellGravity(const std::vector<int>& cells,
                                    const int pos,
                                    const double* gravflux);
        int solveGravityColumn(const std::vector<int>& cells);
        void scToc(const double* x, double* x_c) const;
 
        #ifdef PROFILING
        class Newton_Iter {
        public:
            bool res_s;
            int cell;
            double s;
            double c;
 
            Newton_Iter(bool res_s_val, int cell_val, double s_val, double c_val) {
                res_s = res_s_val;
                cell = cell_val;
                s = s_val;
                c = c_val;
            }
        };
 
        std::list<Newton_Iter> res_counts;
        #endif
 
 
    private:
    const UnstructuredGrid& grid_;
    const double* porosity_;
    const double* porevolume_;  // one volume per cell
    const IncompPropertiesInterface& props_;
    const PolymerProperties& polyprops_;
    std::vector<double> smin_;
    std::vector<double> smax_;
    double tol_;
    double maxit_;
 
    const double* darcyflux_;   // one flux per grid face
    const double* source_;      // one source per cell
    const double* polymer_inflow_c_;
    double dt_;
        std::vector<double> saturation_; // one per cell, only water saturation!
    double* concentration_;
    double* cmax_;
    std::vector<double> fractionalflow_;  // one per cell
    std::vector<double> mc_;  // one per cell
    const double* visc_;
    SingleCellMethod method_;
    double adhoc_safety_;
    
        // For gravity segregation.
        std::vector<double> gravflux_;
        std::vector<double> mob_;
        std::vector<double> cmax0_;
        // For gravity segregation, column variables
        std::vector<double> s0_;
        std::vector<double> c0_;
 
    struct ResidualC;
    struct ResidualS;
 
    class ResidualCGrav;
    class ResidualSGrav;
 
 
    void fracFlow(double s, double c, double cmax, int cell, double& ff) const;
    void fracFlowWithDer(double s, double c, double cmax, int cell, double& ff,
                               double* dff_dsdc) const;
    void fracFlowBoth(double s, double c, double cmax, int cell, double& ff,
                          double* dff_dsdc, bool if_with_der) const;
    void computeMc(double c, double& mc) const;
    void computeMcWithDer(double c, double& mc, double& dmc_dc) const;
        void mobility(double s, double c, int cell, double* mob) const;
    };
 
} // namespace Opm
 
#endif // OPM_TRANSPORTSOLVERTWOPHASEPOLYMER_HEADER_INCLUDED