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

examples/co2injection_pvs_ecfv.cpp

@@ -29,7 +29,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_pvs_ni_ecfv.cpp

@@ -29,7 +29,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_pvs_ni_vcfv.cpp

@@ -29,7 +29,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_pvs_vcfv.cpp

@@ -29,7 +29,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include "problems/co2injectionproblem.hh"

examples/cuvette_pvs.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/cuvetteproblem.hh"
 
 BEGIN_PROPERTIES

examples/diffusion_pvs.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/diffusionproblem.hh"
 
 BEGIN_PROPERTIES

examples/infiltration_pvs.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/infiltrationproblem.hh"
 
 BEGIN_PROPERTIES

examples/obstacle_pvs.cpp

@@ -30,7 +30,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/obstacleproblem.hh"
 
 BEGIN_PROPERTIES

examples/outflow_pvs.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/outflowproblem.hh"
 
 BEGIN_PROPERTIES

examples/problems/cuvetteproblem.hh

@@ -28,7 +28,7 @@
 #ifndef EWOMS_CUVETTE_PROBLEM_HH
 #define EWOMS_CUVETTE_PROBLEM_HH
 
-#include <ewoms/models/pvs/pvsproperties.hh>
+#include <opm/models/pvs/pvsproperties.hh>
 
 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>

examples/problems/infiltrationproblem.hh

@@ -27,7 +27,7 @@
 #ifndef EWOMS_INFILTRATION_PROBLEM_HH
 #define EWOMS_INFILTRATION_PROBLEM_HH
 
-#include <ewoms/models/pvs/pvsproperties.hh>
+#include <opm/models/pvs/pvsproperties.hh>
 
 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
 #include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>

examples/problems/outflowproblem.hh

@@ -27,7 +27,7 @@
 #ifndef EWOMS_OUTFLOW_PROBLEM_HH
 #define EWOMS_OUTFLOW_PROBLEM_HH
 
-#include <ewoms/models/pvs/pvsproperties.hh>
+#include <opm/models/pvs/pvsproperties.hh>
 
 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
 #include <opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp>

examples/problems/waterairproblem.hh

@@ -28,7 +28,7 @@
 #ifndef EWOMS_WATER_AIR_PROBLEM_HH
 #define EWOMS_WATER_AIR_PROBLEM_HH
 
-#include <ewoms/models/pvs/pvsproperties.hh>
+#include <opm/models/pvs/pvsproperties.hh>
 #include <opm/simulators/linalg/parallelistlbackend.hh>
 
 #include <opm/material/fluidsystems/H2OAirFluidSystem.hpp>

examples/waterair_pvs_ni.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/models/pvs/pvsmodel.hh>
+#include <opm/models/pvs/pvsmodel.hh>
 #include "problems/waterairproblem.hh"
 
 BEGIN_PROPERTIES

opm/models/pvs/pvsboundaryratevector.hh

@@ -0,0 +1,204 @@
+// -*- 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::PvsBoundaryRateVector
+ */
+#ifndef EWOMS_PVS_BOUNDARY_RATE_VECTOR_HH
+#define EWOMS_PVS_BOUNDARY_RATE_VECTOR_HH
+
+#include "pvsproperties.hh"
+
+#include <opm/models/common/energymodule.hh>
+#include <opm/material/common/Valgrind.hpp>
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief Implements a rate vector on the boundary for the fully
+ *        implicit compositional multi-phase primary variable
+ *        switching compositional model.
+ */
+template <class TypeTag>
+class PvsBoundaryRateVector : 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, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    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 { conti0EqIdx = Indices::conti0EqIdx };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+
+    typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+
+public:
+    PvsBoundaryRateVector()
+        : ParentType()
+    {}
+
+    /*!
+     * \copydoc
+     * ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar)
+     */
+    PvsBoundaryRateVector(const Evaluation& value)
+        : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const
+     * ImmiscibleBoundaryRateVector& )
+     */
+    PvsBoundaryRateVector(const PvsBoundaryRateVector& value) = default;
+    PvsBoundaryRateVector& operator=(const PvsBoundaryRateVector& 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) = Evaluation(0.0);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            Evaluation density;
+            if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
+                if (focusDofIdx == interiorDofIdx)
+                    density = fluidState.density(phaseIdx);
+                else
+                    density = Opm::getValue(fluidState.density(phaseIdx));
+            }
+            else if (focusDofIdx == interiorDofIdx)
+                density = insideIntQuants.fluidState().density(phaseIdx);
+            else
+                density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx));
+
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                Evaluation molarity;
+                if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
+                    if (focusDofIdx == interiorDofIdx)
+                        molarity = fluidState.molarity(phaseIdx, compIdx);
+                    else
+                        molarity = Opm::getValue(fluidState.molarity(phaseIdx, compIdx));
+                }
+                else if (focusDofIdx == interiorDofIdx)
+                    molarity = insideIntQuants.fluidState().molarity(phaseIdx, compIdx);
+                else
+                    molarity = Opm::getValue(insideIntQuants.fluidState().molarity(phaseIdx, compIdx));
+
+                // add advective flux of current component in current
+                // phase
+                (*this)[conti0EqIdx + compIdx] += extQuants.volumeFlux(phaseIdx)*molarity;
+            }
+
+            if (enableEnergy) {
+                Evaluation specificEnthalpy;
+                if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
+                    if (focusDofIdx == interiorDofIdx)
+                        specificEnthalpy = fluidState.enthalpy(phaseIdx);
+                    else
+                        specificEnthalpy = Opm::getValue(fluidState.enthalpy(phaseIdx));
+                }
+                else if (focusDofIdx == interiorDofIdx)
+                    specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
+                else
+                    specificEnthalpy = Opm::getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
+
+                Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy;
+                EnergyModule::addToEnthalpyRate(*this, enthalpyRate);
+            }
+        }
+
+        if (enableEnergy)
+            // heat conduction
+            EnergyModule::addToEnthalpyRate(*this, EnergyModule::thermalConductionRate(extQuants));
+
+#ifndef NDEBUG
+        for (unsigned i = 0; i < numEq; ++i)
+            Opm::Valgrind::CheckDefined((*this)[i]);
+#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) {
+            Evaluation& val = this->operator[](eqIdx);
+            val = Toolbox::min(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) {
+            Evaluation& val = this->operator[](eqIdx);
+            val = Toolbox::max(0.0, val);
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleBoundaryRateVector::setNoFlow
+     */
+    void setNoFlow()
+    { (*this) = Evaluation(0.0); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsextensivequantities.hh

@@ -0,0 +1,94 @@
+// -*- 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::PvsExtensiveQuantities
+ */
+#ifndef EWOMS_PVS_EXTENSIVE_QUANTITIES_HH
+#define EWOMS_PVS_EXTENSIVE_QUANTITIES_HH
+
+#include "pvsproperties.hh"
+
+#include <opm/models/common/multiphasebaseextensivequantities.hh>
+#include <opm/models/common/energymodule.hh>
+#include <opm/models/common/diffusionmodule.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ * \ingroup ExtensiveQuantities
+ *
+ * \brief Contains all data which is required to calculate all fluxes at a flux
+ *        integration point for the primary variable switching model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+template <class TypeTag>
+class PvsExtensiveQuantities
+    : public MultiPhaseBaseExtensiveQuantities<TypeTag>
+    , public EnergyExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)>
+    , public DiffusionExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion)>
+{
+    typedef MultiPhaseBaseExtensiveQuantities<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
+    typedef Opm::DiffusionExtensiveQuantities<TypeTag, enableDiffusion> DiffusionExtensiveQuantities;
+
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    typedef Opm::EnergyExtensiveQuantities<TypeTag, enableEnergy> EnergyExtensiveQuantities;
+
+public:
+    /*!
+     * \copydoc ParentType::update
+     */
+    void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    {
+        ParentType::update(elemCtx, scvfIdx, timeIdx);
+        DiffusionExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
+        EnergyExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
+    }
+
+    /*!
+     * \copydoc ParentType::updateBoundary
+     */
+    template <class Context, class FluidState>
+    void updateBoundary(const Context& context,
+                        unsigned bfIdx,
+                        unsigned timeIdx,
+                        const FluidState& fluidState)
+    {
+        ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState);
+        DiffusionExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState);
+        EnergyExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState);
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsindices.hh

