Accommodate for two or three phases in compositional simulations.

This commit is contained in:
Svenn Tveit
2025-01-20 12:52:54 +01:00
parent 2f0443947d
commit f59834d3e6
10 changed files with 75 additions and 29 deletions

View File

@@ -77,6 +77,14 @@ struct NumComp<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr int value = 3; static constexpr int value = 3;
}; };
template <class TypeTag, class MyTypeTag>
struct EnableDummyWater { using type = UndefinedProperty; };
template <class TypeTag>
struct EnableDummyWater<TypeTag, TTag::CO2PTBaseProblem> {
static constexpr bool value = true;
};
// Set the grid type: --->2D // Set the grid type: --->2D
template <class TypeTag> template <class TypeTag>
struct Grid<TypeTag, TTag::CO2PTBaseProblem> { using type = Dune::YaspGrid</*dim=*/2>; }; struct Grid<TypeTag, TTag::CO2PTBaseProblem> { using type = Dune::YaspGrid</*dim=*/2>; };
@@ -105,9 +113,10 @@ struct FluidSystem<TypeTag, TTag::CO2PTBaseProblem>
private: private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>; using Scalar = GetPropType<TypeTag, Properties::Scalar>;
static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>(); static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
static constexpr bool enable_water = getPropValue<TypeTag, Properties::EnableDummyWater>();
public: public:
using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp>; using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp, enable_water>;
}; };
// Set the material Law // Set the material Law

View File

@@ -83,6 +83,8 @@ struct EnableDispersion { using type = UndefinedProperty; };
//! Enable convective mixing? //! Enable convective mixing?
template<class TypeTag, class MyTypeTag> template<class TypeTag, class MyTypeTag>
struct EnableConvectiveMixing { using type = UndefinedProperty; }; struct EnableConvectiveMixing { using type = UndefinedProperty; };
template <class TypeTag, class MyTypeTag>
struct EnableWater { using type = UndefinedProperty; };
} // namespace Opm::Properties } // namespace Opm::Properties

View File

@@ -48,19 +48,20 @@ class FlashIndices
{ {
static constexpr int numComponents = getPropValue<TypeTag, Properties::NumComponents>(); static constexpr int numComponents = getPropValue<TypeTag, Properties::NumComponents>();
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() }; enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
enum { enableWater = getPropValue<TypeTag, Properties::EnableWater>() };
using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>; using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
public: public:
//! All phases active (note: immiscible/"dummy" water phase) //! All phases active (note: immiscible/"dummy" water phase)
static constexpr bool waterEnabled = true; static constexpr bool waterEnabled = enableWater;
static constexpr bool gasEnabled = true; static constexpr bool gasEnabled = true;
static constexpr bool oilEnabled = true; static constexpr bool oilEnabled = true;
//! number of active phases //! number of active phases
static constexpr int numPhases = 3; static constexpr int numPhases = enableWater ? 3 : 2;
//! number of equations/primary variables //! number of equations/primary variables
static const int numEq = numComponents + EnergyIndices::numEq_ + 1; static const int numEq = numComponents + EnergyIndices::numEq_ + (enableWater ? 1 : 0);
// Primary variable indices // Primary variable indices
@@ -71,7 +72,7 @@ public:
static constexpr int z0Idx = pressure0Idx + 1; static constexpr int z0Idx = pressure0Idx + 1;
//! Index of water saturation //! Index of water saturation
static constexpr int water0Idx = z0Idx + numComponents - 1; static constexpr int water0Idx = enableWater ? z0Idx + numComponents - 1 : -1000;
// equation indices // equation indices

View File

@@ -78,6 +78,8 @@ class FlashIntensiveQuantities
enum { pressure0Idx = Indices::pressure0Idx }; enum { pressure0Idx = Indices::pressure0Idx };
enum { water0Idx = Indices::water0Idx}; enum { water0Idx = Indices::water0Idx};
static constexpr bool waterEnabled = Indices::waterEnabled;
using Scalar = GetPropType<TypeTag, Properties::Scalar>; using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>; using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@@ -216,8 +218,12 @@ public:
(R * fluidState_.temperature(FluidSystem::gasPhaseIdx)); (R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
// Update oil/gas saturation; water saturation is a primary variable // Update saturation
Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx); Evaluation Sw = 0.0;
if constexpr (waterEnabled) {
Sw = priVars.makeEvaluation(water0Idx, timeIdx);
fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
}
Evaluation L = fluidState_.L(); Evaluation L = fluidState_.L();
Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0); Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
Evaluation Sg = Opm::max(1 - So - Sw, 0.0); Evaluation Sg = Opm::max(1 - So - Sw, 0.0);
@@ -228,7 +234,6 @@ public:
fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So); fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So);
fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg); fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg);
fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L); fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L);
fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V); fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);
@@ -253,7 +258,10 @@ public:
// set the phase viscosity and density // set the phase viscosity and density
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
paramCache.updatePhase(fluidState_, phaseIdx); if (phaseIdx == static_cast<unsigned int>(FluidSystem::oilPhaseIdx)
|| phaseIdx == static_cast<unsigned int>(FluidSystem::gasPhaseIdx)) {
paramCache.updatePhase(fluidState_, phaseIdx);
}
const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx); const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);

