mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-29 12:33:49 -06:00
257 lines
11 KiB
C++
257 lines
11 KiB
C++
/*
|
|
Copyright 2013 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_GEOPROPS_HEADER_INCLUDED
|
|
#define OPM_GEOPROPS_HEADER_INCLUDED
|
|
|
|
#include <opm/core/grid.h>
|
|
#include <opm/autodiff/GridHelpers.hpp>
|
|
#include <opm/core/utility/ErrorMacros.hpp>
|
|
//#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
|
#include <opm/core/pressure/tpfa/TransTpfa.hpp>
|
|
|
|
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
|
|
|
#include <opm/core/utility/platform_dependent/disable_warnings.h>
|
|
|
|
#include <Eigen/Eigen>
|
|
|
|
#ifdef HAVE_DUNE_CORNERPOINT
|
|
#include <dune/common/version.hh>
|
|
#include <dune/grid/CpGrid.hpp>
|
|
#include <dune/grid/common/mcmgmapper.hh>
|
|
#endif
|
|
|
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
|
|
|
#include <cstddef>
|
|
|
|
namespace Opm
|
|
{
|
|
|
|
/// Class containing static geological properties that are
|
|
/// derived from grid and petrophysical properties:
|
|
/// - pore volume
|
|
/// - transmissibilities
|
|
/// - gravity potentials
|
|
class DerivedGeology
|
|
{
|
|
public:
|
|
typedef Eigen::ArrayXd Vector;
|
|
|
|
/// Construct contained derived geological properties
|
|
/// from grid and property information.
|
|
template <class Props, class Grid>
|
|
DerivedGeology(const Grid& grid,
|
|
const Props& props ,
|
|
Opm::EclipseStateConstPtr eclState,
|
|
const double* grav = 0)
|
|
: pvol_ (Opm::AutoDiffGrid::numCells(grid))
|
|
, trans_(Opm::AutoDiffGrid::numFaces(grid))
|
|
, gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
|
|
, z_(Opm::AutoDiffGrid::numCells(grid))
|
|
{
|
|
int numCells = AutoDiffGrid::numCells(grid);
|
|
int numFaces = AutoDiffGrid::numFaces(grid);
|
|
const int *cartDims = AutoDiffGrid::cartDims(grid);
|
|
int numCartesianCells =
|
|
cartDims[0]
|
|
* cartDims[1]
|
|
* cartDims[2];
|
|
|
|
// get the pore volume multipliers from the EclipseState
|
|
std::vector<double> multpv(numCartesianCells, 1.0);
|
|
if (eclState->hasDoubleGridProperty("MULTPV")) {
|
|
multpv = eclState->getDoubleGridProperty("MULTPV")->getData();
|
|
}
|
|
|
|
// get the net-to-gross cell thickness from the EclipseState
|
|
std::vector<double> ntg(numCartesianCells, 1.0);
|
|
if (eclState->hasDoubleGridProperty("NTG")) {
|
|
ntg = eclState->getDoubleGridProperty("NTG")->getData();
|
|
}
|
|
|
|
// Pore volume
|
|
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
|
|
int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
|
|
pvol_[cellIdx] =
|
|
props.porosity()[cellIdx]
|
|
* multpv[cartesianCellIdx]
|
|
* ntg[cartesianCellIdx]
|
|
* AutoDiffGrid::cellVolume(grid, cellIdx);
|
|
}
|
|
|
|
// Transmissibility
|
|
Vector htrans(AutoDiffGrid::numCellFaces(grid));
|
|
Grid* ug = const_cast<Grid*>(& grid);
|
|
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
|
|
|
|
std::vector<double> mult;
|
|
multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
|
|
|
|
// combine the half-face transmissibilites into the final face
|
|
// transmissibilites.
|
|
tpfa_trans_compute(ug, htrans.data(), trans_.data());
|
|
|
|
// multiply the face transmissibilities with their appropriate
|
|
// transmissibility multipliers
|
|
for (int faceIdx = 0; faceIdx < numFaces; faceIdx++) {
|
|
trans_[faceIdx] *= mult[faceIdx];
|
|
}
|
|
|
|
// Compute z coordinates
|
|
for (int c = 0; c<numCells; ++c){
|
|
z_[c] = AutoDiffGrid::cellCentroid(grid, c)[2];
|
|
}
|
|
|
|
|
|
// Gravity potential
|
|
std::fill(gravity_, gravity_ + 3, 0.0);
|
|
if (grav != 0) {
|
|
const typename Vector::Index nd = AutoDiffGrid::dimensions(grid);
|
|
typedef typename AutoDiffGrid::ADCell2FacesTraits<Grid>::Type Cell2Faces;
|
|
Cell2Faces c2f=AutoDiffGrid::cell2Faces(grid);
|
|
|
|
std::size_t i = 0;
|
|
for (typename Vector::Index c = 0; c < numCells; ++c) {
|
|
const double* const cc = AutoDiffGrid::cellCentroid(grid, c);
|
|
|
|
typename Cell2Faces::row_type faces=c2f[c];
|
|
typedef typename Cell2Faces::row_type::iterator Iter;
|
|
|
|
for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f, ++i) {
|
|
const double* const fc = AutoDiffGrid::faceCentroid(grid, *f);
|
|
|
|
for (typename Vector::Index d = 0; d < nd; ++d) {
|
|
gpot_[i] += grav[d] * (fc[d] - cc[d]);
|
|
}
|
|
}
|
|
}
|
|
std::copy(grav, grav + nd, gravity_);
|
|
}
|
|
}
|
|
|
|
const Vector& poreVolume() const { return pvol_ ;}
|
|
const Vector& transmissibility() const { return trans_ ;}
|
|
const Vector& gravityPotential() const { return gpot_ ;}
|
|
const Vector& z() const { return z_ ;}
|
|
const double* gravity() const { return gravity_;}
|
|
|
|
private:
|
|
template <class Grid>
|
|
void multiplyHalfIntersections_(const Grid &grid,
|
|
Opm::EclipseStateConstPtr eclState,
|
|
const std::vector<double> &ntg,
|
|
Vector &halfIntersectTransmissibility,
|
|
std::vector<double> &intersectionTransMult);
|
|
|
|
Vector pvol_ ;
|
|
Vector trans_;
|
|
Vector gpot_ ;
|
|
Vector z_;
|
|
double gravity_[3]; // Size 3 even if grid is 2-dim.
|
|
};
|
|
|
|
template <>
|
|
inline void DerivedGeology::multiplyHalfIntersections_<UnstructuredGrid>(const UnstructuredGrid &grid,
|
|
Opm::EclipseStateConstPtr eclState,
|
|
const std::vector<double> &ntg,
|
|
Vector &halfIntersectTransmissibility,
|
|
std::vector<double> &intersectionTransMult)
|
|
{
|
|
int numCells = grid.number_of_cells;
|
|
|
|
int numIntersections = grid.number_of_faces;
|
|
intersectionTransMult.resize(numIntersections);
|
|
std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
|
|
|
|
std::shared_ptr<const Opm::TransMult> multipliers = eclState->getTransMult();
|
|
|
|
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
|
|
// loop over all logically-Cartesian faces of the current cell
|
|
for (int cellFaceIdx = grid.cell_facepos[cellIdx];
|
|
cellFaceIdx < grid.cell_facepos[cellIdx + 1];
|
|
++ cellFaceIdx)
|
|
{
|
|
// the index of the current cell in arrays for the logically-Cartesian grid
|
|
int cartesianCellIdx = grid.global_cell[cellIdx];
|
|
|
|
// The index of the face in the compressed grid
|
|
int faceIdx = grid.cell_faces[cellFaceIdx];
|
|
|
|
// the logically-Cartesian direction of the face
|
|
int faceTag = grid.cell_facetag[cellFaceIdx];
|
|
|
|
// Translate the C face tag into the enum used by opm-parser's TransMult class
|
|
Opm::FaceDir::DirEnum faceDirection;
|
|
if (faceTag == 0) // left
|
|
faceDirection = Opm::FaceDir::XMinus;
|
|
else if (faceTag == 1) // right
|
|
faceDirection = Opm::FaceDir::XPlus;
|
|
else if (faceTag == 2) // back
|
|
faceDirection = Opm::FaceDir::YMinus;
|
|
else if (faceTag == 3) // front
|
|
faceDirection = Opm::FaceDir::YPlus;
|
|
else if (faceTag == 4) // bottom
|
|
faceDirection = Opm::FaceDir::ZMinus;
|
|
else if (faceTag == 5) // top
|
|
faceDirection = Opm::FaceDir::ZPlus;
|
|
else
|
|
OPM_THROW(std::logic_error, "Unhandled face direction: " << faceTag);
|
|
|
|
// Account for NTG in horizontal one-sided transmissibilities
|
|
switch (faceDirection) {
|
|
case Opm::FaceDir::XMinus:
|
|
case Opm::FaceDir::XPlus:
|
|
case Opm::FaceDir::YMinus:
|
|
case Opm::FaceDir::YPlus:
|
|
halfIntersectTransmissibility[cellFaceIdx] *= ntg[cartesianCellIdx];
|
|
break;
|
|
default:
|
|
// do nothing for the top and bottom faces
|
|
break;
|
|
}
|
|
|
|
// Multiplier contribution on this face
|
|
intersectionTransMult[faceIdx] *=
|
|
multipliers->getMultiplier(cartesianCellIdx, faceDirection);
|
|
}
|
|
}
|
|
}
|
|
|
|
#ifdef HAVE_DUNE_CORNERPOINT
|
|
template <>
|
|
inline void DerivedGeology::multiplyHalfIntersections_<Dune::CpGrid>(const Dune::CpGrid &grid,
|
|
Opm::EclipseStateConstPtr eclState,
|
|
const std::vector<double> &ntg,
|
|
Vector &halfIntersectTransmissibility,
|
|
std::vector<double> &intersectionTransMult)
|
|
{
|
|
#warning "Transmissibility multipliers are not implemented for Dune::CpGrid due to difficulties in mapping intersections to unique indices."
|
|
int numIntersections = grid.numFaces();
|
|
intersectionTransMult.resize(numIntersections);
|
|
std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
|
|
}
|
|
#endif
|
|
}
|
|
|
|
|
|
|
|
#endif // OPM_GEOPROPS_HEADER_INCLUDED
|