Add simple solvent model

assumes:
- solvent is immiscible in the oil phase
- gas pvt and relperms are used for the solvent
- no initial solvent in the model

Solvent is injected using the WSOLVENT keyword
TODO: Make it possible to change WSOLVENT

CMakeLists_files.cmake

@@ -77,9 +77,10 @@ endif()
 list (APPEND EXAMPLE_SOURCE_FILES
 	examples/find_zero.cpp
 	examples/flow.cpp
+	examples/flow_extended.cpp
 	examples/sim_2p_incomp_ad.cpp
 	examples/sim_simple.cpp
-        examples/opm_init_check.cpp
+	examples/opm_init_check.cpp
 	)
 
 # programs listed here will not only be compiled, but also marked for
@@ -105,6 +106,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp
 	opm/autodiff/CPRPreconditioner.hpp
+	opm/autodiff/ExtendedBlackoilModel.hpp
+	opm/autodiff/ExtendedBlackoilModel_impl.hpp
 	opm/autodiff/fastSparseProduct.hpp
 	opm/autodiff/DuneMatrix.hpp
 	opm/autodiff/ExtractParallelGridInformationToISTL.hpp
@@ -124,6 +127,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SimulatorBase.hpp
 	opm/autodiff/SimulatorBase_impl.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
+	opm/autodiff/SimulatorFullyImplicitExtendedBlackoil.hpp
+	opm/autodiff/SimulatorFullyImplicitExtendedBlackoil_impl.hpp
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp

examples/flow_extended.cpp