@@ -0,0 +1,73 @@
+// -*- 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::PvsIndices
+ */
+#ifndef EWOMS_PVS_INDICES_HH
+#define EWOMS_PVS_INDICES_HH
+
+#include "pvsproperties.hh"
+
+#include <opm/models/common/energymodule.hh>
+
+namespace Opm {
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief The indices for the compositional multi-phase primary
+ *        variable switching model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset>
+class PvsIndices
+    : public EnergyIndices<PVOffset + GET_PROP_VALUE(TypeTag, NumComponents),
+                           GET_PROP_VALUE(TypeTag, EnableEnergy)>
+{
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    typedef Opm::EnergyIndices<PVOffset + numComponents, enableEnergy> EnergyIndices;
+
+public:
+    //! Number of partial differential equations or primary variables, respectively
+    static const int numEq = numComponents + EnergyIndices::numEq_;
+
+    // Primary variable indices
+
+    //! Index for the pressure of the first phase
+    static const int pressure0Idx = PVOffset + 0;
+    //! Index of the either the saturation or the mole
+    //! fraction of the phase with the lowest index
+    static const int switch0Idx = PVOffset + 1;
+
+    // equation indices
+
+    //! Index of the mass conservation equation for the first component
+    static const int conti0EqIdx = PVOffset;
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsintensivequantities.hh

@@ -0,0 +1,297 @@
+// -*- 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::PvsIntensiveQuantities
+ */
+#ifndef EWOMS_PVS_INTENSIVE_QUANTITIES_HH
+#define EWOMS_PVS_INTENSIVE_QUANTITIES_HH
+
+#include "pvsproperties.hh"
+
+#include <opm/models/common/energymodule.hh>
+#include <opm/models/common/diffusionmodule.hh>
+
+#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
+#include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
+#include <opm/material/fluidstates/CompositionalFluidState.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <iostream>
+
+namespace Opm {
+/*!
+ * \ingroup PvsModel
+ * \ingroup IntensiveQuantities
+ *
+ * \brief Contains the quantities which are are constant within a
+ *        finite volume in the compositional multi-phase primary
+ *        variable switching model.
+ */
+template <class TypeTag>
+class PvsIntensiveQuantities
+    : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
+    , public DiffusionIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion) >
+    , public EnergyIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy) >
+    , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
+
+    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, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+    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 { switch0Idx = Indices::switch0Idx };
+    enum { pressure0Idx = Indices::pressure0Idx };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem, Evaluation> MiscibleMultiPhaseComposition;
+    typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem, Evaluation> ComputeFromReferencePhase;
+
+    typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
+    typedef Dune::FieldVector<Evaluation, numPhases> EvalPhaseVector;
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
+
+    typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
+    typedef Opm::DiffusionIntensiveQuantities<TypeTag, enableDiffusion> DiffusionIntensiveQuantities;
+    typedef Opm::EnergyIntensiveQuantities<TypeTag, enableEnergy> EnergyIntensiveQuantities;
+
+public:
+    //! The type of the object returned by the fluidState() method
+    typedef Opm::CompositionalFluidState<Evaluation, FluidSystem> FluidState;
+
+    PvsIntensiveQuantities()
+    { }
+
+    PvsIntensiveQuantities(const PvsIntensiveQuantities& other) = default;
+
+    PvsIntensiveQuantities& operator=(const PvsIntensiveQuantities& other) = default;
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::update
+     */
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    {
+        ParentType::update(elemCtx, dofIdx, timeIdx);
+        EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
+
+        const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        const auto& problem = elemCtx.problem();
+
+        /////////////
+        // set the saturations
+        /////////////
+        Evaluation sumSat = 0.0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            fluidState_.setSaturation(phaseIdx, priVars.explicitSaturationValue(phaseIdx, timeIdx));
+            Opm::Valgrind::CheckDefined(fluidState_.saturation(phaseIdx));
+            sumSat += fluidState_.saturation(phaseIdx);
+        }
+        Opm::Valgrind::CheckDefined(priVars.implicitSaturationIdx());
+        Opm::Valgrind::CheckDefined(sumSat);
+        fluidState_.setSaturation(priVars.implicitSaturationIdx(), 1.0 - sumSat);
+
+        /////////////
+        // set the pressures of the fluid phases
+        /////////////
+
+        // calculate capillary pressure
+        const MaterialLawParams& materialParams =
+            problem.materialLawParams(elemCtx, dofIdx, timeIdx);
+        EvalPhaseVector pC;
+        MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
+
+        // set the absolute phase pressures in the fluid state
+        const Evaluation& p0 = priVars.makeEvaluation(pressure0Idx, timeIdx);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            fluidState_.setPressure(phaseIdx, p0 + (pC[phaseIdx] - pC[0]));
+
+        /////////////
+        // calculate the phase compositions
+        /////////////
+
+        typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+        unsigned lowestPresentPhaseIdx = priVars.lowestPresentPhaseIdx();
+        unsigned numNonPresentPhases = 0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!priVars.phaseIsPresent(phaseIdx))
+                ++numNonPresentPhases;
+        }
+
+        // now comes the tricky part: calculate phase compositions
+        if (numNonPresentPhases == numPhases - 1) {
+            // only one phase is present, i.e. the primary variables
+            // contain the complete composition of the phase
+            Evaluation sumx = 0.0;
+            for (unsigned compIdx = 1; compIdx < numComponents; ++compIdx) {
+                const Evaluation& x = priVars.makeEvaluation(switch0Idx + compIdx - 1, timeIdx);
+                fluidState_.setMoleFraction(lowestPresentPhaseIdx, compIdx, x);
+                sumx += x;
+            }
+
+            // set the mole fraction of the first component
+            fluidState_.setMoleFraction(lowestPresentPhaseIdx, 0, 1 - sumx);
+
+            // calculate the composition of the remaining phases (as
+            // well as the densities of all phases). this is the job
+            // of the "ComputeFromReferencePhase" constraint solver
+            ComputeFromReferencePhase::solve(fluidState_, paramCache,
+                                             lowestPresentPhaseIdx,
+                                             /*setViscosity=*/true,
+                                             /*setEnthalpy=*/false);
+        }
+        else {
+            // create the auxiliary constraints
+            unsigned numAuxConstraints = numComponents + numNonPresentPhases - numPhases;
+            Opm::MMPCAuxConstraint<Evaluation> auxConstraints[numComponents];
+
+            unsigned auxIdx = 0;
+            unsigned switchIdx = 0;
+            for (; switchIdx < numPhases - 1; ++switchIdx) {
+                unsigned compIdx = switchIdx + 1;
+                unsigned switchPhaseIdx = switchIdx;
+                if (switchIdx >= lowestPresentPhaseIdx)
+                    switchPhaseIdx += 1;
+
+                if (!priVars.phaseIsPresent(switchPhaseIdx)) {
+                    auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
+                                               priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
+                    ++auxIdx;
+                }
+            }
+
+            for (; auxIdx < numAuxConstraints; ++auxIdx, ++switchIdx) {
+                unsigned compIdx = numPhases - numNonPresentPhases + auxIdx;
+                auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
+                                           priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
+            }
+
+            // both phases are present, i.e. phase compositions are a result of the the
+            // gas <-> liquid equilibrium. This is the job of the
+            // "MiscibleMultiPhaseComposition" constraint solver
+            MiscibleMultiPhaseComposition::solve(fluidState_, paramCache,
+                                                 priVars.phasePresence(),
+                                                 auxConstraints,
+                                                 numAuxConstraints,
+                                                 /*setViscosity=*/true,
+                                                 /*setEnthalpy=*/false);
+        }
+
+#ifndef NDEBUG
+        // make valgrind happy and set the enthalpies to NaN
+        if (!enableEnergy) {
+            Scalar myNan = std::numeric_limits<Scalar>::quiet_NaN();
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                fluidState_.setEnthalpy(phaseIdx, myNan);
+        }
+#endif
+
+        /////////////
+        // calculate the remaining quantities
+        /////////////
+
+        // calculate relative permeabilities
+        MaterialLaw::relativePermeabilities(relativePermeability_,
+                                            materialParams, fluidState_);
+        Opm::Valgrind::CheckDefined(relativePermeability_);
+
+        // mobilities
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            mobility_[phaseIdx] =
+                relativePermeability_[phaseIdx] / fluidState().viscosity(phaseIdx);
+
+        // porosity
+        porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
+        Opm::Valgrind::CheckDefined(porosity_);
+
+        // intrinsic permeability
+        intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
+
+        // update the quantities specific for the velocity model
+        FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
+
+        // energy related quantities
+        EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
+
+        // update the diffusion specific quantities of the intensive quantities
+        DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
+
+        fluidState_.checkDefined();
+    }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::fluidState
+     */
+    const FluidState& fluidState() const
+    { return fluidState_; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::intrinsicPermeability
+     */
+    const DimMatrix& intrinsicPermeability() const
+    { return intrinsicPerm_; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::relativePermeability
+     */
+    const Evaluation& relativePermeability(unsigned phaseIdx) const
+    { return relativePermeability_[phaseIdx]; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::mobility
+     */
+    const Evaluation& mobility(unsigned phaseIdx) const
+    { return mobility_[phaseIdx]; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::porosity
+     */
+    const Evaluation& porosity() const
+    { return porosity_; }
+
+private:
+    FluidState fluidState_;
+    Evaluation porosity_;
+    DimMatrix intrinsicPerm_;
+    Evaluation relativePermeability_[numPhases];
+    Evaluation mobility_[numPhases];
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvslocalresidual.hh

@@ -0,0 +1,200 @@
+// -*- 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::PvsLocalResidual
+ */
+#ifndef EWOMS_PVS_LOCAL_RESIDUAL_HH
+#define EWOMS_PVS_LOCAL_RESIDUAL_HH
+
+#include "pvsproperties.hh"
+
+#include <opm/models/common/diffusionmodule.hh>
+#include <opm/models/common/energymodule.hh>
+
+#include <opm/material/common/Valgrind.hpp>
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief Element-wise calculation of the local residual for the
+ *        compositional multi-phase primary variable switching model.
+ */
+template <class TypeTag>
+class PvsLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
+{
+    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, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+
+    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
+    typedef Opm::DiffusionModule<TypeTag, enableDiffusion> DiffusionModule;
+
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+
+public:
+    /*!
+     * \copydoc ImmiscibleLocalResidual::addPhaseStorage
+     */
+    template <class LhsEval>
+    void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                         const ElementContext& elemCtx,
+                         unsigned dofIdx,
+                         unsigned timeIdx,
+                         unsigned phaseIdx) const
+    {
+        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
+        const auto& fs = intQuants.fluidState();
+
+        // compute storage term of all components within all phases
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+            unsigned eqIdx = conti0EqIdx + compIdx;
+            storage[eqIdx] +=
+                Toolbox::template decay<LhsEval>(fs.molarity(phaseIdx, compIdx))
+                * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
+                * Toolbox::template decay<LhsEval>(intQuants.porosity());
+        }
+
+        EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
+    }
+
+    /*!
+     * \copydoc ImmiscibleLocalResidual::computeStorage
+     */
+    template <class LhsEval>
+    void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                        const ElementContext& elemCtx,
+                        unsigned dofIdx,
+                        unsigned timeIdx) const
+    {
+        storage = 0.0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+            addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx);
+
+        EnergyModule::addSolidEnergyStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx));
+    }
+
+    /*!
+     * \copydoc ImmiscibleLocalResidual::computeFlux
+     */
+    void computeFlux(RateVector& flux,
+                     const ElementContext& elemCtx,
+                     unsigned scvfIdx,
+                     unsigned timeIdx) const
+    {
+        flux = 0.0;
+        addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+        Opm::Valgrind::CheckDefined(flux);
+
+        addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+        Opm::Valgrind::CheckDefined(flux);
+    }
+
+    /*!
+     * \copydoc ImmiscibleLocalResidual::addAdvectiveFlux
+     */
+    void addAdvectiveFlux(RateVector& flux,
+                          const ElementContext& elemCtx,
+                          unsigned scvfIdx,
+                          unsigned timeIdx) const
+    {
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+
+        unsigned focusDofIdx = elemCtx.focusDofIndex();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            // data attached to upstream and the downstream DOFs
+            // of the current phase
+            unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
+            const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+
+            // this is a bit hacky because it is specific to the element-centered
+            // finite volume scheme. (N.B. that if finite differences are used to
+            // linearize the system of equations, it does not matter.)
+            if (upIdx == focusDofIdx) {
+                Evaluation tmp =
+                    up.fluidState().molarDensity(phaseIdx)
+                    * extQuants.volumeFlux(phaseIdx);
+
+                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                    flux[conti0EqIdx + compIdx] +=
+                        tmp*up.fluidState().moleFraction(phaseIdx, compIdx);
+                }
+            }
+            else {
+                Evaluation tmp =
+                    Toolbox::value(up.fluidState().molarDensity(phaseIdx))
+                    * extQuants.volumeFlux(phaseIdx);
+
+                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                    flux[conti0EqIdx + compIdx] +=
+                        tmp*Toolbox::value(up.fluidState().moleFraction(phaseIdx, compIdx));
+                }
+            }
+        }
+
+        EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+    }
+
+    /*!
+     * \copydoc ImmiscibleLocalResidual::addDiffusiveFlux
+     */
+    void addDiffusiveFlux(RateVector& flux,
+                          const ElementContext& elemCtx,
+                          unsigned scvfIdx,
+                          unsigned timeIdx) const
+    {
+        DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+        EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+    }
+
+    /*!
+     * \copydoc ImmiscibleLocalResidual::computeSource
+     */
+    void computeSource(RateVector& source,
+                       const ElementContext& elemCtx,
+                       unsigned dofIdx,
+                       unsigned timeIdx) const
+    {
+        Opm::Valgrind::SetUndefined(source);
+        elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
+        Opm::Valgrind::CheckDefined(source);
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsmodel.hh

@@ -0,0 +1,599 @@
+// -*- 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::PvsModel
+ */
+#ifndef EWOMS_PVS_MODEL_HH
+#define EWOMS_PVS_MODEL_HH
+
+#include <opm/material/densead/Math.hpp>
+
+#include "pvsproperties.hh"
+#include "pvslocalresidual.hh"
+#include "pvsnewtonmethod.hh"
+#include "pvsprimaryvariables.hh"
+#include "pvsratevector.hh"
+#include "pvsboundaryratevector.hh"
+#include "pvsintensivequantities.hh"
+#include "pvsextensivequantities.hh"
+#include "pvsindices.hh"
+
+#include <opm/models/common/multiphasebasemodel.hh>
+#include <opm/models/common/diffusionmodule.hh>
+#include <opm/models/common/energymodule.hh>
+#include <opm/models/io/vtkcompositionmodule.hh>
+#include <opm/models/io/vtkenergymodule.hh>
+#include <opm/models/io/vtkdiffusionmodule.hh>
+
+#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
+#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace Opm {
+template <class TypeTag>
+class PvsModel;
+}
+
+BEGIN_PROPERTIES
+
+//! The type tag for the isothermal single phase problems
+NEW_TYPE_TAG(PvsModel, INHERITS_FROM(MultiPhaseBaseModel,
+                                     VtkPhasePresence,
+                                     VtkComposition,
+                                     VtkEnergy,
+                                     VtkDiffusion));
+
+//! Use the PVS local jacobian operator for the PVS model
+SET_TYPE_PROP(PvsModel,
+              LocalResidual,
+              Opm::PvsLocalResidual<TypeTag>);
+
+//! Use the PVS specific newton method for the PVS model
+SET_TYPE_PROP(PvsModel, NewtonMethod, Opm::PvsNewtonMethod<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(PvsModel, Model, Opm::PvsModel<TypeTag>);
+
+//! the PrimaryVariables property
+SET_TYPE_PROP(PvsModel, PrimaryVariables, Opm::PvsPrimaryVariables<TypeTag>);
+
+//! the RateVector property
+SET_TYPE_PROP(PvsModel, RateVector, Opm::PvsRateVector<TypeTag>);
+
+//! the BoundaryRateVector property
+SET_TYPE_PROP(PvsModel, BoundaryRateVector, Opm::PvsBoundaryRateVector<TypeTag>);
+
+//! the IntensiveQuantities property
+SET_TYPE_PROP(PvsModel, IntensiveQuantities, Opm::PvsIntensiveQuantities<TypeTag>);
+
+//! the ExtensiveQuantities property
+SET_TYPE_PROP(PvsModel, ExtensiveQuantities, Opm::PvsExtensiveQuantities<TypeTag>);
+
+//! The indices required by the isothermal PVS model
+SET_TYPE_PROP(PvsModel, Indices, Opm::PvsIndices<TypeTag, /*PVIdx=*/0>);
+
+// set the model to a medium verbosity
+SET_INT_PROP(PvsModel, PvsVerbosity, 1);
+
+//! Disable the energy equation by default
+SET_BOOL_PROP(PvsModel, EnableEnergy, false);
+
+// disable molecular diffusion by default
+SET_BOOL_PROP(PvsModel, EnableDiffusion, false);
+
+//! The basis value for the weight of the pressure primary variable
+SET_SCALAR_PROP(PvsModel, PvsPressureBaseWeight, 1.0);
+
+//! The basis value for the weight of the saturation primary variables
+SET_SCALAR_PROP(PvsModel, PvsSaturationsBaseWeight, 1.0);
+
+//! The basis value for the weight of the mole fraction primary variables
+SET_SCALAR_PROP(PvsModel, PvsMoleFractionsBaseWeight, 1.0);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief A generic compositional multi-phase model using primary-variable
+ *        switching.
+ *
+ * This model assumes a flow of \f$M \geq 1\f$ fluid phases
+ * \f$\alpha\f$, each of which is assumed to be a mixture \f$N \geq
+ * M\f$ chemical species \f$\kappa\f$.
+ *
+ * By default, the standard multi-phase Darcy approach is used to determine
+ * the velocity, 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 core of the model is the conservation mass of each component by
+ * means of the equation
+ * \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]
+ *
+ * To close the system mathematically, \f$M\f$ model equations are
+ * also required. This model uses the primary variable switching
+ * assumptions, which are given by:
+ * \f[
+ * 0 \stackrel{!}{=}
+ *  f_\alpha = \left\{
+ *  \begin{array}{cl}
+ *    S_\alpha&  \quad \text{if phase }\alpha\text{ is not present} \    \
+ *    1 - \sum_\kappa x_\alpha^\kappa&  \quad \text{else}
+ *  \end{array}
+ *  \right.
+ * \f]
+ *
+ * To make this approach applicable, a pseudo primary variable
+ * <em>phase presence</em> has to be introduced. Its purpose is to
+ * specify for each phase whether it is present or not. It is a
+ * <em>pseudo</em> primary variable because it is not directly considered when
+ * linearizing the system in the Newton method, but after each Newton
+ * iteration, it gets updated like the "real" primary variables.  The
+ * following rules are used for this update procedure:
+ *
+ * <ul>
+ * <li>If phase \f$\alpha\f$ is present according to the pseudo
+ *     primary variable, but \f$S_\alpha < 0\f$ after the Newton
+ *     update, consider the phase \f$\alpha\f$ disappeared for the
+ *     next iteration and use the set of primary variables which
+ *     correspond to the new phase presence.</li>
+ * <li>If phase \f$\alpha\f$ is not present according to the pseudo
+ *     primary variable, but the sum of the component mole fractions
+ *     in the phase is larger than 1, i.e. \f$\sum_\kappa
+ *     x_\alpha^\kappa > 1\f$, consider the phase \f$\alpha\f$ present
+ *     in the the next iteration and update the set of primary
+ *     variables to make it consistent with the new phase
+ *     presence.</li>
+ * <li>In all other cases don't modify the phase presence for phase
+ *     \f$\alpha\f$.</li>
+ *
+ * </ul>
+ *
+ * The model always requires \f$N\f$ primary variables, but their
+ * interpretation is dependent on the phase presence:
+ *
+ * <ul>
+ *
+ * <li>The first primary variable is always interpreted as the
+ *      pressure of the phase with the lowest index \f$PV_0 =
+ *      p_0\f$.</li>
+ *
+ * <li>Then, \f$M - 1\f$ "switching primary variables" follow, which
+ *     are interpreted depending in the presence of the first
+ *     \f$M-1\f$ phases: If phase \f$\alpha\f$ is present, its
+ *     saturation \f$S_\alpha = PV_i\f$ is used as primary variable;
+ *     if it is not present, the mole fraction \f$PV_i =
+ *     x_{\alpha^\star}^\alpha\f$ of the component with index
+ *     \f$\alpha\f$ in the phase with the lowest index that is present
+ *     \f$\alpha^\star\f$ is used instead.</li>
+ *
+ * <li>Finally, the mole fractions of the \f$N-M\f$ components with
+ *     the largest index in the phase with the lowest index that is
+ *     present \f$x_{\alpha^\star}^\kappa\f$ are used as primary
+ *     variables.</li>
+ *
+ * </ul>
+ */
+template <class TypeTag>
+class PvsModel
+    : public MultiPhaseBaseModel<TypeTag>
+{
+    typedef MultiPhaseBaseModel<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+    typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
+
+public:
+    PvsModel(Simulator& simulator)
+        : ParentType(simulator)
+    {
+        verbosity_ = EWOMS_GET_PARAM(TypeTag, int, PvsVerbosity);
+        numSwitched_ = 0;
+    }
+
+    /*!
+     * \brief Register all run-time parameters for the PVS compositional model.
+     */
+    static void registerParameters()
+    {
+        ParentType::registerParameters();
+
+        // register runtime parameters of the VTK output modules
+        Opm::VtkPhasePresenceModule<TypeTag>::registerParameters();
+        Opm::VtkCompositionModule<TypeTag>::registerParameters();
+
+        if (enableDiffusion)
+            Opm::VtkDiffusionModule<TypeTag>::registerParameters();
+
+        if (enableEnergy)
+            Opm::VtkEnergyModule<TypeTag>::registerParameters();
+
+        EWOMS_REGISTER_PARAM(TypeTag, int, PvsVerbosity,
+                             "The verbosity level of the primary variable "
+                             "switching model");
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::name
+     */
+    static std::string name()
+    { return "pvs"; }
+
+    /*!
+     * \copydoc FvBaseDiscretization::primaryVarName
+     */
+    std::string primaryVarName(unsigned pvIdx) const
+    {
+        std::string s;
+        if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
+            return s;
+
+        std::ostringstream oss;
+        if (pvIdx == Indices::pressure0Idx)
+            oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
+        else if (Indices::switch0Idx <= pvIdx
+                 && pvIdx < Indices::switch0Idx + numPhases - 1)
+            oss << "switch_" << pvIdx - Indices::switch0Idx;
+        else if (Indices::switch0Idx + numPhases - 1 <= pvIdx
+                 && pvIdx < Indices::switch0Idx + numComponents - 1)
+            oss << "auxMoleFrac^" << FluidSystem::componentName(pvIdx);
+        else
+            assert(false);
+
+        return oss.str();
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::eqName
+     */
+    std::string eqName(unsigned eqIdx) const
+    {
+        std::string s;
+        if (!(s = EnergyModule::eqName(eqIdx)).empty())
+            return s;
+
+        std::ostringstream oss;
+        if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx
+                                                     + numComponents) {
+            unsigned compIdx = eqIdx - Indices::conti0EqIdx;
+            oss << "continuity^" << FluidSystem::componentName(compIdx);
+        }
+        else
+            assert(false);
+
+        return oss.str();
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::updateFailed
+     */
+    void updateFailed()
+    {
+        ParentType::updateFailed();
+        numSwitched_ = 0;
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::updateBegin
+     */
+    void updateBegin()
+    {
+        ParentType::updateBegin();
+
+        // find the a reference pressure. The first degree of freedom
+        // might correspond to non-interior entities which would lead
+        // to an undefined value, so we have to iterate...
+        size_t nDof = this->numTotalDof();
+        for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
+            if (this->dofTotalVolume(dofIdx) > 0.0) {
+                referencePressure_ =
+                    this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
+                if (referencePressure_ > 0.0)
+                    break;
+            }
+        }
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::primaryVarWeight
+     */
+    Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
+    {
+        Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
+        if (tmp > 0)
+            // energy related quantity
+            return tmp;
+
+        if (Indices::pressure0Idx == pvIdx) {
+            return 10 / referencePressure_;
+        }
+
+        if (Indices::switch0Idx <= pvIdx && pvIdx < Indices::switch0Idx
+                                                    + numPhases - 1) {
+            unsigned phaseIdx = pvIdx - Indices::switch0Idx;
+
+            if (!this->solution(/*timeIdx=*/0)[globalDofIdx].phaseIsPresent(phaseIdx))
+                // for saturations, the weight is always 1
+                return 1;
+
+            // for saturations, the PvsMoleSaturationsBaseWeight
+            // property determines the weight
+            return GET_PROP_VALUE(TypeTag, PvsSaturationsBaseWeight);
+        }
+
+        // for mole fractions, the PvsMoleFractionsBaseWeight
+        // property determines the weight
+        return GET_PROP_VALUE(TypeTag, PvsMoleFractionsBaseWeight);
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::eqWeight
+     */
+    Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
+    {
+        Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
+        if (tmp > 0)
+            // energy related equation
+            return tmp;
+
+        unsigned compIdx = eqIdx - Indices::conti0EqIdx;
+        assert(0 <= compIdx && compIdx <= numComponents);
+
+        // make all kg equal
+        return FluidSystem::molarMass(compIdx);
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::advanceTimeLevel
+     */
+    void advanceTimeLevel()
+    {
+        ParentType::advanceTimeLevel();
+        numSwitched_ = 0;
+    }
+
+    /*!
+     * \brief Return true if the primary variables were switched for
+     *        at least one vertex after the last timestep.
+     */
+    bool switched() const
+    { return numSwitched_ > 0; }
+
+    /*!
+     * \copydoc FvBaseDiscretization::serializeEntity
+     */
+    template <class DofEntity>
+    void serializeEntity(std::ostream& outstream, const DofEntity& dofEntity)
+    {
+        // write primary variables
+        ParentType::serializeEntity(outstream, dofEntity);
+
+        unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
+        if (!outstream.good())
+            throw std::runtime_error("Could not serialize DOF "+std::to_string(dofIdx));
+
+        outstream << this->solution(/*timeIdx=*/0)[dofIdx].phasePresence() << " ";
+    }
+
+    /*!
+     * \copydoc FvBaseDiscretization::deserializeEntity
+     */
+    template <class DofEntity>
+    void deserializeEntity(std::istream& instream, const DofEntity& dofEntity)
+    {
+        // read primary variables
+        ParentType::deserializeEntity(instream, dofEntity);
+
+        // read phase presence
+        unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
+        if (!instream.good())
+            throw std::runtime_error("Could not deserialize DOF "+std::to_string(dofIdx));
+
+        short tmp;
+        instream >> tmp;
+        this->solution(/*timeIdx=*/0)[dofIdx].setPhasePresence(tmp);
+        this->solution(/*timeIdx=*/1)[dofIdx].setPhasePresence(tmp);
+    }
+
+    /*!
+     * \internal
+     * \brief Do the primary variable switching after a Newton iteration.
+     *
+     * This is an internal method that needs to be public because it
+     * gets called by the Newton method after an update.
+     */
+    void switchPrimaryVars_()
+    {
+        numSwitched_ = 0;
+
+        int succeeded;
+        try {
+            std::vector<bool> visited(this->numGridDof(), false);
+            ElementContext elemCtx(this->simulator_);
+
+            ElementIterator elemIt = this->gridView_.template begin<0>();
+            ElementIterator elemEndIt = this->gridView_.template end<0>();
+            for (; elemIt != elemEndIt; ++elemIt) {
+                const Element& elem = *elemIt;
+                if (elem.partitionType() != Dune::InteriorEntity)
+                    continue;
+                elemCtx.updateStencil(elem);
+
+                size_t numLocalDof = elemCtx.stencil(/*timeIdx=*/0).numPrimaryDof();
+                for (unsigned dofIdx = 0; dofIdx < numLocalDof; ++dofIdx) {
+                    unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+
+                    if (visited[globalIdx])
+                        continue;
+                    visited[globalIdx] = true;
+
+                    // compute the intensive quantities of the current degree of freedom
+                    auto& priVars = this->solution(/*timeIdx=*/0)[globalIdx];
+                    elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
+                    const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
+
+                    // evaluate primary variable switch
+                    short oldPhasePresence = priVars.phasePresence();
+
+                    // set the primary variables and the new phase state
+                    // from the current fluid state
+                    priVars.assignNaive(intQuants.fluidState());
+
+                    if (oldPhasePresence != priVars.phasePresence()) {
+                        if (verbosity_ > 1)
+                            printSwitchedPhases_(elemCtx,
+                                                 dofIdx,
+                                                 intQuants.fluidState(),
+                                                 oldPhasePresence,
+                                                 priVars);
+                        ++numSwitched_;
+                    }
+                }
+            }
+
+            succeeded = 1;
+        }
+        catch (...)
+        {
+            std::cout << "rank " << this->simulator_.gridView().comm().rank()
+                      << " caught an exception during primary variable switching"
+                      << "\n"  << std::flush;
+            succeeded = 0;
+        }
+        succeeded = this->simulator_.gridView().comm().min(succeeded);
+
+        if (!succeeded)
+            throw Opm::NumericalIssue("A process did not succeed in adapting the primary variables");
+
+        // make sure that if there was a variable switch in an
+        // other partition we will also set the switch flag
+        // for our partition.
+        numSwitched_ = this->gridView_.comm().sum(numSwitched_);
+
+        if (verbosity_ > 0)
+            this->simulator_.model().newtonMethod().endIterMsg()
+                << ", num switched=" << numSwitched_;
+    }
+
+    template <class FluidState>
+    void printSwitchedPhases_(const ElementContext& elemCtx,
+                              unsigned dofIdx,
+                              const FluidState& fs,
+                              short oldPhasePresence,
+                              const PrimaryVariables& newPv) const
+    {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            bool oldPhasePresent = (oldPhasePresence&  (1 << phaseIdx)) > 0;
+            bool newPhasePresent = newPv.phaseIsPresent(phaseIdx);
+            if (oldPhasePresent == newPhasePresent)
+                continue;
+
+            const auto& pos = elemCtx.pos(dofIdx, /*timeIdx=*/0);
+            if (oldPhasePresent && !newPhasePresent) {
+                std::cout << "'" << FluidSystem::phaseName(phaseIdx)
+                          << "' phase disappears at position " << pos
+                          << ". saturation=" << fs.saturation(phaseIdx)
+                          << std::flush;
+            }
+            else {
+                Scalar sumx = 0;
+                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
+                    sumx += FsToolbox::value(fs.moleFraction(phaseIdx, compIdx));
+
+                std::cout << "'" << FluidSystem::phaseName(phaseIdx)
+                          << "' phase appears at position " << pos
+                          << " sum x = " << sumx  << std::flush;
+            }
+        }
+
+        std::cout << ", new primary variables: ";
+        newPv.print();
+        std::cout << "\n"  << std::flush;
+    }
+
+    void registerOutputModules_()
+    {
+        ParentType::registerOutputModules_();
+
+        // add the VTK output modules which are meaningful for the model
+        this->addOutputModule(new Opm::VtkPhasePresenceModule<TypeTag>(this->simulator_));
+        this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
+        if (enableDiffusion)
+            this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
+        if (enableEnergy)
+            this->addOutputModule(new Opm::VtkEnergyModule<TypeTag>(this->simulator_));
+    }
+
+    mutable Scalar referencePressure_;
+
+    // number of switches of the phase state in the last Newton
+    // iteration
+    unsigned numSwitched_;
+
+    // verbosity of the model
+    int verbosity_;
+};
+}
+
+#endif

opm/models/pvs/pvsnewtonmethod.hh

@@ -0,0 +1,129 @@
+// -*- 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::PvsNewtonMethod
+ */
+#ifndef EWOMS_PVS_NEWTON_METHOD_HH
+#define EWOMS_PVS_NEWTON_METHOD_HH
+
+#include "pvsproperties.hh"
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief A newton solver which is specific to the compositional
+ *        multi-phase PVS model.
+ */
+template <class TypeTag>
+class PvsNewtonMethod : 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, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+    enum { numPhases = FluidSystem::numPhases };
+
+    // primary variable indices
+    enum { pressure0Idx = Indices::pressure0Idx };
+    enum { switch0Idx = Indices::switch0Idx };
+
+public:
+    PvsNewtonMethod(Simulator& simulator) : ParentType(simulator)
+    {}
+
+protected:
+    friend NewtonMethod<TypeTag>;
+    friend ParentType;
+
+    /*!
+     * \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
+     */
+    void updatePrimaryVariables_(unsigned globalDofIdx OPM_UNUSED,
+                                 PrimaryVariables& nextValue,
+                                 const PrimaryVariables& currentValue,
+                                 const EqVector& update,
+                                 const EqVector& currentResidual OPM_UNUSED)
+    {
+        // normal Newton-Raphson update
+        nextValue = currentValue;
+        nextValue -= update;
+
+        ////
+        // put crash barriers along the update path
+        ////
+        // saturations: limit the change of any saturation to at most 20%
+        Scalar sumSatDelta = 0.0;
+        Scalar maxSatDelta = 0.0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
+            if (!currentValue.phaseIsPresent(phaseIdx))
+                continue;
+
+            maxSatDelta = std::max(std::abs(update[switch0Idx + phaseIdx]),
+                                   maxSatDelta);
+            sumSatDelta += update[switch0Idx + phaseIdx];
+        }
+        maxSatDelta = std::max(std::abs(- sumSatDelta), maxSatDelta);
+
+        if (maxSatDelta > 0.2) {
+            Scalar alpha = 0.2/maxSatDelta;
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
+                if (!currentValue.phaseIsPresent(phaseIdx))
+                    continue;
+
+                nextValue[switch0Idx + phaseIdx] =
+                    currentValue[switch0Idx + phaseIdx]
+                    - alpha*update[switch0Idx + phaseIdx];
+            }
+        }
+
+        // limit pressure reference change to 20% of the total value per iteration
+        clampValue_(nextValue[pressure0Idx],
+                    currentValue[pressure0Idx]*0.8,
+                    currentValue[pressure0Idx]*1.2);
+    }
+
+    /*!
+     * \copydoc NewtonMethod::endIteration_
+     */
+    void endIteration_(SolutionVector& uCurrentIter,
+                       const SolutionVector& uLastIter)
+    {
+        ParentType::endIteration_(uCurrentIter, uLastIter);
+        this->problem().model().switchPrimaryVars_();
+    }
+
+    void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
+    { val = std::max(minVal, std::min(val, maxVal)); }
+};
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsprimaryvariables.hh

@@ -0,0 +1,378 @@
+// -*- 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::PvsPrimaryVariables
+ */
+#ifndef EWOMS_PVS_PRIMARY_VARIABLES_HH
+#define EWOMS_PVS_PRIMARY_VARIABLES_HH
+
+#include "pvsindices.hh"
+#include "pvsproperties.hh"
+
+#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
+#include <opm/models/common/energymodule.hh>
+
+#include <opm/material/constraintsolvers/NcpFlash.hpp>
+#include <opm/material/fluidstates/CompositionalFluidState.hpp>
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <iostream>
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief Represents the primary variables used in the primary
+ *        variable switching compositional model.
+ *
+ * This class is basically a Dune::FieldVector which can retrieve its
+ * contents from an aribitatry fluid state.
+ */
+template <class TypeTag>
+class PvsPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
+{
+    typedef FvBasePrimaryVariables<TypeTag> ParentType;
+    typedef PvsPrimaryVariables<TypeTag> ThisType;
+    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, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    // primary variable indices
+    enum { pressure0Idx = Indices::pressure0Idx };
+    enum { switch0Idx = Indices::switch0Idx };
+
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+
+    typedef typename Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
+    typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
+    typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
+
+public:
+    PvsPrimaryVariables() : ParentType()
+    { Opm::Valgrind::SetDefined(*this); }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar)
+     */
+    explicit PvsPrimaryVariables(Scalar value) : ParentType(value)
+    {
+        Opm::Valgrind::CheckDefined(value);
+        Opm::Valgrind::SetDefined(*this);
+
+        phasePresence_ = 0;
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const
+     * ImmisciblePrimaryVariables& )
+     */
+    PvsPrimaryVariables(const PvsPrimaryVariables& value) : ParentType(value)
+    {
+        Opm::Valgrind::SetDefined(*this);
+
+        phasePresence_ = value.phasePresence_;
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::assignMassConservative
+     */
+    template <class FluidState>
+    void assignMassConservative(const FluidState& fluidState,
+                                const MaterialLawParams& matParams,
+                                bool isInEquilibrium = false)
+    {
+#ifndef NDEBUG
+        // make sure the temperature is the same in all fluid phases
+        for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
+            assert(std::abs(fluidState.temperature(0) - fluidState.temperature(phaseIdx)) < 1e-30);
+        }
+#endif // NDEBUG
+
+        // for the equilibrium case, we don't need complicated
+        // computations.
+        if (isInEquilibrium) {
+            assignNaive(fluidState);
+            return;
+        }
+
+        // use a flash calculation to calculate a fluid state in
+        // thermodynamic equilibrium
+        typename FluidSystem::template ParameterCache<Scalar> paramCache;
+        Opm::CompositionalFluidState<Scalar, FluidSystem> fsFlash;
+
+        // use the externally given fluid state as initial value for
+        // the flash calculation
+        fsFlash.assign(fluidState);
+
+        // calculate the phase densities
+        paramCache.updateAll(fsFlash);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            Scalar rho = FluidSystem::density(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) {
+                globalMolarities[compIdx] +=
+                    fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
+            }
+        }
+
+        // run the flash calculation
+        NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
+
+        // use the result to assign the primary variables
+        assignNaive(fsFlash);
+    }
+
+    /*!
+     * \brief Return the fluid phases which are present in a given
+     *        control volume.
+     */
+    short phasePresence() const
+    { return phasePresence_; }
+
+    /*!
+     * \brief Set which fluid phases are present in a given control volume.
+     *
+     * \param value The new phase presence. The phase with index i is
+     *              present if the i-th bit of \c value is 1.
+     */
+    void setPhasePresence(short value)
+    { phasePresence_ = value; }
+
+    /*!
+     * \brief Set whether a given indivividual phase should be present
+     *        or not.
+     *
+     * \param phaseIdx The index of the phase which's presence ought to be set or reset.
+     * \param yesno If true, the presence of the phase is set, else it is reset
+     */
+    void setPhasePresent(unsigned phaseIdx, bool yesno = true)
+    {
+        if (yesno)
+            setPhasePresence(phasePresence_ | (1 << phaseIdx));
+        else
+            setPhasePresence(phasePresence_&  ~(1 << phaseIdx));
+    }
+
+    /*!
+     * \brief Returns the index of the phase with's its saturation is
+     *        determined by the closure condition of saturation.
+     */
+    unsigned implicitSaturationIdx() const
+    { return lowestPresentPhaseIdx(); }
+
+    /*!
+     * \brief Returns true iff a phase is present for a given phase
+     *        presence.
+     *
+     * \param phaseIdx The index of the phase which's presence is
+     *                 queried.
+     * \param phasePresence The bit-map of present phases.
+     */
+    static bool phaseIsPresent(unsigned phaseIdx, short phasePresence)
+    { return phasePresence&  (1 << phaseIdx); }
+
+    /*!
+     * \brief Returns true iff a phase is present for the current
+     *        phase presence.
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     */
+    bool phaseIsPresent(unsigned phaseIdx) const
+    { return phasePresence_&  (1 << phaseIdx); }
+
+    /*!
+     * \brief Returns the phase with the lowest index that is present.
+     */
+    unsigned lowestPresentPhaseIdx() const
+    { return static_cast<unsigned>(ffs(phasePresence_) - 1); }
+
+    /*!
+     * \brief Assignment operator from an other primary variables object
+     */
+    ThisType& operator=(const Implementation& value)
+    {
+        ParentType::operator=(value);
+        phasePresence_ = value.phasePresence_;
+
+        return *this;
+    }
+
+    /*!
+     * \brief Assignment operator from a scalar value
+     */
+    ThisType& operator=(Scalar value)
+    {
+        ParentType::operator=(value);
+
+        phasePresence_ = 0;
+        return *this;
+    }
+
+    /*!
+     * \brief Returns an explcitly stored saturation for a given phase.
+     *
+     * (or 0 if the saturation is not explicitly stored.)
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     */
+    Evaluation explicitSaturationValue(unsigned phaseIdx, unsigned timeIdx) const
+    {
+        if (!phaseIsPresent(phaseIdx) || phaseIdx == lowestPresentPhaseIdx())
+            // non-present phases have saturation 0
+            return 0.0;
+
+        unsigned varIdx = switch0Idx + phaseIdx - 1;
+        if (std::is_same<Evaluation, Scalar>::value)
+            return (*this)[varIdx]; // finite differences
+        else {
+            // automatic differentiation
+            if (timeIdx != 0)
+                Toolbox::createConstant((*this)[varIdx]);
+            return Toolbox::createVariable((*this)[varIdx], varIdx);
+        }
+    }
+
+    /*!
+     * \copydoc ImmisciblePrimaryVariables::assignNaive
+     */
+    template <class FluidState>
+    void assignNaive(const FluidState& fluidState)
+    {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        // assign the phase temperatures. this is out-sourced to
+        // the energy module
+        EnergyModule::setPriVarTemperatures(*this, fluidState);
+
+        // set the pressure of the first phase
+        (*this)[pressure0Idx] = FsToolbox::value(fluidState.pressure(/*phaseIdx=*/0));
+        Opm::Valgrind::CheckDefined((*this)[pressure0Idx]);
+
+        // determine the phase presence.
+        phasePresence_ = 0;
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            // use a NCP condition to determine if the phase is
+            // present or not
+            Scalar a = 1;
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                a -= FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx));
+            }
+            Scalar b = FsToolbox::value(fluidState.saturation(phaseIdx));
+
+            if (b > a)
+                phasePresence_ |= (1 << phaseIdx);
+        }
+
+        // some phase must be present
+        if (phasePresence_ == 0)
+            throw Opm::NumericalIssue("Phase state was 0, i.e., no fluid is present");
+
+        // set the primary variables which correspond to mole
+        // fractions of the present phase which has the lowest index.
+        unsigned lowestPhaseIdx = lowestPresentPhaseIdx();
+        for (unsigned switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
+            unsigned phaseIdx = switchIdx;
+            unsigned compIdx = switchIdx + 1;
+            if (switchIdx >= lowestPhaseIdx)
+                ++phaseIdx;
+
+            if (phaseIsPresent(phaseIdx)) {
+                (*this)[switch0Idx + switchIdx] = FsToolbox::value(fluidState.saturation(phaseIdx));
+                Opm::Valgrind::CheckDefined((*this)[switch0Idx + switchIdx]);
+            }
+            else {
+                (*this)[switch0Idx + switchIdx] =
+                    FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx));
+                Opm::Valgrind::CheckDefined((*this)[switch0Idx + switchIdx]);
+            }
+        }
+
+        // set the mole fractions in of the remaining components in
+        // the phase with the lowest index
+        for (unsigned compIdx = numPhases - 1; compIdx < numComponents - 1; ++compIdx) {
+            (*this)[switch0Idx + compIdx] =
+                FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx + 1));
+            Opm::Valgrind::CheckDefined((*this)[switch0Idx + compIdx]);
+        }
+    }
+
+    /*!
+     * \copydoc FlashPrimaryVariables::print
+     */
+    void print(std::ostream& os = std::cout) const
+    {
+        os << "(p_" << FluidSystem::phaseName(0) << " = "
+           << this->operator[](pressure0Idx);
+        unsigned lowestPhaseIdx = lowestPresentPhaseIdx();
+        for (unsigned switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
+            unsigned phaseIdx = switchIdx;
+            unsigned compIdx = switchIdx + 1;
+            if (phaseIdx >= lowestPhaseIdx)
+                ++phaseIdx; // skip the saturation of the present
+                            // phase with the lowest index
+
+            if (phaseIsPresent(phaseIdx)) {
+                os << ", S_" << FluidSystem::phaseName(phaseIdx) << " = "
+                   << (*this)[switch0Idx + switchIdx];
+            }
+            else {
+                os << ", x_" << FluidSystem::phaseName(lowestPhaseIdx) << "^"
+                   << FluidSystem::componentName(compIdx) << " = "
+                   << (*this)[switch0Idx + switchIdx];
+            }
+        }
+        for (unsigned compIdx = numPhases - 1; compIdx < numComponents - 1;
+             ++compIdx) {
+            os << ", x_" << FluidSystem::phaseName(lowestPhaseIdx) << "^"
+               << FluidSystem::componentName(compIdx + 1) << " = "
+               << (*this)[switch0Idx + compIdx];
+        }
+        os << ")";
+        os << ", phase presence: " << static_cast<int>(phasePresence_) << std::flush;
+    }
+
+private:
+    short phasePresence_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/pvs/pvsproperties.hh

