mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-29 04:23:48 -06:00
25a4bb5fbe
The multiplier contribution from the getRegionMultiplier is added to the face multipliers. The getRegionMultiplier method is called with the cell index on both side of the face in order to return the correct region multiplier across faults.
274 lines
12 KiB
C++
274 lines
12 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 for MULT[XYZ] logical cartesian multipliers
|
|
intersectionTransMult[faceIdx] *=
|
|
multipliers->getMultiplier(cartesianCellIdx, faceDirection);
|
|
|
|
// Multiplier contribution on this fase for region multipliers
|
|
const int cellIdxInside = grid.face_cells[2*faceIdx];
|
|
const int cellIdxOutside = grid.face_cells[2*faceIdx + 1];
|
|
|
|
// Do not apply region multipliers in the case of boundary connections
|
|
if (cellIdxInside < 0 || cellIdxOutside < 0) {
|
|
continue;
|
|
}
|
|
const int cartesianCellIdxInside = grid.global_cell[cellIdxInside];
|
|
const int cartesianCellIdxOutside = grid.global_cell[cellIdxOutside];
|
|
// Only apply the region multipliers from the inside
|
|
if (cartesianCellIdx == cartesianCellIdxInside) {
|
|
intersectionTransMult[faceIdx] *= multipliers->getRegionMultiplier(cartesianCellIdxInside,cartesianCellIdxOutside,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
|