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

dunify and make the initialization parallel

Browse Source
This commit is contained in:
Tor Harald Sandve
2020-11-24 16:24:04 +01:00
parent f924ac7ffe
commit 6c06a72466
3 changed files with 161 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files

228
ebos/equil/initstateequil.hh
View File

@@ -1410,89 +1410,29 @@ invertCapPress(const double pc,
// ===========================================================================
template <typename Grid, typename CellRange>
void verticalExtent(const Grid& grid,
const CellRange& cells,
template <typename CellRange, typename Comm>
void verticalExtent(const CellRange& cells,
const std::vector<std::pair<double, double>>& cellZMinMax,
const Comm& comm,
std::array<double,2>& span)
{
// This code is only supported in three space dimensions
assert (Grid::dimensionworld == 3);
span[0] = std::numeric_limits<double>::max();
span[1] = std::numeric_limits<double>::lowest();
const int nd = Grid::dimensionworld;
// Define vertical span as
//
// [minimum(node depth(cells)), maximum(node depth(cells))]
//
// Note: We use a sledgehammer approach--looping all
// the nodes of all the faces of all the 'cells'--to
// compute those bounds. This necessarily entails
// visiting some nodes (and faces) multiple times.
//
// Note: The implementation of 'RK4IVP<>' implicitly
// imposes the requirement that cell centroids are all
// within this vertical span. That requirement is not
// checked.
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid);
auto faceVertices = Opm::UgGridHelpers::face2Vertices(grid);
for (typename CellRange::const_iterator
ci = cells.begin(), ce = cells.end();
ci != ce; ++ci)
{
for (auto fi = cell2Faces[*ci].begin(),
fe = cell2Faces[*ci].end();
fi != fe; ++fi)
{
for (auto i = faceVertices[*fi].begin(),
e = faceVertices[*fi].end();
i != e; ++i)
{
const double z = Opm::UgGridHelpers::
vertexCoordinates(grid, *i)[nd - 1];
if (z < span[0]) { span[0] = z; }
if (z > span[1]) { span[1] = z; }
}
}
for (const auto& cell : cells) {
if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; }
if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; }
}
}
template <typename Grid, typename CellID>
std::pair<double, double>
horizontalTopBottomDepths(const Grid& grid, const CellID cell)
{
const auto nd = Grid::dimensionworld;
auto c2f = Opm::UgGridHelpers::cell2Faces(grid);
auto top = std::numeric_limits<double>::max();
auto bot = std::numeric_limits<double>::lowest();
const auto topTag = 4; // Top face
const auto botTag = 5; // Bottom face
for (auto f = c2f[cell].begin(), e = c2f[cell].end(); f != e; ++f) {
const auto tag = Opm::UgGridHelpers::faceTag(grid, f);
if ((tag != topTag) && (tag != botTag)) {
// Not top/bottom face. Skip.
continue;
}
const auto depth = Opm::UgGridHelpers::
faceCentroid(grid, *f)[nd - 1];
if (tag == topTag) { // Top face
top = std::min(top, depth);
}
else { // Bottom face (tag == 5)
bot = std::max(bot, depth);
}
}
return std::make_pair(top, bot);
span[0] = comm.min(span[0]);
span[1] = comm.max(span[1]);
}
inline
@@ -1512,17 +1452,15 @@ void subdivisionCentrePoints(const double left,
}
}
template <typename Grid, typename CellID>
template <typename CellID>
std::vector<std::pair<double, double>>
horizontalSubdivision(const Grid& grid,
const CellID cell,
horizontalSubdivision(const CellID cell,
const std::pair<double, double> topbot,
const int numIntervals)
{
auto subdiv = std::vector<std::pair<double, double>>{};
subdiv.reserve(2 * numIntervals);
const auto topbot = horizontalTopBottomDepths(grid, cell);
if (topbot.first > topbot.second) {
throw std::out_of_range {
"Negative thickness (inverted top/bottom faces) in cell "
@@ -1536,6 +1474,59 @@ horizontalSubdivision(const Grid& grid,
return subdiv;
}
template <class Element>
inline double cellCenterDepth(const Element& element)
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
double zz = 0.0;
const Geometry geometry = element.geometry();
const int corners = geometry.corners();
for (int i=0; i < corners; ++i)
zz += geometry.corner(i)[zCoord];
return zz/corners;
}
template <class Element>
inline std::pair<double,double> cellZSpan(const Element& element)
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
double bot = 0.0;
double top = 0.0;
const Geometry geometry = element.geometry();
const int corners = geometry.corners();
assert(corners == 8);
for (int i=0; i < 4; ++i)
bot += geometry.corner(i)[zCoord];
for (int i=4; i < corners; ++i)
top += geometry.corner(i)[zCoord];
return std::make_pair(bot/4, top/4);
}
template <class Element>
inline std::pair<double,double> cellZMinMax(const Element& element)
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
const Geometry geometry = element.geometry();
const int corners = geometry.corners();
assert(corners == 8);
auto min = std::numeric_limits<double>::max();
auto max = std::numeric_limits<double>::lowest();
for (int i=0; i < corners; ++i) {
min = std::min(min, geometry.corner(i)[zCoord]);
max = std::max(max, geometry.corner(i)[zCoord]);
}
return std::make_pair(min, max);
}
} // namespace Details
namespace DeckDependent {
@@ -1571,24 +1562,28 @@ template<class TypeTag>
class InitialStateComputer
{
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
using Grid = GetPropType<TypeTag, Properties::Grid>;
public:
template<class MaterialLawManager>
InitialStateComputer(MaterialLawManager& materialLawManager,
const Opm::EclipseState& eclipseState,
const Grid& grid,
const Vanguard& vanguard,
const double grav = Opm::unit::gravity,
const bool applySwatInit = true)
: temperature_(grid.size(/*codim=*/0)),
saltConcentration_(grid.size(/*codim=*/0)),
: temperature_(vanguard.grid().size(/*codim=*/0)),
saltConcentration_(vanguard.grid().size(/*codim=*/0)),
pp_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))),
std::vector<double>(vanguard.grid().size(/*codim=*/0))),
sat_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))),
rs_(grid.size(/*codim=*/0)),
rv_(grid.size(/*codim=*/0))
std::vector<double>(vanguard.grid().size(/*codim=*/0))),
rs_(vanguard.grid().size(/*codim=*/0)),
rv_(vanguard.grid().size(/*codim=*/0))
{
//Check for presence of kw SWATINIT
if (applySwatInit) {
@@ -1597,6 +1592,10 @@ public:
}
}
//Querry cell depth, cell top-bottom.
updateCellProps_(vanguard);
const Grid& grid = vanguard.grid();
// Get the equilibration records.
const std::vector<Opm::EquilRecord> rec = getEquil(eclipseState);
@@ -1707,10 +1706,11 @@ public:
updateInitialTemperature_(eclipseState);
// EXTRACT the initial salt concentration
updateInitialSaltConcentration_(eclipseState, eqlmap, grid);
updateInitialSaltConcentration_(eclipseState, eqlmap);
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, materialLawManager, grid, grav);
const auto& comm = vanguard.gridView().comm();
calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
@@ -1727,13 +1727,14 @@ public:
const Vec& rv() const { return rv_; }
private:
void updateInitialTemperature_(const Opm::EclipseState& eclState)
{
this->temperature_ = eclState.fieldProps().get_double("TEMPI");
}
template <class RMap>
void updateInitialSaltConcentration_(const Opm::EclipseState& eclState, const RMap& reg, const Grid& grid)
void updateInitialSaltConcentration_(const Opm::EclipseState& eclState, const RMap& reg)
{
const int numEquilReg = rsFunc_.size();
saltVdTable_.resize(numEquilReg);
@@ -1754,7 +1755,7 @@ private:
const auto& cells = reg.cells(i);
for (const auto& cell : cells) {
const double depth = UgGridHelpers::cellCenterDepth(grid, cell);
const double depth = cellCenterDepth_[cell];
this->saltConcentration_[cell] = saltVdTable_[i].eval(depth);
}
}
@@ -1773,6 +1774,30 @@ private:
Vec rs_;
Vec rv_;
Vec swatInit_;
Vec cellCenterDepth_;
std::vector<std::pair<double,double>> cellZSpan_;
std::vector<std::pair<double,double>> cellZMinMax_;
template<class Vanguard>
void updateCellProps_(const Vanguard& vanguard)
{
const auto& gridView = vanguard.gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
int numElements = gridView.size(/*codim=*/0);
cellCenterDepth_.resize(numElements);
cellZSpan_.resize(numElements);
cellZMinMax_.resize(numElements);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
cellZSpan_[elemIdx] = Details::cellZSpan(element);
cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
}
}
template<class RMap>
void setRegionPvtIdx(const Opm::EclipseState& eclState, const RMap& reg)
@@ -1785,11 +1810,11 @@ private:
}
}
template <class RMap, class MaterialLawManager>
template <class RMap, class MaterialLawManager, class Comm>
void calcPressSatRsRv(const RMap& reg,
const std::vector< Opm::EquilRecord >& rec,
MaterialLawManager& materialLawManager,
const Grid& grid,
const Comm& comm,
const double grav)
{
using PhaseSat = Details::PhaseSaturations<
@@ -1800,16 +1825,17 @@ private:
auto psat = PhaseSat { materialLawManager, this->swatInit_ };
auto vspan = std::array<double, 2>{};
for (const auto& r : reg.activeRegions()) {
for (size_t r = 0; r < rec.size(); ++r) {
const auto& cells = reg.cells(r);
Details::verticalExtent(cells, cellZMinMax_, comm, vspan);
if (cells.empty()) {
Opm::OpmLog::warning("Equilibration region " + std::to_string(r + 1)
+ " has no active cells");
continue;
}
Details::verticalExtent(grid, cells, vspan);
const auto eqreg = EquilReg {
rec[r], this->rsFunc_[r], this->rvFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
};
@@ -1824,12 +1850,12 @@ private:
const auto acc = eqreg.equilibrationAccuracy();
if (acc == 0) {
// Centre-point method
this->equilibrateCellCentres(cells, eqreg, grid, ptable, psat);
this->equilibrateCellCentres(cells, eqreg, ptable, psat);
}
else if (acc < 0) {
// Horizontal subdivision
this->equilibrateHorizontal(cells, eqreg, -acc,
grid, ptable, psat);
ptable, psat);
}
}
}
@@ -1876,10 +1902,9 @@ private:
}
}
template <class CellRange, class Grid, class PressTable, class PhaseSat>
template <class CellRange, class PressTable, class PhaseSat>
void equilibrateCellCentres(const CellRange& cells,
const EquilReg& eqreg,
const Grid& grid,
const PressTable& ptable,
PhaseSat& psat)
{
@@ -1887,7 +1912,7 @@ private:
using CellID = std::remove_cv_t<std::remove_reference_t<
decltype(std::declval<CellPos>().cell)>>;
this->cellLoop(cells, [this, &eqreg, &grid, &ptable, &psat]
this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
(const CellID cell,
Details::PhaseQuantityValue& pressures,
Details::PhaseQuantityValue& saturations,
@@ -1895,7 +1920,7 @@ private:
double& Rv) -> void
{
const auto pos = CellPos {
cell, UgGridHelpers::cellCenterDepth(grid, cell)
cell, cellCenterDepth_[cell]
};
saturations = psat.deriveSaturations(pos, eqreg, ptable);
@@ -1911,11 +1936,10 @@ private:
});
}
template <class CellRange, class Grid, class PressTable, class PhaseSat>
template <class CellRange, class PressTable, class PhaseSat>
void equilibrateHorizontal(const CellRange& cells,
const EquilReg& eqreg,
const int acc,
const Grid& grid,
const PressTable& ptable,
PhaseSat& psat)
{
@@ -1923,7 +1947,7 @@ private:
using CellID = std::remove_cv_t<std::remove_reference_t<
decltype(std::declval<CellPos>().cell)>>;
this->cellLoop(cells, [this, acc, &eqreg, &grid, &ptable, &psat]
this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat]
(const CellID cell,
Details::PhaseQuantityValue& pressures,
Details::PhaseQuantityValue& saturations,
@@ -1934,7 +1958,7 @@ private:
saturations.reset();
auto totfrac = 0.0;
for (const auto& [depth, frac] : Details::horizontalSubdivision(grid, cell, acc)) {
for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
const auto pos = CellPos { cell, depth };
saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac);
@@ -1947,7 +1971,7 @@ private:
pressures /= totfrac;
const auto temp = this->temperature_[cell];
const auto cz = UgGridHelpers::cellCenterDepth(grid, cell);
const auto cz = cellCenterDepth_[cell];
Rs = eqreg.dissolutionCalculator()
(cz, pressures.oil, temp, saturations.gas);