OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

implement pore volume and transmissibility multipliers

Browse Source 
this means that the NTG, MULTPV and MULT[XYZ]-? keywords are now
supported.

Actually FAULTS and MULTFAULT are supported too, but that's abstracted
away by opm-parser's TransMult class. (Kudos to [at]joakim-hove for
implementing this.)
This commit is contained in:
Andreas Lauser 2014-07-10 11:53:50 +02:00
parent 63eaecf246
commit e07ef16a45
1 changed files with 145 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

161
opm/autodiff/GeoProps.hpp
View File

@ -33,6 +33,12 @@
#include "reenable_warning_pragmas.h"
#ifdef HAVE_DUNE_CORNERPOINT
#include <dune/common/version.hh>
#include <dune/grid/CpGrid.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#endif
#include <cstddef>
namespace Opm
@ -60,41 +66,76 @@ namespace Opm
, gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
, z_(Opm::AutoDiffGrid::numCells(grid))
{
using namespace Opm::AutoDiffGrid;
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
const typename Vector::Index nc = numCells(grid);
std::transform(beginCellVolumes(grid), endCellVolumes(grid),
props.porosity(), pvol_.data(),
std::multiplies<double>());
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(numCellFaces(grid));
Vector htrans(AutoDiffGrid::numCellFaces(grid));
Grid* ug = const_cast<Grid*>(& grid);
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
tpfa_trans_compute (ug, htrans.data() , trans_.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<nc; ++c){
z_[c] = cellCentroid(grid, c)[2];
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 = dimensions(grid);
typedef typename ADCell2FacesTraits<Grid>::Type Cell2Faces;
Cell2Faces c2f=cell2Faces(grid);
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 < nc; ++c) {
const double* const cc = cellCentroid(grid, c);
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 = faceCentroid(grid, *f);
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]);
@ -112,15 +153,103 @@ namespace Opm
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