@@ -0,0 +1,59 @@
+// -*- 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 PvsModel
+ *
+ * \brief Declares the properties required for the compositional
+ *        multi-phase primary variable switching model.
+ */
+#ifndef EWOMS_PVS_PROPERTIES_HH
+#define EWOMS_PVS_PROPERTIES_HH
+
+#include <opm/models/common/multiphasebaseproperties.hh>
+#include <opm/models/common/diffusionmodule.hh>
+#include <opm/models/common/energymodule.hh>
+#include <opm/models/io/vtkcompositionmodule.hh>
+#include <opm/models/io/vtkphasepresencemodule.hh>
+#include <opm/models/io/vtkdiffusionmodule.hh>
+#include <opm/models/io/vtkenergymodule.hh>
+
+BEGIN_PROPERTIES
+
+//! Specifies whether energy is considered as a conservation quantity or not
+NEW_PROP_TAG(EnableEnergy);
+//! Enable diffusive fluxes?
+NEW_PROP_TAG(EnableDiffusion);
+
+//! The verbosity of the model (0 -> do not print anything, 2 -> spam stdout a lot)
+NEW_PROP_TAG(PvsVerbosity);
+//! The basis value for the weight of the pressure primary variable
+NEW_PROP_TAG(PvsPressureBaseWeight);
+//! The basis value for the weight of the saturation primary variables
+NEW_PROP_TAG(PvsSaturationsBaseWeight);
+//! The basis value for the weight of the mole fraction primary variables
+NEW_PROP_TAG(PvsMoleFractionsBaseWeight);
+
+END_PROPERTIES
+
+#endif

