diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh index 8a9a25083e4..0d2c77d5fa5 100644 --- a/examples/problems/co2ptflashproblem.hh +++ b/examples/problems/co2ptflashproblem.hh @@ -285,6 +285,24 @@ public: return Opm::CompositionalConfig::EOSType::PR; } + template + unsigned pvtRegionIndex(const Context&, unsigned, unsigned) const + { + return 0; + } + + template + Scalar rockCompressibility(const Context&, unsigned, unsigned) const + { + return 0.0; + } + + template + Scalar rockReferencePressure(const Context&, unsigned, unsigned) const + { + return 1.0e5; + } + /*! * \copydoc FvBaseProblem::finishInit */ diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh index 9b97252426d..1f2ed0d409c 100644 --- a/opm/models/ptflash/flashintensivequantities.hh +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -228,6 +228,8 @@ public: // Update phases typename FluidSystem::template ParameterCache paramCache(eos_type); + // the water properties of the cell's PVT region + paramCache.setRegionIndex(problem.pvtRegionIndex(elemCtx, dofIdx, timeIdx)); paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx); paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx); @@ -323,8 +325,15 @@ public: // Compute the remaining quantities ///////////// - // porosity + // porosity, scaled with the pressure by the rock compressibility of the + // cell's region as in the black-oil model porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx); + const Scalar rockCompressibility = problem.rockCompressibility(elemCtx, dofIdx, timeIdx); + if (rockCompressibility > 0.0) { + const Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx); + const Evaluation x = rockCompressibility * (p - rockRefPressure); + porosity_ *= 1.0 + x + 0.5 * x * x; + } Valgrind::CheckDefined(porosity_); // intrinsic permeability diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index f8eb41b4daa..40f49afeb8b 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -437,6 +437,7 @@ class FlowProblemComp : public FlowProblem getEosType(), vanguard.cellCenterDepths(), eqlnum, + this->pvtnum_, vanguard.gridView().comm(), this->gravity()[dimWorld - 1], this->numPressurePointsEquil(), diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp index d5a250e3262..27b184cb159 100644 --- a/opm/simulators/flow/equil/InitStateEquilComp.hpp +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -62,6 +62,7 @@ #include #include #include +#include #include #include #include @@ -148,9 +149,11 @@ class WaterDensityODE WaterDensityODE(const TabulatedFunction& tempVdTable, const CompositionalConfig::EOSType eosType, + const unsigned pvtRegion, const Scalar normGrav) : tempVdTable_(tempVdTable) , eosType_(eosType) + , pvtRegion_(pvtRegion) , g_(normGrav) {} @@ -162,6 +165,7 @@ class WaterDensityODE fs.setPressure(FluidSystem::waterPhaseIdx, press); typename FluidSystem::template ParameterCache paramCache(eosType_); + paramCache.setRegionIndex(pvtRegion_); return FluidSystem::density(fs, paramCache, FluidSystem::waterPhaseIdx) * g_; } @@ -169,6 +173,7 @@ class WaterDensityODE private: const TabulatedFunction& tempVdTable_; CompositionalConfig::EOSType eosType_; + unsigned pvtRegion_; Scalar g_; }; @@ -201,11 +206,18 @@ class WaterDensityODE * * The COMPVD phase column selects the EOS root. If its rows name both phases, * the vapour rows define the gas zone above the gas-oil contact and the liquid - * rows define the liquid zone below it. + * rows define the liquid zone below it. For type 1 a contact outside the gap + * between the last vapour row and the first liquid row moves to the nearest + * edge of that gap. * * Only cell-centre initialization is supported (EQUIL item 9 = 0). * Gas-oil contact capillary pressure must be zero: the downstream flash uses * a single pressure for all phases. + * + * Each PVT region among the cells of an equilibration region has a water column + * of its own, integrated with its water, and each cell takes the column of its + * PVT region. A datum below the water-oil contact gives the pressure of all of + * these waters, so there they must reach the contact at the same pressure. */ template class InitialStateComputer @@ -219,6 +231,8 @@ class InitialStateComputer /// \param[in] eosType Equation of state used by the fluid system. /// \param[in] cellCenterDepth Depth of each cell centre. /// \param[in] eqlnum Zero-based equilibration region of each cell. + /// \param[in] pvtnum Zero-based PVT region of each cell, or empty + /// for a single PVT region. /// \param[in] comm Communicator for parallel runs. /// \param[in] gravity Norm of the gravity vector. /// \param[in] numSamplePoints Sample points in each pressure integration. @@ -230,6 +244,7 @@ class InitialStateComputer const CompositionalConfig::EOSType eosType, const std::vector& cellCenterDepth, const std::vector& eqlnum, + const std::vector& pvtnum, const Parallel::Communication& comm, const Scalar gravity, const int numSamplePoints, @@ -248,6 +263,7 @@ class InitialStateComputer "versus depth from the ZMFVD or the COMPVD keyword."); } + const auto numPvtRegions = tables.getTabdims().getNumPVTTables(); OPM_BEGIN_PARALLEL_TRY_CATCH(); if (eqlnum.size() != cellCenterDepth.size()) { OPM_THROW(std::runtime_error, @@ -263,6 +279,20 @@ class InitialStateComputer cell, region + 1, records.size())); } } + if (!pvtnum.empty() && (pvtnum.size() != cellCenterDepth.size())) { + OPM_THROW(std::runtime_error, + fmt::format("PVTNUM contains {} entries for {} cell depths.", + pvtnum.size(), cellCenterDepth.size())); + } + for (std::size_t cell = 0; cell < pvtnum.size(); ++cell) { + const auto region = pvtnum[cell]; + if (region < 0 || std::cmp_greater_equal(region, numPvtRegions)) { + OPM_THROW(std::runtime_error, + fmt::format("Cell {} has PVTNUM {} outside the {} " + "PVT regions.", + cell, region + 1, numPvtRegions)); + } + } // The endpoint vectors are optional, but a non-empty one is indexed for // every cell. for (const auto& [name, limits] : {std::pair{"connate water", std::cref(connateWater)}, @@ -276,18 +306,22 @@ class InitialStateComputer } OPM_END_PARALLEL_TRY_CATCH("Invalid equilibration input: ", comm); + const auto pvtRegions = + pvtRegionsOfEachRegion(eqlnum, pvtnum, records.size(), numPvtRegions, comm); std::vector regions; regions.reserve(records.size()); for (std::size_t r = 0; r < records.size(); ++r) { - regions.push_back(setupRegion(records.getRecord(r), tables, cellCenterDepth, - eqlnum, comm, gravity, numSamplePoints, r)); + regions.push_back(setupRegion(records.getRecord(r), tables, pvtRegions[r], + cellCenterDepth, eqlnum, comm, gravity, + numSamplePoints, r)); } fluidStates_.resize(cellCenterDepth.size()); referencePressures_.resize(cellCenterDepth.size()); for (std::size_t cell = 0; cell < cellCenterDepth.size(); ++cell) { - referencePressures_[cell] = - assignCell(fluidStates_[cell], regions[eqlnum[cell]], cellCenterDepth[cell], cell); + const auto pvtRegion = pvtnum.empty() ? 0u : static_cast(pvtnum[cell]); + referencePressures_[cell] = assignCell(fluidStates_[cell], regions[eqlnum[cell]], + pvtRegion, cellCenterDepth[cell], cell); } } @@ -325,6 +359,9 @@ class InitialStateComputer /// The equilibrated vertical distributions within one region. struct Region { int initType{1}; // EQUIL item 10 + /// Zero-based PVT regions of the region's cells on all processes. + std::vector pvtRegions; + /// Gas-oil contact. A two-zone COMPVD table may move it, see zoneBoundary(). Scalar zgoc{}; /// Name of the selected composition keyword, used in diagnostics. std::string_view compositionKeyword; @@ -335,6 +372,10 @@ class InitialStateComputer /// COMPVD naming both phases describes a gas zone over a liquid one, /// each with its own composition and its own hydrostatic column. bool twoZone{false}; + /// Depth of the last vapour row of such a table. With EQUIL item 10 = 1 + /// the rows decide the phase, so a cell at this depth stays in the gas + /// zone even when the contact lies on it. + Scalar lastVaporDepth{}; /// Gas-zone composition for a COMPVD table that names both phases. std::vector vaporVdTable; /// Equilibrium vapour at the contact for type 3 without gas-zone rows. @@ -351,7 +392,8 @@ class InitialStateComputer /// anchored at the water-oil contact instead. The fluid there is the one /// just above the contact, which the EOS root has to follow. bool anchoredAtWaterContact{false}; - std::optional waterPressure; + /// The water column of each PVT region in pvtRegions. + std::map waterPressure; }; /// Each pressure column must reach the water-oil contact before another @@ -455,13 +497,17 @@ class InitialStateComputer return rows; } - /// Verify that vapour rows are at or above the gas-oil contact and liquid - /// rows are at or below it. - static void checkZonesStraddleContact(const CompvdTable& compvd, - const std::vector& vaporRows, - const std::vector& liquidRows, - const Scalar zgoc, - const std::size_t regionIdx) + /// Depth at which the gas zone of a COMPVD table naming both phases meets + /// its liquid zone: the gas-oil contact when it lies between the last + /// vapour row and the first liquid row. With EQUIL item 10 = 1 the rows + /// give the composition versus depth, so a contact outside that gap moves + /// to its nearest edge and every row keeps the phase it names. Type 3 uses + /// the contact as its reference depth, so there the rows must agree with it. + static Scalar zoneBoundary(const CompvdTable& compvd, + const std::vector& vaporRows, + const std::vector& liquidRows, + const Region& reg, + const std::size_t regionIdx) { if (vaporRows.empty() || liquidRows.empty()) { OPM_THROW(std::runtime_error, @@ -470,14 +516,39 @@ class InitialStateComputer } const auto& depth = compvd.getDepthColumn(); - if ((depth[vaporRows.back()] > zgoc) || (depth[liquidRows.front()] < zgoc)) { + const Scalar lastVapor = depth[vaporRows.back()]; + const Scalar firstLiquid = depth[liquidRows.front()]; + if (lastVapor > firstLiquid) { + OPM_THROW(std::runtime_error, + fmt::format("The COMPVD table of region {} has a vapour row at {} m " + "below its liquid row at {} m. All vapour rows must lie " + "above the liquid rows.", + regionIdx + 1, lastVapor, firstLiquid)); + } + + const Scalar boundary = std::clamp(reg.zgoc, lastVapor, firstLiquid); + if (boundary == reg.zgoc) { + return boundary; + } + + if (reg.initType != 1) { OPM_THROW(std::runtime_error, fmt::format("The COMPVD table of region {} puts its vapour rows down to " "{} m and its liquid rows from {} m, which do not meet at " - "the gas-oil contact at {} m.", - regionIdx + 1, depth[vaporRows.back()], - depth[liquidRows.front()], zgoc)); + "the gas-oil contact at {} m. EQUIL item 10 = {} takes the " + "contact as the reference depth, so the rows must agree " + "with it.", + regionIdx + 1, lastVapor, firstLiquid, reg.zgoc, + reg.initType)); } + + OpmLog::warning(fmt::format("Equilibration region {}: the gas-oil contact at {} m lies " + "outside the COMPVD phase change between the last vapour " + "row at {} m and the first liquid row at {} m. With EQUIL " + "item 10 = 1 the rows give the composition versus depth, " + "so the gas zone ends at {} m instead.", + regionIdx + 1, reg.zgoc, lastVapor, firstLiquid, boundary)); + return boundary; } /// The COMPVD rows carrying \p phase. @@ -532,8 +603,39 @@ class InitialStateComputer }); } + /// The zero-based PVT regions of the cells of each equilibration region, on + /// all processes. A region without cells has none. + static std::vector> + pvtRegionsOfEachRegion(const std::vector& eqlnum, + const std::vector& pvtnum, + const std::size_t numRegions, + const std::size_t numPvtRegions, + const Parallel::Communication& comm) + { + // Whether equilibration region r holds a cell of PVT region p, at + // r * numPvtRegions + p. + std::vector present(numRegions * numPvtRegions, 0); + for (std::size_t cell = 0; cell < eqlnum.size(); ++cell) { + const auto pvtRegion = pvtnum.empty() ? 0 : pvtnum[cell]; + present[static_cast(eqlnum[cell]) * numPvtRegions + + static_cast(pvtRegion)] = 1; + } + comm.max(present.data(), static_cast(present.size())); + + std::vector> pvtRegions(numRegions); + for (std::size_t r = 0; r < numRegions; ++r) { + for (std::size_t p = 0; p < numPvtRegions; ++p) { + if (present[r * numPvtRegions + p] != 0) { + pvtRegions[r].push_back(static_cast(p)); + } + } + } + return pvtRegions; + } + Region setupRegion(const EquilRecord& record, const TableManager& tables, + const std::vector& pvtRegions, const std::vector& cellCenterDepth, const std::vector& eqlnum, const Parallel::Communication& comm, @@ -543,6 +645,7 @@ class InitialStateComputer { Region reg; + reg.pvtRegions = pvtRegions; reg.initType = record.compositionalInitType(); if (reg.initType != 1 && reg.initType != 3) { OPM_THROW(std::runtime_error, @@ -645,9 +748,9 @@ class InitialStateComputer // the liquid rows the one below the contact. const auto vaporRows = rowsOfPhase(compvd, CompvdTable::Phase::Vapor); const auto liquidRows = rowsOfPhase(compvd, CompvdTable::Phase::Liquid); - checkZonesStraddleContact(compvd, vaporRows, liquidRows, - record.gasOilContactDepth(), regionIdx); + reg.zgoc = zoneBoundary(compvd, vaporRows, liquidRows, reg, regionIdx); reg.twoZone = true; + reg.lastVaporDepth = compvd.getDepthColumn()[vaporRows.back()]; setupComposition(reg.vaporVdTable, compvd, vaporRows); setupComposition(reg.compositionVdTable, compvd, liquidRows); } @@ -715,7 +818,7 @@ class InitialStateComputer integrateWaterPressure(reg, span, gravity, numSamplePoints, record.datumDepth(), record.datumDepthPressure()); hcDatum = reg.zwoc; - hcPressure = reg.waterPressure->value(reg.zwoc) + hcPressure = waterContactPressure(reg, record.datumDepth(), regionIdx) + record.waterOilContactCapillaryPressure(); OpmLog::info(fmt::format("Equilibration region {}: the datum at {} m lies below the " "water-oil contact at {} m, so it gives the water pressure; " @@ -774,7 +877,11 @@ class InitialStateComputer if (!reg.oilPressure.has_value() && !reg.gasPressure.has_value()) { return; } - const auto& hcPressure = reg.oilPressure.has_value() ? reg.oilPressure : reg.gasPressure; + // The water meets the hydrocarbon zone just above its contact: the gas + // zone when the gas-oil contact lies at or below the water-oil contact. + const bool gasAtContact = reg.gasPressure.has_value() + && (!reg.oilPressure.has_value() || (reg.zwoc <= reg.zgoc)); + const auto& hcPressure = gasAtContact ? reg.gasPressure : reg.oilPressure; const Scalar pcow = record.waterOilContactCapillaryPressure(); const Scalar pContact = hcPressure->value(reg.zwoc) - pcow; @@ -784,8 +891,8 @@ class InitialStateComputer "is at {} m.", regionIdx + 1, reg.zwoc)); } - /// Integrates the water pressure over \p span from \p depth, where it is - /// \p pressure. + /// Integrates the water pressure of each PVT region of \p reg over \p span + /// from \p depth, where it is \p pressure. void integrateWaterPressure(Region& reg, const std::array& span, const Scalar gravity, @@ -793,10 +900,43 @@ class InitialStateComputer const Scalar depth, const Scalar pressure) const { - const WaterODE ode(reg.tempVdTable, eosType_, gravity); - reg.waterPressure.emplace(ode, - typename WaterPressFunc::InitCond{depth, pressure}, - numSamplePoints, waterContactSpan(reg, span)); + for (const auto pvtRegion : reg.pvtRegions) { + const WaterODE ode(reg.tempVdTable, eosType_, pvtRegion, gravity); + reg.waterPressure.try_emplace(pvtRegion, ode, + typename WaterPressFunc::InitCond{depth, pressure}, + numSamplePoints, waterContactSpan(reg, span)); + } + } + + /// The water pressure at the water-oil contact when the datum at \p datum, + /// below the contact, gives it. The datum states the pressure of the water + /// of every PVT region in the region, so their columns must reach the + /// contact at one pressure for the hydrocarbon to be anchored there. + static Scalar waterContactPressure(const Region& reg, + const Scalar datum, + const std::size_t regionIdx) + { + // Waters that differ physically part by more than this round-off + // tolerance between the datum and the contact. + constexpr Scalar samePressure{1.0e-10}; + + const auto& [firstPvtRegion, firstColumn] = *reg.waterPressure.begin(); + const Scalar pressure = firstColumn.value(reg.zwoc); + for (const auto& [pvtRegion, column] : reg.waterPressure) { + const Scalar other = column.value(reg.zwoc); + if (std::abs(other - pressure) > samePressure * pressure) { + OPM_THROW(std::runtime_error, + fmt::format("Equilibration region {} places the datum at {} m, below " + "the water-oil contact at {} m, but the water of its PVT " + "regions {} and {} reaches the contact at {:.6g} and " + "{:.6g} bar. Put the datum in the hydrocarbon column.", + regionIdx + 1, datum, reg.zwoc, + firstPvtRegion + 1, pvtRegion + 1, + unit::convert::to(pressure, unit::barsa), + unit::convert::to(other, unit::barsa))); + } + } + return pressure; } /// A per-cell saturation endpoint, or \p fallback when the caller supplied @@ -894,10 +1034,11 @@ class InitialStateComputer const std::array datumSpan{std::min(pressureSpan[0], reg.zgoc), std::max(pressureSpan[1], reg.zgoc)}; if ((reg.zgoc < span[0]) || (reg.zgoc > span[1])) { - OpmLog::warning(fmt::format("Equilibration region {}: the gas-oil contact at {} m " - "lies outside the cells of the region, so the COMPVD " + OpmLog::warning(fmt::format("Equilibration region {}: the gas zone ends at {} m, " + "outside the cells of the region, so the COMPVD " "rows of one phase describe no cell.", - regionIdx + 1, reg.zgoc)); + regionIdx + 1, + reg.zgoc)); } if (datum < reg.zgoc) { @@ -919,9 +1060,11 @@ class InitialStateComputer numSamplePoints, pressureSpan); } - OpmLog::info(fmt::format("Equilibration region {}: COMPVD gives a gas zone above the " - "contact at {} m and a liquid one below it " - "(EQUIL item 10 = 1).", regionIdx + 1, reg.zgoc)); + OpmLog::info(fmt::format("Equilibration region {}: COMPVD gives a gas zone above " + "{} m and a liquid one below it " + "(EQUIL item 10 = 1).", + regionIdx + 1, + reg.zgoc)); } /// EQUIL item 10 type 3: the selected table supplies the liquid composition, @@ -1009,10 +1152,22 @@ class InitialStateComputer numSamplePoints, waterContactSpan(reg, span)); } - Scalar assignCell(FluidState& fs, const Region& reg, const Scalar depth, - const std::size_t cell) const + /// Whether \p depth lies in the gas zone of a region that has one. A depth + /// on the gas-oil contact belongs to the liquid zone, except the last vapour + /// row of a two-zone COMPVD table with EQUIL item 10 = 1, whose rows decide + /// the phase. + static bool isInGasZone(const Region& reg, const Scalar depth) + { + if (reg.twoZone && (reg.initType == 1) && (depth <= reg.lastVaporDepth)) { + return true; + } + return ((reg.initType == 3) || reg.twoZone) && (depth < reg.zgoc); + } + + Scalar assignCell(FluidState& fs, const Region& reg, const unsigned pvtRegion, + const Scalar depth, const std::size_t cell) const { - const bool inGasZone = ((reg.initType == 3) || reg.twoZone) && (depth < reg.zgoc); + const bool inGasZone = isInGasZone(reg, depth); const CompVec z = [®, depth, inGasZone]() { if (!inGasZone) { @@ -1033,8 +1188,12 @@ class InitialStateComputer // saturation leaves some residual hydrocarbon. Keep this separate from // the phase pressures so equilibration retains the capillary offset. const Scalar hydrocarbonPressure = pressFunc->value(depth); - const bool inWaterZone = (depth > reg.zwoc) && reg.waterPressure.has_value(); - const Scalar press = inWaterZone ? reg.waterPressure->value(depth) + // The water column of the cell's own PVT region. + const auto water = reg.waterPressure.find(pvtRegion); + const WaterPressFunc* waterPressure = + (water != reg.waterPressure.end()) ? &water->second : nullptr; + const bool inWaterZone = (depth > reg.zwoc) && (waterPressure != nullptr); + const Scalar press = inWaterZone ? waterPressure->value(depth) : hydrocarbonPressure; fs.setTemperature(Details::evalDepthTable(reg.tempVdTable, depth)); @@ -1055,8 +1214,8 @@ class InitialStateComputer sWat = (depth > reg.zwoc) ? waterLimit(maxWater_, cell, Scalar{1}) : waterLimit(connateWater_, cell, Scalar{0}); fs.setSaturation(FluidSystem::waterPhaseIdx, sWat); - if (reg.waterPressure.has_value()) { - fs.setPressure(FluidSystem::waterPhaseIdx, reg.waterPressure->value(depth)); + if (waterPressure != nullptr) { + fs.setPressure(FluidSystem::waterPhaseIdx, waterPressure->value(depth)); } } diff --git a/regressionTests.cmake b/regressionTests.cmake index 29f184d7f4d..32fbb0c686b 100644 --- a/regressionTests.cmake +++ b/regressionTests.cmake @@ -218,6 +218,23 @@ add_test_compareECLFiles( compositional/equilibration ) +add_test_compareECLFiles( + CASENAME + equil_1d_compvd_water_gascap_contact_mismatch + FILENAME + EQUIL_1D_COMPVD_WATER_GASCAP_CONTACT_MISMATCH + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/equilibration +) + add_test_compareECLFiles( CASENAME equil_1d_compvd_oil @@ -235,6 +252,74 @@ add_test_compareECLFiles( compositional/equilibration ) +add_test_compareECLFiles( + CASENAME + pvtnum_water + FILENAME + PVTNUM_WATER + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + +add_test_compareECLFiles( + CASENAME + pvtnum_equil + FILENAME + PVTNUM_EQUIL + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + +add_test_compareECLFiles( + CASENAME + pvtnum_equil_shared + FILENAME + PVTNUM_EQUIL_SHARED + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + +add_test_compareECLFiles( + CASENAME + pvtnum_rock + FILENAME + PVTNUM_ROCK + SIMULATOR + flow_comp + REFERENCE_SIMULATOR + flow_comp + ABS_TOL + ${abs_tol} + REL_TOL + ${rel_tol} + DIR + compositional/pvt_regions +) + add_test_compareECLFiles( CASENAME spe12 diff --git a/tests/test_compequil.cpp b/tests/test_compequil.cpp index b8d3ac6e788..d029e1a3193 100644 --- a/tests/test_compequil.cpp +++ b/tests/test_compequil.cpp @@ -45,6 +45,9 @@ #include #include #include +#include +#include +#include #include namespace { @@ -131,14 +134,17 @@ struct BasicEquilFixture /// \param connateWater, maxWater Per-cell water saturation endpoints, as /// the simulator reads them from the scaled saturation functions. + /// \param pvtnum Zero-based PVT region of each cell, empty for one region. Computer compute(const std::vector& eqlnum, const std::vector& connateWater = {}, - const std::vector& maxWater = {}) const + const std::vector& maxWater = {}, + const std::vector& pvtnum = {}) const { return Computer(eclState, eclState.compositionalConfig().eosType(0), {depths.begin(), depths.end()}, eqlnum, + pvtnum, Opm::Parallel::Communication{}, gravity, /*numSamplePoints=*/100, @@ -188,6 +194,66 @@ std::string waterDeckString(const std::string& equil, "SOF3\n 0.00 0.0 0.0\n 0.99 1.0 1.0 /\n"); } +// A deck from deckString() declaring two PVT regions, both with its rock. +std::string withTwoPvtRegions(std::string deck) +{ + for (const auto& [from, to] : {std::pair{std::string{"TABDIMS\n/\n"}, + std::string{"TABDIMS\n1* 2 /\n"}}, + std::pair{std::string{"ROCK\n68.9476 0 /\n"}, + std::string{"ROCK\n68.9476 0 /\n68.9476 0 /\n"}}}) { + const auto pos = deck.find(from); + BOOST_REQUIRE(pos != std::string::npos); + deck.replace(pos, from.size(), to); + } + return deck; +} + +// A water zone at the base of each block of ten cells. +constexpr auto twoBlockEquil = "EQUIL\n" + " 2010 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n"; + +// The fresh water of PVT region 1 and the brine of PVT region 2. +constexpr auto freshWaterAndBrine = "PVTW\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + " 75.0 1.01 2.5E-5 0.9 0.0 /\n" + "DENSITY\n" + " 800.0 1000.0 1.0 /\n" + " 800.0 1150.0 1.0 /\n"; + +// Two equilibration regions and two PVT regions over the column. The arguments +// give the EQLNUM and PVTNUM records, the two EQUIL records, and the PVTW and +// DENSITY of the two PVT regions. +std::string twoWaterRegionsDeck(const std::string& eqlnum, const std::string& pvtnum, + const std::string& equil = twoBlockEquil, + const std::string& water = freshWaterAndBrine) +{ + return withTwoPvtRegions(deckString(equil, + "EQLDIMS\n2 /\n", + "REGIONS\n" + "EQLNUM\n" + eqlnum + " /\n" + "PVTNUM\n" + pvtnum + " /\n", + "ZMFVD\n" + " 2000 0 0.7 0.3\n" + " 2100 0 0.3 0.7 /\n" + " 2000 0 0.7 0.3\n" + " 2100 0 0.3 0.7 /\n", + "RTEMP\n100\n/\n", + "OIL\nGAS\nWATER\n", + "SWFN\n 0.01 0.0 0.0\n 1.00 1.0 0.0 /\n" + "SGFN\n 0.00 0.0 0.0\n 0.99 1.0 0.0 /\n" + "SOF3\n 0.00 0.0 0.0\n 0.99 1.0 1.0 /\n" + + water)); +} + +// The water density of a PVT region above at a pressure in bar. +Scalar waterDensity(const Scalar surfaceDensity, const Scalar bw, const Scalar cw, + const Scalar pressure) +{ + const Scalar x = cw * (pressure - 75.0); + return surfaceDensity * (1.0 + x + 0.5 * x * x) / bw; +} + } // Anonymous namespace BOOST_AUTO_TEST_CASE(Type1LiquidRootPressureIntegration) @@ -848,6 +914,111 @@ BOOST_AUTO_TEST_CASE(CompvdTwoZoneClampsCompositionOutsideRows) } } +BOOST_AUTO_TEST_CASE(CompvdTwoZoneContactOutsideRowGap) +{ + // With EQUIL item 10 = 1 the COMPVD rows give the composition versus depth, + // so a gas-oil contact outside the gap between the last vapour row (2049 m) + // and the first liquid row (2051 m) moves to the nearest edge of the gap. + const auto states = [](const std::string& goc) { + const EquilFixture fix(deckString( + "EQUIL\n 2072.5 200 2300 0 " + goc + " 0 /\n", "EQLDIMS\n/\n", "", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2049 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + return fix.compute(std::vector(20, 0)).fluidStates(); + }; + + // The cell at 2057.5 m lies above a contact at 2060 m but below the first + // liquid row, and the one at 2042.5 m below a contact at 2030 m but above + // the last vapour row. + for (const auto& [goc, edge, cell, phaseIdx] : + {std::tuple{"2060", "2051", std::size_t{11}, FluidSystem::oilPhaseIdx}, + std::tuple{"2030", "2049", std::size_t{8}, FluidSystem::gasPhaseIdx}}) { + BOOST_TEST_CONTEXT("Contact at " << goc << " m") { + const auto moved = states(goc); + BOOST_CHECK_CLOSE(moved[cell].saturation(phaseIdx), 1.0, 1e-10); + + const auto expected = states(edge); + BOOST_REQUIRE_EQUAL(moved.size(), expected.size()); + for (std::size_t c = 0; c < moved.size(); ++c) { + BOOST_TEST_CONTEXT("Cell " << c) { + for (const auto p : {FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}) { + BOOST_CHECK_EQUAL(moved[c].pressure(p), expected[c].pressure(p)); + BOOST_CHECK_EQUAL(moved[c].saturation(p), expected[c].saturation(p)); + } + for (int comp = 0; comp < 3; ++comp) { + BOOST_CHECK_EQUAL(moved[c].moleFraction(comp), + expected[c].moleFraction(comp)); + } + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(CompvdTwoZoneKeepsTheLastVapourRowGas) +{ + // A cell on the gas-oil contact belongs to the liquid zone, but with EQUIL + // item 10 = 1 the rows decide the phase. The cell at 2032.5 m sits on the + // last vapour row, both when the contact lies on that row and when it moves + // there from above. + for (const std::string goc : {"2032.5", "2020"}) { + BOOST_TEST_CONTEXT("Contact at " << goc << " m") { + const EquilFixture fix(deckString( + "EQUIL\n 2072.5 200 2300 0 " + goc + " 0 /\n", "EQLDIMS\n/\n", "", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2032.5 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + const auto states = fix.compute(std::vector(20, 0)).fluidStates(); + BOOST_REQUIRE_EQUAL(fix.depths[6], 2032.5); + BOOST_CHECK_CLOSE(states[6].saturation(FluidSystem::gasPhaseIdx), 1.0, 1e-10); + BOOST_CHECK_SMALL(states[6].moleFraction(1) - 0.95, 1e-10); + + // The cells below that row belong to the liquid zone. + BOOST_CHECK_CLOSE(states[7].saturation(FluidSystem::oilPhaseIdx), 1.0, 1e-10); + BOOST_CHECK_SMALL(states[7].moleFraction(1) - 0.60, 1e-10); + } + } +} + +BOOST_AUTO_TEST_CASE(CompvdTwoZoneRejectsRowsAgainstContact) +{ + // Type 3 takes the gas-oil contact as its reference depth, so the rows may + // not put the phase change elsewhere. No item 10 accepts a vapour row below + // a liquid one. + const std::string gapAbove = "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2049 0 0.95 0.05 0 150.0\n" + " 2051 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n"; + const std::string interleaved = "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2040 0 0.60 0.40 1 150.0\n" + " 2060 0 0.95 0.05 0 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n"; + for (const auto& [equil, compvd, reason] : + {std::tuple {"EQUIL\n 2060 200 2300 0 2060 0 3* 3 /\n", + gapAbove, + "so the rows must agree with it"}, + std::tuple {"EQUIL\n 2072.5 200 2300 0 2050 0 /\n", + interleaved, + "must lie above the liquid rows"}}) { + BOOST_TEST_CONTEXT(equil << compvd) { + const EquilFixture fix(deckString(equil, "EQLDIMS\n/\n", "", compvd)); + BOOST_CHECK_EXCEPTION(fix.compute(std::vector(20, 0)), + std::runtime_error, + [reason](const std::runtime_error& e) { + return std::string_view {e.what()}.find(reason) + != std::string_view::npos; + }); + } + } +} + BOOST_AUTO_TEST_CASE(WaterZoneBelowContact) { // The water-oil contact at 2050 m lies inside the column. Above it the @@ -884,6 +1055,182 @@ BOOST_AUTO_TEST_CASE(WaterZoneBelowContact) BOOST_CHECK_GT(rhoWater, rhoHc); } +BOOST_AUTO_TEST_CASE(GasZoneMeetingTheWaterAnchorsIt) +{ + // The vapour rows reach 2060 m, so the gas-oil contact at 2000 m moves down + // to 2060 m, below the water-oil contact at 2052.5 m. The gas zone then + // meets the water, which takes its pressure from the gas at the contact + // rather than from a liquid column extended above its own zone. + const WaterEquilFixture fix(waterDeckString( + "EQUIL\n 2010 150 2052.5 0 2000 0 /\n", + "COMPVD\n" + " 2000 0 0.95 0.05 0 150.0\n" + " 2060 0 0.95 0.05 0 150.0\n" + " 2070 0 0.60 0.40 1 150.0\n" + " 2100 0 0.40 0.60 1 150.0 /\n")); + const auto states = fix.compute(std::vector(20, 0), + std::vector(20, connateSw), + std::vector(20, 1.0)).fluidStates(); + + // Cell 10 is centred on the water-oil contact, where the zero capillary + // pressure leaves the water and the gas at one pressure. + BOOST_REQUIRE_EQUAL(fix.depths[10], 2052.5); + BOOST_CHECK_CLOSE(states[10].pressure(WaterFluidSystem::waterPhaseIdx), + states[10].pressure(WaterFluidSystem::gasPhaseIdx), 1e-8); + BOOST_CHECK_CLOSE(states[10].saturation(WaterFluidSystem::gasPhaseIdx), + 1.0 - connateSw, 1e-10); + BOOST_CHECK_CLOSE(states[11].saturation(WaterFluidSystem::waterPhaseIdx), 1.0, 1e-10); +} + +BOOST_AUTO_TEST_CASE(EachWaterColumnTakesThePvtRegionOfItsCells) +{ + // The lower block holds brine: its water column is steeper than the fresh + // water column of the upper block, and than its own column would be in the + // first PVT region. + std::vector regions(20, 0); + std::fill(regions.begin() + 10, regions.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const WaterEquilFixture fix(twoWaterRegionsDeck("10*1 10*2", "10*1 10*2")); + const auto states = fix.compute(regions, connate, maxWater, regions).fluidStates(); + const auto first = fix.compute(regions, connate, maxWater).fluidStates(); + + const auto pressure = [](const auto& fs) { + return Opm::getValue(fs.pressure(WaterFluidSystem::waterPhaseIdx)); + }; + // Cells 8 and 9 lie below the upper contact at 2035 m, 18 and 19 below the + // lower one at 2085 m. + const Scalar fresh = impliedDensity(pressure(states[8]), pressure(states[9])); + const Scalar brine = impliedDensity(pressure(states[18]), pressure(states[19])); + const Scalar pFresh = 0.5 * (pressure(states[8]) + pressure(states[9])) / barsa; + const Scalar pBrine = 0.5 * (pressure(states[18]) + pressure(states[19])) / barsa; + BOOST_CHECK_CLOSE(fresh, waterDensity(1000.0, 1.03, 4.0e-5, pFresh), 0.01); + BOOST_CHECK_CLOSE(brine, waterDensity(1150.0, 1.01, 2.5e-5, pBrine), 0.01); + BOOST_CHECK_CLOSE(impliedDensity(pressure(first[18]), pressure(first[19])), + waterDensity(1000.0, 1.03, 4.0e-5, pBrine), 0.01); +} + +BOOST_AUTO_TEST_CASE(OneEquilibrationRegionHoldsTheWaterOfEachPvtRegion) +{ + // One equilibration region over both PVT regions, with its contact at 2035 m + // among the fresh-water cells. The hydrocarbon and the fresh water are as + // without the second PVT region, and the brine column starts from the same + // contact. + const std::vector eqlnum(20, 0); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2")); + const auto states = fix.compute(eqlnum, connate, maxWater, pvtnum).fluidStates(); + const auto first = fix.compute(eqlnum, connate, maxWater).fluidStates(); + + const auto pressure = [](const auto& fs, const unsigned phaseIdx) { + return Opm::getValue(fs.pressure(phaseIdx)); + }; + constexpr auto water = WaterFluidSystem::waterPhaseIdx; + for (std::size_t c = 0; c < states.size(); ++c) { + BOOST_TEST_CONTEXT("Cell " << c) { + BOOST_CHECK_EQUAL(pressure(states[c], WaterFluidSystem::oilPhaseIdx), + pressure(first[c], WaterFluidSystem::oilPhaseIdx)); + if (c < 10) { + BOOST_CHECK_EQUAL(pressure(states[c], water), pressure(first[c], water)); + } + } + } + + // Cells 18 and 19 hold brine. Cell 10 lies 17.5 m below the contact, where + // the brine exceeds the fresh water by the weight of the denser column. + const Scalar pBrine = 0.5 * (pressure(states[18], water) + pressure(states[19], water)); + BOOST_CHECK_CLOSE(impliedDensity(pressure(states[18], water), pressure(states[19], water)), + waterDensity(1150.0, 1.01, 2.5e-5, pBrine / barsa), 0.01); + const Scalar pCell10 = 0.5 * (pressure(states[10], water) + pressure(first[10], water)); + BOOST_CHECK_CLOSE(pressure(states[10], water) - pressure(first[10], water), + (waterDensity(1150.0, 1.01, 2.5e-5, pCell10 / barsa) + - waterDensity(1000.0, 1.03, 4.0e-5, pCell10 / barsa)) + * gravity * (fix.depths[10] - 2035.0), + 0.1); +} + +BOOST_AUTO_TEST_CASE(DatumInDifferentWatersIsRejected) +{ + // A datum below the contact gives the pressure of both the fresh water and + // the brine, which cannot both reach the contact at one pressure. + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2", + "EQUIL\n" + " 2045 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n")); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + BOOST_CHECK_EXCEPTION(fix.compute(std::vector(20, 0), + std::vector(20, connateSw), + std::vector(20, 1.0), pvtnum), + std::runtime_error, + [](const std::runtime_error& error) { + return std::string_view{error.what()}.find( + "the water of its PVT regions 1 and 2 reaches the contact") + != std::string_view::npos; + }); +} + +BOOST_AUTO_TEST_CASE(DatumInTheSameWaterIsAccepted) +{ + // PVT regions that share their water, as when they differ in the rock only, + // equilibrate as one. + const WaterEquilFixture fix(twoWaterRegionsDeck("20*1", "10*1 10*2", + "EQUIL\n" + " 2045 150 2035 0 2000 0 /\n" + " 2060 155 2085 0 2050 0 /\n", + "PVTW\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + " 75.0 1.03 4.0E-5 0.3 0.0 /\n" + "DENSITY\n" + " 800.0 1000.0 1.0 /\n" + " 800.0 1000.0 1.0 /\n")); + const std::vector eqlnum(20, 0); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const std::vector connate(20, connateSw); + const std::vector maxWater(20, 1.0); + const auto states = fix.compute(eqlnum, connate, maxWater, pvtnum).fluidStates(); + const auto single = fix.compute(eqlnum, connate, maxWater).fluidStates(); + for (std::size_t c = 0; c < states.size(); ++c) { + for (const auto phaseIdx : {WaterFluidSystem::oilPhaseIdx, + WaterFluidSystem::waterPhaseIdx}) { + BOOST_CHECK_EQUAL(Opm::getValue(states[c].pressure(phaseIdx)), + Opm::getValue(single[c].pressure(phaseIdx))); + } + } +} + +BOOST_AUTO_TEST_CASE(PvtRegionsMayVaryWithoutWater) +{ + // Without water the PVT regions take no part in the equilibration. + const EquilFixture fix(withTwoPvtRegions(deckString( + "EQUIL\n 2010 150 2300 0 2000 0 /\n", "EQLDIMS\n/\n", + "REGIONS\nPVTNUM\n10*1 10*2 /\n"))); + std::vector pvtnum(20, 0); + std::fill(pvtnum.begin() + 10, pvtnum.end(), 1); + const auto states = fix.compute(std::vector(20, 0), {}, {}, pvtnum).fluidStates(); + const auto single = fix.compute(std::vector(20, 0)).fluidStates(); + for (std::size_t c = 0; c < states.size(); ++c) { + BOOST_CHECK_EQUAL(Opm::getValue(states[c].pressure(FluidSystem::oilPhaseIdx)), + Opm::getValue(single[c].pressure(FluidSystem::oilPhaseIdx))); + } +} + +BOOST_AUTO_TEST_CASE(InvalidPvtnumIsRejected) +{ + const WaterEquilFixture fix(twoWaterRegionsDeck("10*1 10*2", "10*1 10*2")); + std::vector regions(20, 0); + std::fill(regions.begin() + 10, regions.end(), 1); + auto pvtnum = regions; + pvtnum.back() = 2; + BOOST_CHECK_THROW(fix.compute(regions, std::vector(20, connateSw), + std::vector(20, 1.0), pvtnum), + std::runtime_error); +} + BOOST_AUTO_TEST_CASE(CoincidentContactsKeepTheGasRoot) { // Contacts that coincide leave no liquid: every hydrocarbon cell sits above diff --git a/tests/test_flashphasepresence.cpp b/tests/test_flashphasepresence.cpp index ec616677f46..85fe8208bd2 100644 --- a/tests/test_flashphasepresence.cpp +++ b/tests/test_flashphasepresence.cpp @@ -59,11 +59,18 @@ struct TestFluidSystem void updatePhase(const FluidState&, unsigned) {} Eval molarVolume(unsigned) const { return Eval{1.0}; } Eval correctedMolarVolume(unsigned) const { return Eval{1.0}; } + void setRegionIndex(unsigned region) { regionIdx = region; } + unsigned regionIndex() const { return regionIdx; } + unsigned regionIdx = 0; }; + // Water is 100 kg/m3 denser in each further PVT region. template - static auto density(const FluidState&, const Cache&, unsigned) - { return typename FluidState::ValueType{1000.0}; } + static auto density(const FluidState&, const Cache& cache, unsigned phase) + { + return typename FluidState::ValueType{ + phase == waterPhaseIdx ? 1000.0 + 100.0 * cache.regionIndex() : 1000.0}; + } template static auto viscosity(const FluidState&, const Cache&, unsigned) @@ -176,6 +183,15 @@ struct TestProblem template double porosity(const Context&, unsigned, unsigned) const { return 0.2; } template + double rockCompressibility(const Context&, unsigned, unsigned) const { return compressibility; } + template + double rockReferencePressure(const Context&, unsigned, unsigned) const { return 5.0e6; } + template + unsigned pvtRegionIndex(const Context&, unsigned, unsigned) const { return pvtRegion; } + + double compressibility = 0.0; // 1/Pa + unsigned pvtRegion = 0; + template Dune::FieldMatrix intrinsicPermeability(const Context&, unsigned, unsigned) const { return Dune::FieldMatrix{1.0}; } }; @@ -283,3 +299,37 @@ BOOST_AUTO_TEST_CASE(PresenceIsRefreshedWhenHydrocarbonEntersAndLeaves) BOOST_CHECK(!iq.phaseIsPresent(1)); BOOST_CHECK_EQUAL(iq.saturationForOutput(2), 1.0); } + +BOOST_AUTO_TEST_CASE(PorosityFollowsTheRockCompressibility) +{ + // The pore volume grows with the pressure over the reference pressure, + // with the second-order expansion of exp(x) the black-oil model uses. + TestContext context; + context.problem_.compressibility = 1.0e-9; + context.priVars.values[2] = 0.5; + IntensiveQuantities iq; + iq.update(context, 0, 0); + + const double x = 1.0e-9 * (1.0e7 - 5.0e6); + BOOST_CHECK_CLOSE(iq.porosity().value(), 0.2 * (1.0 + x + 0.5 * x * x), 1e-12); + BOOST_CHECK_CLOSE(iq.porosity().derivative(0), 0.2 * 1.0e-9 * (1.0 + x), 1e-10); + + context.problem_.compressibility = 0.0; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.porosity().value(), 0.2); + BOOST_CHECK_EQUAL(iq.porosity().derivative(0), 0.0); +} + +BOOST_AUTO_TEST_CASE(WaterPropertiesTakeTheCellsPvtRegion) +{ + TestContext context; + context.priVars.values[2] = 0.5; + IntensiveQuantities iq; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.fluidState().density(2).value(), 1000.0); + + context.problem_.pvtRegion = 1; + iq.update(context, 0, 0); + BOOST_CHECK_EQUAL(iq.fluidState().density(2).value(), 1100.0); + BOOST_CHECK_EQUAL(iq.fluidState().density(0).value(), 1000.0); +}