diff --git a/ebos/eclgenerictracermodel.cc b/ebos/eclgenerictracermodel.cc index 3a1c08f49..290d5d038 100644 --- a/ebos/eclgenerictracermodel.cc +++ b/ebos/eclgenerictracermodel.cc @@ -20,374 +20,10 @@ module for the precise wording of the license and the list of copyright holders. */ - #include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include - -#if HAVE_DUNE_FEM -#include -#include -#include -#endif // HAVE_DUNE_FEM - -#if HAVE_DUNE_ALUGRID -#include -#include -#include "alucartesianindexmapper.hh" -#endif // HAVE_DUNE_ALUGRID - -#include -#include -#include -#include -#include -#include -#include +#include "eclgenerictracermodel_impl.hh" namespace Opm { - -#if HAVE_MPI -template -struct TracerSolverSelector -{ - using Comm = Dune::OwnerOverlapCopyCommunication; - using TracerOperator = Dune::OverlappingSchwarzOperator; - using type = Dune::FlexibleSolver; -}; - -template -std::tuple>>, - std::unique_ptr::type>> -createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&) -{ - OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers."); - return {nullptr, nullptr}; -} - -template -std::tuple>>, - std::unique_ptr::type>> -createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm) -{ - using TracerOperator = Dune::OverlappingSchwarzOperator>; - using TracerSolver = Dune::FlexibleSolver; - const auto& cellComm = grid.cellCommunication(); - auto op = std::make_unique(M, cellComm); - auto dummyWeights = [](){ return Vector();}; - return {std::move(op), std::make_unique(*op, cellComm, prm, dummyWeights, 0)}; -} -#endif - -template -EclGenericTracerModel:: -EclGenericTracerModel(const GridView& gridView, - const EclipseState& eclState, - const CartesianIndexMapper& cartMapper, - const DofMapper& dofMapper, - const std::function(int)> centroids) - : gridView_(gridView) - , eclState_(eclState) - , cartMapper_(cartMapper) - , dofMapper_(dofMapper) - , centroids_(centroids) -{ -} - -template -Scalar EclGenericTracerModel:: -tracerConcentration(int tracerIdx, int globalDofIdx) const -{ - if (tracerConcentration_.empty()) - return 0.0; - - return tracerConcentration_[tracerIdx][globalDofIdx]; -} - -template -void EclGenericTracerModel:: -setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) -{ - this->tracerConcentration_[tracerIdx][globalDofIdx] = value; -} - -template -int EclGenericTracerModel:: -numTracers() const -{ - return this->eclState_.tracer().size(); -} - -template -std::string EclGenericTracerModel:: -fname(int tracerIdx) const -{ - return this->eclState_.tracer()[tracerIdx].fname(); -} - -template -double EclGenericTracerModel:: -currentConcentration_(const Well& eclWell, const std::string& name) const -{ - return eclWell.getTracerProperties().getConcentration(name); -} - -template -const std::string& EclGenericTracerModel:: -name(int tracerIdx) const -{ - return this->eclState_.tracer()[tracerIdx].name; -} - -template -void EclGenericTracerModel:: -doInit(bool rst, size_t numGridDof, - size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx) -{ - const auto& tracers = eclState_.tracer(); - - if (tracers.size() == 0) - return; // tracer treatment is supposed to be disabled - - // retrieve the number of tracers from the deck - const size_t numTracers = tracers.size(); - tracerConcentration_.resize(numTracers); - storageOfTimeIndex1_.resize(numTracers); - - // the phase where the tracer is - tracerPhaseIdx_.resize(numTracers); - for (size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) { - const auto& tracer = tracers[tracerIdx]; - - if (tracer.phase == Phase::WATER) - tracerPhaseIdx_[tracerIdx] = waterPhaseIdx; - else if (tracer.phase == Phase::OIL) - tracerPhaseIdx_[tracerIdx] = oilPhaseIdx; - else if (tracer.phase == Phase::GAS) - tracerPhaseIdx_[tracerIdx] = gasPhaseIdx; - - tracerConcentration_[tracerIdx].resize(numGridDof); - storageOfTimeIndex1_[tracerIdx].resize(numGridDof); - - if (rst) - continue; - - - //TBLK keyword - if (tracer.free_concentration.has_value()){ - const auto& free_concentration = tracer.free_concentration.value(); - int tblkDatasize = free_concentration.size(); - if (tblkDatasize < cartMapper_.cartesianSize()){ - throw std::runtime_error("Wrong size of TBLK for" + tracer.name); - } - for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){ - int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx); - tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx]; - } - } - //TVDPF keyword - else if (tracer.free_tvdp.has_value()) { - const auto& free_tvdp = tracer.free_tvdp.value(); - for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){ - tracerConcentration_[tracerIdx][globalDofIdx] = - free_tvdp.evaluate("TRACER_CONCENTRATION", - centroids_(globalDofIdx)[2]); - } - } else - throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name)); - } - - // allocate matrix for storing the Jacobian of the tracer residual - tracerMatrix_ = std::make_unique(numGridDof, numGridDof, TracerMatrix::random); - - // find the sparsity pattern of the tracer matrix - using NeighborSet = std::set; - std::vector neighbors(numGridDof); - - Stencil stencil(gridView_, dofMapper_); - for (const auto& elem : elements(gridView_)) { - stencil.update(elem); - - for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) { - unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx); - - for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { - unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); - neighbors[myIdx].insert(neighborIdx); - } - } - } - - // allocate space for the rows of the matrix - for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) - tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size()); - tracerMatrix_->endrowsizes(); - - // fill the rows with indices. each degree of freedom talks to - // all of its neighbors. (it also talks to itself since - // degrees of freedom are sometimes quite egocentric.) - for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) { - typename NeighborSet::iterator nIt = neighbors[dofIdx].begin(); - typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end(); - for (; nIt != nEndIt; ++nIt) - tracerMatrix_->addindex(dofIdx, *nIt); - } - tracerMatrix_->endindices(); - - const int sizeCartGrid = cartMapper_.cartesianSize(); - cartToGlobal_.resize(sizeCartGrid); - for (unsigned i = 0; i < numGridDof; ++i) { - int cartIdx = cartMapper_.cartesianIndex(i); - cartToGlobal_[cartIdx] = i; - } -} - -template -bool EclGenericTracerModel:: -linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) -{ -#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7) - Dune::FMatrixPrecision::set_singular_limit(1.e-30); - Dune::FMatrixPrecision::set_absolute_limit(1.e-30); -#endif - x = 0.0; - Scalar tolerance = 1e-2; - int maxIter = 100; - - int verbosity = 0; - PropertyTree prm; - prm.put("maxiter", maxIter); - prm.put("tol", tolerance); - prm.put("verbosity", verbosity); - prm.put("solver", std::string("bicgstab")); - prm.put("preconditioner.type", std::string("ParOverILU0")); - -#if HAVE_MPI - if(gridView_.grid().comm().size() > 1) - { - auto [tracerOperator, solver] = - createParallelFlexibleSolver(gridView_.grid(), M, prm); - (void) tracerOperator; - - Dune::InverseOperatorResult result; - solver->apply(x, b, result); - - // return the result of the solver - return result.converged; - } - else - { -#endif - using TracerSolver = Dune::BiCGSTABSolver; - using TracerOperator = Dune::MatrixAdapter; - using TracerScalarProduct = Dune::SeqScalarProduct; - using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>; - - TracerOperator tracerOperator(M); - TracerScalarProduct tracerScalarProduct; - TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0 - - TracerSolver solver (tracerOperator, tracerScalarProduct, - tracerPreconditioner, tolerance, maxIter, - verbosity); - - Dune::InverseOperatorResult result; - solver.apply(x, b, result); - - // return the result of the solver - return result.converged; -#if HAVE_MPI - } -#endif -} - -template -bool EclGenericTracerModel:: -linearSolveBatchwise_(const TracerMatrix& M, std::vector& x, std::vector& b) -{ -#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7) - Dune::FMatrixPrecision::set_singular_limit(1.e-30); - Dune::FMatrixPrecision::set_absolute_limit(1.e-30); -#endif - Scalar tolerance = 1e-2; - int maxIter = 100; - - int verbosity = 0; - PropertyTree prm; - prm.put("maxiter", maxIter); - prm.put("tol", tolerance); - prm.put("verbosity", verbosity); - prm.put("solver", std::string("bicgstab")); - prm.put("preconditioner.type", std::string("ParOverILU0")); - -#if HAVE_MPI - if(gridView_.grid().comm().size() > 1) - { - auto [tracerOperator, solver] = - createParallelFlexibleSolver(gridView_.grid(), M, prm); - (void) tracerOperator; - bool converged = true; - for (size_t nrhs =0; nrhs < b.size(); ++nrhs) { - x[nrhs] = 0.0; - Dune::InverseOperatorResult result; - solver->apply(x[nrhs], b[nrhs], result); - converged = (converged && result.converged); - } - return converged; - } - else - { -#endif - using TracerSolver = Dune::BiCGSTABSolver; - using TracerOperator = Dune::MatrixAdapter; - using TracerScalarProduct = Dune::SeqScalarProduct; - using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>; - - TracerOperator tracerOperator(M); - TracerScalarProduct tracerScalarProduct; - TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0 - - TracerSolver solver (tracerOperator, tracerScalarProduct, - tracerPreconditioner, tolerance, maxIter, - verbosity); - - bool converged = true; - for (size_t nrhs =0; nrhs < b.size(); ++nrhs) { - x[nrhs] = 0.0; - Dune::InverseOperatorResult result; - solver.apply(x[nrhs], b[nrhs], result); - converged = (converged && result.converged); - } - - // return the result of the solver - return converged; -#if HAVE_MPI - } -#endif -} - #if HAVE_DUNE_FEM template class EclGenericTracerModel>>, diff --git a/ebos/eclgenerictracermodel_impl.hh b/ebos/eclgenerictracermodel_impl.hh new file mode 100644 index 000000000..8c5dde4f6 --- /dev/null +++ b/ebos/eclgenerictracermodel_impl.hh @@ -0,0 +1,399 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/** + * \file + * + * \copydoc Opm::EclTracerModel + */ +#ifndef EWOMS_ECL_GENERIC_TRACER_MODEL_IMPL_HH +#define EWOMS_ECL_GENERIC_TRACER_MODEL_IMPL_HH + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif // HAVE_DUNE_FEM + +#if HAVE_DUNE_ALUGRID +#include +#include +#include "alucartesianindexmapper.hh" +#endif // HAVE_DUNE_ALUGRID + +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +#if HAVE_MPI +template +struct TracerSolverSelector +{ + using Comm = Dune::OwnerOverlapCopyCommunication; + using TracerOperator = Dune::OverlappingSchwarzOperator; + using type = Dune::FlexibleSolver; +}; + +template +std::tuple>>, + std::unique_ptr::type>> +createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&) +{ + OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers."); + return {nullptr, nullptr}; +} + +template +std::tuple>>, + std::unique_ptr::type>> +createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm) +{ + using TracerOperator = Dune::OverlappingSchwarzOperator>; + using TracerSolver = Dune::FlexibleSolver; + const auto& cellComm = grid.cellCommunication(); + auto op = std::make_unique(M, cellComm); + auto dummyWeights = [](){ return Vector();}; + return {std::move(op), std::make_unique(*op, cellComm, prm, dummyWeights, 0)}; +} +#endif + +template +EclGenericTracerModel:: +EclGenericTracerModel(const GridView& gridView, + const EclipseState& eclState, + const CartesianIndexMapper& cartMapper, + const DofMapper& dofMapper, + const std::function(int)> centroids) + : gridView_(gridView) + , eclState_(eclState) + , cartMapper_(cartMapper) + , dofMapper_(dofMapper) + , centroids_(centroids) +{ +} + +template +Scalar EclGenericTracerModel:: +tracerConcentration(int tracerIdx, int globalDofIdx) const +{ + if (tracerConcentration_.empty()) + return 0.0; + + return tracerConcentration_[tracerIdx][globalDofIdx]; +} + +template +void EclGenericTracerModel:: +setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value) +{ + this->tracerConcentration_[tracerIdx][globalDofIdx] = value; +} + +template +int EclGenericTracerModel:: +numTracers() const +{ + return this->eclState_.tracer().size(); +} + +template +std::string EclGenericTracerModel:: +fname(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].fname(); +} + +template +double EclGenericTracerModel:: +currentConcentration_(const Well& eclWell, const std::string& name) const +{ + return eclWell.getTracerProperties().getConcentration(name); +} + +template +const std::string& EclGenericTracerModel:: +name(int tracerIdx) const +{ + return this->eclState_.tracer()[tracerIdx].name; +} + +template +void EclGenericTracerModel:: +doInit(bool rst, size_t numGridDof, + size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx) +{ + const auto& tracers = eclState_.tracer(); + + if (tracers.size() == 0) + return; // tracer treatment is supposed to be disabled + + // retrieve the number of tracers from the deck + const size_t numTracers = tracers.size(); + tracerConcentration_.resize(numTracers); + storageOfTimeIndex1_.resize(numTracers); + + // the phase where the tracer is + tracerPhaseIdx_.resize(numTracers); + for (size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) { + const auto& tracer = tracers[tracerIdx]; + + if (tracer.phase == Phase::WATER) + tracerPhaseIdx_[tracerIdx] = waterPhaseIdx; + else if (tracer.phase == Phase::OIL) + tracerPhaseIdx_[tracerIdx] = oilPhaseIdx; + else if (tracer.phase == Phase::GAS) + tracerPhaseIdx_[tracerIdx] = gasPhaseIdx; + + tracerConcentration_[tracerIdx].resize(numGridDof); + storageOfTimeIndex1_[tracerIdx].resize(numGridDof); + + if (rst) + continue; + + + //TBLK keyword + if (tracer.free_concentration.has_value()){ + const auto& free_concentration = tracer.free_concentration.value(); + int tblkDatasize = free_concentration.size(); + if (tblkDatasize < cartMapper_.cartesianSize()){ + throw std::runtime_error("Wrong size of TBLK for" + tracer.name); + } + for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){ + int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx); + tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx]; + } + } + //TVDPF keyword + else if (tracer.free_tvdp.has_value()) { + const auto& free_tvdp = tracer.free_tvdp.value(); + for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){ + tracerConcentration_[tracerIdx][globalDofIdx] = + free_tvdp.evaluate("TRACER_CONCENTRATION", + centroids_(globalDofIdx)[2]); + } + } else + throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name)); + } + + // allocate matrix for storing the Jacobian of the tracer residual + tracerMatrix_ = std::make_unique(numGridDof, numGridDof, TracerMatrix::random); + + // find the sparsity pattern of the tracer matrix + using NeighborSet = std::set; + std::vector neighbors(numGridDof); + + Stencil stencil(gridView_, dofMapper_); + for (const auto& elem : elements(gridView_)) { + stencil.update(elem); + + for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) { + unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx); + + for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) { + unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx); + neighbors[myIdx].insert(neighborIdx); + } + } + } + + // allocate space for the rows of the matrix + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) + tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size()); + tracerMatrix_->endrowsizes(); + + // fill the rows with indices. each degree of freedom talks to + // all of its neighbors. (it also talks to itself since + // degrees of freedom are sometimes quite egocentric.) + for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) { + typename NeighborSet::iterator nIt = neighbors[dofIdx].begin(); + typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end(); + for (; nIt != nEndIt; ++nIt) + tracerMatrix_->addindex(dofIdx, *nIt); + } + tracerMatrix_->endindices(); + + const int sizeCartGrid = cartMapper_.cartesianSize(); + cartToGlobal_.resize(sizeCartGrid); + for (unsigned i = 0; i < numGridDof; ++i) { + int cartIdx = cartMapper_.cartesianIndex(i); + cartToGlobal_[cartIdx] = i; + } +} + +template +bool EclGenericTracerModel:: +linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b) +{ +#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7) + Dune::FMatrixPrecision::set_singular_limit(1.e-30); + Dune::FMatrixPrecision::set_absolute_limit(1.e-30); +#endif + x = 0.0; + Scalar tolerance = 1e-2; + int maxIter = 100; + + int verbosity = 0; + PropertyTree prm; + prm.put("maxiter", maxIter); + prm.put("tol", tolerance); + prm.put("verbosity", verbosity); + prm.put("solver", std::string("bicgstab")); + prm.put("preconditioner.type", std::string("ParOverILU0")); + +#if HAVE_MPI + if(gridView_.grid().comm().size() > 1) + { + auto [tracerOperator, solver] = + createParallelFlexibleSolver(gridView_.grid(), M, prm); + (void) tracerOperator; + + Dune::InverseOperatorResult result; + solver->apply(x, b, result); + + // return the result of the solver + return result.converged; + } + else + { +#endif + using TracerSolver = Dune::BiCGSTABSolver; + using TracerOperator = Dune::MatrixAdapter; + using TracerScalarProduct = Dune::SeqScalarProduct; + using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>; + + TracerOperator tracerOperator(M); + TracerScalarProduct tracerScalarProduct; + TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0 + + TracerSolver solver (tracerOperator, tracerScalarProduct, + tracerPreconditioner, tolerance, maxIter, + verbosity); + + Dune::InverseOperatorResult result; + solver.apply(x, b, result); + + // return the result of the solver + return result.converged; +#if HAVE_MPI + } +#endif +} + +template +bool EclGenericTracerModel:: +linearSolveBatchwise_(const TracerMatrix& M, std::vector& x, std::vector& b) +{ +#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7) + Dune::FMatrixPrecision::set_singular_limit(1.e-30); + Dune::FMatrixPrecision::set_absolute_limit(1.e-30); +#endif + Scalar tolerance = 1e-2; + int maxIter = 100; + + int verbosity = 0; + PropertyTree prm; + prm.put("maxiter", maxIter); + prm.put("tol", tolerance); + prm.put("verbosity", verbosity); + prm.put("solver", std::string("bicgstab")); + prm.put("preconditioner.type", std::string("ParOverILU0")); + +#if HAVE_MPI + if(gridView_.grid().comm().size() > 1) + { + auto [tracerOperator, solver] = + createParallelFlexibleSolver(gridView_.grid(), M, prm); + (void) tracerOperator; + bool converged = true; + for (size_t nrhs =0; nrhs < b.size(); ++nrhs) { + x[nrhs] = 0.0; + Dune::InverseOperatorResult result; + solver->apply(x[nrhs], b[nrhs], result); + converged = (converged && result.converged); + } + return converged; + } + else + { +#endif + using TracerSolver = Dune::BiCGSTABSolver; + using TracerOperator = Dune::MatrixAdapter; + using TracerScalarProduct = Dune::SeqScalarProduct; + using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>; + + TracerOperator tracerOperator(M); + TracerScalarProduct tracerScalarProduct; + TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0 + + TracerSolver solver (tracerOperator, tracerScalarProduct, + tracerPreconditioner, tolerance, maxIter, + verbosity); + + bool converged = true; + for (size_t nrhs =0; nrhs < b.size(); ++nrhs) { + x[nrhs] = 0.0; + Dune::InverseOperatorResult result; + solver.apply(x[nrhs], b[nrhs], result); + converged = (converged && result.converged); + } + + // return the result of the solver + return converged; +#if HAVE_MPI + } +#endif +} + + +} // namespace Opm +#endif diff --git a/ebos/equil/equilibrationhelpers.cc b/ebos/equil/equilibrationhelpers.cc index 01c19cff2..6ed0d055e 100644 --- a/ebos/equil/equilibrationhelpers.cc +++ b/ebos/equil/equilibrationhelpers.cc @@ -22,597 +22,11 @@ */ #include "config.h" - -#include - -#include - -#include - -#include -#include -#include - -#include +#include "equilibrationhelpers_impl.hh" namespace Opm { namespace EQUIL { - -using FluidSystemSimple = BlackOilFluidSystem; - -// Adjust oil pressure according to gas saturation and cap pressure -using SatOnlyFluidState = SimpleModularFluidState; - -namespace Miscibility { - -template -RsVD::RsVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rs) - : pvtRegionIdx_(pvtRegionIdx) - , rsVsDepth_(depth, rs) -{ -} - -template -double RsVD:: -operator()(const double depth, - const double press, - const double temp, - const double satGas) const -{ - if (satGas > 0.0) { - return satRs(press, temp); - } - else { - if (rsVsDepth_.xMin() > depth) - return rsVsDepth_.valueAt(0); - else if (rsVsDepth_.xMax() < depth) - return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1); - return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false)); - } -} - -template -double RsVD::satRs(const double press, const double temp) const -{ - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); -} - -template -PBVD::PBVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& pbub) - : pvtRegionIdx_(pvtRegionIdx) - , pbubVsDepth_(depth, pbub) -{ -} - -template -double PBVD:: -operator()(const double depth, - const double cellPress, - const double temp, - const double satGas) const -{ - double press = cellPress; - if (satGas <= 0.0) { - if (pbubVsDepth_.xMin() > depth) - press = pbubVsDepth_.valueAt(0); - else if (pbubVsDepth_.xMax() < depth) - press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1); - else - press = pbubVsDepth_.eval(depth, /*extrapolate=*/false); - } - return satRs(std::min(press, cellPress), temp); -} - -template -double PBVD:: -satRs(const double press, const double temp) const -{ - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); -} - -template -PDVD::PDVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& pdew) - : pvtRegionIdx_(pvtRegionIdx) - , pdewVsDepth_(depth, pdew) -{ -} - -template -double PDVD:: -operator()(const double depth, - const double cellPress, - const double temp, - const double satOil) const -{ - double press = cellPress; - if (satOil <= 0.0) { - if (pdewVsDepth_.xMin() > depth) - press = pdewVsDepth_.valueAt(0); - else if (pdewVsDepth_.xMax() < depth) - press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1); - else - press = pdewVsDepth_.eval(depth, /*extrapolate=*/false); - } - return satRv(std::min(press, cellPress), temp); -} - -template -double PDVD:: -satRv(const double press, const double temp) const -{ - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); -} - - -template -RvVD::RvVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rv) - : pvtRegionIdx_(pvtRegionIdx) - , rvVsDepth_(depth, rv) -{ -} - -template -double RvVD:: -operator()(const double depth, - const double press, - const double temp, - const double satOil) const -{ - if (std::abs(satOil) > 1e-16) { - return satRv(press, temp); - } - else { - if (rvVsDepth_.xMin() > depth) - return rvVsDepth_.valueAt(0); - else if (rvVsDepth_.xMax() < depth) - return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1); - return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false)); - } -} - -template -double RvVD:: -satRv(const double press, const double temp) const -{ - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); -} - - -template -RvwVD::RvwVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rvw) - : pvtRegionIdx_(pvtRegionIdx) - , rvwVsDepth_(depth, rvw) -{ -} - -template -double RvwVD:: -operator()(const double depth, - const double press, - const double temp, - const double satWat) const -{ - if (std::abs(satWat) > 1e-16) { - return satRvw(press, temp); //saturated Rvw - } - else { - if (rvwVsDepth_.xMin() > depth) - return rvwVsDepth_.valueAt(0); - else if (rvwVsDepth_.xMax() < depth) - return rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1); - return std::min(satRvw(press, temp), rvwVsDepth_.eval(depth, /*extrapolate=*/false)); - } -} - -template -double RvwVD:: -satRvw(const double press, const double temp) const -{ - return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press); -} - - -template -RsSatAtContact:: -RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) - : pvtRegionIdx_(pvtRegionIdx) -{ - rsSatContact_ = satRs(pContact, T_contact); -} - -template -double RsSatAtContact:: -operator()(const double /* depth */, - const double press, - const double temp, - const double satGas) const -{ - if (satGas > 0.0) { - return satRs(press, temp); - } - else { - return std::min(satRs(press, temp), rsSatContact_); - } -} - -template -double RsSatAtContact:: -satRs(const double press, const double temp) const -{ - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); -} - -template -RvSatAtContact:: -RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) - : pvtRegionIdx_(pvtRegionIdx) -{ - rvSatContact_ = satRv(pContact, T_contact); -} - -template -double RvSatAtContact:: -operator()(const double /*depth*/, - const double press, - const double temp, - const double satOil) const -{ - if (satOil > 0.0) { - return satRv(press, temp); - } - else { - return std::min(satRv(press, temp), rvSatContact_); - } -} - -template -double RvSatAtContact:: -satRv(const double press, const double temp) const -{ - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; -} - -template -RvwSatAtContact:: -RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) - : pvtRegionIdx_(pvtRegionIdx) -{ - rvwSatContact_ = satRvw(pContact, T_contact); -} - -template -double RvwSatAtContact:: -operator()(const double /*depth*/, - const double press, - const double temp, - const double satWat) const -{ - if (satWat > 0.0) { - return satRvw(press, temp); - } - else { - return std::min(satRvw(press, temp), rvwSatContact_); - } -} - -template -double RvwSatAtContact:: -satRvw(const double press, const double temp) const -{ - return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);; -} - -} // namespace Miscibility - -EquilReg::EquilReg(const EquilRecord& rec, - std::shared_ptr rs, - std::shared_ptr rv, - std::shared_ptr rvw, - const TabulatedFunction& saltVdTable, - const int pvtIdx) - : rec_ (rec) - , rs_ (rs) - , rv_ (rv) - , rvw_ (rvw) - , saltVdTable_ (saltVdTable) - , pvtIdx_ (pvtIdx) -{ -} - -double EquilReg::datum() const -{ - return this->rec_.datumDepth(); -} - -double EquilReg::pressure() const -{ - return this->rec_.datumDepthPressure(); -} - -double EquilReg::zwoc() const -{ - return this->rec_.waterOilContactDepth(); -} - -double EquilReg::pcowWoc() const -{ - return this->rec_.waterOilContactCapillaryPressure(); -} - -double EquilReg::zgoc() const -{ - return this->rec_.gasOilContactDepth(); -} - -double EquilReg::pcgoGoc() const -{ - return this->rec_.gasOilContactCapillaryPressure(); -} - -int EquilReg::equilibrationAccuracy() const -{ - return this->rec_.initializationTargetAccuracy(); -} - -const EquilReg::CalcDissolution& -EquilReg::dissolutionCalculator() const -{ - return *this->rs_; -} - -const EquilReg::CalcEvaporation& -EquilReg::evaporationCalculator() const -{ - return *this->rv_; -} - -const EquilReg::CalcWaterEvaporation& -EquilReg::waterEvaporationCalculator() const -{ - return *this->rvw_; -} - -const EquilReg::TabulatedFunction& -EquilReg::saltVdTable() const -{ - return saltVdTable_; -} - -int EquilReg::pvtIdx() const -{ - return this->pvtIdx_; -} - -template -PcEq:: -PcEq(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double targetPc) - : materialLawManager_(materialLawManager), - phase_(phase), - cell_(cell), - targetPc_(targetPc) -{ -} - -template -double PcEq:: -operator()(double s) const -{ - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase_, s); - - std::array pc{0.0}; - using MaterialLaw = typename MaterialLawManager::MaterialLaw; - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; - return pcPhase - targetPc_; -} - -template -PcEqSum:: -PcEqSum(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double targetPc) - : materialLawManager_(materialLawManager), - phase1_(phase1), - phase2_(phase2), - cell_(cell), - targetPc_(targetPc) -{ -} - -template -double PcEqSum:: -operator()(double s) const -{ - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase1_, s); - fluidState.setSaturation(phase2_, 1.0 - s); - - std::array pc {0.0}; - - using MaterialLaw = typename MaterialLawManager::MaterialLaw; - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_]; - double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_]; - return pc1 + pc2 - targetPc_; -} - -template -double minSaturations(const MaterialLawManager& materialLawManager, - const int phase, const int cell) -{ - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx: - return scaledDrainageInfo.Swl; - - case FluidSystem::gasPhaseIdx: - return scaledDrainageInfo.Sgl; - - case FluidSystem::oilPhaseIdx: - throw std::runtime_error("Min saturation not implemented for oil phase."); - - default: - throw std::runtime_error("Unknown phaseIdx ."); - } - return -1.0; -} - -template -double maxSaturations(const MaterialLawManager& materialLawManager, - const int phase, const int cell) -{ - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx: - return scaledDrainageInfo.Swu; - - case FluidSystem::gasPhaseIdx: - return scaledDrainageInfo.Sgu; - - case FluidSystem::oilPhaseIdx: - throw std::runtime_error("Max saturation not implemented for oil phase."); - - default: - throw std::runtime_error("Unknown phaseIdx ."); - } - return -1.0; -} - -template -double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double targetPc, - const bool increasing) -{ - // Find minimum and maximum saturations. - double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - - // Create the equation f(s) = pc(s) - targetPc - const PcEq f(materialLawManager, phase, cell, targetPc); - double f0 = f(s0); - double f1 = f(s1); - if (!std::isfinite(f0 + f1)) - throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1)); - - if (f0 <= 0.0) - return s0; - else if (f1 >= 0.0) - return s1; - - const double tol = 1e-10; - // should at least converge in 2 times bisection but some safety here: - const int maxIter = -2*static_cast(std::log2(tol)) + 10; - int usedIterations = -1; - const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); - return root; -} - -template -double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double targetPc) -{ - // Find minimum and maximum saturations. - double s0 = minSaturations(materialLawManager, phase1, cell); - double s1 = maxSaturations(materialLawManager, phase1, cell); - - // Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc - const PcEqSum f(materialLawManager, phase1, phase2, cell, targetPc); - double f0 = f(s0); - double f1 = f(s1); - if (f0 <= 0.0) - return s0; - else if (f1 >= 0.0) - return s1; - - assert(f0 > 0.0 && f1 < 0.0); - const double tol = 1e-10; - // should at least converge in 2 times bisection but some safety here: - const int maxIter = -2*static_cast(std::log2(tol)) + 10; - int usedIterations = -1; - const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); - return root; -} - -template -double satFromDepth(const MaterialLawManager& materialLawManager, - const double cellDepth, - const double contactDepth, - const int phase, - const int cell, - const bool increasing) -{ - const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - - if (cellDepth < contactDepth) { - return s0; - } - else { - return s1; - } -} - -template -bool isConstPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell) -{ - // Create the equation f(s) = pc(s); - const PcEq f(materialLawManager, phase, cell, 0); - const double f0 = f(minSaturations(materialLawManager, phase, cell)); - const double f1 = f(maxSaturations(materialLawManager, phase, cell)); - return std::abs(f0 - f1) < std::numeric_limits::epsilon(); -} - + using MatLaw = EclMaterialLawManager>; using FS = BlackOilFluidSystem; diff --git a/ebos/equil/equilibrationhelpers_impl.hh b/ebos/equil/equilibrationhelpers_impl.hh new file mode 100644 index 000000000..a242b9b87 --- /dev/null +++ b/ebos/equil/equilibrationhelpers_impl.hh @@ -0,0 +1,615 @@ +// -*- 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. +*/ + +#include + +#include + +#include + +#include +#include +#include + +#include + +namespace Opm { +namespace EQUIL { + +using FluidSystemSimple = BlackOilFluidSystem; + +// Adjust oil pressure according to gas saturation and cap pressure +using SatOnlyFluidState = SimpleModularFluidState; + +namespace Miscibility { + +template +RsVD::RsVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rs) + : pvtRegionIdx_(pvtRegionIdx) + , rsVsDepth_(depth, rs) +{ +} + +template +double RsVD:: +operator()(const double depth, + const double press, + const double temp, + const double satGas) const +{ + if (satGas > 0.0) { + return satRs(press, temp); + } + else { + if (rsVsDepth_.xMin() > depth) + return rsVsDepth_.valueAt(0); + else if (rsVsDepth_.xMax() < depth) + return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1); + return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false)); + } +} + +template +double RsVD::satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +PBVD::PBVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& pbub) + : pvtRegionIdx_(pvtRegionIdx) + , pbubVsDepth_(depth, pbub) +{ +} + +template +double PBVD:: +operator()(const double depth, + const double cellPress, + const double temp, + const double satGas) const +{ + double press = cellPress; + if (satGas <= 0.0) { + if (pbubVsDepth_.xMin() > depth) + press = pbubVsDepth_.valueAt(0); + else if (pbubVsDepth_.xMax() < depth) + press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1); + else + press = pbubVsDepth_.eval(depth, /*extrapolate=*/false); + } + return satRs(std::min(press, cellPress), temp); +} + +template +double PBVD:: +satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +PDVD::PDVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& pdew) + : pvtRegionIdx_(pvtRegionIdx) + , pdewVsDepth_(depth, pdew) +{ +} + +template +double PDVD:: +operator()(const double depth, + const double cellPress, + const double temp, + const double satOil) const +{ + double press = cellPress; + if (satOil <= 0.0) { + if (pdewVsDepth_.xMin() > depth) + press = pdewVsDepth_.valueAt(0); + else if (pdewVsDepth_.xMax() < depth) + press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1); + else + press = pdewVsDepth_.eval(depth, /*extrapolate=*/false); + } + return satRv(std::min(press, cellPress), temp); +} + +template +double PDVD:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); +} + + +template +RvVD::RvVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rv) + : pvtRegionIdx_(pvtRegionIdx) + , rvVsDepth_(depth, rv) +{ +} + +template +double RvVD:: +operator()(const double depth, + const double press, + const double temp, + const double satOil) const +{ + if (std::abs(satOil) > 1e-16) { + return satRv(press, temp); + } + else { + if (rvVsDepth_.xMin() > depth) + return rvVsDepth_.valueAt(0); + else if (rvVsDepth_.xMax() < depth) + return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1); + return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false)); + } +} + +template +double RvVD:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); +} + + +template +RvwVD::RvwVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rvw) + : pvtRegionIdx_(pvtRegionIdx) + , rvwVsDepth_(depth, rvw) +{ +} + +template +double RvwVD:: +operator()(const double depth, + const double press, + const double temp, + const double satWat) const +{ + if (std::abs(satWat) > 1e-16) { + return satRvw(press, temp); //saturated Rvw + } + else { + if (rvwVsDepth_.xMin() > depth) + return rvwVsDepth_.valueAt(0); + else if (rvwVsDepth_.xMax() < depth) + return rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1); + return std::min(satRvw(press, temp), rvwVsDepth_.eval(depth, /*extrapolate=*/false)); + } +} + +template +double RvwVD:: +satRvw(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press); +} + + +template +RsSatAtContact:: +RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rsSatContact_ = satRs(pContact, T_contact); +} + +template +double RsSatAtContact:: +operator()(const double /* depth */, + const double press, + const double temp, + const double satGas) const +{ + if (satGas > 0.0) { + return satRs(press, temp); + } + else { + return std::min(satRs(press, temp), rsSatContact_); + } +} + +template +double RsSatAtContact:: +satRs(const double press, const double temp) const +{ + return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); +} + +template +RvSatAtContact:: +RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rvSatContact_ = satRv(pContact, T_contact); +} + +template +double RvSatAtContact:: +operator()(const double /*depth*/, + const double press, + const double temp, + const double satOil) const +{ + if (satOil > 0.0) { + return satRv(press, temp); + } + else { + return std::min(satRv(press, temp), rvSatContact_); + } +} + +template +double RvSatAtContact:: +satRv(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; +} + +template +RvwSatAtContact:: +RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rvwSatContact_ = satRvw(pContact, T_contact); +} + +template +double RvwSatAtContact:: +operator()(const double /*depth*/, + const double press, + const double temp, + const double satWat) const +{ + if (satWat > 0.0) { + return satRvw(press, temp); + } + else { + return std::min(satRvw(press, temp), rvwSatContact_); + } +} + +template +double RvwSatAtContact:: +satRvw(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);; +} + +} // namespace Miscibility + +EquilReg::EquilReg(const EquilRecord& rec, + std::shared_ptr rs, + std::shared_ptr rv, + std::shared_ptr rvw, + const TabulatedFunction& saltVdTable, + const int pvtIdx) + : rec_ (rec) + , rs_ (rs) + , rv_ (rv) + , rvw_ (rvw) + , saltVdTable_ (saltVdTable) + , pvtIdx_ (pvtIdx) +{ +} + +double EquilReg::datum() const +{ + return this->rec_.datumDepth(); +} + +double EquilReg::pressure() const +{ + return this->rec_.datumDepthPressure(); +} + +double EquilReg::zwoc() const +{ + return this->rec_.waterOilContactDepth(); +} + +double EquilReg::pcowWoc() const +{ + return this->rec_.waterOilContactCapillaryPressure(); +} + +double EquilReg::zgoc() const +{ + return this->rec_.gasOilContactDepth(); +} + +double EquilReg::pcgoGoc() const +{ + return this->rec_.gasOilContactCapillaryPressure(); +} + +int EquilReg::equilibrationAccuracy() const +{ + return this->rec_.initializationTargetAccuracy(); +} + +const EquilReg::CalcDissolution& +EquilReg::dissolutionCalculator() const +{ + return *this->rs_; +} + +const EquilReg::CalcEvaporation& +EquilReg::evaporationCalculator() const +{ + return *this->rv_; +} + +const EquilReg::CalcWaterEvaporation& +EquilReg::waterEvaporationCalculator() const +{ + return *this->rvw_; +} + +const EquilReg::TabulatedFunction& +EquilReg::saltVdTable() const +{ + return saltVdTable_; +} + +int EquilReg::pvtIdx() const +{ + return this->pvtIdx_; +} + +template +PcEq:: +PcEq(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double targetPc) + : materialLawManager_(materialLawManager), + phase_(phase), + cell_(cell), + targetPc_(targetPc) +{ +} + +template +double PcEq:: +operator()(double s) const +{ + const auto& matParams = materialLawManager_.materialLawParams(cell_); + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); + fluidState.setSaturation(phase_, s); + + std::array pc{0.0}; + using MaterialLaw = typename MaterialLawManager::MaterialLaw; + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; + double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; + return pcPhase - targetPc_; +} + +template +PcEqSum:: +PcEqSum(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double targetPc) + : materialLawManager_(materialLawManager), + phase1_(phase1), + phase2_(phase2), + cell_(cell), + targetPc_(targetPc) +{ +} + +template +double PcEqSum:: +operator()(double s) const +{ + const auto& matParams = materialLawManager_.materialLawParams(cell_); + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); + fluidState.setSaturation(phase1_, s); + fluidState.setSaturation(phase2_, 1.0 - s); + + std::array pc {0.0}; + + using MaterialLaw = typename MaterialLawManager::MaterialLaw; + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; + double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_]; + double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; + double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_]; + return pc1 + pc2 - targetPc_; +} + +template +double minSaturations(const MaterialLawManager& materialLawManager, + const int phase, const int cell) +{ + const auto& scaledDrainageInfo = + materialLawManager.oilWaterScaledEpsInfoDrainage(cell); + + // Find minimum and maximum saturations. + switch(phase) { + case FluidSystem::waterPhaseIdx: + return scaledDrainageInfo.Swl; + + case FluidSystem::gasPhaseIdx: + return scaledDrainageInfo.Sgl; + + case FluidSystem::oilPhaseIdx: + throw std::runtime_error("Min saturation not implemented for oil phase."); + + default: + throw std::runtime_error("Unknown phaseIdx ."); + } + return -1.0; +} + +template +double maxSaturations(const MaterialLawManager& materialLawManager, + const int phase, const int cell) +{ + const auto& scaledDrainageInfo = + materialLawManager.oilWaterScaledEpsInfoDrainage(cell); + + // Find minimum and maximum saturations. + switch(phase) { + case FluidSystem::waterPhaseIdx: + return scaledDrainageInfo.Swu; + + case FluidSystem::gasPhaseIdx: + return scaledDrainageInfo.Sgu; + + case FluidSystem::oilPhaseIdx: + throw std::runtime_error("Max saturation not implemented for oil phase."); + + default: + throw std::runtime_error("Unknown phaseIdx ."); + } + return -1.0; +} + +template +double satFromPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double targetPc, + const bool increasing) +{ + // Find minimum and maximum saturations. + double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); + double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); + + // Create the equation f(s) = pc(s) - targetPc + const PcEq f(materialLawManager, phase, cell, targetPc); + double f0 = f(s0); + double f1 = f(s1); + if (!std::isfinite(f0 + f1)) + throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1)); + + if (f0 <= 0.0) + return s0; + else if (f1 >= 0.0) + return s1; + + const double tol = 1e-10; + // should at least converge in 2 times bisection but some safety here: + const int maxIter = -2*static_cast(std::log2(tol)) + 10; + int usedIterations = -1; + const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); + return root; +} + +template +double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double targetPc) +{ + // Find minimum and maximum saturations. + double s0 = minSaturations(materialLawManager, phase1, cell); + double s1 = maxSaturations(materialLawManager, phase1, cell); + + // Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc + const PcEqSum f(materialLawManager, phase1, phase2, cell, targetPc); + double f0 = f(s0); + double f1 = f(s1); + if (f0 <= 0.0) + return s0; + else if (f1 >= 0.0) + return s1; + + assert(f0 > 0.0 && f1 < 0.0); + const double tol = 1e-10; + // should at least converge in 2 times bisection but some safety here: + const int maxIter = -2*static_cast(std::log2(tol)) + 10; + int usedIterations = -1; + const double root = RegulaFalsiBisection::solve(f, s0, s1, maxIter, tol, usedIterations); + return root; +} + +template +double satFromDepth(const MaterialLawManager& materialLawManager, + const double cellDepth, + const double contactDepth, + const int phase, + const int cell, + const bool increasing) +{ + const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); + const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); + + if (cellDepth < contactDepth) { + return s0; + } + else { + return s1; + } +} + +template +bool isConstPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell) +{ + // Create the equation f(s) = pc(s); + const PcEq f(materialLawManager, phase, cell, 0); + const double f0 = f(minSaturations(materialLawManager, phase, cell)); + const double f1 = f(maxSaturations(materialLawManager, phase, cell)); + return std::abs(f0 - f1) < std::numeric_limits::epsilon(); +} + +} +} diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc index 51935c264..0eb5d1091 100644 --- a/ebos/equil/initstateequil.cc +++ b/ebos/equil/initstateequil.cc @@ -22,1927 +22,12 @@ */ #include "config.h" - -#include -#include - -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include - -#include -#include -#include -#include - -#if HAVE_DUNE_FEM -#include -#include -#include -#endif - -#if HAVE_DUNE_ALUGRID -#include -#include -#endif // HAVE_DUNE_ALUGRID +#include "initstateequil.hh" +#include "initstateequil_impl.hh" namespace Opm { namespace EQUIL { - -namespace Details { - -template -void verticalExtent(const CellRange& cells, - const std::vector>& cellZMinMax, - const Comm& comm, - std::array& span) -{ - span[0] = std::numeric_limits::max(); - span[1] = std::numeric_limits::lowest(); - - // Define vertical span as - // - // [minimum(node depth(cells)), maximum(node depth(cells))] - // - // Note: The implementation of 'RK4IVP<>' implicitly - // imposes the requirement that cell centroids are all - // within this vertical span. That requirement is not - // checked. - for (const auto& cell : cells) { - if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; } - if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; } - } - span[0] = comm.min(span[0]); - span[1] = comm.max(span[1]); -} - -void subdivisionCentrePoints(const double left, - const double right, - const int numIntervals, - std::vector>& subdiv) -{ - const auto h = (right - left) / numIntervals; - - auto end = left; - for (auto i = 0*numIntervals; i < numIntervals; ++i) { - const auto start = end; - end = left + (i + 1)*h; - - subdiv.emplace_back((start + end) / 2, h); - } -} - -template -std::vector> -horizontalSubdivision(const CellID cell, - const std::pair topbot, - const int numIntervals) -{ - auto subdiv = std::vector>{}; - subdiv.reserve(2 * numIntervals); - - if (topbot.first > topbot.second) { - throw std::out_of_range { - "Negative thickness (inverted top/bottom faces) in cell " - + std::to_string(cell) - }; - } - - subdivisionCentrePoints(topbot.first, topbot.second, - 2*numIntervals, subdiv); - - return subdiv; -} - -template -double cellCenterDepth(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int zCoord = Element::dimension - 1; - double zz = 0.0; - - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - for (int i=0; i < corners; ++i) - zz += geometry.corner(i)[zCoord]; - - return zz/corners; -} - -template -std::pair cellZSpan(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int zCoord = Element::dimension - 1; - double bot = 0.0; - double top = 0.0; - - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - assert(corners == 8); - for (int i=0; i < 4; ++i) - bot += geometry.corner(i)[zCoord]; - for (int i=4; i < corners; ++i) - top += geometry.corner(i)[zCoord]; - - return std::make_pair(bot/4, top/4); -} - -template -std::pair cellZMinMax(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int zCoord = Element::dimension - 1; - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - assert(corners == 8); - auto min = std::numeric_limits::max(); - auto max = std::numeric_limits::lowest(); - - - for (int i=0; i < corners; ++i) { - min = std::min(min, geometry.corner(i)[zCoord]); - max = std::max(max, geometry.corner(i)[zCoord]); - } - return std::make_pair(min, max); -} - -template -RK4IVP::RK4IVP(const RHS& f, - const std::array& span, - const double y0, - const int N) - : N_(N) - , span_(span) -{ - const double h = stepsize(); - const double h2 = h / 2; - const double 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 double x = span_[0] + i*h; - const double y = y_.back(); - - const double k1 = f_[i]; - const double k2 = f(x + h2, y + h2*k1); - const double k3 = f(x + h2, y + h2*k2); - const double 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() == std::vector::size_type(N + 1)); -} - -template -double RK4IVP:: -operator()(const double x) const -{ - // Dense output (O(h**3)) according to Shampine - // (Hermite interpolation) - const double h = stepsize(); - int i = (x - span_[0]) / h; - const double t = (x - (span_[0] + i*h)) / h; - - // Crude handling of evaluation point outside "span_"; - if (i < 0) { i = 0; } - if (N_ <= i) { i = N_ - 1; } - - const double y0 = y_[i], y1 = y_[i + 1]; - const double f0 = f_[i], f1 = f_[i + 1]; - - double 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 -double RK4IVP:: -stepsize() const -{ - return (span_[1] - span_[0]) / N_; -} - -namespace PhasePressODE { - -template -Water:: -Water(const double temp, - const TabulatedFunction& saltVdTable, - const int pvtRegionIdx, - const double normGrav) - : temp_(temp) - , saltVdTable_(saltVdTable) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -double Water:: -operator()(const double depth, - const double press) const -{ - return this->density(depth, press) * g_; -} - -template -double Water:: -density(const double depth, - const double press) const -{ - // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate - double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); - double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0 /*=Rsw*/, saltConcentration); - rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); - return rho; -} - -template -Oil:: -Oil(const double temp, - const RS& rs, - const int pvtRegionIdx, - const double normGrav) - : temp_(temp) - , rs_(rs) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -double Oil:: -operator()(const double depth, - const double press) const -{ - return this->density(depth, press) * g_; -} - -template -double Oil:: -density(const double depth, - const double press) const -{ - double rs = 0.0; - if(FluidSystem::enableDissolvedGas()) - rs = rs_(depth, press, temp_); - - double bOil = 0.0; - if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press)) { - bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } - else { - bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); - } - double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableDissolvedGas()) { - rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - } - - return rho; -} - -template -Gas:: -Gas(const double temp, - const RV& rv, - const RVW& rvw, - const int pvtRegionIdx, - const double normGrav) - : temp_(temp) - , rv_(rv) - , rvw_(rvw) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -double Gas:: -operator()(const double depth, - const double press) const -{ - return this->density(depth, press) * g_; -} - -template -double Gas:: -density(const double depth, - const double press) const -{ - double rv = 0.0; - if (FluidSystem::enableVaporizedOil()) - rv = rv_(depth, press, temp_); - - double rvw = 0.0; - if (FluidSystem::enableVaporizedWater()) - rvw = rvw_(depth, press, temp_); - - double bGas = 0.0; - - if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { - if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) - && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) - { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, rvw); - } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_) - + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); - return rho; - } - - if (FluidSystem::enableVaporizedOil()){ - if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=rvw*/); - } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - return rho; - } - - if (FluidSystem::enableVaporizedWater()){ - if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } - else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, rvw); - } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); - return rho; - } - - // immiscible gas - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, 0.0/*=rvw*/); - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - - return rho; -} - -} - -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 -double -PressureTable:: -PressureFunction:: -value(const double 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:: -checkPtr(const PressFunc* phasePress, - const std::string& phaseName) const -{ - if (phasePress != nullptr) { return; } - - throw std::invalid_argument { - "Phase pressure function for \"" + phaseName - + "\" most not be null" - }; -} - -template -typename PressureTable::Strategy -PressureTable:: -selectEquilibrationStrategy(const Region& reg) const -{ - if (!this->oilActive()) { - if (reg.datum() > reg.zwoc()) { // Datum in water zone - return &PressureTable::equil_WOG; - } - return &PressureTable::equil_GOW; - } - - if (reg.datum() > reg.zwoc()) { // Datum in water zone - return &PressureTable::equil_WOG; - } - else if (reg.datum() < reg.zgoc()) { // Datum in gas zone - return &PressureTable::equil_GOW; - } - else { // Datum in oil zone - return &PressureTable::equil_OWG; - } -} - -template -void PressureTable:: -copyInPointers(const PressureTable& rhs) -{ - if (rhs.oil_ != nullptr) { - this->oil_ = std::make_unique(*rhs.oil_); - } - - if (rhs.gas_ != nullptr) { - this->gas_ = std::make_unique(*rhs.gas_); - } - - if (rhs.wat_ != nullptr) { - this->wat_ = std::make_unique(*rhs.wat_); - } -} - -template -PhaseSaturations:: -PhaseSaturations(MaterialLawManager& matLawMgr, - const std::vector& swatInit) - : matLawMgr_(matLawMgr) - , swatInit_ (swatInit) -{ -} - -template -PhaseSaturations:: -PhaseSaturations(const PhaseSaturations& rhs) - : matLawMgr_(rhs.matLawMgr_) - , swatInit_ (rhs.swatInit_) - , sat_ (rhs.sat_) - , press_ (rhs.press_) -{ - // Note: We don't need to do anything to the 'fluidState_' here. - this->setEvaluationPoint(*rhs.evalPt_.position, - *rhs.evalPt_.region, - *rhs.evalPt_.ptable); -} - -template -const PhaseQuantityValue& -PhaseSaturations:: -deriveSaturations(const Position& x, - const Region& reg, - const PTable& ptable) -{ - this->setEvaluationPoint(x, reg, ptable); - this->initializePhaseQuantities(); - - if (ptable.waterActive()) { this->deriveWaterSat(); } - if (ptable.gasActive()) { this->deriveGasSat(); } - - if (this->isOverlappingTransition()) { - this->fixUnphysicalTransition(); - } - - if (ptable.oilActive()) { this->deriveOilSat(); } - - this->accountForScaledSaturations(); - - return this->sat_; -} - -template -void PhaseSaturations:: -setEvaluationPoint(const Position& x, - const Region& reg, - const PTable& ptable) -{ - this->evalPt_.position = &x; - this->evalPt_.region = ® - this->evalPt_.ptable = &ptable; -} - -template -void PhaseSaturations:: -initializePhaseQuantities() -{ - this->sat_.reset(); - this->press_.reset(); - - const auto depth = this->evalPt_.position->depth; - const auto& ptable = *this->evalPt_.ptable; - - if (ptable.oilActive()) { - this->press_.oil = ptable.oil(depth); - } - - if (ptable.gasActive()) { - this->press_.gas = ptable.gas(depth); - } - - if (ptable.waterActive()) { - this->press_.water = ptable.water(depth); - } -} - -template -void PhaseSaturations::deriveOilSat() -{ - this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas; -} - -template -void PhaseSaturations::deriveGasSat() -{ - auto& sg = this->sat_.gas; - - const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg. - const auto oilActive = this->evalPt_.ptable->oilActive(); - - if (this->isConstCapPress(this->gasPos())) { - // Sharp interface between phases. Can derive phase saturation - // directly from knowing where 'depth' of evaluation point is - // relative to depth of O/G contact. - const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc(); - sg = this->fromDepthTable(gas_contact, - this->gasPos(), isIncr); - } - else { - // Capillary pressure curve is non-constant, meaning there is a - // transition zone between the gas and oil phases. Invert capillary - // pressure relation - // - // Pcgo(Sg) = Pg - Po - // - // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg). - const auto pw = oilActive? this->press_.oil : this->press_.water; - const auto pcgo = this->press_.gas - pw; - sg = this->invertCapPress(pcgo, this->gasPos(), isIncr); - } -} - -template -void PhaseSaturations::deriveWaterSat() -{ - auto& sw = this->sat_.water; - - const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. - - if (this->isConstCapPress(this->waterPos())) { - // Sharp interface between phases. Can derive phase saturation - // directly from knowing where 'depth' of evaluation point is - // relative to depth of O/W contact. - sw = this->fromDepthTable(this->evalPt_.region->zwoc(), - this->waterPos(), isIncr); - } - else { - // Capillary pressure curve is non-constant, meaning there is a - // transition zone between the oil and water phases. Invert - // capillary pressure relation - // - // Pcow(Sw) = Po - Pw - // - // unless the model uses "SWATINIT". In the latter case, pick the - // saturation directly from the SWATINIT array of the pertinent - // cell. - const auto pcow = this->press_.oil - this->press_.water; - - sw = this->swatInit_.empty() - ? this->invertCapPress(pcow, this->waterPos(), isIncr) - : this->applySwatInit(pcow); - } -} - -template -void PhaseSaturations:: -fixUnphysicalTransition() -{ - auto& sg = this->sat_.gas; - auto& sw = this->sat_.water; - - // Overlapping gas/oil and oil/water transition zones can lead to - // unphysical phase saturations when individual saturations are derived - // directly from inverting O/G and O/W capillary pressure curves. - // - // Recalculate phase saturations using the implied gas/water capillary - // pressure: Pg - Pw. - const auto pcgw = this->press_.gas - this->press_.water; - if (! this->swatInit_.empty()) { - // Re-scale Pc to reflect imposed sw for vanishing oil phase. This - // seems consistent with ECLIPSE, but fails to honour SWATINIT in - // case of non-trivial gas/oil capillary pressure. - sw = this->applySwatInit(pcgw, sw); - } - - sw = satFromSumOfPcs - (this->matLawMgr_, this->waterPos(), this->gasPos(), - this->evalPt_.position->cell, pcgw); - sg = 1.0 - sw; - - this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg); - this->fluidState_.setSaturation(this->gasPos(), sg); - this->fluidState_.setSaturation(this->waterPos(), this->evalPt_ - .ptable->waterActive() ? sw : 0.0); - - // Pcgo = Pg - Po => Po = Pg - Pcgo - this->computeMaterialLawCapPress(); - this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); -} - -template -void PhaseSaturations:: -accountForScaledSaturations() -{ - const auto gasActive = this->evalPt_.ptable->gasActive(); - const auto watActive = this->evalPt_.ptable->waterActive(); - const auto oilActive = this->evalPt_.ptable->oilActive(); - - auto sg = gasActive? this->sat_.gas : 0.0; - auto sw = watActive? this->sat_.water : 0.0; - auto so = oilActive? this->sat_.oil : 0.0; - - this->fluidState_.setSaturation(this->waterPos(), sw); - this->fluidState_.setSaturation(this->oilPos(), so); - this->fluidState_.setSaturation(this->gasPos(), sg); - - const auto& scaledDrainageInfo = this->matLawMgr_ - .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell); - - const auto thresholdSat = 1.0e-6; - if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) { - // Water saturation exceeds maximum possible value. Reset oil phase - // pressure to that which corresponds to maximum possible water - // saturation value. - this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu); - } else if (gasActive) { - this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu); - } - sw = scaledDrainageInfo.Swu; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcow = Po - Pw => Po = Pw + Pcow - this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); - } else { - // Pcgw = Pg - Pw => Pg = Pw + Pcgw - this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); - } - - } - if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { - // Gas saturation exceeds maximum possible value. Reset oil phase - // pressure to that which corresponds to maximum possible gas - // saturation value. - this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu); - } else if (watActive) { - this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu); - } - sg = scaledDrainageInfo.Sgu; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcgo = Pg - Po => Po = Pg - Pcgo - this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); - } else { - // Pcgw = Pg - Pw => Pw = Pg - Pcgw - this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); - } - } - - if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) { - // Water saturation less than minimum possible value in cell. Reset - // water phase pressure to that which corresponds to minimum - // possible water saturation value. - this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl); - } else if (gasActive) { - this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl); - } - sw = scaledDrainageInfo.Swl; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcwo = Po - Pw => Pw = Po - Pcow - this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); - } else { - // Pcgw = Pg - Pw => Pw = Pg - Pcgw - this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); - } - } - - if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { - // Gas saturation less than minimum possible value in cell. Reset - // gas phase pressure to that which corresponds to minimum possible - // gas saturation. - this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl); - } else if (watActive) { - this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl); - } - sg = scaledDrainageInfo.Sgl; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); - } else { - // Pcgw = Pg - Pw => Pg = Pw + Pcgw - this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); - } - } -} - -template -double PhaseSaturations:: -applySwatInit(const double pcow) -{ - return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); -} - -template -double PhaseSaturations:: -applySwatInit(const double pcow, const double sw) -{ - return this->matLawMgr_ - .applySwatinit(this->evalPt_.position->cell, pcow, sw); -} - -template -void PhaseSaturations:: -computeMaterialLawCapPress() -{ - const auto& matParams = this->matLawMgr_ - .materialLawParams(this->evalPt_.position->cell); - - this->matLawCapPress_.fill(0.0); - MaterialLaw::capillaryPressures(this->matLawCapPress_, - matParams, this->fluidState_); -} - -template -double PhaseSaturations:: -materialLawCapPressGasOil() const -{ - return this->matLawCapPress_[this->oilPos()] - + this->matLawCapPress_[this->gasPos()]; -} - -template -double PhaseSaturations:: -materialLawCapPressOilWater() const -{ - return this->matLawCapPress_[this->oilPos()] - - this->matLawCapPress_[this->waterPos()]; -} - -template -double PhaseSaturations:: -materialLawCapPressGasWater() const -{ - return this->matLawCapPress_[this->gasPos()] - - this->matLawCapPress_[this->waterPos()]; -} - -template -bool PhaseSaturations:: -isConstCapPress(const PhaseIdx phaseIdx) const -{ - return isConstPc - (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell); -} - -template -bool PhaseSaturations:: -isOverlappingTransition() const -{ - return this->evalPt_.ptable->gasActive() - && this->evalPt_.ptable->waterActive() - && ((this->sat_.gas + this->sat_.water) > 1.0); -} - -template -double PhaseSaturations:: -fromDepthTable(const double contactdepth, - const PhaseIdx phasePos, - const bool isincr) const -{ - return satFromDepth - (this->matLawMgr_, this->evalPt_.position->depth, - contactdepth, static_cast(phasePos), - this->evalPt_.position->cell, isincr); -} - -template -double PhaseSaturations:: -invertCapPress(const double pc, - const PhaseIdx phasePos, - const bool isincr) const -{ - return satFromPc - (this->matLawMgr_, static_cast(phasePos), - this->evalPt_.position->cell, pc, isincr); -} - -template -PressureTable:: -PressureTable(const double gravity, - const int samplePoints) - : gravity_(gravity) - , nsample_(samplePoints) -{ -} - -template -PressureTable:: -PressureTable(const PressureTable& rhs) - : gravity_(rhs.gravity_) - , nsample_(rhs.nsample_) -{ - this->copyInPointers(rhs); -} - -template -PressureTable:: -PressureTable(PressureTable&& rhs) - : gravity_(rhs.gravity_) - , nsample_(rhs.nsample_) - , oil_ (std::move(rhs.oil_)) - , gas_ (std::move(rhs.gas_)) - , wat_ (std::move(rhs.wat_)) -{ -} - -template -PressureTable& -PressureTable:: -operator=(const PressureTable& rhs) -{ - this->gravity_ = rhs.gravity_; - this->nsample_ = rhs.nsample_; - this->copyInPointers(rhs); - - return *this; -} - -template -PressureTable& -PressureTable:: -operator=(PressureTable&& rhs) -{ - this->gravity_ = rhs.gravity_; - this->nsample_ = rhs.nsample_; - - this->oil_ = std::move(rhs.oil_); - this->gas_ = std::move(rhs.gas_); - this->wat_ = std::move(rhs.wat_); - - return *this; -} - -template -void PressureTable:: -equilibrate(const Region& reg, - const VSpan& span) -{ - // One of the PressureTable::equil_*() member functions. - auto equil = this->selectEquilibrationStrategy(reg); - - (this->*equil)(reg, span); -} - -template -bool PressureTable:: -oilActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); -} - -template -bool PressureTable:: -gasActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); -} - -template -bool PressureTable:: -waterActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); -} - -template -double PressureTable:: -oil(const double depth) const -{ - this->checkPtr(this->oil_.get(), "OIL"); - - return this->oil_->value(depth); -} - -template -double PressureTable:: -gas(const double depth) const -{ - this->checkPtr(this->gas_.get(), "GAS"); - - return this->gas_->value(depth); -} - - -template -double PressureTable:: -water(const double depth) const -{ - this->checkPtr(this->wat_.get(), "WATER"); - - return this->wat_->value(depth); -} - -template -void PressureTable:: -equil_WOG(const Region& reg, const VSpan& span) -{ - // Datum depth in water zone. Calculate phase pressure for water first, - // followed by oil and gas if applicable. - - if (! this->waterActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "WATER zone in model without active water phase" - }; - } - - { - const auto ic = typename WPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeWatPressure(ic, reg, span); - } - - if (this->oilActive()) { - // Pcow = Po - Pw => Po = Pw + Pcow - const auto ic = typename OPress::InitCond { - reg.zwoc(), - this->water(reg.zwoc()) + reg.pcowWoc() - }; - - this->makeOilPressure(ic, reg, span); - } - - if (this->gasActive() && this->oilActive()) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - const auto ic = typename GPress::InitCond { - reg.zgoc(), - this->oil(reg.zgoc()) + reg.pcgoGoc() - }; - - this->makeGasPressure(ic, reg, span); - } else if (this->gasActive() && !this->oilActive()) { - // No oil phase set Pg = Pw + Pcgw - const auto ic = typename GPress::InitCond { - reg.zwoc(), // The WOC is really the GWC for gas/water cases - this->water(reg.zwoc()) + reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases - }; - this->makeGasPressure(ic, reg, span); - } -} - -template -void PressureTable:: -equil_GOW(const Region& reg, const VSpan& span) -{ - // Datum depth in gas zone. Calculate phase pressure for gas first, - // followed by oil and water if applicable. - - if (! this->gasActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "GAS zone in model without active gas phase" - }; - } - - { - const auto ic = typename GPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeGasPressure(ic, reg, span); - } - - if (this->oilActive()) { - // Pcgo = Pg - Po => Po = Pg - Pcgo - const auto ic = typename OPress::InitCond { - reg.zgoc(), - this->gas(reg.zgoc()) - reg.pcgoGoc() - }; - this->makeOilPressure(ic, reg, span); - } - - if (this->waterActive() && this->oilActive()) { - // Pcow = Po - Pw => Pw = Po - Pcow - const auto ic = typename WPress::InitCond { - reg.zwoc(), - this->oil(reg.zwoc()) - reg.pcowWoc() - }; - - this->makeWatPressure(ic, reg, span); - } else if (this->waterActive() && !this->oilActive()) { - // No oil phase set Pw = Pg - Pcgw - const auto ic = typename WPress::InitCond { - reg.zwoc(), // The WOC is really the GWC for gas/water cases - this->gas(reg.zwoc()) - reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases - }; - this->makeWatPressure(ic, reg, span); - } -} - -template -void PressureTable:: -equil_OWG(const Region& reg, const VSpan& span) -{ - // Datum depth in oil zone. Calculate phase pressure for oil first, - // followed by gas and water if applicable. - - if (! this->oilActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "OIL zone in model without active oil phase" - }; - } - - { - const auto ic = typename OPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeOilPressure(ic, reg, span); - } - - if (this->waterActive()) { - // Pcow = Po - Pw => Pw = Po - Pcow - const auto ic = typename WPress::InitCond { - reg.zwoc(), - this->oil(reg.zwoc()) - reg.pcowWoc() - }; - - this->makeWatPressure(ic, reg, span); - } - - if (this->gasActive()) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - const auto ic = typename GPress::InitCond { - reg.zgoc(), - this->oil(reg.zgoc()) + reg.pcgoGoc() - }; - this->makeGasPressure(ic, reg, span); - } -} - -template -void PressureTable:: -makeOilPressure(const typename OPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = OilPressODE { - this->temperature_, reg.dissolutionCalculator(), - reg.pvtIdx(), this->gravity_ - }; - - this->oil_ = std::make_unique(drho, ic, this->nsample_, span); -} - -template -void PressureTable:: -makeGasPressure(const typename GPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = GasPressODE { - this->temperature_, reg.evaporationCalculator(), reg.waterEvaporationCalculator(), - reg.pvtIdx(), this->gravity_ - }; - - this->gas_ = std::make_unique(drho, ic, this->nsample_, span); -} - -template -void PressureTable:: -makeWatPressure(const typename WPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = WatPressODE { - this->temperature_, reg.saltVdTable(), reg.pvtIdx(), this->gravity_ - }; - - this->wat_ = std::make_unique(drho, ic, this->nsample_, span); -} - -} - -namespace DeckDependent { - -std::vector -getEquil(const EclipseState& state) -{ - const auto& init = state.getInitConfig(); - - if(!init.hasEquil()) { - throw std::domain_error("Deck does not provide equilibration data."); - } - - const auto& equil = init.getEquil(); - return { equil.begin(), equil.end() }; -} - -template -std::vector -equilnum(const EclipseState& eclipseState, - const GridView& gridview) -{ - std::vector eqlnum(gridview.size(0), 0); - - if (eclipseState.fieldProps().has_int("EQLNUM")) { - const auto& e = eclipseState.fieldProps().get_int("EQLNUM"); - std::transform(e.begin(), e.end(), eqlnum.begin(), [](int n){ return n - 1;}); - } - - return eqlnum; -} - -template -template -InitialStateComputer:: -InitialStateComputer(MaterialLawManager& materialLawManager, - const EclipseState& eclipseState, - const Grid& grid, - const GridView& gridView, - const CartesianIndexMapper& cartMapper, - const double grav, - const bool applySwatInit) - : temperature_(grid.size(/*codim=*/0)), - saltConcentration_(grid.size(/*codim=*/0)), - saltSaturation_(grid.size(/*codim=*/0)), - pp_(FluidSystem::numPhases, - std::vector(grid.size(/*codim=*/0))), - sat_(FluidSystem::numPhases, - std::vector(grid.size(/*codim=*/0))), - rs_(grid.size(/*codim=*/0)), - rv_(grid.size(/*codim=*/0)), - rvw_(grid.size(/*codim=*/0)), - cartesianIndexMapper_(cartMapper) -{ - //Check for presence of kw SWATINIT - if (applySwatInit) { - if (eclipseState.fieldProps().has_double("SWATINIT")) { - swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); - } - } - - // Querry cell depth, cell top-bottom. - // numerical aquifer cells might be specified with different depths. - const auto& num_aquifers = eclipseState.aquifer().numericalAquifers(); - updateCellProps_(gridView, num_aquifers); - - // Get the equilibration records. - const std::vector rec = getEquil(eclipseState); - const auto& tables = eclipseState.getTableManager(); - // Create (inverse) region mapping. - const RegionMapping<> eqlmap(equilnum(eclipseState, grid)); - const int invalidRegion = -1; - regionPvtIdx_.resize(rec.size(), invalidRegion); - setRegionPvtIdx(eclipseState, eqlmap); - - // Create Rs functions. - rsFunc_.reserve(rec.size()); - if (FluidSystem::enableDissolvedGas()) { - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) { - rsFunc_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].liveOilInitConstantRs()) { - const TableContainer& rsvdTables = tables.getRsvdTables(); - const TableContainer& pbvdTables = tables.getPbvdTables(); - if (rsvdTables.size() > 0) { - - const RsvdTable& rsvdTable = rsvdTables.getTable(i); - std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); - rsFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn, rsColumn)); - } else if (pbvdTables.size() > 0) { - const PbvdTable& pbvdTable = pbvdTables.getTable(i); - std::vector depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy(); - std::vector pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy(); - rsFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn, pbubColumn)); - - } else { - throw std::runtime_error("Cannot initialise: RSVD or PBVD table not available."); - } - - } - else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - throw std::runtime_error("Cannot initialise: when no explicit RSVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); - } - const double pContact = rec[i].datumDepthPressure(); - const double TContact = 273.15 + 20; // standard temperature for now - rsFunc_.push_back(std::make_shared>(pvtIdx, pContact, TContact)); - } - } - } - else { - for (size_t i = 0; i < rec.size(); ++i) { - rsFunc_.push_back(std::make_shared()); - } - } - - rvFunc_.reserve(rec.size()); - if (FluidSystem::enableVaporizedOil()) { - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) { - rvFunc_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].wetGasInitConstantRv()) { - const TableContainer& rvvdTables = tables.getRvvdTables(); - const TableContainer& pdvdTables = tables.getPdvdTables(); - - if (rvvdTables.size() > 0) { - const RvvdTable& rvvdTable = rvvdTables.getTable(i); - std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); - rvFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn, rvColumn)); - } else if (pdvdTables.size() > 0) { - const PdvdTable& pdvdTable = pdvdTables.getTable(i); - std::vector depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy(); - std::vector pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy(); - rvFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn, pdewColumn)); - } else { - throw std::runtime_error("Cannot initialise: RVVD or PDCD table not available."); - } - } - else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - throw std::runtime_error( - "Cannot initialise: when no explicit RVVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); - } - const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double TContact = 273.15 + 20; // standard temperature for now - rvFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); - } - } - } - else { - for (size_t i = 0; i < rec.size(); ++i) { - rvFunc_.push_back(std::make_shared()); - } - } - - rvwFunc_.reserve(rec.size()); - if (FluidSystem::enableVaporizedWater()) { - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) { - rvwFunc_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].humidGasInitConstantRvw()) { - const TableContainer& rvwvdTables = tables.getRvwvdTables(); - - if (rvwvdTables.size() > 0) { - const RvwvdTable& rvwvdTable = rvwvdTables.getTable(i); - std::vector depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy(); - rvwFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn, rvwvdColumn)); - } else { - throw std::runtime_error("Cannot initialise: RVWVD table not available."); - } - } - else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - throw std::runtime_error( - "Cannot initialise: when no explicit RVWVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); - } - const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double TContact = 273.15 + 20; // standard temperature for now - rvwFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); - } - } - } - else { - for (size_t i = 0; i < rec.size(); ++i) { - rvwFunc_.push_back(std::make_shared()); - } - } - - - // EXTRACT the initial temperature - updateInitialTemperature_(eclipseState); - - // EXTRACT the initial salt concentration - updateInitialSaltConcentration_(eclipseState, eqlmap); - - // EXTRACT the initial salt saturation - updateInitialSaltSaturation_(eclipseState, eqlmap); - - // Compute pressures, saturations, rs and rv factors. - const auto& comm = grid.comm(); - calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav); - - // modify the pressure and saturation for numerical aquifer cells - applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage()); - - // Modify oil pressure in no-oil regions so that the pressures of present phases can - // be recovered from the oil pressure and capillary relations. -} - -template -void InitialStateComputer:: -updateInitialTemperature_(const EclipseState& eclState) -{ - this->temperature_ = eclState.fieldProps().get_double("TEMPI"); -} - -template -template -void InitialStateComputer:: -updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) -{ - const int numEquilReg = rsFunc_.size(); - saltVdTable_.resize(numEquilReg); - const auto& tables = eclState.getTableManager(); - const TableContainer& saltvdTables = tables.getSaltvdTables(); - - // If no saltvd table is given, we create a trivial table for the density calculations - if (saltvdTables.empty()) { - std::vector x = {0.0,1.0}; - std::vector y = {0.0,0.0}; - for (auto& table : this->saltVdTable_) { - table.setXYContainers(x, y); - } - } else { - for (size_t i = 0; i < saltvdTables.size(); ++i) { - const SaltvdTable& saltvdTable = saltvdTables.getTable(i); - saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn()); - - const auto& cells = reg.cells(i); - for (const auto& cell : cells) { - const double depth = cellCenterDepth_[cell]; - this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); - } - } - } -} - -template -template -void InitialStateComputer:: -updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg) -{ - const int numEquilReg = rsFunc_.size(); - saltpVdTable_.resize(numEquilReg); - const auto& tables = eclState.getTableManager(); - const TableContainer& saltpvdTables = tables.getSaltpvdTables(); - - for (size_t i = 0; i < saltpvdTables.size(); ++i) { - const SaltpvdTable& saltpvdTable = saltpvdTables.getTable(i); - saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn()); - - const auto& cells = reg.cells(i); - for (const auto& cell : cells) { - const double depth = cellCenterDepth_[cell]; - this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); - } - } -} - - -template -void InitialStateComputer:: -updateCellProps_(const GridView& gridView, - const NumericalAquifers& aquifer) -{ - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - int numElements = gridView.size(/*codim=*/0); - cellCenterDepth_.resize(numElements); - cellZSpan_.resize(numElements); - cellZMinMax_.resize(numElements); - - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - const auto num_aqu_cells = aquifer.allAquiferCells(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element& element = *elemIt; - const unsigned int elemIdx = elemMapper.index(element); - cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); - const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); - cellZSpan_[elemIdx] = Details::cellZSpan(element); - cellZMinMax_[elemIdx] = Details::cellZMinMax(element); - if (!num_aqu_cells.empty()) { - const auto search = num_aqu_cells.find(cartIx); - if (search != num_aqu_cells.end()) { - const auto* aqu_cell = num_aqu_cells.at(cartIx); - const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx]; - cellCenterDepth_[elemIdx] += depth_change_num_aqu; - cellZSpan_[elemIdx].first += depth_change_num_aqu; - cellZSpan_[elemIdx].second += depth_change_num_aqu; - cellZMinMax_[elemIdx].first += depth_change_num_aqu; - cellZMinMax_[elemIdx].second += depth_change_num_aqu; - } - } - } -} - -template -void InitialStateComputer:: -applyNumericalAquifers_(const GridView& gridView, - const NumericalAquifers& aquifer, - const bool co2store) -{ - const auto num_aqu_cells = aquifer.allAquiferCells(); - if (num_aqu_cells.empty()) return; - - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element& element = *elemIt; - const unsigned int elemIdx = elemMapper.index(element); - const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); - const auto search = num_aqu_cells.find(cartIx); - if (search != num_aqu_cells.end()) { - // numerical aquifer cells are filled with water initially - // for co2store the oilphase may be used for the brine - bool co2store_oil_as_brine = co2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const auto watPos = co2store_oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; - if (FluidSystem::phaseIsActive(watPos)) { - this->sat_[watPos][elemIdx] = 1.; - } else { - throw std::logic_error { "Water phase has to be active for numerical aquifer case" }; - } - - const auto oilPos = FluidSystem::oilPhaseIdx; - if (!co2store && FluidSystem::phaseIsActive(oilPos)) { - this->sat_[oilPos][elemIdx] = 0.; - } - - const auto gasPos = FluidSystem::gasPhaseIdx; - if (FluidSystem::phaseIsActive(gasPos)) { - this->sat_[gasPos][elemIdx] = 0.; - } - const auto* aqu_cell = num_aqu_cells.at(cartIx); - const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL " - "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY", - aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id); - OpmLog::info(msg); - - // if pressure is specified for numerical aquifers, we use these pressure values - // for numerical aquifer cells - if (aqu_cell->init_pressure) { - const double pres = *(aqu_cell->init_pressure); - this->pp_[watPos][elemIdx] = pres; - if (FluidSystem::phaseIsActive(gasPos)) { - this->pp_[gasPos][elemIdx] = pres; - } - if (FluidSystem::phaseIsActive(oilPos)) { - this->pp_[oilPos][elemIdx] = pres; - } - } - } - } -} - -template -template -void InitialStateComputer:: -setRegionPvtIdx(const EclipseState& eclState, const RMap& reg) -{ - const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM"); - - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1; - } -} - -template -template -void InitialStateComputer:: -calcPressSatRsRv(const RMap& reg, - const std::vector& rec, - MaterialLawManager& materialLawManager, - const Comm& comm, - const double grav) -{ - using PhaseSat = Details::PhaseSaturations< - MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId - >; - - auto ptable = Details::PressureTable{ grav }; - auto psat = PhaseSat { materialLawManager, this->swatInit_ }; - auto vspan = std::array{}; - - std::vector regionIsEmpty(rec.size(), 0); - for (size_t r = 0; r < rec.size(); ++r) { - const auto& cells = reg.cells(r); - - Details::verticalExtent(cells, cellZMinMax_, comm, vspan); - - const auto acc = rec[r].initializationTargetAccuracy(); - if (acc > 0) { - throw std::runtime_error { - "Cannot initialise model: Positive item 9 is not supported " - "in EQUIL keyword, record " + std::to_string(r + 1) - }; - } - - if (cells.empty()) { - regionIsEmpty[r] = 1; - continue; - } - - const auto eqreg = EquilReg { - rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r] - }; - - // Ensure gas/oil and oil/water contacts are within the span for the - // phase pressure calculation. - vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc())); - vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc())); - - ptable.equilibrate(eqreg, vspan); - - if (acc == 0) { - // Centre-point method - this->equilibrateCellCentres(cells, eqreg, ptable, psat); - } - else if (acc < 0) { - // Horizontal subdivision - this->equilibrateHorizontal(cells, eqreg, -acc, - ptable, psat); - } else { - // Horizontal subdivision with titled fault blocks - // the simulator throw a few line above for the acc > 0 case - // i.e. we should not reach here. - assert(false); - } - } - comm.min(regionIsEmpty.data(),regionIsEmpty.size()); - if (comm.rank() == 0) { - for (size_t r = 0; r < rec.size(); ++r) { - if (regionIsEmpty[r]) //region is empty on all partitions - OpmLog::warning("Equilibration region " + std::to_string(r + 1) - + " has no active cells"); - } - } -} - -template -template -void InitialStateComputer:: -cellLoop(const CellRange& cells, - EquilibrationMethod&& eqmethod) -{ - const auto oilPos = FluidSystem::oilPhaseIdx; - const auto gasPos = FluidSystem::gasPhaseIdx; - const auto watPos = FluidSystem::waterPhaseIdx; - - const auto oilActive = FluidSystem::phaseIsActive(oilPos); - const auto gasActive = FluidSystem::phaseIsActive(gasPos); - const auto watActive = FluidSystem::phaseIsActive(watPos); - - auto pressures = Details::PhaseQuantityValue{}; - auto saturations = Details::PhaseQuantityValue{}; - auto Rs = 0.0; - auto Rv = 0.0; - auto Rvw = 0.0; - - for (const auto& cell : cells) { - eqmethod(cell, pressures, saturations, Rs, Rv, Rvw); - - if (oilActive) { - this->pp_ [oilPos][cell] = pressures.oil; - this->sat_[oilPos][cell] = saturations.oil; - } - - if (gasActive) { - this->pp_ [gasPos][cell] = pressures.gas; - this->sat_[gasPos][cell] = saturations.gas; - } - - if (watActive) { - this->pp_ [watPos][cell] = pressures.water; - this->sat_[watPos][cell] = saturations.water; - } - - if (oilActive && gasActive) { - this->rs_[cell] = Rs; - this->rv_[cell] = Rv; - } - - if (watActive && gasActive) { - this->rvw_[cell] = Rvw; - } - } -} - -template -template -void InitialStateComputer:: -equilibrateCellCentres(const CellRange& cells, - const EquilReg& eqreg, - const PressTable& ptable, - PhaseSat& psat) -{ - using CellPos = typename PhaseSat::Position; - using CellID = std::remove_cv_t().cell)>>; - this->cellLoop(cells, [this, &eqreg, &ptable, &psat] - (const CellID cell, - Details::PhaseQuantityValue& pressures, - Details::PhaseQuantityValue& saturations, - double& Rs, - double& Rv, - double& Rvw) -> void - { - const auto pos = CellPos { - cell, cellCenterDepth_[cell] - }; - - saturations = psat.deriveSaturations(pos, eqreg, ptable); - pressures = psat.correctedPhasePressures(); - - const auto temp = this->temperature_[cell]; - - Rs = eqreg.dissolutionCalculator() - (pos.depth, pressures.oil, temp, saturations.gas); - - Rv = eqreg.evaporationCalculator() - (pos.depth, pressures.gas, temp, saturations.oil); - - Rvw = eqreg.waterEvaporationCalculator() - (pos.depth, pressures.gas, temp, saturations.water); - }); -} - -template -template -void InitialStateComputer:: -equilibrateHorizontal(const CellRange& cells, - const EquilReg& eqreg, - const int acc, - const PressTable& ptable, - PhaseSat& psat) -{ - using CellPos = typename PhaseSat::Position; - using CellID = std::remove_cv_t().cell)>>; - - this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat] - (const CellID cell, - Details::PhaseQuantityValue& pressures, - Details::PhaseQuantityValue& saturations, - double& Rs, - double& Rv, - double& Rvw) -> void - { - pressures .reset(); - saturations.reset(); - - auto totfrac = 0.0; - for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) { - const auto pos = CellPos { cell, depth }; - - saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac); - pressures .axpy(psat.correctedPhasePressures(), frac); - - totfrac += frac; - } - - if (totfrac > 0.) { - saturations /= totfrac; - pressures /= totfrac; - } else { - // Fall back to centre point method for zero-thickness cells. - const auto pos = CellPos { - cell, cellCenterDepth_[cell] - }; - - saturations = psat.deriveSaturations(pos, eqreg, ptable); - pressures = psat.correctedPhasePressures(); - } - - const auto temp = this->temperature_[cell]; - const auto cz = cellCenterDepth_[cell]; - - Rs = eqreg.dissolutionCalculator() - (cz, pressures.oil, temp, saturations.gas); - - Rv = eqreg.evaporationCalculator() - (cz, pressures.gas, temp, saturations.oil); - - Rvw = eqreg.waterEvaporationCalculator() - (cz, pressures.gas, temp, saturations.water); - }); -} - +namespace DeckDependent { #if HAVE_DUNE_FEM using GridView = Dune::Fem::GridPart2GridViewImpl< Dune::Fem::AdaptiveLeafGridPart< @@ -2020,6 +105,6 @@ namespace Details { template std::pair cellZMinMax(const Dune::cpgrid::Entity<0>& element); } - +} } // namespace EQUIL -} // namespace Opm + // namespace Opm diff --git a/ebos/equil/initstateequil_impl.hh b/ebos/equil/initstateequil_impl.hh new file mode 100644 index 000000000..9462069e9 --- /dev/null +++ b/ebos/equil/initstateequil_impl.hh @@ -0,0 +1,1950 @@ +// -*- 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. +*/ +#ifndef EWOMS_INITSTATEEQUIL_IMPL_HH +#define EWOMS_INITSTATEEQUIL_IMPL_HH + + +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +#include +#include +#include +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif + +#if HAVE_DUNE_ALUGRID +#include +#include +#endif // HAVE_DUNE_ALUGRID + +namespace Opm { +namespace EQUIL { + +namespace Details { + +template +void verticalExtent(const CellRange& cells, + const std::vector>& cellZMinMax, + const Comm& comm, + std::array& span) +{ + span[0] = std::numeric_limits::max(); + span[1] = std::numeric_limits::lowest(); + + // Define vertical span as + // + // [minimum(node depth(cells)), maximum(node depth(cells))] + // + // Note: The implementation of 'RK4IVP<>' implicitly + // imposes the requirement that cell centroids are all + // within this vertical span. That requirement is not + // checked. + for (const auto& cell : cells) { + if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; } + if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; } + } + span[0] = comm.min(span[0]); + span[1] = comm.max(span[1]); +} + +void subdivisionCentrePoints(const double left, + const double right, + const int numIntervals, + std::vector>& subdiv) +{ + const auto h = (right - left) / numIntervals; + + auto end = left; + for (auto i = 0*numIntervals; i < numIntervals; ++i) { + const auto start = end; + end = left + (i + 1)*h; + + subdiv.emplace_back((start + end) / 2, h); + } +} + +template +std::vector> +horizontalSubdivision(const CellID cell, + const std::pair topbot, + const int numIntervals) +{ + auto subdiv = std::vector>{}; + subdiv.reserve(2 * numIntervals); + + if (topbot.first > topbot.second) { + throw std::out_of_range { + "Negative thickness (inverted top/bottom faces) in cell " + + std::to_string(cell) + }; + } + + subdivisionCentrePoints(topbot.first, topbot.second, + 2*numIntervals, subdiv); + + return subdiv; +} + +template +double cellCenterDepth(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; + double zz = 0.0; + + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + for (int i=0; i < corners; ++i) + zz += geometry.corner(i)[zCoord]; + + return zz/corners; +} + +template +std::pair cellZSpan(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; + double bot = 0.0; + double top = 0.0; + + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + assert(corners == 8); + for (int i=0; i < 4; ++i) + bot += geometry.corner(i)[zCoord]; + for (int i=4; i < corners; ++i) + top += geometry.corner(i)[zCoord]; + + return std::make_pair(bot/4, top/4); +} + +template +std::pair cellZMinMax(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + assert(corners == 8); + auto min = std::numeric_limits::max(); + auto max = std::numeric_limits::lowest(); + + + for (int i=0; i < corners; ++i) { + min = std::min(min, geometry.corner(i)[zCoord]); + max = std::max(max, geometry.corner(i)[zCoord]); + } + return std::make_pair(min, max); +} + +template +RK4IVP::RK4IVP(const RHS& f, + const std::array& span, + const double y0, + const int N) + : N_(N) + , span_(span) +{ + const double h = stepsize(); + const double h2 = h / 2; + const double 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 double x = span_[0] + i*h; + const double y = y_.back(); + + const double k1 = f_[i]; + const double k2 = f(x + h2, y + h2*k1); + const double k3 = f(x + h2, y + h2*k2); + const double 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() == std::vector::size_type(N + 1)); +} + +template +double RK4IVP:: +operator()(const double x) const +{ + // Dense output (O(h**3)) according to Shampine + // (Hermite interpolation) + const double h = stepsize(); + int i = (x - span_[0]) / h; + const double t = (x - (span_[0] + i*h)) / h; + + // Crude handling of evaluation point outside "span_"; + if (i < 0) { i = 0; } + if (N_ <= i) { i = N_ - 1; } + + const double y0 = y_[i], y1 = y_[i + 1]; + const double f0 = f_[i], f1 = f_[i + 1]; + + double 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 +double RK4IVP:: +stepsize() const +{ + return (span_[1] - span_[0]) / N_; +} + +namespace PhasePressODE { + +template +Water:: +Water(const double temp, + const TabulatedFunction& saltVdTable, + const int pvtRegionIdx, + const double normGrav) + : temp_(temp) + , saltVdTable_(saltVdTable) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +double Water:: +operator()(const double depth, + const double press) const +{ + return this->density(depth, press) * g_; +} + +template +double Water:: +density(const double depth, + const double press) const +{ + // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate + double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); + double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0 /*=Rsw*/, saltConcentration); + rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); + return rho; +} + +template +Oil:: +Oil(const double temp, + const RS& rs, + const int pvtRegionIdx, + const double normGrav) + : temp_(temp) + , rs_(rs) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +double Oil:: +operator()(const double depth, + const double press) const +{ + return this->density(depth, press) * g_; +} + +template +double Oil:: +density(const double depth, + const double press) const +{ + double rs = 0.0; + if(FluidSystem::enableDissolvedGas()) + rs = rs_(depth, press, temp_); + + double bOil = 0.0; + if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press)) { + bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } + else { + bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); + } + double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + if (FluidSystem::enableDissolvedGas()) { + rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + } + + return rho; +} + +template +Gas:: +Gas(const double temp, + const RV& rv, + const RVW& rvw, + const int pvtRegionIdx, + const double normGrav) + : temp_(temp) + , rv_(rv) + , rvw_(rvw) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +double Gas:: +operator()(const double depth, + const double press) const +{ + return this->density(depth, press) * g_; +} + +template +double Gas:: +density(const double depth, + const double press) const +{ + double rv = 0.0; + if (FluidSystem::enableVaporizedOil()) + rv = rv_(depth, press, temp_); + + double rvw = 0.0; + if (FluidSystem::enableVaporizedWater()) + rvw = rvw_(depth, press, temp_); + + double bGas = 0.0; + + if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { + if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) + && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) + { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, rvw); + } + double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_) + + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); + return rho; + } + + if (FluidSystem::enableVaporizedOil()){ + if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=rvw*/); + } + double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + return rho; + } + + if (FluidSystem::enableVaporizedWater()){ + if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } + else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, rvw); + } + double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); + return rho; + } + + // immiscible gas + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, 0.0/*=rvw*/); + double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + + return rho; +} + +} + +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 +double +PressureTable:: +PressureFunction:: +value(const double 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:: +checkPtr(const PressFunc* phasePress, + const std::string& phaseName) const +{ + if (phasePress != nullptr) { return; } + + throw std::invalid_argument { + "Phase pressure function for \"" + phaseName + + "\" most not be null" + }; +} + +template +typename PressureTable::Strategy +PressureTable:: +selectEquilibrationStrategy(const Region& reg) const +{ + if (!this->oilActive()) { + if (reg.datum() > reg.zwoc()) { // Datum in water zone + return &PressureTable::equil_WOG; + } + return &PressureTable::equil_GOW; + } + + if (reg.datum() > reg.zwoc()) { // Datum in water zone + return &PressureTable::equil_WOG; + } + else if (reg.datum() < reg.zgoc()) { // Datum in gas zone + return &PressureTable::equil_GOW; + } + else { // Datum in oil zone + return &PressureTable::equil_OWG; + } +} + +template +void PressureTable:: +copyInPointers(const PressureTable& rhs) +{ + if (rhs.oil_ != nullptr) { + this->oil_ = std::make_unique(*rhs.oil_); + } + + if (rhs.gas_ != nullptr) { + this->gas_ = std::make_unique(*rhs.gas_); + } + + if (rhs.wat_ != nullptr) { + this->wat_ = std::make_unique(*rhs.wat_); + } +} + +template +PhaseSaturations:: +PhaseSaturations(MaterialLawManager& matLawMgr, + const std::vector& swatInit) + : matLawMgr_(matLawMgr) + , swatInit_ (swatInit) +{ +} + +template +PhaseSaturations:: +PhaseSaturations(const PhaseSaturations& rhs) + : matLawMgr_(rhs.matLawMgr_) + , swatInit_ (rhs.swatInit_) + , sat_ (rhs.sat_) + , press_ (rhs.press_) +{ + // Note: We don't need to do anything to the 'fluidState_' here. + this->setEvaluationPoint(*rhs.evalPt_.position, + *rhs.evalPt_.region, + *rhs.evalPt_.ptable); +} + +template +const PhaseQuantityValue& +PhaseSaturations:: +deriveSaturations(const Position& x, + const Region& reg, + const PTable& ptable) +{ + this->setEvaluationPoint(x, reg, ptable); + this->initializePhaseQuantities(); + + if (ptable.waterActive()) { this->deriveWaterSat(); } + if (ptable.gasActive()) { this->deriveGasSat(); } + + if (this->isOverlappingTransition()) { + this->fixUnphysicalTransition(); + } + + if (ptable.oilActive()) { this->deriveOilSat(); } + + this->accountForScaledSaturations(); + + return this->sat_; +} + +template +void PhaseSaturations:: +setEvaluationPoint(const Position& x, + const Region& reg, + const PTable& ptable) +{ + this->evalPt_.position = &x; + this->evalPt_.region = ® + this->evalPt_.ptable = &ptable; +} + +template +void PhaseSaturations:: +initializePhaseQuantities() +{ + this->sat_.reset(); + this->press_.reset(); + + const auto depth = this->evalPt_.position->depth; + const auto& ptable = *this->evalPt_.ptable; + + if (ptable.oilActive()) { + this->press_.oil = ptable.oil(depth); + } + + if (ptable.gasActive()) { + this->press_.gas = ptable.gas(depth); + } + + if (ptable.waterActive()) { + this->press_.water = ptable.water(depth); + } +} + +template +void PhaseSaturations::deriveOilSat() +{ + this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas; +} + +template +void PhaseSaturations::deriveGasSat() +{ + auto& sg = this->sat_.gas; + + const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg. + const auto oilActive = this->evalPt_.ptable->oilActive(); + + if (this->isConstCapPress(this->gasPos())) { + // Sharp interface between phases. Can derive phase saturation + // directly from knowing where 'depth' of evaluation point is + // relative to depth of O/G contact. + const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc(); + sg = this->fromDepthTable(gas_contact, + this->gasPos(), isIncr); + } + else { + // Capillary pressure curve is non-constant, meaning there is a + // transition zone between the gas and oil phases. Invert capillary + // pressure relation + // + // Pcgo(Sg) = Pg - Po + // + // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg). + const auto pw = oilActive? this->press_.oil : this->press_.water; + const auto pcgo = this->press_.gas - pw; + sg = this->invertCapPress(pcgo, this->gasPos(), isIncr); + } +} + +template +void PhaseSaturations::deriveWaterSat() +{ + auto& sw = this->sat_.water; + + const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. + + if (this->isConstCapPress(this->waterPos())) { + // Sharp interface between phases. Can derive phase saturation + // directly from knowing where 'depth' of evaluation point is + // relative to depth of O/W contact. + sw = this->fromDepthTable(this->evalPt_.region->zwoc(), + this->waterPos(), isIncr); + } + else { + // Capillary pressure curve is non-constant, meaning there is a + // transition zone between the oil and water phases. Invert + // capillary pressure relation + // + // Pcow(Sw) = Po - Pw + // + // unless the model uses "SWATINIT". In the latter case, pick the + // saturation directly from the SWATINIT array of the pertinent + // cell. + const auto pcow = this->press_.oil - this->press_.water; + + sw = this->swatInit_.empty() + ? this->invertCapPress(pcow, this->waterPos(), isIncr) + : this->applySwatInit(pcow); + } +} + +template +void PhaseSaturations:: +fixUnphysicalTransition() +{ + auto& sg = this->sat_.gas; + auto& sw = this->sat_.water; + + // Overlapping gas/oil and oil/water transition zones can lead to + // unphysical phase saturations when individual saturations are derived + // directly from inverting O/G and O/W capillary pressure curves. + // + // Recalculate phase saturations using the implied gas/water capillary + // pressure: Pg - Pw. + const auto pcgw = this->press_.gas - this->press_.water; + if (! this->swatInit_.empty()) { + // Re-scale Pc to reflect imposed sw for vanishing oil phase. This + // seems consistent with ECLIPSE, but fails to honour SWATINIT in + // case of non-trivial gas/oil capillary pressure. + sw = this->applySwatInit(pcgw, sw); + } + + sw = satFromSumOfPcs + (this->matLawMgr_, this->waterPos(), this->gasPos(), + this->evalPt_.position->cell, pcgw); + sg = 1.0 - sw; + + this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg); + this->fluidState_.setSaturation(this->gasPos(), sg); + this->fluidState_.setSaturation(this->waterPos(), this->evalPt_ + .ptable->waterActive() ? sw : 0.0); + + // Pcgo = Pg - Po => Po = Pg - Pcgo + this->computeMaterialLawCapPress(); + this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); +} + +template +void PhaseSaturations:: +accountForScaledSaturations() +{ + const auto gasActive = this->evalPt_.ptable->gasActive(); + const auto watActive = this->evalPt_.ptable->waterActive(); + const auto oilActive = this->evalPt_.ptable->oilActive(); + + auto sg = gasActive? this->sat_.gas : 0.0; + auto sw = watActive? this->sat_.water : 0.0; + auto so = oilActive? this->sat_.oil : 0.0; + + this->fluidState_.setSaturation(this->waterPos(), sw); + this->fluidState_.setSaturation(this->oilPos(), so); + this->fluidState_.setSaturation(this->gasPos(), sg); + + const auto& scaledDrainageInfo = this->matLawMgr_ + .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell); + + const auto thresholdSat = 1.0e-6; + if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) { + // Water saturation exceeds maximum possible value. Reset oil phase + // pressure to that which corresponds to maximum possible water + // saturation value. + this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu); + } + sw = scaledDrainageInfo.Swu; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcow = Po - Pw => Po = Pw + Pcow + this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } + + } + if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { + // Gas saturation exceeds maximum possible value. Reset oil phase + // pressure to that which corresponds to maximum possible gas + // saturation value. + this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu); + } + sg = scaledDrainageInfo.Sgu; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcgo = Pg - Po => Po = Pg - Pcgo + this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } + } + + if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) { + // Water saturation less than minimum possible value in cell. Reset + // water phase pressure to that which corresponds to minimum + // possible water saturation value. + this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl); + } + sw = scaledDrainageInfo.Swl; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcwo = Po - Pw => Pw = Po - Pcow + this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } + } + + if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { + // Gas saturation less than minimum possible value in cell. Reset + // gas phase pressure to that which corresponds to minimum possible + // gas saturation. + this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl); + } + sg = scaledDrainageInfo.Sgl; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } + } +} + +template +double PhaseSaturations:: +applySwatInit(const double pcow) +{ + return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); +} + +template +double PhaseSaturations:: +applySwatInit(const double pcow, const double sw) +{ + return this->matLawMgr_ + .applySwatinit(this->evalPt_.position->cell, pcow, sw); +} + +template +void PhaseSaturations:: +computeMaterialLawCapPress() +{ + const auto& matParams = this->matLawMgr_ + .materialLawParams(this->evalPt_.position->cell); + + this->matLawCapPress_.fill(0.0); + MaterialLaw::capillaryPressures(this->matLawCapPress_, + matParams, this->fluidState_); +} + +template +double PhaseSaturations:: +materialLawCapPressGasOil() const +{ + return this->matLawCapPress_[this->oilPos()] + + this->matLawCapPress_[this->gasPos()]; +} + +template +double PhaseSaturations:: +materialLawCapPressOilWater() const +{ + return this->matLawCapPress_[this->oilPos()] + - this->matLawCapPress_[this->waterPos()]; +} + +template +double PhaseSaturations:: +materialLawCapPressGasWater() const +{ + return this->matLawCapPress_[this->gasPos()] + - this->matLawCapPress_[this->waterPos()]; +} + +template +bool PhaseSaturations:: +isConstCapPress(const PhaseIdx phaseIdx) const +{ + return isConstPc + (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell); +} + +template +bool PhaseSaturations:: +isOverlappingTransition() const +{ + return this->evalPt_.ptable->gasActive() + && this->evalPt_.ptable->waterActive() + && ((this->sat_.gas + this->sat_.water) > 1.0); +} + +template +double PhaseSaturations:: +fromDepthTable(const double contactdepth, + const PhaseIdx phasePos, + const bool isincr) const +{ + return satFromDepth + (this->matLawMgr_, this->evalPt_.position->depth, + contactdepth, static_cast(phasePos), + this->evalPt_.position->cell, isincr); +} + +template +double PhaseSaturations:: +invertCapPress(const double pc, + const PhaseIdx phasePos, + const bool isincr) const +{ + return satFromPc + (this->matLawMgr_, static_cast(phasePos), + this->evalPt_.position->cell, pc, isincr); +} + +template +PressureTable:: +PressureTable(const double gravity, + const int samplePoints) + : gravity_(gravity) + , nsample_(samplePoints) +{ +} + +template +PressureTable:: +PressureTable(const PressureTable& rhs) + : gravity_(rhs.gravity_) + , nsample_(rhs.nsample_) +{ + this->copyInPointers(rhs); +} + +template +PressureTable:: +PressureTable(PressureTable&& rhs) + : gravity_(rhs.gravity_) + , nsample_(rhs.nsample_) + , oil_ (std::move(rhs.oil_)) + , gas_ (std::move(rhs.gas_)) + , wat_ (std::move(rhs.wat_)) +{ +} + +template +PressureTable& +PressureTable:: +operator=(const PressureTable& rhs) +{ + this->gravity_ = rhs.gravity_; + this->nsample_ = rhs.nsample_; + this->copyInPointers(rhs); + + return *this; +} + +template +PressureTable& +PressureTable:: +operator=(PressureTable&& rhs) +{ + this->gravity_ = rhs.gravity_; + this->nsample_ = rhs.nsample_; + + this->oil_ = std::move(rhs.oil_); + this->gas_ = std::move(rhs.gas_); + this->wat_ = std::move(rhs.wat_); + + return *this; +} + +template +void PressureTable:: +equilibrate(const Region& reg, + const VSpan& span) +{ + // One of the PressureTable::equil_*() member functions. + auto equil = this->selectEquilibrationStrategy(reg); + + (this->*equil)(reg, span); +} + +template +bool PressureTable:: +oilActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); +} + +template +bool PressureTable:: +gasActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); +} + +template +bool PressureTable:: +waterActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); +} + +template +double PressureTable:: +oil(const double depth) const +{ + this->checkPtr(this->oil_.get(), "OIL"); + + return this->oil_->value(depth); +} + +template +double PressureTable:: +gas(const double depth) const +{ + this->checkPtr(this->gas_.get(), "GAS"); + + return this->gas_->value(depth); +} + + +template +double PressureTable:: +water(const double depth) const +{ + this->checkPtr(this->wat_.get(), "WATER"); + + return this->wat_->value(depth); +} + +template +void PressureTable:: +equil_WOG(const Region& reg, const VSpan& span) +{ + // Datum depth in water zone. Calculate phase pressure for water first, + // followed by oil and gas if applicable. + + if (! this->waterActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "WATER zone in model without active water phase" + }; + } + + { + const auto ic = typename WPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeWatPressure(ic, reg, span); + } + + if (this->oilActive()) { + // Pcow = Po - Pw => Po = Pw + Pcow + const auto ic = typename OPress::InitCond { + reg.zwoc(), + this->water(reg.zwoc()) + reg.pcowWoc() + }; + + this->makeOilPressure(ic, reg, span); + } + + if (this->gasActive() && this->oilActive()) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + const auto ic = typename GPress::InitCond { + reg.zgoc(), + this->oil(reg.zgoc()) + reg.pcgoGoc() + }; + + this->makeGasPressure(ic, reg, span); + } else if (this->gasActive() && !this->oilActive()) { + // No oil phase set Pg = Pw + Pcgw + const auto ic = typename GPress::InitCond { + reg.zwoc(), // The WOC is really the GWC for gas/water cases + this->water(reg.zwoc()) + reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases + }; + this->makeGasPressure(ic, reg, span); + } +} + +template +void PressureTable:: +equil_GOW(const Region& reg, const VSpan& span) +{ + // Datum depth in gas zone. Calculate phase pressure for gas first, + // followed by oil and water if applicable. + + if (! this->gasActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "GAS zone in model without active gas phase" + }; + } + + { + const auto ic = typename GPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeGasPressure(ic, reg, span); + } + + if (this->oilActive()) { + // Pcgo = Pg - Po => Po = Pg - Pcgo + const auto ic = typename OPress::InitCond { + reg.zgoc(), + this->gas(reg.zgoc()) - reg.pcgoGoc() + }; + this->makeOilPressure(ic, reg, span); + } + + if (this->waterActive() && this->oilActive()) { + // Pcow = Po - Pw => Pw = Po - Pcow + const auto ic = typename WPress::InitCond { + reg.zwoc(), + this->oil(reg.zwoc()) - reg.pcowWoc() + }; + + this->makeWatPressure(ic, reg, span); + } else if (this->waterActive() && !this->oilActive()) { + // No oil phase set Pw = Pg - Pcgw + const auto ic = typename WPress::InitCond { + reg.zwoc(), // The WOC is really the GWC for gas/water cases + this->gas(reg.zwoc()) - reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases + }; + this->makeWatPressure(ic, reg, span); + } +} + +template +void PressureTable:: +equil_OWG(const Region& reg, const VSpan& span) +{ + // Datum depth in oil zone. Calculate phase pressure for oil first, + // followed by gas and water if applicable. + + if (! this->oilActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "OIL zone in model without active oil phase" + }; + } + + { + const auto ic = typename OPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeOilPressure(ic, reg, span); + } + + if (this->waterActive()) { + // Pcow = Po - Pw => Pw = Po - Pcow + const auto ic = typename WPress::InitCond { + reg.zwoc(), + this->oil(reg.zwoc()) - reg.pcowWoc() + }; + + this->makeWatPressure(ic, reg, span); + } + + if (this->gasActive()) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + const auto ic = typename GPress::InitCond { + reg.zgoc(), + this->oil(reg.zgoc()) + reg.pcgoGoc() + }; + this->makeGasPressure(ic, reg, span); + } +} + +template +void PressureTable:: +makeOilPressure(const typename OPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = OilPressODE { + this->temperature_, reg.dissolutionCalculator(), + reg.pvtIdx(), this->gravity_ + }; + + this->oil_ = std::make_unique(drho, ic, this->nsample_, span); +} + +template +void PressureTable:: +makeGasPressure(const typename GPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = GasPressODE { + this->temperature_, reg.evaporationCalculator(), reg.waterEvaporationCalculator(), + reg.pvtIdx(), this->gravity_ + }; + + this->gas_ = std::make_unique(drho, ic, this->nsample_, span); +} + +template +void PressureTable:: +makeWatPressure(const typename WPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = WatPressODE { + this->temperature_, reg.saltVdTable(), reg.pvtIdx(), this->gravity_ + }; + + this->wat_ = std::make_unique(drho, ic, this->nsample_, span); +} + +} + +namespace DeckDependent { + +std::vector +getEquil(const EclipseState& state) +{ + const auto& init = state.getInitConfig(); + + if(!init.hasEquil()) { + throw std::domain_error("Deck does not provide equilibration data."); + } + + const auto& equil = init.getEquil(); + return { equil.begin(), equil.end() }; +} + +template +std::vector +equilnum(const EclipseState& eclipseState, + const GridView& gridview) +{ + std::vector eqlnum(gridview.size(0), 0); + + if (eclipseState.fieldProps().has_int("EQLNUM")) { + const auto& e = eclipseState.fieldProps().get_int("EQLNUM"); + std::transform(e.begin(), e.end(), eqlnum.begin(), [](int n){ return n - 1;}); + } + + return eqlnum; +} + +template +template +InitialStateComputer:: +InitialStateComputer(MaterialLawManager& materialLawManager, + const EclipseState& eclipseState, + const Grid& grid, + const GridView& gridView, + const CartesianIndexMapper& cartMapper, + const double grav, + const bool applySwatInit) + : temperature_(grid.size(/*codim=*/0)), + saltConcentration_(grid.size(/*codim=*/0)), + saltSaturation_(grid.size(/*codim=*/0)), + pp_(FluidSystem::numPhases, + std::vector(grid.size(/*codim=*/0))), + sat_(FluidSystem::numPhases, + std::vector(grid.size(/*codim=*/0))), + rs_(grid.size(/*codim=*/0)), + rv_(grid.size(/*codim=*/0)), + rvw_(grid.size(/*codim=*/0)), + cartesianIndexMapper_(cartMapper) +{ + //Check for presence of kw SWATINIT + if (applySwatInit) { + if (eclipseState.fieldProps().has_double("SWATINIT")) { + swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); + } + } + + // Querry cell depth, cell top-bottom. + // numerical aquifer cells might be specified with different depths. + const auto& num_aquifers = eclipseState.aquifer().numericalAquifers(); + updateCellProps_(gridView, num_aquifers); + + // Get the equilibration records. + const std::vector rec = getEquil(eclipseState); + const auto& tables = eclipseState.getTableManager(); + // Create (inverse) region mapping. + const RegionMapping<> eqlmap(equilnum(eclipseState, grid)); + const int invalidRegion = -1; + regionPvtIdx_.resize(rec.size(), invalidRegion); + setRegionPvtIdx(eclipseState, eqlmap); + + // Create Rs functions. + rsFunc_.reserve(rec.size()); + if (FluidSystem::enableDissolvedGas()) { + for (size_t i = 0; i < rec.size(); ++i) { + if (eqlmap.cells(i).empty()) { + rsFunc_.push_back(std::shared_ptr>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].liveOilInitConstantRs()) { + const TableContainer& rsvdTables = tables.getRsvdTables(); + const TableContainer& pbvdTables = tables.getPbvdTables(); + if (rsvdTables.size() > 0) { + + const RsvdTable& rsvdTable = rsvdTables.getTable(i); + std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); + rsFunc_.push_back(std::make_shared>(pvtIdx, + depthColumn, rsColumn)); + } else if (pbvdTables.size() > 0) { + const PbvdTable& pbvdTable = pbvdTables.getTable(i); + std::vector depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy(); + std::vector pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy(); + rsFunc_.push_back(std::make_shared>(pvtIdx, + depthColumn, pbubColumn)); + + } else { + throw std::runtime_error("Cannot initialise: RSVD or PBVD table not available."); + } + + } + else { + if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { + throw std::runtime_error("Cannot initialise: when no explicit RSVD table is given, \n" + "datum depth must be at the gas-oil-contact. " + "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); + } + const double pContact = rec[i].datumDepthPressure(); + const double TContact = 273.15 + 20; // standard temperature for now + rsFunc_.push_back(std::make_shared>(pvtIdx, pContact, TContact)); + } + } + } + else { + for (size_t i = 0; i < rec.size(); ++i) { + rsFunc_.push_back(std::make_shared()); + } + } + + rvFunc_.reserve(rec.size()); + if (FluidSystem::enableVaporizedOil()) { + for (size_t i = 0; i < rec.size(); ++i) { + if (eqlmap.cells(i).empty()) { + rvFunc_.push_back(std::shared_ptr>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].wetGasInitConstantRv()) { + const TableContainer& rvvdTables = tables.getRvvdTables(); + const TableContainer& pdvdTables = tables.getPdvdTables(); + + if (rvvdTables.size() > 0) { + const RvvdTable& rvvdTable = rvvdTables.getTable(i); + std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); + rvFunc_.push_back(std::make_shared>(pvtIdx, + depthColumn, rvColumn)); + } else if (pdvdTables.size() > 0) { + const PdvdTable& pdvdTable = pdvdTables.getTable(i); + std::vector depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy(); + std::vector pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy(); + rvFunc_.push_back(std::make_shared>(pvtIdx, + depthColumn, pdewColumn)); + } else { + throw std::runtime_error("Cannot initialise: RVVD or PDCD table not available."); + } + } + else { + if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { + throw std::runtime_error( + "Cannot initialise: when no explicit RVVD table is given, \n" + "datum depth must be at the gas-oil-contact. " + "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); + } + const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); + const double TContact = 273.15 + 20; // standard temperature for now + rvFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); + } + } + } + else { + for (size_t i = 0; i < rec.size(); ++i) { + rvFunc_.push_back(std::make_shared()); + } + } + + rvwFunc_.reserve(rec.size()); + if (FluidSystem::enableVaporizedWater()) { + for (size_t i = 0; i < rec.size(); ++i) { + if (eqlmap.cells(i).empty()) { + rvwFunc_.push_back(std::shared_ptr>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].humidGasInitConstantRvw()) { + const TableContainer& rvwvdTables = tables.getRvwvdTables(); + + if (rvwvdTables.size() > 0) { + const RvwvdTable& rvwvdTable = rvwvdTables.getTable(i); + std::vector depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy(); + rvwFunc_.push_back(std::make_shared>(pvtIdx, + depthColumn, rvwvdColumn)); + } else { + throw std::runtime_error("Cannot initialise: RVWVD table not available."); + } + } + else { + if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { + throw std::runtime_error( + "Cannot initialise: when no explicit RVWVD table is given, \n" + "datum depth must be at the gas-oil-contact. " + "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); + } + const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); + const double TContact = 273.15 + 20; // standard temperature for now + rvwFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); + } + } + } + else { + for (size_t i = 0; i < rec.size(); ++i) { + rvwFunc_.push_back(std::make_shared()); + } + } + + + // EXTRACT the initial temperature + updateInitialTemperature_(eclipseState); + + // EXTRACT the initial salt concentration + updateInitialSaltConcentration_(eclipseState, eqlmap); + + // EXTRACT the initial salt saturation + updateInitialSaltSaturation_(eclipseState, eqlmap); + + // Compute pressures, saturations, rs and rv factors. + const auto& comm = grid.comm(); + calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav); + + // modify the pressure and saturation for numerical aquifer cells + applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage()); + + // Modify oil pressure in no-oil regions so that the pressures of present phases can + // be recovered from the oil pressure and capillary relations. +} + +template +void InitialStateComputer:: +updateInitialTemperature_(const EclipseState& eclState) +{ + this->temperature_ = eclState.fieldProps().get_double("TEMPI"); +} + +template +template +void InitialStateComputer:: +updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) +{ + const int numEquilReg = rsFunc_.size(); + saltVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + const TableContainer& saltvdTables = tables.getSaltvdTables(); + + // If no saltvd table is given, we create a trivial table for the density calculations + if (saltvdTables.empty()) { + std::vector x = {0.0,1.0}; + std::vector y = {0.0,0.0}; + for (auto& table : this->saltVdTable_) { + table.setXYContainers(x, y); + } + } else { + for (size_t i = 0; i < saltvdTables.size(); ++i) { + const SaltvdTable& saltvdTable = saltvdTables.getTable(i); + saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn()); + + const auto& cells = reg.cells(i); + for (const auto& cell : cells) { + const double depth = cellCenterDepth_[cell]; + this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); + } + } + } +} + +template +template +void InitialStateComputer:: +updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg) +{ + const int numEquilReg = rsFunc_.size(); + saltpVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + const TableContainer& saltpvdTables = tables.getSaltpvdTables(); + + for (size_t i = 0; i < saltpvdTables.size(); ++i) { + const SaltpvdTable& saltpvdTable = saltpvdTables.getTable(i); + saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn()); + + const auto& cells = reg.cells(i); + for (const auto& cell : cells) { + const double depth = cellCenterDepth_[cell]; + this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); + } + } +} + + +template +void InitialStateComputer:: +updateCellProps_(const GridView& gridView, + const NumericalAquifers& aquifer) +{ + ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); + int numElements = gridView.size(/*codim=*/0); + cellCenterDepth_.resize(numElements); + cellZSpan_.resize(numElements); + cellZMinMax_.resize(numElements); + + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + const auto num_aqu_cells = aquifer.allAquiferCells(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& element = *elemIt; + const unsigned int elemIdx = elemMapper.index(element); + cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); + const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); + cellZSpan_[elemIdx] = Details::cellZSpan(element); + cellZMinMax_[elemIdx] = Details::cellZMinMax(element); + if (!num_aqu_cells.empty()) { + const auto search = num_aqu_cells.find(cartIx); + if (search != num_aqu_cells.end()) { + const auto* aqu_cell = num_aqu_cells.at(cartIx); + const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx]; + cellCenterDepth_[elemIdx] += depth_change_num_aqu; + cellZSpan_[elemIdx].first += depth_change_num_aqu; + cellZSpan_[elemIdx].second += depth_change_num_aqu; + cellZMinMax_[elemIdx].first += depth_change_num_aqu; + cellZMinMax_[elemIdx].second += depth_change_num_aqu; + } + } + } +} + +template +void InitialStateComputer:: +applyNumericalAquifers_(const GridView& gridView, + const NumericalAquifers& aquifer, + const bool co2store) +{ + const auto num_aqu_cells = aquifer.allAquiferCells(); + if (num_aqu_cells.empty()) return; + + ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& element = *elemIt; + const unsigned int elemIdx = elemMapper.index(element); + const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); + const auto search = num_aqu_cells.find(cartIx); + if (search != num_aqu_cells.end()) { + // numerical aquifer cells are filled with water initially + // for co2store the oilphase may be used for the brine + bool co2store_oil_as_brine = co2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const auto watPos = co2store_oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; + if (FluidSystem::phaseIsActive(watPos)) { + this->sat_[watPos][elemIdx] = 1.; + } else { + throw std::logic_error { "Water phase has to be active for numerical aquifer case" }; + } + + const auto oilPos = FluidSystem::oilPhaseIdx; + if (!co2store && FluidSystem::phaseIsActive(oilPos)) { + this->sat_[oilPos][elemIdx] = 0.; + } + + const auto gasPos = FluidSystem::gasPhaseIdx; + if (FluidSystem::phaseIsActive(gasPos)) { + this->sat_[gasPos][elemIdx] = 0.; + } + const auto* aqu_cell = num_aqu_cells.at(cartIx); + const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL " + "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY", + aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id); + OpmLog::info(msg); + + // if pressure is specified for numerical aquifers, we use these pressure values + // for numerical aquifer cells + if (aqu_cell->init_pressure) { + const double pres = *(aqu_cell->init_pressure); + this->pp_[watPos][elemIdx] = pres; + if (FluidSystem::phaseIsActive(gasPos)) { + this->pp_[gasPos][elemIdx] = pres; + } + if (FluidSystem::phaseIsActive(oilPos)) { + this->pp_[oilPos][elemIdx] = pres; + } + } + } + } +} + +template +template +void InitialStateComputer:: +setRegionPvtIdx(const EclipseState& eclState, const RMap& reg) +{ + const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM"); + + for (const auto& r : reg.activeRegions()) { + const auto& cells = reg.cells(r); + regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1; + } +} + +template +template +void InitialStateComputer:: +calcPressSatRsRv(const RMap& reg, + const std::vector& rec, + MaterialLawManager& materialLawManager, + const Comm& comm, + const double grav) +{ + using PhaseSat = Details::PhaseSaturations< + MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId + >; + + auto ptable = Details::PressureTable{ grav }; + auto psat = PhaseSat { materialLawManager, this->swatInit_ }; + auto vspan = std::array{}; + + std::vector regionIsEmpty(rec.size(), 0); + for (size_t r = 0; r < rec.size(); ++r) { + const auto& cells = reg.cells(r); + + Details::verticalExtent(cells, cellZMinMax_, comm, vspan); + + const auto acc = rec[r].initializationTargetAccuracy(); + if (acc > 0) { + throw std::runtime_error { + "Cannot initialise model: Positive item 9 is not supported " + "in EQUIL keyword, record " + std::to_string(r + 1) + }; + } + + if (cells.empty()) { + regionIsEmpty[r] = 1; + continue; + } + + const auto eqreg = EquilReg { + rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r] + }; + + // Ensure gas/oil and oil/water contacts are within the span for the + // phase pressure calculation. + vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc())); + vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc())); + + ptable.equilibrate(eqreg, vspan); + + if (acc == 0) { + // Centre-point method + this->equilibrateCellCentres(cells, eqreg, ptable, psat); + } + else if (acc < 0) { + // Horizontal subdivision + this->equilibrateHorizontal(cells, eqreg, -acc, + ptable, psat); + } else { + // Horizontal subdivision with titled fault blocks + // the simulator throw a few line above for the acc > 0 case + // i.e. we should not reach here. + assert(false); + } + } + comm.min(regionIsEmpty.data(),regionIsEmpty.size()); + if (comm.rank() == 0) { + for (size_t r = 0; r < rec.size(); ++r) { + if (regionIsEmpty[r]) //region is empty on all partitions + OpmLog::warning("Equilibration region " + std::to_string(r + 1) + + " has no active cells"); + } + } +} + +template +template +void InitialStateComputer:: +cellLoop(const CellRange& cells, + EquilibrationMethod&& eqmethod) +{ + const auto oilPos = FluidSystem::oilPhaseIdx; + const auto gasPos = FluidSystem::gasPhaseIdx; + const auto watPos = FluidSystem::waterPhaseIdx; + + const auto oilActive = FluidSystem::phaseIsActive(oilPos); + const auto gasActive = FluidSystem::phaseIsActive(gasPos); + const auto watActive = FluidSystem::phaseIsActive(watPos); + + auto pressures = Details::PhaseQuantityValue{}; + auto saturations = Details::PhaseQuantityValue{}; + auto Rs = 0.0; + auto Rv = 0.0; + auto Rvw = 0.0; + + for (const auto& cell : cells) { + eqmethod(cell, pressures, saturations, Rs, Rv, Rvw); + + if (oilActive) { + this->pp_ [oilPos][cell] = pressures.oil; + this->sat_[oilPos][cell] = saturations.oil; + } + + if (gasActive) { + this->pp_ [gasPos][cell] = pressures.gas; + this->sat_[gasPos][cell] = saturations.gas; + } + + if (watActive) { + this->pp_ [watPos][cell] = pressures.water; + this->sat_[watPos][cell] = saturations.water; + } + + if (oilActive && gasActive) { + this->rs_[cell] = Rs; + this->rv_[cell] = Rv; + } + + if (watActive && gasActive) { + this->rvw_[cell] = Rvw; + } + } +} + +template +template +void InitialStateComputer:: +equilibrateCellCentres(const CellRange& cells, + const EquilReg& eqreg, + const PressTable& ptable, + PhaseSat& psat) +{ + using CellPos = typename PhaseSat::Position; + using CellID = std::remove_cv_t().cell)>>; + this->cellLoop(cells, [this, &eqreg, &ptable, &psat] + (const CellID cell, + Details::PhaseQuantityValue& pressures, + Details::PhaseQuantityValue& saturations, + double& Rs, + double& Rv, + double& Rvw) -> void + { + const auto pos = CellPos { + cell, cellCenterDepth_[cell] + }; + + saturations = psat.deriveSaturations(pos, eqreg, ptable); + pressures = psat.correctedPhasePressures(); + + const auto temp = this->temperature_[cell]; + + Rs = eqreg.dissolutionCalculator() + (pos.depth, pressures.oil, temp, saturations.gas); + + Rv = eqreg.evaporationCalculator() + (pos.depth, pressures.gas, temp, saturations.oil); + + Rvw = eqreg.waterEvaporationCalculator() + (pos.depth, pressures.gas, temp, saturations.water); + }); +} + +template +template +void InitialStateComputer:: +equilibrateHorizontal(const CellRange& cells, + const EquilReg& eqreg, + const int acc, + const PressTable& ptable, + PhaseSat& psat) +{ + using CellPos = typename PhaseSat::Position; + using CellID = std::remove_cv_t().cell)>>; + + this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat] + (const CellID cell, + Details::PhaseQuantityValue& pressures, + Details::PhaseQuantityValue& saturations, + double& Rs, + double& Rv, + double& Rvw) -> void + { + pressures .reset(); + saturations.reset(); + + auto totfrac = 0.0; + for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) { + const auto pos = CellPos { cell, depth }; + + saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac); + pressures .axpy(psat.correctedPhasePressures(), frac); + + totfrac += frac; + } + + if (totfrac > 0.) { + saturations /= totfrac; + pressures /= totfrac; + } else { + // Fall back to centre point method for zero-thickness cells. + const auto pos = CellPos { + cell, cellCenterDepth_[cell] + }; + + saturations = psat.deriveSaturations(pos, eqreg, ptable); + pressures = psat.correctedPhasePressures(); + } + + const auto temp = this->temperature_[cell]; + const auto cz = cellCenterDepth_[cell]; + + Rs = eqreg.dissolutionCalculator() + (cz, pressures.oil, temp, saturations.gas); + + Rv = eqreg.evaporationCalculator() + (cz, pressures.gas, temp, saturations.oil); + + Rvw = eqreg.waterEvaporationCalculator() + (cz, pressures.gas, temp, saturations.water); + }); +} + +} +} // namespace EQUIL +} // namespace Opm +#endif