diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 9febd9a1fde..cf07bbbe1cd 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1083,7 +1083,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/equil/EquilibrationHelpers.hpp opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp opm/simulators/flow/equil/InitStateEquil.hpp + opm/simulators/flow/equil/InitStateEquilComp.hpp opm/simulators/flow/equil/InitStateEquil_impl.hpp + opm/simulators/flow/equil/PressureFunction.hpp opm/simulators/flow/rescoup/ReservoirCouplingEnabled.hpp opm/simulators/wells/SegmentState.hpp opm/simulators/wells/WellContainer.hpp diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index f5f52ea57ce..687f7e8b2c0 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -1,7 +1,7 @@ // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- // vi: set et ts=4 sw=4 sts=4: /* - Copyright 2024 SINTEF Digital + Copyright 2024, 2026 SINTEF Digital This file is part of the Open Porous Media project (OPM). @@ -34,6 +34,7 @@ #include #include #include +#include #include @@ -141,17 +142,26 @@ class FlowProblemComp : public FlowProblem [&vg = this->simulator().vanguard()](const unsigned int it) { return vg.gridIdxToEquilGridIdx(it); }); updated = true; }; - // TODO: we might need to do the same with FlowProblemBlackoil for parallel finishTransmissibilities(); if (enableEclOutput_) { - eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities()); + // The output of TRANX, TRANY, TRANZ and NNC is on the whole grid: the + // I/O rank needs the global transmissibilities when running in parallel. + if (simulator.vanguard().grid().comm().size() > 1) { + if (simulator.vanguard().grid().comm().rank() == 0) { + eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility()); + } + } + else { + eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities()); + } std::function equilGridToGrid = [&simulator](unsigned int i) { return simulator.vanguard().gridEquilIdxToGridIdx(i); }; eclWriter_->extractOutputTransAndNNC(equilGridToGrid); } + simulator.vanguard().releaseGlobalTransmissibilities(); const auto& eclState = simulator.vanguard().eclState(); const auto& schedule = simulator.vanguard().schedule(); @@ -450,7 +460,29 @@ class FlowProblemComp : public FlowProblem void readEquilInitialCondition_() override { - throw std::logic_error("Equilibration is not supported by compositional modeling yet"); + const auto& simulator = this->simulator(); + const auto& vanguard = simulator.vanguard(); + const auto& eclState = vanguard.eclState(); + + // Zero-based equilibration region of every cell (EQLNUM, or region 0). + std::vector eqlnum(this->model().numGridDof(), 0); + if (eclState.fieldProps().has_int("EQLNUM")) { + const auto& e = eclState.fieldProps().get_int("EQLNUM"); + std::ranges::transform(e, eqlnum.begin(), [](const int r) { return r - 1; }); + } + + EQUIL::Comp::InitialStateComputer initialState( + eclState, + getEosType(), + vanguard.cellCenterDepths(), + eqlnum, + vanguard.gridView().comm(), + this->gravity()[dimWorld - 1], + this->numPressurePointsEquil()); + + initialFluidStates_ = std::move(initialState.fluidStates()); + // The primary variables are formed from the total composition; see initial(). + zmf_initialization_ = true; } void readEclRestartSolution_() diff --git a/opm/simulators/flow/NonlinearSystemCompositional.hpp b/opm/simulators/flow/NonlinearSystemCompositional.hpp index 2c883278d22..a9932c148b1 100644 --- a/opm/simulators/flow/NonlinearSystemCompositional.hpp +++ b/opm/simulators/flow/NonlinearSystemCompositional.hpp @@ -83,6 +83,8 @@ class NonlinearSystemCompositional : public NonlinearSystem void solveJacobianSystem(BVector& x); + void updateSolution(const BVector& dx); + bool hasNlddSolver() const { return false; } diff --git a/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp b/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp index db941b124db..12f804adba9 100644 --- a/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp +++ b/opm/simulators/flow/NonlinearSystemCompositional_impl.hpp @@ -274,6 +274,32 @@ relativeChange() const return resultDenom > 0.0 ? resultDelta / resultDenom : 0.0; } +template +void +NonlinearSystemCompositional:: +updateSolution(const BVector& dx) +{ + OPM_TIMEBLOCK(updateSolution); + + auto& model = this->simulator_.model(); + auto& solution = model.solution(/*timeIdx=*/0); + + model.newtonMethod().applyUpdate(/*nextSolution=*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx); + + // The linear solver leaves the rows of ghost cells untouched: fetch their + // updated primary variables from the owning processes before the intensive + // quantities are recomputed. + model.syncOverlap(); + + { + OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities); + model.invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } +} + template void NonlinearSystemCompositional:: @@ -306,7 +332,12 @@ reservoirResidualMetrics() const std::vector residualMetrics(numEq, 0.0); - for (unsigned dofIdx = 0; dofIdx < residual.size(); ++dofIdx) { + // Only the interior cells: the residual of a ghost cell misses the flux + // contributions of neighbors outside the overlap layer and does not + // converge, and its converged value lives on the owning process. + const auto& elemMapper = model.elementMapper(); + for (const auto& elem : elements(this->simulator_.gridView(), Dune::Partitions::interior)) { + const unsigned dofIdx = elemMapper.index(elem); if (dofIdx >= model.numGridDof() || model.dofTotalVolume(dofIdx) <= 0.0) { continue; } diff --git a/opm/simulators/flow/equil/InitStateEquil.hpp b/opm/simulators/flow/equil/InitStateEquil.hpp index 36edfd4aab8..346488677cd 100644 --- a/opm/simulators/flow/equil/InitStateEquil.hpp +++ b/opm/simulators/flow/equil/InitStateEquil.hpp @@ -34,6 +34,7 @@ #include #include +#include #include #include @@ -75,25 +76,6 @@ template class EquilReg; namespace Miscibility { template class RsFunction; } namespace Details { -template -class RK4IVP -{ -public: - RK4IVP(const RHS& f, - const std::array& span, - const Scalar y0, - const int N); - - Scalar operator()(const Scalar x) const; - -private: - int N_; - std::array span_; - std::vector y_; - std::vector f_; - - Scalar stepsize() const; -}; namespace PhasePressODE { template @@ -262,38 +244,7 @@ class PressureTable private: template - class PressureFunction - { - public: - struct InitCond { - Scalar depth; - Scalar pressure; - }; - - explicit PressureFunction(const ODE& ode, - const InitCond& ic, - const int nsample, - const VSpan& span); - - PressureFunction(const PressureFunction& rhs); - - PressureFunction(PressureFunction&& rhs) = default; - - PressureFunction& operator=(const PressureFunction& rhs); - - PressureFunction& operator=(PressureFunction&& rhs); - - Scalar value(const Scalar depth) const; - - private: - enum Direction : std::size_t { Up, Down, NumDir }; - - using Distribution = Details::RK4IVP; - using DistrPtr = std::unique_ptr; - - InitCond initial_; - std::array value_; - }; + using PressureFunction = Details::PressureFunction; using OilPressODE = PhasePressODE::Oil< FluidSystem, typename Region::CalcDissolution @@ -773,7 +724,7 @@ class InitialStateComputer PhaseSat& psat); template - void equilibrateTiltedFaultBlock(const CellRange& cells, + void equilibrateTiltedFaultBlock(const CellRange& cells, const EquilReg& eqreg, const GridView& gridView, const int numLevels, const PressTable& ptable, PhaseSat& psat); diff --git a/opm/simulators/flow/equil/InitStateEquilComp.hpp b/opm/simulators/flow/equil/InitStateEquilComp.hpp new file mode 100644 index 00000000000..f029d6f601f --- /dev/null +++ b/opm/simulators/flow/equil/InitStateEquilComp.hpp @@ -0,0 +1,509 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 SINTEF Digital + + 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 3 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 + * + * \brief Hydrostatic equilibration for the compositional simulator (EQUIL + ZMFVD). + */ +#ifndef OPM_INIT_STATE_EQUIL_COMP_HPP +#define OPM_INIT_STATE_EQUIL_COMP_HPP + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace Opm { +namespace EQUIL { +namespace Comp { + +namespace Details { + +/// Right-hand side of the hydrostatic ODE dp/ddepth = rho(depth, p) * g for a +/// fluid whose density follows from the cubic equation of state at the given +/// temperature and composition. The EOS root (liquid or vapour) is selected +/// by the phase index. +template +class EosDensityODE +{ +public: + using Scalar = typename FluidSystem::Scalar; + using CompVec = std::array; + using CompositionFunction = std::function; + using TabulatedFunction = Tabulated1DFunction; + + EosDensityODE(CompositionFunction composition, + const TabulatedFunction& tempVdTable, + const unsigned phaseIdx, + const CompositionalConfig::EOSType eosType, + const Scalar normGrav) + : composition_(std::move(composition)) + , tempVdTable_(tempVdTable) + , phaseIdx_(phaseIdx) + , eosType_(eosType) + , g_(normGrav) + {} + + Scalar operator()(const Scalar depth, + const Scalar press) const + { + const CompVec z = composition_(depth); + const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + + CompositionalFluidState fs; + fs.setTemperature(temp); + fs.setPressure(FluidSystem::oilPhaseIdx, press); + fs.setPressure(FluidSystem::gasPhaseIdx, press); + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { + fs.setMoleFraction(phaseIdx_, compIdx, z[compIdx]); + } + + typename FluidSystem::template ParameterCache paramCache(eosType_); + paramCache.updatePhase(fs, phaseIdx_); + + return FluidSystem::density(fs, paramCache, phaseIdx_) * g_; + } + +private: + CompositionFunction composition_; + const TabulatedFunction& tempVdTable_; + unsigned phaseIdx_; + CompositionalConfig::EOSType eosType_; + Scalar g_; +}; + +} // namespace Details + +/*! + * \brief Computes the initial state of a compositional model from hydrostatic + * equilibrium (the EQUIL and ZMFVD keywords). + * + * The composition versus depth is given by ZMFVD and the temperature by RTEMPVD + * (or the constant RTEMP). The phase pressures are obtained by integrating the + * hydrostatic ODE with the equation-of-state density, reusing the ODE machinery + * of the black-oil equilibration facility. Only the total composition, pressure + * and temperature are needed downstream: the phase split and the saturations are + * recomputed by the flash from these quantities. + * + * The supported initialization procedures (EQUIL item 10) are + * - type 1 (default): ZMFVD provides the total composition and the fluid is + * treated as a single phase throughout the column; + * - type 3: ZMFVD provides the liquid composition below the gas-oil contact. + * The contact acts as the datum, where the pressure is the saturation + * (bubble-point) pressure of the contact liquid unless EQUIL item 11 + * requests the given datum pressure. Above the contact the gas has the + * constant composition of the equilibrium vapour at the contact. + */ +template +class InitialStateComputer +{ +public: + using Scalar = typename FluidSystem::Scalar; + using FluidState = CompositionalFluidState; + + /// \param[in] eclipseState Input state, provides EQUIL, ZMFVD, RTEMP(VD). + /// \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] comm Communicator for parallel runs. + /// \param[in] gravity Norm of the gravity vector. + /// \param[in] numSamplePoints Sample points in each pressure integration. + InitialStateComputer(const EclipseState& eclipseState, + const CompositionalConfig::EOSType eosType, + const std::vector& cellCenterDepth, + const std::vector& eqlnum, + const Parallel::Communication& comm, + const Scalar gravity, + const int numSamplePoints) + : eosType_(eosType) + { + const auto& records = eclipseState.getInitConfig().getEquil(); + const auto& tables = eclipseState.getTableManager(); + + if (!tables.hasTables("ZMFVD")) { + // COMPVD is the other accepted way of giving the composition + // versus depth, but it is not supported here yet; name it so the + // message does not read as if ZMFVD were the only valid input. + const std::string msg = tables.hasTables("COMPVD") + ? "Equilibration of a compositional model with the composition versus " + "depth from COMPVD is not supported; use ZMFVD instead." + : "Equilibration of a compositional model requires the composition " + "versus depth from the ZMFVD keyword."; + OPM_THROW(std::runtime_error, msg); + } + + 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)); + } + + fluidStates_.resize(cellCenterDepth.size()); + for (std::size_t cell = 0; cell < cellCenterDepth.size(); ++cell) { + const auto region = eqlnum[cell]; + if (region < 0 || std::cmp_greater_equal(region, regions.size())) { + OPM_THROW(std::runtime_error, + fmt::format("Cell {} has EQLNUM {} outside the {} " + "equilibration regions.", + cell, region + 1, regions.size())); + } + assignCell(fluidStates_[cell], regions[region], cellCenterDepth[cell]); + } + } + + std::vector& fluidStates() + { return fluidStates_; } + + const std::vector& fluidStates() const + { return fluidStates_; } + +private: + using CompVec = std::array; + using TabulatedFunction = Tabulated1DFunction; + using ODE = Details::EosDensityODE; + using PressFunc = EQUIL::Details::PressureFunction; + + static constexpr int numComponents = FluidSystem::numComponents; + + /// The equilibrated vertical distributions within one region. + struct Region { + int initType{1}; // EQUIL item 10 + Scalar zgoc{}; + CompVec vaporComposition{}; // gas above the contact (type 3) + std::vector zmfVdTable; // per-component ZMFVD + TabulatedFunction tempVdTable; + std::optional oilPressure; + std::optional gasPressure; // type 3 only + }; + + static CompVec composition(const Region& reg, const Scalar depth) + { + CompVec z{}; + Scalar sum = 0.0; + for (int c = 0; c < numComponents; ++c) { + z[c] = std::max(Scalar{0}, reg.zmfVdTable[c].eval(depth, /*extrapolate=*/true)); + sum += z[c]; + } + if (!(sum > 0.0)) { + OPM_THROW(std::runtime_error, + fmt::format("The ZMFVD composition vanishes at depth {} m.", depth)); + } + std::ranges::transform(z, z.begin(), [sum](const Scalar zc) { return zc / sum; }); + return z; + } + + /// Whether the ZMFVD composition differs between the top and the bottom of + /// \p span, i.e. whether the table carries any variation over the region. + static bool compositionVariesAcross(const Region& reg, + const std::array& span) + { + const CompVec top = composition(reg, span[0]); + const CompVec bottom = composition(reg, span[1]); + for (int c = 0; c < numComponents; ++c) { + if (std::abs(top[c] - bottom[c]) > Scalar{1.0e-10}) { + return true; + } + } + return false; + } + + Region setupRegion(const EquilRecord& record, + const TableManager& tables, + const std::vector& cellCenterDepth, + const std::vector& eqlnum, + const Parallel::Communication& comm, + const Scalar gravity, + const int numSamplePoints, + const std::size_t regionIdx) const + { + Region reg; + + reg.initType = record.compositionalInitType(); + if (reg.initType != 1 && reg.initType != 3) { + OPM_THROW(std::runtime_error, + fmt::format("Compositional initialization type {} (EQUIL item 10) is " + "not supported for region {}; only type 1 (total " + "composition) and type 3 (liquid composition) are.", + reg.initType, regionIdx + 1)); + } + + reg.zgoc = record.gasOilContactDepth(); + + const auto& zmfvd = tables.getZmfvdTables().template getTable(regionIdx); + reg.zmfVdTable.resize(numComponents); + std::vector depths(zmfvd.getDepthColumn().begin(), + zmfvd.getDepthColumn().end()); + // A single row means a depth-independent composition; the interpolant + // needs two sample points, so duplicate it onto an arbitrary interval. + const bool constantComposition = (depths.size() == 1); + if (constantComposition) { + depths.push_back(depths.front() + Scalar{1}); + } + for (int c = 0; c < numComponents; ++c) { + const auto& col = zmfvd.getMoleFractionColumn(c); + std::vector values(col.begin(), col.end()); + if (constantComposition) { + values.push_back(values.front()); + } + reg.zmfVdTable[c].setXYContainers(depths, values); + } + + if (tables.hasTables("RTEMPVD")) { + const auto& rtempvd = tables.getRtempvdTables().template getTable(regionIdx); + reg.tempVdTable.setXYContainers(rtempvd.getDepthColumn(), + rtempvd.getTemperatureColumn()); + } + else { + const std::vector x{0.0, 1.0}; + const std::vector y(2, tables.rtemp()); + reg.tempVdTable.setXYContainers(x, y); + } + + // Vertical extent of the region's cells across all processes. + auto span = std::array{std::numeric_limits::max(), + std::numeric_limits::lowest()}; + for (std::size_t cell = 0; cell < cellCenterDepth.size(); ++cell) { + if (std::cmp_equal(eqlnum[cell], regionIdx)) { + span[0] = std::min(span[0], cellCenterDepth[cell]); + span[1] = std::max(span[1], cellCenterDepth[cell]); + } + } + span[0] = comm.min(span[0]); + span[1] = comm.max(span[1]); + if (span[0] > span[1]) { + // No cells anywhere in this region. + return reg; + } + if (span[1] - span[0] < Scalar{1}) { + // Avoid a degenerate integration interval. + span = {span[0] - Scalar{1}, span[1] + Scalar{1}}; + } + + // The equilibration covers the hydrocarbon column only. + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (record.waterOilContactDepth() < span[1]) { + OPM_THROW(std::runtime_error, + fmt::format("Compositional equilibration does not support a water " + "zone: the water-oil contact at {} m is above the " + "deepest cell centre at {} m of region {}.", + record.waterOilContactDepth(), span[1], regionIdx + 1)); + } + OpmLog::info(fmt::format("Equilibration region {}: the water phase is " + "initialized with zero saturation.", regionIdx + 1)); + } + + if (reg.initType == 1) { + setupSinglePhaseRegion(reg, record, span, gravity, numSamplePoints, regionIdx); + } + else { + setupTwoPhaseRegion(reg, record, span, gravity, numSamplePoints, regionIdx); + } + + return reg; + } + + /// EQUIL item 10 type 1: ZMFVD is the total composition and the fluid is a + /// single phase, integrated from the datum with the EOS density. The EOS + /// root is the vapour one if the datum lies in the gas zone (above the + /// gas-oil contact) and the liquid one otherwise. + /// + /// The datum-versus-contact test is enough to pick the root because of the + /// convention on the gas-oil contact (EQUIL item 5): it lies above the top + /// of the reservoir when there is no initial free gas, and below the bottom + /// when the region holds only gas. The defaulted item 5 (0 m, i.e. at the + /// surface) therefore expresses "no free gas" and correctly yields the + /// liquid root. + void setupSinglePhaseRegion(Region& reg, + const EquilRecord& record, + const std::array& span, + const Scalar gravity, + const int numSamplePoints, + const std::size_t regionIdx) const + { + const Scalar datum = record.datumDepth(); + const auto phaseIdx = (datum < reg.zgoc) + ? FluidSystem::gasPhaseIdx : FluidSystem::oilPhaseIdx; + + // Type 1 describes a continuous hydrocarbon phase, i.e. no gas-oil + // contact in the region. A contact placed inside the region is only + // labelled correctly when ZMFVD varies across it, so say so rather + // than initializing a two-phase column as a single phase in silence. + if ((reg.zgoc > span[0]) && (reg.zgoc < span[1]) && + !compositionVariesAcross(reg, span)) + { + OpmLog::warning(fmt::format("Equilibration region {}: EQUIL item 10 is 1 " + "(continuous hydrocarbon phase) but the gas-oil " + "contact at {} m lies inside the region, and ZMFVD " + "gives no compositional variation across it. The " + "phases cannot be labelled reliably; supply a " + "varying ZMFVD or use item 10 = 3.", + regionIdx + 1, reg.zgoc)); + } + + const ODE ode([®](const Scalar depth) { return composition(reg, depth); }, + reg.tempVdTable, phaseIdx, eosType_, gravity); + reg.oilPressure.emplace(ode, + typename PressFunc::InitCond{datum, Scalar(record.datumDepthPressure())}, + numSamplePoints, span); + + OpmLog::info(fmt::format("Equilibration region {}: single phase, total " + "composition specified (EQUIL item 10 is 1).", + regionIdx + 1)); + } + + /// EQUIL item 10 type 3: ZMFVD is the liquid composition. The pressure at + /// the gas-oil contact is the saturation pressure of the contact liquid + /// whenever the given datum pressure disagrees with it by an atmosphere or + /// more (or unconditionally the datum pressure when EQUIL item 11 is 1), + /// and the gas above the contact is the equilibrium vapour of the contact + /// liquid. + void setupTwoPhaseRegion(Region& reg, + const EquilRecord& record, + const std::array& span, + const Scalar gravity, + const int numSamplePoints, + const std::size_t regionIdx) const + { + if (std::abs(record.datumDepth() - reg.zgoc) > 0.0) { + OpmLog::warning(fmt::format("Equilibration region {}: the datum depth {} m " + "must be at the gas-oil contact when EQUIL " + "item 10 is 3; using the contact depth {} m.", + regionIdx + 1, record.datumDepth(), reg.zgoc)); + } + + const CompVec liquid = composition(reg, reg.zgoc); + const Scalar temp = reg.tempVdTable.eval(reg.zgoc, /*extrapolate=*/true); + Scalar psat{}; + CompVec vapor{}; + if (!SaturationPressure::bubblePressure(liquid, temp, eosType_, + psat, vapor)) { + OPM_THROW(std::runtime_error, + fmt::format("The saturation pressure calculation at the gas-oil " + "contact of region {} did not converge.", regionIdx + 1)); + } + reg.vaporComposition = vapor; + + // The datum pressure should already equal the saturation pressure at + // the contact. With EQUIL item 11 defaulted, the two are required to + // agree to within one atmosphere and the datum pressure is reset to the + // computed saturation pressure when they do not; item 11 = 1 keeps the + // given datum pressure whatever the outcome of that test. + constexpr Scalar oneAtmosphere = 101325.0; + const Scalar pDatum = record.datumDepthPressure(); + const bool resetToPsat = record.setToSaturationPressure() + && (std::abs(pDatum - psat) >= oneAtmosphere); + const Scalar pGoc = resetToPsat ? psat : pDatum; + + OpmLog::info(fmt::format("Equilibration region {}: two phases, liquid composition " + "specified (EQUIL item 10 is 3). The saturation pressure " + "at the gas-oil contact ({} m) is {:.6g} bar.", + regionIdx + 1, reg.zgoc, psat / 1e5)); + + if (resetToPsat) { + OpmLog::warning(fmt::format("Equilibration region {}: the datum pressure {:.6g} bar " + "differs from the saturation pressure {:.6g} bar at the " + "gas-oil contact by more than one atmosphere; the " + "saturation pressure is used instead.", + regionIdx + 1, pDatum / 1e5, psat / 1e5)); + } + + const ODE oilOde([®](const Scalar depth) { return composition(reg, depth); }, + reg.tempVdTable, FluidSystem::oilPhaseIdx, eosType_, gravity); + reg.oilPressure.emplace(oilOde, + typename PressFunc::InitCond{reg.zgoc, pGoc}, + numSamplePoints, span); + + const ODE gasOde([vapor](const Scalar) { return vapor; }, + reg.tempVdTable, FluidSystem::gasPhaseIdx, eosType_, gravity); + const Scalar pcgoc = record.gasOilContactCapillaryPressure(); + reg.gasPressure.emplace(gasOde, + typename PressFunc::InitCond{reg.zgoc, pGoc + pcgoc}, + numSamplePoints, span); + } + + void assignCell(FluidState& fs, const Region& reg, const Scalar depth) const + { + const bool inGasZone = (reg.initType == 3) && (depth < reg.zgoc); + + const CompVec z = inGasZone ? reg.vaporComposition : composition(reg, depth); + const auto& pressFunc = inGasZone ? reg.gasPressure : reg.oilPressure; + if (!pressFunc.has_value()) { + OPM_THROW(std::runtime_error, + "Evaluating the equilibrated pressure of a region without cells."); + } + const Scalar press = pressFunc->value(depth); + + fs.setTemperature(reg.tempVdTable.eval(depth, /*extrapolate=*/true)); + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (FluidSystem::phaseIsActive(phaseIdx)) { + fs.setPressure(phaseIdx, press); + fs.setSaturation(phaseIdx, 0.0); + } + } + // Nominal single-phase saturation; the flash recomputes the phase split + // from the total composition, the pressure and the temperature. + fs.setSaturation(inGasZone ? FluidSystem::gasPhaseIdx : FluidSystem::oilPhaseIdx, 1.0); + + for (int c = 0; c < numComponents; ++c) { + fs.setMoleFraction(c, z[c]); + } + } + + CompositionalConfig::EOSType eosType_; + std::vector fluidStates_; +}; + +} // namespace Comp +} // namespace EQUIL +} // namespace Opm + +#endif // OPM_INIT_STATE_EQUIL_COMP_HPP diff --git a/opm/simulators/flow/equil/InitStateEquil_impl.hpp b/opm/simulators/flow/equil/InitStateEquil_impl.hpp index 0ee10d8b119..fcb9a35659f 100644 --- a/opm/simulators/flow/equil/InitStateEquil_impl.hpp +++ b/opm/simulators/flow/equil/InitStateEquil_impl.hpp @@ -315,75 +315,6 @@ Scalar calculateTrueVerticalDepth(Scalar z, Scalar x, Scalar y, return tvd; } -template -RK4IVP::RK4IVP(const RHS& f, - const std::array& span, - const Scalar y0, - const int N) - : N_(N) - , span_(span) -{ - const Scalar h = stepsize(); - const Scalar h2 = h / 2; - const Scalar h6 = h / 6; - - y_.reserve(N + 1); - f_.reserve(N + 1); - - y_.push_back(y0); - f_.push_back(f(span_[0], y0)); - - for (int i = 0; i < N; ++i) { - const Scalar x = span_[0] + i*h; - const Scalar y = y_.back(); - - const Scalar k1 = f_[i]; - const Scalar k2 = f(x + h2, y + h2*k1); - const Scalar k3 = f(x + h2, y + h2*k2); - const Scalar k4 = f(x + h, y + h*k3); - - y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); - f_.push_back(f(x + h, y_.back())); - } - - assert (y_.size() == typename std::vector::size_type(N + 1)); -} - -template -Scalar RK4IVP:: -operator()(const Scalar x) const -{ - // Dense output (O(h**3)) according to Shampine - // (Hermite interpolation) - const Scalar h = stepsize(); - int i = (x - span_[0]) / h; - - // Crude handling of evaluation point outside "span_"; - if (i < 0) { i = 0; } - if (N_ <= i) { i = N_ - 1; } - - // Relative to the interval actually used, so that a point at the end of - // the span lands on t = 1 of the final interval rather than t = 0. - const Scalar t = (x - (span_[0] + i*h)) / h; - - const Scalar y0 = y_[i], y1 = y_[i + 1]; - const Scalar f0 = f_[i], f1 = f_[i + 1]; - - Scalar u = (1 - 2*t) * (y1 - y0); - u += h * ((t - 1)*f0 + t*f1); - u *= t * (t - 1); - u += (1 - t)*y0 + t*y1; - - return u; -} - -template -Scalar RK4IVP:: -stepsize() const -{ - return (span_[1] - span_[0]) / N_; -} - namespace PhasePressODE { template @@ -572,88 +503,6 @@ density(const Scalar depth, } -template -template -PressureTable:: -PressureFunction::PressureFunction(const ODE& ode, - const InitCond& ic, - const int nsample, - const VSpan& span) - : initial_(ic) -{ - this->value_[Direction::Up] = std::make_unique - (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample); - - this->value_[Direction::Down] = std::make_unique - (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample); -} - -template -template -PressureTable:: -PressureFunction::PressureFunction(const PressureFunction& rhs) - : initial_(rhs.initial_) -{ - this->value_[Direction::Up] = - std::make_unique(*rhs.value_[Direction::Up]); - - this->value_[Direction::Down] = - std::make_unique(*rhs.value_[Direction::Down]); -} - -template -template -typename PressureTable::template PressureFunction& -PressureTable:: -PressureFunction:: -operator=(const PressureFunction& rhs) -{ - this->initial_ = rhs.initial_; - - this->value_[Direction::Up] = - std::make_unique(*rhs.value_[Direction::Up]); - - this->value_[Direction::Down] = - std::make_unique(*rhs.value_[Direction::Down]); - - return *this; -} - -template -template -typename PressureTable::template PressureFunction& -PressureTable:: -PressureFunction:: -operator=(PressureFunction&& rhs) -{ - this->initial_ = rhs.initial_; - this->value_ = std::move(rhs.value_); - - return *this; -} - -template -template -typename PressureTable::Scalar -PressureTable:: -PressureFunction:: -value(const Scalar depth) const -{ - if (depth < this->initial_.depth) { - // Value above initial condition depth. - return (*this->value_[Direction::Up])(depth); - } - else if (depth > this->initial_.depth) { - // Value below initial condition depth. - return (*this->value_[Direction::Down])(depth); - } - else { - // Value *at* initial condition depth. - return this->initial_.pressure; - } -} - - template template void PressureTable:: diff --git a/opm/simulators/flow/equil/PressureFunction.hpp b/opm/simulators/flow/equil/PressureFunction.hpp new file mode 100644 index 00000000000..fba664e7c30 --- /dev/null +++ b/opm/simulators/flow/equil/PressureFunction.hpp @@ -0,0 +1,208 @@ +// -*- 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 3 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 + * + * \brief The ODE integrator and phase-pressure function used to solve the + * hydrostatic equilibrium problem, shared by the black-oil and the + * compositional equilibration facilities. + */ +#ifndef OPM_EQUIL_PRESSURE_FUNCTION_HPP +#define OPM_EQUIL_PRESSURE_FUNCTION_HPP + +#include +#include +#include +#include + +namespace Opm { +namespace EQUIL { +namespace Details { + +/// Fourth-order Runge-Kutta integrator for the initial-value problem +/// y'(x) = f(x, y), y(span[0]) = y0 +/// on N equidistant steps covering span, with O(h^3) dense output. +template +class RK4IVP +{ +public: + RK4IVP(const RHS& f, + const std::array& span, + const Scalar y0, + const int N) + : N_(N) + , span_(span) + { + const Scalar h = stepsize(); + const Scalar h2 = h / 2; + const Scalar h6 = h / 6; + + y_.reserve(N + 1); + f_.reserve(N + 1); + + y_.push_back(y0); + f_.push_back(f(span_[0], y0)); + + for (int i = 0; i < N; ++i) { + const Scalar x = span_[0] + i*h; + const Scalar y = y_.back(); + + const Scalar k1 = f_[i]; + const Scalar k2 = f(x + h2, y + h2*k1); + const Scalar k3 = f(x + h2, y + h2*k2); + const Scalar k4 = f(x + h, y + h*k3); + + y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); + f_.push_back(f(x + h, y_.back())); + } + + assert (y_.size() == typename std::vector::size_type(N + 1)); + } + + Scalar operator()(const Scalar x) const + { + // Dense output (O(h**3)) according to Shampine + // (Hermite interpolation) + const Scalar h = stepsize(); + int i = (x - span_[0]) / h; + // Crude handling of evaluation point outside "span_"; + if (i < 0) { i = 0; } + if (N_ <= i) { i = N_ - 1; } + + // Relative to the interval actually used, so that a point at the end of + // the span lands on t = 1 of the final interval rather than t = 0. + const Scalar t = (x - (span_[0] + i*h)) / h; + + const Scalar y0 = y_[i], y1 = y_[i + 1]; + const Scalar f0 = f_[i], f1 = f_[i + 1]; + + Scalar u = (1 - 2*t) * (y1 - y0); + u += h * ((t - 1)*f0 + t*f1); + u *= t * (t - 1); + u += (1 - t)*y0 + t*y1; + + return u; + } + +private: + int N_; + std::array span_; + std::vector y_; + std::vector f_; + + Scalar stepsize() const + { return (span_[1] - span_[0]) / N_; } +}; + +/// Phase pressure as a function of depth, obtained by integrating the +/// hydrostatic ODE +/// dp/ddepth = ODE(depth, p) +/// upwards and downwards from an initial condition (depth, pressure). +template +class PressureFunction +{ +public: + using VSpan = std::array; + + struct InitCond { + Scalar depth; + Scalar pressure; + }; + + explicit PressureFunction(const ODE& ode, + const InitCond& ic, + const int nsample, + const VSpan& span) + : initial_(ic) + { + this->value_[Direction::Up] = std::make_unique + (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample); + + this->value_[Direction::Down] = std::make_unique + (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample); + } + + PressureFunction(const PressureFunction& rhs) + : initial_(rhs.initial_) + { + this->value_[Direction::Up] = + std::make_unique(*rhs.value_[Direction::Up]); + + this->value_[Direction::Down] = + std::make_unique(*rhs.value_[Direction::Down]); + } + + PressureFunction(PressureFunction&& rhs) = default; + + PressureFunction& operator=(const PressureFunction& rhs) + { + this->initial_ = rhs.initial_; + + this->value_[Direction::Up] = + std::make_unique(*rhs.value_[Direction::Up]); + + this->value_[Direction::Down] = + std::make_unique(*rhs.value_[Direction::Down]); + + return *this; + } + + PressureFunction& operator=(PressureFunction&& rhs) + { + this->initial_ = rhs.initial_; + this->value_ = std::move(rhs.value_); + + return *this; + } + + Scalar value(const Scalar depth) const + { + if (depth < this->initial_.depth) { + // Value above initial condition depth. + return (*this->value_[Direction::Up])(depth); + } + else if (depth > this->initial_.depth) { + // Value below initial condition depth. + return (*this->value_[Direction::Down])(depth); + } + else { + // Value *at* initial condition depth. + return this->initial_.pressure; + } + } + +private: + enum Direction : std::size_t { Up, Down, NumDir }; + + using Distribution = RK4IVP; + using DistrPtr = std::unique_ptr; + + InitCond initial_; + std::array value_; +}; + +} // namespace Details +} // namespace EQUIL +} // namespace Opm + +#endif // OPM_EQUIL_PRESSURE_FUNCTION_HPP diff --git a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp index 7193cfc9d11..510d32b65b8 100644 --- a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp @@ -381,8 +381,8 @@ partiallySupported() "EQUIL", { {9,{true, [](int x) { return x >= -20 && x <= 20; }, "EQUIL(OIP_INIT): only values between -20 and 20 are allowed (default is -5)"}}, // OIP_INIT - {10,{false, allow_values {}, "EQUIL(COMP_INIT_TYPE): compositional option not used, should be defaulted"}}, // COMP_INIT_TYPE - {11,{false, allow_values {}, "EQUIL(COMP_NOT_SET_SAT_PRESSURE): compositional option not used, should be defaulted"}}, // COMP_NOT_SET_SAT_PRESSURE + {10,{false, allow_values {1, 3}, "EQUIL(COMP_INIT_TYPE): only the compositional initialization types 1 (continuous hydrocarbon phase) and 3 (liquid composition at the contact) are supported"}}, // COMP_INIT_TYPE + {11,{false, allow_values {0, 1}, "EQUIL(COMP_NOT_SET_SAT_PRESSURE): only values 0 and 1 exist"}}, // COMP_NOT_SET_SAT_PRESSURE }, }, { diff --git a/opm/simulators/utils/PropsDataHandle.hpp b/opm/simulators/utils/PropsDataHandle.hpp index c32ac994d31..488b2ac0d30 100644 --- a/opm/simulators/utils/PropsDataHandle.hpp +++ b/opm/simulators/utils/PropsDataHandle.hpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -72,12 +73,23 @@ class PropsDataHandle m_intKeys = globalProps.keys(); m_doubleKeys = globalProps.keys(); m_distributed_fieldProps.copyTran(globalProps); + + // Some keywords, e.g. the compositional ZMF, carry multiple values + // per cell (stored with the cell index varying fastest). + m_doubleMult.reserve(m_doubleKeys.size()); + for (const auto& doubleKey : m_doubleKeys) + { + const auto& fieldData = globalProps.get_double_field_data(doubleKey, + /* allow_unsupported = */ true); + m_doubleMult.push_back(fieldData.numValuePerCell()); + } } Parallel::MpiSerializer ser(comm); ser.broadcast(Parallel::RootRank{0}, *this); - m_no_data = m_intKeys.size() + m_doubleKeys.size(); + m_no_data = m_intKeys.size() + + std::accumulate(m_doubleMult.begin(), m_doubleMult.end(), std::size_t{0}); if (comm.rank() == 0) { const FieldPropsManager& globalProps = eclState.globalFieldProps(); @@ -86,6 +98,7 @@ class PropsDataHandle using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); + const std::size_t numCells = gridView.size(0); for (const auto &element : elements(gridView, Dune::Partitions::interiorBorder)) { @@ -101,14 +114,18 @@ class PropsDataHandle static_cast(fieldData.value_status[index])); } - for (const auto& doubleKey : m_doubleKeys) + for (std::size_t keyIdx = 0; keyIdx < m_doubleKeys.size(); ++keyIdx) { // We need to allow unsupported keywords to get the data // for TranCalculator, too. - const auto& fieldData = globalProps.get_double_field_data(doubleKey, + const auto& fieldData = globalProps.get_double_field_data(m_doubleKeys[keyIdx], /* allow_unsupported = */ true); - data.emplace_back(fieldData.data[index], - static_cast(fieldData.value_status[index])); + for (std::size_t comp = 0; comp < m_doubleMult[keyIdx]; ++comp) + { + const auto dataIdx = comp * numCells + index; + data.emplace_back(fieldData.data[dataIdx], + static_cast(fieldData.value_status[dataIdx])); + } } } } @@ -117,16 +134,18 @@ class PropsDataHandle ~PropsDataHandle() { // distributed grid is now correctly set up. + const std::size_t numCells = m_grid.size(0); for (const auto& intKey : m_intKeys) { - m_distributed_fieldProps.m_intProps[intKey].data.resize(m_grid.size(0)); - m_distributed_fieldProps.m_intProps[intKey].value_status.resize(m_grid.size(0)); + m_distributed_fieldProps.m_intProps[intKey].data.resize(numCells); + m_distributed_fieldProps.m_intProps[intKey].value_status.resize(numCells); } - for (const auto& doubleKey : m_doubleKeys) + for (std::size_t keyIdx = 0; keyIdx < m_doubleKeys.size(); ++keyIdx) { - m_distributed_fieldProps.m_doubleProps[doubleKey].data.resize(m_grid.size(0)); - m_distributed_fieldProps.m_doubleProps[doubleKey].value_status.resize(m_grid.size(0)); + auto& props = m_distributed_fieldProps.m_doubleProps[m_doubleKeys[keyIdx]]; + props.data.resize(m_doubleMult[keyIdx] * numCells); + props.value_status.resize(m_doubleMult[keyIdx] * numCells); } // copy data for the persistent mao to the field properties @@ -151,11 +170,15 @@ class PropsDataHandle m_distributed_fieldProps.m_intProps[intKey].value_status[index] = static_cast(pair.second); } - for (const auto& doubleKey : m_doubleKeys) + for (std::size_t keyIdx = 0; keyIdx < m_doubleKeys.size(); ++keyIdx) { - const auto& pair = data->second[counter++]; - m_distributed_fieldProps.m_doubleProps[doubleKey].data[index] = pair.first; - m_distributed_fieldProps.m_doubleProps[doubleKey].value_status[index] = static_cast(pair.second); + auto& props = m_distributed_fieldProps.m_doubleProps[m_doubleKeys[keyIdx]]; + for (std::size_t comp = 0; comp < m_doubleMult[keyIdx]; ++comp) + { + const auto& pair = data->second[counter++]; + props.data[comp * numCells + index] = pair.first; + props.value_status[comp * numCells + index] = static_cast(pair.second); + } } } } @@ -208,6 +231,7 @@ class PropsDataHandle { serializer(m_intKeys); serializer(m_doubleKeys); + serializer(m_doubleMult); m_distributed_fieldProps.serializeOp(serializer); } @@ -220,6 +244,8 @@ class PropsDataHandle std::vector m_intKeys; //! \brief The names of the keys of the double fields. std::vector m_doubleKeys; + //! \brief The number of values per cell of each double field. + std::vector m_doubleMult; /// \brief The data per element as a vector mapped from the local id. /// /// each entry is a pair of data and value_status. diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp index a2b3a1cb043..f131ad39c14 100644 --- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp +++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp @@ -730,7 +730,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords() {"WTHPMAX", {true, std::nullopt}}, {"ZIPPY2", {false, std::nullopt}}, {"ZIPP2OFF", {false, std::nullopt}}, - {"ZMFVD", {false, std::nullopt}}, }; return unsupported_keywords;