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

Merge pull request #166 from andlaus/implement_multipliers

Browse Source 
Implement pore volume and permeability multipliers
This commit is contained in:
Bård Skaflestad 2014-07-28 14:38:56 +02:00
commit b98b71e401
12 changed files with 581 additions and 392 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
CMakeLists_files.cmake
View File

@ -49,6 +49,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_span.cpp
tests/test_syntax.cpp
tests/test_scalar_mult.cpp
tests/test_transmissibilitymultipliers.cpp
tests/test_welldensitysegmented.cpp
)
@ -75,8 +76,6 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/sim_2p_comp_ad.cpp
examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp
examples/test_impestpfa_ad.cpp
examples/test_implicit_ad.cpp
)
# programs listed here will not only be compiled, but also marked for

181
examples/sim_2p_comp_ad.cpp
View File

@ -45,6 +45,7 @@
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/SimulatorCompressibleAd.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
@ -85,13 +86,6 @@ try
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
if (!use_deck) {
// This check should be removed when and if this simulator is verified and works without decks.
// The current code for the non-deck case fails for unknown reasons.
OPM_THROW(std::runtime_error, "This simulator cannot run without a deck with wells. Use deck_filename to specify deck.");
}
Opm::DeckConstPtr deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<BlackoilPropertiesInterface> props;
@ -101,84 +95,33 @@ try
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile( deck_filename );
eclipseState.reset(new EclipseState(deck));
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
initBlackoilSurfvol(*grid->c_grid(), *props, state);
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Rock compressibility.
rock_comp.reset(new RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
deck = parser->parseFile( deck_filename );
eclipseState.reset(new EclipseState(deck));
// Grid init
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(deck));
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
initBlackoilSurfvol(*grid->c_grid(), *props, state);
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
// Initialising wells if not from deck.
Wells* simple_wells = 0;
if (!use_deck) {
// Compute pore volumes, in order to enable specifying injection rate
// terms of total pore volume.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/unit::day;
simple_wells = create_wells(2, 2, 2);
const double inj_frac[2] = { 1.0, 0.0 };
const int inj_cell = 0;
const double WI = 1e-8; // This is a completely made-up number.
const double all_fluids[2] = { 1.0, 1.0 };
int ok = add_well(INJECTOR, 0.0, 1, inj_frac, &inj_cell, &WI, "Injector", simple_wells);
ok = ok && append_well_controls(SURFACE_RATE, 0.01*flow_per_sec, all_fluids, 0, simple_wells);
const int prod_cell = grid->c_grid()->number_of_cells - 1;
ok = ok && add_well(PRODUCER, 0.0, 1, NULL, &prod_cell, &WI, "Producer", simple_wells);
ok = ok && append_well_controls(SURFACE_RATE, -0.01*flow_per_sec, all_fluids, 1, simple_wells);
if (!ok) {
OPM_THROW(std::runtime_error, "Simple well init failed.");
}
well_controls_set_current( simple_wells->ctrls[0] , 0);
well_controls_set_current( simple_wells->ctrls[1] , 0);
}
// Linear solver.
LinearSolverFactory linsolver(param);
@ -208,63 +151,47 @@ try
std::cout << "\n\n================ Starting main simulation loop ===============\n";
SimulatorReport rep;
if (!use_deck) {
// Simple simulation without a deck.
WellsManager wells(simple_wells);
// 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);
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of report step.
std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------"
<< "\n (number of steps left: "
<< timeMap->numTimesteps() - reportStepIdx << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
}
simtimer.setCurrentStepNum(reportStepIdx);
// Create and run simulator.
SimulatorCompressibleAd simulator(param,
*grid->c_grid(),
geology,
*props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
grav);
SimulatorTimer simtimer;
simtimer.init(param);
warnIfUnusedParams(param);
WellState well_state;
well_state.init(simple_wells, state);
rep = simulator.run(simtimer, state, well_state);
} else {
// With a deck, we may have more epochs etc.
WellState well_state;
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
SimulatorTimer simtimer;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of report step.
std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------"
<< "\n (number of steps left: "
<< timeMap->numTimesteps() - reportStepIdx << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every report step,
// since number of wells may change etc.
if (reportStepIdx == 0) {
well_state.init(wells.c_wells(), state);
}
simtimer.setCurrentStepNum(reportStepIdx);
// Create and run simulator.
SimulatorCompressibleAd simulator(param,
*grid->c_grid(),
*props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
grav);
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}
// Update total timing report and remember step number.
rep += epoch_rep;
if (reportStepIdx == 0) {
warnIfUnusedParams(param);
}
SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
if (output) {
epoch_rep.reportParam(epoch_os);
}
// Update total timing report and remember step number.
rep += epoch_rep;
}
std::cout << "\n\n================ End of simulation ===============\n\n";

3
examples/sim_fibo_ad.cpp
View File

