OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-26 17:20:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #249 from totto82/locTrans

Browse Source 
Compute half transmissibilities based on local coordinate system
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-17 00:20:11 +01:00
commit 1348658507
5 changed files with 99 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/opm_init_check.cpp
View File

@ -212,7 +212,7 @@ void initOPMTrans(TransGraph& opmTrans , DeckConstPtr deck , std::shared_ptr<con
std::shared_ptr<BlackoilPropsAdInterface> props;
props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
DerivedGeology geology(*grid->c_grid() , *props, eclipseState);
DerivedGeology geology(*grid->c_grid() , *props, eclipseState, false);
const double * opm_trans_data = geology.transmissibility().data();
double SIconversion = Opm::unit::cubic(Opm::unit::meter) * Opm::unit::day * Opm::unit::barsa / (Opm::prefix::centi * Opm::unit::Poise);

2
examples/sim_2p_comp_ad.cpp
View File

@ -154,7 +154,7 @@ try
// With a deck, we may have more epochs etc.
WellState well_state;
Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();
Opm::DerivedGeology geology(*grid->c_grid(), *props, eclipseState);
Opm::DerivedGeology geology(*grid->c_grid(), *props, eclipseState,false);
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {

3
examples/sim_fibo_ad.cpp
View File

@ -207,7 +207,8 @@ try
// initialize variables
simtimer.init(timeMap);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
bool use_local_perm = param.getDefault("use_local_perm", true);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, use_local_perm, grav);
// TODO: also do this in sim_fibo_ad_cp
bool writeTrans = param.getDefault("write_transmissibilities", false);

93
opm/autodiff/GeoProps.hpp
View File

@ -61,7 +61,10 @@ namespace Opm
DerivedGeology(const Grid& grid,
const Props& props ,
Opm::EclipseStateConstPtr eclState,
const double* grav = 0)
const bool use_local_perm,
const double* grav = 0
)
: pvol_ (Opm::AutoDiffGrid::numCells(grid))
, trans_(Opm::AutoDiffGrid::numFaces(grid))
, gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
@ -98,9 +101,24 @@ namespace Opm
}
// Transmissibility
Vector htrans(AutoDiffGrid::numCellFaces(grid));
Grid* ug = const_cast<Grid*>(& grid);
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
#ifdef HAVE_DUNE_CORNERPOINT
if (std::is_same<Grid, Dune::CpGrid>::value) {
if (use_local_perm) {
OPM_THROW(std::runtime_error, "Local coordinate permeability not supported for CpGrid");
}
}
#endif
if (! use_local_perm) {
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
}
else {
tpfa_loc_trans_compute_(grid,props.permeability(),htrans);
}
std::vector<double> mult;
multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
@ -161,13 +179,23 @@ namespace Opm
Vector &halfIntersectTransmissibility,
std::vector<double> &intersectionTransMult);
template <class Grid>
void tpfa_loc_trans_compute_(const Grid &grid,
const double* perm,
Vector &hTrans);
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,
@ -252,6 +280,67 @@ namespace Opm
}
}
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
// hTrans(cellFaceIdx) = K(cellNo,j) * sum( C(:,i) .* N(:,j), 2) / sum(C.*C, 2)
// where K is a diagonal permeability tensor, C is the distance from cell centroid
// 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
int numCells = AutoDiffGrid::numCells(grid);
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 face in the compressed grid
const int faceIdx = grid.cell_faces[cellFaceIdx];
// the logically-Cartesian direction of the face
const int faceTag = grid.cell_facetag[cellFaceIdx];
// d = 0: XPERM d = 4: YPERM d = 8: ZPERM ignores off-diagonal permeability values.
const int d = std::floor(faceTag/2) * 4;
// compute the half transmissibility
double dist = 0.0;
double cn = 0.0;
double sgn = 2.0 * (grid.face_cells[2*faceIdx + 0] == cellIdx) - 1;
const int dim = grid.dimensions;
for (int indx = 0; indx < dim; ++indx) {
const double Ci = grid.face_centroids[faceIdx*dim + indx] - grid.cell_centroids[cellIdx*dim + indx];
dist += Ci*Ci;
cn += sgn * Ci * grid.face_normals[faceIdx*dim + indx];
}
if (cn < 0){
switch (d) {
case 0:
OPM_MESSAGE("Warning: negative X-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
break;
case 4:
OPM_MESSAGE("Warning: negative Y-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
break;
case 8:
OPM_MESSAGE("Warning: negative Z-transmissibility value in cell: " << cellIdx << " replace by absolute value") ;
break;
default:
OPM_THROW(std::logic_error, "Inconsistency in the faceTag in cell: " << cellIdx);
}
cn = -cn;
}
hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist;
}
}
}
#ifdef HAVE_DUNE_CORNERPOINT
template <>
inline void DerivedGeology::multiplyHalfIntersections_<Dune::CpGrid>(const Dune::CpGrid &grid,

8
tests/test_transmissibilitymultipliers.cpp
View File

@ -148,7 +148,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
auto origGridManager = std::make_shared<Opm::GridManager>(origEclipseState->getEclipseGrid());
auto origProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(origDeck, origEclipseState, *(origGridManager->c_grid()));
Opm::DerivedGeology origGeology(*(origGridManager->c_grid()), *origProps, origEclipseState);
Opm::DerivedGeology origGeology(*(origGridManager->c_grid()), *origProps, origEclipseState, false);
/////
/////
@ -159,7 +159,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
auto multGridManager = std::make_shared<Opm::GridManager>(multEclipseState->getEclipseGrid());
auto multProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multDeck, multEclipseState, *(multGridManager->c_grid()));
Opm::DerivedGeology multGeology(*(multGridManager->c_grid()), *multProps, multEclipseState);
Opm::DerivedGeology multGeology(*(multGridManager->c_grid()), *multProps, multEclipseState, false);
/////
/////
@ -171,7 +171,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
auto multMinusGridManager = std::make_shared<Opm::GridManager>(multMinusEclipseState->getEclipseGrid());
auto multMinusProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multMinusDeck, multMinusEclipseState, *(multMinusGridManager->c_grid()));
Opm::DerivedGeology multMinusGeology(*(multMinusGridManager->c_grid()), *multMinusProps, multMinusEclipseState);
Opm::DerivedGeology multMinusGeology(*(multMinusGridManager->c_grid()), *multMinusProps, multMinusEclipseState, false);
/////
/////
@ -182,7 +182,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
auto ntgGridManager = std::make_shared<Opm::GridManager>(ntgEclipseState->getEclipseGrid());
auto ntgProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(ntgDeck, ntgEclipseState, *(ntgGridManager->c_grid()));
Opm::DerivedGeology ntgGeology(*(ntgGridManager->c_grid()), *ntgProps, ntgEclipseState);
Opm::DerivedGeology ntgGeology(*(ntgGridManager->c_grid()), *ntgProps, ntgEclipseState, false);
/////
// compare the transmissibilities (note that for this we assume that the multipliers