OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 10:40:21 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2975 from totto82/dunify_aquifers

Browse Source 
Dunify aquifers
This commit is contained in:
Atgeirr Flø Rasmussen 2020-12-15 09:46:21 +01:00 committed by GitHub
commit 81419aa001
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 198 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
ebos/eclbasevanguard.hh
View File

@ -169,6 +169,7 @@ class EclBaseVanguard : public BaseVanguard<TypeTag>
using Implementation = GetPropType<TypeTag, Properties::Vanguard>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
@ -178,6 +179,8 @@ public:
protected:
static const int dimension = Grid::dimension;
using Element = typename GridView::template Codim<0>::Entity;
public:
/*!
@ -584,6 +587,16 @@ public:
return cartIndex;
}
/*!
* \brief Return compressed index from cartesian index
*
*/
int compressedIndex(int cartesianCellIdx) const
{
int index = cartesianToCompressed_[cartesianCellIdx];
return index;
}
/*!
* \brief Extract Cartesian index triplet (i,j,k) of an active cell.
*
@ -628,6 +641,17 @@ public:
return centroids_;
}
/*!
* \brief Returns the depth of an degree of freedom [m]
*
* For ECL problems this is defined as the average of the depth of an element and is
* thus slightly different from the depth of an element's centroid.
*/
Scalar cellCenterDepth(unsigned globalSpaceIdx) const
{
return cellCenterDepth_[globalSpaceIdx];
}
/*!
* \brief Get the number of cells in the global leaf grid view.
* \warn This is a collective operation that needs to be called
@ -664,6 +688,31 @@ protected:
}
}
}
void updateCartesianToCompressedMapping_()
{
size_t num_cells = asImp_().grid().leafGridView().size(0);
cartesianToCompressed_.resize(cartesianSize(), -1);
for (unsigned i = 0; i < num_cells; ++i) {
unsigned cartesianCellIdx = cartesianIndex(i);
cartesianToCompressed_[cartesianCellIdx] = i;
}
}
void updateCellDepths_()
{
int numCells = this->gridView().size(/*codim=*/0);
cellCenterDepth_.resize(numCells);
ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
auto elemIt = this->gridView().template begin</*codim=*/0>();
const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = cellCenterDepth(element);
}
}
private:
void updateOutputDir_()
{
@ -693,6 +742,20 @@ private:
ioConfig.setEclCompatibleRST(!EWOMS_GET_PARAM(TypeTag, bool, EnableOpmRstFile));
}
Scalar cellCenterDepth(const Element& element) const
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
Scalar 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/Scalar(corners);
}
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
@ -734,12 +797,25 @@ protected:
* Empty otherwise. Used by EclTransmissibilty.
*/
std::vector<double> centroids_;
/*! \brief Mapping between cartesian and compressed cells.
* It is initialized the first time it is called
*/
std::vector<int> cartesianToCompressed_;
/*! \brief Cell center depths computed
* from averaging cell corner depths
*/
std::vector<Scalar> cellCenterDepth_;
/*! \brief information about wells in parallel
*
* For each well in the model there is an entry with its name
* and a boolean indicating whether it perforates local cells.
*/
std::vector<std::pair<std::string,bool>> parallelWells_;
};
template <class TypeTag>

4
ebos/eclcpgridvanguard.hh
View File

@ -244,6 +244,10 @@ public:
cartesianIndexMapper_.reset(new CartesianIndexMapper(*grid_));
this->updateGridView_();
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
#if HAVE_MPI
if (mpiSize > 1) {
try

2
ebos/eclpolyhedralgridvanguard.hh
View File

@ -185,6 +185,8 @@ protected:
{
grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true));
cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
}
void filterConnections_()

41
ebos/eclproblem.hh
View File

@ -898,7 +898,6 @@ public:
maxOilSaturation_.resize(numDof, 0.0);
}
updateElementDepths_();
readRockParameters_();
readMaterialParameters_();
readThermalParameters_();
@ -1468,9 +1467,10 @@ public:
Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return elementCenterDepth_[globalSpaceIdx];
return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
}
/*!
* \copydoc BlackoilProblem::rockCompressibility
*/
@ -2241,40 +2241,6 @@ private:
}
Scalar cellCenterDepth(const Element& element) const
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
Scalar 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/Scalar(corners);
}
void updateElementDepths_()
{
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& gridView = vanguard.gridView();
const auto& elemMapper = this->elementMapper();;
int numElements = gridView.size(/*codim=*/0);
elementCenterDepth_.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);
elementCenterDepth_[elemIdx] = cellCenterDepth(element);
}
}
// update the parameters needed for DRSDT and DRVDT
void updateCompositionChangeLimits_()
{
@ -2475,7 +2441,7 @@ private:
if (!rockTableIdx_.empty()) {
tableIdx = rockTableIdx_[elemIdx];
}
overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(elementCenterDepth_[elemIdx], /*extrapolation=*/true);
overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(vanguard.cellCenterDepth(elemIdx), /*extrapolation=*/true);
}
}
@ -3313,7 +3279,6 @@ private:
static std::string briefDescription_;
std::array<std::vector<Scalar>, 2> referencePorosity_;
std::vector<Scalar> elementCenterDepth_;
EclTransmissibility<TypeTag> transmissibilities_;
std::shared_ptr<EclMaterialLawManager> materialLawManager_;

