changed: ewoms/models/blackoil -> opm/models/blackoil

examples/problems/reservoirproblem.hh

@@ -28,7 +28,7 @@
 #ifndef EWOMS_RESERVOIR_PROBLEM_HH
 #define EWOMS_RESERVOIR_PROBLEM_HH
 
-#include <ewoms/models/blackoil/blackoilproperties.hh>
+#include <opm/models/blackoil/blackoilproperties.hh>
 
 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>

examples/reservoir_blackoil_ecfv.cpp

@@ -29,7 +29,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilmodel.hh>
 #include <ewoms/disc/ecfv/ecfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 

examples/reservoir_blackoil_vcfv.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilmodel.hh>
 #include <ewoms/disc/vcfv/vcfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 

opm/models/blackoil/blackoilboundaryratevector.hh

@@ -0,0 +1,271 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilBoundaryRateVector
+ */
+#ifndef EWOMS_BLACK_OIL_BOUNDARY_RATE_VECTOR_HH
+#define EWOMS_BLACK_OIL_BOUNDARY_RATE_VECTOR_HH
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/constraintsolvers/NcpFlash.hpp>
+
+#include "blackoilintensivequantities.hh"
+#include "blackoilenergymodules.hh"
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief Implements a boundary vector for the fully implicit black-oil model.
+ */
+template <class TypeTag>
+class BlackOilBoundaryRateVector : public GET_PROP_TYPE(TypeTag, RateVector)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { contiEnergyEqIdx = Indices::contiEnergyEqIdx };
+    enum { enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam) };
+
+    static constexpr bool blackoilConserveSurfaceVolume = GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume);
+
+    typedef Opm::BlackOilEnergyModule<TypeTag, enableEnergy> EnergyModule;
+
+public:
+    /*!
+     * \brief Default constructor
+     */
+    BlackOilBoundaryRateVector() : ParentType()
+    {}
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar)
+     */
+    BlackOilBoundaryRateVector(Scalar value) : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const ImmiscibleBoundaryRateVector& )
+     */
+    BlackOilBoundaryRateVector(const BlackOilBoundaryRateVector& value) = default;
+    BlackOilBoundaryRateVector& operator=(const BlackOilBoundaryRateVector& value) = default;
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::setFreeFlow
+     */
+    template <class Context, class FluidState>
+    void setFreeFlow(const Context& context,
+                     unsigned bfIdx,
+                     unsigned timeIdx,
+                     const FluidState& fluidState)
+    {
+        ExtensiveQuantities extQuants;
+        extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState);
+        const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
+        unsigned focusDofIdx = context.focusDofIndex();
+        unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx);
+
+        ////////
+        // advective fluxes of all components in all phases
+        ////////
+        (*this) = 0.0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            const auto& pBoundary = fluidState.pressure(phaseIdx);
+            const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
+
+            RateVector tmp;
+
+            // mass conservation
+            if (pBoundary < pInside)
+                // outflux
+                LocalResidual::template evalPhaseFluxes_<Evaluation>(tmp,
+                                                                     phaseIdx,
+                                                                     insideIntQuants.pvtRegionIndex(),
+                                                                     extQuants,
+                                                                     insideIntQuants.fluidState());
+            else if (pBoundary > pInside) {
+                typedef typename std::conditional<std::is_same<typename FluidState::Scalar, Evaluation>::value,
+                                                  Evaluation, Scalar>::type RhsEval;
+                // influx
+                LocalResidual::template evalPhaseFluxes_<RhsEval>(tmp,
+                                                                  phaseIdx,
+                                                                  insideIntQuants.pvtRegionIndex(),
+                                                                  extQuants,
+                                                                  fluidState);
+            }
+
+            for (unsigned i = 0; i < tmp.size(); ++i)
+                (*this)[i] += tmp[i];
+
+            // energy conservation
+            if (enableEnergy) {
+                Evaluation density;
+                Evaluation specificEnthalpy;
+                if (pBoundary > pInside) {
+                    if (focusDofIdx == interiorDofIdx) {
+                        density = fluidState.density(phaseIdx);
+                        specificEnthalpy = fluidState.enthalpy(phaseIdx);
+                    }
+                    else {
+                        density = Opm::getValue(fluidState.density(phaseIdx));
+                        specificEnthalpy = Opm::getValue(fluidState.enthalpy(phaseIdx));
+                    }
+                }
+                else if (focusDofIdx == interiorDofIdx) {
+                    density = insideIntQuants.fluidState().density(phaseIdx);
+                    specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
+                }
+                else {
+                    density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx));
+                    specificEnthalpy = Opm::getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
+                }
+
+                Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy;
+                EnergyModule::addToEnthalpyRate(*this, enthalpyRate);
+            }
+        }
+
+        if (enableSolvent) {
+            (*this)[Indices::contiSolventEqIdx] = extQuants.solventVolumeFlux();
+            if (blackoilConserveSurfaceVolume)
+                (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventInverseFormationVolumeFactor();
+            else
+                (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventDensity();
+
+        }
+
+        if (enablePolymer) {
+            (*this)[Indices::contiPolymerEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.polymerConcentration();
+        }
+
+        // make sure that the right mass conservation quantities are used
+        LocalResidual::adaptMassConservationQuantities_(*this, insideIntQuants.pvtRegionIndex());
+
+        // heat conduction
+        if (enableEnergy)
+            EnergyModule::addToEnthalpyRate(*this, extQuants.energyFlux());
+
+#ifndef NDEBUG
+        for (unsigned i = 0; i < numEq; ++i) {
+            Opm::Valgrind::CheckDefined((*this)[i]);
+        }
+        Opm::Valgrind::CheckDefined(*this);
+#endif
+    }
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::setInFlow
+     */
+    template <class Context, class FluidState>
+    void setInFlow(const Context& context,
+                   unsigned bfIdx,
+                   unsigned timeIdx,
+                   const FluidState& fluidState)
+    {
+        this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
+
+        // we only allow fluxes in the direction opposite to the outer
+        // unit normal
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+            Scalar& val = this->operator[](eqIdx);
+            val = std::min<Scalar>(0.0, val);
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::setOutFlow
+     */
+    template <class Context, class FluidState>
+    void setOutFlow(const Context& context,
+                    unsigned bfIdx,
+                    unsigned timeIdx,
+                    const FluidState& fluidState)
+    {
+        this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
+
+        // we only allow fluxes in the same direction as the outer
+        // unit normal
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+            Scalar& val = this->operator[](eqIdx);
+            val = std::max( Scalar(0), val);
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::setNoFlow
+     */
+    void setNoFlow()
+    { (*this) = Scalar(0); }
+
+    /*!
+     * \copydoc Specify an energy flux that corresponds to the thermal conduction from
+     *          the domain boundary
+     *
+     * This means that a "thermal flow" boundary is a no-flow condition for mass and thermal
+     * conduction for energy.
+     */
+    template <class Context, class FluidState>
+    void setThermalFlow(const Context& context,
+                        unsigned bfIdx,
+                        unsigned timeIdx,
+                        const FluidState& boundaryFluidState)
+    {
+        // set the mass no-flow condition
+        setNoFlow();
+
+        if (!enableEnergy)
+            // if we do not conserve energy there is nothing we should do in addition
+            return;
+
+        ExtensiveQuantities extQuants;
+        extQuants.updateBoundary(context, bfIdx, timeIdx, boundaryFluidState);
+
+        (*this)[contiEnergyEqIdx] += extQuants.energyFlux();
+
+#ifndef NDEBUG
+        for (unsigned i = 0; i < numEq; ++i)
+            Opm::Valgrind::CheckDefined((*this)[i]);
+        Opm::Valgrind::CheckDefined(*this);
+#endif
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoildarcyfluxmodule.hh

@@ -0,0 +1,100 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief This file contains the default flux module of the blackoil model.
+ *
+ * It is neccessary to accomodate the extensions of the black-oil model.
+ */
+#ifndef EWOMS_BLACK_OIL_DARCY_FLUX_MODULE_HH
+#define EWOMS_BLACK_OIL_DARCY_FLUX_MODULE_HH
+
+#include <opm/models/blackoil/blackoilproperties.hh>
+
+#include <opm/models/common/darcyfluxmodule.hh>
+#include <opm/models/utils/propertysystem.hh>
+
+namespace Opm {
+
+template <class TypeTag>
+class BlackOilDarcyExtensiveQuantities;
+
+/*!
+ * \ingroup FluxModules
+ * \brief Provides a Darcy flux module for the blackoil model
+ */
+template <class TypeTag>
+struct BlackOilDarcyFluxModule
+{
+    typedef DarcyIntensiveQuantities<TypeTag> FluxIntensiveQuantities;
+    typedef BlackOilDarcyExtensiveQuantities<TypeTag> FluxExtensiveQuantities;
+    typedef DarcyBaseProblem<TypeTag> FluxBaseProblem;
+
+    /*!
+     * \brief Register all run-time parameters for the flux module.
+     */
+    static void registerParameters()
+    { }
+};
+
+
+/*!
+ * \ingroup FluxModules
+ * \brief Specifies the extensive quantities for the black-oil model if using Darcy relation.
+ *
+ * This class basically forwards everything to the default Darcy flux module and adds a
+ * few methods needed by the extensions of the black-oil model. (i.e. the solvent and the
+ * polymer extensions.)
+ */
+template <class TypeTag>
+class BlackOilDarcyExtensiveQuantities : public DarcyExtensiveQuantities<TypeTag>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+public:
+    /*!
+     * \brief Update the extensive quantities which are specific to the solvent extension
+     * of the black-oil model.
+     */
+    void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    {
+        asImp_().updateVolumeFluxPerm(elemCtx,
+                                      scvfIdx,
+                                      timeIdx);
+
+    }
+
+    void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    { asImp_().updateShearMultipliersPerm(elemCtx, scvfIdx, timeIdx); }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilenergymodules.hh

@@ -0,0 +1,624 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief Contains the classes required to extend the black-oil model by energy.
+ */
+#ifndef EWOMS_BLACK_OIL_ENERGY_MODULE_HH
+#define EWOMS_BLACK_OIL_ENERGY_MODULE_HH
+
+#include "blackoilproperties.hh"
+#include <ewoms/io/vtkblackoilenergymodule.hh>
+#include <opm/models/common/quantitycallbacks.hh>
+
+#include <opm/material/common/Tabulated1DFunction.hpp>
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <string>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOil
+ * \brief Contains the high level supplements required to extend the black oil
+ *        model by energy.
+ */
+template <class TypeTag, bool enableEnergyV = GET_PROP_VALUE(TypeTag, EnableEnergy)>
+class BlackOilEnergyModule
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
+    static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
+
+    static constexpr unsigned enableEnergy = enableEnergyV;
+    static constexpr unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
+    static constexpr unsigned numPhases = FluidSystem::numPhases;
+
+public:
+    /*!
+     * \brief Register all run-time parameters for the black-oil energy module.
+     */
+    static void registerParameters()
+    {
+        if (!enableEnergy)
+            // energys have been disabled at compile time
+            return;
+
+        Opm::VtkBlackOilEnergyModule<TypeTag>::registerParameters();
+    }
+
+    /*!
+     * \brief Register all energy specific VTK and ECL output modules.
+     */
+    static void registerOutputModules(Model& model,
+                                      Simulator& simulator)
+    {
+        if (!enableEnergy)
+            // energys have been disabled at compile time
+            return;
+
+        model.addOutputModule(new Opm::VtkBlackOilEnergyModule<TypeTag>(simulator));
+    }
+
+    static bool primaryVarApplies(unsigned pvIdx)
+    {
+        if (!enableEnergy)
+            // energys have been disabled at compile time
+            return false;
+
+        return pvIdx == temperatureIdx;
+    }
+
+    static std::string primaryVarName(unsigned pvIdx OPM_OPTIM_UNUSED)
+    {
+        assert(primaryVarApplies(pvIdx));
+
+        return "temperature";
+    }
+
+    static Scalar primaryVarWeight(unsigned pvIdx OPM_OPTIM_UNUSED)
+    {
+        assert(primaryVarApplies(pvIdx));
+
+        // TODO: it may be beneficial to chose this differently.
+        return static_cast<Scalar>(1.0);
+    }
+
+    static bool eqApplies(unsigned eqIdx)
+    {
+        if (!enableEnergy)
+            return false;
+
+        return eqIdx == contiEnergyEqIdx;
+    }
+
+    static std::string eqName(unsigned eqIdx OPM_OPTIM_UNUSED)
+    {
+        assert(eqApplies(eqIdx));
+
+        return "conti^energy";
+    }
+
+    static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
+    {
+        assert(eqApplies(eqIdx));
+
+        return 1.0;
+    }
+
+    // must be called after water storage is computed
+    template <class LhsEval>
+    static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                           const IntensiveQuantities& intQuants)
+    {
+        if (!enableEnergy)
+            return;
+
+        const auto& poro = Opm::decay<LhsEval>(intQuants.porosity());
+
+        // accumulate the internal energy of the fluids
+        const auto& fs = intQuants.fluidState();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            const auto& u = Opm::decay<LhsEval>(fs.internalEnergy(phaseIdx));
+            const auto& S = Opm::decay<LhsEval>(fs.saturation(phaseIdx));
+            const auto& rho = Opm::decay<LhsEval>(fs.density(phaseIdx));
+
+            storage[contiEnergyEqIdx] += poro*S*u*rho;
+        }
+
+        // add the internal energy of the rock
+        Scalar refPoro = intQuants.referencePorosity();
+        const auto& uRock = Opm::decay<LhsEval>(intQuants.rockInternalEnergy());
+        storage[contiEnergyEqIdx] += (1.0 - refPoro)*uRock;
+        storage[contiEnergyEqIdx] *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor);
+    }
+
+    static void computeFlux(RateVector& flux,
+                            const ElementContext& elemCtx,
+                            unsigned scvfIdx,
+                            unsigned timeIdx)
+    {
+        if (!enableEnergy)
+            return;
+
+        flux[contiEnergyEqIdx] = 0.0;
+
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned focusIdx = elemCtx.focusDofIndex();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
+            if (upIdx == focusIdx)
+                addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
+            else
+                addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
+        }
+
+        // diffusive energy flux
+        flux[contiEnergyEqIdx] += extQuants.energyFlux();
+        flux[contiEnergyEqIdx] *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor);
+    }
+
+    template <class UpstreamEval>
+    static void addPhaseEnthalpyFlux_(RateVector& flux,
+                                      unsigned phaseIdx,
+                                      const ElementContext& elemCtx,
+                                      unsigned scvfIdx,
+                                      unsigned timeIdx)
+    {
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
+        const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+        const auto& fs = up.fluidState();
+
+        const auto& volFlux = extQuants.volumeFlux(phaseIdx);
+        flux[contiEnergyEqIdx] +=
+            Opm::decay<UpstreamEval>(fs.enthalpy(phaseIdx))
+            * Opm::decay<UpstreamEval>(fs.density(phaseIdx))
+            * volFlux;
+    }
+
+    static void addToEnthalpyRate(RateVector& flux,
+                                  const Evaluation& hRate)
+    {
+        if (!enableEnergy)
+            return;
+
+        flux[contiEnergyEqIdx] += hRate;
+    }
+
+    /*!
+     * \brief Assign the energy specific primary variables to a PrimaryVariables object
+     */
+    static void assignPrimaryVars(PrimaryVariables& priVars,
+                                  Scalar temperature OPM_UNUSED)
+    {
+        if (!enableEnergy)
+            return;
+
+        priVars[temperatureIdx] = temperatureIdx;
+    }
+
+    /*!
+     * \brief Assign the energy specific primary variables to a PrimaryVariables object
+     */
+    template <class FluidState>
+    static void assignPrimaryVars(PrimaryVariables& priVars,
+                                  const FluidState& fluidState)
+    {
+        if (!enableEnergy)
+            return;
+
+        priVars[temperatureIdx] = fluidState.temperature(/*phaseIdx=*/0);
+    }
+
+    /*!
+     * \brief Do a Newton-Raphson update the primary variables of the energys.
+     */
+    static void updatePrimaryVars(PrimaryVariables& newPv,
+                                  const PrimaryVariables& oldPv,
+                                  const EqVector& delta)
+    {
+        if (!enableEnergy)
+            return;
+
+        // do a plain unchopped Newton update
+        newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
+    }
+
+    /*!
+     * \brief Return how much a Newton-Raphson update is considered an error
+     */
+    static Scalar computeUpdateError(const PrimaryVariables& oldPv OPM_UNUSED,
+                                     const EqVector& delta OPM_UNUSED)
+    {
+        // do not consider consider the cange of energy primary variables for
+        // convergence
+        // TODO: maybe this should be changed
+        return static_cast<Scalar>(0.0);
+    }
+
+    /*!
+     * \brief Return how much a residual is considered an error
+     */
+    static Scalar computeResidualError(const EqVector& resid)
+    {
+        // do not weight the residual of energy when it comes to convergence
+        return std::abs(Opm::scalarValue(resid[contiEnergyEqIdx]));
+    }
+
+    template <class DofEntity>
+    static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
+    {
+        if (!enableEnergy)
+            return;
+
+        unsigned dofIdx = model.dofMapper().index(dof);
+        const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
+        outstream << priVars[temperatureIdx];
+    }
+
+    template <class DofEntity>
+    static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
+    {
+        if (!enableEnergy)
+            return;
+
+        unsigned dofIdx = model.dofMapper().index(dof);
+        PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
+        PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
+
+        instream >> priVars0[temperatureIdx];
+
+        // set the primary variables for the beginning of the current time step.
+        priVars1 = priVars0[temperatureIdx];
+    }
+};
+
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilEnergyIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the equations needed by the
+ *        energys extension of the black-oil model.
+ */
+template <class TypeTag, bool enableEnergyV = GET_PROP_VALUE(TypeTag, EnableEnergy)>
+class BlackOilEnergyIntensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, SolidEnergyLaw) SolidEnergyLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductionLaw) ThermalConductionLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    typedef BlackOilEnergyModule<TypeTag> EnergyModule;
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    static constexpr int temperatureIdx = Indices::temperatureIdx;
+    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+
+public:
+    /*!
+     * \brief Update the temperature of the intensive quantity's fluid state
+     *
+     */
+    void updateTemperature_(const ElementContext& elemCtx,
+                            unsigned dofIdx,
+                            unsigned timeIdx)
+    {
+        auto& fs = asImp_().fluidState_;
+        const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+
+        // set temperature
+        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx));
+    }
+
+    /*!
+     * \brief Compute the intensive quantities needed to handle energy conservation
+     *
+     */
+    void updateEnergyQuantities_(const ElementContext& elemCtx,
+                                 unsigned dofIdx,
+                                 unsigned timeIdx,
+                                 const typename FluidSystem::template ParameterCache<Evaluation>& paramCache)
+    {
+        auto& fs = asImp_().fluidState_;
+
+        // compute the specific enthalpy of the fluids, the specific enthalpy of the rock
+        // and the thermal condictivity coefficients
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                fs.setEnthalpy(phaseIdx, 0.0);
+                continue;
+            }
+
+            const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
+            fs.setEnthalpy(phaseIdx, h);
+        }
+
+        const auto& solidEnergyLawParams = elemCtx.problem().solidEnergyLawParams(elemCtx, dofIdx, timeIdx);
+        rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
+
+        const auto& thermalConductionLawParams = elemCtx.problem().thermalConductionLawParams(elemCtx, dofIdx, timeIdx);
+        totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
+    }
+
+    const Evaluation& rockInternalEnergy() const
+    { return rockInternalEnergy_; }
+
+    const Evaluation& totalThermalConductivity() const
+    { return totalThermalConductivity_; }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Evaluation rockInternalEnergy_;
+    Evaluation totalThermalConductivity_;
+};
+
+template <class TypeTag>
+class BlackOilEnergyIntensiveQuantities<TypeTag, false>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+    static constexpr bool enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature);
+
+public:
+    void updateTemperature_(const ElementContext& elemCtx,
+                            unsigned dofIdx,
+                            unsigned timeIdx)
+    {
+        if (enableTemperature) {
+            // even if energy is conserved, the temperature can vary over the spatial
+            // domain if the EnableTemperature property is set to true
+            auto& fs = asImp_().fluidState_;
+            Scalar T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
+            fs.setTemperature(T);
+        }
+    }
+
+    void updateEnergyQuantities_(const ElementContext& elemCtx OPM_UNUSED,
+                                 unsigned dofIdx OPM_UNUSED,
+                                 unsigned timeIdx OPM_UNUSED,
+                                 const typename FluidSystem::template ParameterCache<Evaluation>& paramCache OPM_UNUSED)
+    { }
+
+    const Evaluation& rockInternalEnergy() const
+    { throw std::logic_error("Requested the rock internal energy, which is "
+                             "unavailable because energy is not conserved"); }
+
+    const Evaluation& totalThermalConductivity() const
+    { throw std::logic_error("Requested the total thermal conductivity, which is "
+                             "unavailable because energy is not conserved"); }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+};
+
+
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilEnergyExtensiveQuantities
+ *
+ * \brief Provides the energy specific extensive quantities to the generic black-oil
+ *        module's extensive quantities.
+ */
+template <class TypeTag, bool enableEnergyV = GET_PROP_VALUE(TypeTag, EnableEnergy)>
+class BlackOilEnergyExtensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+
+    typedef BlackOilEnergyModule<TypeTag> EnergyModule;
+
+    static const int dimWorld = GridView::dimensionworld;
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+    typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
+
+public:
+    void updateEnergy(const ElementContext& elemCtx,
+                      unsigned scvfIdx,
+                      unsigned timeIdx)
+    {
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        Scalar faceArea = scvf.area();
+        unsigned inIdx = scvf.interiorIndex();
+        unsigned exIdx = scvf.exteriorIndex();
+        const auto& inIq = elemCtx.intensiveQuantities(inIdx, timeIdx);
+        const auto& exIq = elemCtx.intensiveQuantities(exIdx, timeIdx);
+        const auto& inFs = inIq.fluidState();
+        const auto& exFs = exIq.fluidState();
+
+        Evaluation deltaT;
+        if (elemCtx.focusDofIndex() == inIdx)
+            deltaT =
+                Opm::decay<Scalar>(exFs.temperature(/*phaseIdx=*/0))
+                - inFs.temperature(/*phaseIdx=*/0);
+        else if (elemCtx.focusDofIndex() == exIdx)
+            deltaT =
+                exFs.temperature(/*phaseIdx=*/0)
+                - Opm::decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
+        else
+            deltaT =
+                Opm::decay<Scalar>(exFs.temperature(/*phaseIdx=*/0))
+                - Opm::decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
+
+        Evaluation inLambda;
+        if (elemCtx.focusDofIndex() == inIdx)
+            inLambda = inIq.totalThermalConductivity();
+        else
+            inLambda = Opm::decay<Scalar>(inIq.totalThermalConductivity());
+
+        Evaluation exLambda;
+        if (elemCtx.focusDofIndex() == exIdx)
+            exLambda = exIq.totalThermalConductivity();
+        else
+            exLambda = Opm::decay<Scalar>(exIq.totalThermalConductivity());
+
+        auto distVec = elemCtx.pos(exIdx, timeIdx);
+        distVec -= elemCtx.pos(inIdx, timeIdx);
+
+        Evaluation H;
+        if (inLambda > 0.0 && exLambda > 0.0) {
+            // compute the "thermal transmissibility". In contrast to the normal
+            // transmissibility this cannot be done as a preprocessing step because the
+            // average thermal thermal conductivity is analogous to the permeability but
+            // depends on the solution.
+            Scalar alpha = elemCtx.problem().thermalHalfTransmissibility(elemCtx, scvfIdx, timeIdx);
+            const Evaluation& inH = inLambda*alpha;
+            const Evaluation& exH = exLambda*alpha;
+            H = 1.0/(1.0/inH + 1.0/exH);
+        }
+        else
+            H = 0.0;
+
+        energyFlux_ = deltaT * (-H/faceArea);
+    }
+
+    template <class Context, class BoundaryFluidState>
+    void updateEnergyBoundary(const Context& ctx,
+                              unsigned scvfIdx,
+                              unsigned timeIdx,
+                              const BoundaryFluidState& boundaryFs)
+    {
+        const auto& stencil = ctx.stencil(timeIdx);
+        const auto& scvf = stencil.boundaryFace(scvfIdx);
+
+        unsigned inIdx = scvf.interiorIndex();
+        const auto& inIq = ctx.intensiveQuantities(inIdx, timeIdx);
+        const auto& inFs = inIq.fluidState();
+
+        Evaluation deltaT;
+        if (ctx.focusDofIndex() == inIdx)
+            deltaT =
+                boundaryFs.temperature(/*phaseIdx=*/0)
+                - inFs.temperature(/*phaseIdx=*/0);
+        else
+            deltaT =
+                Opm::decay<Scalar>(boundaryFs.temperature(/*phaseIdx=*/0))
+                - Opm::decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
+
+        Evaluation lambda;
+        if (ctx.focusDofIndex() == inIdx)
+            lambda = inIq.totalThermalConductivity();
+        else
+            lambda = Opm::decay<Scalar>(inIq.totalThermalConductivity());
+
+        auto distVec = scvf.integrationPos();
+        distVec -= ctx.pos(inIdx, timeIdx);
+
+        if (lambda > 0.0) {
+            // compute the "thermal transmissibility". In contrast to the normal
+            // transmissibility this cannot be done as a preprocessing step because the
+            // average thermal conductivity is analogous to the permeability but depends
+            // on the solution.
+            Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
+            energyFlux_ = deltaT*lambda*(-alpha);
+        }
+        else
+            energyFlux_ = 0.0;
+    }
+
+    const Evaluation& energyFlux()  const
+    { return energyFlux_; }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Evaluation energyFlux_;
+};
+
+template <class TypeTag>
+class BlackOilEnergyExtensiveQuantities<TypeTag, false>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+
+public:
+    void updateEnergy(const ElementContext& elemCtx OPM_UNUSED,
+                      unsigned scvfIdx OPM_UNUSED,
+                      unsigned timeIdx OPM_UNUSED)
+    {}
+
+    template <class Context, class BoundaryFluidState>
+    void updateEnergyBoundary(const Context& ctx OPM_UNUSED,
+                              unsigned scvfIdx OPM_UNUSED,
+                              unsigned timeIdx OPM_UNUSED,
+                              const BoundaryFluidState& boundaryFs OPM_UNUSED)
+    {}
+
+    const Evaluation& energyFlux()  const
+    { throw std::logic_error("Requested the energy flux, but energy is not conserved"); }
+};
+
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilextensivequantities.hh