View File

@@ -65,6 +65,8 @@ class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResid
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() }; enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>; using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
static const bool waterEnabled = Indices::waterEnabled;
using Toolbox = Opm::MathToolbox<Evaluation>; using Toolbox = Opm::MathToolbox<Evaluation>;
public: public:
@@ -82,7 +84,7 @@ public:
const auto& fs = intQuants.fluidState(); const auto& fs = intQuants.fluidState();
// compute water storage term // compute water storage term
if (phaseIdx == waterPhaseIdx) { if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
unsigned eqIdx = conti0EqIdx + numComponents; unsigned eqIdx = conti0EqIdx + numComponents;
storage[eqIdx] = storage[eqIdx] =
Toolbox::template decay<LhsEval>(fs.density(phaseIdx)) Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
@@ -160,7 +162,7 @@ public:
up.fluidState().density(phaseIdx) up.fluidState().density(phaseIdx)
* extQuants.volumeFlux(phaseIdx); * extQuants.volumeFlux(phaseIdx);
if (phaseIdx == waterPhaseIdx) { if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
unsigned eqIdx = conti0EqIdx + numComponents; unsigned eqIdx = conti0EqIdx + numComponents;
flux[eqIdx] = tmp; flux[eqIdx] = tmp;
} }
@@ -176,7 +178,7 @@ public:
Toolbox::value(up.fluidState().density(phaseIdx)) Toolbox::value(up.fluidState().density(phaseIdx))
* extQuants.volumeFlux(phaseIdx); * extQuants.volumeFlux(phaseIdx);
if (phaseIdx == waterPhaseIdx) { if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
unsigned eqIdx = conti0EqIdx + numComponents; unsigned eqIdx = conti0EqIdx + numComponents;
flux[eqIdx] = tmp; flux[eqIdx] = tmp;
} }

View File

