mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Removes the dependency of FullyImpliciteBlackoilSolver onto UnstructuredGrid.
With these changes it will be possible to use CpGrid with FIBOS except for the output routines.
This commit is contained in:
parent
0a5262b7c3
commit
5112b8af26
@ -21,6 +21,7 @@
|
|||||||
#define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
|
#define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
|
||||||
|
|
||||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||||
|
#include <opm/autodiff/GridHelpers.hpp>
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
|
|
||||||
@ -60,28 +61,17 @@ struct HelperOps
|
|||||||
/// Constructs all helper vectors and matrices.
|
/// Constructs all helper vectors and matrices.
|
||||||
HelperOps(const UnstructuredGrid& grid)
|
HelperOps(const UnstructuredGrid& grid)
|
||||||
{
|
{
|
||||||
const int nc = grid.number_of_cells;
|
using namespace AutoDiffGrid;
|
||||||
const int nf = grid.number_of_faces;
|
const int nc = numCells(grid);
|
||||||
|
const int nf = numFaces(grid);
|
||||||
// Define some neighbourhood-derived helper arrays.
|
// Define some neighbourhood-derived helper arrays.
|
||||||
typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
|
typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
|
||||||
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
|
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
|
||||||
typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
|
typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
|
||||||
TwoColInt nb = Eigen::Map<TwoColInt>(grid.face_cells, nf, 2);
|
TwoColInt nbi;
|
||||||
// std::cout << "nb = \n" << nb << std::endl;
|
extractInternalFaces(grid, internal_faces, nbi);
|
||||||
TwoColBool nbib = nb >= 0;
|
int num_internal=internal_faces.size();
|
||||||
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;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// std::cout << "nbi = \n" << nbi << std::endl;
|
// std::cout << "nbi = \n" << nbi << std::endl;
|
||||||
// Create matrices.
|
// Create matrices.
|
||||||
ngrad.resize(num_internal, nc);
|
ngrad.resize(num_internal, nc);
|
||||||
@ -103,6 +93,7 @@ struct HelperOps
|
|||||||
div = ngrad.transpose();
|
div = ngrad.transpose();
|
||||||
std::vector<Tri> fullngrad_tri;
|
std::vector<Tri> fullngrad_tri;
|
||||||
fullngrad_tri.reserve(2*nf);
|
fullngrad_tri.reserve(2*nf);
|
||||||
|
ADFaceCellTraits<UnstructuredGrid>::Type nb=faceCells(grid);
|
||||||
for (int i = 0; i < nf; ++i) {
|
for (int i = 0; i < nf; ++i) {
|
||||||
if (nb(i,0) >= 0) {
|
if (nb(i,0) >= 0) {
|
||||||
fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
|
fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
|
||||||
@ -131,8 +122,11 @@ struct HelperOps
|
|||||||
const HelperOps& h,
|
const HelperOps& h,
|
||||||
const typename ADB::V& ifaceflux)
|
const typename ADB::V& ifaceflux)
|
||||||
{
|
{
|
||||||
|
using namespace AutoDiffGrid;
|
||||||
typedef HelperOps::IFaces::Index IFIndex;
|
typedef HelperOps::IFaces::Index IFIndex;
|
||||||
const IFIndex nif = h.internal_faces.size();
|
const IFIndex nif = h.internal_faces.size();
|
||||||
|
ADFaceCellTraits<UnstructuredGrid>::Type
|
||||||
|
face_cells = faceCells(g);
|
||||||
assert(nif == ifaceflux.size());
|
assert(nif == ifaceflux.size());
|
||||||
|
|
||||||
// Define selector structure.
|
// Define selector structure.
|
||||||
@ -140,8 +134,8 @@ struct HelperOps
|
|||||||
std::vector<Triplet> s; s.reserve(nif);
|
std::vector<Triplet> s; s.reserve(nif);
|
||||||
for (IFIndex iface = 0; iface < nif; ++iface) {
|
for (IFIndex iface = 0; iface < nif; ++iface) {
|
||||||
const int f = h.internal_faces[iface];
|
const int f = h.internal_faces[iface];
|
||||||
const int c1 = g.face_cells[2*f + 0];
|
const int c1 = face_cells(f,0);
|
||||||
const int c2 = g.face_cells[2*f + 1];
|
const int c2 = face_cells(f,1);
|
||||||
|
|
||||||
assert ((c1 >= 0) && (c2 >= 0));
|
assert ((c1 >= 0) && (c2 >= 0));
|
||||||
|
|
||||||
@ -152,7 +146,7 @@ struct HelperOps
|
|||||||
}
|
}
|
||||||
|
|
||||||
// Assemble explicit selector operator.
|
// Assemble explicit selector operator.
|
||||||
select_.resize(nif, g.number_of_cells);
|
select_.resize(nif, numCells(g));
|
||||||
select_.setFromTriplets(s.begin(), s.end());
|
select_.setFromTriplets(s.begin(), s.end());
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -51,9 +51,68 @@ namespace Opm
|
|||||||
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const EclipseGridParser& deck,
|
||||||
const UnstructuredGrid& grid,
|
const UnstructuredGrid& grid,
|
||||||
const bool init_rock)
|
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){
|
if (init_rock){
|
||||||
rock_.init(deck, grid);
|
rock_.init(deck, number_of_cells, global_cell, cart_dims);
|
||||||
}
|
}
|
||||||
const int samples = 0;
|
const int samples = 0;
|
||||||
const int region_number = 0;
|
const int region_number = 0;
|
||||||
@ -122,7 +181,7 @@ namespace Opm
|
|||||||
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
||||||
satprops_.reset(ptr);
|
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()) {
|
if (phase_usage_.num_phases != satprops_->numPhases()) {
|
||||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
||||||
@ -131,14 +190,17 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template<class T>
|
||||||
/// Constructor wrapping an opm-core black oil interface.
|
void BlackoilPropsAdFromDeck::init(Opm::DeckConstPtr newParserDeck,
|
||||||
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(Opm::DeckConstPtr newParserDeck,
|
int number_of_cells,
|
||||||
const UnstructuredGrid& grid,
|
const int* global_cell,
|
||||||
|
const int* cart_dims,
|
||||||
|
T begin_cell_centroids,
|
||||||
|
int dimensions,
|
||||||
const bool init_rock)
|
const bool init_rock)
|
||||||
{
|
{
|
||||||
if (init_rock){
|
if (init_rock){
|
||||||
rock_.init(newParserDeck, grid);
|
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
|
||||||
}
|
}
|
||||||
|
|
||||||
const int samples = 0;
|
const int samples = 0;
|
||||||
@ -218,7 +280,7 @@ namespace Opm
|
|||||||
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
|
||||||
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
|
||||||
satprops_.reset(ptr);
|
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()) {
|
if (phase_usage_.num_phases != satprops_->numPhases()) {
|
||||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
||||||
|
@ -33,6 +33,10 @@
|
|||||||
#include <boost/shared_ptr.hpp>
|
#include <boost/shared_ptr.hpp>
|
||||||
#include <boost/scoped_ptr.hpp>
|
#include <boost/scoped_ptr.hpp>
|
||||||
|
|
||||||
|
#ifdef HAVE_DUNE_CORNERPOINT
|
||||||
|
#include <dune/grid/CpGrid.hpp>
|
||||||
|
#endif
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
|
|
||||||
@ -60,6 +64,41 @@ namespace Opm
|
|||||||
const UnstructuredGrid& grid,
|
const UnstructuredGrid& grid,
|
||||||
const bool init_rock = true );
|
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 //
|
// Rock interface //
|
||||||
////////////////////////////
|
////////////////////////////
|
||||||
@ -329,6 +368,24 @@ namespace Opm
|
|||||||
const std::vector<int>& cells);
|
const std::vector<int>& cells);
|
||||||
|
|
||||||
private:
|
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_;
|
RockFromDeck rock_;
|
||||||
boost::scoped_ptr<SaturationPropsInterface> satprops_;
|
boost::scoped_ptr<SaturationPropsInterface> satprops_;
|
||||||
PhaseUsage phase_usage_;
|
PhaseUsage phase_usage_;
|
||||||
|
@ -22,6 +22,7 @@
|
|||||||
|
|
||||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||||
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
||||||
|
#include <opm/autodiff/GridHelpers.hpp>
|
||||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||||
#include <opm/autodiff/GeoProps.hpp>
|
#include <opm/autodiff/GeoProps.hpp>
|
||||||
|
|
||||||
@ -79,15 +80,20 @@ namespace {
|
|||||||
const HelperOps& ops ,
|
const HelperOps& ops ,
|
||||||
const GeoProps& geo )
|
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 (int c = 0, i = 0; c < nc; ++c) {
|
||||||
for (; i < grid.cell_facepos[c + 1]; ++i) {
|
typename SparseTableView::row_type faces=c2f[c];
|
||||||
const int f = grid.cell_faces[ i ];
|
typedef typename SparseTableView::row_type::iterator Iter;
|
||||||
const int p = 0 + (grid.face_cells[2*f + 0] != c);
|
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);
|
std::vector<Tri> grav; grav.reserve(2 * ni);
|
||||||
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
|
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
|
||||||
const int f = ops.internal_faces[ i ];
|
const int f = ops.internal_faces[ i ];
|
||||||
const int c1 = grid.face_cells[2*f + 0];
|
const int c1 = faceCells(grid)(f, 0);
|
||||||
const int c2 = grid.face_cells[2*f + 1];
|
const int c2 = faceCells(grid)(f, 1);
|
||||||
|
|
||||||
assert ((c1 >= 0) && (c2 >= 0));
|
assert ((c1 >= 0) && (c2 >= 0));
|
||||||
|
|
||||||
@ -125,9 +131,10 @@ namespace {
|
|||||||
|
|
||||||
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
|
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 nw = wells.number_of_wells;
|
||||||
const int nperf = wells.well_connpos[nw];
|
const int nperf = wells.well_connpos[nw];
|
||||||
const int dim = grid.dimensions;
|
const int dim = dimensions(grid);
|
||||||
V wdp = V::Zero(nperf,1);
|
V wdp = V::Zero(nperf,1);
|
||||||
assert(wdp.size() == rho.size());
|
assert(wdp.size() == rho.size());
|
||||||
|
|
||||||
@ -141,7 +148,7 @@ namespace {
|
|||||||
const double ref_depth = wells.depth_ref[w];
|
const double ref_depth = wells.depth_ref[w];
|
||||||
for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
|
for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
|
||||||
const int cell = wells.well_cells[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);
|
wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -203,7 +210,7 @@ namespace {
|
|||||||
, linsolver_ (linsolver)
|
, linsolver_ (linsolver)
|
||||||
, active_(activePhases(fluid.phaseUsage()))
|
, active_(activePhases(fluid.phaseUsage()))
|
||||||
, canph_ (active2Canonical(fluid.phaseUsage()))
|
, canph_ (active2Canonical(fluid.phaseUsage()))
|
||||||
, cells_ (buildAllCells(grid.number_of_cells))
|
, cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
|
||||||
, ops_ (grid)
|
, ops_ (grid)
|
||||||
, wops_ (wells)
|
, wops_ (wells)
|
||||||
, grav_ (gravityOperator(grid_, ops_, geo_))
|
, grav_ (gravityOperator(grid_, ops_, geo_))
|
||||||
@ -333,7 +340,8 @@ namespace {
|
|||||||
FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
|
FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
|
||||||
const WellState& xw)
|
const WellState& xw)
|
||||||
{
|
{
|
||||||
const int nc = grid_.number_of_cells;
|
using namespace Opm::AutoDiffGrid;
|
||||||
|
const int nc = numCells(grid_);
|
||||||
const int np = x.numPhases();
|
const int np = x.numPhases();
|
||||||
|
|
||||||
// The block pattern assumes the following primary variables:
|
// The block pattern assumes the following primary variables:
|
||||||
@ -429,7 +437,8 @@ namespace {
|
|||||||
FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
|
FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
|
||||||
const WellState& xw)
|
const WellState& xw)
|
||||||
{
|
{
|
||||||
const int nc = grid_.number_of_cells;
|
using namespace Opm::AutoDiffGrid;
|
||||||
|
const int nc = numCells(grid_);
|
||||||
const int np = x.numPhases();
|
const int np = x.numPhases();
|
||||||
|
|
||||||
std::vector<V> vars0;
|
std::vector<V> vars0;
|
||||||
@ -618,6 +627,7 @@ namespace {
|
|||||||
const BlackoilState& x ,
|
const BlackoilState& x ,
|
||||||
const WellState& xw )
|
const WellState& xw )
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
// Create the primary variables.
|
// Create the primary variables.
|
||||||
const SolutionState state = variableState(x, xw);
|
const SolutionState state = variableState(x, xw);
|
||||||
|
|
||||||
@ -680,7 +690,7 @@ namespace {
|
|||||||
|
|
||||||
// Contribution to mass balance will have to wait.
|
// 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 np = wells_.number_of_phases;
|
||||||
const int nw = wells_.number_of_wells;
|
const int nw = wells_.number_of_wells;
|
||||||
const int nperf = wells_.well_connpos[nw];
|
const int nperf = wells_.well_connpos[nw];
|
||||||
@ -700,7 +710,7 @@ namespace {
|
|||||||
// Compute well pressure differentials.
|
// Compute well pressure differentials.
|
||||||
// Construct pressure difference vector for wells.
|
// Construct pressure difference vector for wells.
|
||||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||||
const int dim = grid_.dimensions;
|
const int dim = dimensions(grid_);
|
||||||
const double* g = geo_.gravity();
|
const double* g = geo_.gravity();
|
||||||
if (g) {
|
if (g) {
|
||||||
// Guard against gravity in anything but last dimension.
|
// Guard against gravity in anything but last dimension.
|
||||||
@ -880,8 +890,9 @@ namespace {
|
|||||||
BlackoilState& state,
|
BlackoilState& state,
|
||||||
WellState& well_state)
|
WellState& well_state)
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
const int np = fluid_.numPhases();
|
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 int nw = wells_.number_of_wells;
|
||||||
const V null;
|
const V null;
|
||||||
assert(null.size() == 0);
|
assert(null.size() == 0);
|
||||||
@ -1128,7 +1139,8 @@ namespace {
|
|||||||
std::vector<ADB>
|
std::vector<ADB>
|
||||||
FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state) const
|
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 std::vector<int>& bpat = state.pressure.blockPattern();
|
||||||
|
|
||||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||||
@ -1153,7 +1165,8 @@ namespace {
|
|||||||
std::vector<ADB>
|
std::vector<ADB>
|
||||||
FullyImplicitBlackoilSolver::computePressures(const SolutionState& state) const
|
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 std::vector<int>& bpat = state.pressure.blockPattern();
|
||||||
|
|
||||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||||
|
@ -22,6 +22,7 @@
|
|||||||
|
|
||||||
#include <opm/core/grid.h>
|
#include <opm/core/grid.h>
|
||||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||||
|
#include <opm/autodiff/GridHelpers.hpp>
|
||||||
#include <Eigen/Eigen>
|
#include <Eigen/Eigen>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
@ -43,14 +44,15 @@ namespace Opm
|
|||||||
DerivedGeology(const UnstructuredGrid& grid,
|
DerivedGeology(const UnstructuredGrid& grid,
|
||||||
const Props& props ,
|
const Props& props ,
|
||||||
const double* grav = 0)
|
const double* grav = 0)
|
||||||
: pvol_ (grid.number_of_cells)
|
: pvol_ (Opm::AutoDiffGrid::numCells(grid))
|
||||||
, trans_(grid.number_of_faces)
|
, trans_(Opm::AutoDiffGrid::numFaces(grid))
|
||||||
, gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
|
, gpot_ (Vector::Zero(Opm::AutoDiffGrid::cell2Faces(grid).noEntries(), 1))
|
||||||
, z_(grid.number_of_cells)
|
, z_(Opm::AutoDiffGrid::numCells(grid))
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
// Pore volume
|
// Pore volume
|
||||||
const typename Vector::Index nc = grid.number_of_cells;
|
const typename Vector::Index nc = numCells(grid);
|
||||||
std::transform(grid.cell_volumes, grid.cell_volumes + nc,
|
std::transform(beginCellVolumes(grid), endCellVolumes(grid),
|
||||||
props.porosity(), pvol_.data(),
|
props.porosity(), pvol_.data(),
|
||||||
std::multiplies<double>());
|
std::multiplies<double>());
|
||||||
|
|
||||||
@ -62,26 +64,27 @@ namespace Opm
|
|||||||
|
|
||||||
// Compute z coordinates
|
// Compute z coordinates
|
||||||
for (int c = 0; c<nc; ++c){
|
for (int c = 0; c<nc; ++c){
|
||||||
z_[c] = grid.cell_centroids[c*3 + 2];
|
z_[c] = cellCentroid(grid, c)[2];
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
// Gravity potential
|
// Gravity potential
|
||||||
std::fill(gravity_, gravity_ + 3, 0.0);
|
std::fill(gravity_, gravity_ + 3, 0.0);
|
||||||
if (grav != 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) {
|
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;
|
typename SparseTableView::row_type faces=c2f[c];
|
||||||
for (int i = p[c]; i < p[c + 1]; ++i) {
|
typedef SparseTableView::row_type::iterator Iter;
|
||||||
const int f = grid.cell_faces[i];
|
|
||||||
|
|
||||||
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) {
|
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]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -19,6 +19,7 @@
|
|||||||
|
|
||||||
#include <opm/autodiff/ImpesTPFAAD.hpp>
|
#include <opm/autodiff/ImpesTPFAAD.hpp>
|
||||||
#include <opm/autodiff/GeoProps.hpp>
|
#include <opm/autodiff/GeoProps.hpp>
|
||||||
|
#include <opm/autodiff/GridHelpers.hpp>
|
||||||
|
|
||||||
#include <opm/core/simulator/BlackoilState.hpp>
|
#include <opm/core/simulator/BlackoilState.hpp>
|
||||||
#include <opm/core/simulator/WellState.hpp>
|
#include <opm/core/simulator/WellState.hpp>
|
||||||
@ -54,15 +55,21 @@ namespace {
|
|||||||
const HelperOps& ops ,
|
const HelperOps& ops ,
|
||||||
const GeoProps& geo )
|
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);
|
std::vector<int> f2hf(2 * numFaces(grid), -1);
|
||||||
for (int c = 0, i = 0; c < nc; ++c) {
|
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
|
||||||
for (; i < grid.cell_facepos[c + 1]; ++i) {
|
face_cells;
|
||||||
const int f = grid.cell_faces[ i ];
|
SparseTableView c2f=cell2Faces(grid);
|
||||||
const int p = 0 + (grid.face_cells[2*f + 0] != c);
|
for (int c = 0; c < nc; ++c) {
|
||||||
|
typename SparseTableView::row_type
|
||||||
f2hf[2*f + p] = i;
|
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);
|
std::vector<Tri> grav; grav.reserve(2 * ni);
|
||||||
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
|
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
|
||||||
const int f = ops.internal_faces[ i ];
|
const int f = ops.internal_faces[ i ];
|
||||||
const int c1 = grid.face_cells[2*f + 0];
|
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
|
||||||
const int c2 = grid.face_cells[2*f + 1];
|
face_cells=faceCells(grid);
|
||||||
|
const int c1 = face_cells(f,0);
|
||||||
|
const int c2 = face_cells(f,1);
|
||||||
|
|
||||||
assert ((c1 >= 0) && (c2 >= 0));
|
assert ((c1 >= 0) && (c2 >= 0));
|
||||||
|
|
||||||
@ -98,9 +107,10 @@ namespace {
|
|||||||
|
|
||||||
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
|
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 nw = wells.number_of_wells;
|
||||||
const int nperf = wells.well_connpos[nw];
|
const int nperf = wells.well_connpos[nw];
|
||||||
const int dim = grid.dimensions;
|
const int dim = dimensions(grid);
|
||||||
V wdp = V::Zero(nperf,1);
|
V wdp = V::Zero(nperf,1);
|
||||||
assert(wdp.size() == rho.size());
|
assert(wdp.size() == rho.size());
|
||||||
|
|
||||||
@ -157,7 +167,8 @@ namespace {
|
|||||||
BlackoilState& state,
|
BlackoilState& state,
|
||||||
WellState& well_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();
|
const int np = state.numPhases();
|
||||||
|
|
||||||
well_flow_residual_.resize(np, ADB::null());
|
well_flow_residual_.resize(np, ADB::null());
|
||||||
@ -226,15 +237,16 @@ namespace {
|
|||||||
const BlackoilState& state,
|
const BlackoilState& state,
|
||||||
const WellState& well_state)
|
const WellState& well_state)
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
// Suppress warnings about "unused parameters".
|
// Suppress warnings about "unused parameters".
|
||||||
static_cast<void>(dt);
|
static_cast<void>(dt);
|
||||||
static_cast<void>(well_state);
|
static_cast<void>(well_state);
|
||||||
|
|
||||||
const int nc = grid_.number_of_cells;
|
const int nc = numCells(grid_);
|
||||||
const int np = state.numPhases();
|
const int np = state.numPhases();
|
||||||
const int nw = wells_.number_of_wells;
|
const int nw = wells_.number_of_wells;
|
||||||
const int nperf = wells_.well_connpos[nw];
|
const int nperf = wells_.well_connpos[nw];
|
||||||
const int dim = grid_.dimensions;
|
const int dim = dimensions(grid_);
|
||||||
|
|
||||||
const std::vector<int> cells = buildAllCells(nc);
|
const std::vector<int> cells = buildAllCells(nc);
|
||||||
|
|
||||||
@ -284,9 +296,9 @@ namespace {
|
|||||||
const BlackoilState& state,
|
const BlackoilState& state,
|
||||||
const WellState& well_state)
|
const WellState& well_state)
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
const V& pv = geo_.poreVolume();
|
const V& pv = geo_.poreVolume();
|
||||||
const int nc = grid_.number_of_cells;
|
const int nc = numCells(grid_); ;
|
||||||
const int np = state.numPhases();
|
const int np = state.numPhases();
|
||||||
const int nw = wells_.number_of_wells;
|
const int nw = wells_.number_of_wells;
|
||||||
const int nperf = wells_.well_connpos[nw];
|
const int nperf = wells_.well_connpos[nw];
|
||||||
@ -415,7 +427,8 @@ namespace {
|
|||||||
ImpesTPFAAD::solveJacobianSystem(BlackoilState& state,
|
ImpesTPFAAD::solveJacobianSystem(BlackoilState& state,
|
||||||
WellState& well_state) const
|
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 nw = wells_.number_of_wells;
|
||||||
// const int np = state.numPhases();
|
// const int np = state.numPhases();
|
||||||
|
|
||||||
@ -458,9 +471,10 @@ namespace {
|
|||||||
ImpesTPFAAD::computeFluxes(BlackoilState& state,
|
ImpesTPFAAD::computeFluxes(BlackoilState& state,
|
||||||
WellState& well_state) const
|
WellState& well_state) const
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
// This method computes state.faceflux(),
|
// This method computes state.faceflux(),
|
||||||
// well_state.perfRates() and well_state.perfPress().
|
// 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 np = state.numPhases();
|
||||||
const int nw = wells_.number_of_wells;
|
const int nw = wells_.number_of_wells;
|
||||||
const int nperf = wells_.well_connpos[nw];
|
const int nperf = wells_.well_connpos[nw];
|
||||||
@ -515,8 +529,8 @@ namespace {
|
|||||||
flux += face_mob * head;
|
flux += face_mob * head;
|
||||||
}
|
}
|
||||||
|
|
||||||
V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces);
|
V all_flux = superset(flux, ops_.internal_faces, numFaces(grid_));
|
||||||
std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().begin());
|
std::copy(all_flux.data(), all_flux.data() + numFaces(grid_), state.faceflux().begin());
|
||||||
|
|
||||||
perf_flux = -perf_flux; // well_state.perfRates() assumed to be inflows.
|
perf_flux = -perf_flux; // well_state.perfRates() assumed to be inflows.
|
||||||
std::copy(perf_flux.data(), perf_flux.data() + nperf, well_state.perfRates().begin());
|
std::copy(perf_flux.data(), perf_flux.data() + nperf, well_state.perfRates().begin());
|
||||||
|
@ -158,7 +158,14 @@ namespace Opm
|
|||||||
dm["saturation"] = &state.saturation();
|
dm["saturation"] = &state.saturation();
|
||||||
dm["pressure"] = &state.pressure();
|
dm["pressure"] = &state.pressure();
|
||||||
std::vector<double> cell_velocity;
|
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;
|
dm["velocity"] = &cell_velocity;
|
||||||
Opm::writeVtkData(grid, dm, vtkfile);
|
Opm::writeVtkData(grid, dm, vtkfile);
|
||||||
}
|
}
|
||||||
@ -174,7 +181,14 @@ namespace Opm
|
|||||||
dm["pressure"] = &state.pressure();
|
dm["pressure"] = &state.pressure();
|
||||||
dm["surfvolume"] = &state.surfacevol();
|
dm["surfvolume"] = &state.surfacevol();
|
||||||
std::vector<double> cell_velocity;
|
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;
|
dm["velocity"] = &cell_velocity;
|
||||||
|
|
||||||
// Write data (not grid) in Matlab format
|
// Write data (not grid) in Matlab format
|
||||||
@ -300,7 +314,7 @@ namespace Opm
|
|||||||
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
|
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
|
||||||
|
|
||||||
// Misc init.
|
// Misc init.
|
||||||
const int num_cells = grid.number_of_cells;
|
const int num_cells = AutoDiffGrid::numCells(grid);
|
||||||
allcells_.resize(num_cells);
|
allcells_.resize(num_cells);
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
allcells_[cell] = cell;
|
allcells_[cell] = cell;
|
||||||
|
@ -23,6 +23,7 @@
|
|||||||
#endif // HAVE_CONFIG_H
|
#endif // HAVE_CONFIG_H
|
||||||
|
|
||||||
#include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
|
#include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
|
||||||
|
#include <opm/autodiff/GridHelpers.hpp>
|
||||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
|
|
||||||
@ -400,7 +401,7 @@ namespace Opm
|
|||||||
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
|
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
|
||||||
|
|
||||||
// Misc init.
|
// Misc init.
|
||||||
const int num_cells = grid.number_of_cells;
|
const int num_cells = Opm::AutoDiffGrid::numCells(grid);
|
||||||
allcells_.resize(num_cells);
|
allcells_.resize(num_cells);
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
allcells_[cell] = cell;
|
allcells_[cell] = cell;
|
||||||
@ -498,11 +499,13 @@ namespace Opm
|
|||||||
double av_prev_press = 0.0;
|
double av_prev_press = 0.0;
|
||||||
double av_press = 0.0;
|
double av_press = 0.0;
|
||||||
double tot_vol = 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) {
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
|
av_prev_press += initial_pressure[cell]*
|
||||||
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
|
Opm::AutoDiffGrid::cellVolume(grid_, cell);
|
||||||
tot_vol += grid_.cell_volumes[cell];
|
av_press += state.pressure()[cell]*
|
||||||
|
Opm::AutoDiffGrid::cellVolume(grid_, cell);
|
||||||
|
tot_vol += Opm::AutoDiffGrid::cellVolume(grid_, cell);
|
||||||
}
|
}
|
||||||
// Renormalization constant
|
// Renormalization constant
|
||||||
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
||||||
|
@ -51,21 +51,22 @@ namespace Opm
|
|||||||
tol_(param.getDefault("nl_tolerance", 1e-9)),
|
tol_(param.getDefault("nl_tolerance", 1e-9)),
|
||||||
maxit_(param.getDefault("nl_maxiter", 30))
|
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);
|
allcells_.resize(nc);
|
||||||
for (int i = 0; i < nc; ++i) {
|
for (int i = 0; i < nc; ++i) {
|
||||||
allcells_[i] = i;
|
allcells_[i] = i;
|
||||||
}
|
}
|
||||||
if (gravity && gravity[grid_.dimensions - 1] != 0.0) {
|
if (gravity && gravity[dimensions(grid_) - 1] != 0.0) {
|
||||||
gravity_ = gravity[grid_.dimensions - 1];
|
gravity_ = gravity[dimensions(grid_) - 1];
|
||||||
for (int dd = 0; dd < grid_.dimensions - 1; ++dd) {
|
for (int dd = 0; dd < dimensions(grid_) - 1; ++dd) {
|
||||||
if (gravity[dd] != 0.0) {
|
if (gravity[dd] != 0.0) {
|
||||||
OPM_THROW(std::runtime_error, "TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
|
OPM_THROW(std::runtime_error, "TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
V htrans(grid.cell_facepos[grid.number_of_cells]);
|
V htrans(grid.cell_facepos[grid.number_of_cells]);
|
||||||
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
|
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());
|
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans.data());
|
||||||
transi_ = subset(trans, ops_.internal_faces);
|
transi_ = subset(trans, ops_.internal_faces);
|
||||||
}
|
}
|
||||||
@ -161,16 +162,17 @@ namespace Opm
|
|||||||
const double dt,
|
const double dt,
|
||||||
TwophaseState& state)
|
TwophaseState& state)
|
||||||
{
|
{
|
||||||
|
using namespace Opm::AutoDiffGrid;
|
||||||
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
|
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
|
||||||
typedef Eigen::Map<const V> Vec;
|
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);
|
const TwoCol s0 = Eigen::Map<const TwoCol>(state.saturation().data(), nc, 2);
|
||||||
double res_norm = 1e100;
|
double res_norm = 1e100;
|
||||||
const V sw0 = s0.leftCols<1>();
|
const V sw0 = s0.leftCols<1>();
|
||||||
// sw1 is the object that will be changed every Newton iteration.
|
// sw1 is the object that will be changed every Newton iteration.
|
||||||
// V sw1 = sw0;
|
// V sw1 = sw0;
|
||||||
V sw1 = 0.5*V::Ones(nc,1);
|
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();
|
const int num_internal = ops_.internal_faces.size();
|
||||||
V dflux = subset(dflux_all, ops_.internal_faces);
|
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 p1 = Vec(state.pressure().data(), nc, 1);
|
||||||
const V ndp = (ops_.ngrad * p1.matrix()).array();
|
const V ndp = (ops_.ngrad * p1.matrix()).array();
|
||||||
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DynArr;
|
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();
|
const V ndz = (ops_.ngrad * z.matrix()).array();
|
||||||
assert(num_internal == ndp.size());
|
assert(num_internal == ndp.size());
|
||||||
const double* density = props_.density();
|
const double* density = props_.density();
|
||||||
|
Loading…
Reference in New Issue
Block a user