OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-30 13:03:49 -06:00
Code Issues Packages Projects Releases Wiki Activity

Ported initStateEquil to using the GridHelpers.

Browse Source 
Currently the keyword EQUIL is not supported by the fully
implicit blackoil simulator when using CpGrid. This
commit is a first step towards this as it makes the
implementation of initStateEquil generic.
This commit is contained in:
Markus Blatt 2015-02-20 09:26:03 +01:00 committed by Andreas Lauser
parent 733061c943
commit dec7a93918
2 changed files with 87 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
opm/core/simulator/initStateEquil.hpp
View File

@ -1,5 +1,7 @@
/* /*
Copyright 2014 SINTEF ICT, Applied Mathematics. Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -20,6 +22,7 @@
#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
#define OPM_INITSTATEEQUIL_HEADER_INCLUDED #define OPM_INITSTATEEQUIL_HEADER_INCLUDED
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/simulator/EquilibrationHelpers.hpp> #include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp> #include <opm/core/props/BlackoilPropertiesInterface.hpp>
@ -59,7 +62,8 @@ namespace Opm
* \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/ */
void initStateEquil(const UnstructuredGrid& grid, template<class Grid>
void initStateEquil(const Grid& grid,
const BlackoilPropertiesInterface& props, const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr deck, const Opm::DeckConstPtr deck,
const Opm::EclipseStateConstPtr eclipseState, const Opm::EclipseStateConstPtr eclipseState,
@ -114,9 +118,9 @@ namespace Opm
* of pressure values in each cell in the current * of pressure values in each cell in the current
* equilibration region. * equilibration region.
*/ */
template <class Region, class CellRange> template <class Grid, class Region, class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
phasePressures(const UnstructuredGrid& G, phasePressures(const Grid& G,
const Region& reg, const Region& reg,
const CellRange& cells, const CellRange& cells,
const double grav = unit::gravity); const double grav = unit::gravity);
@ -147,9 +151,10 @@ namespace Opm
* \return Phase saturations, one vector for each phase, each containing * \return Phase saturations, one vector for each phase, each containing
* one saturation value per cell in the region. * one saturation value per cell in the region.
*/ */
template <class Region, class CellRange> template <class Grid, class Region, class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
phaseSaturations(const Region& reg, phaseSaturations(const Grid& grid,
const Region& reg,
const CellRange& cells, const CellRange& cells,
BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
const std::vector<double> swat_init, const std::vector<double> swat_init,
@ -175,8 +180,8 @@ namespace Opm
* \param[in] rs_func Rs as function of pressure and depth. * \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range. * \return Rs values, one for each cell in the 'cells' range.
*/ */
template <class CellRangeType> template <class Grid, class CellRangeType>
std::vector<double> computeRs(const UnstructuredGrid& grid, std::vector<double> computeRs(const Grid& grid,
const CellRangeType& cells, const CellRangeType& cells,
const std::vector<double> oil_pressure, const std::vector<double> oil_pressure,
const std::vector<double>& temperature, const std::vector<double>& temperature,
@ -230,19 +235,21 @@ namespace Opm
} }
} }
template<class Grid>
inline inline
std::vector<int> std::vector<int>
equilnum(const Opm::DeckConstPtr deck, equilnum(const Opm::DeckConstPtr deck,
const Opm::EclipseStateConstPtr eclipseState, const Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& G ) const Grid& G )
{ {
std::vector<int> eqlnum; std::vector<int> eqlnum;
if (deck->hasKeyword("EQLNUM")) { if (deck->hasKeyword("EQLNUM")) {
eqlnum.resize(G.number_of_cells); const int nc = UgGridHelpers::numCells(G);
eqlnum.resize(nc);
const std::vector<int>& e = const std::vector<int>& e =
eclipseState->getIntGridProperty("EQLNUM")->getData(); eclipseState->getIntGridProperty("EQLNUM")->getData();
const int* gc = G.global_cell; const int* gc = UgGridHelpers::globalCell(G);
for (int cell = 0; cell < G.number_of_cells; ++cell) { for (int cell = 0; cell < nc; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell]; const int deck_pos = (gc == NULL) ? cell : gc[cell];
eqlnum[cell] = e[deck_pos] - 1; eqlnum[cell] = e[deck_pos] - 1;
} }
@ -250,7 +257,7 @@ namespace Opm
else { else {
// No explicit equilibration region. // No explicit equilibration region.
// All cells in region zero. // All cells in region zero.
eqlnum.assign(G.number_of_cells, 0); eqlnum.assign(UgGridHelpers::numCells(G), 0);
} }
return eqlnum; return eqlnum;
@ -259,17 +266,18 @@ namespace Opm
class InitialStateComputer { class InitialStateComputer {
public: public:
template<class Grid>
InitialStateComputer(BlackoilPropertiesInterface& props, InitialStateComputer(BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr deck, const Opm::DeckConstPtr deck,
const Opm::EclipseStateConstPtr eclipseState, const Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& G , const Grid& G ,
const double grav = unit::gravity) const double grav = unit::gravity)
: pp_(props.numPhases(), : pp_(props.numPhases(),
std::vector<double>(G.number_of_cells)), std::vector<double>(UgGridHelpers::numCells(G))),
sat_(props.numPhases(), sat_(props.numPhases(),
std::vector<double>(G.number_of_cells)), std::vector<double>(UgGridHelpers::numCells(G))),
rs_(G.number_of_cells), rs_(UgGridHelpers::numCells(G)),
rv_(G.number_of_cells) rv_(UgGridHelpers::numCells(G))
{ {
// Get the equilibration records. // Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(deck); const std::vector<EquilRecord> rec = getEquil(deck);
@ -346,9 +354,10 @@ namespace Opm
// Check for presence of kw SWATINIT // Check for presence of kw SWATINIT
if (deck->hasKeyword("SWATINIT")) { if (deck->hasKeyword("SWATINIT")) {
const std::vector<double>& swat_init = eclipseState->getDoubleGridProperty("SWATINIT")->getData(); const std::vector<double>& swat_init = eclipseState->getDoubleGridProperty("SWATINIT")->getData();
swat_init_.resize(G.number_of_cells); const int nc = UgGridHelpers::numCells(G);
const int* gc = G.global_cell; swat_init_.resize(nc);
for (int c = 0; c < G.number_of_cells; ++c) { const int* gc = UgGridHelpers::globalCell(G);
for (int c = 0; c < nc; ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c]; const int deck_pos = (gc == NULL) ? c : gc[c];
swat_init_[c] = swat_init[deck_pos]; swat_init_[c] = swat_init[deck_pos];
} }
@ -382,12 +391,12 @@ namespace Opm
Vec rv_; Vec rv_;
Vec swat_init_; Vec swat_init_;
template <class RMap> template <class RMap, class Grid>
void void
calcPressSatRsRv(const RMap& reg , calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec , const std::vector< EquilRecord >& rec ,
Opm::BlackoilPropertiesInterface& props, Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G , const Grid& G ,
const double grav) const double grav)
{ {
for (typename RMap::RegionId for (typename RMap::RegionId

105
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -1,5 +1,7 @@
/* /*
Copyright 2014 SINTEF ICT, Applied Mathematics. Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -21,6 +23,7 @@
#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED #define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/initState.hpp> #include <opm/core/simulator/initState.hpp>
@ -300,16 +303,17 @@ namespace Opm
} // namespace PhaseIndex } // namespace PhaseIndex
namespace PhasePressure { namespace PhasePressure {
template <class PressFunction, template <class Grid,
class PressFunction,
class CellRange> class CellRange>
void void
assign(const UnstructuredGrid& G , assign(const Grid& G ,
const std::array<PressFunction, 2>& f , const std::array<PressFunction, 2>& f ,
const double split, const double split,
const CellRange& cells, const CellRange& cells,
std::vector<double>& p ) std::vector<double>& p )
{ {
const int nd = G.dimensions; const int nd = UgGridHelpers::dimensions(G);
enum { up = 0, down = 1 }; enum { up = 0, down = 1 };
@ -320,15 +324,16 @@ namespace Opm
{ {
assert (c < p.size()); assert (c < p.size());
const double z = G.cell_centroids[(*ci)*nd + (nd - 1)]; const double z = UgGridHelpers::cellCentroidCoordinate(G, *ci, nd-1);
p[c] = (z < split) ? f[up](z) : f[down](z); p[c] = (z < split) ? f[up](z) : f[down](z);
} }
} }
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
void void
water(const UnstructuredGrid& G , water(const Grid& G ,
const Region& reg , const Region& reg ,
const std::array<double,2>& span , const std::array<double,2>& span ,
const double grav , const double grav ,
@ -363,10 +368,11 @@ namespace Opm
assign(G, wpress, z0, cells, press); assign(G, wpress, z0, cells, press);
} }
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
void void
oil(const UnstructuredGrid& G , oil(const Grid& G ,
const Region& reg , const Region& reg ,
const std::array<double,2>& span , const std::array<double,2>& span ,
const double grav , const double grav ,
@ -418,10 +424,11 @@ namespace Opm
} }
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
void void
gas(const UnstructuredGrid& G , gas(const Grid& G ,
const Region& reg , const Region& reg ,
const std::array<double,2>& span , const std::array<double,2>& span ,
const double grav , const double grav ,
@ -463,10 +470,11 @@ namespace Opm
} }
} // namespace PhasePressure } // namespace PhasePressure
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
void void
equilibrateOWG(const UnstructuredGrid& G, equilibrateOWG(const Grid& G,
const Region& reg, const Region& reg,
const double grav, const double grav,
const std::array<double,2>& span, const std::array<double,2>& span,
@ -500,10 +508,11 @@ namespace Opm
namespace Equil { namespace Equil {
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
phasePressures(const UnstructuredGrid& G, phasePressures(const Grid& G,
const Region& reg, const Region& reg,
const CellRange& cells, const CellRange& cells,
const double grav) const double grav)
@ -515,18 +524,9 @@ namespace Opm
int ncell = 0; int ncell = 0;
{ {
// This code is only supported in three space dimensions // This code is only supported in three space dimensions
assert (G.dimensions == 3); assert (UgGridHelpers::dimensions(G) == 3);
const int nd = G.dimensions; const int nd = UgGridHelpers::dimensions(G);
// Define short-name aliases to reduce visual clutter.
const double* const nc = & G.node_coordinates[0];
const int* const cfp = & G.cell_facepos[0];
const int* const cf = & G.cell_faces[0];
const int* const fnp = & G.face_nodepos[0];
const int* const fn = & G.face_nodes[0];
// Define vertical span as // Define vertical span as
// //
@ -541,21 +541,24 @@ namespace Opm
// imposes the requirement that cell centroids are all // imposes the requirement that cell centroids are all
// within this vertical span. That requirement is not // within this vertical span. That requirement is not
// checked. // checked.
typename UgGridHelpers::Cell2FacesTraits<Grid>::Type
cell2Faces = UgGridHelpers::cell2Faces(G);
typename UgGridHelpers::Face2VerticesTraits<Grid>::Type
faceVertices = UgGridHelpers::face2Vertices(G);
for (typename CellRange::const_iterator for (typename CellRange::const_iterator
ci = cells.begin(), ce = cells.end(); ci = cells.begin(), ce = cells.end();
ci != ce; ++ci, ++ncell) ci != ce; ++ci, ++ncell)
{ {
for (const int for (auto fi=cell2Faces[*ci].begin(),
*fi = & cf[ cfp[*ci + 0] ], fe=cell2Faces[*ci].end();
*fe = & cf[ cfp[*ci + 1] ]; fi != fe;
fi != fe; ++fi) ++fi)
{ {
for (const int for (auto i = faceVertices[*fi].begin(), e=faceVertices[*fi].end();
*i = & fn[ fnp[*fi + 0] ],
*e = & fn[ fnp[*fi + 1] ];
i != e; ++i) i != e; ++i)
{ {
const double z = nc[(*i)*nd + (nd - 1)]; const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1];
if (z < span[0]) { span[0] = z; } if (z < span[0]) { span[0] = z; }
if (z > span[1]) { span[1] = z; } if (z > span[1]) { span[1] = z; }
@ -587,10 +590,11 @@ namespace Opm
return press; return press;
} }
template <class Region, template <class Grid,
class Region,
class CellRange> class CellRange>
std::vector<double> std::vector<double>
temperature(const UnstructuredGrid& /* G */, temperature(const Grid& /* G */,
const Region& /* reg */, const Region& /* reg */,
const CellRange& cells) const CellRange& cells)
{ {
@ -598,9 +602,9 @@ namespace Opm
return std::vector<double>(cells.size(), 273.15 + 20.0); return std::vector<double>(cells.size(), 273.15 + 20.0);
} }
template <class Region, class CellRange> template <class Grid, class Region, class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
phaseSaturations(const UnstructuredGrid& G, phaseSaturations(const Grid& G,
const Region& reg, const Region& reg,
const CellRange& cells, const CellRange& cells,
BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
@ -621,8 +625,6 @@ namespace Opm
double smin[BlackoilPhases::MaxNumPhases] = { 0.0 }; double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
double smax[BlackoilPhases::MaxNumPhases] = { 0.0 }; double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
const double* const cc = & G.cell_centroids[0];
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua]; const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour]; const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid]; const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
@ -637,8 +639,10 @@ namespace Opm
double sw = 0.0; double sw = 0.0;
if (water) { if (water) {
if (isConstPc(props,waterpos,cell)){ if (isConstPc(props,waterpos,cell)){
const int nd = G.dimensions; const int nd = UgGridHelpers::dimensions(G);
const double cellDepth = cc[nd * cell + nd-1]; const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
cell,
nd-1);
sw = satFromDepth(props,cellDepth,zwoc,waterpos,cell,false); sw = satFromDepth(props,cellDepth,zwoc,waterpos,cell,false);
phase_saturations[waterpos][local_index] = sw; phase_saturations[waterpos][local_index] = sw;
} }
@ -657,8 +661,10 @@ namespace Opm
double sg = 0.0; double sg = 0.0;
if (gas) { if (gas) {
if (isConstPc(props,gaspos,cell)){ if (isConstPc(props,gaspos,cell)){
const int nd = G.dimensions; const int nd = UgGridHelpers::dimensions(G);
const double cellDepth = cc[nd * cell + nd-1]; const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
cell,
nd-1);
sg = satFromDepth(props,cellDepth,zgoc,gaspos,cell,true); sg = satFromDepth(props,cellDepth,zgoc,gaspos,cell,true);
phase_saturations[gaspos][local_index] = sg; phase_saturations[gaspos][local_index] = sg;
} }
@ -742,19 +748,19 @@ namespace Opm
* \param[in] rs_func Rs as function of pressure and depth. * \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range. * \return Rs values, one for each cell in the 'cells' range.
*/ */
template <class CellRangeType> template <class Grid, class CellRangeType>
std::vector<double> computeRs(const UnstructuredGrid& grid, std::vector<double> computeRs(const Grid& grid,
const CellRangeType& cells, const CellRangeType& cells,
const std::vector<double> oil_pressure, const std::vector<double> oil_pressure,
const std::vector<double>& temperature, const std::vector<double>& temperature,
const Miscibility::RsFunction& rs_func, const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation) const std::vector<double> gas_saturation)
{ {
assert(grid.dimensions == 3); assert(UgGridHelpers::dimensions(grid) == 3);
std::vector<double> rs(cells.size()); std::vector<double> rs(cells.size());
int count = 0; int count = 0;
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) { for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
const double depth = grid.cell_centroids[3*(*it) + 2]; const double depth = UgGridHelpers::cellCentroidCoordinate(grid, *it, 2);
rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]); rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
} }
return rs; return rs;
@ -798,7 +804,8 @@ namespace Opm
* \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/ */
void initStateEquil(const UnstructuredGrid& grid, template<class Grid>
void initStateEquil(const Grid& grid,
BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr deck, const Opm::DeckConstPtr deck,
const Opm::EclipseStateConstPtr eclipseState, const Opm::EclipseStateConstPtr eclipseState,
@ -815,7 +822,7 @@ namespace Opm
state.saturation() = convertSats(isc.saturation()); state.saturation() = convertSats(isc.saturation());
state.gasoilratio() = isc.rs(); state.gasoilratio() = isc.rs();
state.rv() = isc.rv(); state.rv() = isc.rv();
initBlackoilSurfvolUsingRSorRV(grid, props, state); initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
} }