mirror of
https://github.com/OPM/opm-simulators.git
synced 2026-11-02 20:30:55 -06:00
Markus Blatt
55fb6fa7ed
The constructor now needs and additional boolean to flag whether to local permeabilities. The old one used is not there anymore. Therefor this patch moves the code to the new constructor interface.
353 lines
13 KiB
C++
353 lines
13 KiB
C++
/*
|
|
Copyright 2014 Andreas Lauser
|
|
|
|
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/>.
|
|
*/
|
|
#include <config.h>
|
|
|
|
#if HAVE_DYNAMIC_BOOST_TEST
|
|
#define BOOST_TEST_DYN_LINK
|
|
#endif
|
|
|
|
#define BOOST_TEST_MODULE TransmissibilityMultipliers
|
|
|
|
#include <opm/autodiff/GeoProps.hpp>
|
|
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
|
|
|
|
#include <opm/core/grid/GridManager.hpp>
|
|
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
|
|
|
#include <opm/parser/eclipse/Parser/Parser.hpp>
|
|
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
|
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
|
|
|
#if HAVE_DUNE_CORNERPOINT
|
|
#include <dune/common/version.hh>
|
|
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,3)
|
|
#include <dune/common/parallel/mpihelper.hh>
|
|
#else
|
|
#include <dune/common/mpihelper.hh>
|
|
#endif
|
|
#include <dune/grid/CpGrid.hpp>
|
|
#endif
|
|
|
|
#include <boost/test/unit_test.hpp>
|
|
|
|
#include <string>
|
|
|
|
#include <stdlib.h>
|
|
|
|
// as surprising as it seems, this is a minimal deck required to get to the point where
|
|
// the transmissibilities are calculated. The problem here is that the OPM property
|
|
// objects mix fluid and rock properties, so that all properties need to be defined even
|
|
// if they are not of interest :/
|
|
std::string deckPreMult =
|
|
"RUNSPEC\n"
|
|
"TABDIMS\n"
|
|
"/\n"
|
|
"OIL\n"
|
|
"GAS\n"
|
|
"WATER\n"
|
|
"METRIC\n"
|
|
"DIMENS\n"
|
|
"2 2 2/\n"
|
|
"GRID\n"
|
|
"DXV\n"
|
|
"1.0 2.0 /\n"
|
|
"DYV\n"
|
|
"3.0 4.0 /\n"
|
|
"DZV\n"
|
|
"5.0 6.0/\n"
|
|
"TOPS\n"
|
|
"4*100 /\n";
|
|
std::string deckPostMult =
|
|
"PROPS\n"
|
|
"DENSITY\n"
|
|
"100 200 300 /\n"
|
|
"PVTW\n"
|
|
" 100 1 1e-6 1.0 0 /\n"
|
|
"PVDG\n"
|
|
"1 1 1e-2\n"
|
|
"100 0.25 2e-2 /\n"
|
|
"PVTO\n"
|
|
"1e-3 1.0 1.0 1.0\n"
|
|
" 100.0 1.0 1.0\n"
|
|
"/\n"
|
|
"1.0 10.0 1.1 0.9\n"
|
|
" 100.0 1.1 0.9\n"
|
|
"/\n"
|
|
"/\n"
|
|
"SWOF\n"
|
|
"0.0 0.0 1.0 0.0\n"
|
|
"1.0 1.0 0.0 1.0/\n"
|
|
"SGOF\n"
|
|
"0.0 0.0 1.0 0.0\n"
|
|
"1.0 1.0 0.0 1.0/\n"
|
|
"PORO\n"
|
|
"8*0.3 /\n"
|
|
"PERMX\n"
|
|
"8*1 /\n"
|
|
"SCHEDULE\n"
|
|
"TSTEP\n"
|
|
"1.0 2.0 3.0 4.0 /\n"
|
|
"/\n";
|
|
|
|
std::string origDeckString = deckPreMult + deckPostMult;
|
|
std::string multDeckString =
|
|
deckPreMult +
|
|
"MULTX\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
"MULTY\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
"MULTZ\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
deckPostMult;
|
|
|
|
std::string multMinusDeckString =
|
|
deckPreMult +
|
|
"MULTX-\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
"MULTY-\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
"MULTZ-\n" +
|
|
"1 2 3 4 5 6 7 8 /\n" +
|
|
deckPostMult;
|
|
|
|
// the NTG values get harmonically averaged for the transmissibilites. If the value
|
|
// is the same on both sides, the averaging buils down to a no-op, though...
|
|
std::string ntgDeckString =
|
|
deckPreMult +
|
|
"NTG\n" +
|
|
"8*0.5 /\n" +
|
|
deckPostMult;
|
|
|
|
|
|
BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
|
|
{
|
|
Opm::parameter::ParameterGroup param;
|
|
Opm::ParserPtr parser(new Opm::Parser() );
|
|
|
|
/////
|
|
// create a DerivedGeology object without any multipliers involved
|
|
Opm::DeckConstPtr origDeck = parser->parseString(origDeckString);
|
|
Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck));
|
|
|
|
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, false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object _with_ transmissibility multipliers involved
|
|
Opm::DeckConstPtr multDeck = parser->parseString(multDeckString);
|
|
Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck));
|
|
|
|
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, false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object _with_ transmissibility multipliers involved for
|
|
// the negative faces
|
|
Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString);
|
|
Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck));
|
|
|
|
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, false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object with the NTG keyword involved
|
|
Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString);
|
|
Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck));
|
|
|
|
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, false);
|
|
/////
|
|
|
|
// compare the transmissibilities (note that for this we assume that the multipliers
|
|
// do not change the grid topology)
|
|
const UnstructuredGrid cGrid = *(origGridManager->c_grid());
|
|
int numFaces = cGrid.number_of_faces;
|
|
for (int faceIdx = 0; faceIdx < numFaces; ++ faceIdx) {
|
|
// in DUNE-speak, a face here is more like an intersection which is not specific
|
|
// to a codim-0 entity (i.e., cell)
|
|
|
|
// get the cell indices of the compressed grid for the face's interior and
|
|
// exterior cell
|
|
int insideCellIdx = cGrid.face_cells[2*faceIdx + 0];
|
|
int outsideCellIdx = cGrid.face_cells[2*faceIdx + 1];
|
|
|
|
if (insideCellIdx < 0 || outsideCellIdx < 0) {
|
|
// do not consider cells at the domain boundary: Their would only be used for
|
|
// Dirichlet-like boundary conditions which have not been implemented so
|
|
// far...
|
|
continue;
|
|
}
|
|
|
|
// translate these to canonical indices (i.e., the logically Cartesian ones used by the deck)
|
|
if (cGrid.global_cell) {
|
|
insideCellIdx = cGrid.global_cell[insideCellIdx];
|
|
outsideCellIdx = cGrid.global_cell[outsideCellIdx];
|
|
}
|
|
|
|
double origTrans = origGeology.transmissibility()[faceIdx];
|
|
double multTrans = multGeology.transmissibility()[faceIdx];
|
|
double multMinusTrans = multMinusGeology.transmissibility()[faceIdx];
|
|
double ntgTrans = ntgGeology.transmissibility()[faceIdx];
|
|
BOOST_CHECK_CLOSE(origTrans*(insideCellIdx + 1), multTrans, 1e-6);
|
|
BOOST_CHECK_CLOSE(origTrans*(outsideCellIdx + 1), multMinusTrans, 1e-6);
|
|
|
|
int insideCellKIdx = insideCellIdx/(cGrid.cartdims[0]*cGrid.cartdims[1]);
|
|
int outsideCellKIdx = outsideCellIdx/(cGrid.cartdims[0]*cGrid.cartdims[1]);
|
|
if (insideCellKIdx == outsideCellKIdx)
|
|
// NTG only reduces the permebility of the X-Y plane
|
|
BOOST_CHECK_CLOSE(origTrans*0.5, ntgTrans, 1e-6);
|
|
}
|
|
}
|
|
|
|
#if HAVE_DUNE_CORNERPOINT
|
|
BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
|
|
{
|
|
int argc = 1;
|
|
char **argv;
|
|
argv = new (char*);
|
|
argv[0] = strdup("footest");
|
|
|
|
Dune::MPIHelper::instance(argc, argv);
|
|
|
|
Opm::parameter::ParameterGroup param;
|
|
Opm::ParserPtr parser(new Opm::Parser() );
|
|
|
|
#warning TODO: there seems to be some index mess-up in DerivedGeology in the case of Dune::CpGrid
|
|
#if 0
|
|
|
|
/////
|
|
// create a DerivedGeology object without any multipliers involved
|
|
Opm::DeckConstPtr origDeck = parser->parseString(origDeckString);
|
|
Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck));
|
|
|
|
auto origGrid = std::make_shared<Dune::CpGrid>();
|
|
origGrid->processEclipseFormat(origEclipseState->getEclipseGrid(), 0.0, false);
|
|
|
|
auto origProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(origDeck,
|
|
origEclipseState,
|
|
*origGrid);
|
|
|
|
Opm::DerivedGeology origGeology(*origGrid, *origProps, origEclipseState, false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object _with_ transmissibility multipliers involved
|
|
Opm::DeckConstPtr multDeck = parser->parseString(multDeckString);
|
|
Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck));
|
|
|
|
auto multGrid = std::make_shared<Dune::CpGrid>();
|
|
multGrid->processEclipseFormat(multEclipseState->getEclipseGrid(), 0.0, false);
|
|
|
|
auto multProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multDeck, multEclipseState, *multGrid);
|
|
|
|
Opm::DerivedGeology multGeology(*multGrid, *multProps, multEclipseState, false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object _with_ transmissibility multipliers involved for
|
|
// the negative faces
|
|
Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString);
|
|
Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck));
|
|
|
|
auto multMinusGrid = std::make_shared<Dune::CpGrid>();
|
|
multMinusGrid->processEclipseFormat(multMinusEclipseState->getEclipseGrid(), 0.0, false);
|
|
|
|
auto multMinusProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multMinusDeck, multMinusEclipseState, *multMinusGrid);
|
|
|
|
Opm::DerivedGeology multMinusGeology(*multMinusGrid, *multMinusProps, multMinusEclipseState,
|
|
false);
|
|
/////
|
|
|
|
/////
|
|
// create a DerivedGeology object with the NTG keyword involved
|
|
Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString);
|
|
Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck));
|
|
|
|
auto ntgGrid = std::make_shared<Dune::CpGrid>();
|
|
ntgGrid->processEclipseFormat(ntgEclipseState->getEclipseGrid(), 0.0, false);
|
|
|
|
auto ntgProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(ntgDeck, ntgEclipseState, *ntgGrid);
|
|
|
|
Opm::DerivedGeology ntgGeology(*ntgGrid, *ntgProps, ntgEclipseState);
|
|
/////
|
|
|
|
// compare the transmissibilities (note that for this we assume that the multipliers
|
|
// do not change the grid topology)
|
|
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,3)
|
|
auto gridView = origGrid->leafGridView();
|
|
#else
|
|
auto gridView = origGrid->leafView();
|
|
#endif
|
|
typedef Dune::MultipleCodimMultipleGeomTypeMapper<Dune::CpGrid::LeafGridView,
|
|
Dune::MCMGElementLayout> ElementMapper;
|
|
ElementMapper elementMapper(gridView);
|
|
auto eIt = gridView.begin<0>();
|
|
const auto& eEndIt = gridView.end<0>();
|
|
for (; eIt < eEndIt; ++ eIt) {
|
|
// loop over the intersections of the current element
|
|
auto isIt = gridView.ibegin(*eIt);
|
|
const auto& isEndIt = gridView.iend(*eIt);
|
|
for (; isIt != isEndIt; ++isIt) {
|
|
if (isIt->boundary())
|
|
continue; // ignore domain the boundaries
|
|
|
|
// get the cell indices of the compressed grid for the face's interior and
|
|
// exterior cell
|
|
int insideCellIdx = elementMapper.map(*isIt->inside());
|
|
int outsideCellIdx = elementMapper.map(*isIt->outside());
|
|
|
|
// translate these to canonical indices (i.e., the logically Cartesian ones used by the deck)
|
|
insideCellIdx = origGrid->globalCell()[insideCellIdx];
|
|
outsideCellIdx = origGrid->globalCell()[outsideCellIdx];
|
|
|
|
#warning TODO: how to get the intersection index for the compressed grid??
|
|
#if 0
|
|
int globalIsIdx = 0; // <- how to get this??
|
|
double origTrans = origGeology.transmissibility()[globalIsIdx];
|
|
double multTrans = multGeology.transmissibility()[globalIsIdx];
|
|
double multMinusTrans = multMinusGeology.transmissibility()[globalIsIdx];
|
|
double ntgTrans = ntgGeology.transmissibility()[globalIsIdx];
|
|
BOOST_CHECK_CLOSE(origTrans*(insideCellIdx + 1), multTrans, 1e-6);
|
|
BOOST_CHECK_CLOSE(origTrans*(outsideCellIdx + 1), multMinusTrans, 1e-6);
|
|
|
|
std::array<int, 3> ijkInside, ijkOutside;
|
|
origGrid->getIJK(insideCellIdx, ijkInside);
|
|
origGrid->getIJK(outsideCellIdx, ijkOutside);
|
|
if (ijkInside[2] == ijkOutside[2])
|
|
// NTG only reduces the permebility of the X-Y plane
|
|
BOOST_CHECK_CLOSE(origTrans*0.5, ntgTrans, 1e-6);
|
|
#endif
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
#endif
|