@@ -0,0 +1,426 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 IRIS AS
+
+  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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+
+#include <dune/common/version.hh>
+
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
+#include <dune/common/parallel/mpihelper.hh>
+#else
+#include <dune/common/mpihelper.hh>
+#endif
+
+#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID
+#define USE_DUNE_CORNERPOINTGRID 1
+#include <dune/grid/CpGrid.hpp>
+#include <dune/grid/common/GridAdapter.hpp>
+#else
+#undef USE_DUNE_CORNERPOINTGRID
+#endif
+
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cornerpoint_grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/autodiff/GridHelpers.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/initStateEquil.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/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
+
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+#include <opm/autodiff/SimulatorFullyImplicitExtendedBlackoil.hpp>
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/RedistributeDataHandles.hpp>
+
+#include <opm/core/utility/share_obj.hpp>
+
+#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
+#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
+#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <memory>
+#include <algorithm>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <cstdlib>
+
+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;
+#if USE_DUNE_CORNERPOINTGRID
+    // Must ensure an instance of the helper is created to initialise MPI.
+    const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
+    const int mpi_rank = mpi_helper.rank();
+    const int mpi_size = mpi_helper.size();
+#else
+    // default values for serial run
+    const int mpi_rank = 0;
+    const int mpi_size = 1;
+#endif
+
+    // Write parameters used for later reference. (only if rank is zero)
+    const bool output_cout = ( mpi_rank == 0 );
+
+    if(output_cout)
+    {
+        std::cout << "**********************************************************************\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "*                   This is Flow-Solvent (version XXXX.XX)           *\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "*     Flow-Solvent is a simulator for fully implicit three-phase,    *\n";
+        std::cout << "*    forth component (black-oil + solvent) flow, and is part of OPM. *\n";
+        std::cout << "*           For more information see http://opm-project.org          *\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "**********************************************************************\n\n";
+    }
+
+    // Read parameters, see if a deck was specified on the command line.
+    if ( output_cout )
+    {
+        std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+    }
+
+    parameter::ParameterGroup param(argc, argv, false, output_cout);
+    if( !output_cout )
+    {
+        param.disableOutput();
+    }
+
+    if (!param.unhandledArguments().empty()) {
+        if (param.unhandledArguments().size() != 1) {
+            std::cerr << "You can only specify a single input deck on the command line.\n";
+            return EXIT_FAILURE;
+        } else {
+            param.insertParameter("deck_filename", param.unhandledArguments()[0]);
+        }
+    }
+
+    // We must have an input deck. Grid and props will be read from that.
+    if (!param.has("deck_filename")) {
+        std::cerr << "This program must be run with an input deck.\n"
+            "Specify the deck filename either\n"
+            "    a) as a command line argument by itself\n"
+            "    b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
+            "    c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
+        return EXIT_FAILURE;
+    }
+
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    std::string deck_filename = param.get<std::string>("deck_filename");
+
+    // Write parameters used for later reference. (only if rank is zero)
+    bool output = ( mpi_rank == 0 ) && param.getDefault("output", true);
+    std::string output_dir;
+    if (output) {
+        // Create output directory if needed.
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            std::cerr << "Creating directories failed: " << fpath << std::endl;
+            return EXIT_FAILURE;
+        }
+        // Write simulation parameters.
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    std::string logFile = output_dir + "/LOGFILE.txt";
+    Opm::ParserPtr parser(new Opm::Parser());
+    {
+        std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
+        std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
+
+        Opm::OpmLog::addBackend( "STREAM" , streamLog );
+        Opm::OpmLog::addBackend( "COUNTER" , counterLog );
+    }
+
+    Opm::DeckConstPtr deck;
+    std::shared_ptr<EclipseState> eclipseState;
+    try {
+        deck = parser->parseFile(deck_filename);
+        Opm::checkDeck(deck);
+        eclipseState.reset(new Opm::EclipseState(deck));
+    }
+    catch (const std::invalid_argument& e) {
+        std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
+        std::cerr << "Exception caught: " << e.what() << std::endl;
+        return EXIT_FAILURE;
+    }
+
+
+    std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
+#if USE_DUNE_CORNERPOINTGRID
+    // Dune::CpGrid as grid manager
+    typedef Dune::CpGrid  Grid;
+    // Grid init
+    Grid grid;
+    grid.processEclipseFormat(deck, false, false, false, porv);
+#else
+    // UnstructuredGrid as grid manager
+    typedef UnstructuredGrid  Grid;
+    GridManager gridManager( eclipseState->getEclipseGrid(), porv );
+    const Grid& grid = *(gridManager.c_grid());
+#endif
+
+    // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
+    if (param.has("output_interval")) {
+        int output_interval = param.get<int>("output_interval");
+        IOConfigPtr ioConfig = eclipseState->getIOConfig();
+        ioConfig->overrideRestartWriteInterval((size_t)output_interval);
+    }
+
+    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
+    Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu );
+
+    // Rock and fluid init
+    BlackoilPropertiesFromDeck props( deck, eclipseState,
+                                      Opm::UgGridHelpers::numCells(grid),
+                                      Opm::UgGridHelpers::globalCell(grid),
+                                      Opm::UgGridHelpers::cartDims(grid),
+                                      Opm::UgGridHelpers::beginCellCentroids(grid),
+                                      Opm::UgGridHelpers::dimensions(grid), param);
+
+    BlackoilPropsAdFromDeck new_props( deck, eclipseState, grid );
+    // check_well_controls = param.getDefault("check_well_controls", false);
+    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+    // Rock compressibility.
+    RockCompressibility rock_comp(deck, eclipseState);
+
+    // Gravity.
+    gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+
+    ExtendedBlackoilState state;
+    // Init state variables (saturation and pressure).
+    if (param.has("init_saturation")) {
+        initStateBasic(Opm::UgGridHelpers::numCells(grid),
+                       Opm::UgGridHelpers::globalCell(grid),
+                       Opm::UgGridHelpers::cartDims(grid),
+                       Opm::UgGridHelpers::numFaces(grid),
+                       Opm::UgGridHelpers::faceCells(grid),
+                       Opm::UgGridHelpers::beginFaceCentroids(grid),
+                       Opm::UgGridHelpers::beginCellCentroids(grid),
+                       Opm::UgGridHelpers::dimensions(grid),
+                       props, param, gravity[2], state);
+
+        initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state);
+
+        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
+        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
+            const int numPhases = props.numPhases();
+            const int numCells  = Opm::UgGridHelpers::numCells(grid);
+            for (int c = 0; c < numCells; ++c) {
+                state.gasoilratio()[c] = state.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
+                    / state.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
+            }
+        }
+    } else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) {
+        state.init(Opm::UgGridHelpers::numCells(grid),
+                   Opm::UgGridHelpers::numFaces(grid),
+                   props.numPhases());
+        const double grav = param.getDefault("gravity", unit::gravity);
+        initStateEquil(grid, props, deck, eclipseState, grav, state);
+        state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
+    } else {
+        initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
+                                  Opm::UgGridHelpers::globalCell(grid),
+                                  Opm::UgGridHelpers::numFaces(grid),
+                                  Opm::UgGridHelpers::faceCells(grid),
+                                  Opm::UgGridHelpers::beginFaceCentroids(grid),
+                                  Opm::UgGridHelpers::beginCellCentroids(grid),
+                                  Opm::UgGridHelpers::dimensions(grid),
+                                  props, deck, gravity[2], state);
+    }
+
+
+    // The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
+    if (deck->hasKeyword("SWATINIT")) {
+        const int numCells = Opm::UgGridHelpers::numCells(grid);
+        std::vector<int> cells(numCells);
+        for (int c = 0; c < numCells; ++c) { cells[c] = c; }
+        std::vector<double> pc = state.saturation();
+        props.capPress(numCells, state.saturation().data(), cells.data(), pc.data(),NULL);
+        new_props.setSwatInitScaling(state.saturation(),pc);
+    }
+
+    //state.solvent_saturation()[49] = 0.1;
+    //state.saturation()[49*3 + 2] -= 0.1;
+
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+#if USE_DUNE_CORNERPOINTGRID
+    if(output_cout)
+    {
+        std::cout << std::endl << "Warning: use of local perm is not yet implemented for CpGrid!" << std::endl << std::endl;
+    }
+    const bool use_local_perm = false;
+#else
+    const bool use_local_perm = param.getDefault("use_local_perm", true);
+#endif
+
+    DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav);
+    boost::any parallel_information;
+
+    // At this point all properties and state variables are correctly initialized
+    // If there are more than one processors involved, we now repartition the grid
+    // and initilialize new properties and states for it.
+    if( mpi_size > 1 )
+    {
+        if( param.getDefault("output_matlab", false) || param.getDefault("output_ecl", true) )
+        {
+        	std::cerr << "We only support vtk output during parallel runs. \n"
+                      << "Please use \"output_matlab=false output_ecl=false\" to deactivate the \n"
+                      << "other outputs!" << std::endl;
+        	return EXIT_FAILURE;
+        }
+
+        Opm::distributeGridAndData( grid, eclipseState, state, new_props, geoprops, parallel_information, use_local_perm );
+    }
+
+    // Solver for Newton iterations.
+    std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
+    if (param.getDefault("use_interleaved", false)) {
+        fis_solver.reset(new NewtonIterationBlackoilInterleaved(param));
+    } else if (param.getDefault("use_cpr", true)) {
+        fis_solver.reset(new NewtonIterationBlackoilCPR(param));
+    } else {
+        fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
+    }
+
+    Opm::ScheduleConstPtr schedule = eclipseState->getSchedule();
+    Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
+    SimulatorTimer simtimer;
+
+    // initialize variables
+    simtimer.init(timeMap);
+
+    std::vector<double> threshold_pressures = thresholdPressures(eclipseState, grid);
+
+    SimulatorFullyImplicitExtendedBlackoil< Grid >  simulator(param,
+                                                      grid,
+                                                      geoprops,
+                                                      new_props,
+                                                      rock_comp.isActive() ? &rock_comp : 0,
+                                                      *fis_solver,
+                                                      grav,
+                                                      deck->hasKeyword("DISGAS"),
+                                                      deck->hasKeyword("VAPOIL"),
+                                                      eclipseState,
+                                                      outputWriter,
+                                                      deck,
+                                                      threshold_pressures,
+                                                      deck->hasKeyword("SOLVENT") );
+
+    if (!schedule->initOnly()){
+        if( output_cout )
+        {
+            std::cout << "\n\n================ Starting main simulation loop ===============\n"
+                      << std::flush;
+        }
+
+        SimulatorReport fullReport = simulator.run(simtimer, state);
+
+        if( output_cout )
+        {
+            std::cout << "\n\n================    End of simulation     ===============\n\n";
+            fullReport.reportFullyImplicit(std::cout);
+        }
+
+        if (output) {
+            std::string filename = output_dir + "/walltime.txt";
+            std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+            fullReport.reportParam(tot_os);
+            warnIfUnusedParams(param);
+        }
+    } else {
+        outputWriter.writeInit( simtimer );
+        if ( output_cout )
+        {
+            std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
+        }
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    return EXIT_FAILURE;
+}
+

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -843,7 +843,7 @@ namespace detail {
         V trans_all(transi.size() + trans_nnc.size());
         trans_all << transi, trans_nnc;
 
-        const std::vector<ADB> kr = computeRelPerm(state);
+        const std::vector<ADB> kr = asImpl().computeRelPerm(state);
         for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
             asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
 
@@ -1569,10 +1569,10 @@ namespace detail {
             }
         }
 
