mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge remote-tracking branch 'upstream/master' into master-refactor-for-cpgrid-support
Manually resolved conflicts in: opm/core/io/eclipse/EclipseWriter.cpp opm/core/io/eclipse/EclipseWriter.hpp opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
This commit is contained in:
Markus Blatt
@@ -97,9 +97,10 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const
|
||||
double* solution,
|
||||
const boost::any& add) const
|
||||
{
|
||||
return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution);
|
||||
return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution, add);
|
||||
}
|
||||
|
||||
void LinearSolverFactory::setTolerance(const double tol)
|
||||
|
||||
@@ -74,7 +74,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const;
|
||||
double* solution,
|
||||
const boost::any& add=boost::any()) const;
|
||||
|
||||
/// Set tolerance for the linear solver.
|
||||
/// \param[in] tol tolerance value
|
||||
|
||||
@@ -20,6 +20,7 @@
|
||||
#ifndef OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED
|
||||
#define OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED
|
||||
|
||||
#include<boost/any.hpp>
|
||||
|
||||
struct CSRMatrix;
|
||||
|
||||
@@ -69,7 +70,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const = 0;
|
||||
double* solution,
|
||||
const boost::any& add=boost::any()) const = 0;
|
||||
|
||||
/// Set tolerance for the linear solver.
|
||||
/// \param[in] tol tolerance value
|
||||
|
||||
@@ -23,6 +23,8 @@
|
||||
#endif
|
||||
|
||||
#include <opm/core/linalg/LinearSolverIstl.hpp>
|
||||
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
// Silence compatibility warning from DUNE headers since we don't use
|
||||
// the deprecated member anyway (in this compilation unit)
|
||||
@@ -35,7 +37,9 @@
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/operators.hh>
|
||||
#include <dune/istl/io.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
#include <dune/istl/paamg/amg.hh>
|
||||
#include <dune/istl/paamg/kamg.hh>
|
||||
@@ -47,7 +51,7 @@
|
||||
|
||||
#include <stdexcept>
|
||||
#include <iostream>
|
||||
|
||||
#include <type_traits>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
@@ -134,7 +138,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const
|
||||
double* solution,
|
||||
const boost::any& comm) const
|
||||
{
|
||||
// Build Istl structures from input.
|
||||
// System matrix
|
||||
@@ -155,10 +160,26 @@ namespace Opm
|
||||
if (maxit == 0) {
|
||||
maxit = 5000;
|
||||
}
|
||||
Dune::SeqScalarProduct<Vector> sp;
|
||||
Dune::Amg::SequentialInformation comm;
|
||||
Operator opA(A);
|
||||
return solveSystem(opA, solution, rhs, sp, comm, maxit);
|
||||
#if HAVE_MPI
|
||||
if(comm.type()==typeid(ParallelISTLInformation))
|
||||
{
|
||||
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
|
||||
const ParallelISTLInformation& info = boost::any_cast<const ParallelISTLInformation&>(comm);
|
||||
Comm istlComm(info.communicator());
|
||||
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
|
||||
Dune::OverlappingSchwarzOperator<Mat,Vector,Vector, Comm>
|
||||
opA(A, istlComm);
|
||||
Dune::OverlappingSchwarzScalarProduct<Vector,Comm> sp(istlComm);
|
||||
return solveSystem(opA, solution, rhs, sp, istlComm, maxit);
|
||||
}
|
||||
else
|
||||
#endif
|
||||
{
|
||||
Dune::SeqScalarProduct<Vector> sp;
|
||||
Dune::Amg::SequentialInformation comm;
|
||||
Operator opA(A);
|
||||
return solveSystem(opA, solution, rhs, sp, comm, maxit);
|
||||
}
|
||||
}
|
||||
|
||||
template<class O, class S, class C>
|
||||
@@ -204,6 +225,14 @@ namespace Opm
|
||||
break;
|
||||
case FastAMG:
|
||||
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
|
||||
|
||||
#if HAVE_MPI
|
||||
if(std::is_same<C,Dune::OwnerOverlapCopyCommunication<int,int> >::value)
|
||||
{
|
||||
OPM_THROW(std::runtime_error, "Trying to use sequential FastAMG solver for a parallel problem!");
|
||||
}
|
||||
#endif // HAVE_MPI
|
||||
|
||||
res = solveFastAMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
|
||||
linsolver_prolongate_factor_);
|
||||
#else
|
||||
@@ -234,28 +263,47 @@ namespace Opm
|
||||
return linsolver_residual_tolerance_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
template<class P, class O, class C>
|
||||
struct SmootherChooser
|
||||
{
|
||||
typedef P Type;
|
||||
};
|
||||
|
||||
#if HAVE_MPI
|
||||
template<class P, class O>
|
||||
struct SmootherChooser<P, O, Dune::OwnerOverlapCopyCommunication<int,int> >
|
||||
{
|
||||
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
|
||||
typedef Dune::BlockPreconditioner<typename O::domain_type, typename O::range_type,
|
||||
Comm, P>
|
||||
Type;
|
||||
};
|
||||
#endif
|
||||
|
||||
|
||||
template<class P, class O, class C>
|
||||
struct PreconditionerTraits
|
||||
{
|
||||
typedef std::shared_ptr<P> PointerType;
|
||||
typedef typename SmootherChooser<P,O,C>::Type SmootherType;
|
||||
typedef std::shared_ptr<SmootherType> PointerType;
|
||||
};
|
||||
|
||||
template<class P, class O, class C>
|
||||
typename PreconditionerTraits<P,O>::PointerType
|
||||
typename PreconditionerTraits<P,O,C>::PointerType
|
||||
makePreconditioner(O& opA, double relax, const C& comm, int iterations=1)
|
||||
{
|
||||
typename Dune::Amg::SmootherTraits<P>::Arguments args;
|
||||
typename Dune::Amg::ConstructionTraits<P>::Arguments cargs;
|
||||
typedef typename SmootherChooser<P,O,C>::Type SmootherType;
|
||||
typedef typename PreconditionerTraits<P,O,C>::PointerType PointerType;
|
||||
typename Dune::Amg::SmootherTraits<SmootherType>::Arguments args;
|
||||
typename Dune::Amg::ConstructionTraits<SmootherType>::Arguments cargs;
|
||||
cargs.setMatrix(opA.getmat());
|
||||
args.iterations=iterations;
|
||||
args.relaxationFactor=relax;
|
||||
cargs.setArgs(args);
|
||||
cargs.setComm(comm);
|
||||
return std::shared_ptr<P>(Dune::Amg::ConstructionTraits<P>::construct(cargs));
|
||||
return PointerType(Dune::Amg::ConstructionTraits<SmootherType>::construct(cargs));
|
||||
}
|
||||
|
||||
template<class O, class S, class C>
|
||||
@@ -323,19 +371,20 @@ namespace Opm
|
||||
#endif
|
||||
|
||||
#if SMOOTHER_ILU
|
||||
typedef Dune::SeqILU0<Mat,Vector,Vector> Smoother;
|
||||
typedef Dune::SeqILU0<Mat,Vector,Vector> SeqSmoother;
|
||||
#else
|
||||
typedef Dune::SeqSOR<Mat,Vector,Vector> Smoother;
|
||||
typedef Dune::SeqSOR<Mat,Vector,Vector> SeqSmoother;
|
||||
#endif
|
||||
typedef typename SmootherChooser<SeqSmoother, O, C>::Type Smoother;
|
||||
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
|
||||
typedef Dune::Amg::AMG<Operator,Vector,Smoother> Precond;
|
||||
typedef Dune::Amg::AMG<O,Vector,Smoother,C> Precond;
|
||||
|
||||
// Construct preconditioner.
|
||||
Criterion criterion;
|
||||
Precond::SmootherArgs smootherArgs;
|
||||
typename Precond::SmootherArgs smootherArgs;
|
||||
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
|
||||
linsolver_smooth_steps);
|
||||
Precond precond(opA, criterion, smootherArgs);
|
||||
Precond precond(opA, criterion, smootherArgs, comm);
|
||||
|
||||
// Construct linear solver.
|
||||
Dune::CGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
|
||||
@@ -360,6 +409,7 @@ namespace Opm
|
||||
double linsolver_prolongate_factor, int linsolver_smooth_steps)
|
||||
{
|
||||
// Solve with AMG solver.
|
||||
Dune::MatrixAdapter<typename O::matrix_type,Vector,Vector> sOpA(opA.getmat());
|
||||
|
||||
#if FIRST_DIAGONAL
|
||||
typedef Dune::Amg::FirstDiagonal CouplingMetric;
|
||||
@@ -386,10 +436,10 @@ namespace Opm
|
||||
Criterion criterion;
|
||||
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
|
||||
linsolver_smooth_steps);
|
||||
Precond precond(opA, criterion, smootherArgs);
|
||||
Precond precond(sOpA, criterion, smootherArgs);
|
||||
|
||||
// Construct linear solver.
|
||||
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
|
||||
Dune::GeneralizedPCGSolver<Vector> linsolve(sOpA, precond, tolerance, maxit, verbosity);
|
||||
|
||||
// Solve system.
|
||||
Dune::InverseOperatorResult result;
|
||||
@@ -409,6 +459,8 @@ namespace Opm
|
||||
double linsolver_prolongate_factor)
|
||||
{
|
||||
// Solve with AMG solver.
|
||||
typedef Dune::MatrixAdapter<typename O::matrix_type, Vector, Vector> Operator;
|
||||
Operator sOpA(opA.getmat());
|
||||
|
||||
#if FIRST_DIAGONAL
|
||||
typedef Dune::Amg::FirstDiagonal CouplingMetric;
|
||||
@@ -434,10 +486,10 @@ namespace Opm
|
||||
parms.setNoPreSmoothSteps(smooth_steps);
|
||||
parms.setNoPostSmoothSteps(smooth_steps);
|
||||
parms.setProlongationDampingFactor(linsolver_prolongate_factor);
|
||||
Precond precond(opA, criterion, parms);
|
||||
Precond precond(sOpA, criterion, parms);
|
||||
|
||||
// Construct linear solver.
|
||||
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
|
||||
Dune::GeneralizedPCGSolver<Vector> linsolve(sOpA, precond, tolerance, maxit, verbosity);
|
||||
|
||||
// Solve system.
|
||||
Dune::InverseOperatorResult result;
|
||||
|
||||
@@ -24,12 +24,11 @@
|
||||
#include <opm/core/linalg/LinearSolverInterface.hpp>
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <string>
|
||||
|
||||
#include <boost/any.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
/// Concrete class encapsulating some dune-istl linear solvers.
|
||||
class LinearSolverIstl : public LinearSolverInterface
|
||||
{
|
||||
@@ -75,7 +74,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const;
|
||||
double* solution,
|
||||
const boost::any& comm=boost::any()) const;
|
||||
|
||||
/// Set tolerance for the residual in dune istl linear solver.
|
||||
/// \param[in] tol tolerance value
|
||||
|
||||
@@ -46,7 +46,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const
|
||||
double* solution,
|
||||
const boost::any&) const
|
||||
{
|
||||
CSRMatrix A = {
|
||||
(size_t)size,
|
||||
|
||||
@@ -55,7 +55,8 @@ namespace Opm
|
||||
const int* ja,
|
||||
const double* sa,
|
||||
const double* rhs,
|
||||
double* solution) const;
|
||||
double* solution,
|
||||
const boost::any& add=boost::any()) const;
|
||||
|
||||
/// Set tolerance for the linear solver.
|
||||
/// \param[in] tol tolerance value
|
||||
|
||||
160
opm/core/linalg/ParallelIstlInformation.hpp
Normal file
160
opm/core/linalg/ParallelIstlInformation.hpp
Normal file
@@ -0,0 +1,160 @@
|
||||
/*
|
||||
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
|
||||
#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
|
||||
|
||||
#if HAVE_MPI && HAVE_DUNE_ISTL
|
||||
|
||||
#include "mpi.h"
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
#include <dune/common/parallel/interface.hh>
|
||||
#include <dune/common/parallel/communicator.hh>
|
||||
#include <dune/common/enumset.hh>
|
||||
|
||||
#include<algorithm>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/// \brief Class that encapsulates the parallelization information needed by the
|
||||
/// ISTL solvers.
|
||||
class ParallelISTLInformation
|
||||
{
|
||||
public:
|
||||
/// \brief The type of the parallel index set used.
|
||||
typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
|
||||
/// \brief The type of the remote indices information used.
|
||||
typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
|
||||
|
||||
/// \brief Constructs an empty parallel information object using MPI_COMM_WORLD
|
||||
ParallelISTLInformation()
|
||||
: indexSet_(new ParallelIndexSet),
|
||||
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
|
||||
communicator_(MPI_COMM_WORLD)
|
||||
{}
|
||||
/// \brief Constructs an empty parallel information object using a communicator.
|
||||
/// \param communicator The communicator to use.
|
||||
ParallelISTLInformation(MPI_Comm communicator)
|
||||
: indexSet_(new ParallelIndexSet),
|
||||
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
|
||||
communicator_(communicator)
|
||||
{}
|
||||
/// \brief Constructs a parallel information object from the specified information.
|
||||
/// \param indexSet The parallel index set to use.
|
||||
/// \param remoteIndices The remote indices information to use.
|
||||
/// \param communicator The communicator to use.
|
||||
ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
|
||||
const std::shared_ptr<RemoteIndices>& remoteIndices,
|
||||
MPI_Comm communicator)
|
||||
: indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
|
||||
{}
|
||||
/// \brief Copy constructor.
|
||||
///
|
||||
/// The information will be shared by the the two objects.
|
||||
ParallelISTLInformation(const ParallelISTLInformation& other)
|
||||
: indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
|
||||
communicator_(other.communicator_)
|
||||
{}
|
||||
/// \brief Get a pointer to the underlying index set.
|
||||
std::shared_ptr<ParallelIndexSet> indexSet() const
|
||||
{
|
||||
return indexSet_;
|
||||
}
|
||||
/// \brief Get a pointer to the remote indices information.
|
||||
std::shared_ptr<RemoteIndices> remoteIndices() const
|
||||
{
|
||||
return remoteIndices_;
|
||||
}
|
||||
/// \brief Get the MPI communicator that we use.
|
||||
MPI_Comm communicator() const
|
||||
{
|
||||
return communicator_;
|
||||
}
|
||||
/// \brief Copy the information stored to the specified objects.
|
||||
/// \param[out] indexSet The object to store the index set in.
|
||||
/// \param[out] remoteIndices The object to store the remote indices information in.
|
||||
void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices) const
|
||||
{
|
||||
indexSet.beginResize();
|
||||
IndexSetInserter<ParallelIndexSet> inserter(indexSet);
|
||||
std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
|
||||
indexSet.endResize();
|
||||
remoteIndices.rebuild<false>();
|
||||
}
|
||||
/// \brief Communcate the dofs owned by us to the other process.
|
||||
///
|
||||
/// Afterwards all associated dofs will contain the same data.
|
||||
template<class T>
|
||||
void copyOwnerToAll (const T& source, T& dest) const
|
||||
{
|
||||
typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
|
||||
typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
|
||||
OwnerOverlapSet sourceFlags;
|
||||
AllSet destFlags;
|
||||
Dune::Interface interface(communicator_);
|
||||
if(!remoteIndices_->isSynced())
|
||||
remoteIndices_->rebuild<false>();
|
||||
interface.build(*remoteIndices_,sourceFlags,destFlags);
|
||||
Dune::BufferedCommunicator communicator;
|
||||
communicator.template build<T>(interface);
|
||||
communicator.template forward<CopyGatherScatter<T> >(source,dest);
|
||||
communicator.free();
|
||||
}
|
||||
private:
|
||||
/** \brief gather/scatter callback for communcation */
|
||||
template<typename T>
|
||||
struct CopyGatherScatter
|
||||
{
|
||||
typedef typename Dune::CommPolicy<T>::IndexedType V;
|
||||
|
||||
static V gather(const T& a, std::size_t i)
|
||||
{
|
||||
return a[i];
|
||||
}
|
||||
|
||||
static void scatter(T& a, V v, std::size_t i)
|
||||
{
|
||||
a[i] = v;
|
||||
}
|
||||
};
|
||||
template<class T>
|
||||
class IndexSetInserter
|
||||
{
|
||||
public:
|
||||
typedef T ParallelIndexSet;
|
||||
|
||||
IndexSetInserter(ParallelIndexSet& indexSet)
|
||||
: indexSet_(&indexSet)
|
||||
{}
|
||||
void operator()(const typename ParallelIndexSet::IndexPair& pair)
|
||||
{
|
||||
indexSet_->add(pair.global(), pair.local());
|
||||
}
|
||||
|
||||
private:
|
||||
ParallelIndexSet* indexSet_;
|
||||
};
|
||||
|
||||
std::shared_ptr<ParallelIndexSet> indexSet_;
|
||||
std::shared_ptr<RemoteIndices> remoteIndices_;
|
||||
MPI_Comm communicator_;
|
||||
};
|
||||
} // end namespace Opm
|
||||
#endif
|
||||
#endif
|
||||
@@ -471,8 +471,7 @@ namespace Opm
|
||||
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
|
||||
}
|
||||
|
||||
|
||||
// Initialize saturation scaling parameter
|
||||
// Initialize saturation scaling parameters
|
||||
template <class SatFuncSet>
|
||||
template <class T>
|
||||
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
|
||||
@@ -1003,8 +1002,9 @@ namespace Opm
|
||||
bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD");
|
||||
bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD");
|
||||
int itab = 0;
|
||||
std::vector<std::vector<std::vector<double> > > table_dummy;
|
||||
std::vector<std::vector<std::vector<double> > >& table = table_dummy;
|
||||
std::vector<std::vector<double> > param_col;
|
||||
std::vector<std::vector<double> > depth_col;
|
||||
std::vector<std::string> col_names;
|
||||
|
||||
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
|
||||
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
|
||||
@@ -1070,22 +1070,14 @@ namespace Opm
|
||||
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
|
||||
}
|
||||
if (!useKeyword && itab > 0) {
|
||||
Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD"));
|
||||
table.resize(1); // only one region supported so far
|
||||
for (unsigned i = 0; i < table.size(); ++i) {
|
||||
table[i].resize(9);
|
||||
|
||||
for (int k = 0; k < enptvd.numRows(); ++k) {
|
||||
table[i][0] = enptvd.getDepthColumn();
|
||||
table[i][1] = enptvd.getSwcoColumn();
|
||||
table[i][2] = enptvd.getSwcritColumn();
|
||||
table[i][3] = enptvd.getSwmaxColumn();
|
||||
table[i][4] = enptvd.getSgcoColumn();
|
||||
table[i][5] = enptvd.getSgcritColumn();
|
||||
table[i][6] = enptvd.getSgmaxColumn();
|
||||
table[i][7] = enptvd.getSowcritColumn();
|
||||
table[i][8] = enptvd.getSogcritColumn();
|
||||
}
|
||||
int num_tables = newParserDeck->getKeyword("ENPTVD")->size();
|
||||
param_col.resize(num_tables);
|
||||
depth_col.resize(num_tables);
|
||||
col_names.resize(9);
|
||||
for (int table_num=0; table_num<num_tables; ++table_num) {
|
||||
Opm::SimpleTable enptvd(newParserDeck->getKeyword("ENPTVD"), col_names, table_num);
|
||||
depth_col[table_num] = enptvd.getColumn(0); // depth
|
||||
param_col[table_num] = enptvd.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
|
||||
}
|
||||
}
|
||||
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
|
||||
@@ -1141,21 +1133,14 @@ namespace Opm
|
||||
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
|
||||
}
|
||||
if (!useKeyword && itab > 0) {
|
||||
Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD"));
|
||||
table.resize(1); // only one region supported so far
|
||||
for (unsigned i = 0; i < table.size(); ++i) {
|
||||
table[i].resize(8);
|
||||
|
||||
for (int k = 0; k < enkrvd.numRows(); ++k) {
|
||||
table[i][0] = enkrvd.getDepthColumn();
|
||||
table[i][1] = enkrvd.getKrwmaxColumn();
|
||||
table[i][2] = enkrvd.getKrgmaxColumn();
|
||||
table[i][3] = enkrvd.getKromaxColumn();
|
||||
table[i][4] = enkrvd.getKrwcritColumn();
|
||||
table[i][5] = enkrvd.getKrgcritColumn();
|
||||
table[i][6] = enkrvd.getKrocritgColumn();
|
||||
table[i][7] = enkrvd.getKrocritwColumn();
|
||||
}
|
||||
int num_tables = newParserDeck->getKeyword("ENKRVD")->size();
|
||||
param_col.resize(num_tables);
|
||||
depth_col.resize(num_tables);
|
||||
col_names.resize(8);
|
||||
for (int table_num=0; table_num<num_tables; ++table_num) {
|
||||
Opm::SimpleTable enkrvd(newParserDeck->getKeyword("ENKRVD"), col_names, table_num);
|
||||
depth_col[table_num] = enkrvd.getColumn(0); // depth
|
||||
param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1172,25 +1157,15 @@ namespace Opm
|
||||
scaleparam[c] = val[deck_pos];
|
||||
}
|
||||
} else {
|
||||
// TODO for new parser
|
||||
/*
|
||||
std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
|
||||
if (keyword[0] == 'S')
|
||||
newParserDeck.getENPTVD().write(std::cout);
|
||||
else
|
||||
newParserDeck.getENKRVD().write(std::cout);
|
||||
*/
|
||||
const int dim = dimensions;
|
||||
for (int cell = 0; cell < number_of_cells; ++cell) {
|
||||
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
|
||||
if (table[itab][jtab][0] != -1.0) {
|
||||
std::vector<double>& depth = table[0][jtab];
|
||||
std::vector<double>& val = table[itab][jtab];
|
||||
if (param_col[jtab][0] >= 0.0) {
|
||||
double zc = UgGridHelpers
|
||||
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
|
||||
dim-1);
|
||||
if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
|
||||
scaleparam[cell] = linearInterpolation(depth, val, zc);
|
||||
if (zc >= depth_col[jtab].front() && zc <= depth_col[jtab].back()) { //don't want extrap outside depth interval
|
||||
scaleparam[cell] = linearInterpolation(depth_col[jtab], param_col[jtab], zc);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
829
opm/core/simulator/EquilibrationHelpers.hpp
Normal file
829
opm/core/simulator/EquilibrationHelpers.hpp
Normal file
@@ -0,0 +1,829 @@
|
||||
/*
|
||||
Copyright 2014 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
#include <opm/core/utility/RegionMapping.hpp>
|
||||
#include <opm/core/utility/RootFinders.hpp>
|
||||
|
||||
#include <memory>
|
||||
|
||||
|
||||
/*
|
||||
---- synopsis of EquilibrationHelpers.hpp ----
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace Equil {
|
||||
|
||||
template <class Props>
|
||||
class DensityCalculator;
|
||||
|
||||
template <>
|
||||
class DensityCalculator< BlackoilPropertiesInterface >;
|
||||
|
||||
namespace Miscibility {
|
||||
class RsFunction;
|
||||
class NoMixing;
|
||||
class RsVD;
|
||||
class RsSatAtContact;
|
||||
}
|
||||
|
||||
struct EquilRecord;
|
||||
|
||||
template <class DensCalc>
|
||||
class EquilReg;
|
||||
|
||||
|
||||
struct PcEq;
|
||||
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc,
|
||||
const bool increasing = false);
|
||||
struct PcEqSum
|
||||
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc);
|
||||
} // namespace Equil
|
||||
} // namespace Opm
|
||||
|
||||
---- end of synopsis of EquilibrationHelpers.hpp ----
|
||||
*/
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
/**
|
||||
* Types and routines that collectively implement a basic
|
||||
* ECLIPSE-style equilibration-based initialisation scheme.
|
||||
*
|
||||
* This namespace is intentionally nested to avoid name clashes
|
||||
* with other parts of OPM.
|
||||
*/
|
||||
namespace Equil {
|
||||
|
||||
|
||||
template <class Props>
|
||||
class DensityCalculator;
|
||||
|
||||
/**
|
||||
* Facility for calculating phase densities based on the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* Implements the crucial <CODE>operator()(p,svol)</CODE>
|
||||
* function that is expected by class EquilReg.
|
||||
*/
|
||||
template <>
|
||||
class DensityCalculator< BlackoilPropertiesInterface > {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props Implementation of the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* \param[in] c Single cell used as a representative cell
|
||||
* in a PVT region.
|
||||
*/
|
||||
DensityCalculator(const BlackoilPropertiesInterface& props,
|
||||
const int c)
|
||||
: props_(props)
|
||||
, c_(1, c)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute phase densities of all phases at phase point
|
||||
* given by (pressure, surface volume) tuple.
|
||||
*
|
||||
* \param[in] p Fluid pressure.
|
||||
*
|
||||
* \param[in] z Surface volumes of all phases.
|
||||
*
|
||||
* \return Phase densities at phase point.
|
||||
*/
|
||||
std::vector<double>
|
||||
operator()(const double p,
|
||||
const std::vector<double>& z) const
|
||||
{
|
||||
const int np = props_.numPhases();
|
||||
std::vector<double> A(np * np, 0);
|
||||
|
||||
assert (z.size() == std::vector<double>::size_type(np));
|
||||
|
||||
double* dAdp = 0;
|
||||
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
|
||||
|
||||
std::vector<double> rho(np, 0.0);
|
||||
props_.density(1, &A[0], &rho[0]);
|
||||
|
||||
return rho;
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const std::vector<int> c_;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Types and routines relating to phase mixing in
|
||||
* equilibration calculations.
|
||||
*/
|
||||
namespace Miscibility {
|
||||
|
||||
/**
|
||||
* Base class for phase mixing functions.
|
||||
*/
|
||||
class RsFunction
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Function call operator.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
virtual double operator()(const double depth,
|
||||
const double press,
|
||||
const double sat = 0.0) const = 0;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Type that implements "no phase mixing" policy.
|
||||
*/
|
||||
class NoMixing : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press. In "no mixing
|
||||
* policy", this is identically zero.
|
||||
*/
|
||||
double
|
||||
operator()(const double /* depth */,
|
||||
const double /* press */,
|
||||
const double sat = 0.0) const
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Type that implements "dissolved gas-oil ratio"
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RSVD'.
|
||||
*/
|
||||
class RsVD : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rs Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RsVD(const BlackoilPropertiesInterface& props,
|
||||
const int cell,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: props_(props)
|
||||
, cell_(cell)
|
||||
, depth_(depth)
|
||||
, rs_(rs)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double sat_gas = 0.0) const
|
||||
{
|
||||
if (sat_gas > 0.0) {
|
||||
return satRs(press);
|
||||
} else {
|
||||
return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rs/Bo is in the gas row and oil column of A_.
|
||||
// 1/Bo is in the oil row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*opos + gpos] / A_[np*opos + opos];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Type that implements "vaporized oil-gas ratio"
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RVVD'.
|
||||
*/
|
||||
class RvVD : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rv Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RvVD(const BlackoilPropertiesInterface& props,
|
||||
const int cell,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rv)
|
||||
: props_(props)
|
||||
, cell_(cell)
|
||||
, depth_(depth)
|
||||
, rv_(rv)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RV
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RV
|
||||
* value.
|
||||
*
|
||||
* \return Vaporized oil-gas ratio (RV) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double sat_oil = 0.0 ) const
|
||||
{
|
||||
if (sat_oil > 0.0) {
|
||||
return satRv(press);
|
||||
} else {
|
||||
return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rv_; /**< Vaporized oil-gas ratio */
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRv(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rv/Bg is in the oil row and gas column of A_.
|
||||
// 1/Bg is in the gas row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*gpos + opos] / A_[np*gpos + gpos];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Class that implements "dissolved gas-oil ratio" (Rs)
|
||||
* as function of depth and pressure as follows:
|
||||
*
|
||||
* 1. The Rs at the gas-oil contact is equal to the
|
||||
* saturated Rs value, Rs_sat_contact.
|
||||
*
|
||||
* 2. The Rs elsewhere is equal to Rs_sat_contact, but
|
||||
* constrained to the saturated value as given by the
|
||||
* local pressure.
|
||||
*
|
||||
* This should yield Rs-values that are constant below the
|
||||
* contact, and decreasing above the contact.
|
||||
*/
|
||||
class RsSatAtContact : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
*/
|
||||
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
|
||||
: props_(props), cell_(cell)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||
rs_sat_contact_ = satRs(p_contact);
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double /* depth */,
|
||||
const double press,
|
||||
const double sat_gas = 0.0) const
|
||||
{
|
||||
if (sat_gas > 0.0) {
|
||||
return satRs(press);
|
||||
} else {
|
||||
return std::min(satRs(press), rs_sat_contact_);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
double rs_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rs/Bo is in the gas row and oil column of A_.
|
||||
// 1/Bo is in the oil row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*opos + gpos] / A_[np*opos + opos];
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Class that implements "vaporized oil-gas ratio" (Rv)
|
||||
* as function of depth and pressure as follows:
|
||||
*
|
||||
* 1. The Rv at the gas-oil contact is equal to the
|
||||
* saturated Rv value, Rv_sat_contact.
|
||||
*
|
||||
* 2. The Rv elsewhere is equal to Rv_sat_contact, but
|
||||
* constrained to the saturated value as given by the
|
||||
* local pressure.
|
||||
*
|
||||
* This should yield Rv-values that are constant below the
|
||||
* contact, and decreasing above the contact.
|
||||
*/
|
||||
class RvSatAtContact : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
*/
|
||||
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
|
||||
: props_(props), cell_(cell)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
|
||||
rv_sat_contact_ = satRv(p_contact);
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RV
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure at which to calculate RV
|
||||
* value.
|
||||
*
|
||||
* \return Dissolved oil-gas ratio (RV) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double /*depth*/,
|
||||
const double press,
|
||||
const double sat_oil = 0.0) const
|
||||
{
|
||||
if (sat_oil > 0.0) {
|
||||
return satRv(press);
|
||||
} else {
|
||||
return std::min(satRv(press), rv_sat_contact_);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
double rv_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRv(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rv/Bg is in the oil row and gas column of A_.
|
||||
// 1/Bg is in the gas row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*gpos + opos] / A_[np*gpos + gpos];
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Miscibility
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Equilibration record.
|
||||
*
|
||||
* Layout and contents inspired by first six items of
|
||||
* ECLIPSE's 'EQUIL' records. This is the minimum amount of
|
||||
* input data needed to define phase pressures in an
|
||||
* equilibration region.
|
||||
*
|
||||
* Data consists of three pairs of depth and pressure values:
|
||||
* 1. main
|
||||
* - @c depth Main datum depth.
|
||||
* - @c press Pressure at datum depth.
|
||||
*
|
||||
* 2. woc
|
||||
* - @c depth Depth of water-oil contact
|
||||
* - @c press water-oil capillary pressure at water-oil contact.
|
||||
* Capillary pressure defined as "P_oil - P_water".
|
||||
*
|
||||
* 3. goc
|
||||
* - @c depth Depth of gas-oil contact
|
||||
* - @c press Gas-oil capillary pressure at gas-oil contact.
|
||||
* Capillary pressure defined as "P_gas - P_oil".
|
||||
*/
|
||||
struct EquilRecord {
|
||||
struct {
|
||||
double depth;
|
||||
double press;
|
||||
} main, woc, goc;
|
||||
int live_oil_table_index;
|
||||
int wet_gas_table_index;
|
||||
int N;
|
||||
};
|
||||
|
||||
/**
|
||||
* Aggregate information base of an equilibration region.
|
||||
*
|
||||
* Provides inquiry methods for retrieving depths of contacs
|
||||
* and pressure values as well as a means of calculating fluid
|
||||
* densities, dissolved gas-oil ratio and vapourised oil-gas
|
||||
* ratios.
|
||||
*
|
||||
* \tparam DensCalc Type that provides access to a phase
|
||||
* density calculation facility. Must implement an operator()
|
||||
* declared as
|
||||
* <CODE>
|
||||
* std::vector<double>
|
||||
* operator()(const double press,
|
||||
* const std::vector<double>& svol )
|
||||
* </CODE>
|
||||
* that calculates the phase densities of all phases in @c
|
||||
* svol at fluid pressure @c press.
|
||||
*/
|
||||
template <class DensCalc>
|
||||
class EquilReg {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] rec Equilibration data of current region.
|
||||
* \param[in] density Density calculator of current region.
|
||||
* \param[in] rs Calculator of dissolved gas-oil ratio.
|
||||
* \param[in] rv Calculator of vapourised oil-gas ratio.
|
||||
* \param[in] pu Summary of current active phases.
|
||||
*/
|
||||
EquilReg(const EquilRecord& rec,
|
||||
const DensCalc& density,
|
||||
std::shared_ptr<Miscibility::RsFunction> rs,
|
||||
std::shared_ptr<Miscibility::RsFunction> rv,
|
||||
const PhaseUsage& pu)
|
||||
: rec_ (rec)
|
||||
, density_(density)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, pu_ (pu)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Type of density calculator.
|
||||
*/
|
||||
typedef DensCalc CalcDensity;
|
||||
|
||||
/**
|
||||
* Type of dissolved gas-oil ratio calculator.
|
||||
*/
|
||||
typedef Miscibility::RsFunction CalcDissolution;
|
||||
|
||||
/**
|
||||
* Type of vapourised oil-gas ratio calculator.
|
||||
*/
|
||||
typedef Miscibility::RsFunction CalcEvaporation;
|
||||
|
||||
/**
|
||||
* Datum depth in current region
|
||||
*/
|
||||
double datum() const { return this->rec_.main.depth; }
|
||||
|
||||
/**
|
||||
* Pressure at datum depth in current region.
|
||||
*/
|
||||
double pressure() const { return this->rec_.main.press; }
|
||||
|
||||
/**
|
||||
* Depth of water-oil contact.
|
||||
*/
|
||||
double zwoc() const { return this->rec_.woc .depth; }
|
||||
|
||||
/**
|
||||
* water-oil capillary pressure at water-oil contact.
|
||||
*
|
||||
* \return P_o - P_w at WOC.
|
||||
*/
|
||||
double pcow_woc() const { return this->rec_.woc .press; }
|
||||
|
||||
/**
|
||||
* Depth of gas-oil contact.
|
||||
*/
|
||||
double zgoc() const { return this->rec_.goc .depth; }
|
||||
|
||||
/**
|
||||
* Gas-oil capillary pressure at gas-oil contact.
|
||||
*
|
||||
* \return P_g - P_o at GOC.
|
||||
*/
|
||||
double pcgo_goc() const { return this->rec_.goc .press; }
|
||||
|
||||
/**
|
||||
* Retrieve phase density calculator of current region.
|
||||
*/
|
||||
const CalcDensity&
|
||||
densityCalculator() const { return this->density_; }
|
||||
|
||||
/**
|
||||
* Retrieve dissolved gas-oil ratio calculator of current
|
||||
* region.
|
||||
*/
|
||||
const CalcDissolution&
|
||||
dissolutionCalculator() const { return *this->rs_; }
|
||||
|
||||
/**
|
||||
* Retrieve vapourised oil-gas ratio calculator of current
|
||||
* region.
|
||||
*/
|
||||
const CalcEvaporation&
|
||||
evaporationCalculator() const { return *this->rv_; }
|
||||
|
||||
/**
|
||||
* Retrieve active fluid phase summary.
|
||||
*/
|
||||
const PhaseUsage&
|
||||
phaseUsage() const { return this->pu_; }
|
||||
|
||||
private:
|
||||
EquilRecord rec_; /**< Equilibration data */
|
||||
DensCalc density_; /**< Density calculator */
|
||||
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
|
||||
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
|
||||
PhaseUsage pu_; /**< Active phase summary */
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Functor for inverting capillary pressure function.
|
||||
/// Function represented is
|
||||
/// f(s) = pc(s) - target_pc
|
||||
struct PcEq
|
||||
{
|
||||
PcEq(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase_] = s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Compute saturation of some phase corresponding to a given
|
||||
/// capillary pressure.
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc,
|
||||
const bool increasing = false)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
|
||||
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
|
||||
|
||||
// Create the equation f(s) = pc(s) - target_pc
|
||||
const PcEq f(props, phase, cell, target_pc);
|
||||
const double f0 = f(s0);
|
||||
const double f1 = f(s1);
|
||||
if (f0 <= 0.0) {
|
||||
return s0;
|
||||
} else if (f1 > 0.0) {
|
||||
return s1;
|
||||
} else {
|
||||
const int max_iter = 30;
|
||||
const double tol = 1e-6;
|
||||
int iter_used = -1;
|
||||
typedef RegulaFalsi<ThrowOnError> ScalarSolver;
|
||||
const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
|
||||
return sol;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Functor for inverting a sum of capillary pressure functions.
|
||||
/// Function represented is
|
||||
/// f(s) = pc1(s) + pc2(1 - s) - target_pc
|
||||
struct PcEqSum
|
||||
{
|
||||
PcEqSum(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase1_] = s;
|
||||
s_[phase2_] = 1.0 - s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase1_] + pc_[phase2_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase1_;
|
||||
const int phase2_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/// Compute saturation of some phase corresponding to a given
|
||||
/// capillary pressure, where the capillary pressure function
|
||||
/// is given as a sum of two other functions.
|
||||
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double smin = sminarr[phase1];
|
||||
const double smax = smaxarr[phase1];
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
|
||||
const PcEqSum f(props, phase1, phase2, cell, target_pc);
|
||||
const double f0 = f(smin);
|
||||
const double f1 = f(smax);
|
||||
if (f0 <= 0.0) {
|
||||
return smin;
|
||||
} else if (f1 > 0.0) {
|
||||
return smax;
|
||||
} else {
|
||||
const int max_iter = 30;
|
||||
const double tol = 1e-6;
|
||||
int iter_used = -1;
|
||||
typedef RegulaFalsi<ThrowOnError> ScalarSolver;
|
||||
const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
|
||||
return sol;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Equil
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
673
opm/core/simulator/initStateEquil.hpp
Normal file
673
opm/core/simulator/initStateEquil.hpp
Normal file
@@ -0,0 +1,673 @@
|
||||
/*
|
||||
Copyright 2014 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
||||
#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/simulator/EquilibrationHelpers.hpp>
|
||||
#include <opm/core/simulator/BlackoilState.hpp>
|
||||
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/utility/RegionMapping.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
|
||||
#include <opm/parser/eclipse/Utility/SimpleTable.hpp>
|
||||
|
||||
#include <array>
|
||||
#include <cassert>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Facilities for an ECLIPSE-style equilibration-based
|
||||
* initialisation scheme (keyword 'EQUIL').
|
||||
*/
|
||||
struct UnstructuredGrid;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/**
|
||||
* Compute initial state by an equilibration procedure.
|
||||
*
|
||||
* The following state fields are modified:
|
||||
* pressure(),
|
||||
* saturation(),
|
||||
* surfacevol(),
|
||||
* gasoilratio(),
|
||||
* rv().
|
||||
*
|
||||
* \param[in] grid Grid.
|
||||
* \param[in] props Property object, pvt and capillary properties are used.
|
||||
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
|
||||
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
|
||||
*/
|
||||
void initStateEquil(const UnstructuredGrid& grid,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const EclipseGridParser& deck,
|
||||
const double gravity,
|
||||
BlackoilState& state);
|
||||
|
||||
|
||||
void initStateEquil(const UnstructuredGrid& grid,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const Opm::DeckConstPtr newParserDeck,
|
||||
const double gravity,
|
||||
BlackoilState& state);
|
||||
|
||||
|
||||
/**
|
||||
* Types and routines that collectively implement a basic
|
||||
* ECLIPSE-style equilibration-based initialisation scheme.
|
||||
*
|
||||
* This namespace is intentionally nested to avoid name clashes
|
||||
* with other parts of OPM.
|
||||
*/
|
||||
namespace Equil {
|
||||
|
||||
/**
|
||||
* Compute initial phase pressures by means of equilibration.
|
||||
*
|
||||
* This function uses the information contained in an
|
||||
* equilibration record (i.e., depths and pressurs) as well as
|
||||
* a density calculator and related data to vertically
|
||||
* integrate the phase pressure ODE
|
||||
* \f[
|
||||
* \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} =
|
||||
* \rho_{\alpha}(z,p_{\alpha})\cdot g
|
||||
* \f]
|
||||
* in which \f$\rho_{\alpha}$ denotes the fluid density of
|
||||
* fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the
|
||||
* corresponding phase pressure, \f$z\f$ is the depth and
|
||||
* \f$g\f$ is the acceleration due to gravity (assumed
|
||||
* directed downwords, in the positive \f$z\f$ direction).
|
||||
*
|
||||
* \tparam Region Type of an equilibration region information
|
||||
* base. Typically an instance of the EquilReg
|
||||
* class template.
|
||||
*
|
||||
* \tparam CellRange Type of cell range that demarcates the
|
||||
* cells pertaining to the current
|
||||
* equilibration region. Must implement
|
||||
* methods begin() and end() to bound the range
|
||||
* as well as provide an inner type,
|
||||
* const_iterator, to traverse the range.
|
||||
*
|
||||
* \param[in] G Grid.
|
||||
* \param[in] reg Current equilibration region.
|
||||
* \param[in] cells Range that spans the cells of the current
|
||||
* equilibration region.
|
||||
* \param[in] grav Acceleration of gravity.
|
||||
*
|
||||
* \return Phase pressures, one vector for each active phase,
|
||||
* of pressure values in each cell in the current
|
||||
* equilibration region.
|
||||
*/
|
||||
template <class Region, class CellRange>
|
||||
std::vector< std::vector<double> >
|
||||
phasePressures(const UnstructuredGrid& G,
|
||||
const Region& reg,
|
||||
const CellRange& cells,
|
||||
const double grav = unit::gravity);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial phase saturations by means of equilibration.
|
||||
*
|
||||
* \tparam Region Type of an equilibration region information
|
||||
* base. Typically an instance of the EquilReg
|
||||
* class template.
|
||||
*
|
||||
* \tparam CellRange Type of cell range that demarcates the
|
||||
* cells pertaining to the current
|
||||
* equilibration region. Must implement
|
||||
* methods begin() and end() to bound the range
|
||||
* as well as provide an inner type,
|
||||
* const_iterator, to traverse the range.
|
||||
*
|
||||
* \param[in] reg Current equilibration region.
|
||||
* \param[in] cells Range that spans the cells of the current
|
||||
* equilibration region.
|
||||
* \param[in] props Property object, needed for capillary functions.
|
||||
* \param[in] phase_pressures Phase pressures, one vector for each active phase,
|
||||
* of pressure values in each cell in the current
|
||||
* equilibration region.
|
||||
* \return Phase saturations, one vector for each phase, each containing
|
||||
* one saturation value per cell in the region.
|
||||
*/
|
||||
template <class Region, class CellRange>
|
||||
std::vector< std::vector<double> >
|
||||
phaseSaturations(const Region& reg,
|
||||
const CellRange& cells,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const std::vector< std::vector<double> >& phase_pressures);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial Rs values.
|
||||
*
|
||||
* \tparam CellRangeType Type of cell range that demarcates the
|
||||
* cells pertaining to the current
|
||||
* equilibration region. Must implement
|
||||
* methods begin() and end() to bound the range
|
||||
* as well as provide an inner type,
|
||||
* const_iterator, to traverse the range.
|
||||
*
|
||||
* \param[in] grid Grid.
|
||||
* \param[in] cells Range that spans the cells of the current
|
||||
* equilibration region.
|
||||
* \param[in] oil_pressure Oil pressure for each cell in range.
|
||||
* \param[in] rs_func Rs as function of pressure and depth.
|
||||
* \return Rs values, one for each cell in the 'cells' range.
|
||||
*/
|
||||
template <class CellRangeType>
|
||||
std::vector<double> computeRs(const UnstructuredGrid& grid,
|
||||
const CellRangeType& cells,
|
||||
const std::vector<double> oil_pressure,
|
||||
const Miscibility::RsFunction& rs_func,
|
||||
const std::vector<double> gas_saturation);
|
||||
|
||||
namespace DeckDependent {
|
||||
|
||||
inline
|
||||
std::vector<EquilRecord>
|
||||
getEquil(const EclipseGridParser& deck)
|
||||
{
|
||||
if (deck.hasField("EQUIL")) {
|
||||
const EQUIL& eql = deck.getEQUIL();
|
||||
|
||||
typedef std::vector<EquilLine>::size_type sz_t;
|
||||
const sz_t nrec = eql.equil.size();
|
||||
|
||||
std::vector<EquilRecord> ret;
|
||||
ret.reserve(nrec);
|
||||
for (sz_t r = 0; r < nrec; ++r) {
|
||||
const EquilLine& rec = eql.equil[r];
|
||||
|
||||
EquilRecord record =
|
||||
{
|
||||
{ rec.datum_depth_ ,
|
||||
rec.datum_depth_pressure_ }
|
||||
,
|
||||
{ rec.water_oil_contact_depth_ ,
|
||||
rec.oil_water_cap_pressure_ }
|
||||
,
|
||||
{ rec.gas_oil_contact_depth_ ,
|
||||
rec.gas_oil_cap_pressure_ }
|
||||
,
|
||||
rec.live_oil_table_index_
|
||||
,
|
||||
rec.wet_gas_table_index_
|
||||
,
|
||||
rec.N_
|
||||
};
|
||||
if (record.N != 0) {
|
||||
OPM_THROW(std::domain_error,
|
||||
"kw EQUIL, item 9: Only N=0 supported.");
|
||||
}
|
||||
ret.push_back(record);
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
else {
|
||||
OPM_THROW(std::domain_error,
|
||||
"Deck does not provide equilibration data.");
|
||||
}
|
||||
}
|
||||
|
||||
inline
|
||||
std::vector<EquilRecord>
|
||||
getEquil(const Opm::DeckConstPtr newParserDeck)
|
||||
{
|
||||
if (newParserDeck->hasKeyword("EQUIL")) {
|
||||
|
||||
Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
|
||||
|
||||
const int nrec = eql.numRegions();
|
||||
|
||||
std::vector<EquilRecord> ret;
|
||||
ret.reserve(nrec);
|
||||
for (int r = 0; r < nrec; ++r) {
|
||||
|
||||
EquilRecord record =
|
||||
{
|
||||
{ eql.datumDepth(r) ,
|
||||
eql.datumDepthPressure(r) }
|
||||
,
|
||||
{ eql.waterOilContactDepth(r) ,
|
||||
eql.waterOilContactCapillaryPressure(r) }
|
||||
,
|
||||
{ eql.gasOilContactDepth(r) ,
|
||||
eql.gasOilContactCapillaryPressure(r) }
|
||||
,
|
||||
eql.liveOilInitProceedure(r)
|
||||
,
|
||||
eql.wetGasInitProceedure(r)
|
||||
,
|
||||
eql.initializationTargetAccuracy(r)
|
||||
};
|
||||
if (record.N != 0) {
|
||||
OPM_THROW(std::domain_error,
|
||||
"kw EQUIL, item 9: Only N=0 supported.");
|
||||
}
|
||||
ret.push_back(record);
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
else {
|
||||
OPM_THROW(std::domain_error,
|
||||
"Deck does not provide equilibration data.");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline
|
||||
std::vector<int>
|
||||
equilnum(const EclipseGridParser& deck,
|
||||
const UnstructuredGrid& G )
|
||||
{
|
||||
std::vector<int> eqlnum;
|
||||
if (deck.hasField("EQLNUM")) {
|
||||
const std::vector<int>& e = deck.getIntegerValue("EQLNUM");
|
||||
eqlnum.reserve(e.size());
|
||||
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
|
||||
std::bind2nd(std::minus<int>(), 1));
|
||||
}
|
||||
else {
|
||||
// No explicit equilibration region.
|
||||
// All cells in region zero.
|
||||
eqlnum.assign(G.number_of_cells, 0);
|
||||
}
|
||||
|
||||
return eqlnum;
|
||||
}
|
||||
|
||||
|
||||
inline
|
||||
std::vector<int>
|
||||
equilnum(const Opm::DeckConstPtr newParserDeck,
|
||||
const UnstructuredGrid& G )
|
||||
{
|
||||
std::vector<int> eqlnum;
|
||||
if (newParserDeck->hasKeyword("EQLNUM")) {
|
||||
const std::vector<int>& e = newParserDeck->getKeyword("EQLNUM")->getIntData();
|
||||
eqlnum.reserve(e.size());
|
||||
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
|
||||
std::bind2nd(std::minus<int>(), 1));
|
||||
}
|
||||
else {
|
||||
// No explicit equilibration region.
|
||||
// All cells in region zero.
|
||||
eqlnum.assign(G.number_of_cells, 0);
|
||||
}
|
||||
|
||||
return eqlnum;
|
||||
}
|
||||
|
||||
|
||||
template <class InputDeck>
|
||||
class InitialStateComputer;
|
||||
|
||||
template <>
|
||||
class InitialStateComputer<Opm::EclipseGridParser> {
|
||||
public:
|
||||
InitialStateComputer(const BlackoilPropertiesInterface& props,
|
||||
const EclipseGridParser& deck ,
|
||||
const UnstructuredGrid& G ,
|
||||
const double grav = unit::gravity)
|
||||
: pp_(props.numPhases(),
|
||||
std::vector<double>(G.number_of_cells)),
|
||||
sat_(props.numPhases(),
|
||||
std::vector<double>(G.number_of_cells)),
|
||||
rs_(G.number_of_cells),
|
||||
rv_(G.number_of_cells)
|
||||
{
|
||||
// Get the equilibration records.
|
||||
const std::vector<EquilRecord> rec = getEquil(deck);
|
||||
|
||||
// Create (inverse) region mapping.
|
||||
const RegionMapping<> eqlmap(equilnum(deck, G));
|
||||
|
||||
// Create Rs functions.
|
||||
rs_func_.reserve(rec.size());
|
||||
if (deck.hasField("DISGAS")) {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
if (rec[i].live_oil_table_index > 0) {
|
||||
if (deck.hasField("RSVD")) {
|
||||
// TODO When this kw is actually parsed, also check for proper number of available tables
|
||||
// For now, just use dummy ...
|
||||
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
|
||||
std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
|
||||
}
|
||||
} else {
|
||||
if (rec[i].goc.depth != rec[i].main.depth) {
|
||||
OPM_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 " << (i + 1) << " (counting from 1), this does not hold.");
|
||||
}
|
||||
const double p_contact = rec[i].main.press;
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
|
||||
}
|
||||
}
|
||||
|
||||
rv_func_.reserve(rec.size());
|
||||
if (deck.hasField("VAPOIL")) {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
if (rec[i].wet_gas_table_index > 0) {
|
||||
if (deck.hasField("RVVD")) {
|
||||
// TODO When this kw is actually parsed, also check for proper number of available tables
|
||||
// For now, just use dummy ...
|
||||
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
|
||||
std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
|
||||
}
|
||||
} else {
|
||||
if (rec[i].goc.depth != rec[i].main.depth) {
|
||||
OPM_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 " << (i + 1) << " (counting from 1), this does not hold.");
|
||||
}
|
||||
const double p_contact = rec[i].main.press + rec[i].goc.press;
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
|
||||
}
|
||||
}
|
||||
|
||||
// Compute pressures, saturations, rs and rv factors.
|
||||
calcPressSatRsRv(eqlmap, rec, props, G, grav);
|
||||
|
||||
// Modify oil pressure in no-oil regions so that the pressures of present phases can
|
||||
// be recovered from the oil pressure and capillary relations.
|
||||
}
|
||||
|
||||
typedef std::vector<double> Vec;
|
||||
typedef std::vector<Vec> PVec; // One per phase.
|
||||
|
||||
const PVec& press() const { return pp_; }
|
||||
const PVec& saturation() const { return sat_; }
|
||||
const Vec& rs() const { return rs_; }
|
||||
const Vec& rv() const { return rv_; }
|
||||
|
||||
private:
|
||||
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
|
||||
typedef EquilReg<RhoCalc> EqReg;
|
||||
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
|
||||
|
||||
PVec pp_;
|
||||
PVec sat_;
|
||||
Vec rs_;
|
||||
Vec rv_;
|
||||
|
||||
template <class RMap>
|
||||
void
|
||||
calcPressSatRsRv(const RMap& reg ,
|
||||
const std::vector< EquilRecord >& rec ,
|
||||
const Opm::BlackoilPropertiesInterface& props,
|
||||
const UnstructuredGrid& G ,
|
||||
const double grav)
|
||||
{
|
||||
typedef Miscibility::NoMixing NoMix;
|
||||
|
||||
for (typename RMap::RegionId
|
||||
r = 0, nr = reg.numRegions();
|
||||
r < nr; ++r)
|
||||
{
|
||||
const typename RMap::CellRange cells = reg.cells(r);
|
||||
|
||||
const int repcell = *cells.begin();
|
||||
const RhoCalc calc(props, repcell);
|
||||
const EqReg eqreg(rec[r], calc,
|
||||
rs_func_[r], rv_func_[r],
|
||||
props.phaseUsage());
|
||||
|
||||
PVec press = phasePressures(G, eqreg, cells, grav);
|
||||
|
||||
const PVec sat = phaseSaturations(eqreg, cells, props, press);
|
||||
|
||||
const int np = props.numPhases();
|
||||
for (int p = 0; p < np; ++p) {
|
||||
copyFromRegion(press[p], cells, pp_[p]);
|
||||
copyFromRegion(sat[p], cells, sat_[p]);
|
||||
}
|
||||
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
|
||||
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
|
||||
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
|
||||
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
|
||||
copyFromRegion(rs, cells, rs_);
|
||||
copyFromRegion(rv, cells, rv_);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class CellRangeType>
|
||||
void copyFromRegion(const Vec& source,
|
||||
const CellRangeType& cells,
|
||||
Vec& destination)
|
||||
{
|
||||
auto s = source.begin();
|
||||
auto c = cells.begin();
|
||||
const auto e = cells.end();
|
||||
for (; c != e; ++c, ++s) {
|
||||
destination[*c] = *s;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
template <>
|
||||
class InitialStateComputer<Opm::DeckConstPtr> {
|
||||
public:
|
||||
InitialStateComputer(const BlackoilPropertiesInterface& props,
|
||||
const Opm::DeckConstPtr newParserDeck,
|
||||
const UnstructuredGrid& G ,
|
||||
const double grav = unit::gravity)
|
||||
: pp_(props.numPhases(),
|
||||
std::vector<double>(G.number_of_cells)),
|
||||
sat_(props.numPhases(),
|
||||
std::vector<double>(G.number_of_cells)),
|
||||
rs_(G.number_of_cells),
|
||||
rv_(G.number_of_cells)
|
||||
{
|
||||
// Get the equilibration records.
|
||||
const std::vector<EquilRecord> rec = getEquil(newParserDeck);
|
||||
|
||||
// Create (inverse) region mapping.
|
||||
const RegionMapping<> eqlmap(equilnum(newParserDeck, G));
|
||||
|
||||
// Create Rs functions.
|
||||
rs_func_.reserve(rec.size());
|
||||
if (newParserDeck->hasKeyword("DISGAS")) {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
if (rec[i].live_oil_table_index > 0) {
|
||||
if (newParserDeck->hasKeyword("RSVD") && rec[i].live_oil_table_index <= newParserDeck->getKeyword("RSVD")->size()) {
|
||||
Opm::SimpleTable rsvd(newParserDeck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
|
||||
std::vector<double> vd(rsvd.getColumn("vd"));
|
||||
std::vector<double> rs(rsvd.getColumn("rs"));
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
|
||||
}
|
||||
} else {
|
||||
if (rec[i].goc.depth != rec[i].main.depth) {
|
||||
OPM_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 " << (i + 1) << " (counting from 1), this does not hold.");
|
||||
}
|
||||
const double p_contact = rec[i].main.press;
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
|
||||
}
|
||||
}
|
||||
|
||||
rv_func_.reserve(rec.size());
|
||||
if (newParserDeck->hasKeyword("VAPOIL")) {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
if (rec[i].wet_gas_table_index > 0) {
|
||||
if (newParserDeck->hasKeyword("RVVD") && rec[i].wet_gas_table_index <= newParserDeck->getKeyword("RVVD")->size()) {
|
||||
Opm::SimpleTable rvvd(newParserDeck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);
|
||||
std::vector<double> vd(rvvd.getColumn("vd"));
|
||||
std::vector<double> rv(rvvd.getColumn("rv"));
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, vd, rv));
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
|
||||
}
|
||||
} else {
|
||||
if (rec[i].goc.depth != rec[i].main.depth) {
|
||||
OPM_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 " << (i + 1) << " (counting from 1), this does not hold.");
|
||||
}
|
||||
const double p_contact = rec[i].main.press + rec[i].goc.press;
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
|
||||
}
|
||||
}
|
||||
|
||||
// Compute pressures, saturations, rs and rv factors.
|
||||
calcPressSatRsRv(eqlmap, rec, props, G, grav);
|
||||
|
||||
// Modify oil pressure in no-oil regions so that the pressures of present phases can
|
||||
// be recovered from the oil pressure and capillary relations.
|
||||
}
|
||||
|
||||
typedef std::vector<double> Vec;
|
||||
typedef std::vector<Vec> PVec; // One per phase.
|
||||
|
||||
const PVec& press() const { return pp_; }
|
||||
const PVec& saturation() const { return sat_; }
|
||||
const Vec& rs() const { return rs_; }
|
||||
const Vec& rv() const { return rv_; }
|
||||
|
||||
private:
|
||||
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
|
||||
typedef EquilReg<RhoCalc> EqReg;
|
||||
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
|
||||
|
||||
PVec pp_;
|
||||
PVec sat_;
|
||||
Vec rs_;
|
||||
Vec rv_;
|
||||
|
||||
template <class RMap>
|
||||
void
|
||||
calcPressSatRsRv(const RMap& reg ,
|
||||
const std::vector< EquilRecord >& rec ,
|
||||
const Opm::BlackoilPropertiesInterface& props,
|
||||
const UnstructuredGrid& G ,
|
||||
const double grav)
|
||||
{
|
||||
typedef Miscibility::NoMixing NoMix;
|
||||
|
||||
for (typename RMap::RegionId
|
||||
r = 0, nr = reg.numRegions();
|
||||
r < nr; ++r)
|
||||
{
|
||||
const typename RMap::CellRange cells = reg.cells(r);
|
||||
|
||||
const int repcell = *cells.begin();
|
||||
const RhoCalc calc(props, repcell);
|
||||
const EqReg eqreg(rec[r], calc,
|
||||
rs_func_[r], rv_func_[r],
|
||||
props.phaseUsage());
|
||||
|
||||
PVec press = phasePressures(G, eqreg, cells, grav);
|
||||
|
||||
const PVec sat = phaseSaturations(eqreg, cells, props, press);
|
||||
|
||||
const int np = props.numPhases();
|
||||
for (int p = 0; p < np; ++p) {
|
||||
copyFromRegion(press[p], cells, pp_[p]);
|
||||
copyFromRegion(sat[p], cells, sat_[p]);
|
||||
}
|
||||
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
|
||||
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
|
||||
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
|
||||
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
|
||||
copyFromRegion(rs, cells, rs_);
|
||||
copyFromRegion(rv, cells, rv_);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class CellRangeType>
|
||||
void copyFromRegion(const Vec& source,
|
||||
const CellRangeType& cells,
|
||||
Vec& destination)
|
||||
{
|
||||
auto s = source.begin();
|
||||
auto c = cells.begin();
|
||||
const auto e = cells.end();
|
||||
for (; c != e; ++c, ++s) {
|
||||
destination[*c] = *s;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
} // namespace DeckDependent
|
||||
} // namespace Equil
|
||||
} // namespace Opm
|
||||
|
||||
#include <opm/core/simulator/initStateEquil_impl.hpp>
|
||||
|
||||
#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
||||
781
opm/core/simulator/initStateEquil_impl.hpp
Normal file
781
opm/core/simulator/initStateEquil_impl.hpp
Normal file
@@ -0,0 +1,781 @@
|
||||
/*
|
||||
Copyright 2014 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
|
||||
#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
#include <functional>
|
||||
#include <vector>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace Details {
|
||||
template <class RHS>
|
||||
class RK4IVP {
|
||||
public:
|
||||
RK4IVP(const RHS& f ,
|
||||
const std::array<double,2>& 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<double>::size_type(N + 1));
|
||||
}
|
||||
|
||||
double
|
||||
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;
|
||||
}
|
||||
|
||||
private:
|
||||
int N_;
|
||||
std::array<double,2> span_;
|
||||
std::vector<double> y_;
|
||||
std::vector<double> f_;
|
||||
|
||||
double
|
||||
stepsize() const { return (span_[1] - span_[0]) / N_; }
|
||||
};
|
||||
|
||||
namespace PhasePressODE {
|
||||
template <class Density>
|
||||
class Water {
|
||||
public:
|
||||
Water(const Density& rho,
|
||||
const int np,
|
||||
const int ix,
|
||||
const double norm_grav)
|
||||
: rho_(rho)
|
||||
, svol_(np, 0)
|
||||
, ix_(ix)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[ix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
operator()(const double /* depth */,
|
||||
const double press) const
|
||||
{
|
||||
return this->density(press) * g_;
|
||||
}
|
||||
|
||||
private:
|
||||
const Density& rho_;
|
||||
std::vector<double> svol_;
|
||||
const int ix_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double press) const
|
||||
{
|
||||
const std::vector<double>& rho = rho_(press, svol_);
|
||||
|
||||
return rho[ix_];
|
||||
}
|
||||
};
|
||||
|
||||
template <class Density, class RS>
|
||||
class Oil {
|
||||
public:
|
||||
Oil(const Density& rho,
|
||||
const RS& rs,
|
||||
const int np,
|
||||
const int oix,
|
||||
const int gix,
|
||||
const double norm_grav)
|
||||
: rho_(rho)
|
||||
, rs_(rs)
|
||||
, svol_(np, 0)
|
||||
, oix_(oix)
|
||||
, gix_(gix)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[oix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
return this->density(depth, press) * g_;
|
||||
}
|
||||
|
||||
private:
|
||||
const Density& rho_;
|
||||
const RS& rs_;
|
||||
mutable std::vector<double> svol_;
|
||||
const int oix_;
|
||||
const int gix_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
if (gix_ >= 0) {
|
||||
svol_[gix_] = rs_(depth, press);
|
||||
}
|
||||
|
||||
const std::vector<double>& rho = rho_(press, svol_);
|
||||
return rho[oix_];
|
||||
}
|
||||
};
|
||||
|
||||
template <class Density, class RV>
|
||||
class Gas {
|
||||
public:
|
||||
Gas(const Density& rho,
|
||||
const RV& rv,
|
||||
const int np,
|
||||
const int gix,
|
||||
const int oix,
|
||||
const double norm_grav)
|
||||
: rho_(rho)
|
||||
, rv_(rv)
|
||||
, svol_(np, 0)
|
||||
, gix_(gix)
|
||||
, oix_(oix)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[gix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
return this->density(depth, press) * g_;
|
||||
}
|
||||
|
||||
private:
|
||||
const Density& rho_;
|
||||
const RV& rv_;
|
||||
mutable std::vector<double> svol_;
|
||||
const int gix_;
|
||||
const int oix_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
if (oix_ >= 0) {
|
||||
svol_[oix_] = rv_(depth, press);
|
||||
}
|
||||
|
||||
const std::vector<double>& rho = rho_(press, svol_);
|
||||
return rho[gix_];
|
||||
}
|
||||
};
|
||||
} // namespace PhasePressODE
|
||||
|
||||
namespace PhaseUsed {
|
||||
inline bool
|
||||
water(const PhaseUsage& pu)
|
||||
{
|
||||
return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
|
||||
}
|
||||
|
||||
inline bool
|
||||
oil(const PhaseUsage& pu)
|
||||
{
|
||||
return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
|
||||
}
|
||||
|
||||
inline bool
|
||||
gas(const PhaseUsage& pu)
|
||||
{
|
||||
return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
|
||||
}
|
||||
} // namespace PhaseUsed
|
||||
|
||||
namespace PhaseIndex {
|
||||
inline int
|
||||
water(const PhaseUsage& pu)
|
||||
{
|
||||
int i = -1;
|
||||
if (PhaseUsed::water(pu)) {
|
||||
i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
|
||||
}
|
||||
|
||||
return i;
|
||||
}
|
||||
|
||||
inline int
|
||||
oil(const PhaseUsage& pu)
|
||||
{
|
||||
int i = -1;
|
||||
if (PhaseUsed::oil(pu)) {
|
||||
i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
|
||||
}
|
||||
|
||||
return i;
|
||||
}
|
||||
|
||||
inline int
|
||||
gas(const PhaseUsage& pu)
|
||||
{
|
||||
int i = -1;
|
||||
if (PhaseUsed::gas(pu)) {
|
||||
i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
|
||||
}
|
||||
|
||||
return i;
|
||||
}
|
||||
} // namespace PhaseIndex
|
||||
|
||||
namespace PhasePressure {
|
||||
template <class PressFunction,
|
||||
class CellRange>
|
||||
void
|
||||
assign(const UnstructuredGrid& G ,
|
||||
const std::array<PressFunction, 2>& f ,
|
||||
const double split,
|
||||
const CellRange& cells,
|
||||
std::vector<double>& p )
|
||||
{
|
||||
const int nd = G.dimensions;
|
||||
|
||||
enum { up = 0, down = 1 };
|
||||
|
||||
std::vector<double>::size_type c = 0;
|
||||
for (typename CellRange::const_iterator
|
||||
ci = cells.begin(), ce = cells.end();
|
||||
ci != ce; ++ci, ++c)
|
||||
{
|
||||
assert (c < p.size());
|
||||
|
||||
const double z = G.cell_centroids[(*ci)*nd + (nd - 1)];
|
||||
p[c] = (z < split) ? f[up](z) : f[down](z);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Region,
|
||||
class CellRange>
|
||||
void
|
||||
water(const UnstructuredGrid& G ,
|
||||
const Region& reg ,
|
||||
const std::array<double,2>& span ,
|
||||
const double grav ,
|
||||
const double po_woc,
|
||||
const CellRange& cells ,
|
||||
std::vector<double>& press )
|
||||
{
|
||||
using PhasePressODE::Water;
|
||||
typedef Water<typename Region::CalcDensity> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav);
|
||||
|
||||
const double z0 = reg.zwoc();
|
||||
const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
|
||||
|
||||
std::array<double,2> up = {{ z0, span[0] }};
|
||||
std::array<double,2> down = {{ z0, span[1] }};
|
||||
|
||||
typedef Details::RK4IVP<ODE> WPress;
|
||||
std::array<WPress,2> wpress = {
|
||||
{
|
||||
WPress(drho, up , p0, 100)
|
||||
,
|
||||
WPress(drho, down, p0, 100)
|
||||
}
|
||||
};
|
||||
|
||||
assign(G, wpress, z0, cells, press);
|
||||
}
|
||||
|
||||
template <class Region,
|
||||
class CellRange>
|
||||
void
|
||||
oil(const UnstructuredGrid& G ,
|
||||
const Region& reg ,
|
||||
const std::array<double,2>& span ,
|
||||
const double grav ,
|
||||
const CellRange& cells ,
|
||||
std::vector<double>& press ,
|
||||
double& po_woc,
|
||||
double& po_goc)
|
||||
{
|
||||
using PhasePressODE::Oil;
|
||||
typedef Oil<typename Region::CalcDensity,
|
||||
typename Region::CalcDissolution> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
ODE drho(reg.densityCalculator(),
|
||||
reg.dissolutionCalculator(),
|
||||
pu.num_phases, oix, gix, grav);
|
||||
|
||||
const double z0 = reg.datum();
|
||||
const double p0 = reg.pressure();
|
||||
|
||||
std::array<double,2> up = {{ z0, span[0] }};
|
||||
std::array<double,2> down = {{ z0, span[1] }};
|
||||
|
||||
typedef Details::RK4IVP<ODE> OPress;
|
||||
std::array<OPress,2> opress = {
|
||||
{
|
||||
OPress(drho, up , p0, 100)
|
||||
,
|
||||
OPress(drho, down, p0, 100)
|
||||
}
|
||||
};
|
||||
|
||||
assign(G, opress, z0, cells, press);
|
||||
|
||||
po_woc = -std::numeric_limits<double>::max();
|
||||
po_goc = -std::numeric_limits<double>::max();
|
||||
|
||||
const double woc = reg.zwoc();
|
||||
// Compute Oil pressure at WOC
|
||||
if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
|
||||
else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
|
||||
else { po_woc = p0; } // WOC *at* datum
|
||||
|
||||
const double goc = reg.zgoc();
|
||||
// Compute Oil pressure at GOC
|
||||
if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
|
||||
else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
|
||||
else { po_goc = p0; } // GOC *at* datum
|
||||
}
|
||||
|
||||
template <class Region,
|
||||
class CellRange>
|
||||
void
|
||||
gas(const UnstructuredGrid& G ,
|
||||
const Region& reg ,
|
||||
const std::array<double,2>& span ,
|
||||
const double grav ,
|
||||
const double po_goc,
|
||||
const CellRange& cells ,
|
||||
std::vector<double>& press )
|
||||
{
|
||||
using PhasePressODE::Gas;
|
||||
typedef Gas<typename Region::CalcDensity,
|
||||
typename Region::CalcEvaporation> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
|
||||
ODE drho(reg.densityCalculator(),
|
||||
reg.evaporationCalculator(),
|
||||
pu.num_phases, gix, oix, grav);
|
||||
|
||||
const double z0 = reg.zgoc();
|
||||
const double p0 = po_goc + reg.pcgo_goc(); // Pcog = Pg - Po
|
||||
|
||||
std::array<double,2> up = {{ z0, span[0] }};
|
||||
std::array<double,2> down = {{ z0, span[1] }};
|
||||
|
||||
typedef Details::RK4IVP<ODE> GPress;
|
||||
std::array<GPress,2> gpress = {
|
||||
{
|
||||
GPress(drho, up , p0, 100)
|
||||
,
|
||||
GPress(drho, down, p0, 100)
|
||||
}
|
||||
};
|
||||
|
||||
assign(G, gpress, z0, cells, press);
|
||||
}
|
||||
} // namespace PhasePressure
|
||||
|
||||
template <class Region,
|
||||
class CellRange>
|
||||
void
|
||||
equilibrateOWG(const UnstructuredGrid& G,
|
||||
const Region& reg,
|
||||
const double grav,
|
||||
const std::array<double,2>& span,
|
||||
const CellRange& cells,
|
||||
std::vector< std::vector<double> >& press)
|
||||
{
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
double po_woc = -1, po_goc = -1;
|
||||
if (PhaseUsed::oil(pu)) {
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
PhasePressure::oil(G, reg, span, grav, cells,
|
||||
press[ oix ], po_woc, po_goc);
|
||||
}
|
||||
|
||||
if (PhaseUsed::water(pu)) {
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
PhasePressure::water(G, reg, span, grav, po_woc,
|
||||
cells, press[ wix ]);
|
||||
}
|
||||
|
||||
if (PhaseUsed::gas(pu)) {
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
PhasePressure::gas(G, reg, span, grav, po_goc,
|
||||
cells, press[ gix ]);
|
||||
}
|
||||
}
|
||||
} // namespace Details
|
||||
|
||||
|
||||
namespace Equil {
|
||||
|
||||
|
||||
template <class Region,
|
||||
class CellRange>
|
||||
std::vector< std::vector<double> >
|
||||
phasePressures(const UnstructuredGrid& G,
|
||||
const Region& reg,
|
||||
const CellRange& cells,
|
||||
const double grav)
|
||||
{
|
||||
std::array<double,2> span =
|
||||
{{ std::numeric_limits<double>::max() ,
|
||||
-std::numeric_limits<double>::max() }}; // Symm. about 0.
|
||||
|
||||
int ncell = 0;
|
||||
{
|
||||
// This code is only supported in three space dimensions
|
||||
assert (G.dimensions == 3);
|
||||
|
||||
const int nd = G.dimensions;
|
||||
|
||||
// Define short-name aliases to reduce visual clutter.
|
||||
const double* const nc = & G.node_coordinates[0];
|
||||
|
||||
const int* const cfp = & G.cell_facepos[0];
|
||||
const int* const cf = & G.cell_faces[0];
|
||||
|
||||
const int* const fnp = & G.face_nodepos[0];
|
||||
const int* const fn = & G.face_nodes[0];
|
||||
|
||||
// Define vertical span as
|
||||
//
|
||||
// [minimum(node depth(cells)), maximum(node depth(cells))]
|
||||
//
|
||||
// Note: We use a sledgehammer approach--looping all
|
||||
// the nodes of all the faces of all the 'cells'--to
|
||||
// compute those bounds. This necessarily entails
|
||||
// visiting some nodes (and faces) multiple times.
|
||||
//
|
||||
// Note: The implementation of 'RK4IVP<>' implicitly
|
||||
// imposes the requirement that cell centroids are all
|
||||
// within this vertical span. That requirement is not
|
||||
// checked.
|
||||
for (typename CellRange::const_iterator
|
||||
ci = cells.begin(), ce = cells.end();
|
||||
ci != ce; ++ci, ++ncell)
|
||||
{
|
||||
for (const int
|
||||
*fi = & cf[ cfp[*ci + 0] ],
|
||||
*fe = & cf[ cfp[*ci + 1] ];
|
||||
fi != fe; ++fi)
|
||||
{
|
||||
for (const int
|
||||
*i = & fn[ fnp[*fi + 0] ],
|
||||
*e = & fn[ fnp[*fi + 1] ];
|
||||
i != e; ++i)
|
||||
{
|
||||
const double z = nc[(*i)*nd + (nd - 1)];
|
||||
|
||||
if (z < span[0]) { span[0] = z; }
|
||||
if (z > span[1]) { span[1] = z; }
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const int np = reg.phaseUsage().num_phases;
|
||||
|
||||
typedef std::vector<double> pval;
|
||||
std::vector<pval> press(np, pval(ncell, 0.0));
|
||||
|
||||
const double z0 = reg.datum();
|
||||
const double zwoc = reg.zwoc ();
|
||||
const double zgoc = reg.zgoc ();
|
||||
|
||||
if (! ((zgoc > z0) || (z0 > zwoc))) {
|
||||
// Datum depth in oil zone (zgoc <= z0 <= zwoc)
|
||||
Details::equilibrateOWG(G, reg, grav, span, cells, press);
|
||||
} else {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
|
||||
}
|
||||
|
||||
return press;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class Region, class CellRange>
|
||||
std::vector< std::vector<double> >
|
||||
phaseSaturations(const Region& reg,
|
||||
const CellRange& cells,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
std::vector< std::vector<double> >& phase_pressures)
|
||||
{
|
||||
const double z0 = reg.datum();
|
||||
const double zwoc = reg.zwoc ();
|
||||
const double zgoc = reg.zgoc ();
|
||||
if ((zgoc > z0) || (z0 > zwoc)) {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
|
||||
}
|
||||
if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
|
||||
}
|
||||
|
||||
std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
|
||||
double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
|
||||
double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
|
||||
|
||||
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
|
||||
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
|
||||
const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
|
||||
const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
std::vector<double>::size_type local_index = 0;
|
||||
for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
|
||||
const int cell = *ci;
|
||||
props.satRange(1, &cell, smin, smax);
|
||||
// Find saturations from pressure differences by
|
||||
// inverting capillary pressure functions.
|
||||
double sw = 0.0;
|
||||
if (water) {
|
||||
const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
|
||||
sw = satFromPc(props, waterpos, cell, pcov);
|
||||
phase_saturations[waterpos][local_index] = sw;
|
||||
}
|
||||
double sg = 0.0;
|
||||
if (gas) {
|
||||
// Note that pcog is defined to be (pg - po), not (po - pg).
|
||||
const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
|
||||
const double increasing = true; // pcog(sg) expected to be increasing function
|
||||
sg = satFromPc(props, gaspos, cell, pcog, increasing);
|
||||
phase_saturations[gaspos][local_index] = sg;
|
||||
}
|
||||
bool overlap = false;
|
||||
if (gas && water && (sg + sw > 1.0)) {
|
||||
// Overlapping gas-oil and oil-water transition
|
||||
// zones can lead to unphysical saturations when
|
||||
// treated as above. Must recalculate using gas-water
|
||||
// capillary pressure.
|
||||
const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
|
||||
sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
|
||||
sg = 1.0 - sw;
|
||||
phase_saturations[waterpos][local_index] = sw;
|
||||
phase_saturations[gaspos][local_index] = sg;
|
||||
overlap = true;
|
||||
}
|
||||
phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
|
||||
|
||||
// Adjust phase pressures for max and min saturation ...
|
||||
double pc[BlackoilPhases::MaxNumPhases];
|
||||
double sat[BlackoilPhases::MaxNumPhases];
|
||||
if (sw > smax[waterpos]-1.0e-6) {
|
||||
sat[waterpos] = smax[waterpos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
|
||||
} else if (overlap || sg > smax[gaspos]-1.0e-6) {
|
||||
sat[gaspos] = smax[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
|
||||
}
|
||||
if (sg < smin[gaspos]+1.0e-6) {
|
||||
sat[gaspos] = smin[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
|
||||
}
|
||||
if (sw < smin[waterpos]+1.0e-6) {
|
||||
sat[waterpos] = smin[waterpos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
|
||||
}
|
||||
}
|
||||
return phase_saturations;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial Rs values.
|
||||
*
|
||||
* \tparam CellRangeType Type of cell range that demarcates the
|
||||
* cells pertaining to the current
|
||||
* equilibration region. Must implement
|
||||
* methods begin() and end() to bound the range
|
||||
* as well as provide an inner type,
|
||||
* const_iterator, to traverse the range.
|
||||
*
|
||||
* \param[in] grid Grid.
|
||||
* \param[in] cells Range that spans the cells of the current
|
||||
* equilibration region.
|
||||
* \param[in] oil_pressure Oil pressure for each cell in range.
|
||||
* \param[in] rs_func Rs as function of pressure and depth.
|
||||
* \return Rs values, one for each cell in the 'cells' range.
|
||||
*/
|
||||
template <class CellRangeType>
|
||||
std::vector<double> computeRs(const UnstructuredGrid& grid,
|
||||
const CellRangeType& cells,
|
||||
const std::vector<double> oil_pressure,
|
||||
const Miscibility::RsFunction& rs_func,
|
||||
const std::vector<double> gas_saturation)
|
||||
{
|
||||
assert(grid.dimensions == 3);
|
||||
std::vector<double> rs(cells.size());
|
||||
int count = 0;
|
||||
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
|
||||
const double depth = grid.cell_centroids[3*(*it) + 2];
|
||||
rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]);
|
||||
}
|
||||
return rs;
|
||||
}
|
||||
|
||||
} // namespace Equil
|
||||
|
||||
|
||||
namespace
|
||||
{
|
||||
/// Convert saturations from a vector of individual phase saturation vectors
|
||||
/// to an interleaved format where all values for a given cell come before all
|
||||
/// values for the next cell, all in a single vector.
|
||||
std::vector<double> convertSats(const std::vector< std::vector<double> >& sat)
|
||||
{
|
||||
const int np = sat.size();
|
||||
const int nc = sat[0].size();
|
||||
std::vector<double> s(np * nc);
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
s[np*c + p] = sat[p][c];
|
||||
}
|
||||
}
|
||||
return s;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial state by an equilibration procedure.
|
||||
*
|
||||
* The following state fields are modified:
|
||||
* pressure(),
|
||||
* saturation(),
|
||||
* surfacevol(),
|
||||
* gasoilratio(),
|
||||
* rv().
|
||||
*
|
||||
* \param[in] grid Grid.
|
||||
* \param[in] props Property object, pvt and capillary properties are used.
|
||||
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
|
||||
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
|
||||
*/
|
||||
void initStateEquil(const UnstructuredGrid& grid,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const EclipseGridParser& deck,
|
||||
const double gravity,
|
||||
BlackoilState& state)
|
||||
{
|
||||
typedef Equil::DeckDependent::InitialStateComputer<EclipseGridParser> ISC;
|
||||
ISC isc(props, deck, grid, gravity);
|
||||
const auto pu = props.phaseUsage();
|
||||
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
|
||||
? pu.phase_pos[BlackoilPhases::Liquid]
|
||||
: pu.phase_pos[BlackoilPhases::Aqua];
|
||||
state.pressure() = isc.press()[ref_phase];
|
||||
state.saturation() = convertSats(isc.saturation());
|
||||
state.gasoilratio() = isc.rs();
|
||||
state.rv() = isc.rv();
|
||||
// TODO: state.surfacevol() must be computed from s, rs, rv.
|
||||
}
|
||||
|
||||
|
||||
void initStateEquil(const UnstructuredGrid& grid,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const Opm::DeckConstPtr newParserDeck,
|
||||
const double gravity,
|
||||
BlackoilState& state)
|
||||
{
|
||||
typedef Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> ISC;
|
||||
ISC isc(props, newParserDeck, grid, gravity);
|
||||
const auto pu = props.phaseUsage();
|
||||
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
|
||||
? pu.phase_pos[BlackoilPhases::Liquid]
|
||||
: pu.phase_pos[BlackoilPhases::Aqua];
|
||||
state.pressure() = isc.press()[ref_phase];
|
||||
state.saturation() = convertSats(isc.saturation());
|
||||
state.gasoilratio() = isc.rs();
|
||||
state.rv() = isc.rv();
|
||||
// TODO: state.surfacevol() must be computed from s, rs, rv.
|
||||
}
|
||||
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
|
||||
@@ -65,8 +65,10 @@ struct Wells
|
||||
/**
|
||||
* Component fractions for each well. Array of size
|
||||
* <CODE>number_of_wells * number_of_phases</CODE>.
|
||||
* This is intended to be used for injection wells. For production wells
|
||||
* the component fractions will vary and cannot be specified a priori.
|
||||
* For injection wells, this gives the injected component mix.
|
||||
* For production wells the component fractions of the wellbore
|
||||
* will vary and cannot be specified a priori, the component mix
|
||||
* given here should be considered a default or preferred mix.
|
||||
*/
|
||||
double *comp_frac;
|
||||
|
||||
|
||||
@@ -550,6 +550,30 @@ namespace Opm
|
||||
OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
|
||||
}
|
||||
set_current_control(well_index, cpos, w_);
|
||||
|
||||
// Set well component fraction to match preferred phase for the well.
|
||||
double cf[3] = { 0.0, 0.0, 0.0 };
|
||||
{
|
||||
switch (well->getPreferredPhase()) {
|
||||
case Phase::WATER:
|
||||
if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
|
||||
OPM_THROW(std::runtime_error, "Water phase not used, yet found water-preferring well.");
|
||||
}
|
||||
cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
|
||||
break;
|
||||
case Phase::OIL:
|
||||
if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
|
||||
OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-preferring well.");
|
||||
}
|
||||
cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||
case Phase::GAS:
|
||||
if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
|
||||
OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-preferring well.");
|
||||
}
|
||||
cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||
}
|
||||
std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
|
||||
}
|
||||
}
|
||||
well_index++;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user