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 transmissibilities using the generic GridHelpers interface.

Browse Source
This commit is contained in:
Markus Blatt 2015-01-13 20:17:30 +01:00
parent 656e688692
commit a3acd22f24
1 changed files with 41 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

94
opm/autodiff/GeoProps.hpp
View File

@ -196,35 +196,40 @@ namespace Opm
}; };
template <> template <class GridType>
inline void DerivedGeology::multiplyHalfIntersections_<UnstructuredGrid>(const UnstructuredGrid &grid, inline void DerivedGeology::multiplyHalfIntersections_(const GridType &grid,
Opm::EclipseStateConstPtr eclState, Opm::EclipseStateConstPtr eclState,
const std::vector<double> &ntg, const std::vector<double> &ntg,
Vector &halfIntersectTransmissibility, Vector &halfIntersectTransmissibility,
std::vector<double> &intersectionTransMult) std::vector<double> &intersectionTransMult)
{ {
int numCells = grid.number_of_cells; int numCells = Opm::AutoDiffGrid::numCells(grid);
int numIntersections = grid.number_of_faces; int numIntersections = Opm::AutoDiffGrid::numFaces(grid);
intersectionTransMult.resize(numIntersections); intersectionTransMult.resize(numIntersections);
std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0); std::fill(intersectionTransMult.begin(), intersectionTransMult.end(), 1.0);
std::shared_ptr<const Opm::TransMult> multipliers = eclState->getTransMult(); std::shared_ptr<const Opm::TransMult> multipliers = eclState->getTransMult();
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
auto faceCells = Opm::AutoDiffGrid::faceCells(grid);
const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
int cellFaceIdx = 0;
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
// loop over all logically-Cartesian faces of the current cell // loop over all logically-Cartesian faces of the current cell
for (int cellFaceIdx = grid.cell_facepos[cellIdx]; auto cellFacesRange = cell2Faces[cellIdx];
cellFaceIdx < grid.cell_facepos[cellIdx + 1];
++ cellFaceIdx) for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
{ {
// the index of the current cell in arrays for the logically-Cartesian grid // the index of the current cell in arrays for the logically-Cartesian grid
int cartesianCellIdx = grid.global_cell[cellIdx]; int cartesianCellIdx = global_cell[cellIdx];
// The index of the face in the compressed grid // The index of the face in the compressed grid
int faceIdx = grid.cell_faces[cellFaceIdx]; int faceIdx = *cellFaceIter;
// the logically-Cartesian direction of the face // the logically-Cartesian direction of the face
int faceTag = grid.cell_facetag[cellFaceIdx]; int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
// Translate the C face tag into the enum used by opm-parser's TransMult class // Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection; Opm::FaceDir::DirEnum faceDirection;
@ -261,15 +266,15 @@ namespace Opm
multipliers->getMultiplier(cartesianCellIdx, faceDirection); multipliers->getMultiplier(cartesianCellIdx, faceDirection);
// Multiplier contribution on this fase for region multipliers // Multiplier contribution on this fase for region multipliers
const int cellIdxInside = grid.face_cells[2*faceIdx]; const int cellIdxInside = faceCells(faceIdx, 0);
const int cellIdxOutside = grid.face_cells[2*faceIdx + 1]; const int cellIdxOutside = faceCells(faceIdx, 1);
// Do not apply region multipliers in the case of boundary connections // Do not apply region multipliers in the case of boundary connections
if (cellIdxInside < 0 || cellIdxOutside < 0) { if (cellIdxInside < 0 || cellIdxOutside < 0) {
continue; continue;
} }
const int cartesianCellIdxInside = grid.global_cell[cellIdxInside]; const int cartesianCellIdxInside = global_cell[cellIdxInside];
const int cartesianCellIdxOutside = grid.global_cell[cellIdxOutside]; const int cartesianCellIdxOutside = global_cell[cellIdxOutside];
// Only apply the region multipliers from the inside // Only apply the region multipliers from the inside
if (cartesianCellIdx == cartesianCellIdxInside) { if (cartesianCellIdx == cartesianCellIdxInside) {
intersectionTransMult[faceIdx] *= multipliers->getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection); intersectionTransMult[faceIdx] *= multipliers->getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,faceDirection);
@ -281,17 +286,9 @@ namespace Opm
} }
template <class GridType> template <class GridType>
inline void DerivedGeology::tpfa_loc_trans_compute_(const GridType& /* grid */, inline void DerivedGeology::tpfa_loc_trans_compute_(const GridType& grid,
const double* /* perm */, const double* perm,
Vector& /* hTrans */) Vector& hTrans){
{
OPM_THROW(std::logic_error, "No generic implementation exists for tpfa_loc_trans_compute_().");
}
template <>
inline void DerivedGeology::tpfa_loc_trans_compute_<UnstructuredGrid>(const UnstructuredGrid& grid,
const double* perm,
Vector& hTrans){
// Using Local coordinate system for the transmissibility calculations // Using Local coordinate system for the transmissibility calculations
// hTrans(cellFaceIdx) = K(cellNo,j) * sum( C(:,i) .* N(:,j), 2) / sum(C.*C, 2) // hTrans(cellFaceIdx) = K(cellNo,j) * sum( C(:,i) .* N(:,j), 2) / sum(C.*C, 2)
@ -299,17 +296,22 @@ namespace Opm
// to face centroid and N is the normal vector pointing outwards with norm equal to the face area. // to face centroid and N is the normal vector pointing outwards with norm equal to the face area.
// Off-diagonal permeability values are ignored without warning // Off-diagonal permeability values are ignored without warning
int numCells = AutoDiffGrid::numCells(grid); int numCells = AutoDiffGrid::numCells(grid);
int cellFaceIdx = 0;
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
auto faceCells = Opm::UgGridHelpers::faceCells(grid);
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
// loop over all logically-Cartesian faces of the current cell // loop over all logically-Cartesian faces of the current cell
for (int cellFaceIdx = grid.cell_facepos[cellIdx]; auto cellFacesRange = cell2Faces[cellIdx];
cellFaceIdx < grid.cell_facepos[cellIdx + 1];
++ cellFaceIdx) for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
{ {
// The index of the face in the compressed grid // The index of the face in the compressed grid
const int faceIdx = grid.cell_faces[cellFaceIdx]; const int faceIdx = *cellFaceIter;
// the logically-Cartesian direction of the face // the logically-Cartesian direction of the face
const int faceTag = grid.cell_facetag[cellFaceIdx]; const int faceTag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
// d = 0: XPERM d = 4: YPERM d = 8: ZPERM ignores off-diagonal permeability values. // d = 0: XPERM d = 4: YPERM d = 8: ZPERM ignores off-diagonal permeability values.
const int d = std::floor(faceTag/2) * 4; const int d = std::floor(faceTag/2) * 4;
@ -317,12 +319,13 @@ namespace Opm
// compute the half transmissibility // compute the half transmissibility
double dist = 0.0; double dist = 0.0;
double cn = 0.0; double cn = 0.0;
double sgn = 2.0 * (grid.face_cells[2*faceIdx + 0] == cellIdx) - 1; double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1;
const int dim = grid.dimensions; const int dim = Opm::UgGridHelpers::dimensions(grid);
for (int indx = 0; indx < dim; ++indx) { for (int indx = 0; indx < dim; ++indx) {
const double Ci = grid.face_centroids[faceIdx*dim + indx] - grid.cell_centroids[cellIdx*dim + indx]; const double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx] -
Opm::UgGridHelpers::cellCentroidCoordinate(grid, cellIdx, indx);
dist += Ci*Ci; dist += Ci*Ci;
cn += sgn * Ci * grid.face_normals[faceIdx*dim + indx]; cn += sgn * Ci * Opm::UgGridHelpers::faceNormal(grid, faceIdx)[indx];
} }
if (cn < 0){ if (cn < 0){
@ -348,21 +351,6 @@ namespace Opm
} }
#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
} }