-        const V sumSat = sw + so + sg;
-        sw = sw / sumSat;
-        so = so / sumSat;
-        sg = sg / sumSat;
+        //const V sumSat = sw + so + sg;
+        //sw = sw / sumSat;
+        //so = so / sumSat;
+        //sg = sg / sumSat;
 
         // Update the reservoir_state
         for (int c = 0; c < nc; ++c) {
@@ -1607,7 +1607,6 @@ namespace detail {
             rv = rv_old - drv_limited;
         }
 
-
         // Sg is used as primal variable for water only cells.
         const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
         auto watOnly = sw >  (1 - epsilon);

opm/autodiff/ExtendedBlackoilModel.hpp

@@ -0,0 +1,287 @@
+/*
+  Copyright 2015 IRIS AS
+
+  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_EXTENDEDBLACKOILMODEL_HEADER_INCLUDED
+#define OPM_EXTENDEDBLACKOILMODEL_HEADER_INCLUDED
+
+#include <opm/autodiff/BlackoilModelBase.hpp>
+#include <opm/autodiff/BlackoilModelParameters.hpp>
+#include <opm/autodiff/ExtendedBlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
+
+namespace Opm {
+
+    /// A model implementation for three-phase black oil
+    /// with one extra component.
+    ///
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    template<class Grid>
+    class ExtendedBlackoilModel : public BlackoilModelBase<Grid, ExtendedBlackoilModel<Grid> >
+    {
+    public:
+
+        // ---------  Types and enums  ---------
+
+        typedef BlackoilModelBase<Grid, ExtendedBlackoilModel<Grid> > Base;
+        typedef typename Base::ReservoirState ReservoirState;
+        typedef typename Base::WellState WellState;
+        // The next line requires C++11 support available in g++ 4.7.
+        // friend Base;
+        friend class BlackoilModelBase<Grid, ExtendedBlackoilModel<Grid> >;
+
+        /// Construct the model. It will retain references to the
+        /// arguments of this functions, and they are expected to
+        /// remain in scope for the lifetime of the solver.
+        /// \param[in] param               parameters
+        /// \param[in] grid                grid data structure
+        /// \param[in] fluid               fluid properties
+        /// \param[in] geo                 rock properties
+        /// \param[in] rock_comp_props     if non-null, rock compressibility properties
+        /// \param[in] wells               well structure
+        /// \param[in] linsolver           linear solver
+        /// \param[in] has_disgas          turn on dissolved gas
+        /// \param[in] has_vapoil          turn on vaporized oil feature
+        /// \param[in] has_polymer         turn on polymer feature
+        /// \param[in] has_plyshlog        true when PLYSHLOG keyword available
+        /// \param[in] has_shrate          true when PLYSHLOG keyword available
+        /// \param[in] wells_rep_radius    representative radius of well perforations during shear effects calculation
+        /// \param[in] wells_perf_length   perforation length for well perforations
+        /// \param[in] wells_bore_diameter wellbore diameters for well performations
+        /// \param[in] terminal_output  request output to cout/cerr
+        ExtendedBlackoilModel(const typename Base::ModelParameters&   param,
+                             const Grid&                             grid,
+                             const BlackoilPropsAdInterface&         fluid,
+                             const DerivedGeology&                   geo,
+                             const RockCompressibility*              rock_comp_props,
+                             const Wells*                            wells,
+                             const NewtonIterationBlackoilInterface& linsolver,
+                             const EclipseStateConstPtr              eclState,
+                             const bool                              has_disgas,
+                             const bool                              has_vapoil,
+                             const bool                              terminal_output,
+                             const bool                              has_solvent);
+
+        /// Called once before each time step.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void prepareStep(const double dt,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Called once after each time step.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void afterStep(const double dt,
+                       ReservoirState& reservoir_state,
+                       WellState& well_state);
+
+        /// Apply an update to the primary variables, chopped if appropriate.
+        /// \param[in]      dx                updates to apply to primary variables
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void updateState(const V& dx,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Compute convergence based on total mass balance (tol_mb) and maximum
+        /// residual mass balance (tol_cnv).
+        /// \param[in]   dt          timestep length
+        /// \param[in]   iteration   current iteration number
+        bool getConvergence(const double dt, const int iteration);
+
+        /// Assemble the residual and Jacobian of the nonlinear system.
+        /// \param[in]      reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
+        void assemble(const ReservoirState& reservoir_state,
+                      WellState& well_state,
+                      const bool initial_assembly);
+
+
+    protected:
+
+        // ---------  Types and enums  ---------
+
+        typedef typename Base::SolutionState SolutionState;
+        typedef typename Base::DataBlock DataBlock;
+        enum { Solvent = CanonicalVariablePositions::Next };
+
+        // ---------  Data members  ---------
+        const bool has_solvent_;
+        const int solvent_pos_;
+
+        // Need to declare Base members we want to use here.
+        using Base::grid_;
+        using Base::fluid_;
+        using Base::geo_;
+        using Base::rock_comp_props_;
+        using Base::wells_;
+        using Base::linsolver_;
+        using Base::active_;
+        using Base::canph_;
+        using Base::cells_;
+        using Base::ops_;
+        using Base::wops_;
+        using Base::has_disgas_;
+        using Base::has_vapoil_;
+        using Base::param_;
+        using Base::use_threshold_pressure_;
+        using Base::threshold_pressures_by_interior_face_;
+        using Base::rq_;
+        using Base::phaseCondition_;
+        using Base::well_perforation_pressure_diffs_;
+        using Base::residual_;
+        using Base::terminal_output_;
+        using Base::primalVariable_;
+        using Base::pvdt_;
+
+        // ---------  Protected methods  ---------
+
+        // Need to declare Base members we want to use here.
+        using Base::wellsActive;
+        using Base::wells;
+        using Base::variableState;
+        using Base::computePressures;
+        using Base::computeGasPressure;
+        using Base::applyThresholdPressures;
+        using Base::fluidViscosity;
+        using Base::fluidReciprocFVF;
+        using Base::fluidDensity;
+        using Base::fluidRsSat;
+        using Base::fluidRvSat;
+        using Base::poroMult;
+        using Base::transMult;
+        using Base::updatePrimalVariableFromState;
+        using Base::updatePhaseCondFromPrimalVariable;
+        using Base::dpMaxRel;
+        using Base::dsMax;
+        using Base::drMaxRel;
+        using Base::maxResidualAllowed;
+
+        using Base::updateWellControls;
+        using Base::computeWellConnectionPressures;
+        using Base::addWellControlEq;
+
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+
+        void
+        makeConstantState(SolutionState& state) const;
+
+        std::vector<V>
+        variableStateInitials(const ReservoirState& x,
+                              const WellState& xw) const;
+
+        std::vector<int>
+        variableStateIndices() const;
+
+        SolutionState
+        variableStateExtractVars(const ReservoirState& x,
+                                 const std::vector<int>& indices,
+                                 std::vector<ADB>& vars) const;
+
+        void
+        computeAccum(const SolutionState& state,
+                     const int            aix  );
+
+        void
+        assembleMassBalanceEq(const SolutionState& state);
+
+        void
+        addWellContributionToMassBalanceEq(std::vector<ADB>& cq_s,
+                                           const SolutionState& state,
+                                           WellState& xw);
+
+        void
+        computeMassFlux(const int               actph ,
+                        const V&                transi,
+                        const ADB&              kr    ,
+                        const ADB&              p     ,
+                        const SolutionState&    state );
+
+        const std::vector<PhasePresence>
+        phaseCondition() const {return this->phaseCondition_;}
+
+        /// \brief Compute the reduction within the convergence check.
+        /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. tempV.col(i) contains the
+        ///                   values
+        ///                  for phase i.
+        /// \param[in] R     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+        ///                   of R for the phase i.
+        /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+        ///                   maximum of tempV for the phase i.
+        /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+        ///                   of B for the phase i.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \return The total pore volume over all cells.
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
+                             std::array<double,MaxNumPhases+1>& R_sum,
+                             std::array<double,MaxNumPhases+1>& maxCoeff,
+                             std::array<double,MaxNumPhases+1>& B_avg,
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int nw) const;
+
+    };
+
+
+
+    /// Need to include concentration in our state variables, otherwise all is as
+    /// the default blackoil model.
+    struct ExtendedBlackoilSolutionState : public DefaultBlackoilSolutionState
+    {
+        explicit ExtendedBlackoilSolutionState(const int np)
+            : DefaultBlackoilSolutionState(np),
+              solvent_saturation( ADB::null())
+        {
+        }
+        ADB solvent_saturation;
+    };
+
+
+
+    /// Providing types by template specialisation of ModelTraits for BlackoilPolymerModel.
+    template <class Grid>
+    struct ModelTraits< ExtendedBlackoilModel<Grid> >
+    {
+        typedef ExtendedBlackoilState ReservoirState;
+        typedef WellStateFullyImplicitBlackoilSolvent WellState;
+        typedef BlackoilModelParameters ModelParameters;
+        typedef ExtendedBlackoilSolutionState SolutionState;
+    };
+
+} // namespace Opm
+
+#include "ExtendedBlackoilModel_impl.hpp"
+
+#endif // OPM_EXTENDEDBLACKOILMODEL_HEADER_INCLUDED

opm/autodiff/ExtendedBlackoilModel_impl.hpp

@@ -0,0 +1,701 @@
+/*
+  Copyright 2015 IRIS AS
+
+  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_EXTENDEDBLACKOILMODEL_IMPL_HEADER_INCLUDED
+#define OPM_EXTENDEDBLACKOILMODEL_IMPL_HEADER_INCLUDED
+
+#include <opm/autodiff/ExtendedBlackoilModel.hpp>
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/WellDensitySegmented.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/Exceptions.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/well_controls.h>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+#include <limits>
+
+namespace Opm {
+
+
+
+    namespace detail {
+
+        template <class PU>
+        int solventPos(const PU& pu)
+        {
+            const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+            int pos = 0;
+            for (int phase = 0; phase < maxnp; ++phase) {
+                if (pu.phase_used[phase]) {
+                    pos++;
+                }
+            }
+
+            return pos;
+        }
+
+    } // namespace detail
+
+
+
+    template <class Grid>
+    ExtendedBlackoilModel<Grid>::ExtendedBlackoilModel(const typename Base::ModelParameters&   param,
+                                                     const Grid&                             grid,
+                                                     const BlackoilPropsAdInterface&         fluid,
+                                                     const DerivedGeology&                   geo,
+                                                     const RockCompressibility*              rock_comp_props,
+                                                     const Wells*                            wells,
+                                                     const NewtonIterationBlackoilInterface& linsolver,
+                                                     const EclipseStateConstPtr              eclState,
+                                                     const bool                              has_disgas,
+                                                     const bool                              has_vapoil,
+                                                     const bool                              terminal_output,
+                                                     const bool                              has_solvent)
+        : Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
+               eclState, has_disgas, has_vapoil, terminal_output),
+          has_solvent_(has_solvent),
+          solvent_pos_(detail::solventPos(fluid.phaseUsage()))
+    {
+        if (has_solvent_) {
+
+            // If deck has solvent, residual_ should contain solvent equation.
+            rq_.resize(fluid_.numPhases() + 1);
+            residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
+            assert(solvent_pos_ == fluid_.numPhases());
+            if (has_vapoil_) {
+                OPM_THROW(std::runtime_error, "Solvent option only works with dead gas\n");
+            }
+        }
+    }
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::
+    prepareStep(const double dt,
+                ReservoirState& reservoir_state,
+                WellState& well_state)
+    {
+        Base::prepareStep(dt, reservoir_state, well_state);
+    }
+
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::
+    afterStep(const double /* dt */,
+              ReservoirState& /* reservoir_state */,
+              WellState& /* well_state */)
+    {
+
+    }
+
+
+
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::makeConstantState(SolutionState& state) const
+    {
+        Base::makeConstantState(state);
+        state.solvent_saturation = ADB::constant(state.solvent_saturation.value());
+    }
+
+
+
+
+
+    template <class Grid>
+    std::vector<V>
+    ExtendedBlackoilModel<Grid>::variableStateInitials(const ReservoirState& x,
+                                                      const WellState& xw) const
+    {
+        std::vector<V> vars0 = Base::variableStateInitials(x, xw);
+        assert(int(vars0.size()) == fluid_.numPhases() + 2);
+
+        // Initial polymer concentration.
+        if (has_solvent_) {
+            assert (not x.solvent_saturation().empty());
+            const int nc = x.solvent_saturation().size();
+            const V ss = Eigen::Map<const V>(&x.solvent_saturation()[0], nc);
+            // Solvent belongs after other reservoir vars but before well vars.
+            auto solvent_pos = vars0.begin() + fluid_.numPhases();
+            assert(solvent_pos == vars0.end() - 2);
+            vars0.insert(solvent_pos, ss);
+        }
+        return vars0;
+    }
+
+
+
+
+
+    template <class Grid>
+    std::vector<int>
+    ExtendedBlackoilModel<Grid>::variableStateIndices() const
+    {
+        std::vector<int> ind = Base::variableStateIndices();
+        assert(ind.size() == 5);
+        if (has_solvent_) {
+            ind.resize(6);
+            // Solvent belongs after other reservoir vars but before well vars.
+            ind[Solvent] = fluid_.numPhases();
+            // Solvent is pushing back the well vars.
+            ++ind[Qs];
+            ++ind[Bhp];
+        }
+        return ind;
+    }
+
+
+
+
+    template <class Grid>
+    typename ExtendedBlackoilModel<Grid>::SolutionState
+    ExtendedBlackoilModel<Grid>::variableStateExtractVars(const ReservoirState& x,
+                                                         const std::vector<int>& indices,
+                                                         std::vector<ADB>& vars) const
+    {
+        SolutionState state = Base::variableStateExtractVars(x, indices, vars);
+        if (has_solvent_) {
+            state.solvent_saturation = std::move(vars[indices[Solvent]]);
+            if (active_[ Oil ]) {
+                // Note that so is never a primary variable.
+                const Opm::PhaseUsage pu = fluid_.phaseUsage();
+                state.saturation[pu.phase_pos[ Oil ]] -= state.solvent_saturation;
+            }
+        }
+        return state;
+    }
+
+
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::computeAccum(const SolutionState& state,
+                                             const int            aix  )
+    {
+        Base::computeAccum(state, aix);
+
+        // Compute accumulation of the solvent
+        if (has_solvent_) {
+            const ADB& press = state.pressure;
+            const ADB& ss = state.solvent_saturation;
+            const ADB pv_mult = poroMult(press); // also computed in Base::computeAccum, could be optimized.
+            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+            // Compute solvent accumulation term.
+            // WARNING: Currently the solvent uses the same pvt as gas
+            // TODO: Add support for PVDS, SDENSITY
+            rq_[solvent_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Gas]].b * ss;
+
+        }
+    }
+
+
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::
+    assembleMassBalanceEq(const SolutionState& state)
+    {
+
+        Base::assembleMassBalanceEq(state);
+
+        if (has_solvent_) {
+            residual_.material_balance_eq[ solvent_pos_ ] =
+                pvdt_ * (rq_[solvent_pos_].accum[1] - rq_[solvent_pos_].accum[0])
+                + ops_.div*rq_[solvent_pos_].mflux;
+        }
+
+    }
+
+
+
+
+
+    template <class Grid>
+    void ExtendedBlackoilModel<Grid>::addWellContributionToMassBalanceEq(std::vector<ADB>& cq_s,
+                                                                        const SolutionState& state,
+                                                                        WellState& xw)
+
+    {
+
+        // Add well contributions to solvent mass balance equation
+
+        Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
+
+        if (has_solvent_) {
+            const int nperf = wells().well_connpos[wells().number_of_wells];
+            const int nc = Opm::AutoDiffGrid::numCells(grid_);
+
+            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+            const ADB zero = ADB::constant(V::Zero(nc));
+            const ADB& ss = state.solvent_saturation;
+            const ADB& sg = (active_[ Gas ]
+                             ? state.saturation[ pu.phase_pos[ Gas ] ]
+                             : zero);
+
+            Selector<double> zero_selector(ss.value(), Selector<double>::Zero);
+            V F_solvent = zero_selector.select(ss, ss / (ss + sg)).value();
+
+            const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
+            const int nw = wells().number_of_wells;
+            V wellSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
+
+            for (int w = 0; w < nw; ++w) {
+                if(wells().type[w] == PRODUCER) {
+                    for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
+                        wellSolventFraction[perf] = F_solvent[well_cells[perf]];
+                    }
+                }
+            }
+
+            const ADB& rs_perfcells = subset(state.rs, well_cells);
+            int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
+            int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
+            // remove contribution from the dissolved gas.
+            const ADB cq_s_solvent = wellSolventFraction * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]);
+            //cq_s[gas_pos] = ( ones - wellSolventFraction ) *  cq_s[gas_pos];;
+            residual_.material_balance_eq[solvent_pos_] -= superset(cq_s_solvent, well_cells, nc);
+            residual_.material_balance_eq[gas_pos] += superset(cq_s_solvent, well_cells, nc);
+
+        }
+    }
+
+
+
+
+
+
+    template <class Grid>
+    void ExtendedBlackoilModel<Grid>::updateState(const V& dx,
+                                                  ReservoirState& reservoir_state,
+                                                  WellState& well_state)
+    {
+
+        if (has_solvent_) {
+            // Extract solvent change.
+            const int np = fluid_.numPhases();
+            const int nc = Opm::AutoDiffGrid::numCells(grid_);
+            const V zero = V::Zero(nc);
+            const int solvent_start = nc * np;
+            const V dss = subset(dx, Span(nc, 1, solvent_start));
+
+            // Create new dx with the dss part deleted.
+            V modified_dx = V::Zero(dx.size() - nc);
+            modified_dx.head(solvent_start) = dx.head(solvent_start);
+            const int tail_len = dx.size() - solvent_start - nc;
+            modified_dx.tail(tail_len) = dx.tail(tail_len);
+
+            // Call base version.
+            Base::updateState(modified_dx, reservoir_state, well_state);
+
+            // Update solvent.
+            const V ss_old = Eigen::Map<const V>(&reservoir_state.solvent_saturation()[0], nc, 1);
+            const V ss = (ss_old - dss).max(zero);
+            std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin());
+
+            // adjust oil saturation
+            const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+            const int pos = pu.phase_pos[ Oil ];
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pos] = 1 - ss[c] - reservoir_state.saturation()[c*np + 0] - reservoir_state.saturation()[c*np + 2] ;
+            }
+
+        } else {
+            // Just forward call to base version.
+            Base::updateState(dx, reservoir_state, well_state);
+        }
+    }
+
+
+
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::computeMassFlux(const int               actph ,
+                                                const V&                transi,
+                                                const ADB&              kr    ,
+                                                const ADB&              phasePressure,
+                                                const SolutionState&    state)
+    {
+        Base::computeMassFlux(actph, transi, kr, phasePressure, state);
+
+        const int canonicalPhaseIdx = canph_[ actph ];
+        if (canonicalPhaseIdx == Gas) {
+            if (has_solvent_) {
+                const int  nc   = Opm::UgGridHelpers::numCells(grid_);
+
+                const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+                const ADB zero = ADB::constant(V::Zero(nc));
+                const ADB& ss = state.solvent_saturation;
+                const ADB& sg = (active_[ Gas ]
+                                 ? state.saturation[ pu.phase_pos[ Gas ] ]
+                                 : zero);
+
+                Selector<double> zero_selector(ss.value(), Selector<double>::Zero);
+                ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
+                V ones = V::Constant(nc, 1.0);
+
+                // WARNING: The solvent mobility is simply given as a fraction of the gas mobility
+                // TODO: Add support for gas/solvent relative permeability functions (SSFN)
+                rq_[solvent_pos_].mob = F_solvent * rq_[actph].mob; //tr_mult * F_solvent * kr / mu;
+                rq_[actph].mob = (ones - F_solvent) * rq_[actph].mob;
+
+                // WARNING: use gas values for the solvent
+                rq_[solvent_pos_].b = rq_[actph].b;
+                rq_[solvent_pos_].dh = rq_[actph].dh;
+                UpwindSelector<double> upwind(grid_, ops_, rq_[solvent_pos_].dh.value());
+                // Compute solvent flux.
+                rq_[solvent_pos_].mflux = upwind.select(rq_[solvent_pos_].b * rq_[solvent_pos_].mob) * (transi * rq_[solvent_pos_].dh);
+
+            }
+        }
+
+    }
+
+    template <class Grid>
+    std::vector<ADB>
+    ExtendedBlackoilModel<Grid>::computeRelPerm(const SolutionState& state) const
+    {
+        using namespace Opm::AutoDiffGrid;
+        const int               nc   = numCells(grid_);
+
+        const ADB zero = ADB::constant(V::Zero(nc));
+
+        const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+        const ADB& sw = (active_[ Water ]
+                         ? state.saturation[ pu.phase_pos[ Water ] ]
+                         : zero);
+
+        const ADB& so = (active_[ Oil ]
+                         ? state.saturation[ pu.phase_pos[ Oil ] ]
+                         : zero);
+
+        const ADB& sg = (active_[ Gas ]
+                         ? state.saturation[ pu.phase_pos[ Gas ] ]
+                         : zero);
+
+        if (has_solvent_) {
+            const ADB& ss = state.solvent_saturation;
+            return fluid_.relperm(sw, so, sg+ss, cells_);
+        } else {
+            return fluid_.relperm(sw, so, sg, cells_);
+        }
+
+    }
+
+
+    template <class Grid>
+    void
+    ExtendedBlackoilModel<Grid>::assemble(const ReservoirState& reservoir_state,
+                                         WellState& well_state,
+                                         const bool initial_assembly)
+    {
+
+        using namespace Opm::AutoDiffGrid;
+
+        // Possibly switch well controls and updating well state to
+        // get reasonable initial conditions for the wells
+        updateWellControls(well_state);
+
+        // Create the primary variables.
+        SolutionState state = variableState(reservoir_state, well_state);
+
+        if (initial_assembly) {
+            // Create the (constant, derivativeless) initial state.
+            SolutionState state0 = state;
+            makeConstantState(state0);
+            // Compute initial accumulation contributions
+            // and well connection pressures.
+            computeAccum(state0, 0);
+            computeWellConnectionPressures(state0, well_state);
+        }
+
+        // -------- Mass balance equations --------
+        assembleMassBalanceEq(state);
+
+        // -------- Well equations ----------
+        if ( ! wellsActive() ) {
+            return;
+        }
+
+        V aliveWells;
+
+        const int np = wells().number_of_phases;
+        std::vector<ADB> cq_s(np, ADB::null());
+
+        const int nw = wells().number_of_wells;
+        const int nperf = wells().well_connpos[nw];
+        const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
+
+        std::vector<ADB> mob_perfcells(np, ADB::null());
+        std::vector<ADB> b_perfcells(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+            mob_perfcells[phase] = subset(rq_[phase].mob, well_cells);
+            b_perfcells[phase] = subset(rq_[phase].b, well_cells);
+        }
+
+        if (has_solvent_) {
+            int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
+            mob_perfcells[gas_pos] += subset(rq_[solvent_pos_].mob, well_cells);
+        }
+        if (param_.solve_welleq_initially_ && initial_assembly) {
+            // solve the well equations as a pre-processing step
+            Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
+        }
+
+        Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
+        Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
+        Base::addWellFluxEq(cq_s, state);
+        addWellContributionToMassBalanceEq(cq_s, state, well_state);
+        Base::addWellControlEq(state, well_state, aliveWells);
+
+    }
+
+
+    template <class Grid>
+    double
+    ExtendedBlackoilModel<Grid>::convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
+                                                     const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
+                                                     const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
+                                                     std::array<double,MaxNumPhases+1>& R_sum,
+                                                     std::array<double,MaxNumPhases+1>& maxCoeff,
+                                                     std::array<double,MaxNumPhases+1>& B_avg,
+                                                     std::vector<double>& maxNormWell,
+                                                     int nc,
+                                                     int nw) const
+    {
+        // Do the global reductions
+#if HAVE_MPI
+        if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+        {
+            const ParallelISTLInformation& info =
+                boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+
+            // Compute the global number of cells and porevolume
+            std::vector<int> v(nc, 1);
+            auto nc_and_pv = std::tuple<int, double>(0, 0.0);
+            auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
+                                                        Opm::Reduction::makeGlobalSumFunctor<double>());
+            auto nc_and_pv_containers  = std::make_tuple(v, geo_.poreVolume());
+            info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
+
+            for ( int idx=0; idx<MaxNumPhases+1; ++idx )
+            {
+                if ((idx == MaxNumPhases && has_solvent_) || active_[idx]) { // Dealing with solvent *or* an active phase.
+                    auto values     = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
+                    auto containers = std::make_tuple(B.col(idx),
+                                                      tempV.col(idx),
+                                                      R.col(idx));
+                    auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalMaxFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalSumFunctor<double>());
+                    info.computeReduction(containers, operators, values);
+                    B_avg[idx]       = std::get<0>(values)/std::get<0>(nc_and_pv);
+                    maxCoeff[idx]    = std::get<1>(values);
+                    R_sum[idx]       = std::get<2>(values);
+                    if (idx != MaxNumPhases) { // We do not compute a well flux residual for solvent.
+                        maxNormWell[idx] = 0.0;
+                        for ( int w=0; w<nw; ++w ) {
+                            maxNormWell[idx]  = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
+                        }
+                    }
+                }
+                else
+                {
+                    maxNormWell[idx] = R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
+                }
+            }
+            info.communicator().max(&maxNormWell[0], MaxNumPhases+1);
+            // Compute pore volume
+            return std::get<1>(nc_and_pv);
+        }
+        else
+#endif
+        {
+            for ( int idx=0; idx<MaxNumPhases+1; ++idx )
+            {
+                if (((idx == MaxNumPhases && has_solvent_) || active_[idx]) ) { // Dealing with solvent *or* an active phase.
+                    B_avg[idx] = B.col(idx).sum()/nc;
+                    maxCoeff[idx] = tempV.col(idx).maxCoeff();
+                    R_sum[idx] = R.col(idx).sum();
+                }
+                else
+                {
+                    R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
+                }
+                if (idx != MaxNumPhases) { // We do not compute a well flux residual for polymer.
+                    maxNormWell[idx] = 0.0;
+                    for ( int w=0; w<nw; ++w ) {
+                        maxNormWell[idx]  = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
+                    }
+                }
+            }
+            // Compute total pore volume
+            return geo_.poreVolume().sum();
+        }
+    }
+
+
+
+
+
+
+
+    template <class Grid>
+    bool
+    ExtendedBlackoilModel<Grid>::getConvergence(const double dt, const int iteration)
+    {
+        const double tol_mb    = param_.tolerance_mb_;
+        const double tol_cnv   = param_.tolerance_cnv_;
+        const double tol_wells = param_.tolerance_wells_;
+
+        const int nc = Opm::AutoDiffGrid::numCells(grid_);
+        const int nw = wellsActive() ? wells().number_of_wells : 0;
+        const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+        const V pv = geo_.poreVolume();
+
+        const std::vector<PhasePresence> cond = phaseCondition();
+
+        std::array<double,MaxNumPhases+1> CNV                   = {{0., 0., 0., 0.}};
+        std::array<double,MaxNumPhases+1> R_sum                 = {{0., 0., 0., 0.}};
+        std::array<double,MaxNumPhases+1> B_avg                 = {{0., 0., 0., 0.}};
+        std::array<double,MaxNumPhases+1> maxCoeff              = {{0., 0., 0., 0.}};
+        std::array<double,MaxNumPhases+1> mass_balance_residual = {{0., 0., 0., 0.}};
+        std::array<double,MaxNumPhases> well_flux_residual    = {{0., 0., 0.}};
+        std::size_t cols = MaxNumPhases+1; // needed to pass the correct type to Eigen
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> B(nc, cols);
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> R(nc, cols);
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases+1> tempV(nc, cols);
+        std::vector<double> maxNormWell(MaxNumPhases);
+
+        for ( int idx=0; idx<MaxNumPhases; ++idx )
+        {
+            if (active_[idx]) {
+                const int pos    = pu.phase_pos[idx];
+                const ADB& tempB = rq_[pos].b;
+                B.col(idx)       = 1./tempB.value();
+                R.col(idx)       = residual_.material_balance_eq[idx].value();
+                tempV.col(idx)   = R.col(idx).abs()/pv;
+            }
+        }
+        if (has_solvent_) {
+            const ADB& tempB = rq_[solvent_pos_].b;
+            B.col(MaxNumPhases) = 1. / tempB.value();
+            R.col(MaxNumPhases) = residual_.material_balance_eq[solvent_pos_].value();
+            tempV.col(MaxNumPhases) = R.col(MaxNumPhases).abs()/pv;
+        }
+
+        const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg,
+                                                  maxNormWell, nc, nw);
+
+        bool converged_MB = true;
+        bool converged_CNV = true;
+        bool converged_Well = true;
+        // Finish computation
+        for ( int idx = 0; idx < (MaxNumPhases + 1) ; ++idx )
+        {
+            CNV[idx]                   = B_avg[idx] * dt * maxCoeff[idx];
+            mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
+            converged_MB               = converged_MB && (mass_balance_residual[idx] < tol_mb);
+            converged_CNV              = converged_CNV && (CNV[idx] < tol_cnv);
+            if (idx != MaxNumPhases) { // No well flux residual for polymer.
+                well_flux_residual[idx]    = B_avg[idx] * maxNormWell[idx];
+                converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
+            }
+        }
+
+        const double residualWell     = detail::infinityNormWell(residual_.well_eq,
+                                                                 linsolver_.parallelInformation());
+        converged_Well   = converged_Well && (residualWell < Opm::unit::barsa);
+        const bool   converged        = converged_MB && converged_CNV && converged_Well;
+
+        // if one of the residuals is NaN, throw exception, so that the solver can be restarted
+        if (std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
+            std::isnan(mass_balance_residual[Oil])   || mass_balance_residual[Oil]   > maxResidualAllowed() ||
+            std::isnan(mass_balance_residual[Gas])   || mass_balance_residual[Gas]   > maxResidualAllowed() ||
+            std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() ||
+            std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() ||
+            std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() ||
+            std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() ||
+            std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() ||
+            std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ||
+            std::isnan(residualWell)     || residualWell     > maxResidualAllowed() )
+        {
+            OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!");
+        }
+
+        if ( terminal_output_ )
+        {
+            // Only rank 0 does print to std::cout
+            if (iteration == 0) {
+                std::cout << "\nIter  MB(WATER)   MB(OIL)    MB(GAS)    MB(SOLVENT)    CNVW       CNVO       CNVG       CNVS   W-FLUX(W)  W-FLUX(O)  W-FLUX(G)\n";
+            }
+            const std::streamsize oprec = std::cout.precision(3);
+            const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
+            std::cout << std::setw(4) << iteration
+                      << std::setw(11) << mass_balance_residual[Water]
+                      << std::setw(11) << mass_balance_residual[Oil]
+                      << std::setw(11) << mass_balance_residual[Gas]
+                      << std::setw(11) << mass_balance_residual[solvent_pos_]
+                      << std::setw(11) << CNV[Water]
+                      << std::setw(11) << CNV[Oil]
+                      << std::setw(11) << CNV[Gas]
+                      << std::setw(11) << CNV[solvent_pos_]
+                      << std::setw(11) << well_flux_residual[Water]
+                      << std::setw(11) << well_flux_residual[Oil]
+                      << std::setw(11) << well_flux_residual[Gas]
+                      << std::endl;
+            std::cout.precision(oprec);
+            std::cout.flags(oflags);
+        }
+        return converged;
+    }
+}
+
+
+#endif // OPM_EXTENDEDBLACKOIL_IMPL_HEADER_INCLUDED

