diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh index 6cae6c4f3..c40d0462c 100644 --- a/opm/models/blackoil/blackoilenergymodules.hh +++ b/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; } diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 15b3b819c..c8454d0fd 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/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) - Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx); - + 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 = - Opm::max(fluidState_.saturation(oilPhaseIdx), - elemCtx.problem().maxOilSaturation(globalSpaceIdx)); + Evaluation SoMax = 0.0; + 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 { - assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv); - + else if (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 paramCache; paramCache.setRegionIndex(pvtRegionIdx); - paramCache.setMaxOilSat(SoMax); + 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; } diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh new file mode 100644 index 000000000..5cc866507 --- /dev/null +++ b/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 . + + 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 + +namespace Ewoms { + +/*! + * \ingroup BlackOilModel + * + * \brief The primary variable and equation indices for the black-oil model. + */ +template +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 diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index b16683278..9acc4a467 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/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]; diff --git a/opm/models/blackoil/blackoilratevector.hh b/opm/models/blackoil/blackoilratevector.hh index 14db92874..5e14e2ed7 100644 --- a/opm/models/blackoil/blackoilratevector.hh +++ b/opm/models/blackoil/blackoilratevector.hh @@ -83,17 +83,6 @@ public: BlackOilRateVector(Scalar value) : ParentType(Toolbox::createConstant(value)) {} - template - BlackOilRateVector(const typename std::enable_if::value, Evaluation>::type& value) : ParentType(value) - {} - - /*! - * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const - * ImmiscibleRateVector& ) - */ - BlackOilRateVector(const BlackOilRateVector& value) : ParentType(value) - {} - /*! * \copydoc ImmiscibleRateVector::setMassRate */ diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index 129517a58..0346c6c88 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/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 = @@ -167,9 +167,9 @@ struct BlackOilTwoPhaseIndices static const int contiPolymerMWEqIdx = numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; - //! Index of the continuity equation for the foam component + //! 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 = diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index 3f6dfa990..1b302b594 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/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]; } diff --git a/opm/models/utils/simulator.hh b/opm/models/utils/simulator.hh index edd5d8132..a0d013038 100644 --- a/opm/models/utils/simulator.hh +++ b/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();