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

Removes the dependency of FullyImpliciteBlackoilSolver onto UnstructuredGrid.

Browse Source 
With these changes it will be possible to use CpGrid with FIBOS except for the
output routines.
This commit is contained in:
Markus Blatt 2014-02-20 13:15:02 +01:00
parent 0a5262b7c3
commit 5112b8af26
9 changed files with 260 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
opm/autodiff/AutoDiffHelpers.hpp
View File

@ -21,6 +21,7 @@
#define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
@ -60,28 +61,17 @@ struct HelperOps
/// Constructs all helper vectors and matrices.
HelperOps(const UnstructuredGrid& grid)
{
const int nc = grid.number_of_cells;
const int nf = grid.number_of_faces;
using namespace AutoDiffGrid;
const int nc = numCells(grid);
const int nf = numFaces(grid);
// Define some neighbourhood-derived helper arrays.
typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
TwoColInt nb = Eigen::Map<TwoColInt>(grid.face_cells, nf, 2);
// std::cout << "nb = \n" << nb << std::endl;
TwoColBool nbib = nb >= 0;
OneColBool ifaces = nbib.rowwise().all();
const int num_internal = ifaces.cast<int>().sum();
// std::cout << num_internal << " internal faces." << std::endl;
TwoColInt nbi(num_internal, 2);
internal_faces.resize(num_internal);
int fi = 0;
for (int f = 0; f < nf; ++f) {
if (ifaces[f]) {
internal_faces[fi] = f;
nbi.row(fi) = nb.row(f);
++fi;
}
}
TwoColInt nbi;
extractInternalFaces(grid, internal_faces, nbi);
int num_internal=internal_faces.size();
// std::cout << "nbi = \n" << nbi << std::endl;
// Create matrices.
ngrad.resize(num_internal, nc);
@ -103,6 +93,7 @@ struct HelperOps
div = ngrad.transpose();
std::vector<Tri> fullngrad_tri;
fullngrad_tri.reserve(2*nf);
ADFaceCellTraits<UnstructuredGrid>::Type nb=faceCells(grid);
for (int i = 0; i < nf; ++i) {
if (nb(i,0) >= 0) {
fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
@ -131,8 +122,11 @@ struct HelperOps
const HelperOps& h,
const typename ADB::V& ifaceflux)
{
using namespace AutoDiffGrid;
typedef HelperOps::IFaces::Index IFIndex;
const IFIndex nif = h.internal_faces.size();
ADFaceCellTraits<UnstructuredGrid>::Type
face_cells = faceCells(g);
assert(nif == ifaceflux.size());
// Define selector structure.
@ -140,8 +134,8 @@ struct HelperOps
std::vector<Triplet> s; s.reserve(nif);
for (IFIndex iface = 0; iface < nif; ++iface) {
const int f = h.internal_faces[iface];
const int c1 = g.face_cells[2*f + 0];
const int c2 = g.face_cells[2*f + 1];
const int c1 = face_cells(f,0);
const int c2 = face_cells(f,1);
assert ((c1 >= 0) && (c2 >= 0));
@ -152,7 +146,7 @@ struct HelperOps
}
// Assemble explicit selector operator.
select_.resize(nif, g.number_of_cells);
select_.resize(nif, numCells(g));
select_.setFromTriplets(s.begin(), s.end());
}

80
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -51,9 +51,68 @@ namespace Opm
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const bool init_rock)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.cell_centroids, grid.dimensions, init_rock);
}
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const bool init_rock)
{
init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.cell_centroids, grid.dimensions, init_rock);
}
#ifdef HAVE_DUNE_CORNERPOINT
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
const Dune::CpGrid& grid,
const bool init_rock )
{
init(deck, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
static_cast<const int*>(&grid.logicalCartesianSize()[0]),
grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
}
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
const Dune::CpGrid& grid,
const bool init_rock )
{
init(newParserDeck, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
static_cast<const int*>(&grid.logicalCartesianSize()[0]),
grid.beginCellCentroids(), Dune::CpGrid::dimension, init_rock);
}
#endif
/// Constructor wrapping an opm-core black oil interface.
template<class T>
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimensions,
const bool init_rock)
{
init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids,
dimensions, init_rock);
}
template<class T>
void BlackoilPropsAdFromDeck::init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimensions,
const bool init_rock)
{
if (init_rock){
rock_.init(deck, grid);
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
const int samples = 0;
const int region_number = 0;
@ -122,7 +181,7 @@ namespace Opm
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, grid, -1);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions, -1);
if (phase_usage_.num_phases != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
@ -131,14 +190,17 @@ namespace Opm
}
}
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const bool init_rock)
template<class T>
void BlackoilPropsAdFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimensions,
const bool init_rock)
{
if (init_rock){
rock_.init(newParserDeck, grid);
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
}
const int samples = 0;
@ -218,7 +280,7 @@ namespace Opm
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, -1);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, dimensions, -1);
if (phase_usage_.num_phases != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "

57
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@ -33,6 +33,10 @@
#include <boost/shared_ptr.hpp>
#include <boost/scoped_ptr.hpp>
#ifdef HAVE_DUNE_CORNERPOINT
#include <dune/grid/CpGrid.hpp>
#endif
namespace Opm
{
@ -60,6 +64,41 @@ namespace Opm
const UnstructuredGrid& grid,
const bool init_rock = true );
#ifdef HAVE_DUNE_CORNERPOINT
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
const Dune::CpGrid& grid,
const bool init_rock = true );
/// Constructor wrapping an opm-core black oil interface.
BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
const Dune::CpGrid& grid,
const bool init_rock = true );
#endif
/// Constructor taking not a grid but only the needed information
template<class T>
BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimensions,
const bool init_rock);
/// Constructor taking not a grid but only the needed information
template<class T>
BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimensions,
const bool init_rock);
////////////////////////////
// Rock interface //
////////////////////////////
@ -329,6 +368,24 @@ namespace Opm
const std::vector<int>& cells);
private:
/// Initializes the properties.
template<class T>
void init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const bool init_rock);
/// Initializes the properties.
template<class T>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const bool init_rock);
RockFromDeck rock_;
boost::scoped_ptr<SaturationPropsInterface> satprops_;
PhaseUsage phase_usage_;

