diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index 6b7f24a99..ca588b5bf 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -33,6 +33,12 @@ #include "reenable_warning_pragmas.h" +#ifdef HAVE_DUNE_CORNERPOINT +#include +#include +#include +#endif + #include 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 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 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()); + 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); tpfa_htrans_compute(ug, props.permeability(), htrans.data()); - tpfa_trans_compute (ug, htrans.data() , trans_.data()); + + std::vector 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::Type Cell2Faces; - Cell2Faces c2f=cell2Faces(grid); + const typename Vector::Index nd = AutoDiffGrid::dimensions(grid); + typedef typename AutoDiffGrid::ADCell2FacesTraits::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 + void multiplyHalfIntersections_(const Grid &grid, + Opm::EclipseStateConstPtr eclState, + const std::vector &ntg, + Vector &halfIntersectTransmissibility, + std::vector &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_(const UnstructuredGrid &grid, + Opm::EclipseStateConstPtr eclState, + const std::vector &ntg, + Vector &halfIntersectTransmissibility, + std::vector &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 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_(const Dune::CpGrid &grid, + Opm::EclipseStateConstPtr eclState, + const std::vector &ntg, + Vector &halfIntersectTransmissibility, + std::vector &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