@@ -0,0 +1,103 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilExtensiveQuantities
+ */
+#ifndef EWOMS_BLACK_OIL_EXTENSIVE_QUANTITIES_HH
+#define EWOMS_BLACK_OIL_EXTENSIVE_QUANTITIES_HH
+
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilenergymodules.hh"
+
+#include <opm/models/common/multiphasebaseextensivequantities.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ * \ingroup ExtensiveQuantities
+ *
+ * \brief This template class contains the data which is required to
+ *        calculate the fluxes of the fluid phases over a face of a
+ *        finite volume for the black-oil model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the intergration point, etc.
+ */
+template <class TypeTag>
+class BlackOilExtensiveQuantities
+    : public MultiPhaseBaseExtensiveQuantities<TypeTag>
+    , public BlackOilSolventExtensiveQuantities<TypeTag>
+    , public BlackOilPolymerExtensiveQuantities<TypeTag>
+    , public BlackOilEnergyExtensiveQuantities<TypeTag>
+{
+    typedef MultiPhaseBaseExtensiveQuantities<TypeTag> MultiPhaseParent;
+
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+    /*!
+     * \brief Update the extensive quantities for a given sub-control-volume-face.
+     *
+     * \param elemCtx Reference to the current element context.
+     * \param scvfIdx The local index of the sub-control-volume face for
+     *                which the extensive quantities should be calculated.
+     * \param timeIdx The index used by the time discretization.
+     */
+    void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    {
+        MultiPhaseParent::update(elemCtx, scvfIdx, timeIdx);
+
+        asImp_().updateSolvent(elemCtx, scvfIdx, timeIdx);
+        asImp_().updatePolymer(elemCtx, scvfIdx, timeIdx);
+        asImp_().updateEnergy(elemCtx, scvfIdx, timeIdx);
+    }
+
+    template <class Context, class FluidState>
+    void updateBoundary(const Context& ctx,
+                        unsigned bfIdx,
+                        unsigned timeIdx,
+                        const FluidState& fluidState)
+    {
+        MultiPhaseParent::updateBoundary(ctx, bfIdx, timeIdx, fluidState);
+
+        asImp_().updateEnergyBoundary(ctx, bfIdx, timeIdx, fluidState);
+    }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilfoammodules.hh

@@ -0,0 +1,619 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief Contains the classes required to extend the black-oil model to include the effects of foam.
+ */
+#ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH
+#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
+
+#include "blackoilproperties.hh"
+//#include <ewoms/io/vtkblackoilfoammodule.hh>
+#include <opm/models/common/quantitycallbacks.hh>
+
+#include <opm/material/common/Tabulated1DFunction.hpp>
+//#include <opm/material/common/IntervalTabulated2DFunction.hpp>
+
+#if HAVE_ECL_INPUT
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/FoamadsTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/FoammobTable.hpp>
+#endif
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <string>
+#include <math.h>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOil
+ * \brief Contains the high level supplements required to extend the black oil
+ *        model to include the effects of foam.
+ */
+template <class TypeTag, bool enableFoamV = GET_PROP_VALUE(TypeTag, EnableFoam)>
+class BlackOilFoamModule
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+
+    typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
+
+    static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
+    static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
+    static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
+
+    static constexpr unsigned enableFoam = enableFoamV;
+    static constexpr bool enableVtkOutput = GET_PROP_VALUE(TypeTag, EnableVtkOutput);
+
+    static constexpr unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
+    static constexpr unsigned numPhases = FluidSystem::numPhases;
+
+public:
+    // a struct containing constants to calculate change to relative permeability,
+    // based on model (1-9) in Table 1 of
+    // Kun Ma, Guangwei Ren, Khalid Mateen, Danielle Morel, and Philippe Cordelier:
+    // "Modeling techniques for foam flow in porous media", SPE Journal, 20(03):453–470, jun 2015.
+    // The constants are provided by various deck keywords as shown in the comments below.
+    struct FoamCoefficients {
+        Scalar fm_min = 1e-20;   // FOAMFSC
+        Scalar fm_mob = 1.0;     // FOAMFRM
+
+        Scalar fm_surf = 1.0;    // FOAMFSC
+        Scalar ep_surf = 1.0;    // FOAMFSC
+
+        Scalar fm_oil = 1.0;     // FOAMFSO
+        Scalar fl_oil = 0.0;     // FOAMFSO
+        Scalar ep_oil = 0.0;     // FOAMFSO
+
+        Scalar fm_cap = 1.0;     // FOAMFCN
+        Scalar ep_cap = 0.0;     // FOAMFCN
+
+        Scalar fm_dry = 1.0;     // FOAMFSW
+        Scalar ep_dry = 0.0;     // FOAMFSW
+    };
+
+#if HAVE_ECL_INPUT
+    /*!
+     * \brief Initialize all internal data structures needed by the foam module
+     */
+    static void initFromDeck(const Opm::Deck& deck, const Opm::EclipseState& eclState)
+    {
+        // some sanity checks: if foam is enabled, the FOAM keyword must be
+        // present, if foam is disabled the keyword must not be present.
+        if (enableFoam && !deck.hasKeyword("FOAM")) {
+            throw std::runtime_error("Non-trivial foam treatment requested at compile time, but "
+                                     "the deck does not contain the FOAM keyword");
+        }
+        else if (!enableFoam && deck.hasKeyword("FOAM")) {
+            throw std::runtime_error("Foam treatment disabled at compile time, but the deck "
+                                     "contains the FOAM keyword");
+        }
+
+        if (!deck.hasKeyword("FOAM")) {
+            return; // foam treatment is supposed to be disabled
+        }
+
+        // Check that only implemented options are used.
+        // We only support the default values of FOAMOPTS (GAS, TAB).
+        if (deck.hasKeyword("FOAMOPTS")) {
+            const auto kw = deck.getKeyword("FOAMOPTS");
+            if (kw.getRecord(0).getItem("TRANSPORT_PHASE").get<std::string>(0) != "GAS") {
+                throw std::runtime_error("In FOAMOPTS, only GAS is allowed for the transport phase.");
+            }
+            if (kw.getRecord(0).getItem("MODEL").get<std::string>(0) != "TAB") {
+                throw std::runtime_error("In FOAMOPTS, only TAB is allowed for the gas mobility factor reduction model.");
+            }
+        }
+
+        const auto& tableManager = eclState.getTableManager();
+        const unsigned int numSatRegions = tableManager.getTabdims().getNumSatTables();
+        setNumSatRegions(numSatRegions);
+        const unsigned int numPvtRegions = tableManager.getTabdims().getNumPVTTables();
+        setNumPvtRegions(numPvtRegions);
+
+        // Get and check FOAMROCK data.
+        const Opm::FoamConfig& foamConf = eclState.getInitConfig().getFoamConfig();
+        if (numSatRegions != foamConf.size()) {
+            throw std::runtime_error("Inconsistent sizes, number of saturation regions differ from the number of elements "
+                                     "in FoamConfig, which typically corresponds to the number of records in FOAMROCK.");
+        }
+
+        // Get and check FOAMADS data.
+        const auto& foamadsTables = tableManager.getFoamadsTables();
+        if (foamadsTables.empty()) {
+            throw std::runtime_error("FOAMADS must be specified in FOAM runs");
+        }
+        if (numSatRegions != foamadsTables.size()) {
+            throw std::runtime_error("Inconsistent sizes, number of saturation regions differ from the "
+                                     "number of FOAMADS tables.");
+        }
+
+        // Set data that vary with saturation region.
+        for (std::size_t satReg = 0; satReg < numSatRegions; ++satReg) {
+            const auto& rec = foamConf.getRecord(satReg);
+            foamCoefficients_[satReg] = FoamCoefficients();
+            foamCoefficients_[satReg].fm_min = rec.minimumSurfactantConcentration();
+            foamCoefficients_[satReg].fm_surf = rec.referenceSurfactantConcentration();
+            foamCoefficients_[satReg].ep_surf = rec.exponent();
+            foamRockDensity_[satReg] = rec.rockDensity();
+            foamAllowDesorption_[satReg] = rec.allowDesorption();
+            const auto& foamadsTable = foamadsTables.template getTable<Opm::FoamadsTable>(satReg);
+            const auto& conc = foamadsTable.getFoamConcentrationColumn();
+            const auto& ads = foamadsTable.getAdsorbedFoamColumn();
+            adsorbedFoamTable_[satReg].setXYContainers(conc, ads);
+        }
+
+        // Get and check FOAMMOB data.
+        const auto& foammobTables = tableManager.getFoammobTables();
+        if (foammobTables.empty()) {
+            // When in the future adding support for the functional
+            // model, FOAMMOB will not be required anymore (functional
+            // family of keywords can be used instead, FOAMFSC etc.).
+            throw std::runtime_error("FOAMMOB must be specified in FOAM runs");
+        }
+        if (numPvtRegions != foammobTables.size()) {
+            throw std::runtime_error("Inconsistent sizes, number of PVT regions differ from the "
+                                     "number of FOAMMOB tables.");
+        }
+
+        // Set data that vary with PVT region.
+        for (std::size_t pvtReg = 0; pvtReg < numPvtRegions; ++pvtReg) {
+            const auto& foammobTable = foammobTables.template getTable<Opm::FoammobTable>(pvtReg);
+            const auto& conc = foammobTable.getFoamConcentrationColumn();
+            const auto& mobMult = foammobTable.getMobilityMultiplierColumn();
+            gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult);
+        }
+    }
+#endif
+
+    /*!
+     * \brief Specify the number of saturation regions.
+     */
+    static void setNumSatRegions(unsigned numRegions)
+    {
+        foamCoefficients_.resize(numRegions);
+        foamRockDensity_.resize(numRegions);
+        foamAllowDesorption_.resize(numRegions);
+        adsorbedFoamTable_.resize(numRegions);
+    }
+
+
+    /*!
+     * \brief Specify the number of PVT regions.
+     */
+    static void setNumPvtRegions(unsigned numRegions)
+    {
+        gasMobilityMultiplierTable_.resize(numRegions);
+    }
+
+
+    /*!
+     * \brief Register all run-time parameters for the black-oil foam module.
+     */
+    static void registerParameters()
+    {
+        if (!enableFoam)
+            // foam has been disabled at compile time
+            return;
+
+        //Opm::VtkBlackOilFoamModule<TypeTag>::registerParameters();
+    }
+
+    /*!
+     * \brief Register all foam specific VTK and ECL output modules.
+     */
+    static void registerOutputModules(Model& model,
+                                      Simulator& simulator)
+    {
+        if (!enableFoam)
+            // foam have been disabled at compile time
+            return;
+
+        if (enableVtkOutput) {
+            Opm::OpmLog::warning("VTK output requested, currently unsupported by the foam module.");
+        }
+        //model.addOutputModule(new Opm::VtkBlackOilFoamModule<TypeTag>(simulator));
+    }
+
+    static bool primaryVarApplies(unsigned pvIdx)
+    {
+        if (!enableFoam) {
+            return false;
+        } else {
+            return pvIdx == foamConcentrationIdx;
+        }
+    }
+
+    static std::string primaryVarName(unsigned pvIdx OPM_OPTIM_UNUSED)
+    {
+        assert(primaryVarApplies(pvIdx));
+        return "foam_concentration";
+    }
+
+    static Scalar primaryVarWeight(unsigned pvIdx OPM_OPTIM_UNUSED)
+    {
+       assert(primaryVarApplies(pvIdx));
+
+       // TODO: it may be beneficial to chose this differently.
+       return static_cast<Scalar>(1.0);
+    }
+
+    static bool eqApplies(unsigned eqIdx)
+    {
+        if (!enableFoam)
+            return false;
+
+        return eqIdx == contiFoamEqIdx;
+
+    }
+
+    static std::string eqName(unsigned eqIdx)
+    {
+        assert(eqApplies(eqIdx));
+
+        return "conti^foam";
+    }
+
+    static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
+    {
+       assert(eqApplies(eqIdx));
+
+       // TODO: it may be beneficial to chose this differently.
+       return static_cast<Scalar>(1.0);
+    }
+
+    // must be called after water storage is computed
+    template <class LhsEval>
+    static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                           const IntensiveQuantities& intQuants)
+    {
+        if (!enableFoam)
+            return;
+
+        const auto& fs = intQuants.fluidState();
+
+        LhsEval surfaceVolumeFreeGas =
+            Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
+            * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx))
+            * Toolbox::template decay<LhsEval>(intQuants.porosity());
+
+        // Avoid singular matrix if no gas is present.
+        surfaceVolumeFreeGas = Opm::max(surfaceVolumeFreeGas, 1e-10);
+
+        // Foam/surfactant in gas phase.
+        const LhsEval gasFoam = surfaceVolumeFreeGas
+            * Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
+
+        // Adsorbed foam/surfactant.
+        const LhsEval adsorbedFoam =
+            Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
+            * Toolbox::template decay<LhsEval>(intQuants.foamRockDensity())
+            * Toolbox::template decay<LhsEval>(intQuants.foamAdsorbed());
+
+        LhsEval accumulationFoam = gasFoam + adsorbedFoam;
+        storage[contiFoamEqIdx] += accumulationFoam;
+    }
+
+    static void computeFlux(RateVector& flux,
+                            const ElementContext& elemCtx,
+                            unsigned scvfIdx,
+                            unsigned timeIdx)
+
+    {
+        if (!enableFoam) {
+            return;
+        }
+
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
+        const unsigned inIdx = extQuants.interiorIndex();
+        const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+
+        // The effect of the gas mobility reduction factor is
+        // incorporated in the mobility, so the oil (if vaporized oil
+        // is active) and gas fluxes do not need modification here.
+        if (upIdx == inIdx) {
+            flux[contiFoamEqIdx] =
+                extQuants.volumeFlux(gasPhaseIdx)
+                *up.fluidState().invB(gasPhaseIdx)
+                *up.foamConcentration();
+        } else {
+            flux[contiFoamEqIdx] =
+                extQuants.volumeFlux(gasPhaseIdx)
+                *Opm::decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
+                *Opm::decay<Scalar>(up.foamConcentration());
+        }
+
+    }
+
+    /*!
+     * \brief Return how much a Newton-Raphson update is considered an error
+     */
+    static Scalar computeUpdateError(const PrimaryVariables& oldPv OPM_UNUSED,
+                                     const EqVector& delta OPM_UNUSED)
+    {
+        // do not consider the change of foam primary variables for convergence
+        // TODO: maybe this should be changed
+        return static_cast<Scalar>(0.0);
+    }
+
+    template <class DofEntity>
+    static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
+    {
+        if (!enableFoam)
+            return;
+
+        unsigned dofIdx = model.dofMapper().index(dof);
+        const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
+        outstream << priVars[foamConcentrationIdx];
+    }
+
+    template <class DofEntity>
+    static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
+    {
+        if (!enableFoam)
+            return;
+
+        unsigned dofIdx = model.dofMapper().index(dof);
+        PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
+        PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
+
+        instream >> priVars0[foamConcentrationIdx];
+
+        // set the primary variables for the beginning of the current time step.
+        priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
+    }
+
+    static const Scalar foamRockDensity(const ElementContext& elemCtx,
+                                        unsigned scvIdx,
+                                        unsigned timeIdx)
+    {
+        unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
+        return foamRockDensity_[satnumRegionIdx];
+    }
+
+    static bool foamAllowDesorption(const ElementContext& elemCtx,
+                                    unsigned scvIdx,
+                                    unsigned timeIdx)
+    {
+        unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
+        return foamAllowDesorption_[satnumRegionIdx];
+    }
+
+    static const TabulatedFunction& adsorbedFoamTable(const ElementContext& elemCtx,
+                                                      unsigned scvIdx,
+                                                      unsigned timeIdx)
+    {
+       unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
+       return adsorbedFoamTable_[satnumRegionIdx];
+    }
+
+    static const TabulatedFunction& gasMobilityMultiplierTable(const ElementContext& elemCtx,
+                                                               unsigned scvIdx,
+                                                               unsigned timeIdx)
+    {
+       unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
+       return gasMobilityMultiplierTable_[pvtnumRegionIdx];
+    }
+
+    static const FoamCoefficients& foamCoefficients(const ElementContext& elemCtx,
+                                                    const unsigned scvIdx,
+                                                    const unsigned timeIdx)
+    {
+        unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
+        return foamCoefficients_[satnumRegionIdx];
+    }
+
+private:
+    static std::vector<Scalar> foamRockDensity_;
+    static std::vector<bool> foamAllowDesorption_;
+    static std::vector<FoamCoefficients> foamCoefficients_;
+    static std::vector<TabulatedFunction> adsorbedFoamTable_;
+    static std::vector<TabulatedFunction> gasMobilityMultiplierTable_;
+};
+
+
+
+template <class TypeTag, bool enableFoam>
+std::vector<typename BlackOilFoamModule<TypeTag, enableFoam>::Scalar>
+BlackOilFoamModule<TypeTag, enableFoam>::foamRockDensity_;
+
+template <class TypeTag, bool enableFoam>
+std::vector<bool>
+BlackOilFoamModule<TypeTag, enableFoam>::foamAllowDesorption_;
+
+template <class TypeTag, bool enableFoam>
+std::vector<typename BlackOilFoamModule<TypeTag, enableFoam>::FoamCoefficients>
+BlackOilFoamModule<TypeTag, enableFoam>::foamCoefficients_;
+
+template <class TypeTag, bool enableFoam>
+std::vector<typename BlackOilFoamModule<TypeTag, enableFoam>::TabulatedFunction>
+BlackOilFoamModule<TypeTag, enableFoam>::adsorbedFoamTable_;
+
+template <class TypeTag, bool enableFoam>
+std::vector<typename BlackOilFoamModule<TypeTag, enableFoam>::TabulatedFunction>
+BlackOilFoamModule<TypeTag, enableFoam>::gasMobilityMultiplierTable_;
+
+
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilFoamIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the equations needed by the
+ *        polymers extension of the black-oil model.
+ */
+template <class TypeTag, bool enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam)>
+class BlackOilFoamIntensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    typedef BlackOilFoamModule<TypeTag> FoamModule;
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
+    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
+
+public:
+
+    /*!
+     * \brief Update the intensive properties needed to handle polymers from the
+     *        primary variables
+     *
+     */
+    void foamPropertiesUpdate_(const ElementContext& elemCtx,
+                               unsigned dofIdx,
+                               unsigned timeIdx)
+    {
+        const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        foamConcentration_ = priVars.makeEvaluation(foamConcentrationIdx, timeIdx);
+        const auto& fs = asImp_().fluidState_;
+
+        // Compute gas mobility reduction factor
+        Evaluation mobilityReductionFactor = 1.0;
+        if (false) {
+            // The functional model is used.
+            // TODO: allow this model.
+            // In order to do this we must allow transport to be in the water phase, not just the gas phase.
+            const auto& foamCoefficients = FoamModule::foamCoefficients(elemCtx, dofIdx, timeIdx);
+
+            const Scalar fm_mob = foamCoefficients.fm_mob;
+
+            const Scalar fm_surf = foamCoefficients.fm_surf;
+            const Scalar ep_surf = foamCoefficients.ep_surf;
+
+            const Scalar fm_oil = foamCoefficients.fm_oil;
+            const Scalar fl_oil = foamCoefficients.fl_oil;
+            const Scalar ep_oil = foamCoefficients.ep_oil;
+
+            const Scalar fm_dry = foamCoefficients.fm_dry;
+            const Scalar ep_dry = foamCoefficients.ep_dry;
+
+            const Scalar fm_cap = foamCoefficients.fm_cap;
+            const Scalar ep_cap = foamCoefficients.ep_cap;
+
+            const Evaluation C_surf = foamConcentration_;
+            const Evaluation Ca = 1e10; // TODO: replace with proper capillary number.
+            const Evaluation S_o = fs.saturation(oilPhaseIdx);
+            const Evaluation S_w = fs.saturation(waterPhaseIdx);
+
+            Evaluation F1 = pow(C_surf/fm_surf, ep_surf);
+            Evaluation F2 = pow((fm_oil-S_o)/(fm_oil-fl_oil), ep_oil);
+            Evaluation F3 = pow(fm_cap/Ca, ep_cap);
+            Evaluation F7 = 0.5 + atan(ep_dry*(S_w-fm_dry))/M_PI;
+
+            mobilityReductionFactor = 1./(1. + fm_mob*F1*F2*F3*F7);
+        } else {
+            // The tabular model is used.
+            // Note that the current implementation only includes the effect of foam concentration (FOAMMOB),
+            // and not the optional pressure dependence (FOAMMOBP) or shear dependence (FOAMMOBS).
+            const auto& gasMobilityMultiplier = FoamModule::gasMobilityMultiplierTable(elemCtx, dofIdx, timeIdx);
+            mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_, /* extrapolate = */ true);
+        }
+
+        // adjust gas mobility
+        asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
+
+        foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
+
+        const auto& adsorbedFoamTable = FoamModule::adsorbedFoamTable(elemCtx, dofIdx, timeIdx);
+        foamAdsorbed_ = adsorbedFoamTable.eval(foamConcentration_, /*extrapolate=*/true);
+        if (!FoamModule::foamAllowDesorption(elemCtx, dofIdx, timeIdx)) {
+            throw std::runtime_error("Foam module does not support the 'no desorption' option.");
+        }
+    }
+
+    const Evaluation& foamConcentration() const
+    { return foamConcentration_; }
+
+    Scalar foamRockDensity() const
+    { return foamRockDensity_; }
+
+    const Evaluation& foamAdsorbed() const
+    { return foamAdsorbed_; }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Evaluation foamConcentration_;
+    Scalar foamRockDensity_;
+    Evaluation foamAdsorbed_;
+};
+
+template <class TypeTag>
+class BlackOilFoamIntensiveQuantities<TypeTag, false>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+public:
+    void foamPropertiesUpdate_(const ElementContext& elemCtx OPM_UNUSED,
+                                  unsigned scvIdx OPM_UNUSED,
+                                  unsigned timeIdx OPM_UNUSED)
+    { }
+
+
+    const Evaluation& foamConcentration() const
+    { throw std::runtime_error("foamConcentration() called but foam is disabled"); }
+
+    Scalar foamRockDensity() const
+    { throw std::runtime_error("foamRockDensity() called but foam is disabled"); }
+
+    Scalar foamAdsorbed() const
+    { throw std::runtime_error("foamAdsorbed() called but foam is disabled"); }
+};
+
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilindices.hh