49
opm/autodiff/FullyImplicitBlackoilSolver.cpp
View File

@ -22,6 +22,7 @@
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp>
@ -79,15 +80,20 @@ namespace {
const HelperOps& ops ,
const GeoProps& geo )
{
const int nc = grid.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid);
SparseTableView c2f = cell2Faces(grid);
std::vector<int> f2hf(2 * grid.number_of_faces, -1);
std::vector<int> f2hf(2 * numFaces(grid), -1);
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
face_cells = faceCells(grid);
for (int c = 0, i = 0; c < nc; ++c) {
for (; i < grid.cell_facepos[c + 1]; ++i) {
const int f = grid.cell_faces[ i ];
const int p = 0 + (grid.face_cells[2*f + 0] != c);
typename SparseTableView::row_type faces=c2f[c];
typedef typename SparseTableView::row_type::iterator Iter;
for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
const int p = 0 + (face_cells(*f, 0) != c);
f2hf[2*f + p] = i;
f2hf[2*(*f) + p] = i;
}
}
@ -103,8 +109,8 @@ namespace {
std::vector<Tri> grav; grav.reserve(2 * ni);
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
const int f = ops.internal_faces[ i ];
const int c1 = grid.face_cells[2*f + 0];
const int c2 = grid.face_cells[2*f + 1];
const int c1 = faceCells(grid)(f, 0);
const int c2 = faceCells(grid)(f, 1);
assert ((c1 >= 0) && (c2 >= 0));
@ -125,9 +131,10 @@ namespace {
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
{
using namespace Opm::AutoDiffGrid;
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
const int dim = grid.dimensions;
const int dim = dimensions(grid);
V wdp = V::Zero(nperf,1);
assert(wdp.size() == rho.size());
@ -141,7 +148,7 @@ namespace {
const double ref_depth = wells.depth_ref[w];
for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
const int cell = wells.well_cells[j];
const double cell_depth = grid.cell_centroids[dim * cell + dim - 1];
const double cell_depth = cellCentroid(grid, cell)[dim - 1];
wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
}
}
@ -203,7 +210,7 @@ namespace {
, linsolver_ (linsolver)
, active_(activePhases(fluid.phaseUsage()))
, canph_ (active2Canonical(fluid.phaseUsage()))
, cells_ (buildAllCells(grid.number_of_cells))
, cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
, ops_ (grid)
, wops_ (wells)
, grav_ (gravityOperator(grid_, ops_, geo_))
@ -333,7 +340,8 @@ namespace {
FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
const WellState& xw)
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = x.numPhases();
// The block pattern assumes the following primary variables:
@ -429,7 +437,8 @@ namespace {
FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
const WellState& xw)
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = x.numPhases();
std::vector<V> vars0;
@ -618,6 +627,7 @@ namespace {
const BlackoilState& x ,
const WellState& xw )
{
using namespace Opm::AutoDiffGrid;
// Create the primary variables.
const SolutionState state = variableState(x, xw);
@ -680,7 +690,7 @@ namespace {
// Contribution to mass balance will have to wait.
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_);
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
@ -700,7 +710,7 @@ namespace {
// Compute well pressure differentials.
// Construct pressure difference vector for wells.
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const int dim = grid_.dimensions;
const int dim = dimensions(grid_);
const double* g = geo_.gravity();
if (g) {
// Guard against gravity in anything but last dimension.
@ -880,8 +890,9 @@ namespace {
BlackoilState& state,
WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
const int np = fluid_.numPhases();
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_);
const int nw = wells_.number_of_wells;
const V null;
assert(null.size() == 0);
@ -1128,7 +1139,8 @@ namespace {
std::vector<ADB>
FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state) const
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const std::vector<int>& bpat = state.pressure.blockPattern();
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
@ -1153,7 +1165,8 @@ namespace {
std::vector<ADB>
FullyImplicitBlackoilSolver::computePressures(const SolutionState& state) const
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const std::vector<int>& bpat = state.pressure.blockPattern();
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);

31
opm/autodiff/GeoProps.hpp
View File

@ -22,6 +22,7 @@
#include <opm/core/grid.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/autodiff/GridHelpers.hpp>
#include <Eigen/Eigen>
namespace Opm
@ -43,14 +44,15 @@ namespace Opm
DerivedGeology(const UnstructuredGrid& grid,
const Props& props ,
const double* grav = 0)
: pvol_ (grid.number_of_cells)
, trans_(grid.number_of_faces)
, gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
, z_(grid.number_of_cells)
: 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;
// Pore volume
const typename Vector::Index nc = grid.number_of_cells;
std::transform(grid.cell_volumes, grid.cell_volumes + nc,
const typename Vector::Index nc = numCells(grid);
std::transform(beginCellVolumes(grid), endCellVolumes(grid),
props.porosity(), pvol_.data(),
std::multiplies<double>());
@ -62,26 +64,27 @@ namespace Opm
// Compute z coordinates
for (int c = 0; c<nc; ++c){
z_[c] = grid.cell_centroids[c*3 + 2];
z_[c] = cellCentroid(grid, c)[2];
}
// Gravity potential
std::fill(gravity_, gravity_ + 3, 0.0);
if (grav != 0) {
const typename Vector::Index nd = grid.dimensions;
const typename Vector::Index nd = dimensions(grid);
SparseTableView c2f=cell2Faces(grid);
for (typename Vector::Index c = 0; c < nc; ++c) {
const double* const cc = & grid.cell_centroids[c*nd + 0];
const double* const cc = cellCentroid(grid, c);
const int* const p = grid.cell_facepos;
for (int i = p[c]; i < p[c + 1]; ++i) {
const int f = grid.cell_faces[i];
typename SparseTableView::row_type faces=c2f[c];
typedef SparseTableView::row_type::iterator Iter;
const double* const fc = & grid.face_centroids[f*nd + 0];
for (Iter f=faces.begin(), end=faces.end(); f!=end; ++f) {
const double* const fc = faceCentroid(grid, *f);
for (typename Vector::Index d = 0; d < nd; ++d) {
gpot_[i] += grav[d] * (fc[d] - cc[d]);
gpot_[f-faces.begin()] += grav[d] * (fc[d] - cc[d]);
}
}
}

54
opm/autodiff/ImpesTPFAAD.cpp
View File

@ -19,6 +19,7 @@
#include <opm/autodiff/ImpesTPFAAD.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
@ -54,15 +55,21 @@ namespace {
const HelperOps& ops ,
const GeoProps& geo )
{
const int nc = grid.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid);
std::vector<int> f2hf(2 * grid.number_of_faces, -1);
for (int c = 0, i = 0; c < nc; ++c) {
for (; i < grid.cell_facepos[c + 1]; ++i) {
const int f = grid.cell_faces[ i ];
const int p = 0 + (grid.face_cells[2*f + 0] != c);
f2hf[2*f + p] = i;
std::vector<int> f2hf(2 * numFaces(grid), -1);
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
face_cells;
SparseTableView c2f=cell2Faces(grid);
for (int c = 0; c < nc; ++c) {
typename SparseTableView::row_type
cell_faces = c2f[c];
typedef typename SparseTableView::row_type::iterator Iter;
for (Iter f=cell_faces.begin(), end=cell_faces.end();
f!=end; ++end) {
const int p = 0 + (face_cells(*f,0) != c);
f2hf[2*(*f) + p] = f-c2f[0].begin();
}
}
@ -78,8 +85,10 @@ namespace {
std::vector<Tri> grav; grav.reserve(2 * ni);
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
const int f = ops.internal_faces[ i ];
const int c1 = grid.face_cells[2*f + 0];
const int c2 = grid.face_cells[2*f + 1];
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
face_cells=faceCells(grid);
const int c1 = face_cells(f,0);
const int c2 = face_cells(f,1);
assert ((c1 >= 0) && (c2 >= 0));
@ -98,9 +107,10 @@ namespace {
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
{
using namespace Opm::AutoDiffGrid;
const int nw = wells.number_of_wells;
const int nperf = wells.well_connpos[nw];
const int dim = grid.dimensions;
const int dim = dimensions(grid);
V wdp = V::Zero(nperf,1);
assert(wdp.size() == rho.size());
@ -157,7 +167,8 @@ namespace {
BlackoilState& state,
WellState& well_state)
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = state.numPhases();
well_flow_residual_.resize(np, ADB::null());
@ -226,15 +237,16 @@ namespace {
const BlackoilState& state,
const WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
// Suppress warnings about "unused parameters".
static_cast<void>(dt);
static_cast<void>(well_state);
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_);
const int np = state.numPhases();
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
const int dim = grid_.dimensions;
const int dim = dimensions(grid_);
const std::vector<int> cells = buildAllCells(nc);
@ -284,9 +296,9 @@ namespace {
const BlackoilState& state,
const WellState& well_state)
{
using namespace Opm::AutoDiffGrid;
const V& pv = geo_.poreVolume();
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_); ;
const int np = state.numPhases();
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
@ -415,7 +427,8 @@ namespace {
ImpesTPFAAD::solveJacobianSystem(BlackoilState& state,
WellState& well_state) const
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int nw = wells_.number_of_wells;
// const int np = state.numPhases();
@ -458,9 +471,10 @@ namespace {
ImpesTPFAAD::computeFluxes(BlackoilState& state,
WellState& well_state) const
{
using namespace Opm::AutoDiffGrid;
// This method computes state.faceflux(),
// well_state.perfRates() and well_state.perfPress().
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_);
const int np = state.numPhases();
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
@ -515,8 +529,8 @@ namespace {
flux += face_mob * head;
}
V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces);
std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().begin());
V all_flux = superset(flux, ops_.internal_faces, numFaces(grid_));
std::copy(all_flux.data(), all_flux.data() + numFaces(grid_), state.faceflux().begin());
perf_flux = -perf_flux; // well_state.perfRates() assumed to be inflows.
std::copy(perf_flux.data(), perf_flux.data() + nperf, well_state.perfRates().begin());

20
opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
View File

@ -158,7 +158,14 @@ namespace Opm
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
AutoDiffGrid::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(grid),
state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
}
@ -174,7 +181,14 @@ namespace Opm
dm["pressure"] = &state.pressure();
dm["surfvolume"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
AutoDiffGrid::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(grid),
state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
// Write data (not grid) in Matlab format
@ -300,7 +314,7 @@ namespace Opm
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Misc init.
const int num_cells = grid.number_of_cells;
const int num_cells = AutoDiffGrid::numCells(grid);
allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells_[cell] = cell;

13
opm/autodiff/SimulatorIncompTwophaseAd.cpp
View File

@ -23,6 +23,7 @@
#endif // HAVE_CONFIG_H
#include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@ -400,7 +401,7 @@ namespace Opm
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
// Misc init.
const int num_cells = grid.number_of_cells;
const int num_cells = Opm::AutoDiffGrid::numCells(grid);
allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells_[cell] = cell;
@ -498,11 +499,13 @@ namespace Opm
double av_prev_press = 0.0;
double av_press = 0.0;
double tot_vol = 0.0;
const int num_cells = grid_.number_of_cells;
const int num_cells = Opm::AutoDiffGrid::numCells(grid_);
for (int cell = 0; cell < num_cells; ++cell) {
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
tot_vol += grid_.cell_volumes[cell];
av_prev_press += initial_pressure[cell]*
Opm::AutoDiffGrid::cellVolume(grid_, cell);
av_press += state.pressure()[cell]*
Opm::AutoDiffGrid::cellVolume(grid_, cell);
tot_vol += Opm::AutoDiffGrid::cellVolume(grid_, cell);
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;

18
opm/autodiff/TransportSolverTwophaseAd.cpp
View File

@ -51,21 +51,22 @@ namespace Opm
tol_(param.getDefault("nl_tolerance", 1e-9)),
maxit_(param.getDefault("nl_maxiter", 30))
{
const int nc = grid_.number_of_cells;
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
allcells_.resize(nc);
for (int i = 0; i < nc; ++i) {
allcells_[i] = i;
}
if (gravity && gravity[grid_.dimensions - 1] != 0.0) {
gravity_ = gravity[grid_.dimensions - 1];
for (int dd = 0; dd < grid_.dimensions - 1; ++dd) {
if (gravity && gravity[dimensions(grid_) - 1] != 0.0) {
gravity_ = gravity[dimensions(grid_) - 1];
for (int dd = 0; dd < dimensions(grid_) - 1; ++dd) {
if (gravity[dd] != 0.0) {
OPM_THROW(std::runtime_error, "TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
}
}
V htrans(grid.cell_facepos[grid.number_of_cells]);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
V trans(grid_.number_of_faces);
V trans(numFaces(grid_));
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans.data());
transi_ = subset(trans, ops_.internal_faces);
}
@ -161,16 +162,17 @@ namespace Opm
const double dt,
TwophaseState& state)
{
using namespace Opm::AutoDiffGrid;
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
typedef Eigen::Map<const V> Vec;
const int nc = grid_.number_of_cells;
const int nc = numCells(grid_);
const TwoCol s0 = Eigen::Map<const TwoCol>(state.saturation().data(), nc, 2);
double res_norm = 1e100;
const V sw0 = s0.leftCols<1>();
// sw1 is the object that will be changed every Newton iteration.
// V sw1 = sw0;
V sw1 = 0.5*V::Ones(nc,1);
const V dflux_all = Vec(state.faceflux().data(), grid_.number_of_faces, 1);
const V dflux_all = Vec(state.faceflux().data(), numFaces(grid_), 1);
const int num_internal = ops_.internal_faces.size();
V dflux = subset(dflux_all, ops_.internal_faces);
@ -184,7 +186,7 @@ namespace Opm
const V p1 = Vec(state.pressure().data(), nc, 1);
const V ndp = (ops_.ngrad * p1.matrix()).array();
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DynArr;
const V z = Eigen::Map<DynArr>(grid_.cell_centroids, nc, grid_.dimensions).rightCols<1>();
const V z = cellCentroidsZ(grid_);
const V ndz = (ops_.ngrad * z.matrix()).array();
assert(num_internal == ndp.size());
const double* density = props_.density();