opm/autodiff/ExtendedBlackoilState.hpp

@@ -0,0 +1,60 @@
+/*
+  Copyright 2015 IRIS AS, 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_EXTENDEDBLACKOILSTATE_HEADER_INCLUDED
+#define OPM_EXTENDEDBLACKOILSTATE_HEADER_INCLUDED
+
+#include <opm/autodiff/ExtendedBlackoilState.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/grid.h>
+#include <vector>
+
+
+namespace Opm
+{
+
+    /// Simulator state for blackoil simulator with solvent.
+    /// We use the Blackoil state parameters.
+    class ExtendedBlackoilState : public BlackoilState
+    {
+    public:
+        void init(const UnstructuredGrid& g, int num_phases)
+        {
+            this->init(g.number_of_cells, g.number_of_faces, num_phases);
+        }
+
+        void init(int number_of_cells, int number_of_faces, int num_phases)
+        {
+            BlackoilState::init(number_of_cells, number_of_faces, num_phases);
+            solventId_  = SimulatorState::registerCellData( "SSOL", 1 );
+        }
+
+        std::vector<double>& solvent_saturation()    { return cellData()[ solventId_ ]; }
+        const std::vector<double>& solvent_saturation() const    { return cellData()[ solventId_ ]; }
+
+    private:
+        int solventId_;
+    };
+
+} // namespace Opm
+
+
+
+
+#endif // OPM_EXTENDEDBLACKOILSTATE_HEADER_INCLUDED