@@ -0,0 +1,153 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilIndices
+ */
+#ifndef EWOMS_BLACK_OIL_INDICES_HH
+#define EWOMS_BLACK_OIL_INDICES_HH
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief The primary variable and equation indices for the black-oil model.
+ */
+template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, unsigned PVOffset>
+struct BlackOilIndices
+{
+    //! Number of phases active at all times
+    static const int numPhases = 3;
+
+    //! All phases are enabled
+    static const bool oilEnabled = true;
+    static const bool waterEnabled = true;
+    static const bool gasEnabled = true;
+
+    //! Are solvents involved?
+    static const bool enableSolvent = numSolventsV > 0;
+
+    //! Are polymers involved?
+    static const bool enablePolymer = numPolymersV > 0;
+
+    //! Shall energy be conserved?
+    static const bool enableEnergy = numEnergyV > 0;
+
+    //! Number of solvent components to be considered
+    static const int numSolvents = enableSolvent ? numSolventsV : 0;
+
+    //! Number of polymer components to be considered
+    static const int numPolymers = enablePolymer ? numPolymersV : 0;
+
+    //! Number of energy equations to be considered
+    static const int numEnergy = enableEnergy ? numEnergyV : 0;
+
+    //! Number of foam equations to be considered
+    static const int numFoam = enableFoam? 1 : 0;
+
+    //! The number of equations
+    static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam;
+
+    //! \brief returns the index of "active" component
+    static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx)
+    { return compIdx; }
+
+    static constexpr unsigned activeToCanonicalComponentIndex(unsigned compIdx)
+    { return compIdx; }
+
+    ////////
+    // Primary variable indices
+    ////////
+
+    //! The index of the water saturation
+    static const int waterSaturationIdx = PVOffset + 0;
+
+    //! Index of the oil pressure in a vector of primary variables
+    static const int pressureSwitchIdx = PVOffset + 1;
+
+    /*!
+     * \brief Index of the switching variable which determines the composition of the
+     *        hydrocarbon phases.
+     *
+     * Depending on the phases present, this variable is either interpreted as the
+     * saturation of the gas phase, as the mole fraction of the gas component in the oil
+     * phase or as the mole fraction of the oil component in the gas phase.
+     */
+    static const int compositionSwitchIdx = PVOffset + 2;
+
+    //! Index of the primary variable for the first solvent
+    static const int solventSaturationIdx =
+        enableSolvent ? PVOffset + numPhases : -1000;
+
+    //! Index of the primary variable for the first polymer
+    static const int polymerConcentrationIdx =
+        enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the primary variable for the second polymer primary variable (molecular weight)
+    static const int polymerMoleWeightIdx =
+        numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the foam
+    static const int foamConcentrationIdx =
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+
+    //! Index of the primary variable for temperature
+    static const int temperatureIdx  =
+        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : - 1000;
+
+
+    ////////
+    // Equation indices
+    ////////
+
+    //! Index of the continuity equation of the first phase
+    static const int conti0EqIdx = PVOffset + 0;
+    // two continuity equations follow
+
+    //! Index of the continuity equation for the first solvent component
+    static const int contiSolventEqIdx =
+        enableSolvent ? PVOffset + numPhases : -1000;
+
+    //! Index of the continuity equation for the first polymer component
+    static const int contiPolymerEqIdx =
+        enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the continuity equation for the second polymer component (molecular weight)
+    static const int contiPolymerMWEqIdx =
+        numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the foam component
+    static const int contiFoamEqIdx =
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+
+
+    //! Index of the continuity equation for energy
+    static const int contiEnergyEqIdx =
+        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilintensivequantities.hh

