OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-26 01:01:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

[refactor] Move cell centroid lookup magic to vanguard for reuse.

Browse Source 
In addition to transmissibilities We will need this for the tracers,
too, and should prevent code duplication.
This commit is contained in:
Markus Blatt 2021-09-30 11:50:22 +02:00
parent 2168b8c296
commit 457d270bdf
7 changed files with 83 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
ebos/eclalugridvanguard.hh
View File

@ -102,7 +102,6 @@ public:
~EclAluGridVanguard()
{
delete cartesianIndexMapper_;
delete equilCartesianIndexMapper_;
delete grid_;
delete equilGrid_;
@ -178,6 +177,18 @@ public:
const EquilCartesianIndexMapper& equilCartesianIndexMapper() const
{ return *equilCartesianIndexMapper_; }
/*!
* \brief Get function to query cell centroids for a distributed grid.
*
* Currently this only non-empty for a loadbalanced CpGrid.
* It is a function return the centroid for the given element
* index.
*/
std::function<std::array<double,dimensionworld>(int)>
cellCentroids() const
{
return this->cellCentroids_(cartesianIndexMapper_.get());
}
protected:
void createGrids_()
{
@ -210,7 +221,8 @@ protected:
grid_ = factory.convert(*equilGrid_, cartesianCellId_);
cartesianIndexMapper_ =
new CartesianIndexMapper(*grid_, cartesianDimension_, cartesianCellId_);
std::make_unique<CartesianIndexMapper>(*grid_, cartesianDimension_,
cartesianCellId_);
}
void filterConnections_()
@ -222,7 +234,7 @@ protected:
EquilGrid* equilGrid_;
std::vector<int> cartesianCellId_;
std::array<int,dimension> cartesianDimension_;
CartesianIndexMapper* cartesianIndexMapper_;
std::unique_ptr<CartesianIndexMapper> cartesianIndexMapper_;
EquilCartesianIndexMapper* equilCartesianIndexMapper_;
};

42
ebos/eclbasevanguard.hh
View File

@ -36,6 +36,7 @@
#include <opm/grid/common/GridEnums.hpp>
#include <opm/grid/common/CartesianIndexMapper.hpp>
#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquiferCell.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
#include <array>
@ -185,6 +186,7 @@ public:
protected:
static const int dimension = Grid::dimension;
static const int dimensionworld = Grid::dimensionworld;
using Element = typename GridView::template Codim<0>::Entity;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
@ -329,16 +331,6 @@ public:
{ return asImp_().equilCartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); }
/*!
* \brief Get the cell centroids for a distributed grid.
*
* Currently this only non-empty for a loadbalanced CpGrid.
*/
const std::vector<double>& cellCentroids() const
{
return centroids_;
}
/*!
* \brief Returns the depth of a degree of freedom [m]
*
@ -389,6 +381,36 @@ public:
}
protected:
/*!
* \brief Get function to query cell centroids for a distributed grid.
*
* Currently this only non-empty for a loadbalanced CpGrid.
* It is a function return the centroid for the given element
* index.
* \param cartMapper The cartesian index mapper for lookup of
* cartesian indices
*/
template<class CartMapper>
std::function<std::array<double,dimensionworld>(int)>
cellCentroids_(const CartMapper& cartMapper) const
{
return [this, cartMapper](int elemIdx) {
const auto& centroids = this->centroids_;
auto rank = this->gridView().comm().rank();
std::array<double,dimensionworld> centroid;
if (rank == 0) {
unsigned cartesianCellIdx = cartMapper.cartesianIndex(elemIdx);
centroid = this->eclState().getInputGrid().getCellCenter(cartesianCellIdx);
} else
{
std::copy(centroids.begin() + elemIdx * dimensionworld,
centroids.begin() + (elemIdx + 1) * dimensionworld,
centroid.begin());
}
return centroid;
};
}
void callImplementationInit()
{
asImp_().createGrids_();

14
ebos/eclcpgridvanguard.hh
View File

@ -90,6 +90,7 @@ public:
using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper, Scalar>;
static const int dimensionworld = Grid::dimensionworld;
private:
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
@ -138,7 +139,6 @@ public:
this->eclState(), this->parallelWells_);
#endif
this->allocCartMapper();
this->updateGridView_();
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
@ -149,6 +149,18 @@ public:
#endif
}
/*!
* \brief Get function to query cell centroids for a distributed grid.
*
* Currently this only non-empty for a loadbalanced CpGrid.
* It is a function return the centroid for the given element
* index.
*/
std::function<std::array<double,dimensionworld>(int)>
cellCentroids() const
{
return this->cellCentroids_(this->cartesianIndexMapper());
}
protected:
void createGrids_()

14
ebos/eclgenericcpgridvanguard.cc
View File