55
opm/simulators/aquifers/AquiferCarterTracy.hpp
View File

@ -46,14 +46,14 @@ public:
using typename Base::RateVector;
using typename Base::Scalar;
using typename Base::Simulator;
using typename Base::ElementMapper;
using Base::waterCompIdx;
using Base::waterPhaseIdx;
AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator,
const AquiferCT::AQUCT_data& aquct_data)
: Base(aquct_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
: Base(aquct_data.aquiferID, connections, ebosSimulator)
, aquct_data_(aquct_data)
{
}
@ -91,18 +91,10 @@ protected:
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
inline void initializeConnections() override
{
const auto& eclState = this->ebos_simulator_.vanguard().eclState();
const auto& ugrid = this->ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
// We hack the cell depth values for now. We can actually get it from elementcontext pos
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
@ -110,18 +102,41 @@ protected:
Scalar denom_face_areas = 0.;
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < this->size(); ++idx) {
const int cell_index = this->cartesian_to_compressed_.at(this->connections_[idx].global_index);
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
if (cell_index < 0) //the global_index is not part of this grid
continue;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get default areas for all intersections
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
const auto cellFacesRange = cell2Faces[cell_index];
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
// The index of the face in the compressed grid
const int faceIdx = *cellFaceIter;
// only deal with connections given by the aquifer
if( idx < 0)
continue;
// the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
switch (faceTag) {
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = Opm::FaceDir::XMinus;
break;
@ -146,12 +161,10 @@ protected:
}
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_.at(idx) = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
}
}
auto cellCenter = grid.getCellCenter(this->connections_[idx].global_index);
this->cell_depth_.at(idx) = cellCenter[2];
}
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());

56
opm/simulators/aquifers/AquiferFetkovich.hpp
View File

@ -47,15 +47,15 @@ public:
using typename Base::RateVector;
using typename Base::Scalar;
using typename Base::Simulator;
using typename Base::ElementMapper;
using Base::waterCompIdx;
using Base::waterPhaseIdx;
AquiferFetkovich(const std::vector<Aquancon::AquancCell>& connections,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator,
const Aquifetp::AQUFETP_data& aqufetp_data)
: Base(aqufetp_data.aquiferID, connections, cartesian_to_compressed, ebosSimulator)
: Base(aqufetp_data.aquiferID, connections, ebosSimulator)
, aqufetp_data_(aqufetp_data)
{
}
@ -94,18 +94,10 @@ protected:
inline void initializeConnections() override
{
const auto& eclState = this->ebos_simulator_.vanguard().eclState();
const auto& ugrid = this->ebos_simulator_.vanguard().grid();
const auto& grid = eclState.getInputGrid();
// We hack the cell depth values for now. We can actually get it from elementcontext pos
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
@ -114,22 +106,40 @@ protected:
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->cartesian_to_compressed_.at(global_index);
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
if (cell_index < 0) //the global_index is not part of this grid
continue;
this->cellToConnectionIdx_[cell_index] = idx;
const auto cellCenter = grid.getCellCenter(global_index);
this->cell_depth_.at(idx) = cellCenter[2];
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get areas for all connections
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
if (!this->connections_[idx].influx_coeff.first) { // influx_coeff is defaulted
const auto cellFacesRange = cell2Faces[cell_index];
for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
// The index of the face in the compressed grid
const int faceIdx = *cellFaceIter;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// the logically-Cartesian direction of the face
const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
// only deal with grid boundaries
if (!intersection.boundary())
continue;
switch (faceTag) {
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = Opm::FaceDir::XMinus;
break;
@ -150,11 +160,11 @@ protected:
break;
default:
OPM_THROW(Opm::NumericalIssue,
"Initialization of Aquifer problem. Make sure faceTag is correctly defined");
}
"Initialization of Aquifer Fetkovich problem. Make sure faceTag is correctly defined"); }
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
break;
}
}

32
opm/simulators/aquifers/AquiferInterface.hpp
View File

