Merge pull request #530 from atgeirr/enable-onephase-revised

Enable onephase revised

opm/models/blackoil/blackoilenergymodules.hh

@@ -376,7 +376,6 @@ public:
         // and the thermal condictivity coefficients
         for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
-                fs.setEnthalpy(phaseIdx, 0.0);
                 continue;
             }
 

opm/models/blackoil/blackoilintensivequantities.hh

@@ -128,16 +128,21 @@ public:
 
         // extract the water and the gas saturations for convenience
         Evaluation Sw = 0.0;
-        if (waterEnabled)
+        if (waterEnabled) {
+            if (priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) {
+                Sw = 1.0;
+            } else {
                 Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
-
+            }
+        }
         Evaluation Sg = 0.0;
         if (compositionSwitchEnabled)
         {
-            if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
+            if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
                 // -> threephase case
+                assert( priVars.primaryVarsMeaning() != PrimaryVariables::OnePhase_p );
                 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
-            else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
+            } else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
                 // -> gas-water case
                 Sg = 1.0 - Sw;
 
@@ -201,9 +206,11 @@ public:
         // update the Saturation functions for the blackoil solvent module.
         asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
 
-        const Evaluation& SoMax =
+        Evaluation SoMax = 0.0;
-            Opm::max(fluidState_.saturation(oilPhaseIdx),
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+            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
@@ -257,9 +264,7 @@ public:
             else
                 fluidState_.setRv(0.0);
         }
-        else {
+        else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
-            assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
-
             const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
             fluidState_.setRv(Rv);
 
@@ -274,14 +279,18 @@ public:
                                                             SoMax);
 
                 fluidState_.setRs(Opm::min(RsMax, RsSat));
-            }
+            } else {
-            else
                 fluidState_.setRs(0.0);
             }
+        } else {
+            assert(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p);
+        }
 
         typename FluidSystem::template ParameterCache<Evaluation> paramCache;
         paramCache.setRegionIndex(pvtRegionIdx);
+        if(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
             paramCache.setMaxOilSat(SoMax);
+        }
         paramCache.updateAll(fluidState_);
 
         // compute the phase densities and transform the phase permeabilities into mobilities
@@ -338,7 +347,14 @@ public:
         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);
+            Evaluation x;
+            if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
+            } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
+                x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
+            } else {
+                x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
+            }
             porosity_ *= 1.0 + x + 0.5*x*x;
         }
 

opm/models/blackoil/blackoilonephaseindices.hh

@@ -0,0 +1,163 @@
+// -*- 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 Ewoms::BlackOilTwoPhaseIndices
+ */
+#ifndef EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH
+#define EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH
+
+#include <cassert>
+
+namespace Ewoms {
+
+/*!
+ * \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 canonicalCompIdx>
+struct BlackOilOnePhaseIndices
+{
+    //! Is phase enabled or not
+    static const bool oilEnabled = (canonicalCompIdx == 0);
+    static const bool waterEnabled = (canonicalCompIdx == 1);
+    static const bool gasEnabled = (canonicalCompIdx == 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 = 1;
+
+    //! 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  = -10000;
+
+    //! Index of the oil pressure in a vector of primary variables
+    static const int pressureSwitchIdx  = 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 = -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 ? 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
+    //////////////////////
+
+    //! \brief returns the index of "active" component
+    static unsigned canonicalToActiveComponentIndex(unsigned compIdx)
+    {
+        return 0;
+    }
+
+    static unsigned activeToCanonicalComponentIndex(unsigned compIdx)
+    {
+        // assumes canonical oil = 0, water = 1, gas = 2;
+        assert(compIdx == 0);
+        if(gasEnabled) {
+            return 2;
+        } else if (waterEnabled) {
+            return 1;
+        } else {
+            assert(oilEnabled);
+        }
+        return 0;
+    }
+
+    //! Index of the continuity equation of the first (and only) phase
+    static const int conti0EqIdx = PVOffset + 0;
+
+    //! 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 Ewoms
+
+#endif

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -112,6 +112,7 @@ public:
         Sw_po_Sg, // threephase case
         Sw_po_Rs, // water + oil case
         Sw_pg_Rv, // water + gas case
+        OnePhase_p, // onephase case
     };
 
     BlackOilPrimaryVariables()
@@ -264,10 +265,14 @@ public:
         EnergyModule::assignPrimaryVars(*this, fluidState);
 
         // determine the meaning of the primary variables
-        if ((gasPresent && oilPresent) || onlyWater)
+        if (FluidSystem::numActivePhases() == 1) {
+            primaryVarsMeaning_ = OnePhase_p;
+        }
+        else if ((gasPresent && oilPresent) || (onlyWater && FluidSystem::phaseIsActive(oilPhaseIdx))) {
             // 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
@@ -287,7 +292,16 @@ public:
         }
 
         // assign the actual primary variables
-        if (primaryVarsMeaning() == Sw_po_Sg) {
+        if (primaryVarsMeaning() == OnePhase_p) {
+            if (waterEnabled) {
+                (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+                (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx));
+            } else {
+                throw std::logic_error("For single-phase runs, only pure water is presently allowed.");
+            }
+            
+        }
+        else if (primaryVarsMeaning() == Sw_po_Sg) {
             if (waterEnabled)
                 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
             (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
@@ -340,6 +354,9 @@ public:
         // 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.
+        if (primaryVarsMeaning() == OnePhase_p){
+            return false;
+        }
         Scalar Sw = 0.0;
         if (waterEnabled)
             Sw = (*this)[Indices::waterSaturationIdx];

opm/models/blackoil/blackoilratevector.hh

@@ -83,17 +83,6 @@ public:
     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
      */

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -104,7 +104,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the primary variable for the foam
     static const int foamConcentrationIdx =
-        enableFoam ? polymerMoleWeightIdx + 1 : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
@@ -169,7 +169,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the continuity equation for the foam component
     static const int contiFoamEqIdx =
-        enableFoam ? contiPolymerMWEqIdx + 1 : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =

opm/models/discretization/common/fvbaselocalresidual.hh

@@ -579,8 +579,10 @@ protected:
 
             // Use the implicit Euler time discretization
             for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                double dt = elemCtx.simulator().timeStepSize();
+                assert(dt > 0);
                 tmp[eqIdx] -= tmp2[eqIdx];
-                tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+                tmp[eqIdx] *= scvVolume / dt;
 
                 residual[dofIdx][eqIdx] += tmp[eqIdx];
             }

opm/models/utils/simulator.hh

@@ -125,7 +125,7 @@ public:
         time_ = 0.0;
         endTime_ = EWOMS_GET_PARAM(TypeTag, Scalar, EndTime);
         timeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize);
-
+        assert(timeStepSize_ > 0);
         const std::string& predetTimeStepFile =
             EWOMS_GET_PARAM(TypeTag, std::string, PredeterminedTimeStepsFile);
         if (!predetTimeStepFile.empty()) {
@@ -400,7 +400,9 @@ public:
      * \param timeStepSize The new value for the time step size \f$\mathrm{[s]}\f$
      */
     void setTimeStepSize(Scalar value)
-    { timeStepSize_ = value; }
+    {
+         timeStepSize_ = value;
+    }
 
     /*!
      * \brief Set the current time step index to a given value.
@@ -760,7 +762,7 @@ public:
                 else
                     // ask the problem to provide the next time step size
                     dt = std::min(maxTimeStepSize(), problem_->nextTimeStepSize());
-
+                assert(dt > 0);
                 setTimeStepSize(dt);
             }
             prePostProcessTimer_.stop();