@ -92,7 +92,6 @@ void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::doLoadBalance_(Dun
// its edge weights. since this is (kind of) a layering violation and
// transmissibilities are relatively expensive to compute, we only do it if
// more than a single process is involved in the simulation.
cartesianIndexMapper_.reset(new CartesianIndexMapper(*grid_));
if (grid_->size(0))
{
this->allocTrans();
@ -175,8 +174,6 @@ void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::doLoadBalance_(Dun
}
grid_->switchToDistributedView();
cartesianIndexMapper_.reset();
// Calling Schedule::filterConnections would remove any perforated
// cells that exist only on other ranks even in the case of distributed wells
// But we need all connections to figure out the first cell of a well (e.g. for
@ -210,12 +207,6 @@ void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::distributeFieldPro
}
#endif
template<class ElementMapper, class GridView, class Scalar>
void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::allocCartMapper()
{
this->cartesianIndexMapper_.reset(new CartesianIndexMapper(this->grid()));
}
template<class ElementMapper, class GridView, class Scalar>
void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(EclipseState& eclState)
{
@ -258,6 +249,9 @@ void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Ecl
}
}
cartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*grid_);
#if HAVE_MPI
{
const bool has_numerical_aquifer = eclState.aquifer().hasNumericalAquifer();
@ -285,7 +279,7 @@ void EclGenericCpGridVanguard<ElementMapper,GridView,Scalar>::doCreateGrids_(Ecl
if (mpiRank == 0)
{
equilGrid_.reset(new Dune::CpGrid(*grid_));
equilCartesianIndexMapper_.reset(new CartesianIndexMapper(*equilGrid_));
equilCartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*equilGrid_);
std::vector<int> actnum = UgGridHelpers::createACTNUM(*grid_);
auto &field_props = eclState.fieldProps();

17
ebos/eclpolyhedralgridvanguard.hh
View File

@ -92,7 +92,7 @@ private:
typedef Grid* GridPointer;
typedef EquilGrid* EquilGridPointer;
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
typedef CartesianIndexMapper* CartesianIndexMapperPointer;
typedef std::unique_ptr<CartesianIndexMapper> CartesianIndexMapperPointer;
public:
EclPolyhedralGridVanguard(Simulator& simulator)
@ -104,7 +104,6 @@ public:
~EclPolyhedralGridVanguard()
{
delete cartesianIndexMapper_;
delete grid_;
}
@ -184,11 +183,23 @@ public:
return simulator_.problem().eclTransmissibilities();
}
/*!
* \brief Get function to query cell centroids for a distributed grid.
*
* Currently this only non-empty for a loadbalanced CpGrid.
* It is a function return the centroid for the given element
* index.
*/
std::function<std::array<double,dimensionworld>(int)>
cellCentroids() const
{
return this->cellCentroids_(this->cartesianIndexMapper());
}
protected:
void createGrids_()
{
grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true));
cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
cartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*grid_);
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
}

12
ebos/ecltransmissibility.cc
View File

@ -88,7 +88,7 @@ EclTransmissibility(const EclipseState& eclState,
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const Grid& grid,
const std::vector<double>& centroids,
std::function<std::array<double,dimWorld>(int)> centroids,
bool enableEnergy,
bool enableDiffusivity)
: eclState_(eclState)
@ -174,15 +174,7 @@ update(bool global)
// compute the axis specific "centroids" used for the transmissibilities. for
// consistency with the flow simulator, we use the element centers as
// computed by opm-parser's Opm::EclipseGrid class for all axes.
std::array<double, 3> centroid;
if (gridView_.comm().rank() == 0) {
const auto& eclGrid = eclState_.getInputGrid();
unsigned cartesianCellIdx = cartMapper_.cartesianIndex(elemIdx);
centroid = eclGrid.getCellCenter(cartesianCellIdx);
} else
std::copy(centroids_.begin() + centroidIdx * dimWorld,
centroids_.begin() + (centroidIdx + 1) * dimWorld,
centroid.begin());
std::array<double, dimWorld> centroid = centroids_(elemIdx);
for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx)
for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx)

5
ebos/ecltransmissibility.hh
View File

@ -38,6 +38,7 @@
#include <tuple>
#include <vector>
#include <unordered_map>
#include <functional>
namespace Opm {
@ -58,7 +59,7 @@ public:
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const Grid& grid,
const std::vector<double>& centroids,
std::function<std::array<double,dimWorld>(int)> centroids,
bool enableEnergy,
bool enableDiffusivity);
@ -236,7 +237,7 @@ protected:
const GridView& gridView_;
const Dune::CartesianIndexMapper<Grid>& cartMapper_;
const Grid& grid_;
const std::vector<double>& centroids_;
std::function<std::array<double,dimWorld>(int)> centroids_;
Scalar transmissibilityThreshold_;
std::map<std::pair<unsigned, unsigned>, Scalar> transBoundary_;
std::map<std::pair<unsigned, unsigned>, Scalar> thermalHalfTransBoundary_;