opm/models/pvs/pvsratevector.hh

@@ -0,0 +1,148 @@
+// -*- 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::PvsRateVector
+ */
+#ifndef EWOMS_PVS_RATE_VECTOR_HH
+#define EWOMS_PVS_RATE_VECTOR_HH
+
+#include "pvsindices.hh"
+
+#include <opm/models/common/energymodule.hh>
+#include <opm/material/constraintsolvers/NcpFlash.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup PvsModel
+ *
+ * \brief Implements a vector representing molar, mass or volumetric rates.
+ *
+ * This class is basically a Dune::FieldVector which can be set using either mass, molar
+ * or volumetric rates.
+ */
+template <class TypeTag>
+class PvsRateVector
+    : 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;
+
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+
+    typedef Dune::FieldVector<Evaluation, numEq> ParentType;
+    typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
+
+public:
+    PvsRateVector() : ParentType()
+    { Opm::Valgrind::SetUndefined(*this); }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar)
+     */
+    PvsRateVector(const Evaluation& value)
+        : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const ImmiscibleRateVector&)
+     */
+    PvsRateVector(const PvsRateVector& value) : ParentType(value)
+    {}
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setMassRate
+     */
+    void setMassRate(const ParentType& value)
+    {
+        // convert to molar rates
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+            (*this)[conti0EqIdx + compIdx] = value[conti0EqIdx + compIdx];
+            (*this)[conti0EqIdx + compIdx] /= FluidSystem::molarMass(compIdx);
+        }
+    }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setMolarRate
+     */
+    void setMolarRate(const ParentType& value)
+    { ParentType::operator=(value); }
+
+    /*!
+     * \copydoc ImmiscibleRateVector::setEnthalpyRate
+     */
+    template <class RhsEval>
+    void setEnthalpyRate(const RhsEval& rate)
+    { EnergyModule::setEnthalpyRate(*this, rate); }
+
+    /*!
+     * \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, compIdx)
+                * fluidState.moleFraction(phaseIdx, compIdx)
+                * volume;
+
+        EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume);
+    }
+
+    /*!
+     * \brief Assignment operator from a scalar or a function evaluation
+     */
+    template <class RhsEval>
+    PvsRateVector& operator=(const RhsEval& value)
+    {
+        for (unsigned i=0; i < this->size(); ++i)
+            (*this)[i] = value;
+        return *this;
+    }
+
+    /*!
+     * \brief Assignment operator from another rate vector
+     */
+    PvsRateVector& operator=(const PvsRateVector& other)
+    {
+        for (unsigned i=0; i < this->size(); ++i)
+            (*this)[i] = other[i];
+        return *this;
+    }
+};
+
+} // namespace Opm
+
+#endif