@ -193,6 +193,8 @@ try
// initialize variables
simtimer.init(timeMap);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState);
SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of a report step.
@ -224,6 +226,7 @@ try
// Create and run simulator.
SimulatorFullyImplicitBlackoil<UnstructuredGrid> simulator(param,
*grid->c_grid(),
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,

5
examples/sim_fibo_ad_cp.cpp
View File

@ -230,12 +230,14 @@ try
<< std::flush;
WellStateFullyImplicitBlackoil well_state;
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
simtimer.init(timeMap);
Opm::DerivedGeology geology(*grid, *new_props, eclipseState);
SimulatorReport fullReport;
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
// Report on start of a report step.
@ -271,6 +273,7 @@ try
SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
*grid,
geology,
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,

110
examples/test_impestpfa_ad.cpp
View File

@ -1,110 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
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>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/ImpesTPFAAD.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/wells.h>
#include <algorithm>
#include <iostream>
int
main(int argc, char* argv[])
try
{
const Opm::parameter::ParameterGroup param(argc, argv, false);
const Opm::GridManager gm(5, 5);
const UnstructuredGrid* g = gm.c_grid();
const int nc = g->number_of_cells;
const Opm::BlackoilPropertiesBasic oldprops(param, 2, nc);
const Opm::BlackoilPropsAd props(oldprops);
Wells* wells = create_wells(2, 2, 5);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int num_inj = 3;
const int inj_cells[num_inj] = { 0, 1, 2 };
const int num_prod = 2;
const int prod_cells[num_prod] = { 20, 21 };
const double WI[3] = { 1e-12, 1e-12, 1e-12 };
bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells);
ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells);
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells);
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells);
if (!ok) {
OPM_THROW(std::runtime_error, "Something went wrong with well init.");
}
set_current_control(0, 0, wells);
set_current_control(1, 0, wells);
double grav[] = { /*1.0*/ 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav);
Opm::LinearSolverFactory linsolver(param);
Opm::ImpesTPFAAD ps(*g, props, geo, *wells, linsolver);
Opm::BlackoilState state;
initStateBasic(*g, oldprops, param, 0.0, state);
initBlackoilSurfvol(*g, oldprops, state);
Opm::WellState well_state;
well_state.init(wells, state);
ps.solve(1.0, state, well_state);
std::cout << "Cell pressure:" << std::endl;
std::copy(state.pressure().begin(), state.pressure().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Face flux:" << std::endl;
std::copy(state.faceflux().begin(), state.faceflux().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::cout << "Well bhp pressure:" << std::endl;
std::copy(well_state.bhp().begin(), well_state.bhp().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
return 0;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

127
examples/test_implicit_ad.cpp
View File

@ -1,127 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
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>
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <memory>
#include <algorithm>
#include <iostream>
namespace {
std::shared_ptr<Wells>
createWellConfig()
{
std::shared_ptr<Wells> wells(create_wells(2, 2, 2),
destroy_wells);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int num_inj = 1;
const int inj_cells[num_inj] = { 0 };
const int num_prod = 1;
const int prod_cells[num_prod] = { 19 };
const double WI[3] = { 1e-12, 1e-12, 1e-12 };
bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells.get());
ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells.get());
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells.get());
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells.get());
if (!ok) {
OPM_THROW(std::runtime_error, "Something went wrong with well init.");
}
set_current_control(0, 0, wells.get());
set_current_control(1, 0, wells.get());
return wells;
}
template <class Ostream, typename T, class A>
Ostream&
operator<<(Ostream& os, const std::vector<T,A>& v)
{
std::copy(v.begin(), v.end(), std::ostream_iterator<T>(os, " "));
return os;
}
}
int
main(int argc, char* argv[])
try
{
const Opm::parameter::ParameterGroup param(argc, argv, false);
const Opm::GridManager gm(20, 1);
const UnstructuredGrid* g = gm.c_grid();
const int nc = g->number_of_cells;
const Opm::BlackoilPropertiesBasic props0(param, 2, nc);
const Opm::BlackoilPropsAd props(props0);
std::shared_ptr<Wells> wells = createWellConfig();
double grav[] = { 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav);
Opm::NewtonIterationBlackoilSimple fis_solver(param);
Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver, /*hasDisgas*/ true, /*hasVapoil=*/false);
Opm::BlackoilState state;
initStateBasic(*g, props0, param, 0.0, state);
initBlackoilSurfvol(*g, props0, state);
Opm::WellStateFullyImplicitBlackoil well_state;
well_state.init(wells.get(), state);
solver.step(1.0, state, well_state);
std::cout << state.pressure() << '\n'
<< well_state.bhp() << '\n';
return 0;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

170
opm/autodiff/GeoProps.hpp
View File