@@ -135,6 +135,10 @@ template<class TypeTag>
struct EnableEnergy<TypeTag, TTag::FlashModel> struct EnableEnergy<TypeTag, TTag::FlashModel>
{ static constexpr bool value = false; }; { static constexpr bool value = false; };
template<class TypeTag>
struct EnableWater<TypeTag, TTag::MultiPhaseBaseModel>
{ static constexpr int value = GetPropType<TypeTag, Properties::FluidSystem>::waterEnabled; };
} // namespace Opm::Properties } // namespace Opm::Properties
namespace Opm { namespace Opm {

View File

@@ -63,6 +63,8 @@ class FlashNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMeth
enum { z0Idx = Indices::z0Idx }; enum { z0Idx = Indices::z0Idx };
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() }; enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
static constexpr bool waterEnabled = Indices::waterEnabled;
public: public:
/*! /*!
* \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& ) * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& )
@@ -126,10 +128,12 @@ protected:
clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol); clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol);
} }
// limit change in water saturation to 0.2 if constexpr (waterEnabled) {
constexpr Scalar dSwMax = 0.2; // limit change in water saturation to 0.2
if (update[Indices::water0Idx] > dSwMax) { constexpr Scalar dSwMax = 0.2;
nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax; if (update[Indices::water0Idx] > dSwMax) {
nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax;
}
} }
} }
private: private:

View File

@@ -63,6 +63,8 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
enum { pressure0Idx = Indices::pressure0Idx }; enum { pressure0Idx = Indices::pressure0Idx };
enum { water0Idx = Indices::water0Idx }; enum { water0Idx = Indices::water0Idx };
static constexpr bool waterEnabled = Indices::waterEnabled;
// phase indices // phase indices
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@@ -122,7 +124,9 @@ public:
(*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx)); (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx));
// assign water saturation // assign water saturation
(*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx)); if constexpr (waterEnabled) {
(*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
}
} }
/*! /*!
@@ -138,7 +142,9 @@ public:
os << ", z_" << FluidSystem::componentName(compIdx) << " = " os << ", z_" << FluidSystem::componentName(compIdx) << " = "
<< this->operator[](z0Idx + compIdx); << this->operator[](z0Idx + compIdx);
} }
os << ", S_w = " << this->operator[](water0Idx); if constexpr (waterEnabled) {
os << ", S_w = " << this->operator[](water0Idx);
}
os << ")" << std::flush; os << ")" << std::flush;
} }
}; };

View File

@@ -362,7 +362,7 @@ public:
Dune::FieldVector<Scalar, numComponents> z(0.0); Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0; Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (phaseIdx == waterPhaseIdx){ if (Indices::waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)){
continue; continue;
} }
const auto saturation = fs.saturation(phaseIdx); const auto saturation = fs.saturation(phaseIdx);
@@ -494,7 +494,7 @@ protected:
const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx); const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx); const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
if (water_active && Indices::numPhases > 1) if (water_active && Indices::numPhases > 2)
waterSaturationData = fp.get_double("SWAT"); waterSaturationData = fp.get_double("SWAT");
else else
waterSaturationData.resize(numDof); waterSaturationData.resize(numDof);

View File

@@ -1676,17 +1676,27 @@ INSTANTIATE_TYPE(double)
INSTANTIATE_TYPE(float) INSTANTIATE_TYPE(float)
#endif #endif
#define INSTANTIATE_COMP(NUM) \ #define INSTANTIATE_COMP_THREEPHASE(NUM) \
template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM>; \ template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
template class GenericOutputBlackoilModule<FS##NUM<double>>; template class GenericOutputBlackoilModule<FS##NUM<double>>;
INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput
INSTANTIATE_COMP_THREEPHASE(2)
INSTANTIATE_COMP_THREEPHASE(3)
INSTANTIATE_COMP_THREEPHASE(4)
INSTANTIATE_COMP_THREEPHASE(5)
INSTANTIATE_COMP_THREEPHASE(6)
INSTANTIATE_COMP_THREEPHASE(7)
INSTANTIATE_COMP(2) #define INSTANTIATE_COMP_TWOPHASE(NUM) \
INSTANTIATE_COMP(3) template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
INSTANTIATE_COMP(4) template class GenericOutputBlackoilModule<GFS##NUM<double>>;
INSTANTIATE_COMP(5)
INSTANTIATE_COMP(6) INSTANTIATE_COMP_TWOPHASE(2)
INSTANTIATE_COMP(7) INSTANTIATE_COMP_TWOPHASE(3)
INSTANTIATE_COMP_TWOPHASE(4)
INSTANTIATE_COMP_TWOPHASE(5)
INSTANTIATE_COMP_TWOPHASE(6)
INSTANTIATE_COMP_TWOPHASE(7)
} // namespace Opm } // namespace Opm