mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-16 20:24:48 -06:00
Refactoring to be able to use template class/methods even if they are not initiated
This commit is contained in:
hnil
parent
10a22c1da9
commit
3e1fe57e60
@ -20,374 +20,10 @@
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
#include <ebos/eclgenerictracermodel.hh>
|
||||
|
||||
#include <opm/simulators/linalg/ilufirstelement.hh>
|
||||
#include <opm/simulators/linalg/PropertyTree.hpp>
|
||||
#include <opm/simulators/linalg/FlexibleSolver.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
#include <opm/grid/CpGrid.hpp>
|
||||
#include <opm/grid/polyhedralgrid.hh>
|
||||
#include <opm/models/discretization/ecfv/ecfvstencil.hh>
|
||||
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Phase.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
|
||||
|
||||
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
|
||||
#include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
|
||||
|
||||
#include <dune/istl/operators.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
|
||||
#if HAVE_DUNE_FEM
|
||||
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
|
||||
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
|
||||
#include <ebos/femcpgridcompat.hh>
|
||||
#endif // HAVE_DUNE_FEM
|
||||
|
||||
#if HAVE_DUNE_ALUGRID
|
||||
#include <dune/alugrid/grid.hh>
|
||||
#include <dune/alugrid/3d/gridview.hh>
|
||||
#include "alucartesianindexmapper.hh"
|
||||
#endif // HAVE_DUNE_ALUGRID
|
||||
|
||||
#include <fmt/format.h>
|
||||
#include <iostream>
|
||||
#include <set>
|
||||
#include <stdexcept>
|
||||
#include <functional>
|
||||
#include <array>
|
||||
#include <string>
|
||||
#include "eclgenerictracermodel_impl.hh"
|
||||
|
||||
namespace Opm {
|
||||
|
||||
#if HAVE_MPI
|
||||
template<class M, class V>
|
||||
struct TracerSolverSelector
|
||||
{
|
||||
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
|
||||
using TracerOperator = Dune::OverlappingSchwarzOperator<M, V, V, Comm>;
|
||||
using type = Dune::FlexibleSolver<TracerOperator>;
|
||||
};
|
||||
|
||||
template<class Vector, class Grid, class Matrix>
|
||||
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>>,
|
||||
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
|
||||
createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&)
|
||||
{
|
||||
OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers.");
|
||||
return {nullptr, nullptr};
|
||||
}
|
||||
|
||||
template<class Vector, class Matrix>
|
||||
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>>,
|
||||
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
|
||||
createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm)
|
||||
{
|
||||
using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>;
|
||||
using TracerSolver = Dune::FlexibleSolver<TracerOperator>;
|
||||
const auto& cellComm = grid.cellCommunication();
|
||||
auto op = std::make_unique<TracerOperator>(M, cellComm);
|
||||
auto dummyWeights = [](){ return Vector();};
|
||||
return {std::move(op), std::make_unique<TracerSolver>(*op, cellComm, prm, dummyWeights, 0)};
|
||||
}
|
||||
#endif
|
||||
|
||||
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
EclGenericTracerModel(const GridView& gridView,
|
||||
const EclipseState& eclState,
|
||||
const CartesianIndexMapper& cartMapper,
|
||||
const DofMapper& dofMapper,
|
||||
const std::function<std::array<double,dimWorld>(int)> centroids)
|
||||
: gridView_(gridView)
|
||||
, eclState_(eclState)
|
||||
, cartMapper_(cartMapper)
|
||||
, dofMapper_(dofMapper)
|
||||
, centroids_(centroids)
|
||||
{
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
Scalar EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
tracerConcentration(int tracerIdx, int globalDofIdx) const
|
||||
{
|
||||
if (tracerConcentration_.empty())
|
||||
return 0.0;
|
||||
|
||||
return tracerConcentration_[tracerIdx][globalDofIdx];
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
|
||||
{
|
||||
this->tracerConcentration_[tracerIdx][globalDofIdx] = value;
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
int EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
numTracers() const
|
||||
{
|
||||
return this->eclState_.tracer().size();
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
std::string EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
fname(int tracerIdx) const
|
||||
{
|
||||
return this->eclState_.tracer()[tracerIdx].fname();
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
double EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
currentConcentration_(const Well& eclWell, const std::string& name) const
|
||||
{
|
||||
return eclWell.getTracerProperties().getConcentration(name);
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
const std::string& EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
name(int tracerIdx) const
|
||||
{
|
||||
return this->eclState_.tracer()[tracerIdx].name;
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
doInit(bool rst, size_t numGridDof,
|
||||
size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx)
|
||||
{
|
||||
const auto& tracers = eclState_.tracer();
|
||||
|
||||
if (tracers.size() == 0)
|
||||
return; // tracer treatment is supposed to be disabled
|
||||
|
||||
// retrieve the number of tracers from the deck
|
||||
const size_t numTracers = tracers.size();
|
||||
tracerConcentration_.resize(numTracers);
|
||||
storageOfTimeIndex1_.resize(numTracers);
|
||||
|
||||
// the phase where the tracer is
|
||||
tracerPhaseIdx_.resize(numTracers);
|
||||
for (size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) {
|
||||
const auto& tracer = tracers[tracerIdx];
|
||||
|
||||
if (tracer.phase == Phase::WATER)
|
||||
tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
|
||||
else if (tracer.phase == Phase::OIL)
|
||||
tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
|
||||
else if (tracer.phase == Phase::GAS)
|
||||
tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
|
||||
|
||||
tracerConcentration_[tracerIdx].resize(numGridDof);
|
||||
storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
|
||||
|
||||
if (rst)
|
||||
continue;
|
||||
|
||||
|
||||
//TBLK keyword
|
||||
if (tracer.free_concentration.has_value()){
|
||||
const auto& free_concentration = tracer.free_concentration.value();
|
||||
int tblkDatasize = free_concentration.size();
|
||||
if (tblkDatasize < cartMapper_.cartesianSize()){
|
||||
throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
|
||||
}
|
||||
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
|
||||
int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
|
||||
tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
|
||||
}
|
||||
}
|
||||
//TVDPF keyword
|
||||
else if (tracer.free_tvdp.has_value()) {
|
||||
const auto& free_tvdp = tracer.free_tvdp.value();
|
||||
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
|
||||
tracerConcentration_[tracerIdx][globalDofIdx] =
|
||||
free_tvdp.evaluate("TRACER_CONCENTRATION",
|
||||
centroids_(globalDofIdx)[2]);
|
||||
}
|
||||
} else
|
||||
throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name));
|
||||
}
|
||||
|
||||
// allocate matrix for storing the Jacobian of the tracer residual
|
||||
tracerMatrix_ = std::make_unique<TracerMatrix>(numGridDof, numGridDof, TracerMatrix::random);
|
||||
|
||||
// find the sparsity pattern of the tracer matrix
|
||||
using NeighborSet = std::set<unsigned>;
|
||||
std::vector<NeighborSet> neighbors(numGridDof);
|
||||
|
||||
Stencil stencil(gridView_, dofMapper_);
|
||||
for (const auto& elem : elements(gridView_)) {
|
||||
stencil.update(elem);
|
||||
|
||||
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
|
||||
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
|
||||
|
||||
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
|
||||
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
|
||||
neighbors[myIdx].insert(neighborIdx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// allocate space for the rows of the matrix
|
||||
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
|
||||
tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
|
||||
tracerMatrix_->endrowsizes();
|
||||
|
||||
// fill the rows with indices. each degree of freedom talks to
|
||||
// all of its neighbors. (it also talks to itself since
|
||||
// degrees of freedom are sometimes quite egocentric.)
|
||||
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
|
||||
typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
|
||||
typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
|
||||
for (; nIt != nEndIt; ++nIt)
|
||||
tracerMatrix_->addindex(dofIdx, *nIt);
|
||||
}
|
||||
tracerMatrix_->endindices();
|
||||
|
||||
const int sizeCartGrid = cartMapper_.cartesianSize();
|
||||
cartToGlobal_.resize(sizeCartGrid);
|
||||
for (unsigned i = 0; i < numGridDof; ++i) {
|
||||
int cartIdx = cartMapper_.cartesianIndex(i);
|
||||
cartToGlobal_[cartIdx] = i;
|
||||
}
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
bool EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
|
||||
{
|
||||
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
|
||||
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
|
||||
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
|
||||
#endif
|
||||
x = 0.0;
|
||||
Scalar tolerance = 1e-2;
|
||||
int maxIter = 100;
|
||||
|
||||
int verbosity = 0;
|
||||
PropertyTree prm;
|
||||
prm.put("maxiter", maxIter);
|
||||
prm.put("tol", tolerance);
|
||||
prm.put("verbosity", verbosity);
|
||||
prm.put("solver", std::string("bicgstab"));
|
||||
prm.put("preconditioner.type", std::string("ParOverILU0"));
|
||||
|
||||
#if HAVE_MPI
|
||||
if(gridView_.grid().comm().size() > 1)
|
||||
{
|
||||
auto [tracerOperator, solver] =
|
||||
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
|
||||
(void) tracerOperator;
|
||||
|
||||
Dune::InverseOperatorResult result;
|
||||
solver->apply(x, b, result);
|
||||
|
||||
// return the result of the solver
|
||||
return result.converged;
|
||||
}
|
||||
else
|
||||
{
|
||||
#endif
|
||||
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
|
||||
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
|
||||
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
|
||||
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
|
||||
|
||||
TracerOperator tracerOperator(M);
|
||||
TracerScalarProduct tracerScalarProduct;
|
||||
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
|
||||
|
||||
TracerSolver solver (tracerOperator, tracerScalarProduct,
|
||||
tracerPreconditioner, tolerance, maxIter,
|
||||
verbosity);
|
||||
|
||||
Dune::InverseOperatorResult result;
|
||||
solver.apply(x, b, result);
|
||||
|
||||
// return the result of the solver
|
||||
return result.converged;
|
||||
#if HAVE_MPI
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
bool EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
|
||||
{
|
||||
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
|
||||
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
|
||||
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
|
||||
#endif
|
||||
Scalar tolerance = 1e-2;
|
||||
int maxIter = 100;
|
||||
|
||||
int verbosity = 0;
|
||||
PropertyTree prm;
|
||||
prm.put("maxiter", maxIter);
|
||||
prm.put("tol", tolerance);
|
||||
prm.put("verbosity", verbosity);
|
||||
prm.put("solver", std::string("bicgstab"));
|
||||
prm.put("preconditioner.type", std::string("ParOverILU0"));
|
||||
|
||||
#if HAVE_MPI
|
||||
if(gridView_.grid().comm().size() > 1)
|
||||
{
|
||||
auto [tracerOperator, solver] =
|
||||
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
|
||||
(void) tracerOperator;
|
||||
bool converged = true;
|
||||
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
|
||||
x[nrhs] = 0.0;
|
||||
Dune::InverseOperatorResult result;
|
||||
solver->apply(x[nrhs], b[nrhs], result);
|
||||
converged = (converged && result.converged);
|
||||
}
|
||||
return converged;
|
||||
}
|
||||
else
|
||||
{
|
||||
#endif
|
||||
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
|
||||
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
|
||||
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
|
||||
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
|
||||
|
||||
TracerOperator tracerOperator(M);
|
||||
TracerScalarProduct tracerScalarProduct;
|
||||
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
|
||||
|
||||
TracerSolver solver (tracerOperator, tracerScalarProduct,
|
||||
tracerPreconditioner, tolerance, maxIter,
|
||||
verbosity);
|
||||
|
||||
bool converged = true;
|
||||
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
|
||||
x[nrhs] = 0.0;
|
||||
Dune::InverseOperatorResult result;
|
||||
solver.apply(x[nrhs], b[nrhs], result);
|
||||
converged = (converged && result.converged);
|
||||
}
|
||||
|
||||
// return the result of the solver
|
||||
return converged;
|
||||
#if HAVE_MPI
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
#if HAVE_DUNE_FEM
|
||||
template class EclGenericTracerModel<Dune::CpGrid,
|
||||
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
|
||||
|
||||
399
ebos/eclgenerictracermodel_impl.hh
Normal file
399
ebos/eclgenerictracermodel_impl.hh
Normal file
@ -0,0 +1,399 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
/**
|
||||
* \file
|
||||
*
|
||||
* \copydoc Opm::EclTracerModel
|
||||
*/
|
||||
#ifndef EWOMS_ECL_GENERIC_TRACER_MODEL_IMPL_HH
|
||||
#define EWOMS_ECL_GENERIC_TRACER_MODEL_IMPL_HH
|
||||
|
||||
#include <ebos/eclgenerictracermodel.hh>
|
||||
|
||||
#include <opm/simulators/linalg/ilufirstelement.hh>
|
||||
#include <opm/simulators/linalg/PropertyTree.hpp>
|
||||
#include <opm/simulators/linalg/FlexibleSolver.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
#include <opm/grid/CpGrid.hpp>
|
||||
#include <opm/grid/polyhedralgrid.hh>
|
||||
#include <opm/models/discretization/ecfv/ecfvstencil.hh>
|
||||
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Phase.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Tables/TracerVdTable.hpp>
|
||||
|
||||
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
|
||||
#include <opm/input/eclipse/Schedule/Well/WellTracerProperties.hpp>
|
||||
|
||||
#include <dune/istl/operators.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
|
||||
#if HAVE_DUNE_FEM
|
||||
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
|
||||
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
|
||||
#include <ebos/femcpgridcompat.hh>
|
||||
#endif // HAVE_DUNE_FEM
|
||||
|
||||
#if HAVE_DUNE_ALUGRID
|
||||
#include <dune/alugrid/grid.hh>
|
||||
#include <dune/alugrid/3d/gridview.hh>
|
||||
#include "alucartesianindexmapper.hh"
|
||||
#endif // HAVE_DUNE_ALUGRID
|
||||
|
||||
#include <fmt/format.h>
|
||||
#include <iostream>
|
||||
#include <set>
|
||||
#include <stdexcept>
|
||||
#include <functional>
|
||||
#include <array>
|
||||
#include <string>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
#if HAVE_MPI
|
||||
template<class M, class V>
|
||||
struct TracerSolverSelector
|
||||
{
|
||||
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
|
||||
using TracerOperator = Dune::OverlappingSchwarzOperator<M, V, V, Comm>;
|
||||
using type = Dune::FlexibleSolver<TracerOperator>;
|
||||
};
|
||||
|
||||
template<class Vector, class Grid, class Matrix>
|
||||
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>>,
|
||||
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
|
||||
createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&)
|
||||
{
|
||||
OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers.");
|
||||
return {nullptr, nullptr};
|
||||
}
|
||||
|
||||
template<class Vector, class Matrix>
|
||||
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>>,
|
||||
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
|
||||
createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm)
|
||||
{
|
||||
using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
|
||||
Dune::OwnerOverlapCopyCommunication<int,int>>;
|
||||
using TracerSolver = Dune::FlexibleSolver<TracerOperator>;
|
||||
const auto& cellComm = grid.cellCommunication();
|
||||
auto op = std::make_unique<TracerOperator>(M, cellComm);
|
||||
auto dummyWeights = [](){ return Vector();};
|
||||
return {std::move(op), std::make_unique<TracerSolver>(*op, cellComm, prm, dummyWeights, 0)};
|
||||
}
|
||||
#endif
|
||||
|
||||
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
EclGenericTracerModel(const GridView& gridView,
|
||||
const EclipseState& eclState,
|
||||
const CartesianIndexMapper& cartMapper,
|
||||
const DofMapper& dofMapper,
|
||||
const std::function<std::array<double,dimWorld>(int)> centroids)
|
||||
: gridView_(gridView)
|
||||
, eclState_(eclState)
|
||||
, cartMapper_(cartMapper)
|
||||
, dofMapper_(dofMapper)
|
||||
, centroids_(centroids)
|
||||
{
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
Scalar EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
tracerConcentration(int tracerIdx, int globalDofIdx) const
|
||||
{
|
||||
if (tracerConcentration_.empty())
|
||||
return 0.0;
|
||||
|
||||
return tracerConcentration_[tracerIdx][globalDofIdx];
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
|
||||
{
|
||||
this->tracerConcentration_[tracerIdx][globalDofIdx] = value;
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
int EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
numTracers() const
|
||||
{
|
||||
return this->eclState_.tracer().size();
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
std::string EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
fname(int tracerIdx) const
|
||||
{
|
||||
return this->eclState_.tracer()[tracerIdx].fname();
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
double EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
currentConcentration_(const Well& eclWell, const std::string& name) const
|
||||
{
|
||||
return eclWell.getTracerProperties().getConcentration(name);
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
const std::string& EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
name(int tracerIdx) const
|
||||
{
|
||||
return this->eclState_.tracer()[tracerIdx].name;
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
void EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
doInit(bool rst, size_t numGridDof,
|
||||
size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx)
|
||||
{
|
||||
const auto& tracers = eclState_.tracer();
|
||||
|
||||
if (tracers.size() == 0)
|
||||
return; // tracer treatment is supposed to be disabled
|
||||
|
||||
// retrieve the number of tracers from the deck
|
||||
const size_t numTracers = tracers.size();
|
||||
tracerConcentration_.resize(numTracers);
|
||||
storageOfTimeIndex1_.resize(numTracers);
|
||||
|
||||
// the phase where the tracer is
|
||||
tracerPhaseIdx_.resize(numTracers);
|
||||
for (size_t tracerIdx = 0; tracerIdx < numTracers; tracerIdx++) {
|
||||
const auto& tracer = tracers[tracerIdx];
|
||||
|
||||
if (tracer.phase == Phase::WATER)
|
||||
tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
|
||||
else if (tracer.phase == Phase::OIL)
|
||||
tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
|
||||
else if (tracer.phase == Phase::GAS)
|
||||
tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
|
||||
|
||||
tracerConcentration_[tracerIdx].resize(numGridDof);
|
||||
storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
|
||||
|
||||
if (rst)
|
||||
continue;
|
||||
|
||||
|
||||
//TBLK keyword
|
||||
if (tracer.free_concentration.has_value()){
|
||||
const auto& free_concentration = tracer.free_concentration.value();
|
||||
int tblkDatasize = free_concentration.size();
|
||||
if (tblkDatasize < cartMapper_.cartesianSize()){
|
||||
throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
|
||||
}
|
||||
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
|
||||
int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
|
||||
tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
|
||||
}
|
||||
}
|
||||
//TVDPF keyword
|
||||
else if (tracer.free_tvdp.has_value()) {
|
||||
const auto& free_tvdp = tracer.free_tvdp.value();
|
||||
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
|
||||
tracerConcentration_[tracerIdx][globalDofIdx] =
|
||||
free_tvdp.evaluate("TRACER_CONCENTRATION",
|
||||
centroids_(globalDofIdx)[2]);
|
||||
}
|
||||
} else
|
||||
throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name));
|
||||
}
|
||||
|
||||
// allocate matrix for storing the Jacobian of the tracer residual
|
||||
tracerMatrix_ = std::make_unique<TracerMatrix>(numGridDof, numGridDof, TracerMatrix::random);
|
||||
|
||||
// find the sparsity pattern of the tracer matrix
|
||||
using NeighborSet = std::set<unsigned>;
|
||||
std::vector<NeighborSet> neighbors(numGridDof);
|
||||
|
||||
Stencil stencil(gridView_, dofMapper_);
|
||||
for (const auto& elem : elements(gridView_)) {
|
||||
stencil.update(elem);
|
||||
|
||||
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
|
||||
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
|
||||
|
||||
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
|
||||
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
|
||||
neighbors[myIdx].insert(neighborIdx);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// allocate space for the rows of the matrix
|
||||
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
|
||||
tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
|
||||
tracerMatrix_->endrowsizes();
|
||||
|
||||
// fill the rows with indices. each degree of freedom talks to
|
||||
// all of its neighbors. (it also talks to itself since
|
||||
// degrees of freedom are sometimes quite egocentric.)
|
||||
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
|
||||
typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
|
||||
typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
|
||||
for (; nIt != nEndIt; ++nIt)
|
||||
tracerMatrix_->addindex(dofIdx, *nIt);
|
||||
}
|
||||
tracerMatrix_->endindices();
|
||||
|
||||
const int sizeCartGrid = cartMapper_.cartesianSize();
|
||||
cartToGlobal_.resize(sizeCartGrid);
|
||||
for (unsigned i = 0; i < numGridDof; ++i) {
|
||||
int cartIdx = cartMapper_.cartesianIndex(i);
|
||||
cartToGlobal_[cartIdx] = i;
|
||||
}
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
bool EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
|
||||
{
|
||||
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
|
||||
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
|
||||
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
|
||||
#endif
|
||||
x = 0.0;
|
||||
Scalar tolerance = 1e-2;
|
||||
int maxIter = 100;
|
||||
|
||||
int verbosity = 0;
|
||||
PropertyTree prm;
|
||||
prm.put("maxiter", maxIter);
|
||||
prm.put("tol", tolerance);
|
||||
prm.put("verbosity", verbosity);
|
||||
prm.put("solver", std::string("bicgstab"));
|
||||
prm.put("preconditioner.type", std::string("ParOverILU0"));
|
||||
|
||||
#if HAVE_MPI
|
||||
if(gridView_.grid().comm().size() > 1)
|
||||
{
|
||||
auto [tracerOperator, solver] =
|
||||
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
|
||||
(void) tracerOperator;
|
||||
|
||||
Dune::InverseOperatorResult result;
|
||||
solver->apply(x, b, result);
|
||||
|
||||
// return the result of the solver
|
||||
return result.converged;
|
||||
}
|
||||
else
|
||||
{
|
||||
#endif
|
||||
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
|
||||
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
|
||||
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
|
||||
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
|
||||
|
||||
TracerOperator tracerOperator(M);
|
||||
TracerScalarProduct tracerScalarProduct;
|
||||
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
|
||||
|
||||
TracerSolver solver (tracerOperator, tracerScalarProduct,
|
||||
tracerPreconditioner, tolerance, maxIter,
|
||||
verbosity);
|
||||
|
||||
Dune::InverseOperatorResult result;
|
||||
solver.apply(x, b, result);
|
||||
|
||||
// return the result of the solver
|
||||
return result.converged;
|
||||
#if HAVE_MPI
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
|
||||
bool EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
|
||||
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
|
||||
{
|
||||
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
|
||||
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
|
||||
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
|
||||
#endif
|
||||
Scalar tolerance = 1e-2;
|
||||
int maxIter = 100;
|
||||
|
||||
int verbosity = 0;
|
||||
PropertyTree prm;
|
||||
prm.put("maxiter", maxIter);
|
||||
prm.put("tol", tolerance);
|
||||
prm.put("verbosity", verbosity);
|
||||
prm.put("solver", std::string("bicgstab"));
|
||||
prm.put("preconditioner.type", std::string("ParOverILU0"));
|
||||
|
||||
#if HAVE_MPI
|
||||
if(gridView_.grid().comm().size() > 1)
|
||||
{
|
||||
auto [tracerOperator, solver] =
|
||||
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
|
||||
(void) tracerOperator;
|
||||
bool converged = true;
|
||||
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
|
||||
x[nrhs] = 0.0;
|
||||
Dune::InverseOperatorResult result;
|
||||
solver->apply(x[nrhs], b[nrhs], result);
|
||||
converged = (converged && result.converged);
|
||||
}
|
||||
return converged;
|
||||
}
|
||||
else
|
||||
{
|
||||
#endif
|
||||
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
|
||||
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
|
||||
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
|
||||
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
|
||||
|
||||
TracerOperator tracerOperator(M);
|
||||
TracerScalarProduct tracerScalarProduct;
|
||||
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
|
||||
|
||||
TracerSolver solver (tracerOperator, tracerScalarProduct,
|
||||
tracerPreconditioner, tolerance, maxIter,
|
||||
verbosity);
|
||||
|
||||
bool converged = true;
|
||||
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
|
||||
x[nrhs] = 0.0;
|
||||
Dune::InverseOperatorResult result;
|
||||
solver.apply(x[nrhs], b[nrhs], result);
|
||||
converged = (converged && result.converged);
|
||||
}
|
||||
|
||||
// return the result of the solver
|
||||
return converged;
|
||||
#if HAVE_MPI
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
#endif
|
||||
@ -22,597 +22,11 @@
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
|
||||
#include <ebos/equil/equilibrationhelpers.hh>
|
||||
|
||||
#include <opm/common/TimingMacros.hpp>
|
||||
|
||||
#include <opm/common/utility/numeric/RootFinders.hpp>
|
||||
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
|
||||
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
|
||||
|
||||
#include <fmt/format.h>
|
||||
#include "equilibrationhelpers_impl.hh"
|
||||
|
||||
namespace Opm {
|
||||
namespace EQUIL {
|
||||
|
||||
using FluidSystemSimple = BlackOilFluidSystem<double>;
|
||||
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
using SatOnlyFluidState = SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystemSimple,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false>;
|
||||
|
||||
namespace Miscibility {
|
||||
|
||||
template<class FluidSystem>
|
||||
RsVD<FluidSystem>::RsVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rsVsDepth_(depth, rs)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
if (satGas > 0.0) {
|
||||
return satRs(press, temp);
|
||||
}
|
||||
else {
|
||||
if (rsVsDepth_.xMin() > depth)
|
||||
return rsVsDepth_.valueAt(0);
|
||||
else if (rsVsDepth_.xMax() < depth)
|
||||
return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1);
|
||||
return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsVD<FluidSystem>::satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
PBVD<FluidSystem>::PBVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& pbub)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, pbubVsDepth_(depth, pbub)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PBVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double cellPress,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
double press = cellPress;
|
||||
if (satGas <= 0.0) {
|
||||
if (pbubVsDepth_.xMin() > depth)
|
||||
press = pbubVsDepth_.valueAt(0);
|
||||
else if (pbubVsDepth_.xMax() < depth)
|
||||
press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
|
||||
else
|
||||
press = pbubVsDepth_.eval(depth, /*extrapolate=*/false);
|
||||
}
|
||||
return satRs(std::min(press, cellPress), temp);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PBVD<FluidSystem>::
|
||||
satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
PDVD<FluidSystem>::PDVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& pdew)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, pdewVsDepth_(depth, pdew)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PDVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double cellPress,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
double press = cellPress;
|
||||
if (satOil <= 0.0) {
|
||||
if (pdewVsDepth_.xMin() > depth)
|
||||
press = pdewVsDepth_.valueAt(0);
|
||||
else if (pdewVsDepth_.xMax() < depth)
|
||||
press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
|
||||
else
|
||||
press = pdewVsDepth_.eval(depth, /*extrapolate=*/false);
|
||||
}
|
||||
return satRv(std::min(press, cellPress), temp);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PDVD<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RvVD<FluidSystem>::RvVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rv)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rvVsDepth_(depth, rv)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
if (std::abs(satOil) > 1e-16) {
|
||||
return satRv(press, temp);
|
||||
}
|
||||
else {
|
||||
if (rvVsDepth_.xMin() > depth)
|
||||
return rvVsDepth_.valueAt(0);
|
||||
else if (rvVsDepth_.xMax() < depth)
|
||||
return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1);
|
||||
return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvVD<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RvwVD<FluidSystem>::RvwVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rvw)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rvwVsDepth_(depth, rvw)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satWat) const
|
||||
{
|
||||
if (std::abs(satWat) > 1e-16) {
|
||||
return satRvw(press, temp); //saturated Rvw
|
||||
}
|
||||
else {
|
||||
if (rvwVsDepth_.xMin() > depth)
|
||||
return rvwVsDepth_.valueAt(0);
|
||||
else if (rvwVsDepth_.xMax() < depth)
|
||||
return rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1);
|
||||
return std::min(satRvw(press, temp), rvwVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwVD<FluidSystem>::
|
||||
satRvw(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RsSatAtContact<FluidSystem>::
|
||||
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rsSatContact_ = satRs(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsSatAtContact<FluidSystem>::
|
||||
operator()(const double /* depth */,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
if (satGas > 0.0) {
|
||||
return satRs(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRs(press, temp), rsSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsSatAtContact<FluidSystem>::
|
||||
satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
RvSatAtContact<FluidSystem>::
|
||||
RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rvSatContact_ = satRv(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvSatAtContact<FluidSystem>::
|
||||
operator()(const double /*depth*/,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
if (satOil > 0.0) {
|
||||
return satRv(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRv(press, temp), rvSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvSatAtContact<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
RvwSatAtContact<FluidSystem>::
|
||||
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rvwSatContact_ = satRvw(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwSatAtContact<FluidSystem>::
|
||||
operator()(const double /*depth*/,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satWat) const
|
||||
{
|
||||
if (satWat > 0.0) {
|
||||
return satRvw(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRvw(press, temp), rvwSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwSatAtContact<FluidSystem>::
|
||||
satRvw(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
|
||||
}
|
||||
|
||||
} // namespace Miscibility
|
||||
|
||||
EquilReg::EquilReg(const EquilRecord& rec,
|
||||
std::shared_ptr<Miscibility::RsFunction> rs,
|
||||
std::shared_ptr<Miscibility::RsFunction> rv,
|
||||
std::shared_ptr<Miscibility::RsFunction> rvw,
|
||||
const TabulatedFunction& saltVdTable,
|
||||
const int pvtIdx)
|
||||
: rec_ (rec)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, rvw_ (rvw)
|
||||
, saltVdTable_ (saltVdTable)
|
||||
, pvtIdx_ (pvtIdx)
|
||||
{
|
||||
}
|
||||
|
||||
double EquilReg::datum() const
|
||||
{
|
||||
return this->rec_.datumDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pressure() const
|
||||
{
|
||||
return this->rec_.datumDepthPressure();
|
||||
}
|
||||
|
||||
double EquilReg::zwoc() const
|
||||
{
|
||||
return this->rec_.waterOilContactDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pcowWoc() const
|
||||
{
|
||||
return this->rec_.waterOilContactCapillaryPressure();
|
||||
}
|
||||
|
||||
double EquilReg::zgoc() const
|
||||
{
|
||||
return this->rec_.gasOilContactDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pcgoGoc() const
|
||||
{
|
||||
return this->rec_.gasOilContactCapillaryPressure();
|
||||
}
|
||||
|
||||
int EquilReg::equilibrationAccuracy() const
|
||||
{
|
||||
return this->rec_.initializationTargetAccuracy();
|
||||
}
|
||||
|
||||
const EquilReg::CalcDissolution&
|
||||
EquilReg::dissolutionCalculator() const
|
||||
{
|
||||
return *this->rs_;
|
||||
}
|
||||
|
||||
const EquilReg::CalcEvaporation&
|
||||
EquilReg::evaporationCalculator() const
|
||||
{
|
||||
return *this->rv_;
|
||||
}
|
||||
|
||||
const EquilReg::CalcWaterEvaporation&
|
||||
EquilReg::waterEvaporationCalculator() const
|
||||
{
|
||||
return *this->rvw_;
|
||||
}
|
||||
|
||||
const EquilReg::TabulatedFunction&
|
||||
EquilReg::saltVdTable() const
|
||||
{
|
||||
return saltVdTable_;
|
||||
}
|
||||
|
||||
int EquilReg::pvtIdx() const
|
||||
{
|
||||
return this->pvtIdx_;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
PcEq<FluidSystem,MaterialLawManager>::
|
||||
PcEq(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
targetPc_(targetPc)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double PcEq<FluidSystem,MaterialLawManager>::
|
||||
operator()(double s) const
|
||||
{
|
||||
const auto& matParams = materialLawManager_.materialLawParams(cell_);
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(phase_, s);
|
||||
|
||||
std::array<double, FluidSystem::numPhases> pc{0.0};
|
||||
using MaterialLaw = typename MaterialLawManager::MaterialLaw;
|
||||
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
|
||||
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
|
||||
return pcPhase - targetPc_;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
PcEqSum<FluidSystem,MaterialLawManager>::
|
||||
PcEqSum(const MaterialLawManager& materialLawManager,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
targetPc_(targetPc)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double PcEqSum<FluidSystem,MaterialLawManager>::
|
||||
operator()(double s) const
|
||||
{
|
||||
const auto& matParams = materialLawManager_.materialLawParams(cell_);
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(phase1_, s);
|
||||
fluidState.setSaturation(phase2_, 1.0 - s);
|
||||
|
||||
std::array<double, FluidSystem::numPhases> pc {0.0};
|
||||
|
||||
using MaterialLaw = typename MaterialLawManager::MaterialLaw;
|
||||
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
|
||||
double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
|
||||
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
|
||||
return pc1 + pc2 - targetPc_;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double minSaturations(const MaterialLawManager& materialLawManager,
|
||||
const int phase, const int cell)
|
||||
{
|
||||
const auto& scaledDrainageInfo =
|
||||
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
|
||||
|
||||
// Find minimum and maximum saturations.
|
||||
switch(phase) {
|
||||
case FluidSystem::waterPhaseIdx:
|
||||
return scaledDrainageInfo.Swl;
|
||||
|
||||
case FluidSystem::gasPhaseIdx:
|
||||
return scaledDrainageInfo.Sgl;
|
||||
|
||||
case FluidSystem::oilPhaseIdx:
|
||||
throw std::runtime_error("Min saturation not implemented for oil phase.");
|
||||
|
||||
default:
|
||||
throw std::runtime_error("Unknown phaseIdx .");
|
||||
}
|
||||
return -1.0;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double maxSaturations(const MaterialLawManager& materialLawManager,
|
||||
const int phase, const int cell)
|
||||
{
|
||||
const auto& scaledDrainageInfo =
|
||||
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
|
||||
|
||||
// Find minimum and maximum saturations.
|
||||
switch(phase) {
|
||||
case FluidSystem::waterPhaseIdx:
|
||||
return scaledDrainageInfo.Swu;
|
||||
|
||||
case FluidSystem::gasPhaseIdx:
|
||||
return scaledDrainageInfo.Sgu;
|
||||
|
||||
case FluidSystem::oilPhaseIdx:
|
||||
throw std::runtime_error("Max saturation not implemented for oil phase.");
|
||||
|
||||
default:
|
||||
throw std::runtime_error("Unknown phaseIdx .");
|
||||
}
|
||||
return -1.0;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double satFromPc(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double targetPc,
|
||||
const bool increasing)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
// Create the equation f(s) = pc(s) - targetPc
|
||||
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, targetPc);
|
||||
double f0 = f(s0);
|
||||
double f1 = f(s1);
|
||||
if (!std::isfinite(f0 + f1))
|
||||
throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1));
|
||||
|
||||
if (f0 <= 0.0)
|
||||
return s0;
|
||||
else if (f1 >= 0.0)
|
||||
return s1;
|
||||
|
||||
const double tol = 1e-10;
|
||||
// should at least converge in 2 times bisection but some safety here:
|
||||
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
|
||||
int usedIterations = -1;
|
||||
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
|
||||
return root;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
double s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc
|
||||
const PcEqSum<FluidSystem, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, targetPc);
|
||||
double f0 = f(s0);
|
||||
double f1 = f(s1);
|
||||
if (f0 <= 0.0)
|
||||
return s0;
|
||||
else if (f1 >= 0.0)
|
||||
return s1;
|
||||
|
||||
assert(f0 > 0.0 && f1 < 0.0);
|
||||
const double tol = 1e-10;
|
||||
// should at least converge in 2 times bisection but some safety here:
|
||||
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
|
||||
int usedIterations = -1;
|
||||
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
|
||||
return root;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double satFromDepth(const MaterialLawManager& materialLawManager,
|
||||
const double cellDepth,
|
||||
const double contactDepth,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const bool increasing)
|
||||
{
|
||||
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
if (cellDepth < contactDepth) {
|
||||
return s0;
|
||||
}
|
||||
else {
|
||||
return s1;
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
bool isConstPc(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell)
|
||||
{
|
||||
// Create the equation f(s) = pc(s);
|
||||
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, 0);
|
||||
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
|
||||
}
|
||||
|
||||
|
||||
using MatLaw = EclMaterialLawManager<ThreePhaseMaterialTraits<double,0,1,2>>;
|
||||
using FS = BlackOilFluidSystem<double>;
|
||||
|
||||
|
||||
615
ebos/equil/equilibrationhelpers_impl.hh
Normal file
615
ebos/equil/equilibrationhelpers_impl.hh
Normal file
@ -0,0 +1,615 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Consult the COPYING file in the top-level source directory of this
|
||||
module for the precise wording of the license and the list of
|
||||
copyright holders.
|
||||
*/
|
||||
|
||||
#include <ebos/equil/equilibrationhelpers.hh>
|
||||
|
||||
#include <opm/common/TimingMacros.hpp>
|
||||
|
||||
#include <opm/common/utility/numeric/RootFinders.hpp>
|
||||
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
|
||||
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
|
||||
|
||||
#include <fmt/format.h>
|
||||
|
||||
namespace Opm {
|
||||
namespace EQUIL {
|
||||
|
||||
using FluidSystemSimple = BlackOilFluidSystem<double>;
|
||||
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
using SatOnlyFluidState = SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystemSimple,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false>;
|
||||
|
||||
namespace Miscibility {
|
||||
|
||||
template<class FluidSystem>
|
||||
RsVD<FluidSystem>::RsVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rsVsDepth_(depth, rs)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
if (satGas > 0.0) {
|
||||
return satRs(press, temp);
|
||||
}
|
||||
else {
|
||||
if (rsVsDepth_.xMin() > depth)
|
||||
return rsVsDepth_.valueAt(0);
|
||||
else if (rsVsDepth_.xMax() < depth)
|
||||
return rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1);
|
||||
return std::min(satRs(press, temp), rsVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsVD<FluidSystem>::satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
PBVD<FluidSystem>::PBVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& pbub)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, pbubVsDepth_(depth, pbub)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PBVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double cellPress,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
double press = cellPress;
|
||||
if (satGas <= 0.0) {
|
||||
if (pbubVsDepth_.xMin() > depth)
|
||||
press = pbubVsDepth_.valueAt(0);
|
||||
else if (pbubVsDepth_.xMax() < depth)
|
||||
press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
|
||||
else
|
||||
press = pbubVsDepth_.eval(depth, /*extrapolate=*/false);
|
||||
}
|
||||
return satRs(std::min(press, cellPress), temp);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PBVD<FluidSystem>::
|
||||
satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
PDVD<FluidSystem>::PDVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& pdew)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, pdewVsDepth_(depth, pdew)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PDVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double cellPress,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
double press = cellPress;
|
||||
if (satOil <= 0.0) {
|
||||
if (pdewVsDepth_.xMin() > depth)
|
||||
press = pdewVsDepth_.valueAt(0);
|
||||
else if (pdewVsDepth_.xMax() < depth)
|
||||
press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
|
||||
else
|
||||
press = pdewVsDepth_.eval(depth, /*extrapolate=*/false);
|
||||
}
|
||||
return satRv(std::min(press, cellPress), temp);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double PDVD<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RvVD<FluidSystem>::RvVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rv)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rvVsDepth_(depth, rv)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
if (std::abs(satOil) > 1e-16) {
|
||||
return satRv(press, temp);
|
||||
}
|
||||
else {
|
||||
if (rvVsDepth_.xMin() > depth)
|
||||
return rvVsDepth_.valueAt(0);
|
||||
else if (rvVsDepth_.xMax() < depth)
|
||||
return rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1);
|
||||
return std::min(satRv(press, temp), rvVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvVD<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RvwVD<FluidSystem>::RvwVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rvw)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, rvwVsDepth_(depth, rvw)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwVD<FluidSystem>::
|
||||
operator()(const double depth,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satWat) const
|
||||
{
|
||||
if (std::abs(satWat) > 1e-16) {
|
||||
return satRvw(press, temp); //saturated Rvw
|
||||
}
|
||||
else {
|
||||
if (rvwVsDepth_.xMin() > depth)
|
||||
return rvwVsDepth_.valueAt(0);
|
||||
else if (rvwVsDepth_.xMax() < depth)
|
||||
return rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1);
|
||||
return std::min(satRvw(press, temp), rvwVsDepth_.eval(depth, /*extrapolate=*/false));
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwVD<FluidSystem>::
|
||||
satRvw(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
|
||||
template<class FluidSystem>
|
||||
RsSatAtContact<FluidSystem>::
|
||||
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rsSatContact_ = satRs(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsSatAtContact<FluidSystem>::
|
||||
operator()(const double /* depth */,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satGas) const
|
||||
{
|
||||
if (satGas > 0.0) {
|
||||
return satRs(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRs(press, temp), rsSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RsSatAtContact<FluidSystem>::
|
||||
satRs(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
RvSatAtContact<FluidSystem>::
|
||||
RvSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rvSatContact_ = satRv(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvSatAtContact<FluidSystem>::
|
||||
operator()(const double /*depth*/,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satOil) const
|
||||
{
|
||||
if (satOil > 0.0) {
|
||||
return satRv(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRv(press, temp), rvSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvSatAtContact<FluidSystem>::
|
||||
satRv(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
RvwSatAtContact<FluidSystem>::
|
||||
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
rvwSatContact_ = satRvw(pContact, T_contact);
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwSatAtContact<FluidSystem>::
|
||||
operator()(const double /*depth*/,
|
||||
const double press,
|
||||
const double temp,
|
||||
const double satWat) const
|
||||
{
|
||||
if (satWat > 0.0) {
|
||||
return satRvw(press, temp);
|
||||
}
|
||||
else {
|
||||
return std::min(satRvw(press, temp), rvwSatContact_);
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem>
|
||||
double RvwSatAtContact<FluidSystem>::
|
||||
satRvw(const double press, const double temp) const
|
||||
{
|
||||
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
|
||||
}
|
||||
|
||||
} // namespace Miscibility
|
||||
|
||||
EquilReg::EquilReg(const EquilRecord& rec,
|
||||
std::shared_ptr<Miscibility::RsFunction> rs,
|
||||
std::shared_ptr<Miscibility::RsFunction> rv,
|
||||
std::shared_ptr<Miscibility::RsFunction> rvw,
|
||||
const TabulatedFunction& saltVdTable,
|
||||
const int pvtIdx)
|
||||
: rec_ (rec)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, rvw_ (rvw)
|
||||
, saltVdTable_ (saltVdTable)
|
||||
, pvtIdx_ (pvtIdx)
|
||||
{
|
||||
}
|
||||
|
||||
double EquilReg::datum() const
|
||||
{
|
||||
return this->rec_.datumDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pressure() const
|
||||
{
|
||||
return this->rec_.datumDepthPressure();
|
||||
}
|
||||
|
||||
double EquilReg::zwoc() const
|
||||
{
|
||||
return this->rec_.waterOilContactDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pcowWoc() const
|
||||
{
|
||||
return this->rec_.waterOilContactCapillaryPressure();
|
||||
}
|
||||
|
||||
double EquilReg::zgoc() const
|
||||
{
|
||||
return this->rec_.gasOilContactDepth();
|
||||
}
|
||||
|
||||
double EquilReg::pcgoGoc() const
|
||||
{
|
||||
return this->rec_.gasOilContactCapillaryPressure();
|
||||
}
|
||||
|
||||
int EquilReg::equilibrationAccuracy() const
|
||||
{
|
||||
return this->rec_.initializationTargetAccuracy();
|
||||
}
|
||||
|
||||
const EquilReg::CalcDissolution&
|
||||
EquilReg::dissolutionCalculator() const
|
||||
{
|
||||
return *this->rs_;
|
||||
}
|
||||
|
||||
const EquilReg::CalcEvaporation&
|
||||
EquilReg::evaporationCalculator() const
|
||||
{
|
||||
return *this->rv_;
|
||||
}
|
||||
|
||||
const EquilReg::CalcWaterEvaporation&
|
||||
EquilReg::waterEvaporationCalculator() const
|
||||
{
|
||||
return *this->rvw_;
|
||||
}
|
||||
|
||||
const EquilReg::TabulatedFunction&
|
||||
EquilReg::saltVdTable() const
|
||||
{
|
||||
return saltVdTable_;
|
||||
}
|
||||
|
||||
int EquilReg::pvtIdx() const
|
||||
{
|
||||
return this->pvtIdx_;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
PcEq<FluidSystem,MaterialLawManager>::
|
||||
PcEq(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
targetPc_(targetPc)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double PcEq<FluidSystem,MaterialLawManager>::
|
||||
operator()(double s) const
|
||||
{
|
||||
const auto& matParams = materialLawManager_.materialLawParams(cell_);
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(phase_, s);
|
||||
|
||||
std::array<double, FluidSystem::numPhases> pc{0.0};
|
||||
using MaterialLaw = typename MaterialLawManager::MaterialLaw;
|
||||
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
|
||||
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
|
||||
return pcPhase - targetPc_;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
PcEqSum<FluidSystem,MaterialLawManager>::
|
||||
PcEqSum(const MaterialLawManager& materialLawManager,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
targetPc_(targetPc)
|
||||
{
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double PcEqSum<FluidSystem,MaterialLawManager>::
|
||||
operator()(double s) const
|
||||
{
|
||||
const auto& matParams = materialLawManager_.materialLawParams(cell_);
|
||||
SatOnlyFluidState fluidState;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
|
||||
fluidState.setSaturation(phase1_, s);
|
||||
fluidState.setSaturation(phase2_, 1.0 - s);
|
||||
|
||||
std::array<double, FluidSystem::numPhases> pc {0.0};
|
||||
|
||||
using MaterialLaw = typename MaterialLawManager::MaterialLaw;
|
||||
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
|
||||
double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
|
||||
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
|
||||
return pc1 + pc2 - targetPc_;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double minSaturations(const MaterialLawManager& materialLawManager,
|
||||
const int phase, const int cell)
|
||||
{
|
||||
const auto& scaledDrainageInfo =
|
||||
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
|
||||
|
||||
// Find minimum and maximum saturations.
|
||||
switch(phase) {
|
||||
case FluidSystem::waterPhaseIdx:
|
||||
return scaledDrainageInfo.Swl;
|
||||
|
||||
case FluidSystem::gasPhaseIdx:
|
||||
return scaledDrainageInfo.Sgl;
|
||||
|
||||
case FluidSystem::oilPhaseIdx:
|
||||
throw std::runtime_error("Min saturation not implemented for oil phase.");
|
||||
|
||||
default:
|
||||
throw std::runtime_error("Unknown phaseIdx .");
|
||||
}
|
||||
return -1.0;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double maxSaturations(const MaterialLawManager& materialLawManager,
|
||||
const int phase, const int cell)
|
||||
{
|
||||
const auto& scaledDrainageInfo =
|
||||
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
|
||||
|
||||
// Find minimum and maximum saturations.
|
||||
switch(phase) {
|
||||
case FluidSystem::waterPhaseIdx:
|
||||
return scaledDrainageInfo.Swu;
|
||||
|
||||
case FluidSystem::gasPhaseIdx:
|
||||
return scaledDrainageInfo.Sgu;
|
||||
|
||||
case FluidSystem::oilPhaseIdx:
|
||||
throw std::runtime_error("Max saturation not implemented for oil phase.");
|
||||
|
||||
default:
|
||||
throw std::runtime_error("Unknown phaseIdx .");
|
||||
}
|
||||
return -1.0;
|
||||
}
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double satFromPc(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double targetPc,
|
||||
const bool increasing)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
// Create the equation f(s) = pc(s) - targetPc
|
||||
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, targetPc);
|
||||
double f0 = f(s0);
|
||||
double f1 = f(s1);
|
||||
if (!std::isfinite(f0 + f1))
|
||||
throw std::logic_error(fmt::format("The capillary pressure values {} and {} are not finite", f0, f1));
|
||||
|
||||
if (f0 <= 0.0)
|
||||
return s0;
|
||||
else if (f1 >= 0.0)
|
||||
return s1;
|
||||
|
||||
const double tol = 1e-10;
|
||||
// should at least converge in 2 times bisection but some safety here:
|
||||
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
|
||||
int usedIterations = -1;
|
||||
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
|
||||
return root;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double targetPc)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
double s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - targetPc
|
||||
const PcEqSum<FluidSystem, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, targetPc);
|
||||
double f0 = f(s0);
|
||||
double f1 = f(s1);
|
||||
if (f0 <= 0.0)
|
||||
return s0;
|
||||
else if (f1 >= 0.0)
|
||||
return s1;
|
||||
|
||||
assert(f0 > 0.0 && f1 < 0.0);
|
||||
const double tol = 1e-10;
|
||||
// should at least converge in 2 times bisection but some safety here:
|
||||
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
|
||||
int usedIterations = -1;
|
||||
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
|
||||
return root;
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
double satFromDepth(const MaterialLawManager& materialLawManager,
|
||||
const double cellDepth,
|
||||
const double contactDepth,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const bool increasing)
|
||||
{
|
||||
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
if (cellDepth < contactDepth) {
|
||||
return s0;
|
||||
}
|
||||
else {
|
||||
return s1;
|
||||
}
|
||||
}
|
||||
|
||||
template<class FluidSystem, class MaterialLawManager>
|
||||
bool isConstPc(const MaterialLawManager& materialLawManager,
|
||||
const int phase,
|
||||
const int cell)
|
||||
{
|
||||
// Create the equation f(s) = pc(s);
|
||||
const PcEq<FluidSystem, MaterialLawManager> f(materialLawManager, phase, cell, 0);
|
||||
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
1950
ebos/equil/initstateequil_impl.hh
Normal file
1950
ebos/equil/initstateequil_impl.hh
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user