@ -50,6 +50,7 @@ public:
using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
@ -75,12 +76,10 @@ public:
// Constructor
AquiferInterface(int aqID,
const std::vector<Aquancon::AquancCell>& connections,
const std::unordered_map<int, int>& cartesian_to_compressed,
const Simulator& ebosSimulator)
: aquiferID(aqID)
, connections_(connections)
, ebos_simulator_(ebosSimulator)
, cartesian_to_compressed_(cartesian_to_compressed)
{
}
@ -201,29 +200,13 @@ protected:
rhow_.at(idx) = fs.density(waterPhaseIdx);
}
template <class faceCellType, class ugridType>
inline double getFaceArea(const faceCellType& faceCells,
const ugridType& ugrid,
const int faceIdx,
const int idx) const
template <class Intersection>
inline double getFaceArea(const Intersection& intersection,
unsigned idx) const
{
// Check now if the face is outside of the reservoir, or if it adjoins an inactive cell
// Do not make the connection if the product of the two cellIdx > 0. This is because the
// face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell
// adjoining)
double faceArea = 0.;
const auto cellNeighbour0 = faceCells(faceIdx, 0);
const auto cellNeighbour1 = faceCells(faceIdx, 1);
const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
const auto calculatedFaceArea
= (!this->connections_[idx].influx_coeff.first) ? defaultFaceArea : this->connections_[idx].influx_coeff.second;
faceArea = (cellNeighbour0 * cellNeighbour1 > 0) ? 0. : calculatedFaceArea;
if (cellNeighbour1 == 0) {
faceArea = (cellNeighbour0 < 0) ? faceArea : 0.;
} else if (cellNeighbour0 == 0) {
faceArea = (cellNeighbour1 < 0) ? faceArea : 0.;
}
return faceArea;
const auto& geometry = intersection.geometry();
const auto defaultFaceArea = geometry.volume();
return (!this->connections_[idx].influx_coeff.first) ? defaultFaceArea : this->connections_[idx].influx_coeff.second;
}
virtual void endTimeStep() = 0;
@ -231,7 +214,6 @@ protected:
const int aquiferID;
const std::vector<Aquancon::AquancCell> connections_;
const Simulator& ebos_simulator_;
const std::unordered_map<int, int> cartesian_to_compressed_;
// Grid variables
std::vector<Scalar> faceArea_connected_;

34
opm/simulators/aquifers/BlackoilAquiferModel.hpp
View File

@ -35,13 +35,45 @@
#include <opm/simulators/aquifers/AquiferCarterTracy.hpp>
#include <opm/simulators/aquifers/AquiferFetkovich.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
#include <opm/material/densead/Math.hpp>
#include <vector>
#include <type_traits>
namespace Opm
{
template<class Grid>
class SupportsFaceTag
: public std::bool_constant<false>
{};
template<>
class SupportsFaceTag<Dune::CpGrid>
: public std::bool_constant<true>
{};
template<>
class SupportsFaceTag<Dune::PolyhedralGrid<3, 3>>
: public std::bool_constant<true>
{};
#if HAVE_DUNE_ALUGRID
template<>
class SupportsFaceTag<Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>
: public std::bool_constant<true>
{};
#endif
/// Class for handling the blackoil well model.
template <typename TypeTag>
class BlackoilAquiferModel
@ -49,6 +81,7 @@ class BlackoilAquiferModel
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
public:
explicit BlackoilAquiferModel(Simulator& simulator);
@ -83,7 +116,6 @@ protected:
Simulator& simulator_;
std::unordered_map<int, int> cartesian_to_compressed_;
// TODO: probably we can use one variable to store both types of aquifers, because
// they share the base class
mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy;

15
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
View File

@ -18,7 +18,6 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/grid/utility/cartesianToCompressed.hpp>
namespace Opm
{
@ -26,6 +25,10 @@ template <typename TypeTag>
BlackoilAquiferModel<TypeTag>::BlackoilAquiferModel(Simulator& simulator)
: simulator_(simulator)
{
// Grid needs to support Facetag
using Grid = std::remove_const_t<std::remove_reference_t<decltype(simulator.vanguard().grid())>>;
static_assert(SupportsFaceTag<Grid>::value, "Grid has to support assumptions about face tag.");
init();
}
@ -170,21 +173,15 @@ BlackoilAquiferModel<TypeTag>::init()
throw std::runtime_error("Aquifers currently do not work in parallel.");
// Get all the carter tracy aquifer properties data and put it in aquifers vector
const auto& ugrid = simulator_.vanguard().grid();
const int number_of_cells = simulator_.gridView().size(0);
cartesian_to_compressed_ = cartesianToCompressed(number_of_cells,
Opm::UgGridHelpers::globalCell(ugrid));
const auto& connections = aquifer.connections();
for (const auto& aq : aquifer.ct()) {
aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
cartesian_to_compressed_, this->simulator_, aq);
this->simulator_, aq);
}
for (const auto& aq : aquifer.fetp()) {
aquifers_Fetkovich.emplace_back(connections[aq.aquiferID],
cartesian_to_compressed_, this->simulator_, aq);
this->simulator_, aq);
}
}
template <typename TypeTag>