@ -25,12 +25,20 @@
//#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/tpfa/TransTpfa.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include "disable_warning_pragmas.h"
#include <Eigen/Eigen>
#include "reenable_warning_pragmas.h"
#ifdef HAVE_DUNE_CORNERPOINT
#include <dune/common/version.hh>
#include <dune/grid/CpGrid.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#endif
#include <cstddef>
namespace Opm
@ -49,49 +57,85 @@ namespace Opm
/// Construct contained derived geological properties
/// from grid and property information.
template <class Props, class Grid>
DerivedGeology(const Grid& grid,
const Props& props ,
const double* grav = 0)
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))
{
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<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
const typename Vector::Index nc = numCells(grid);
std::transform(beginCellVolumes(grid), endCellVolumes(grid),
props.porosity(), pvol_.data(),
std::multiplies<double>());
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*>(& grid);
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
tpfa_trans_compute (ug, htrans.data() , trans_.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<nc; ++c){
z_[c] = cellCentroid(grid, c)[2];
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 = dimensions(grid);
typedef typename ADCell2FacesTraits<Grid>::Type Cell2Faces;
Cell2Faces c2f=cell2Faces(grid);
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 < 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]);
@ -109,15 +153,103 @@ namespace Opm
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
intersectionTransMult[faceIdx] *=
multipliers->getMultiplier(cartesianCellIdx, 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

9
opm/autodiff/SimulatorCompressibleAd.cpp
View File

@ -67,6 +67,7 @@ namespace Opm
public:
Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
@ -100,7 +101,7 @@ namespace Opm
const double* gravity_;
// Solvers
BlackoilPropsAd fluid_;
DerivedGeology geo_;
const DerivedGeology& geo_;
ImpesTPFAAD psolver_;
TransportSolverCompressibleTwophaseReorder tsolver_;
// Needed by column-based gravity segregation solver.
@ -114,13 +115,14 @@ namespace Opm
SimulatorCompressibleAd::SimulatorCompressibleAd(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity));
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, wells_manager, linsolver, gravity));
}
@ -232,6 +234,7 @@ namespace Opm
// \TODO: make CompressibleTpfa take bcs.
SimulatorCompressibleAd::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
@ -244,7 +247,7 @@ namespace Opm
wells_(wells_manager.c_wells()),
gravity_(gravity),
fluid_(props_),
geo_(grid_, fluid_, gravity_),
geo_(geo),
psolver_(grid_, fluid_, geo_, *wells_manager.c_wells(), linsolver),
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),

3
opm/autodiff/SimulatorCompressibleAd.hpp
View File

@ -31,6 +31,7 @@ namespace Opm
{
namespace parameter { class ParameterGroup; }
class BlackoilPropertiesInterface;
class DerivedGeology;
class RockCompressibility;
class WellsManager;
class LinearSolverInterface;
@ -60,6 +61,7 @@ namespace Opm
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] geo the "ready to use" geological properties of the reservoir
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] well_manager well manager
@ -67,6 +69,7 @@ namespace Opm
/// \param[in] gravity if non-null, gravity vector
SimulatorCompressibleAd(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,

2
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@ -32,6 +32,7 @@ namespace Opm
namespace parameter { class ParameterGroup; }
class BlackoilPropsAdInterface;
class RockCompressibility;
class DerivedGeology;
class WellsManager;
class NewtonIterationBlackoilInterface;
class SimulatorTimer;
@ -70,6 +71,7 @@ namespace Opm
/// \param[in] gravity if non-null, gravity vector
SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,

9
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -62,6 +62,7 @@ namespace Opm
public:
Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
@ -93,7 +94,7 @@ namespace Opm
const Wells* wells_;
const double* gravity_;
// Solvers
DerivedGeology geo_;
const DerivedGeology &geo_;
FullyImplicitBlackoilSolver<Grid> solver_;
// Misc. data
std::vector<int> allcells_;
@ -105,6 +106,7 @@ namespace Opm
template<class T>
SimulatorFullyImplicitBlackoil<T>::SimulatorFullyImplicitBlackoil(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
@ -114,7 +116,7 @@ namespace Opm
const bool has_vapoil )
{
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, has_disgas, has_vapoil));
pimpl_.reset(new Impl(param, grid, geo, props, rock_comp_props, wells_manager, linsolver, gravity, has_disgas, has_vapoil));
}
@ -192,6 +194,7 @@ namespace Opm
template<class T>
SimulatorFullyImplicitBlackoil<T>::Impl::Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
@ -205,7 +208,7 @@ namespace Opm
wells_manager_(wells_manager),
wells_(wells_manager.c_wells()),
gravity_(gravity),
geo_(grid_, props_, gravity_),
geo_(geo),
solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver, has_disgas, has_vapoil)
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),

351
tests/test_transmissibilitymultipliers.cpp Normal file
View File

@ -0,0 +1,351 @@
/*
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);
/////
/////
// 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);
/////
/////
// 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);
/////
/////
// 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);
/////
// 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);
/////
/////
// 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);
/////
/////
// 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);
/////
/////
// 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