@@ -0,0 +1,445 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilIntensiveQuantities
+ */
+#ifndef EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
+#define EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
+
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilfoammodules.hh"
+#include "blackoilenergymodules.hh"
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/common/fmatrix.hh>
+
+#include <cstring>
+#include <utility>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOilModel
+ * \ingroup IntensiveQuantities
+ *
+ * \brief Contains the quantities which are are constant within a
+ *        finite volume in the black-oil model.
+ */
+template <class TypeTag>
+class BlackOilIntensiveQuantities
+    : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
+    , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
+    , public BlackOilSolventIntensiveQuantities<TypeTag>
+    , public BlackOilPolymerIntensiveQuantities<TypeTag>
+    , public BlackOilFoamIntensiveQuantities<TypeTag>
+    , public BlackOilEnergyIntensiveQuantities<TypeTag>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam) };
+    enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { enableExperiments = GET_PROP_VALUE(TypeTag, EnableExperiments) };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { dimWorld = GridView::dimensionworld };
+    enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+
+    static const bool compositionSwitchEnabled = Indices::gasEnabled;
+    static const bool waterEnabled = Indices::waterEnabled;
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
+    typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
+    typedef Opm::BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled,  Indices::numPhases > FluidState;
+
+public:
+    BlackOilIntensiveQuantities()
+    {
+        if (compositionSwitchEnabled) {
+            fluidState_.setRs(0.0);
+            fluidState_.setRv(0.0);
+        }
+    }
+
+    BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) = default;
+
+    BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
+
+    /*!
+     * \copydoc IntensiveQuantities::update
+     */
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    {
+        ParentType::update(elemCtx, dofIdx, timeIdx);
+
+        const auto& problem = elemCtx.problem();
+        const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+
+        asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
+
+        unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+        unsigned pvtRegionIdx = priVars.pvtRegionIndex();
+        fluidState_.setPvtRegionIndex(pvtRegionIdx);
+
+        // extract the water and the gas saturations for convenience
+        Evaluation Sw = 0.0;
+        if (waterEnabled)
+            Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
+
+        Evaluation Sg = 0.0;
+        if (compositionSwitchEnabled)
+        {
+            if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
+                // -> threephase case
+                Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
+                // -> gas-water case
+                Sg = 1.0 - Sw;
+
+                // deal with solvent
+                if (enableSolvent)
+                    Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+            }
+            else
+            {
+                assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs);
+                // -> oil-water case
+                Sg = 0.0;
+            }
+        }
+
+        Opm::Valgrind::CheckDefined(Sg);
+        Opm::Valgrind::CheckDefined(Sw);
+
+        Evaluation So = 1.0 - Sw - Sg;
+
+        // deal with solvent
+        if (enableSolvent)
+            So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx))
+            fluidState_.setSaturation(waterPhaseIdx, Sw);
+
+        if (FluidSystem::phaseIsActive(gasPhaseIdx))
+            fluidState_.setSaturation(gasPhaseIdx, Sg);
+
+        if (FluidSystem::phaseIsActive(oilPhaseIdx))
+            fluidState_.setSaturation(oilPhaseIdx, So);
+
+        asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+
+        // now we compute all phase pressures
+        Evaluation pC[numPhases];
+        const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
+        MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
+
+        //oil is the reference phase for pressure
+        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
+            const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                if (FluidSystem::phaseIsActive(phaseIdx))
+                    fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
+        }
+
+        else {
+            const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                if (FluidSystem::phaseIsActive(phaseIdx))
+                    fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
+        }
+
+        // calculate relative permeabilities. note that we store the result into the
+        // mobility_ class attribute. the division by the phase viscosity happens later.
+        MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_);
+        Opm::Valgrind::CheckDefined(mobility_);
+
+        // update the Saturation functions for the blackoil solvent module.
+        asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+
+        const Evaluation& SoMax =
+            Opm::max(fluidState_.saturation(oilPhaseIdx),
+                     elemCtx.problem().maxOilSaturation(globalSpaceIdx));
+
+        // take the meaning of the switiching primary variable into account for the gas
+        // and oil phase compositions
+        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
+            // in the threephase case, gas and oil phases are potentially present, i.e.,
+            // we use the compositions of the gas-saturated oil and oil-saturated gas.
+            if (FluidSystem::enableDissolvedGas()) {
+                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
+                const Evaluation& RsSat =
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            oilPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+                fluidState_.setRs(Opm::min(RsMax, RsSat));
+            }
+            else if (compositionSwitchEnabled)
+                fluidState_.setRs(0.0);
+
+            if (FluidSystem::enableVaporizedOil()) {
+                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
+                const Evaluation& RvSat =
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+                fluidState_.setRv(Opm::min(RvMax, RvSat));
+            }
+            else if (compositionSwitchEnabled)
+                fluidState_.setRv(0.0);
+        }
+        else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
+            // if the switching variable is the mole fraction of the gas component in the
+            Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
+
+            // oil phase, we can directly set the composition of the oil phase
+            const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            fluidState_.setRs(Opm::min(RsMax, Rs));
+
+            if (FluidSystem::enableVaporizedOil()) {
+                // the gas phase is not present, but we need to compute its "composition"
+                // for the gravity correction anyway
+                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
+                const auto& RvSat =
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+
+                fluidState_.setRv(Opm::min(RvMax, RvSat));
+            }
+            else
+                fluidState_.setRv(0.0);
+        }
+        else {
+            assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
+
+            const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            fluidState_.setRv(Rv);
+
+            if (FluidSystem::enableDissolvedGas()) {
+                // the oil phase is not present, but we need to compute its "composition" for
+                // the gravity correction anyway
+                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
+                const auto& RsSat =
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            oilPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+
+                fluidState_.setRs(Opm::min(RsMax, RsSat));
+            }
+            else
+                fluidState_.setRs(0.0);
+        }
+
+        typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+        paramCache.setRegionIndex(pvtRegionIdx);
+        paramCache.setMaxOilSat(SoMax);
+        paramCache.updateAll(fluidState_);
+
+        // compute the phase densities and transform the phase permeabilities into mobilities
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
+            fluidState_.setInvB(phaseIdx, b);
+
+            const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
+            mobility_[phaseIdx] /= mu;
+        }
+        Opm::Valgrind::CheckDefined(mobility_);
+
+        // calculate the phase densities
+        Evaluation rho;
+        if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
+            rho = fluidState_.invB(waterPhaseIdx);
+            rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            fluidState_.setDensity(waterPhaseIdx, rho);
+        }
+
+        if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+            rho = fluidState_.invB(gasPhaseIdx);
+            rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            if (FluidSystem::enableVaporizedOil()) {
+                rho +=
+                    fluidState_.invB(gasPhaseIdx) *
+                    fluidState_.Rv() *
+                    FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            }
+            fluidState_.setDensity(gasPhaseIdx, rho);
+        }
+
+        if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+            rho = fluidState_.invB(oilPhaseIdx);
+            rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            if (FluidSystem::enableDissolvedGas()) {
+                rho +=
+                    fluidState_.invB(oilPhaseIdx) *
+                    fluidState_.Rs() *
+                    FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            }
+            fluidState_.setDensity(oilPhaseIdx, rho);
+        }
+
+        // retrieve the porosity from the problem
+        referencePorosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
+        porosity_ = referencePorosity_;
+
+        // the porosity must be modified by the compressibility of the
+        // rock...
+        Scalar rockCompressibility = problem.rockCompressibility(elemCtx, dofIdx, timeIdx);
+        if (rockCompressibility > 0.0) {
+            Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx);
+            Evaluation x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
+            porosity_ *= 1.0 + x + 0.5*x*x;
+        }
+
+        if (enableExperiments)
+            // deal with water induced rock compaction
+            porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
+
+        asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
+        asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+
+        // update the quantities which are required by the chosen
+        // velocity model
+        FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
+
+#ifndef NDEBUG
+        // some safety checks in debug mode
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            assert(Opm::isfinite(fluidState_.density(phaseIdx)));
+            assert(Opm::isfinite(fluidState_.saturation(phaseIdx)));
+            assert(Opm::isfinite(fluidState_.temperature(phaseIdx)));
+            assert(Opm::isfinite(fluidState_.pressure(phaseIdx)));
+            assert(Opm::isfinite(fluidState_.invB(phaseIdx)));
+        }
+        assert(Opm::isfinite(fluidState_.Rs()));
+        assert(Opm::isfinite(fluidState_.Rv()));
+#endif
+    }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::fluidState
+     */
+    const FluidState& fluidState() const
+    { return fluidState_; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::mobility
+     */
+    const Evaluation& mobility(unsigned phaseIdx) const
+    { return mobility_[phaseIdx]; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::porosity
+     */
+    const Evaluation& porosity() const
+    { return porosity_; }
+
+    /*!
+     * \brief Returns the index of the PVT region used to calculate the thermodynamic
+     *        quantities.
+     *
+     * This allows to specify different Pressure-Volume-Temperature (PVT) relations in
+     * different parts of the spatial domain. Note that this concept should be seen as a
+     * work-around of the fact that the black-oil model does not capture the
+     * thermodynamics well enough. (Because there is, err, only a single real world with
+     * in which all substances follow the same physical laws and hence the same
+     * thermodynamics.) Anyway: Since the ECL file format uses multiple PVT regions, we
+     * support it as well in our black-oil model. (Note that, if it is not explicitly
+     * specified, the PVT region index is 0.)
+     */
+    auto pvtRegionIndex() const
+        -> decltype(std::declval<FluidState>().pvtRegionIndex())
+    { return fluidState_.pvtRegionIndex(); }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::relativePermeability
+     */
+    Evaluation relativePermeability(unsigned phaseIdx) const
+    {
+        // warning: slow
+        return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
+    }
+
+    /*!
+     * \brief Returns the porosity of the rock at reference conditions.
+     *
+     * I.e., the porosity of rock which is not perturbed by pressure and temperature
+     * changes.
+     */
+    Scalar referencePorosity() const
+    { return referencePorosity_; }
+
+private:
+    friend BlackOilSolventIntensiveQuantities<TypeTag>;
+    friend BlackOilPolymerIntensiveQuantities<TypeTag>;
+    friend BlackOilEnergyIntensiveQuantities<TypeTag>;
+    friend BlackOilFoamIntensiveQuantities<TypeTag>;
+
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    FluidState fluidState_;
+    Scalar referencePorosity_;
+    Evaluation porosity_;
+    Evaluation mobility_[numPhases];
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoillocalresidual.hh

@@ -0,0 +1,299 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilLocalResidual
+ */
+#ifndef EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH
+#define EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH
+
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilenergymodules.hh"
+#include "blackoilfoammodules.hh"
+
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief Calculates the local residual of the black oil model.
+ */
+template <class TypeTag>
+class BlackOilLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+
+    static const bool waterEnabled = Indices::waterEnabled;
+    static const bool gasEnabled = Indices::gasEnabled;
+    static const bool oilEnabled = Indices::oilEnabled;
+    static const bool compositionSwitchEnabled = (compositionSwitchIdx >= 0);
+
+    static constexpr bool blackoilConserveSurfaceVolume = GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume);
+    static constexpr bool enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy);
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef BlackOilSolventModule<TypeTag> SolventModule;
+    typedef BlackOilPolymerModule<TypeTag> PolymerModule;
+    typedef BlackOilEnergyModule<TypeTag> EnergyModule;
+    typedef BlackOilFoamModule<TypeTag> FoamModule;
+
+public:
+    /*!
+     * \copydoc FvBaseLocalResidual::computeStorage
+     */
+    template <class LhsEval>
+    void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                        const ElementContext& elemCtx,
+                        unsigned dofIdx,
+                        unsigned timeIdx) const
+    {
+        // retrieve the intensive quantities for the SCV at the specified point in time
+        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
+        const auto& fs = intQuants.fluidState();
+
+        storage = 0.0;
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                if (Indices::numPhases == 3) { // add trivial equation for the pseudo phase
+                    unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+                    if (timeIdx == 0)
+                        storage[conti0EqIdx + activeCompIdx] = Opm::variable<LhsEval>(0.0, conti0EqIdx + activeCompIdx);
+                    else
+                        storage[conti0EqIdx + activeCompIdx] = 0.0;
+                }
+                continue;
+            }
+
+            unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+            LhsEval surfaceVolume =
+                Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
+                * Toolbox::template decay<LhsEval>(fs.invB(phaseIdx))
+                * Toolbox::template decay<LhsEval>(intQuants.porosity());
+
+            storage[conti0EqIdx + activeCompIdx] += surfaceVolume;
+
+            // account for dissolved gas
+            if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
+                unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+                storage[conti0EqIdx + activeGasCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
+                    * surfaceVolume;
+            }
+
+            // account for vaporized oil
+            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
+                unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+                storage[conti0EqIdx + activeOilCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv())
+                    * surfaceVolume;
+            }
+        }
+
+        adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());
+
+        // deal with solvents (if present)
+        SolventModule::addStorage(storage, intQuants);
+
+        // deal with polymer (if present)
+        PolymerModule::addStorage(storage, intQuants);
+
+        // deal with energy (if present)
+        EnergyModule::addStorage(storage, intQuants);
+
+        // deal with foam (if present)
+        FoamModule::addStorage(storage, intQuants);
+    }
+
+    /*!
+     * \copydoc FvBaseLocalResidual::computeFlux
+     */
+    void computeFlux(RateVector& flux,
+                     const ElementContext& elemCtx,
+                     unsigned scvfIdx,
+                     unsigned timeIdx) const
+    {
+        assert(timeIdx == 0);
+
+        flux = 0.0;
+
+        const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned focusDofIdx = elemCtx.focusDofIndex();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
+            const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+            unsigned pvtRegionIdx = up.pvtRegionIndex();
+            if (upIdx == focusDofIdx)
+                evalPhaseFluxes_<Evaluation>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
+            else
+                evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
+        }
+
+        // deal with solvents (if present)
+        SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with polymer (if present)
+        PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with energy (if present)
+        EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with foam (if present)
+        FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+    }
+
+    /*!
+     * \copydoc FvBaseLocalResidual::computeSource
+     */
+    void computeSource(RateVector& source,
+                       const ElementContext& elemCtx,
+                       unsigned dofIdx,
+                       unsigned timeIdx) const
+    {
+        // retrieve the source term intrinsic to the problem
+        elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
+
+        // scale the source term of the energy equation
+        if (enableEnergy)
+            source[Indices::contiEnergyEqIdx] *= GET_PROP_VALUE(TypeTag, BlackOilEnergyScalingFactor);
+    }
+
+    /*!
+     * \brief Helper function to calculate the flux of mass in terms of conservation
+     *        quantities via specific fluid phase over a face.
+     */
+    template <class UpEval, class FluidState>
+    static void evalPhaseFluxes_(RateVector& flux,
+                                 unsigned phaseIdx,
+                                 unsigned pvtRegionIdx,
+                                 const ExtensiveQuantities& extQuants,
+                                 const FluidState& upFs)
+    {
+        const auto& invB = Opm::getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
+        const auto& surfaceVolumeFlux = invB*extQuants.volumeFlux(phaseIdx);
+        unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+
+        if (blackoilConserveSurfaceVolume)
+            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
+        else
+            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux*FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
+
+        if (phaseIdx == oilPhaseIdx) {
+            // dissolved gas (in the oil phase).
+            if (FluidSystem::enableDissolvedGas()) {
+                const auto& Rs = Opm::BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            }
+        }
+        else if (phaseIdx == gasPhaseIdx) {
+            // vaporized oil (in the gas phase).
+            if (FluidSystem::enableVaporizedOil()) {
+                const auto& Rv = Opm::BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux*FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            }
+        }
+    }
+
+    /*!
+     * \brief Helper function to convert the mass-related parts of a Dune::FieldVector
+     *        that stores conservation quantities in terms of "surface-volume" to the
+     *        conservation quantities used by the model.
+     *
+     * Depending on the value of the BlackoilConserveSurfaceVolume property, the model
+     * either conserves mass by means of "surface volume" of the components or mass
+     * directly. In the former case, this method is a no-op; in the latter, the values
+     * passed are multiplied by their respective pure component's density at surface
+     * conditions.
+     */
+    template <class Scalar>
+    static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container, unsigned pvtRegionIdx)
+    {
+        if (blackoilConserveSurfaceVolume)
+            return;
+
+        // convert "surface volume" to mass. this is complicated a bit by the fact that
+        // not all phases are necessarily enabled. (we here assume that if a fluid phase
+        // is disabled, its respective "main" component is not considered as well.)
+
+        if (waterEnabled) {
+            unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+            container[conti0EqIdx + activeWaterCompIdx] *=
+                FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+        }
+
+        if (gasEnabled) {
+            unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+            container[conti0EqIdx + activeGasCompIdx] *=
+                FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+        }
+
+        if (oilEnabled) {
+            unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+            container[conti0EqIdx + activeOilCompIdx] *=
+                FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+        }
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilmodel.hh

@@ -0,0 +1,566 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilModel
+ */
+#ifndef EWOMS_BLACK_OIL_MODEL_HH
+#define EWOMS_BLACK_OIL_MODEL_HH
+
+#include <opm/material/densead/Math.hpp>
+
+#include "blackoilproblem.hh"
+#include "blackoilindices.hh"
+#include "blackoiltwophaseindices.hh"
+#include "blackoilextensivequantities.hh"
+#include "blackoilprimaryvariables.hh"
+#include "blackoilintensivequantities.hh"
+#include "blackoilratevector.hh"
+#include "blackoilboundaryratevector.hh"
+#include "blackoillocalresidual.hh"
+#include "blackoilnewtonmethod.hh"
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilfoammodules.hh"
+#include "blackoildarcyfluxmodule.hh"
+
+#include <opm/models/common/multiphasebasemodel.hh>
+#include <ewoms/io/vtkcompositionmodule.hh>
+#include <ewoms/io/vtkblackoilmodule.hh>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <sstream>
+#include <string>
+
+namespace Opm {
+template <class TypeTag>
+class BlackOilModel;
+
+template <class TypeTag>
+class EclVanguard;
+}
+
+BEGIN_PROPERTIES
+
+//! The type tag for the black-oil problems
+NEW_TYPE_TAG(BlackOilModel, INHERITS_FROM(MultiPhaseBaseModel,
+                                          VtkBlackOil,
+                                          VtkBlackOilSolvent,
+                                          VtkBlackOilPolymer,
+                                          VtkBlackOilEnergy,
+                                          VtkComposition));
+
+//! Set the local residual function
+SET_TYPE_PROP(BlackOilModel, LocalResidual,
+              Opm::BlackOilLocalResidual<TypeTag>);
+
+//! Use the black-oil specific newton method
+SET_TYPE_PROP(BlackOilModel, NewtonMethod, Opm::BlackOilNewtonMethod<TypeTag>);
+
+//! The Model property
+SET_TYPE_PROP(BlackOilModel, Model, Opm::BlackOilModel<TypeTag>);
+
+//! The Problem property
+SET_TYPE_PROP(BlackOilModel, BaseProblem, Opm::BlackOilProblem<TypeTag>);
+
+//! the RateVector property
+SET_TYPE_PROP(BlackOilModel, RateVector, Opm::BlackOilRateVector<TypeTag>);
+
+//! the BoundaryRateVector property
+SET_TYPE_PROP(BlackOilModel, BoundaryRateVector, Opm::BlackOilBoundaryRateVector<TypeTag>);
+
+//! the PrimaryVariables property
+SET_TYPE_PROP(BlackOilModel, PrimaryVariables, Opm::BlackOilPrimaryVariables<TypeTag>);
+
+//! the IntensiveQuantities property
+SET_TYPE_PROP(BlackOilModel, IntensiveQuantities, Opm::BlackOilIntensiveQuantities<TypeTag>);
+
+//! the ExtensiveQuantities property
+SET_TYPE_PROP(BlackOilModel, ExtensiveQuantities, Opm::BlackOilExtensiveQuantities<TypeTag>);
+
+//! Use the the velocity module which is aware of the black-oil specific model extensions
+//! (i.e., the polymer and solvent extensions)
+SET_TYPE_PROP(BlackOilModel, FluxModule, Opm::BlackOilDarcyFluxModule<TypeTag>);
+
+//! The indices required by the model
+SET_TYPE_PROP(BlackOilModel, Indices,
+              Opm::BlackOilIndices<GET_PROP_VALUE(TypeTag, EnableSolvent),
+                                   GET_PROP_VALUE(TypeTag, EnablePolymer),
+                                   GET_PROP_VALUE(TypeTag, EnableEnergy),
+                                   GET_PROP_VALUE(TypeTag, EnableFoam),
+                                   /*PVOffset=*/0>);
+
+//! Set the fluid system to the black-oil fluid system by default
+SET_PROP(BlackOilModel, FluidSystem)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+
+public:
+    typedef Opm::BlackOilFluidSystem<Scalar> type;
+};
+
+// by default, all ECL extension modules are disabled
+SET_BOOL_PROP(BlackOilModel, EnableSolvent, false);
+SET_BOOL_PROP(BlackOilModel, EnablePolymer, false);
+SET_BOOL_PROP(BlackOilModel, EnablePolymerMW, false);
+SET_BOOL_PROP(BlackOilModel, EnableFoam, false);
+
+//! By default, the blackoil model is isothermal and does not conserve energy
+SET_BOOL_PROP(BlackOilModel, EnableTemperature, false);
+SET_BOOL_PROP(BlackOilModel, EnableEnergy, false);
+
+//! by default, scale the energy equation by the inverse of the energy required to heat
+//! up one kg of water by 30 Kelvin. If we conserve surface volumes, this must be divided
+//! by the weight of one cubic meter of water. This is required to make the "dumb" linear
+//! solvers that do not weight the components of the solutions do the right thing.
+//! by default, don't scale the energy equation, i.e. assume that a reasonable linear
+//! solver is used. (Not scaling it makes debugging quite a bit easier.)
+SET_PROP(BlackOilModel, BlackOilEnergyScalingFactor)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    static constexpr Scalar alpha = GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume) ? 1000.0 : 1.0;
+
+public:
+    typedef Scalar type;
+    static const Scalar value;
+};
+
+PROP_STATIC_CONST_MEMBER_DEFINITION_PREFIX_(BlackOilModel, BlackOilEnergyScalingFactor)
+    ::value = 1.0/(30*4184.0*alpha);
+
+// by default, ebos formulates the conservation equations in terms of mass not surface
+// volumes
+SET_BOOL_PROP(BlackOilModel, BlackoilConserveSurfaceVolume, false);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ * \brief A fully-implicit black-oil flow model.
+ *
+ * The black-oil model is a three-phase, three-component model widely
+ * used for oil reservoir simulation.  The phases are denoted by lower
+ * index \f$\alpha \in \{ w, g, o \}\f$ ("water", "gas" and "oil") and
+ * the components by upper index \f$\kappa \in \{ W, G, O \}\f$
+ * ("Water", "Gas" and "Oil"). The model assumes partial miscibility:
+ *
+ * - Water and the gas phases are immisicible and are assumed to be
+ *   only composed of the water and gas components respectively-
+ * - The oil phase is assumed to be a mixture of the gas and the oil
+ *  components.
+ *
+ * The densities of the phases are determined by so-called
+ * <i>formation volume factors</i>:
+ *
+ * \f[
+ * B_\alpha := \frac{\varrho_\alpha(1\,\text{bar})}{\varrho_\alpha(p_\alpha)}
+ * \f]
+ *
+ * Since the gas and water phases are assumed to be immiscible, this
+ * is sufficint to calculate their density. For the formation volume
+ * factor of the the oil phase \f$B_o\f$ determines the density of
+ * *saturated* oil, i.e. the density of the oil phase if some gas
+ * phase is present.
+ *
+ * The composition of the oil phase is given by the <i>gas dissolution factor</i>
+ * \f$R_s\f$, which defined as the volume of gas at atmospheric pressure that is
+ * dissolved in a given amount of oil at reservoir pressure:
+ *
+ * \f[
+ * R_s := \frac{\varrho_{o}^G}{\varrho_o^O}\;.
+ * \f]
+ *
+ * This allows to calculate all quantities required for the
+ * mass-conservation equations for each component, i.e.
+ *
+ * \f[
+ * \sum_\alpha \frac{\partial\;\phi c_\alpha^\kappa S_\alpha }{\partial t}
+ * - \sum_\alpha \mathrm{div} \left\{ c_\alpha^\kappa \mathbf{v}_\alpha \right\}
+ * - q^\kappa = 0 \;,
+ * \f]
+ * where \f$\mathrm{v}_\alpha\f$ is the filter velocity of the phase
+ * \f$\alpha\f$.
+ *
+ * By default \f$\mathrm{v}_\alpha\f$ is determined by using the
+ * standard multi-phase Darcy approach, i.e.
+ * \f[ \mathbf{v}_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
+ *\left(\mathbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right) \;, \f]
+ * although the actual approach which is used can be specified via the
+ * \c FluxModule property. For example, the velocity model can by
+ * changed to the Forchheimer approach by
+ * \code
+ * SET_TYPE_PROP(MyProblemTypeTag, FluxModule, Opm::ForchheimerFluxModule<TypeTag>);
+ * \endcode
+ *
+ * The primary variables used by this model are:
+ * - The pressure of the phase with the lowest index
+ * - The two saturations of the phases with the lowest indices
+ */
+template<class TypeTag >
+class BlackOilModel
+    : public MultiPhaseBaseModel<TypeTag>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
+    typedef MultiPhaseBaseModel<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = FluidSystem::numComponents };
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+    static const bool compositionSwitchEnabled = Indices::gasEnabled;
+    static const bool waterEnabled = Indices::waterEnabled;
+
+    typedef BlackOilSolventModule<TypeTag> SolventModule;
+    typedef BlackOilPolymerModule<TypeTag> PolymerModule;
+    typedef BlackOilEnergyModule<TypeTag> EnergyModule;
+public:
+    BlackOilModel(Simulator& simulator)
+        : ParentType(simulator)
+    {}
+
+    /*!
+     * \brief Register all run-time parameters for the immiscible model.
+     */
+    static void registerParameters()
+    {
+        ParentType::registerParameters();
+
+        SolventModule::registerParameters();
+        PolymerModule::registerParameters();
+        EnergyModule::registerParameters();
+
+        // register runtime parameters of the VTK output modules
+        Opm::VtkBlackOilModule<TypeTag>::registerParameters();
+        Opm::VtkCompositionModule<TypeTag>::registerParameters();
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::name
+     */
+    static std::string name()
+    { return "blackoil"; }
+
+    /*!
+     * \copydoc FvBaseDiscretization::primaryVarName
+     */
+    std::string primaryVarName(int pvIdx) const
+    {
+        std::ostringstream oss;
+
+        if (pvIdx == Indices::waterSaturationIdx)
+            oss << "saturation_" << FluidSystem::phaseName(FluidSystem::waterPhaseIdx);
+        else if (pvIdx == Indices::pressureSwitchIdx)
+            oss << "pressure_switching";
+        else if (static_cast<int>(pvIdx) == Indices::compositionSwitchIdx)
+            oss << "composition_switching";
+        else if (SolventModule::primaryVarApplies(pvIdx))
+            return SolventModule::primaryVarName(pvIdx);
+        else if (PolymerModule::primaryVarApplies(pvIdx))
+            return PolymerModule::primaryVarName(pvIdx);
+        else if (EnergyModule::primaryVarApplies(pvIdx))
+            return EnergyModule::primaryVarName(pvIdx);
+        else
+            assert(false);
+
+        return oss.str();
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::eqName
+     */
+    std::string eqName(int eqIdx) const
+    {
+        std::ostringstream oss;
+
+        if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx + numComponents)
+            oss << "conti_" << FluidSystem::phaseName(eqIdx - Indices::conti0EqIdx);
+        else if (SolventModule::eqApplies(eqIdx))
+            return SolventModule::eqName(eqIdx);
+        else if (PolymerModule::eqApplies(eqIdx))
+            return PolymerModule::eqName(eqIdx);
+        else if (EnergyModule::eqApplies(eqIdx))
+            return EnergyModule::eqName(eqIdx);
+        else
+            assert(false);
+
+        return oss.str();
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::primaryVarWeight
+     */
+    Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
+    {
+        // do not care about the auxiliary equations as they are supposed to scale
+        // themselves
+        if (globalDofIdx >= this->numGridDof())
+            return 1.0;
+
+        // saturations are always in the range [0, 1]!
+        if (int(Indices::waterSaturationIdx) == int(pvIdx))
+            return 1.0;
+
+        // oil pressures usually are in the range of 100 to 500 bars for typical oil
+        // reservoirs (which is the only relevant application for the black-oil model).
+        else if (int(Indices::pressureSwitchIdx) == int(pvIdx))
+            return 1.0/300e5;
+
+        // deal with primary variables stemming from the solvent module
+        else if (SolventModule::primaryVarApplies(pvIdx))
+            return SolventModule::primaryVarWeight(pvIdx);
+
+        // deal with primary variables stemming from the polymer module
+        else if (PolymerModule::primaryVarApplies(pvIdx))
+            return PolymerModule::primaryVarWeight(pvIdx);
+
+        // deal with primary variables stemming from the energy module
+        else if (EnergyModule::primaryVarApplies(pvIdx))
+            return EnergyModule::primaryVarWeight(pvIdx);
+
+        // if the primary variable is either the gas saturation, Rs or Rv
+        assert(int(Indices::compositionSwitchIdx) == int(pvIdx));
+
+        auto pvMeaning = this->solution(0)[globalDofIdx].primaryVarsMeaning();
+        if (pvMeaning == PrimaryVariables::Sw_po_Sg)
+            return 1.0; // gas saturation
+        else if (pvMeaning == PrimaryVariables::Sw_po_Rs)
+            return 1.0/250.; // gas dissolution factor
+        else {
+            assert(pvMeaning == PrimaryVariables::Sw_pg_Rv);
+            return 1.0/0.025; // oil vaporization factor
+        }
+
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::eqWeight
+     */
+    Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
+    {
+        // do not care about the auxiliary equations as they are supposed to scale
+        // themselves
+        if (globalDofIdx >= this->numGridDof())
+            return 1.0;
+
+        // we do not care much about water, so it gets de-prioritized by a factor of 100
+        static constexpr Scalar waterPriority = 1e-2;
+
+        if (GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume)) {
+            // Roughly convert the surface volume of the fluids from m^3 to kg. (in this
+            // context, it does not really matter if the actual densities are off by a
+            // factor of two or three.)
+            switch (eqIdx) {
+            case Indices::conti0EqIdx + FluidSystem::waterCompIdx:
+                return 1000.0*waterPriority;
+
+            case Indices::conti0EqIdx + FluidSystem::gasCompIdx:
+                return 1.0;
+
+            case Indices::conti0EqIdx + FluidSystem::oilCompIdx:
+                return 650.0;
+            }
+        }
+
+        if (SolventModule::eqApplies(eqIdx))
+            return SolventModule::eqWeight(eqIdx);
+
+        else if (PolymerModule::eqApplies(eqIdx))
+            return PolymerModule::eqWeight(eqIdx);
+
+        else if (EnergyModule::eqApplies(eqIdx))
+            return EnergyModule::eqWeight(eqIdx);
+
+        // it is said that all kilograms are born equal (except water)!
+        if (eqIdx == Indices::conti0EqIdx + FluidSystem::waterCompIdx)
+            return waterPriority;
+        return 1.0;
+    }
+
+    /*!
+     * \brief Write the current solution for a degree of freedom to a
+     *        restart file.
+     *
+     * \param outstream The stream into which the vertex data should
+     *                  be serialized to
+     * \param dof The Dune entity which's data should be serialized
+     */
+    template <class DofEntity>
+    void serializeEntity(std::ostream& outstream, const DofEntity& dof)
+    {
+        unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
+
+        // write phase state
+        if (!outstream.good())
+            throw std::runtime_error("Could not serialize degree of freedom "+std::to_string(dofIdx));
+
+        // write the primary variables
+        const auto& priVars = this->solution(/*timeIdx=*/0)[dofIdx];
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+            outstream << priVars[eqIdx] << " ";
+
+        // write the pseudo primary variables
+        outstream << priVars.primaryVarsMeaning() << " ";
+        outstream << priVars.pvtRegionIndex() << " ";
+
+        SolventModule::serializeEntity(*this, outstream, dof);
+        PolymerModule::serializeEntity(*this, outstream, dof);
+        EnergyModule::serializeEntity(*this, outstream, dof);
+    }
+
+    /*!
+     * \brief Reads the current solution variables for a degree of
+     *        freedom from a restart file.
+     *
+     * \param instream The stream from which the vertex data should
+     *                  be deserialized from
+     * \param dof The Dune entity which's data should be deserialized
+     */
+    template <class DofEntity>
+    void deserializeEntity(std::istream& instream,
+                           const DofEntity& dof)
+    {
+        unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
+
+        // read in the "real" primary variables of the DOF
+        auto& priVars = this->solution(/*timeIdx=*/0)[dofIdx];
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+            if (!instream.good())
+                throw std::runtime_error("Could not deserialize degree of freedom "+std::to_string(dofIdx));
+            instream >> priVars[eqIdx];
+        }
+
+        // read the pseudo primary variables
+        unsigned primaryVarsMeaning;
+        instream >> primaryVarsMeaning;
+
+        unsigned pvtRegionIdx;
+        instream >> pvtRegionIdx;
+
+        if (!instream.good())
+            throw std::runtime_error("Could not deserialize degree of freedom "+std::to_string(dofIdx));
+
+        SolventModule::deserializeEntity(*this, instream, dof);
+        PolymerModule::deserializeEntity(*this, instream, dof);
+        EnergyModule::deserializeEntity(*this, instream, dof);
+
+        typedef typename PrimaryVariables::PrimaryVarsMeaning PVM;
+        priVars.setPrimaryVarsMeaning(static_cast<PVM>(primaryVarsMeaning));
+        priVars.setPvtRegionIndex(pvtRegionIdx);
+    }
+
+    /*!
+     * \brief Deserializes the state of the model.
+     *
+     * \tparam Restarter The type of the serializer class
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void deserialize(Restarter& res)
+    {
+        ParentType::deserialize(res);
+
+        // set the PVT indices of the primary variables. This is also done by writing
+        // them into the restart file and re-reading them, but it is better to calculate
+        // them from scratch because the input could have been changed in this regard...
+        ElementContext elemCtx(this->simulator_);
+        auto elemIt = this->gridView().template begin</*codim=*/0>();
+        auto elemEndIt = this->gridView().template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++ elemIt) {
+            elemCtx.updateStencil(*elemIt);
+            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timIdx=*/0); ++dofIdx) {
+                unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timIdx=*/0);
+                updatePvtRegionIndex_(this->solution(/*timeIdx=*/0)[globalDofIdx],
+                                      elemCtx,
+                                      dofIdx,
+                                      /*timeIdx=*/0);
+            }
+        }
+
+        this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);
+    }
+
+/*
+    // hack: this interferes with the static polymorphism trick
+protected:
+    friend ParentType;
+    friend Discretization;
+*/
+
+    template <class Context>
+    void supplementInitialSolution_(PrimaryVariables& priVars,
+                                    const Context& context,
+                                    unsigned dofIdx,
+                                    unsigned timeIdx)
+    { updatePvtRegionIndex_(priVars, context, dofIdx, timeIdx); }
+
+    void registerOutputModules_()
+    {
+        ParentType::registerOutputModules_();
+
+        // add the VTK output modules which make sense for the blackoil model
+        SolventModule::registerOutputModules(*this, this->simulator_);
+        PolymerModule::registerOutputModules(*this, this->simulator_);
+        EnergyModule::registerOutputModules(*this, this->simulator_);
+
+        this->addOutputModule(new Opm::VtkBlackOilModule<TypeTag>(this->simulator_));
+        this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
+    }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+    template <class Context>
+    void updatePvtRegionIndex_(PrimaryVariables& priVars,
+                               const Context& context,
+                               unsigned dofIdx,
+                               unsigned timeIdx)
+    {
+        unsigned regionIdx = context.problem().pvtRegionIndex(context, dofIdx, timeIdx);
+        priVars.setPvtRegionIndex(regionIdx);
+    }
+};
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -0,0 +1,324 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilNewtonMethod
+ */
+#ifndef EWOMS_BLACK_OIL_NEWTON_METHOD_HH
+#define EWOMS_BLACK_OIL_NEWTON_METHOD_HH
+
+#include "blackoilproperties.hh"
+
+#include <opm/models/utils/signum.hh>
+
+#include <opm/material/common/Unused.hpp>
+
+BEGIN_PROPERTIES
+
+NEW_PROP_TAG(DpMaxRel);
+NEW_PROP_TAG(DsMax);
+NEW_PROP_TAG(PriVarOscilationThreshold);
+
+SET_SCALAR_PROP(NewtonMethod, DpMaxRel, 0.3);
+SET_SCALAR_PROP(NewtonMethod, DsMax, 0.2);
+SET_SCALAR_PROP(NewtonMethod, PriVarOscilationThreshold, 1e-5);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief A newton solver which is specific to the black oil model.
+ */
+template <class TypeTag>
+class BlackOilNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
+
+    static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
+
+public:
+    BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator)
+    {
+        priVarOscilationThreshold_ = EWOMS_GET_PARAM(TypeTag, Scalar, PriVarOscilationThreshold);
+        dpMaxRel_ = EWOMS_GET_PARAM(TypeTag, Scalar, DpMaxRel);
+        dsMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, DsMax);
+    }
+
+    /*!
+     * \copydoc NewtonMethod::finishInit()
+     */
+    void finishInit()
+    {
+        ParentType::finishInit();
+
+        wasSwitched_.resize(this->model().numTotalDof());
+        std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
+    }
+
+    /*!
+     * \brief Register all run-time parameters for the immiscible model.
+     */
+    static void registerParameters()
+    {
+        ParentType::registerParameters();
+
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, DpMaxRel, "Maximum relative change of pressure in a single iteration");
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, DsMax, "Maximum absolute change of any saturation in a single iteration");
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, PriVarOscilationThreshold,
+                             "The threshold value for the primary variable switching conditions after its meaning has switched to hinder oscilations");
+    }
+
+    /*!
+     * \brief Returns the number of degrees of freedom for which the
+     *        interpretation has changed for the most recent iteration.
+     */
+    unsigned numPriVarsSwitched() const
+    { return numPriVarsSwitched_; }
+
+protected:
+    friend NewtonMethod<TypeTag>;
+    friend ParentType;
+
+    /*!
+     * \copydoc FvBaseNewtonMethod::beginIteration_
+     */
+    void beginIteration_()
+    {
+        numPriVarsSwitched_ = 0;
+        ParentType::beginIteration_();
+    }
+
+    /*!
+     * \copydoc FvBaseNewtonMethod::endIteration_
+     */
+    void endIteration_(SolutionVector& uCurrentIter,
+                       const SolutionVector& uLastIter)
+    {
+#if HAVE_MPI
+        // in the MPI enabled case we need to add up the number of DOF
+        // for which the interpretation changed over all processes.
+        int localSwitched = numPriVarsSwitched_;
+        MPI_Allreduce(&localSwitched,
+                      &numPriVarsSwitched_,
+                      /*num=*/1,
+                      MPI_INT,
+                      MPI_SUM,
+                      MPI_COMM_WORLD);
+#endif // HAVE_MPI
+
+        this->simulator_.model().newtonMethod().endIterMsg()
+            << ", num switched=" << numPriVarsSwitched_;
+
+        ParentType::endIteration_(uCurrentIter, uLastIter);
+    }
+
+public:
+    void update_(SolutionVector& nextSolution,
+                 const SolutionVector& currentSolution,
+                 const GlobalEqVector& solutionUpdate,
+                 const GlobalEqVector& currentResidual)
+    {
+        const auto& comm = this->simulator_.gridView().comm();
+
+        int succeeded;
+        try {
+            ParentType::update_(nextSolution,
+                                currentSolution,
+                                solutionUpdate,
+                                currentResidual);
+            succeeded = 1;
+        }
+        catch (...) {
+            std::cout << "Newton update threw an exception on rank "
+                      << comm.rank() << "\n";
+            succeeded = 0;
+        }
+        succeeded = comm.min(succeeded);
+
+        if (!succeeded)
+            throw Opm::NumericalIssue("A process did not succeed in adapting the primary variables");
+
+        numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
+    }
+
+protected:
+    /*!
+     * \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
+     */
+    void updatePrimaryVariables_(unsigned globalDofIdx,
+                                 PrimaryVariables& nextValue,
+                                 const PrimaryVariables& currentValue,
+                                 const EqVector& update,
+                                 const EqVector& currentResidual)
+    {
+        static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
+        static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
+        static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
+        static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
+        static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
+
+        currentValue.checkDefined();
+        Opm::Valgrind::CheckDefined(update);
+        Opm::Valgrind::CheckDefined(currentResidual);
+
+        // saturation delta for each phase
+        Scalar deltaSw = 0.0;
+        Scalar deltaSo = 0.0;
+        Scalar deltaSg = 0.0;
+        Scalar deltaSs = 0.0;
+
+        if (Indices::waterEnabled) {
+            deltaSw = update[Indices::waterSaturationIdx];
+            deltaSo = -deltaSw;
+        }
+
+        if (Indices::gasEnabled && currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
+            deltaSg = update[Indices::compositionSwitchIdx];
+            deltaSo -= deltaSg;
+        }
+
+        if (enableSolvent) {
+            deltaSs = update[Indices::solventSaturationIdx];
+            deltaSo -= deltaSs;
+        }
+
+        // maximum saturation delta
+        Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
+        maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
+        maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
+
+        // scaling factor for saturation deltas to make sure that none of them exceeds
+        // the specified threshold value.
+        Scalar satAlpha = 1.0;
+        if (maxSatDelta > dsMax_)
+            satAlpha = dsMax_/maxSatDelta;
+
+        for (int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
+            // calculate the update of the current primary variable. For the black-oil
+            // model we limit the pressure delta relative to the pressure's current
+            // absolute value (Default: 30%) and saturation deltas to an absolute change
+            // (Default: 20%). Further, we ensure that the R factors, solvent
+            // "saturation" and polymer concentration do not become negative after the
+            // update.
+            Scalar delta = update[pvIdx];
+
+            // limit pressure delta
+            if (pvIdx == Indices::pressureSwitchIdx) {
+                if (std::abs(delta) > dpMaxRel_*currentValue[pvIdx])
+                    delta = Opm::signum(delta)*dpMaxRel_*currentValue[pvIdx];
+            }
+            // water saturation delta
+            else if (pvIdx == Indices::waterSaturationIdx)
+                delta *= satAlpha;
+            else if (pvIdx == Indices::compositionSwitchIdx) {
+                // the switching primary variable for composition is tricky because the
+                // "reasonable" value ranges it exhibits vary widely depending on its
+                // interpretation since it can represent Sg, Rs or Rv. For now, we only
+                // limit saturation deltas and ensure that the R factors do not become
+                // negative.
+                if (currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
+                    delta *= satAlpha;
+                else {
+                    if (delta > currentValue[Indices::compositionSwitchIdx])
+                        delta = currentValue[Indices::compositionSwitchIdx];
+                }
+            }
+            else if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
+                // solvent saturation updates are also subject to the Appleyard chop
+                delta *= satAlpha;
+            else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
+                const double sign = delta >= 0. ? 1. : -1.;
+                // maximum change of polymer molecular weight, the unit is MDa.
+                // applying this limit to stabilize the simulation. The value itself is still experimental.
+                const double maxMolarWeightChange = 100.0;
+                delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
+                delta *= satAlpha;
+            }
+
+            // do the actual update
+            nextValue[pvIdx] = currentValue[pvIdx] - delta;
+
+            // keep the solvent saturation between 0 and 1
+            if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
+                nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
+
+            // keep the polymer concentration above 0
+            if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx)
+                nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
+
+            if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
+                nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
+                const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
+                if (polymerConcentration < 1.e-10)
+                    nextValue[pvIdx] = 0.0;
+            }
+
+            // keep the foam concentration above 0
+            if (enableFoam && pvIdx == Indices::foamConcentrationIdx)
+                nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
+
+            // keep the temperature above 100 and below 1000 Kelvin
+            if (enableEnergy && pvIdx == Indices::temperatureIdx)
+                nextValue[pvIdx] = std::max(std::min(nextValue[pvIdx], 1000.0), 100.0);
+        }
+
+        // switch the new primary variables to something which is physically meaningful.
+        // use a threshold value after a switch to make it harder to switch back
+        // immediately.
+        if (wasSwitched_[globalDofIdx])
+            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, priVarOscilationThreshold_);
+        else
+            wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx);
+
+        if (wasSwitched_[globalDofIdx])
+            ++ numPriVarsSwitched_;
+
+        nextValue.checkDefined();
+    }
+
+private:
+    int numPriVarsSwitched_;
+
+    Scalar priVarOscilationThreshold_;
+    Scalar dpMaxRel_;
+    Scalar dsMax_;
+
+    // keep track of cells where the primary variable meaning has changed
+    // to detect and hinder oscillations
+    std::vector<bool> wasSwitched_;
+};
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -0,0 +1,647 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilPrimaryVariables
+ */
+#ifndef EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
+#define EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
+
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilenergymodules.hh"
+#include "blackoilfoammodules.hh"
+
+#include <ewoms/disc/common/fvbaseprimaryvariables.hh>
+
+#include <dune/common/fvector.hh>
+
+#include <opm/material/constraintsolvers/NcpFlash.hpp>
+#include <opm/material/fluidstates/CompositionalFluidState.hpp>
+#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+namespace Opm {
+template <class TypeTag, bool enableSolvent>
+class BlackOilSolventModule;
+
+template <class TypeTag, bool enablePolymer>
+class BlackOilPolymerModule;
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief Represents the primary variables used by the black-oil model.
+ */
+template <class TypeTag>
+class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
+{
+    typedef FvBasePrimaryVariables<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+
+    // number of equations
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+    // primary variable indices
+    enum { waterSaturationIdx = Indices::waterSaturationIdx };
+    enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
+    enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+
+    static const bool compositionSwitchEnabled = Indices::gasEnabled;
+    static const bool waterEnabled = Indices::waterEnabled;
+
+    // phase indices from the fluid system
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+
+    // component indices from the fluid system
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+
+    typedef typename Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+    typedef BlackOilSolventModule<TypeTag, enableSolvent> SolventModule;
+    typedef BlackOilPolymerModule<TypeTag, enablePolymer> PolymerModule;
+    typedef BlackOilEnergyModule<TypeTag, enableEnergy> EnergyModule;
+    typedef BlackOilFoamModule<TypeTag, enableFoam> FoamModule;
+
+    static_assert(numPhases == 3, "The black-oil model assumes three phases!");
+    static_assert(numComponents == 3, "The black-oil model assumes three components!");
+
+public:
+    enum PrimaryVarsMeaning {
+        Sw_po_Sg, // threephase case
+        Sw_po_Rs, // water + oil case
+        Sw_pg_Rv, // water + gas case
+    };
+
+    BlackOilPrimaryVariables()
+        : ParentType()
+    {
+        Opm::Valgrind::SetUndefined(*this);
+        pvtRegionIdx_ = 0;
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar)
+     */
+    BlackOilPrimaryVariables(Scalar value)
+        : ParentType(value)
+    {
+        Opm::Valgrind::SetUndefined(primaryVarsMeaning_);
+        pvtRegionIdx_ = 0;
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const ImmisciblePrimaryVariables& )
+     */
+    BlackOilPrimaryVariables(const BlackOilPrimaryVariables& value) = default;
+
+    /*!
+     * \brief Set the index of the region which should be used for PVT properties.
+     *
+     * The concept of PVT regions is a hack to work around the fact that the
+     * pseudo-components used by the black oil model (i.e., oil, gas and water) change
+     * their composition within the spatial domain. We implement them because, the ECL
+     * file format mandates them.
+     */
+    void setPvtRegionIndex(unsigned value)
+    { pvtRegionIdx_ = static_cast<unsigned short>(value); }
+
+    /*!
+     * \brief Return the index of the region which should be used for PVT properties.
+     */
+    unsigned pvtRegionIndex() const
+    { return pvtRegionIdx_; }
+
+    /*!
+     * \brief Return the interpretation which should be applied to the switching primary
+     *        variables.
+     */
+    PrimaryVarsMeaning primaryVarsMeaning() const
+    { return primaryVarsMeaning_; }
+
+    /*!
+     * \brief Set the interpretation which should be applied to the switching primary
+     *        variables.
+     */
+    void setPrimaryVarsMeaning(PrimaryVarsMeaning newMeaning)
+    { primaryVarsMeaning_ = newMeaning; }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::assignMassConservative
+     */
+    template <class FluidState>
+    void assignMassConservative(const FluidState& fluidState,
+                                const MaterialLawParams& matParams,
+                                bool isInEquilibrium = false)
+    {
+        typedef typename std::remove_reference<typename FluidState::Scalar>::type ConstEvaluation;
+        typedef typename std::remove_const<ConstEvaluation>::type FsEvaluation;
+        typedef typename Opm::MathToolbox<FsEvaluation> FsToolbox;
+
+#ifndef NDEBUG
+        // make sure the temperature is the same in all fluid phases
+        for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
+            Opm::Valgrind::CheckDefined(fluidState.temperature(0));
+            Opm::Valgrind::CheckDefined(fluidState.temperature(phaseIdx));
+
+            assert(fluidState.temperature(0) == fluidState.temperature(phaseIdx));
+        }
+#endif // NDEBUG
+
+        // for the equilibrium case, we don't need complicated
+        // computations.
+        if (isInEquilibrium) {
+            assignNaive(fluidState);
+            return;
+        }
+
+        // If your compiler bails out here, you're probably not using a suitable black
+        // oil fluid system.
+        typename FluidSystem::template ParameterCache<Scalar> paramCache;
+        paramCache.setRegionIndex(pvtRegionIdx_);
+        paramCache.setMaxOilSat(FsToolbox::value(fluidState.saturation(oilPhaseIdx)));
+
+        // create a mutable fluid state with well defined densities based on the input
+        typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
+        typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FlashFluidState;
+        FlashFluidState fsFlash;
+        fsFlash.setTemperature(FsToolbox::value(fluidState.temperature(/*phaseIdx=*/0)));
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            fsFlash.setPressure(phaseIdx, FsToolbox::value(fluidState.pressure(phaseIdx)));
+            fsFlash.setSaturation(phaseIdx, FsToolbox::value(fluidState.saturation(phaseIdx)));
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
+                fsFlash.setMoleFraction(phaseIdx, compIdx, FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx)));
+        }
+
+        paramCache.updateAll(fsFlash);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            Scalar rho = FluidSystem::template density<FlashFluidState, Scalar>(fsFlash, paramCache, phaseIdx);
+            fsFlash.setDensity(phaseIdx, rho);
+        }
+
+        // calculate the "global molarities"
+        ComponentVector globalMolarities(0.0);
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                if (!FluidSystem::phaseIsActive(phaseIdx))
+                    continue;
+
+                globalMolarities[compIdx] +=
+                    fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
+            }
+        }
+
+        // use a flash calculation to calculate a fluid state in
+        // thermodynamic equilibrium
+
+        // run the flash calculation
+        NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
+
+        // use the result to assign the primary variables
+        assignNaive(fsFlash);
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::assignNaive
+     */
+    template <class FluidState>
+    void assignNaive(const FluidState& fluidState)
+    {
+        typedef typename std::remove_reference<typename FluidState::Scalar>::type ConstEvaluation;
+        typedef typename std::remove_const<ConstEvaluation>::type FsEvaluation;
+        typedef typename Opm::MathToolbox<FsEvaluation> FsToolbox;
+
+        bool gasPresent = FluidSystem::phaseIsActive(gasPhaseIdx)?(fluidState.saturation(gasPhaseIdx) > 0.0):false;
+        bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 0.0):false;
+        static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
+        bool onlyWater = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > thresholdWaterFilledCell):false;
+
+        // deal with the primary variables for the energy extension
+        EnergyModule::assignPrimaryVars(*this, fluidState);
+
+        // determine the meaning of the primary variables
+        if ((gasPresent && oilPresent) || onlyWater)
+            // gas and oil: both hydrocarbon phases are in equilibrium (i.e., saturated
+            // with the "protagonist" component of the other phase.)
+            primaryVarsMeaning_ = Sw_po_Sg;
+        else if (oilPresent) {
+            // only oil: if dissolved gas is enabled, we need to consider the oil phase
+            // composition, if it is disabled, the gas component must stick to its phase
+            if (FluidSystem::enableDissolvedGas())
+                primaryVarsMeaning_ = Sw_po_Rs;
+            else
+                primaryVarsMeaning_ = Sw_po_Sg;
+        }
+        else {
+            assert(gasPresent);
+            // only gas: if vaporized oil is enabled, we need to consider the gas phase
+            // composition, if it is disabled, the oil component must stick to its phase
+            if (FluidSystem::enableVaporizedOil())
+                primaryVarsMeaning_ = Sw_pg_Rv;
+            else
+                primaryVarsMeaning_ = Sw_po_Sg;
+        }
+
+        // assign the actual primary variables
+        if (primaryVarsMeaning() == Sw_po_Sg) {
+            if (waterEnabled)
+                (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+            (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
+            if( compositionSwitchEnabled )
+                (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
+        }
+        else if (primaryVarsMeaning() == Sw_po_Rs) {
+            const auto& Rs = Opm::BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
+
+            if (waterEnabled)
+                (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+            (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
+            if( compositionSwitchEnabled )
+                (*this)[compositionSwitchIdx] = Rs;
+        }
+        else {
+            assert(primaryVarsMeaning() == Sw_pg_Rv);
+
+            const auto& Rv = Opm::BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
+            if (waterEnabled)
+                (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+
+            (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
+            if( compositionSwitchEnabled )
+                (*this)[compositionSwitchIdx] = Rv;
+        }
+
+        checkDefined();
+    }
+
+    /*!
+     * \brief Adapt the interpretation of the switching variables to be physically
+     *        meaningful.
+     *
+     * If the meaning of the primary variables changes, their values are also adapted in a
+     * meaningful manner. (e.g. if the gas phase appears and the composition switching
+     * variable changes its meaning from the gas dissolution factor Rs to the gas
+     * saturation Sg, the value for this variable is set to zero.)
+     * A Scalar eps can be passed to make the switching condition more strict.
+     * Useful for avoiding ocsilation in the primaryVarsMeaning.
+     *
+     * \return true Iff the interpretation of one of the switching variables was changed
+     */
+    bool adaptPrimaryVariables(const Problem& problem, unsigned globalDofIdx, Scalar eps = 0.0)
+    {
+        static const Scalar thresholdWaterFilledCell = 1.0 - eps;
+
+        // this function accesses quite a few black-oil specific low-level functions
+        // directly for better performance (instead of going the canonical way through
+        // the IntensiveQuantities). The reason is that most intensive quantities are not
+        // required to be able to decide if the primary variables needs to be switched or
+        // not, so it would be a waste to compute them.
+        Scalar Sw = 0.0;
+        if (waterEnabled)
+            Sw = (*this)[Indices::waterSaturationIdx];
+
+        if (primaryVarsMeaning() == Sw_po_Sg) {
+
+            // special case for cells with almost only water
+            if (Sw >= thresholdWaterFilledCell) {
+
+                // make sure water saturations does not exceed 1.0
+                if (waterEnabled)
+                    (*this)[Indices::waterSaturationIdx] = 1.0;
+                // the hydrocarbon gas saturation is set to 0.0
+                if (compositionSwitchEnabled)
+                    (*this)[Indices::compositionSwitchIdx] = 0.0;
+
+                return false;
+            }
+
+            // phase equilibrium, i.e., both hydrocarbon phases are present.
+            Scalar Sg = 0.0;
+            if (compositionSwitchEnabled)
+                Sg = (*this)[Indices::compositionSwitchIdx];
+
+            Scalar So = 1.0 - Sw - Sg - solventSaturation_();
+
+            Scalar So2 = 1.0 - Sw - solventSaturation_();
+            if (Sg < -eps && So2 > 0.0 && FluidSystem::enableDissolvedGas()) {
+                // the hydrocarbon gas phase disappeared and some oil phase is left,
+                // i.e., switch the primary variables to { Sw, po, Rs }.
+                //
+                // by a "lucky" coincidence the pressure switching variable already
+                // represents the oil phase pressure, so we do not need to change
+                // this. For the gas dissolution factor, we use the low-level blackoil
+                // PVT objects to calculate the mole fraction of gas saturated oil.
+                Scalar po = (*this)[Indices::pressureSwitchIdx];
+                Scalar T = asImp_().temperature_();
+                Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+                Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
+                Scalar RsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
+                                                                                   T,
+                                                                                   po,
+                                                                                   So2,
+                                                                                   SoMax);
+
+                setPrimaryVarsMeaning(Sw_po_Rs);
+                if (compositionSwitchEnabled)
+                    (*this)[Indices::compositionSwitchIdx] =
+                        std::min(RsMax, RsSat);
+                return true;
+            }
+
+            Scalar Sg2 = 1.0 - Sw - solventSaturation_();
+            if (So < -eps && Sg2 > 0.0 && FluidSystem::enableVaporizedOil()) {
+                // the oil phase disappeared and some hydrocarbon gas phase is still
+                // present, i.e., switch the primary variables to { Sw, pg, Rv }.
+                Scalar po = (*this)[Indices::pressureSwitchIdx];
+
+                // we only have the oil pressure readily available, but we need the gas
+                // pressure, i.e. we must determine capillary pressure
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC, /*So=*/0.0, Sg2 + solventSaturation_(), Sw, matParams);
+                Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+
+                // we start at the Rv value that corresponds to that of oil-saturated
+                // hydrocarbon gas
+                Scalar T = asImp_().temperature_();
+                Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+                Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
+                Scalar RvSat =
+                    FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                                         T,
+                                                                         pg,
+                                                                         Scalar(0),
+                                                                         SoMax);
+                setPrimaryVarsMeaning(Sw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                if (compositionSwitchEnabled)
+                    (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+
+                return true;
+            }
+
+            return false;
+        }
+        else if (primaryVarsMeaning() == Sw_po_Rs) {
+            assert(compositionSwitchEnabled);
+
+            // special case for cells with almost only water
+            if (Sw >= thresholdWaterFilledCell) {
+                // switch back to phase equilibrium mode if the oil phase vanishes (i.e.,
+                // the water-only case)
+                setPrimaryVarsMeaning(Sw_po_Sg);
+                if (waterEnabled)
+                    (*this)[Indices::waterSaturationIdx] = 1.0; // water saturation
+
+                (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
+
+                return true;
+            }
+
+            // Only the oil and the water phases are present. The hydrocarbon gas phase
+            // appears as soon as more of the gas component is present in the oil phase
+            // than what saturated oil can hold.
+            Scalar T = asImp_().temperature_();
+            Scalar po = (*this)[Indices::pressureSwitchIdx];
+            Scalar So = 1.0 - Sw - solventSaturation_();
+            Scalar SoMax = std::max(So, problem.maxOilSaturation(globalDofIdx));
+            Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
+            Scalar RsSat =
+                FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
+                                                                    T,
+                                                                    po,
+                                                                    So,
+                                                                    SoMax);
+
+            Scalar Rs = (*this)[Indices::compositionSwitchIdx];
+            if (Rs > std::min(RsMax, RsSat*(1.0 + eps))) {
+                // the gas phase appears, i.e., switch the primary variables to { Sw, po,
+                // Sg }.
+                setPrimaryVarsMeaning(Sw_po_Sg);
+                (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
+
+                return true;
+            }
+
+            return false;
+        }
+        else {
+            assert(primaryVarsMeaning() == Sw_pg_Rv);
+            assert(compositionSwitchEnabled);
+
+            Scalar pg = (*this)[Indices::pressureSwitchIdx];
+            Scalar Sg = 1.0 - Sw - solventSaturation_();
+
+            // special case for cells with almost only water
+            if (Sw >= thresholdWaterFilledCell) {
+                // switch to phase equilibrium mode because the hydrocarbon gas phase
+                // disappears. here we need the capillary pressures to calculate the oil
+                // phase pressure using the gas phase pressure
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC,
+                                           /*So=*/0.0,
+                                           /*Sg=*/Sg + solventSaturation_(),
+                                           Sw,
+                                           matParams);
+                Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
+
+                setPrimaryVarsMeaning(Sw_po_Sg);
+                if (waterEnabled)
+                    (*this)[Indices::waterSaturationIdx] = 1.0;
+
+                (*this)[Indices::pressureSwitchIdx] = po;
+                (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
+
+                return true;
+            }
+
+            // Only the gas and the water phases are present. The oil phase appears as
+            // soon as more of the oil component is present in the hydrocarbon gas phase
+            // than what saturated gas contains. Note that we use the blackoil specific
+            // low-level PVT objects here for performance reasons.
+            Scalar T = asImp_().temperature_();
+            Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+            Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
+            Scalar RvSat =
+                FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                                     T,
+                                                                     pg,
+                                                                     /*So=*/Scalar(0.0),
+                                                                     SoMax);
+
+            Scalar Rv = (*this)[Indices::compositionSwitchIdx];
+            if (Rv > std::min(RvMax, RvSat*(1.0 + eps))) {
+                // switch to phase equilibrium mode because the oil phase appears. here
+                // we also need the capillary pressures to calculate the oil phase
+                // pressure using the gas phase pressure
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC,
+                                           /*So=*/0.0,
+                                           /*Sg=*/Sg + solventSaturation_(),
+                                           Sw,
+                                           matParams);
+                Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
+
+                setPrimaryVarsMeaning(Sw_po_Sg);
+                (*this)[Indices::pressureSwitchIdx] = po;
+                (*this)[Indices::compositionSwitchIdx] = Sg; // hydrocarbon gas saturation
+
+                return true;
+            }
+
+            return false;
+        }
+
+        assert(false);
+        return false;
+    }
+
+    BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;
+    BlackOilPrimaryVariables& operator=(Scalar value)
+    {
+        for (unsigned i = 0; i < numEq; ++i)
+            (*this)[i] = value;
+
+        return *this;
+    }
+
+    /*!
+     * \brief Instruct valgrind to check the definedness of all attributes of this class.
+     *
+     * We cannot simply check the definedness of the whole object because there might be
+     * "alignedness holes" in the memory layout which are caused by the pseudo primary
+     * variables.
+     */
+    void checkDefined() const
+    {
+#ifndef NDEBUG
+        // check the "real" primary variables
+        for (unsigned i = 0; i < this->size(); ++i)
+            Opm::Valgrind::CheckDefined((*this)[i]);
+
+        // check the "pseudo" primary variables
+        Opm::Valgrind::CheckDefined(primaryVarsMeaning_);
+        Opm::Valgrind::CheckDefined(pvtRegionIdx_);
+#endif // NDEBUG
+    }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+    Scalar solventSaturation_() const
+    {
+        if (!enableSolvent)
+            return 0.0;
+
+        return (*this)[Indices::solventSaturationIdx];
+    }
+
+    Scalar polymerConcentration_() const
+    {
+        if (!enablePolymer)
+            return 0.0;
+
+        return (*this)[Indices::polymerConcentrationIdx];
+    }
+
+    Scalar foamConcentration_() const
+    {
+        if (!enableFoam)
+            return 0.0;
+
+        return (*this)[Indices::foamConcentrationIdx];
+    }
+
+    Scalar temperature_() const
+    {
+        if (!enableEnergy)
+            return FluidSystem::reservoirTemperature();
+
+        return (*this)[Indices::temperatureIdx];
+    }
+
+    template <class Container>
+    void computeCapillaryPressures_(Container& result,
+                                    Scalar So,
+                                    Scalar Sg,
+                                    Scalar Sw,
+                                    const MaterialLawParams& matParams) const
+    {
+        typedef Opm::SimpleModularFluidState<Scalar,
+                                             numPhases,
+                                             numComponents,
+                                             FluidSystem,
+                                             /*storePressure=*/false,
+                                             /*storeTemperature=*/false,
+                                             /*storeComposition=*/false,
+                                             /*storeFugacity=*/false,
+                                             /*storeSaturation=*/true,
+                                             /*storeDensity=*/false,
+                                             /*storeViscosity=*/false,
+                                             /*storeEnthalpy=*/false> SatOnlyFluidState;
+        SatOnlyFluidState fluidState;
+        fluidState.setSaturation(waterPhaseIdx, Sw);
+        fluidState.setSaturation(oilPhaseIdx, So);
+        fluidState.setSaturation(gasPhaseIdx, Sg);
+
+        MaterialLaw::capillaryPressures(result, matParams, fluidState);
+    }
+
+    PrimaryVarsMeaning primaryVarsMeaning_;
+    unsigned short pvtRegionIdx_;
+};
+
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilproblem.hh

@@ -0,0 +1,176 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilProblem
+ */
+#ifndef EWOMS_BLACKOIL_PROBLEM_HH
+#define EWOMS_BLACKOIL_PROBLEM_HH
+
+#include "blackoilproperties.hh"
+
+#include <opm/models/common/multiphasebaseproblem.hh>
+
+#include <opm/material/common/Unused.hpp>
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ * \brief Base class for all problems which use the black-oil model.
+ */
+template<class TypeTag>
+class BlackOilProblem : public MultiPhaseBaseProblem<TypeTag>
+{
+private:
+    typedef MultiPhaseBaseProblem<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+
+public:
+    /*!
+     * \copydoc Doxygen::defaultProblemConstructor
+     *
+     * \param simulator The manager object of the simulation
+     */
+    BlackOilProblem(Simulator& simulator)
+        : ParentType(simulator)
+    {}
+
+    /*!
+     * \brief Returns the maximum value of the gas dissolution factor at the current time
+     *        for a given degree of freedom.
+     *
+     * This is required for the DRSDT keyword.
+     */
+    Scalar maxGasDissolutionFactor(unsigned timeIdx OPM_UNUSED, unsigned globalDofIdx OPM_UNUSED) const
+    { return std::numeric_limits<Scalar>::max()/2; }
+
+    /*!
+     * \brief Returns the maximum value of the oil vaporization factor at the current
+     *        time for a given degree of freedom.
+     *
+     * This is required for the DRVDT keyword.
+     */
+    Scalar maxOilVaporizationFactor(unsigned timeIdx OPM_UNUSED, unsigned globalDofIdx OPM_UNUSED) const
+    { return std::numeric_limits<Scalar>::max()/2; }
+
+    /*!
+     * \brief Returns the maximum value of the oil saturation seen at the current time
+     *        for a given degree of freedom.
+     *
+     * This is required for the VAPPARS keyword.
+     */
+    Scalar maxOilSaturation(unsigned globalDofIdx OPM_UNUSED) const
+    { return 1.0; }
+
+    /*!
+     * \brief Returns the index of the relevant region for thermodynmic properties
+     */
+    template <class Context>
+    unsigned pvtRegionIndex(const Context& context OPM_UNUSED,
+                            unsigned spaceIdx OPM_UNUSED,
+                            unsigned timeIdx OPM_UNUSED) const
+    { return 0; }
+
+    /*!
+     * \brief Returns the index of the relevant region for saturation functions
+     */
+    template <class Context>
+    unsigned satnumRegionIndex(const Context& context OPM_UNUSED,
+                               unsigned spaceIdx OPM_UNUSED,
+                               unsigned timeIdx OPM_UNUSED) const
+    { return 0; }
+
+    /*!
+     * \brief Returns the index of the relevant region for solvent mixing functions
+     */
+    template <class Context>
+    unsigned miscnumRegionIndex(const Context& context OPM_UNUSED,
+                                unsigned spaceIdx OPM_UNUSED,
+                                unsigned timeIdx OPM_UNUSED) const
+    { return 0; }
+
+    /*!
+     * \brief Returns the index of the relevant region for polymer mixing functions
+     */
+    template <class Context>
+    unsigned plmixnumRegionIndex(const Context& context OPM_UNUSED,
+                                 unsigned spaceIdx OPM_UNUSED,
+                                 unsigned timeIdx OPM_UNUSED) const
+    { return 0; }
+
+    /*!
+     * \brief Returns the compressibility of the porous medium of a cell
+     */
+    template <class Context>
+    Scalar rockCompressibility(const Context& context OPM_UNUSED,
+                               unsigned spaceIdx OPM_UNUSED,
+                               unsigned timeIdx OPM_UNUSED) const
+    { return 0.0; }
+
+    /*!
+     * \brief Returns the reference pressure for rock the compressibility of a cell
+     */
+    template <class Context>
+    Scalar rockReferencePressure(const Context& context OPM_UNUSED,
+                                 unsigned spaceIdx OPM_UNUSED,
+                                 unsigned timeIdx OPM_UNUSED) const
+    { return 1e5; }
+
+    /*!
+     * \brief Returns the reference temperature
+     *
+     * This is only relevant for temperature dependent quantities, in particular those
+     * needed by the module for energy conservation.
+     */
+    Scalar referenceTemperature() const
+    { return 273.15 + 15.56; /* [K] */ }
+
+    /*!
+     * \brief Returns the porosity multiplier due to water-induced rock compaction
+     *
+     * This is a somewhat exotic feature. Most likely you will not need to touch this
+     * method.
+     */
+    template <class Evaluation>
+    Scalar rockCompPoroMultiplier(const IntensiveQuantities& intQuants OPM_UNUSED,
+                                  unsigned globalSpaceIdx OPM_UNUSED) const
+    { return 1.0; }
+
+private:
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation& asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilproperties.hh

@@ -0,0 +1,79 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ * \ingroup BlackOilModel
+ *
+ * \brief Declares the properties required by the black oil model.
+ */
+#ifndef EWOMS_BLACK_OIL_PROPERTIES_HH
+#define EWOMS_BLACK_OIL_PROPERTIES_HH
+
+#include <opm/models/common/multiphasebaseproperties.hh>
+
+BEGIN_PROPERTIES
+
+//! Specifies if the simulation should write output files that are
+//! compatible with those produced by the commercial Eclipse simulator
+NEW_PROP_TAG(EnableEclipseOutput);
+//! The material law for thermal conduction
+NEW_PROP_TAG(ThermalConductionLaw);
+//! The parameters of the material law for thermal conduction
+NEW_PROP_TAG(ThermalConductionLawParams);
+//! The material law for energy storage of the rock
+NEW_PROP_TAG(SolidEnergyLaw);
+//! The parameters for material law for energy storage of the rock
+NEW_PROP_TAG(SolidEnergyLawParams);
+//! Enable the ECL-blackoil extension for solvents. ("Second gas")
+NEW_PROP_TAG(EnableSolvent);
+//! Enable the ECL-blackoil extension for polymer.
+NEW_PROP_TAG(EnablePolymer);
+//! Enable the tracking polymer molecular weight tracking and related functionalities
+NEW_PROP_TAG(EnablePolymerMW);
+//! Enable surface volume scaling
+NEW_PROP_TAG(BlackoilConserveSurfaceVolume);
+//! Enable the ECL-blackoil extension for foam
+NEW_PROP_TAG(EnableFoam);
+
+//! Allow the spatial and temporal domains to exhibit non-constant temperature
+//! in the black-oil model
+NEW_PROP_TAG(EnableTemperature);
+
+//! Enable the ECL-blackoil extension for energy conservation
+//!
+//! Setting this property to true implies EnableTemperature.
+NEW_PROP_TAG(EnableEnergy);
+
+//! The relative weight of the residual of the energy equation compared to the mass
+//! residuals
+//!
+//! this is basically a hack to work around the limitation that the convergence criterion
+//! of unmodified dune-istl linear solvers cannot weight the individual equations. if the
+//! energy equation is not scaled, its absolute value is normally several orders of
+//! magnitude larger than that of the mass balance equations
+NEW_PROP_TAG(BlackOilEnergyScalingFactor);
+
+
+END_PROPERTIES
+
+#endif

opm/models/blackoil/blackoilratevector.hh

@@ -0,0 +1,213 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilRateVector
+ */
+#ifndef EWOMS_BLACK_OIL_RATE_VECTOR_HH
+#define EWOMS_BLACK_OIL_RATE_VECTOR_HH
+
+#include <dune/common/fvector.hh>
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/constraintsolvers/NcpFlash.hpp>
+
+#include "blackoilintensivequantities.hh"
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief Implements a vector representing mass, molar or volumetric rates for
+ *        the black oil model.
+ *
+ * This class is basically a Dune::FieldVector which can be set using
+ * either mass, molar or volumetric rates.
+ */
+template <class TypeTag>
+class BlackOilRateVector
+    : public Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Evaluation),
+                               GET_PROP_VALUE(TypeTag, NumEq)>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    typedef BlackOilSolventModule<TypeTag> SolventModule;
+    typedef BlackOilPolymerModule<TypeTag> PolymerModule;
+    typedef BlackOilFoamModule<TypeTag> FoamModule;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { contiEnergyEqIdx = Indices::contiEnergyEqIdx };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enablePolymerMolarWeight = GET_PROP_VALUE(TypeTag, EnablePolymerMW) };
+    enum { enableFoam = GET_PROP_VALUE(TypeTag, EnableFoam) };
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldVector<Evaluation, numEq> ParentType;
+
+public:
+    BlackOilRateVector() : ParentType()
+    { Opm::Valgrind::SetUndefined(*this); }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar)
+     */
+    BlackOilRateVector(Scalar value) : ParentType(Toolbox::createConstant(value))
+    {}
+
+    template <class Eval = Evaluation>
+    BlackOilRateVector(const typename std::enable_if<std::is_same<Eval, Evaluation>::value, Evaluation>::type& value) : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const
+     * ImmiscibleRateVector& )
+     */
+    BlackOilRateVector(const BlackOilRateVector& value) : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setMassRate
+     */
+    void setMassRate(const ParentType& value, unsigned pvtRegionIdx = 0)
+    {
+        ParentType::operator=(value);
+
+        // convert to "surface volume" if requested
+        if (GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume)) {
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                (*this)[FluidSystem::gasCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx);
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                (*this)[FluidSystem::oilCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx);
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                (*this)[FluidSystem::waterCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx);
+            }
+            if (enableSolvent) {
+                const auto& solventPvt = SolventModule::solventPvt();
+                (*this)[Indices::contiSolventEqIdx] /=
+                        solventPvt.referenceDensity(pvtRegionIdx);
+            }
+
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setMolarRate
+     */
+    void setMolarRate(const ParentType& value, unsigned pvtRegionIdx = 0)
+    {
+        // first, assign molar rates
+        ParentType::operator=(value);
+
+        // then, convert them to mass rates
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
+            (*this)[conti0EqIdx + compIdx] *= FluidSystem::molarMass(compIdx, pvtRegionIdx);
+
+        const auto& solventPvt = SolventModule::solventPvt();
+        (*this)[Indices::contiSolventEqIdx] *= solventPvt.molarMass(pvtRegionIdx);
+
+        if ( enablePolymer ) {
+            if (enablePolymerMolarWeight )
+                throw std::logic_error("Set molar rate with polymer weight tracking not implemented");
+
+            (*this)[Indices::contiPolymerEqIdx] *= PolymerModule::molarMass(pvtRegionIdx);
+        }
+
+        if ( enableFoam ) {
+            throw std::logic_error("setMolarRate() not implemented for foam");
+        }
+
+        // convert to "surface volume" if requested
+        if (GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume)) {
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                (*this)[FluidSystem::gasCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx);
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                (*this)[FluidSystem::oilCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx);
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                (*this)[FluidSystem::waterCompIdx] /=
+                        FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx);
+            }
+            if (enableSolvent) {
+                (*this)[Indices::contiSolventEqIdx] /=
+                        solventPvt.referenceDensity(pvtRegionIdx);
+            }
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setVolumetricRate
+     */
+    template <class FluidState, class RhsEval>
+    void setVolumetricRate(const FluidState& fluidState,
+                           unsigned phaseIdx,
+                           const RhsEval& volume)
+    {
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
+            (*this)[conti0EqIdx + compIdx] =
+                fluidState.density(phaseIdx)
+                * fluidState.massFraction(phaseIdx, compIdx)
+                * volume;
+    }
+
+    /*!
+     * \brief Assignment operator from a scalar or a function evaluation
+     */
+    template <class RhsEval>
+    BlackOilRateVector& operator=(const RhsEval& value)
+    {
+        for (unsigned i=0; i < this->size(); ++i)
+            (*this)[i] = value;
+        return *this;
+    }
+
+    /*!
+     * \brief Assignment operator from another rate vector
+     */
+    BlackOilRateVector& operator=(const BlackOilRateVector& other)
+    {
+        for (unsigned i=0; i < this->size(); ++i)
+            (*this)[i] = other[i];
+        return *this;
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -0,0 +1,181 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilTwoPhaseIndices
+ */
+#ifndef EWOMS_BLACK_OIL_TWO_PHASE_INDICES_HH
+#define EWOMS_BLACK_OIL_TWO_PHASE_INDICES_HH
+
+#include <cassert>
+
+namespace Opm {
+
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief The primary variable and equation indices for the black-oil model.
+ */
+template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
+struct BlackOilTwoPhaseIndices
+{
+    //! Is phase enabled or not
+    static const bool oilEnabled = (disabledCanonicalCompIdx != 0);
+    static const bool waterEnabled = (disabledCanonicalCompIdx != 1);
+    static const bool gasEnabled = (disabledCanonicalCompIdx != 2);
+
+    //! Are solvents involved?
+    static const bool enableSolvent = numSolventsV > 0;
+
+    //! Are polymers involved?
+    static const bool enablePolymer = numPolymersV > 0;
+
+    //! Shall energy be conserved?
+    static const bool enableEnergy = numEnergyV > 0;
+
+    //! Number of solvent components to be considered
+    static const int numSolvents = enableSolvent ? numSolventsV : 0;
+
+    //! Number of polymer components to be considered
+    static const int numPolymers = enablePolymer ? numPolymersV : 0;
+
+    //! Number of energy equations to be considered
+    static const int numEnergy = enableEnergy ? numEnergyV : 0;
+
+    //! Number of foam equations to be considered
+    static const int numFoam = enableFoam? 1 : 0;
+
+    //! The number of fluid phases
+    static const int numPhases = 2;
+
+    //! The number of equations
+    static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam;
+
+    //////////////////////////////
+    // Primary variable indices
+    //////////////////////////////
+
+    //! The index of the water saturation. For two-phase oil gas models this is disabled.
+    static const int waterSaturationIdx  = waterEnabled ? PVOffset + 0 : -10000;
+
+    //! Index of the oil pressure in a vector of primary variables
+    static const int pressureSwitchIdx  = waterEnabled ? PVOffset + 1 : PVOffset + 0;
+
+    /*!
+     * \brief Index of the switching variable which determines the composition of the
+     *        hydrocarbon phases.
+     *
+     * \note For two-phase water oil models this is disabled.
+     */
+    static const int compositionSwitchIdx = gasEnabled ? PVOffset + 1 : -10000;
+
+    //! Index of the primary variable for the first solvent
+    static const int solventSaturationIdx =
+        enableSolvent ? PVOffset + numPhases : -1000;
+
+    //! Index of the primary variable for the first polymer
+    static const int polymerConcentrationIdx =
+        enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the primary variable for the second polymer primary variable (molecular weight)
+    static const int polymerMoleWeightIdx =
+        numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the foam
+    static const int foamConcentrationIdx =
+        enableFoam ? polymerMoleWeightIdx + 1 : -1000;
+
+    //! Index of the primary variable for temperature
+    static const int temperatureIdx  =
+        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : - 1000;
+
+    //////////////////////
+    // Equation indices
+    //////////////////////
+
+    //! \brief returns the index of "active" component
+    static unsigned canonicalToActiveComponentIndex(unsigned compIdx)
+    {
+        // assumes canonical oil = 0, water = 1, gas = 2;
+        if(!gasEnabled) {
+            assert(compIdx != 2);
+            // oil = 0, water = 1
+            return compIdx;
+        } else if (!waterEnabled) {
+            assert(compIdx != 1);
+            // oil = 0, gas = 1
+            return compIdx / 2;
+        } else {
+            assert(!oilEnabled);
+            assert(compIdx != 0);
+        }
+        // water = 0, gas = 1;
+        return compIdx-1;
+    }
+
+    static unsigned activeToCanonicalComponentIndex(unsigned compIdx)
+    {
+        // assumes canonical oil = 0, water = 1, gas = 2;
+        assert(compIdx < 2);
+        if(!gasEnabled) {
+            // oil = 0, water = 1
+            return compIdx;
+        } else if (!waterEnabled) {
+            // oil = 0, gas = 1
+            return compIdx * 2;
+        } else {
+            assert(!oilEnabled);
+        }
+        // water = 0, gas = 1;
+        return compIdx+1;
+    }
+
+    //! Index of the continuity equation of the first phase
+    static const int conti0EqIdx = PVOffset + 0;
+    // one continuity equation follows
+
+    //! Index of the continuity equation for the first solvent component
+    static const int contiSolventEqIdx =
+        enableSolvent ? PVOffset + numPhases : -1000;
+
+    //! Index of the continuity equation for the first polymer component
+    static const int contiPolymerEqIdx =
+        enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the continuity equation for the second polymer component (molecular weight)
+    static const int contiPolymerMWEqIdx =
+        numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the foam  component
+    static const int contiFoamEqIdx =
+        enableFoam ? contiPolymerMWEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for energy
+    static const int contiEnergyEqIdx =
+        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+};
+
+} // namespace Opm
+
+#endif