diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 8e86c4ef9..14deaf677 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -1,5 +1,7 @@ /* 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). @@ -20,6 +22,7 @@ #ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED #define OPM_INITSTATEEQUIL_HEADER_INCLUDED +#include #include #include #include @@ -59,7 +62,8 @@ namespace Opm * \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ - void initStateEquil(const UnstructuredGrid& grid, + template + void initStateEquil(const Grid& grid, const BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, @@ -114,9 +118,9 @@ namespace Opm * of pressure values in each cell in the current * equilibration region. */ - template + template std::vector< std::vector > - phasePressures(const UnstructuredGrid& G, + phasePressures(const Grid& G, const Region& reg, const CellRange& cells, const double grav = unit::gravity); @@ -147,9 +151,10 @@ namespace Opm * \return Phase saturations, one vector for each phase, each containing * one saturation value per cell in the region. */ - template + template std::vector< std::vector > - phaseSaturations(const Region& reg, + phaseSaturations(const Grid& grid, + const Region& reg, const CellRange& cells, BlackoilPropertiesInterface& props, const std::vector swat_init, @@ -175,8 +180,8 @@ namespace Opm * \param[in] rs_func Rs as function of pressure and depth. * \return Rs values, one for each cell in the 'cells' range. */ - template - std::vector computeRs(const UnstructuredGrid& grid, + template + std::vector computeRs(const Grid& grid, const CellRangeType& cells, const std::vector oil_pressure, const std::vector& temperature, @@ -230,19 +235,21 @@ namespace Opm } } + template inline std::vector equilnum(const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, - const UnstructuredGrid& G ) + const Grid& G ) { std::vector eqlnum; if (deck->hasKeyword("EQLNUM")) { - eqlnum.resize(G.number_of_cells); + const int nc = UgGridHelpers::numCells(G); + eqlnum.resize(nc); const std::vector& e = eclipseState->getIntGridProperty("EQLNUM")->getData(); - const int* gc = G.global_cell; - for (int cell = 0; cell < G.number_of_cells; ++cell) { + const int* gc = UgGridHelpers::globalCell(G); + for (int cell = 0; cell < nc; ++cell) { const int deck_pos = (gc == NULL) ? cell : gc[cell]; eqlnum[cell] = e[deck_pos] - 1; } @@ -250,7 +257,7 @@ namespace Opm else { // No explicit equilibration region. // All cells in region zero. - eqlnum.assign(G.number_of_cells, 0); + eqlnum.assign(UgGridHelpers::numCells(G), 0); } return eqlnum; @@ -259,17 +266,18 @@ namespace Opm class InitialStateComputer { public: + template InitialStateComputer(BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, - const UnstructuredGrid& G , + const Grid& G , const double grav = unit::gravity) : pp_(props.numPhases(), - std::vector(G.number_of_cells)), + std::vector(UgGridHelpers::numCells(G))), sat_(props.numPhases(), - std::vector(G.number_of_cells)), - rs_(G.number_of_cells), - rv_(G.number_of_cells) + std::vector(UgGridHelpers::numCells(G))), + rs_(UgGridHelpers::numCells(G)), + rv_(UgGridHelpers::numCells(G)) { // Get the equilibration records. const std::vector rec = getEquil(deck); @@ -346,9 +354,10 @@ namespace Opm // Check for presence of kw SWATINIT if (deck->hasKeyword("SWATINIT")) { const std::vector& swat_init = eclipseState->getDoubleGridProperty("SWATINIT")->getData(); - swat_init_.resize(G.number_of_cells); - const int* gc = G.global_cell; - for (int c = 0; c < G.number_of_cells; ++c) { + const int nc = UgGridHelpers::numCells(G); + swat_init_.resize(nc); + const int* gc = UgGridHelpers::globalCell(G); + for (int c = 0; c < nc; ++c) { const int deck_pos = (gc == NULL) ? c : gc[c]; swat_init_[c] = swat_init[deck_pos]; } @@ -382,12 +391,12 @@ namespace Opm Vec rv_; Vec swat_init_; - template + template void calcPressSatRsRv(const RMap& reg , const std::vector< EquilRecord >& rec , Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , + const Grid& G , const double grav) { for (typename RMap::RegionId diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index e85a67042..9d0c6e7dd 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -1,5 +1,7 @@ /* 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). @@ -21,6 +23,7 @@ #define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED #include +#include #include #include @@ -300,16 +303,17 @@ namespace Opm } // namespace PhaseIndex namespace PhasePressure { - template void - assign(const UnstructuredGrid& G , + assign(const Grid& G , const std::array& f , const double split, const CellRange& cells, std::vector& p ) { - const int nd = G.dimensions; + const int nd = UgGridHelpers::dimensions(G); enum { up = 0, down = 1 }; @@ -320,15 +324,16 @@ namespace Opm { 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); } } - template void - water(const UnstructuredGrid& G , + water(const Grid& G , const Region& reg , const std::array& span , const double grav , @@ -363,10 +368,11 @@ namespace Opm assign(G, wpress, z0, cells, press); } - template void - oil(const UnstructuredGrid& G , + oil(const Grid& G , const Region& reg , const std::array& span , const double grav , @@ -418,10 +424,11 @@ namespace Opm } - template void - gas(const UnstructuredGrid& G , + gas(const Grid& G , const Region& reg , const std::array& span , const double grav , @@ -463,10 +470,11 @@ namespace Opm } } // namespace PhasePressure - template void - equilibrateOWG(const UnstructuredGrid& G, + equilibrateOWG(const Grid& G, const Region& reg, const double grav, const std::array& span, @@ -500,10 +508,11 @@ namespace Opm namespace Equil { - template std::vector< std::vector > - phasePressures(const UnstructuredGrid& G, + phasePressures(const Grid& G, const Region& reg, const CellRange& cells, const double grav) @@ -515,18 +524,9 @@ namespace Opm int ncell = 0; { // This code is only supported in three space dimensions - assert (G.dimensions == 3); + assert (UgGridHelpers::dimensions(G) == 3); - const int nd = G.dimensions; - - // 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]; + const int nd = UgGridHelpers::dimensions(G); // Define vertical span as // @@ -541,21 +541,24 @@ namespace Opm // imposes the requirement that cell centroids are all // within this vertical span. That requirement is not // checked. + typename UgGridHelpers::Cell2FacesTraits::Type + cell2Faces = UgGridHelpers::cell2Faces(G); + typename UgGridHelpers::Face2VerticesTraits::Type + faceVertices = UgGridHelpers::face2Vertices(G); + for (typename CellRange::const_iterator ci = cells.begin(), ce = cells.end(); ci != ce; ++ci, ++ncell) { - for (const int - *fi = & cf[ cfp[*ci + 0] ], - *fe = & cf[ cfp[*ci + 1] ]; - fi != fe; ++fi) + for (auto fi=cell2Faces[*ci].begin(), + fe=cell2Faces[*ci].end(); + fi != fe; + ++fi) { - for (const int - *i = & fn[ fnp[*fi + 0] ], - *e = & fn[ fnp[*fi + 1] ]; + for (auto i = faceVertices[*fi].begin(), e=faceVertices[*fi].end(); 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[1]) { span[1] = z; } @@ -587,10 +590,11 @@ namespace Opm return press; } - template std::vector - temperature(const UnstructuredGrid& /* G */, + temperature(const Grid& /* G */, const Region& /* reg */, const CellRange& cells) { @@ -598,9 +602,9 @@ namespace Opm return std::vector(cells.size(), 273.15 + 20.0); } - template + template std::vector< std::vector > - phaseSaturations(const UnstructuredGrid& G, + phaseSaturations(const Grid& G, const Region& reg, const CellRange& cells, BlackoilPropertiesInterface& props, @@ -621,8 +625,6 @@ namespace Opm double smin[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 gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour]; const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid]; @@ -637,8 +639,10 @@ namespace Opm double sw = 0.0; if (water) { if (isConstPc(props,waterpos,cell)){ - const int nd = G.dimensions; - const double cellDepth = cc[nd * cell + nd-1]; + const int nd = UgGridHelpers::dimensions(G); + const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G, + cell, + nd-1); sw = satFromDepth(props,cellDepth,zwoc,waterpos,cell,false); phase_saturations[waterpos][local_index] = sw; } @@ -657,8 +661,10 @@ namespace Opm double sg = 0.0; if (gas) { if (isConstPc(props,gaspos,cell)){ - const int nd = G.dimensions; - const double cellDepth = cc[nd * cell + nd-1]; + const int nd = UgGridHelpers::dimensions(G); + const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G, + cell, + nd-1); sg = satFromDepth(props,cellDepth,zgoc,gaspos,cell,true); phase_saturations[gaspos][local_index] = sg; } @@ -742,19 +748,19 @@ namespace Opm * \param[in] rs_func Rs as function of pressure and depth. * \return Rs values, one for each cell in the 'cells' range. */ - template - std::vector computeRs(const UnstructuredGrid& grid, + template + std::vector computeRs(const Grid& grid, const CellRangeType& cells, const std::vector oil_pressure, const std::vector& temperature, const Miscibility::RsFunction& rs_func, const std::vector gas_saturation) { - assert(grid.dimensions == 3); + assert(UgGridHelpers::dimensions(grid) == 3); std::vector rs(cells.size()); int count = 0; 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]); } return rs; @@ -798,7 +804,8 @@ namespace Opm * \param[in] deck Simulation deck, used to obtain EQUIL and related data. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ - void initStateEquil(const UnstructuredGrid& grid, + template + void initStateEquil(const Grid& grid, BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, @@ -815,7 +822,7 @@ namespace Opm state.saturation() = convertSats(isc.saturation()); state.gasoilratio() = isc.rs(); state.rv() = isc.rv(); - initBlackoilSurfvolUsingRSorRV(grid, props, state); + initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state); }