From 17e9a73f8154d025c87ed11868a58f71ba47ae2a Mon Sep 17 00:00:00 2001 From: hnil Date: Sun, 19 May 2019 20:42:41 +0200 Subject: [PATCH 01/14] tried to enable onephase with water as phase --- .../blackoil/blackoilintensivequantities.hh | 23 ++++++++++++++----- .../blackoil/blackoilprimaryvariables.hh | 17 ++++++++++++-- opm/models/blackoil/blackoilratevector.hh | 10 ++++---- 3 files changed, 37 insertions(+), 13 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 15b3b819c..e1223b4a6 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -128,9 +128,13 @@ 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::px){ + Sw = 1.0; + }else{ + Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx); + } + } Evaluation Sg = 0.0; if (compositionSwitchEnabled) { @@ -257,8 +261,8 @@ public: else fluidState_.setRv(0.0); } - else { - assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv); + 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); @@ -277,6 +281,8 @@ public: } else fluidState_.setRs(0.0); + }else{ + assert(priVars.primaryVarsMeaning() == PrimaryVariables::px); } typename FluidSystem::template ParameterCache paramCache; @@ -338,7 +344,12 @@ 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{ + x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); + } porosity_ *= 1.0 + x + 0.5*x*x; } diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index b16683278..e4a779d42 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 + px , // onephase case }; BlackOilPrimaryVariables() @@ -264,10 +265,13 @@ public: EnergyModule::assignPrimaryVars(*this, fluidState); // determine the meaning of the primary variables - if ((gasPresent && oilPresent) || onlyWater) + if ( FluidSystem::numActivePhases() == 1 ){ + primaryVarsMeaning_ = px; + }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 +291,16 @@ public: } // assign the actual primary variables - if (primaryVarsMeaning() == Sw_po_Sg) { + if (primaryVarsMeaning() == px ) { + if (waterEnabled){ + (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); + (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx)); + }else{ + throw std::logic_error("Only pure ware 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)); diff --git a/opm/models/blackoil/blackoilratevector.hh b/opm/models/blackoil/blackoilratevector.hh index 14db92874..0a99c10c7 100644 --- a/opm/models/blackoil/blackoilratevector.hh +++ b/opm/models/blackoil/blackoilratevector.hh @@ -83,16 +83,16 @@ public: BlackOilRateVector(Scalar value) : ParentType(Toolbox::createConstant(value)) {} - template - BlackOilRateVector(const typename std::enable_if::value, Evaluation>::type& value) : ParentType(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) - {} + // BlackOilRateVector(const BlackOilRateVector& value) : ParentType(value) + // {} /*! * \copydoc ImmiscibleRateVector::setMassRate From 9c0520aeeb1864eff6a0cc401f839e69118b87a7 Mon Sep 17 00:00:00 2001 From: hnil Date: Sun, 19 May 2019 20:46:38 +0200 Subject: [PATCH 02/14] added onephase indices --- .../blackoil/blackoilonephaseindices.hh | 153 ++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 opm/models/blackoil/blackoilonephaseindices.hh diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh new file mode 100644 index 000000000..571725ba3 --- /dev/null +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -0,0 +1,153 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + 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 = (enableCanonicalCompIdx == 0); + static const bool waterEnabled = (enableCanonicalCompIdx == 1); + static const bool gasEnabled = (enableCanonicalCompIdx == 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; + + //! The number of fluid phases + static const int numPhases = 1; + + //! The number of equations + static const int numEq = numPhases + numSolvents + numPolymers + numEnergy; + + ////////////////////////////// + // 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 temperature + static const int temperatureIdx = + enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers : - 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 < 2); + if(gasEnabled) { + return 2; + } else if (waterEnabled) { + return 1; + } else { + assert(oilEnabled); + } + return 0; + } + + //! Index of the continuity equation of the first phase + static const int conti0EqIdx = PVOffset + 0; + // one continuity equation follows + + //! Index of the continuity equation for the first solvent component + static const int contiSolventEqIdx = + enableSolvent ? PVOffset + numPhases : -1000; + + //! Index of the continuity equation for the first polymer component + static const int contiPolymerEqIdx = + enablePolymer ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the continuity equation for the second polymer component (molecular weight) + static const int contiPolymerMWEqIdx = + numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; + + //! Index of the continuity equation for energy + static const int contiEnergyEqIdx = + enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers : -1000; +}; + +} // namespace Ewoms + +#endif From 8572e1d1089e215fdb817787f13905cd1a6ed5db Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 20 May 2019 20:30:40 +0200 Subject: [PATCH 03/14] make it possible with only one phase in genneral? --- opm/models/blackoil/blackoilintensivequantities.hh | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index e1223b4a6..bd59dbec3 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -348,7 +348,11 @@ public: if(FluidSystem::phaseIsActive(oilPhaseIdx)){ x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure); }else{ - x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); + 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; } From efe3313f85d8073e2adf26d856d5b0de700fb9d0 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 21 May 2019 14:32:25 +0200 Subject: [PATCH 04/14] not adapting any variables in onephase --- opm/models/blackoil/blackoilprimaryvariables.hh | 3 +++ 1 file changed, 3 insertions(+) diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index e4a779d42..8ddf27773 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -353,6 +353,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() == px){ + return false; + } Scalar Sw = 0.0; if (waterEnabled) Sw = (*this)[Indices::waterSaturationIdx]; From c1640aaf3e91972761bcdc86c7b451f50e23aeca Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 21 May 2019 21:52:00 +0200 Subject: [PATCH 05/14] avoid setting entalpy in non active phases (assert will else be broken) --- opm/models/blackoil/blackoilenergymodules.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh index 6cae6c4f3..2021b74ef 100644 --- a/opm/models/blackoil/blackoilenergymodules.hh +++ b/opm/models/blackoil/blackoilenergymodules.hh @@ -376,7 +376,7 @@ public: // and the thermal condictivity coefficients for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { - fs.setEnthalpy(phaseIdx, 0.0); + //fs.setEnthalpy(phaseIdx, 0.0); continue; } From d9b981e0591cfc60dcf9801007c6cef995de6a50 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 21 May 2019 22:33:14 +0200 Subject: [PATCH 06/14] made onephase water runs possible --- opm/models/blackoil/blackoilintensivequantities.hh | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index bd59dbec3..1f1ca165f 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -205,9 +205,13 @@ 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; + //const Evaluation& + 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 @@ -287,7 +291,9 @@ public: 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 From 3a56213ed8af83453830229b22bd660f147822d5 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 22 May 2019 10:07:15 +0200 Subject: [PATCH 07/14] Added several checks to ensure timestep>0 --- .../discretization/common/fvbaselocalresidual.hh | 4 +++- opm/models/utils/simulator.hh | 11 ++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index 3f6dfa990..d4efde481 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..f52ac7604 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,12 @@ public: * \param timeStepSize The new value for the time step size \f$\mathrm{[s]}\f$ */ void setTimeStepSize(Scalar value) - { timeStepSize_ = value; } + { + if( not(value>0)){ + std::cout << "Time step set to zero" << std::endl; + } + timeStepSize_ = value; + } /*! * \brief Set the current time step index to a given value. @@ -443,7 +448,7 @@ public: */ bool finished() const { - assert(timeStepSize_ >= 0.0); + assert(timeStepSize_ > 0.0); Scalar eps = std::max(Scalar(std::abs(this->time())), timeStepSize()) *std::numeric_limits::epsilon()*1e3; From ee482e87a6c7933ec3bb30d2e54163973471157f Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 23 May 2019 09:13:45 +0200 Subject: [PATCH 08/14] change according to comments in pullrequest --- opm/models/blackoil/blackoilenergymodules.hh | 1 - .../blackoil/blackoilintensivequantities.hh | 22 ++++++++----------- .../blackoil/blackoilonephaseindices.hh | 10 ++++----- opm/models/utils/simulator.hh | 9 +++----- 4 files changed, 17 insertions(+), 25 deletions(-) diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh index 2021b74ef..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 1f1ca165f..984f5a3b5 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -140,6 +140,7 @@ public: { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) // -> threephase case + assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::px) ) Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { // -> gas-water case @@ -206,7 +207,6 @@ public: asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx); Evaluation SoMax=0; - //const Evaluation& if(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ SoMax = Opm::max(fluidState_.saturation(oilPhaseIdx), @@ -266,11 +266,9 @@ public: fluidState_.setRv(0.0); } 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); - + if (FluidSystem::enableDissolvedGas()) { // the oil phase is not present, but we need to compute its "composition" for // the gravity correction anyway @@ -280,12 +278,12 @@ public: oilPhaseIdx, pvtRegionIdx, SoMax); - + fluidState_.setRs(Opm::min(RsMax, RsSat)); - } - else + } else { fluidState_.setRs(0.0); - }else{ + } + } else { assert(priVars.primaryVarsMeaning() == PrimaryVariables::px); } @@ -353,12 +351,10 @@ public: Evaluation x; if(FluidSystem::phaseIsActive(oilPhaseIdx)){ x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure); + }else if( FluidSystem::phaseIsActive(waterPhaseIdx) ){ + x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); }else{ - if(FluidSystem::phaseIsActive(waterPhaseIdx)){ - x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); - }else{ - x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure); - } + 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 index 571725ba3..8be4f4622 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -37,13 +37,13 @@ namespace Ewoms { * * \brief The primary variable and equation indices for the black-oil model. */ -template +template struct BlackOilOnePhaseIndices { //! Is phase enabled or not - static const bool oilEnabled = (enableCanonicalCompIdx == 0); - static const bool waterEnabled = (enableCanonicalCompIdx == 1); - static const bool gasEnabled = (enableCanonicalCompIdx == 2); + 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; @@ -116,7 +116,7 @@ struct BlackOilOnePhaseIndices static unsigned activeToCanonicalComponentIndex(unsigned compIdx) { // assumes canonical oil = 0, water = 1, gas = 2; - assert(compIdx < 2); + assert(compIdx == 0); if(gasEnabled) { return 2; } else if (waterEnabled) { diff --git a/opm/models/utils/simulator.hh b/opm/models/utils/simulator.hh index f52ac7604..f8a9a16c6 100644 --- a/opm/models/utils/simulator.hh +++ b/opm/models/utils/simulator.hh @@ -401,10 +401,7 @@ public: */ void setTimeStepSize(Scalar value) { - if( not(value>0)){ - std::cout << "Time step set to zero" << std::endl; - } - timeStepSize_ = value; + timeStepSize_ = value; } /*! @@ -448,7 +445,7 @@ public: */ bool finished() const { - assert(timeStepSize_ > 0.0); + assert(timeStepSize_ >= 0.0); Scalar eps = std::max(Scalar(std::abs(this->time())), timeStepSize()) *std::numeric_limits::epsilon()*1e3; @@ -765,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(); From b4bd3d66e4f27df861179769887695bba0fa9dbd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 8 Aug 2019 09:21:35 +0200 Subject: [PATCH 09/14] Add braces where necessary (not just a style change). --- opm/models/blackoil/blackoilintensivequantities.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 984f5a3b5..7eacd742f 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -138,11 +138,11 @@ public: Evaluation Sg = 0.0; if (compositionSwitchEnabled) { - if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) + if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { // -> threephase case - assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::px) ) + assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::px) ); 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; From 544aeea40ba66b55f4160a1b5c559389e5964b8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 8 Aug 2019 10:27:17 +0200 Subject: [PATCH 10/14] Minor modifications requested in review. --- opm/models/blackoil/blackoilintensivequantities.hh | 6 +++--- opm/models/blackoil/blackoilonephaseindices.hh | 3 +-- opm/models/blackoil/blackoilprimaryvariables.hh | 8 ++++---- opm/models/blackoil/blackoilratevector.hh | 11 ----------- opm/models/utils/simulator.hh | 2 +- 5 files changed, 9 insertions(+), 21 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 7eacd742f..50dfa82f6 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -129,7 +129,7 @@ public: // extract the water and the gas saturations for convenience Evaluation Sw = 0.0; if (waterEnabled){ - if(priVars.primaryVarsMeaning() == PrimaryVariables::px){ + if(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p){ Sw = 1.0; }else{ Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx); @@ -140,7 +140,7 @@ public: { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { // -> threephase case - assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::px) ); + assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) ); Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); } else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { // -> gas-water case @@ -284,7 +284,7 @@ public: fluidState_.setRs(0.0); } } else { - assert(priVars.primaryVarsMeaning() == PrimaryVariables::px); + assert(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p); } typename FluidSystem::template ParameterCache paramCache; diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index 8be4f4622..d857f0d3d 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -127,9 +127,8 @@ struct BlackOilOnePhaseIndices return 0; } - //! Index of the continuity equation of the first phase + //! Index of the continuity equation of the first (and only) phase static const int conti0EqIdx = PVOffset + 0; - // one continuity equation follows //! Index of the continuity equation for the first solvent component static const int contiSolventEqIdx = diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index 8ddf27773..b21efa69e 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -112,7 +112,7 @@ public: Sw_po_Sg, // threephase case Sw_po_Rs, // water + oil case Sw_pg_Rv, // water + gas case - px , // onephase case + OnePhase_p, // onephase case }; BlackOilPrimaryVariables() @@ -266,7 +266,7 @@ public: // determine the meaning of the primary variables if ( FluidSystem::numActivePhases() == 1 ){ - primaryVarsMeaning_ = px; + 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.) @@ -291,7 +291,7 @@ public: } // assign the actual primary variables - if (primaryVarsMeaning() == px ) { + if (primaryVarsMeaning() == OnePhase_p) { if (waterEnabled){ (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx)); @@ -353,7 +353,7 @@ 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() == px){ + if (primaryVarsMeaning() == OnePhase_p){ return false; } Scalar Sw = 0.0; diff --git a/opm/models/blackoil/blackoilratevector.hh b/opm/models/blackoil/blackoilratevector.hh index 0a99c10c7..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/utils/simulator.hh b/opm/models/utils/simulator.hh index f8a9a16c6..c1db35ee4 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); + assert(timeStepSize_ > 0); const std::string& predetTimeStepFile = EWOMS_GET_PARAM(TypeTag, std::string, PredeterminedTimeStepsFile); if (!predetTimeStepFile.empty()) { From 8977d64bc6ef28504f470fcea3ba4c2960498175 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 8 Aug 2019 11:38:50 +0200 Subject: [PATCH 11/14] Whitespace fixes and style consistency. --- .../blackoil/blackoilintensivequantities.hh | 27 +++++++++---------- .../blackoil/blackoilprimaryvariables.hh | 11 ++++---- .../common/fvbaselocalresidual.hh | 4 +-- opm/models/utils/simulator.hh | 2 +- 4 files changed, 22 insertions(+), 22 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 50dfa82f6..c8454d0fd 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -128,10 +128,10 @@ public: // extract the water and the gas saturations for convenience Evaluation Sw = 0.0; - if (waterEnabled){ - if(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p){ + if (waterEnabled) { + if (priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) { Sw = 1.0; - }else{ + } else { Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx); } } @@ -140,7 +140,7 @@ public: { if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { // -> threephase case - assert( not(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) ); + assert( priVars.primaryVarsMeaning() != PrimaryVariables::OnePhase_p ); Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); } else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { // -> gas-water case @@ -206,11 +206,10 @@ public: // update the Saturation functions for the blackoil solvent module. asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx); - Evaluation SoMax=0; - if(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ - 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 @@ -268,7 +267,7 @@ public: else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) { const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx); fluidState_.setRv(Rv); - + if (FluidSystem::enableDissolvedGas()) { // the oil phase is not present, but we need to compute its "composition" for // the gravity correction anyway @@ -278,7 +277,7 @@ public: oilPhaseIdx, pvtRegionIdx, SoMax); - + fluidState_.setRs(Opm::min(RsMax, RsSat)); } else { fluidState_.setRs(0.0); @@ -349,11 +348,11 @@ public: if (rockCompressibility > 0.0) { Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx); Evaluation x; - if(FluidSystem::phaseIsActive(oilPhaseIdx)){ + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure); - }else if( FluidSystem::phaseIsActive(waterPhaseIdx) ){ + } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){ x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure); - }else{ + } else { x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure); } porosity_ *= 1.0 + x + 0.5*x*x; diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index b21efa69e..6807d76c7 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -265,13 +265,14 @@ public: EnergyModule::assignPrimaryVars(*this, fluidState); // determine the meaning of the primary variables - if ( FluidSystem::numActivePhases() == 1 ){ + if (FluidSystem::numActivePhases() == 1) { primaryVarsMeaning_ = OnePhase_p; - }else if ((gasPresent && oilPresent) || (onlyWater && FluidSystem::phaseIsActive(oilPhaseIdx)) ){ + } + 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 @@ -292,10 +293,10 @@ public: // assign the actual primary variables if (primaryVarsMeaning() == OnePhase_p) { - if (waterEnabled){ + if (waterEnabled) { (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx)); - }else{ + } else { throw std::logic_error("Only pure ware is presently allowed"); } diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index d4efde481..1b302b594 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/opm/models/discretization/common/fvbaselocalresidual.hh @@ -580,9 +580,9 @@ protected: // Use the implicit Euler time discretization for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { double dt = elemCtx.simulator().timeStepSize(); - assert(dt>0); + assert(dt > 0); tmp[eqIdx] -= tmp2[eqIdx]; - tmp[eqIdx] *= scvVolume / dt; + tmp[eqIdx] *= scvVolume / dt; residual[dofIdx][eqIdx] += tmp[eqIdx]; } diff --git a/opm/models/utils/simulator.hh b/opm/models/utils/simulator.hh index c1db35ee4..a0d013038 100644 --- a/opm/models/utils/simulator.hh +++ b/opm/models/utils/simulator.hh @@ -762,7 +762,7 @@ public: else // ask the problem to provide the next time step size dt = std::min(maxTimeStepSize(), problem_->nextTimeStepSize()); - assert(dt>0); + assert(dt > 0); setTimeStepSize(dt); } prePostProcessTimer_.stop(); From 7fe1ad5c695fe38c5187b3126b48423298d0b9d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 09:54:43 +0200 Subject: [PATCH 12/14] Fixed error message and reinstated dummy enthalpy update. --- opm/models/blackoil/blackoilenergymodules.hh | 1 + opm/models/blackoil/blackoilprimaryvariables.hh | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh index c40d0462c..6cae6c4f3 100644 --- a/opm/models/blackoil/blackoilenergymodules.hh +++ b/opm/models/blackoil/blackoilenergymodules.hh @@ -376,6 +376,7 @@ 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/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index 6807d76c7..9acc4a467 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -297,7 +297,7 @@ public: (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx)); } else { - throw std::logic_error("Only pure ware is presently allowed"); + throw std::logic_error("For single-phase runs, only pure water is presently allowed."); } } From bc0eaca9432f12985958848f07c8eec9b95820e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 15:55:04 +0200 Subject: [PATCH 13/14] Remove buggy call to setEnthalpy(). --- opm/models/blackoil/blackoilenergymodules.hh | 1 - 1 file changed, 1 deletion(-) 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; } From b9509cdf78663a5d1a9e254540adefa07fc06050 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 16:08:04 +0200 Subject: [PATCH 14/14] Add or fix foam-related indices for one- and two-phase cases. --- .../blackoil/blackoilonephaseindices.hh | 21 ++++++++++++++----- .../blackoil/blackoiltwophaseindices.hh | 6 +++--- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index d857f0d3d..5cc866507 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -26,7 +26,7 @@ * \copydoc Ewoms::BlackOilTwoPhaseIndices */ #ifndef EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH -#define EWOMS_BLACK_OIL_ONe_PHASE_INDICES_HH +#define EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH #include @@ -37,7 +37,7 @@ namespace Ewoms { * * \brief The primary variable and equation indices for the black-oil model. */ -template +template struct BlackOilOnePhaseIndices { //! Is phase enabled or not @@ -63,11 +63,14 @@ struct BlackOilOnePhaseIndices //! 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; + static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam; ////////////////////////////// // Primary variable indices @@ -99,9 +102,13 @@ struct BlackOilOnePhaseIndices 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 : - 1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam: - 1000; ////////////////////// // Equation indices @@ -142,9 +149,13 @@ struct BlackOilOnePhaseIndices 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 : -1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000; }; } // namespace Ewoms 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 =