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

rebasing

Browse Source
This commit is contained in:
Elyes Ahmed 2021-12-01 14:00:21 +01:00
parent 2e11d7a13a
commit 5c2b60bcd0
35 changed files with 898 additions and 243 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
CMakeLists.txt
View File

@ -30,6 +30,7 @@ option(ENABLE_FPGA "Enable FPGA kernels integration?" OFF)
option(USE_CHOW_PATEL_ILU "Use the iterative ILU by Chow and Patel?" OFF)
option(USE_CHOW_PATEL_ILU_GPU "Run iterative ILU decomposition on GPU? Requires USE_CHOW_PATEL_ILU" OFF)
option(USE_CHOW_PATEL_ILU_GPU_PARALLEL "Try to use more parallelism on the GPU during the iterative ILU decomposition? Requires USE_CHOW_PATEL_ILU_GPU" OFF)
option(BUILD_FLOW_ALU_GRID "Build flow blackoil with alu grid" OFF)
# The following was copied from CMakeLists.txt in opm-common.
# TODO: factor out the common parts in opm-common and opm-simulator as a cmake module
@ -453,7 +454,6 @@ if(TARGET fmt::fmt)
endif()
set_property(TARGET flow_libblackoil_poly PROPERTY EXCLUDE_FROM_ALL ${FLOW_POLY_ONLY_DEFAULT_ENABLE_IF})
# the production oriented general-purpose ECL simulator
opm_add_test(flow
ONLY_COMPILE
ALWAYS_ENABLE
@ -489,6 +489,23 @@ opm_add_test(flow_distribute_z
${FLOW_TGTS}
$<TARGET_OBJECTS:moduleVersion>
)
if (NOT BUILD_FLOW_ALU_GRID)
set(FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF "FALSE")
else()
set(FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF "TRUE")
endif()
# the production oriented general-purpose ECL simulator
opm_add_test(flow_alugrid
ONLY_COMPILE
ALWAYS_ENABLE
DEFAULT_ENABLE_IF ${FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF}
DEPENDS opmsimulators
LIBRARIES opmsimulators
SOURCES
flow/flow_blackoil_alugrid.cpp
$<TARGET_OBJECTS:moduleVersion>)
target_compile_definitions(flow_alugrid PRIVATE USE_ALUGRID)
if (BUILD_FLOW)
install(TARGETS flow DESTINATION bin)

35
ebos/alucartesianindexmapper.hh
View File

@ -30,6 +30,11 @@
#include <dune/grid/common/datahandleif.hh>
#include <dune/grid/utility/persistentcontainer.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#endif // HAVE_DUNE_ALUGRID
#include <opm/grid/common/CartesianIndexMapper.hpp>
#include <array>
#include <memory>
@ -37,16 +42,24 @@
#include <vector>
#include <cassert>
namespace Opm {
namespace Dune {
/*!
* \brief Interface class to access the logical Cartesian grid as used in industry
* standard simulator decks.
*/
template <class Grid>
class AluCartesianIndexMapper
//#if HAVE_MPI
template <>
class CartesianIndexMapper<Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>
{
public:
#if HAVE_MPI
using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
// data handle for communicating global ids during load balance and communication
template <class GridView>
class GlobalIndexDataHandle : public Dune::CommDataHandleIF<GlobalIndexDataHandle<GridView>, int>
@ -93,10 +106,10 @@ public:
~GlobalIndexDataHandle()
{ finalize(); }
bool contains(int dim, int codim) const
bool contains(int /* dim */, int codim) const
{ return codim == 0; }
bool fixedsize(int dim, int codim) const
bool fixedsize(int /* dim */, int /* codim */) const
{ return true; }
//! \brief loop over all internal data handlers and call gather for
@ -111,7 +124,7 @@ public:
//! \brief loop over all internal data handlers and call scatter for
//! given entity
template<class MessageBufferImp, class EntityType>
void scatter(MessageBufferImp& buff, const EntityType& element, size_t n)
void scatter(MessageBufferImp& buff, const EntityType& element, size_t /* n */)
{
int globalIdx = -1;
buff.read(globalIdx);
@ -125,7 +138,7 @@ public:
//! \brief loop over all internal data handlers and return sum of data
//! size of given entity
template<class EntityType>
size_t size(const EntityType& en) const
size_t size(const EntityType& /* en */) const
{ return 1; }
protected:
@ -164,9 +177,9 @@ public:
static const int dimension = Grid::dimension ;
/** \brief constructor taking grid */
AluCartesianIndexMapper(const Grid& grid,
const std::array<int, dimension>& cartDims,
const std::vector<int>& cartesianIndex)
CartesianIndexMapper(const Grid& grid,
const std::array<int, dimension>& cartDims,
const std::vector<int>& cartesianIndex)
: grid_(grid)
, cartesianDimensions_(cartDims)
, cartesianIndex_(cartesianIndex)
@ -248,6 +261,6 @@ protected:
const int cartesianSize_ ;
};
} // end namespace Opm
} // end namespace Dune
#endif

72
ebos/collecttoiorank.cc
View File

@ -35,7 +35,14 @@
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif // HAVE_DUNE_FEM
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <algorithm>
#include <cassert>
@ -817,7 +824,8 @@ CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid,
{
std::set<int> send, recv;
using EquilGridView = typename EquilGrid::LeafGridView;
typename std::is_same<Grid, EquilGrid>::type isSameGrid;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
ElementMapper elemMapper(localGridView, Dune::mcmgElementLayout());
sortedCartesianIdx_.reserve(localGridView.size(0));
@ -843,10 +851,12 @@ CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid,
using EquilElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView>;
EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
// Scatter the global index to local index for lookup during restart
ElementIndexScatterHandle<EquilElementMapper,ElementMapper> handle(equilElemMapper, elemMapper, localIdxToGlobalIdx_);
grid.scatterData(handle);
// Scatter the global index to local index for lookup during restart
if constexpr (isSameGrid) {
ElementIndexScatterHandle<EquilElementMapper,ElementMapper> handle(equilElemMapper, elemMapper, localIdxToGlobalIdx_);
grid.scatterData(handle);
}
// loop over all elements (global grid) and store Cartesian index
auto elemIt = equilGrid->leafGridView().template begin<0>();
const auto& elemEndIt = equilGrid->leafGridView().template end<0>();
@ -869,15 +879,19 @@ CollectDataToIORank(const Grid& grid, const EquilGrid* equilGrid,
// Scatter the global index to local index for lookup during restart
// This is a bit hacky since the type differs from the iorank.
// But should work since we only receive, i.e. use the second parameter.
ElementIndexScatterHandle<ElementMapper, ElementMapper> handle(elemMapper, elemMapper,
localIdxToGlobalIdx_);
grid.scatterData(handle);
if constexpr (isSameGrid) {
ElementIndexScatterHandle<ElementMapper, ElementMapper> handle(elemMapper, elemMapper, localIdxToGlobalIdx_);
grid.scatterData(handle);
}
}
// Sync the global element indices
ElementIndexHandle<ElementMapper> handle(elemMapper, localIdxToGlobalIdx_);
grid.communicate(handle, Dune::InteriorBorder_All_Interface,
if constexpr (isSameGrid) {
ElementIndexHandle<ElementMapper> handle(elemMapper, localIdxToGlobalIdx_);
grid.communicate(handle, Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication);
}
localIndexMap_.clear();
const size_t gridSize = grid.size(0);
localIndexMap_.reserve(gridSize);
@ -1056,11 +1070,45 @@ template class CollectDataToIORank<Dune::CpGrid,
Dune::CpGrid,
Dune::PartitionIteratorType(4),
false> > >;
#ifdef HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class CollectDataToIORank<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>;
template class CollectDataToIORank<ALUGrid3CN,
Dune::CpGrid,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>>;
#endif //HAVE_DUNE_ALUGRID
#else // ! HAVE_DUNE_FEM
template class CollectDataToIORank<Dune::CpGrid,
Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>;
#endif
#ifdef HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class CollectDataToIORank<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>>;
#endif //HAVE_DUNE_ALUGRID
#endif // ! HAVE_DUNE_FEM
template class CollectDataToIORank<Dune::PolyhedralGrid<3,3,double>,
Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>;

5
ebos/collecttoiorank.hh
View File

@ -28,9 +28,7 @@
#include <opm/output/data/Solution.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/output/data/Groups.hpp>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
#include <opm/grid/common/p2pcommunicator.hh>
#include <ebos/eclinterregflows.hh>
@ -112,6 +110,9 @@ public:
{ return toIORankComm_.size() > 1; }
int localIdxToGlobalIdx(unsigned localIdx) const;
bool doesNeedReordering() const
{ return needsReordering;}
size_t numCells () const
{ return globalCartesianIndex_.size(); }

171
ebos/eclalugridvanguard.hh
View File

@ -28,11 +28,15 @@
#define EWOMS_ECL_ALU_GRID_VANGUARD_HH
#include "eclbasevanguard.hh"
#include "ecltransmissibility.hh"
#include "alucartesianindexmapper.hh"
#include <opm/models/common/multiphasebaseproperties.hh>
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/common/fromtogridfactory.hh>
#include <dune/alugrid/dgf.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/utils/ParallelEclipseState.hpp>
namespace Opm {
template <class TypeTag>
@ -55,7 +59,11 @@ struct Vanguard<TypeTag, TTag::EclAluGridVanguard> {
};
template<class TypeTag>
struct Grid<TypeTag, TTag::EclAluGridVanguard> {
using type = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>;
#if HAVE_MPI
using type = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using type = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
};
template<class TypeTag>
struct EquilGrid<TypeTag, TTag::EclAluGridVanguard> {
@ -79,33 +87,30 @@ class EclAluGridVanguard : public EclBaseVanguard<TypeTag>
friend class EclBaseVanguard<TypeTag>;
typedef EclBaseVanguard<TypeTag> ParentType;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
public:
using Grid = GetPropType<TypeTag, Properties::Grid>;
using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
private:
typedef Opm::AluCartesianIndexMapper<Grid> CartesianIndexMapper;
using GridView = GetPropType<TypeTag, Properties::GridView>;
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
typedef Dune::CartesianIndexMapper<EquilGrid> EquilCartesianIndexMapper;
using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper, CartesianIndexMapper, Scalar>;
typedef Dune::FromToGridFactory<Grid> Factory;
static const int dimension = Grid::dimension;
static const int dimensionworld = Grid::dimensionworld;
public:
EclAluGridVanguard(Simulator& simulator)
: EclBaseVanguard<TypeTag>(simulator)
{
this->callImplementationInit();
{
this->mpiRank = EclGenericVanguard::comm().rank();
this->callImplementationInit();
}
~EclAluGridVanguard()
{
delete equilCartesianIndexMapper_;
delete grid_;
delete equilGrid_;
}
~EclAluGridVanguard() = default;
/*!
* \brief Return a reference to the simulation grid.
@ -162,6 +167,50 @@ public:
grid().communicate(*dataHandle,
Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication );
if (grid().size(0))
{
globalTrans_.reset(new TransmissibilityType(this->eclState(),
this->gridView(),
this->cartesianIndexMapper(),
this->grid(),
this->cellCentroids(),
getPropValue<TypeTag,
Properties::EnableEnergy>(),
getPropValue<TypeTag,
Properties::EnableDiffusion>()));
// Re-ordering for ALUGrid
globalTrans_->update(false, [&](unsigned int i) { return gridEquilIdxToGridIdx(i);});
}
}
template<class DataHandle>
void scatterData(DataHandle& handle) const
{
// not existing for this type of grid yet
}
template<class DataHandle>
void gatherData(DataHandle& handle) const
{
// not existing for this type of grid yet
}
template<class DataHandle, class InterfaceType, class CommunicationDirection>
void communicate (DataHandle& data, InterfaceType iftype, CommunicationDirection dir) const
{
// not existing for this type of grid yet
}
/*!
* \brief Free the memory occupied by the global transmissibility object.
*
* After writing the initial solution, this array should not be necessary anymore.
*/
void releaseGlobalTransmissibilities()
{
globalTrans_.reset();
}
/*!
@ -187,8 +236,39 @@ public:
std::function<std::array<double,dimensionworld>(int)>
cellCentroids() const
{
return this->cellCentroids_(cartesianIndexMapper_.get());
return this->cellCentroids_(this->cartesianIndexMapper());
}
const TransmissibilityType& globalTransmissibility() const
{
assert( globalTrans_ != nullptr );
return *globalTrans_;
}
void releaseGlobalTransmissibility()
{
globalTrans_.reset();
}
const std::vector<int>& globalCell()
{
return cartesianCellId_;
}
std::vector<int> cellPartition() const
{
// not required for this type of grid yet
return {};
}
unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const {
return equilGridToGrid_[elemIndex];
}
unsigned int gridIdxToEquilGridIdx(unsigned int elemIndex) const {
return ordering_[elemIndex];
}
protected:
void createGrids_()
{
@ -200,29 +280,51 @@ protected:
/////
// create the EQUIL grid
/////
equilGrid_ = new EquilGrid();
equilGrid_->processEclipseFormat(&(this->eclState().getInputGrid()),
&(this->eclState()),
/*isPeriodic=*/false,
/*flipNormals=*/false,
/*clipZ=*/false);
const EclipseGrid* input_grid = nullptr;
std::vector<double> global_porv;
// At this stage the ParallelEclipseState instance is still in global
// view; on rank 0 we have undistributed data for the entire grid, on
// the other ranks the EclipseState is empty.
if (mpiRank == 0) {
// Processing grid
input_grid = &this->eclState().getInputGrid();
global_porv = this->eclState().fieldProps().porv(true);
OpmLog::info("\nProcessing grid");
}
cartesianCellId_ = equilGrid_->globalCell();
this->equilGrid_ = std::make_unique<Dune::CpGrid>(EclGenericVanguard::comm());
// Note: removed_cells is guaranteed to be empty on ranks other than 0.
auto removed_cells =
this->equilGrid_->processEclipseFormat(input_grid,
&this->eclState(),
/*isPeriodic=*/false,
/*flipNormals=*/false,
/*clipZ=*/false);
cartesianCellId_ = this->equilGrid_->globalCell();
for (unsigned i = 0; i < dimension; ++i)
cartesianDimension_[i] = equilGrid_->logicalCartesianSize()[i];
cartesianDimension_[i] = this->equilGrid_->logicalCartesianSize()[i];
equilCartesianIndexMapper_ = new EquilCartesianIndexMapper(*equilGrid_);
equilCartesianIndexMapper_ = std::make_unique<EquilCartesianIndexMapper>(*equilGrid_);
/////
// create the simulation grid
/////
Dune::FromToGridFactory<Grid> factory;
grid_ = factory.convert(*equilGrid_, cartesianCellId_);
cartesianIndexMapper_ =
std::make_unique<CartesianIndexMapper>(*grid_, cartesianDimension_,
cartesianCellId_);
factory_ = std::make_unique<Factory>();
grid_ = factory_->convert(*equilGrid_, cartesianCellId_, ordering_);
OpmLog::warning("Space Filling Curve Ordering is not yet supported: DISABLE_ALUGRID_SFC_ORDERING is enabled");
equilGridToGrid_.resize(ordering_.size());
for (size_t index = 0; index<ordering_.size(); ++index) {
equilGridToGrid_[ordering_[index]] = index;
}
cartesianIndexMapper_ = std::make_unique<CartesianIndexMapper>(*grid_, cartesianDimension_, cartesianCellId_);
this->updateGridView_();
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
this->updateCellThickness_();
}
void filterConnections_()
@ -230,12 +332,17 @@ protected:
// not handling the removal of completions for this type of grid yet.
}
Grid* grid_;
EquilGrid* equilGrid_;
std::unique_ptr<Grid> grid_;
std::unique_ptr<EquilGrid> equilGrid_;
std::vector<int> cartesianCellId_;
std::vector<unsigned int> ordering_;
std::vector<unsigned int> equilGridToGrid_;
std::array<int,dimension> cartesianDimension_;
std::unique_ptr<CartesianIndexMapper> cartesianIndexMapper_;
EquilCartesianIndexMapper* equilCartesianIndexMapper_;
std::unique_ptr<EquilCartesianIndexMapper> equilCartesianIndexMapper_;
std::unique_ptr<Factory> factory_;
std::unique_ptr<TransmissibilityType> globalTrans_;
int mpiRank;
};
} // namespace Opm

40
ebos/eclbasevanguard.hh
View File

@ -428,6 +428,10 @@ public:
return grid.comm().sum(local_cells);
}
void setupCartesianToCompressed_() {
this->updateCartesianToCompressedMapping_();
}
protected:
/*!
* \brief Get function to query cell centroids for a distributed grid.
@ -502,15 +506,14 @@ protected:
const auto num_aqu_cells = this->allAquiferCells();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
for(const auto& element : elements(this->gridView())) {
const unsigned int elemIdx = elemMapper.index(element);
cellCenterDepth_[elemIdx] = cellCenterDepth(element);
if (!num_aqu_cells.empty()) {
const unsigned int global_index = cartesianIndex(elemIdx);
const auto search = num_aqu_cells.find(global_index);
if (search != num_aqu_cells.end()) {
const unsigned int global_index = cartesianIndex(elemIdx);
const auto search = num_aqu_cells.find(global_index);
if (search != num_aqu_cells.end()) {
// updating the cell depth using aquifer cell depth
cellCenterDepth_[elemIdx] = search->second->depth;
}
@ -530,9 +533,9 @@ protected:
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 auto& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellThickness_[elemIdx] = asImp_().computeCellThickness(element);
cellThickness_[elemIdx] = computeCellThickness(element);
}
}
@ -551,6 +554,29 @@ private:
return zz/Scalar(corners);
}
Scalar computeCellThickness(const typename GridView::template Codim<0>::Entity& element) const
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
Scalar zz1 = 0.0;
Scalar zz2 = 0.0;
const Geometry& geometry = element.geometry();
// This code only works with CP-grid where the
// number of corners are 8 and
// also assumes that the first
// 4 corners are the top surface and
// the 4 next are the bottomn.
assert(geometry.corners() == 8);
for (int i=0; i < 4; ++i){
zz1 += geometry.corner(i)[zCoord];
zz2 += geometry.corner(i+4)[zCoord];
}
zz1 /=4;
zz2 /=4;
return zz2-zz1;
}
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }

16
ebos/eclcpgridvanguard.hh
View File

@ -87,13 +87,13 @@ class EclCpGridVanguard : public EclBaseVanguard<TypeTag>
public:
using Grid = GetPropType<TypeTag, Properties::Grid>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper, Scalar>;
using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper, CartesianIndexMapper, Scalar>;
static const int dimensionworld = Grid::dimensionworld;
private:
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
using Element = typename GridView::template Codim<0>::Entity;
public:
@ -149,6 +149,13 @@ public:
#endif
}
unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const {
return elemIndex;
}
unsigned int gridIdxToEquilGridIdx(unsigned int elemIndex) const {
return elemIndex;
}
/*!
* \brief Get function to query cell centroids for a distributed grid.
*
@ -162,6 +169,11 @@ public:
return this->cellCentroids_(this->cartesianIndexMapper());
}
const std::vector<int>& globalCell()
{
return this->grid().globalCell();
}
protected:
void createGrids_()
{

2
ebos/eclequilinitializer.hh
View File

@ -103,7 +103,7 @@ public:
EQUIL::DeckDependent::InitialStateComputer<TypeTag> initialState(materialLawManager,
eclState,
vanguard.gridView(),
vanguard,
vanguard.cartesianMapper(),
simulator.problem().gravity()[dimWorld - 1]);

2
ebos/eclgenericcpgridvanguard.cc
View File

@ -48,7 +48,7 @@
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif //HAVE_DUNE_FEM
#include <cassert>
#include <numeric>

43
ebos/eclgenericproblem.cc
View File

@ -34,12 +34,17 @@
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif // HAVE_DUNE_FEM
#include <boost/date_time.hpp>
@ -729,12 +734,46 @@ template class EclGenericProblem<Dune::Fem::GridPart2GridViewImpl<
Opm::BlackOilFluidSystem<
double,
Opm::BlackOilDefaultIndexTraits>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericProblem<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
template class EclGenericProblem<Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>,
Opm::BlackOilFluidSystem<
double,
Opm::BlackOilDefaultIndexTraits>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericProblem<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
#endif
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = const Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericProblem<Dune::GridView<Dune::ALU3dLeafGridViewTraits<ALUGrid3CN, Dune::PartitionIteratorType(4)>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif //HAVE_DUNE_FEM
template class EclGenericProblem<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,

48
ebos/eclgenericthresholdpressure.cc
View File

@ -37,12 +37,17 @@
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif // HAVE_DUNE_FEM
#include <algorithm>
#include <cassert>
@ -263,11 +268,52 @@ template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
double>;
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false> >,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericThresholdPressure<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericThresholdPressure<ALUGrid3CN,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif
template class EclGenericThresholdPressure<Dune::PolyhedralGrid<3,3,double>,

46
ebos/eclgenerictracermodel.cc
View File

@ -45,7 +45,13 @@
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif // HAVE_DUNE_FEM
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <fmt/format.h>
#include <iostream>
@ -395,13 +401,49 @@ template class EclGenericTracerModel<Dune::CpGrid,
Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >,
false, false>,
double>;
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericTracerModel<ALUGrid3CN, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>, Opm::EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,false,false>,
double>;
template class EclGenericTracerModel<ALUGrid3CN,
Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> > >,
Opm::EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
false, false>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else
template class EclGenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>,
double>;
#endif
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericTracerModel<ALUGrid3CN,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
Dune::PartitionIteratorType(4)>>>,
Opm::EcfvStencil<double,Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
Dune::PartitionIteratorType(4)>>,false,false>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif //HAVE_DUNE_FEM
template class EclGenericTracerModel<Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,

89
ebos/eclgenericwriter.cc
View File

@ -28,6 +28,11 @@
#include <opm/grid/cpgrid/GridHelpers.hpp>
#include <opm/grid/polyhedralgrid.hh>
#include <opm/grid/utility/cartesianToCompressed.hpp>
#if HAVE_DUNE_ALUGRID
#include "eclalugridvanguard.hh"
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#endif // HAVE_DUNE_ALUGRID
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
@ -47,7 +52,7 @@
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif
#endif // HAVE_DUNE_FEM
#if HAVE_MPI
#include <mpi.h>
@ -242,21 +247,20 @@ eclIO() const
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
writeInit()
writeInit(const std::function<unsigned int(unsigned int)>& map)
{
if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int> > integerVectors;
if (collectToIORank_.isParallel())
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
auto cartMap = cartesianToCompressed(equilGrid_->size(0),
UgGridHelpers::globalCell(*equilGrid_));
eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap));
auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
eclIO_->writeInitial(computeTrans_(cartMap, map), integerVectors, exportNncStructure_(cartMap, map));
}
}
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const
computeTrans_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
{
const auto& cartMapper = *equilCartMapper_;
const auto& cartDims = cartMapper.cartesianDimensions();
@ -304,6 +308,12 @@ computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const
int gc1 = std::min(cartIdx1, cartIdx2);
int gc2 = std::max(cartIdx1, cartIdx2);
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
tranx.data[gc1] = globalTrans().transmissibility(c1, c2);
continue; // skip other if clauses as they are false, last one needs some computation
@ -327,7 +337,7 @@ computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const
template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const
exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
{
std::size_t nx = eclState_.getInputGrid().getNX();
std::size_t ny = eclState_.getInputGrid().getNY();
@ -378,14 +388,20 @@ exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const
if (c1 > c2)
continue; // we only need to handle each connection once, thank you.
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1];
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2];
std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
if ( cc2 < cc1 )
std::swap(cc1, cc2);
auto cellDiff = cc2 - cc1;
// Re-ordering in case of non-empty mapping between equilGrid to grid
if (map) {
c1 = map(c1); // equilGridToGrid map
c2 = map(c2);
}
if (cellDiff != 1 &&
cellDiff != nx &&
cellDiff != nx*ny &&
@ -430,9 +446,10 @@ doWriteOutput(const int reportStepNum,
bool doublePrecision)
{
const auto isParallel = this->collectToIORank_.isParallel();
const bool needsReordering = this->collectToIORank_.doesNeedReordering();
RestartValue restartValue {
isParallel ? this->collectToIORank_.globalCellData()
(isParallel || needsReordering) ? this->collectToIORank_.globalCellData()
: std::move(localCellData),
isParallel ? this->collectToIORank_.globalWellData()
@ -561,8 +578,7 @@ globalTrans() const
#if HAVE_DUNE_FEM
template class EclGenericWriter<Dune::CpGrid,
Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
double>;
template class EclGenericWriter<Dune::CpGrid,
Dune::CpGrid,
@ -578,18 +594,59 @@ template class EclGenericWriter<Dune::CpGrid,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#ifdef HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
double>;
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false>>>,
double>;
#endif // HAVE_DUNE_ALUGRID
#else // !HAVE_DUNE_FEM
template class EclGenericWriter<Dune::CpGrid,
Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
double>;
#endif
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclGenericWriter<ALUGrid3CN,
Dune::CpGrid,
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>>,
double>;
#endif // HAVE_DUNE_ALUGRID
#endif // !HAVE_DUNE_FEM
template class EclGenericWriter<Dune::PolyhedralGrid<3,3,double>,
Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
double>;
} // namespace Opm

21
ebos/eclgenericwriter.hh
View File

@ -30,7 +30,10 @@
#include "collecttoiorank.hh"
#include <ebos/ecltransmissibility.hh>
#if HAVE_DUNE_ALUGRID
#include "alucartesianindexmapper.hh"
#include <dune/alugrid/common/fromtogridfactory.hh>
#endif // HAVE_DUNE_ALUGRID
#include <opm/models/parallel/tasklets.hh>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
@ -62,9 +65,11 @@ namespace Opm {
template <class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
class EclGenericWriter
{
{
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
typedef Dune::CartesianIndexMapper<EquilGrid> EquilCartesianIndexMapper;
using CollectDataToIORankType = CollectDataToIORank<Grid,EquilGrid,GridView>;
using TransmissibilityType = EclTransmissibility<Grid,GridView,ElementMapper,Scalar>;
using TransmissibilityType = EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>;
public:
// The Simulator object should preferably have been const - the
@ -83,7 +88,7 @@ public:
const EclipseIO& eclIO() const;
void writeInit();
void writeInit(const std::function<unsigned int(unsigned int)>& map);
void setTransmissibilities(const TransmissibilityType* globalTrans)
{
@ -98,10 +103,10 @@ public:
{
simulation_report_ = report;
}
protected:
const TransmissibilityType& globalTrans() const;
unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const;
void doWriteOutput(const int reportStepNum,
const bool isSubStep,
@ -151,8 +156,8 @@ protected:
SimulatorReport simulation_report_;
private:
data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const;
std::vector<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const;
data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const;
std::vector<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const;
};
} // namespace Opm

11
ebos/ecloutputblackoilmodule.hh
View File

@ -417,7 +417,7 @@ public:
this->bubblePointPressure_[globalDofIdx] = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
}
catch (const NumericalIssue&) {
const auto cartesianIdx = elemCtx.simulator().vanguard().grid().globalCell()[globalDofIdx];
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
this->failedCellsPb_.push_back(cartesianIdx);
}
}
@ -426,7 +426,7 @@ public:
this->dewPointPressure_[globalDofIdx] = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
}
catch (const NumericalIssue&) {
const auto cartesianIdx = elemCtx.simulator().vanguard().grid().globalCell()[globalDofIdx];
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
this->failedCellsPd_.push_back(cartesianIdx);
}
}
@ -528,10 +528,11 @@ public:
updateFluidInPlace_(elemCtx, dofIdx);
// Adding block data
const auto cartesianIdx = elemCtx.simulator().vanguard().grid().globalCell()[globalDofIdx];
for (auto& val : this->blockData_) {
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
for (auto& val: this->blockData_) {
const auto& key = val.first;
int cartesianIdxBlock = key.second - 1;
assert(key.second > 0);
unsigned int cartesianIdxBlock = key.second - 1;
if (cartesianIdx == cartesianIdxBlock) {
if ((key.first == "BWSAT") || (key.first == "BSWAT"))
val.second = getValue(fs.saturation(waterPhaseIdx));

1
ebos/eclpolyhedralgridvanguard.hh
View File

@ -195,6 +195,7 @@ public:
{
return this->cellCentroids_(this->cartesianIndexMapper());
}
protected:
void createGrids_()
{

98
ebos/eclproblem.hh
View File

@ -28,24 +28,19 @@
#ifndef EWOMS_ECL_PROBLEM_HH
#define EWOMS_ECL_PROBLEM_HH
//#define DISABLE_ALUGRID_SFC_ORDERING 1
//#define EBOS_USE_ALUGRID 1
// make sure that the EBOS_USE_ALUGRID macro. using the preprocessor for this is slightly
// hacky...
#if EBOS_USE_ALUGRID
//#define DISABLE_ALUGRID_SFC_ORDERING 1
#if USE_ALUGRID
#define DISABLE_ALUGRID_SFC_ORDERING 1
#if !HAVE_DUNE_ALUGRID
#warning "ALUGrid was indicated to be used for the ECL black oil simulator, but this "
#warning "requires the presence of dune-alugrid >= 2.4. Falling back to Dune::CpGrid"
#undef EBOS_USE_ALUGRID
#define EBOS_USE_ALUGRID 0
#undef USE_ALUGRID
#define USE_ALUGRID 0
#endif
#else
#define EBOS_USE_ALUGRID 0
#define USE_ALUGRID 0
#endif
#if EBOS_USE_ALUGRID
#if USE_ALUGRID
#include "eclalugridvanguard.hh"
#elif USE_POLYHEDRALGRID
#include "eclpolyhedralgridvanguard.hh"
@ -110,6 +105,7 @@
#include <vector>
#include <string>
#include <algorithm>
#include<functional>
namespace Opm {
template <class TypeTag>
@ -120,7 +116,7 @@ namespace Opm::Properties {
namespace TTag {
#if EBOS_USE_ALUGRID
#if USE_ALUGRID
struct EclBaseProblem {
using InheritsFrom = std::tuple<VtkEclTracer, EclOutputBlackOil, EclAluGridVanguard>;
};
@ -867,7 +863,19 @@ public:
this->readRockParameters_(simulator.vanguard().cellCenterDepths());
readMaterialParameters_();
readThermalParameters_();
transmissibilities_.finishInit();
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> gridToEquilGrid;
#ifdef HAVE_DUNE_ALUGRID
using Grid = GetPropType<TypeTag, Properties::Grid>;
typename std::is_same<Grid, Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>::type isAlugrid;
if (isAlugrid) {
gridToEquilGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridIdxToEquilGridIdx(i);
};
}
#endif // HAVE_DUNE_ALUGRID
transmissibilities_.finishInit(gridToEquilGrid);
const auto& initconfig = eclState.getInitConfig();
tracerModel_.init(initconfig.restartRequested());
@ -937,7 +945,16 @@ public:
} else
eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
eclWriter_->writeInit();
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> equilGridToGrid;
#ifdef HAVE_DUNE_ALUGRID
if (isAlugrid) {
equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
}
#endif // HAVE_DUNE_ALUGRID
eclWriter_->writeInit(equilGridToGrid);
}
simulator.vanguard().releaseGlobalTransmissibilities();
@ -1023,8 +1040,21 @@ public:
eclState.apply_schedule_keywords( miniDeck );
eclBroadcast(cc, eclState.getTransMult() );
// Re-ordering in case of ALUGrid
std::function<unsigned int(unsigned int)> equilGridToGrid;
#ifdef HAVE_DUNE_ALUGRID
using Grid = GetPropType<TypeTag, Properties::Grid>;
typename std::is_same<Grid, Dune::ALUGrid<3, 3, Dune::cube,
Dune::nonconforming>>::type isAlugrid;
if (isAlugrid) {
equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
}
#endif // HAVE_DUNE_ALUGRID
// re-compute all quantities which may possibly be affected.
transmissibilities_.update(true);
transmissibilities_.update(true, equilGridToGrid);
this->referencePorosity_[1] = this->referencePorosity_[0];
updateReferencePorosity_();
updatePffDofData_();
@ -1161,13 +1191,27 @@ public:
auto& schedule = simulator.vanguard().schedule();
auto& ecl_state = simulator.vanguard().eclState();
int episodeIdx = this->episodeIndex();
// Re-ordering in case of Alugrid
std::function<unsigned int(unsigned int)> gridToEquilGrid;
#ifdef HAVE_DUNE_ALUGRID
using Grid = GetPropType<TypeTag, Properties::Grid>;
typename std::is_same<Grid, Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>::type isAlugrid;
if (isAlugrid) {
gridToEquilGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridIdxToEquilGridIdx(i);
};
}
#endif // HAVE_DUNE_ALUGRID
this->applyActions(episodeIdx,
simulator.time() + simulator.timeStepSize(),
simulator.vanguard().grid().comm(),
ecl_state,
schedule,
simulator.vanguard().actionState(),
simulator.vanguard().summaryState());
simulator.vanguard().summaryState(),
gridToEquilGrid);
// deal with "clogging" for the MICP model
if constexpr (enableMICP){
@ -1279,7 +1323,7 @@ public:
the simulator properties which need to be updated. This functionality is
probably not complete.
*/
void applySimulatorUpdate(int report_step, Parallel::Communication comm, const SimulatorUpdate& sim_update, EclipseState& ecl_state, Schedule& schedule, SummaryState& summary_state, bool& commit_wellstate) {
void applySimulatorUpdate(int report_step, Parallel::Communication comm, const SimulatorUpdate& sim_update, EclipseState& ecl_state, Schedule& schedule, SummaryState& summary_state, bool& commit_wellstate, const std::function<unsigned int(unsigned int)>& map = {}) {
this->wellModel_.updateEclWells(report_step, sim_update.affected_wells, summary_state);
if (!sim_update.affected_wells.empty())
commit_wellstate = true;
@ -1290,7 +1334,7 @@ public:
eclBroadcast(comm, ecl_state.getTransMult() );
// re-compute transmissibility
transmissibilities_.update(true);
transmissibilities_.update(true,map);
}
}
@ -1302,7 +1346,8 @@ public:
EclipseState& ecl_state,
Schedule& schedule,
Action::State& actionState,
SummaryState& summaryState) {
SummaryState& summaryState,
const std::function<unsigned int(unsigned int)>& map = {}) {
const auto& actions = schedule[reportStep].actions();
if (actions.empty())
return;
@ -1322,7 +1367,7 @@ public:
bool commit_wellstate = false;
for (const auto& pyaction : actions.pending_python(actionState)) {
auto sim_update = schedule.runPyAction(reportStep, *pyaction, actionState, ecl_state, summaryState);
this->applySimulatorUpdate(reportStep, comm, sim_update, ecl_state, schedule, summaryState, commit_wellstate);
this->applySimulatorUpdate(reportStep, comm, sim_update, ecl_state, schedule, summaryState, commit_wellstate, map);
}
auto simTime = asTimeT(now);
@ -1342,7 +1387,7 @@ public:
const auto& wellpi = this->fetchWellPI(reportStep, *action, schedule, matching_wells);
auto sim_update = schedule.applyAction(reportStep, *action, actionResult.wells(), wellpi);
this->applySimulatorUpdate(reportStep, comm, sim_update, ecl_state, schedule, summaryState, commit_wellstate);
this->applySimulatorUpdate(reportStep, comm, sim_update, ecl_state, schedule, summaryState, commit_wellstate, map);
actionState.add_run(*action, simTime, std::move(actionResult));
} else {
std::string msg = "The action: " + action->name() + " evaluated to false at " + ts;
@ -1502,6 +1547,17 @@ public:
Scalar dofCenterDepth(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
{
unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
return dofCenterDepth(globalSpaceIdx);
}
/*!
* \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 dofCenterDepth(unsigned globalSpaceIdx) const
{
return this->simulator().vanguard().cellCenterDepth(globalSpaceIdx);
}

211
ebos/ecltransmissibility.cc
View File

@ -30,6 +30,12 @@
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
@ -82,11 +88,11 @@ std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
namespace Opm {
template<class Grid, class GridView, class ElementMapper, class Scalar>
EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
EclTransmissibility(const EclipseState& eclState,
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const CartesianIndexMapper& cartMapper,
const Grid& grid,
std::function<std::array<double,dimWorld>(int)> centroids,
bool enableEnergy,
@ -103,36 +109,36 @@ EclTransmissibility(const EclipseState& eclState,
transmissibilityThreshold_ = unitSystem.parse("Transmissibility").getSIScaling() * 1e-6;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
{
return trans_.at(isId(elemIdx1, elemIdx2));
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
transmissibilityBoundary(unsigned elemIdx, unsigned boundaryFaceIdx) const
{
return transBoundary_.at(std::make_pair(elemIdx, boundaryFaceIdx));
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
{
return thermalHalfTrans_.at(directionalIsId(insideElemIdx, outsideElemIdx));
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
{
return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx));
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
{
if (diffusivity_.empty())
@ -142,9 +148,9 @@ diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
update(bool global)
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
update(bool global, const std::function<unsigned int(unsigned int)>& map)
{
const auto& cartDims = cartMapper_.cartesianDimensions();
auto& transMult = eclState_.getTransMult();
@ -155,8 +161,11 @@ update(bool global)
const std::vector<double>& ntg = eclState_.fieldProps().get_double("NTG");
const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_;
unsigned numElements = elemMapper.size();
extractPermeability_();
if (map)
extractPermeability_(map);
else
extractPermeability_();
// calculate the axis specific centroids of all elements
std::array<std::vector<DimVector>, dimWorld> axisCentroids;
@ -481,8 +490,8 @@ update(bool global)
removeSmallNonCartesianTransmissibilities_();
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
extractPermeability_()
{
unsigned numElem = gridView_.size(/*codim=*/0);
@ -523,8 +532,51 @@ extractPermeability_()
"(The PERM{X,Y,Z} keywords are missing)");
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
extractPermeability_(const std::function<unsigned int(unsigned int)>& map)
{
unsigned numElem = gridView_.size(/*codim=*/0);
permeability_.resize(numElem);
// read the intrinsic permeabilities from the eclState. Note that all arrays
// provided by eclState are one-per-cell of "uncompressed" grid, whereas the
// simulation grid might remove a few elements. (e.g. because it is distributed
// over several processes.)
const auto& fp = eclState_.fieldProps();
if (fp.has_double("PERMX")) {
const std::vector<double>& permxData = fp.get_double("PERMX");
std::vector<double> permyData;
if (fp.has_double("PERMY"))
permyData = fp.get_double("PERMY");
else
permyData = permxData;
std::vector<double> permzData;
if (fp.has_double("PERMZ"))
permzData = fp.get_double("PERMZ");
else
permzData = permxData;
for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
permeability_[dofIdx] = 0.0;
size_t inputDofIdx = map(dofIdx);
permeability_[dofIdx][0][0] = permxData[inputDofIdx];
permeability_[dofIdx][1][1] = permyData[inputDofIdx];
permeability_[dofIdx][2][2] = permzData[inputDofIdx];
}
// for now we don't care about non-diagonal entries
}
else
throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
"(The PERM{X,Y,Z} keywords are missing)");
}
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
extractPorosity_()
{
// read the intrinsic porosity from the eclState. Note that all arrays
@ -540,8 +592,8 @@ extractPorosity_()
"(The PORO keywords are missing)");
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
removeSmallNonCartesianTransmissibilities_()
{
const auto& cartDims = cartMapper_.cartesianDimensions();
@ -562,8 +614,8 @@ removeSmallNonCartesianTransmissibilities_()
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper, Scalar>::
applyAllZMultipliers_(Scalar& trans,
unsigned insideFaceIdx,
unsigned outsideFaceIdx,
@ -601,8 +653,8 @@ applyAllZMultipliers_(Scalar& trans,
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
updateFromEclState_(bool global)
{
const FieldPropsManager* fp =
@ -633,9 +685,9 @@ updateFromEclState_(bool global)
resetTransmissibilityFromArrays_(is_tran, trans);
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
std::array<std::vector<double>,3>
EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
createTransmissibilityArrays_(const std::array<bool,3>& is_tran)
{
const auto& cartDims = cartMapper_.cartesianDimensions();
@ -700,8 +752,8 @@ createTransmissibilityArrays_(const std::array<bool,3>& is_tran)
return trans;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
resetTransmissibilityFromArrays_(const std::array<bool,3>& is_tran,
const std::array<std::vector<double>,3>& trans)
{
@ -758,9 +810,9 @@ resetTransmissibilityFromArrays_(const std::array<bool,3>& is_tran,
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
template<class Intersection>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
computeFaceProperties(const Intersection& intersection,
const int,
const int,
@ -781,9 +833,9 @@ computeFaceProperties(const Intersection& intersection,
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
template<class Intersection>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
computeFaceProperties(const Intersection& intersection,
const int insideElemIdx,
const int insideFaceIdx,
@ -800,9 +852,9 @@ computeFaceProperties(const Intersection& intersection,
faceAreaNormal = grid_.faceAreaNormalEcl(faceIdx);
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
std::tuple<std::vector<NNCdata>, std::vector<NNCdata>>
EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompressed)
{
// First scale NNCs with EDITNNC.
@ -853,8 +905,8 @@ applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompr
return std::make_tuple(processedNnc, unprocessedNnc);
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLocal)
{
const auto& nnc_input = eclState_.getInputNNC();
@ -926,8 +978,8 @@ applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLoca
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
computeHalfTrans_(Scalar& halfTrans,
const DimVector& areaNormal,
int faceIdx, // in the reference element that contains the intersection
@ -947,8 +999,8 @@ computeHalfTrans_(Scalar& halfTrans,
halfTrans /= distance.two_norm2();
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
computeHalfDiffusivity_(Scalar& halfDiff,
const DimVector& areaNormal,
const DimVector& distance,
@ -963,9 +1015,9 @@ computeHalfDiffusivity_(Scalar& halfDiff,
halfDiff /= distance.two_norm2();
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
typename EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::DimVector
EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
typename EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::DimVector
EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
distanceVector_(const DimVector& center,
int faceIdx, // in the reference element that contains the intersection
unsigned elemIdx,
@ -980,8 +1032,8 @@ distanceVector_(const DimVector& center,
return x;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
applyMultipliers_(Scalar& trans,
unsigned faceIdx,
unsigned cartElemIdx,
@ -1014,8 +1066,8 @@ applyMultipliers_(Scalar& trans,
}
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,Scalar>::
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
applyNtg_(Scalar& trans,
unsigned faceIdx,
unsigned elemIdx,
@ -1047,6 +1099,7 @@ applyNtg_(Scalar& trans,
template class EclTransmissibility<Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
Dune::CartesianIndexMapper<Dune::CpGrid>,
double>;
template class EclTransmissibility<Dune::CpGrid,
Dune::Fem::GridPart2GridViewImpl<
@ -1060,17 +1113,65 @@ template class EclTransmissibility<Dune::CpGrid,
Dune::CpGrid,
Dune::PartitionIteratorType(4),
false> > >,
Dune::CartesianIndexMapper<Dune::CpGrid>,
double>;
#else
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class EclTransmissibility<ALUGrid3CN, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
Dune::CartesianIndexMapper<ALUGrid3CN>,
double>;
template class EclTransmissibility<ALUGrid3CN,
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false> >,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::Fem::GridPart2GridViewImpl<
Dune::Fem::AdaptiveLeafGridPart<
ALUGrid3CN,
Dune::PartitionIteratorType(4),
false> > >,
Dune::CartesianIndexMapper<ALUGrid3CN>,
double>;
#endif //HAVE_DUNE_ALUGRID
#else // !DUNE_FEM
template class EclTransmissibility<Dune::CpGrid,
Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
Dune::CartesianIndexMapper<Dune::CpGrid>,
double>;
#endif
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
#else
using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
#endif //HAVE_MPI
template class
EclTransmissibility<ALUGrid3CN,
Dune::GridView<
Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>
>,
Dune::MultipleCodimMultipleGeomTypeMapper<
Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>
>,
Dune::CartesianIndexMapper<ALUGrid3CN>,
double>;
#endif //HAVE_DUNE_ALUGRID
#endif //HAVE_DUNE_FEM
template class EclTransmissibility<Dune::PolyhedralGrid<3,3,double>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>,
Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
Dune::CartesianIndexMapper<Dune::PolyhedralGrid<3,3,double>>,
double>;
} // namespace Opm

22
ebos/ecltransmissibility.hh
View File

@ -33,6 +33,12 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#include "alucartesianindexmapper.hh"
#endif // HAVE_DUNE_ALUGRID
#include <array>
#include <map>
#include <tuple>
@ -46,7 +52,7 @@ class EclipseState;
struct NNCdata;
class TransMult;
template<class Grid, class GridView, class ElementMapper, class Scalar>
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
class EclTransmissibility {
// Grid and world dimension
enum { dimWorld = GridView::dimensionworld };
@ -57,7 +63,7 @@ public:
EclTransmissibility(const EclipseState& eclState,
const GridView& gridView,
const Dune::CartesianIndexMapper<Grid>& cartMapper,
const CartesianIndexMapper& cartMapper,
const Grid& grid,
std::function<std::array<double,dimWorld>(int)> centroids,
bool enableEnergy,
@ -116,8 +122,8 @@ public:
* permeabilities. This approach is probably not always correct
* either but at least it seems to be much better.
*/
void finishInit()
{ update(true); }
void finishInit(const std::function<unsigned int(unsigned int)>& map = {})
{ update(true,map); }
/*!
* \brief Compute all transmissibilities
@ -125,7 +131,7 @@ public:
* \param global If true, update is called on all processes
* Also, this updates the "thermal half transmissibilities" if energy is enabled.
*/
void update(bool global);
void update(bool global, const std::function<unsigned int(unsigned int)>& map = {});
protected:
void updateFromEclState_(bool global);
@ -201,7 +207,9 @@ protected:
void applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLocal);
void extractPermeability_();
void extractPermeability_(const std::function<unsigned int(unsigned int)>& map);
void extractPorosity_();
void computeHalfTrans_(Scalar& halfTrans,
@ -235,7 +243,7 @@ protected:
std::unordered_map<std::uint64_t, Scalar> trans_;
const EclipseState& eclState_;
const GridView& gridView_;
const Dune::CartesianIndexMapper<Grid>& cartMapper_;
const CartesianIndexMapper& cartMapper_;
const Grid& grid_;
std::function<std::array<double,dimWorld>(int)> centroids_;
Scalar transmissibilityThreshold_;

6
ebos/eclwriter.hh
View File

@ -300,8 +300,8 @@ public:
// add cell data to perforations for Rft output
this->eclOutputModule_->addRftDataToWells(localWellData, reportStepNum);
}
if (this->collectToIORank_.isParallel()) {
if (this->collectToIORank_.isParallel()|| this->collectToIORank_.doesNeedReordering()) {
this->collectToIORank_.collect(localCellData,
eclOutputModule_->getBlockData(),
eclOutputModule_->getWBPData(),
@ -461,7 +461,7 @@ private:
eclOutputModule_->processElement(elemCtx);
}
OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ", simulator_.vanguard().grid().comm())
OPM_END_PARALLEL_TRY_CATCH("EclWriter::prepareLocalCellData() failed: ", simulator_.vanguard().grid().comm());
}
void captureLocalFluxData()

39
ebos/equil/initstateequil.hh
View File

@ -1588,26 +1588,27 @@ class InitialStateComputer
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>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
public:
template<class MaterialLawManager>
InitialStateComputer(MaterialLawManager& materialLawManager,
const EclipseState& eclipseState,
const GridView& gridView,
const Opm::EclipseState& eclipseState,
const Vanguard& vanguard,
const CartesianIndexMapper& cartMapper,
const double grav = unit::gravity,
const double grav = Opm::unit::gravity,
const bool applySwatInit = true)
: temperature_(gridView.size(/*codim=*/0)),
saltConcentration_(gridView.size(/*codim=*/0)),
saltSaturation_(gridView.size(/*codim=*/0)),
: temperature_(vanguard.grid().size(/*codim=*/0)),
saltConcentration_(vanguard.grid().size(/*codim=*/0)),
saltSaturation_(vanguard.grid().size(/*codim=*/0)),
pp_(FluidSystem::numPhases,
std::vector<double>(gridView.size(/*codim=*/0))),
std::vector<double>(vanguard.grid().size(/*codim=*/0))),
sat_(FluidSystem::numPhases,
std::vector<double>(gridView.size(/*codim=*/0))),
rs_(gridView.size(/*codim=*/0)),
rv_(gridView.size(/*codim=*/0)),
std::vector<double>(vanguard.grid().size(/*codim=*/0))),
rs_(vanguard.grid().size(/*codim=*/0)),
rv_(vanguard.grid().size(/*codim=*/0)),
cartesianIndexMapper_(cartMapper)
{
//Check for presence of kw SWATINIT
@ -1620,13 +1621,14 @@ public:
// Querry cell depth, cell top-bottom.
// numerical aquifer cells might be specified with different depths.
const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
updateCellProps_(gridView, num_aquifers);
updateCellProps_(vanguard, num_aquifers);
const Grid& grid = vanguard.grid();
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(eclipseState);
const auto& tables = eclipseState.getTableManager();
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(eclipseState, gridView));
const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
const int invalidRegion = -1;
regionPvtIdx_.resize(rec.size(), invalidRegion);
setRegionPvtIdx(eclipseState, eqlmap);
@ -1737,11 +1739,11 @@ public:
updateInitialSaltSaturation_(eclipseState, eqlmap);
// Compute pressures, saturations, rs and rv factors.
const auto& comm = gridView.comm();
const auto& comm = grid.comm();
calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
// modify the pressure and saturation for numerical aquifer cells
applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage());
applyNumericalAquifers_(vanguard, num_aquifers, eclipseState.runspec().co2Storage());
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
@ -1835,10 +1837,12 @@ private:
std::vector<std::pair<double,double>> cellZSpan_;
std::vector<std::pair<double,double>> cellZMinMax_;
void updateCellProps_(const GridView& gridView,
void updateCellProps_(const Vanguard& vanguard,
const NumericalAquifers& aquifer)
{
const auto& gridView = vanguard.gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
//const auto& elemMapper = gridView.elementMapper();
int numElements = gridView.size(/*codim=*/0);
cellCenterDepth_.resize(numElements);
cellZSpan_.resize(numElements);
@ -1869,14 +1873,15 @@ private:
}
}
void applyNumericalAquifers_(const GridView& gridView,
void applyNumericalAquifers_(const Vanguard& vanguard,
const NumericalAquifers& aquifer,
const bool co2store)
{
const auto num_aqu_cells = aquifer.allAquiferCells();
if (num_aqu_cells.empty()) return;
const auto& gridView = vanguard.gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
//const auto& elemMapper = vanguard.gridView().elementMapper();;
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {

47
flow/flow_blackoil_alugrid.cpp Normal file
View File

@ -0,0 +1,47 @@
/*
Copyright 2020, NORCE AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/simulators/flow/Main.hpp>
namespace Opm::Properties {
namespace TTag {
struct EclFlowProblemAlugrid {
using InheritsFrom = std::tuple<EclFlowProblem>;
};
}
// by default use the dummy aquifer "model"
template<class TypeTag>
struct EclAquiferModel<TypeTag, TTag::EclFlowProblemAlugrid> {
using type = Opm::EclBaseAquiferModel<TypeTag>;
};
// Enable aquifers by default in experimental mode
template<class TypeTag>
struct EclEnableAquifers<TypeTag, TTag::EclFlowProblemAlugrid> {
static constexpr bool value = false;
};
}
int main(int argc, char** argv)
{
using TypeTag = Opm::Properties::TTag::EclFlowProblemAlugrid;
auto mainObject = Opm::Main(argc, argv);
return mainObject.runStatic<TypeTag>();
}

12
opm/core/props/satfunc/RelpermDiagnostics.cpp
View File

@ -38,7 +38,10 @@
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
#ifdef HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/3d/gridview.hh>
#endif // HAVE_DUNE_ALUGRID
namespace Opm{
bool RelpermDiagnostics::phaseCheck_(const EclipseState& es)
@ -809,4 +812,11 @@ namespace Opm{
INSTANCE_DIAGNOSIS(Dune::CpGrid)
INSTANCE_DIAGNOSIS(Dune::PolyhedralGrid<3,3>)
#if HAVE_DUNE_ALUGRID
#if HAVE_MPI
INSTANCE_DIAGNOSIS(Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>)
#else
INSTANCE_DIAGNOSIS(Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>)
#endif //HAVE_MPI
#endif //HAVE_DUNE_ALUGRID
} //namespace Opm

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

@ -187,7 +187,6 @@ BlackoilAquiferModel<TypeTag>::init()
return;
}
const auto& connections = aquifer.connections();
for (const auto& aq : aquifer.ct()) {
aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
@ -201,9 +200,8 @@ BlackoilAquiferModel<TypeTag>::init()
if (aquifer.hasNumericalAquifer()) {
const auto& num_aquifers = aquifer.numericalAquifers().aquifers();
const auto& ugrid = simulator_.vanguard().grid();
const int number_of_cells = simulator_.gridView().size(0);
const int* global_cell = UgGridHelpers::globalCell(ugrid);
const int* global_cell = this->simulator_.vanguard().globalCell().data();
const std::unordered_map<int, int> cartesian_to_compressed = cartesianToCompressed(number_of_cells,
global_cell);
for ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {

4
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@ -229,7 +229,7 @@ namespace Opm {
, well_model_ (well_model)
, terminal_output_ (terminal_output)
, current_relaxation_(1.0)
, dx_old_(UgGridHelpers::numCells(grid_))
, dx_old_(ebosSimulator_.model().numGridDof())
{
// compute global sum of number of cells
global_nc_ = detail::countGlobalCells(grid_);
@ -350,7 +350,7 @@ namespace Opm {
//residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ;
// Compute the nonlinear update.
const int nc = UgGridHelpers::numCells(grid_);
unsigned nc = ebosSimulator_.model().numGridDof();
BVector x(nc);
// Solve the linear system.

2
opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -80,6 +80,7 @@ public:
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using AquiferModel = GetPropType<TypeTag, Properties::EclAquiferModel>;
typedef AdaptiveTimeSteppingEbos<TypeTag> TimeStepper;
typedef BlackOilPolymerModule<TypeTag> PolymerModule;
@ -90,7 +91,6 @@ public:
typedef typename Model::ModelParameters ModelParameters;
typedef typename Solver::SolverParameters SolverParameters;
typedef BlackoilWellModel<TypeTag> WellModel;
typedef BlackoilAquiferModel<TypeTag> AquiferModel;
/// Initialise from parameters and objects to observe.

5
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -159,7 +159,7 @@ namespace Opm
extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_);
// For some reason simulator_.model().elementMapper() is not initialized at this stage
// Hence const auto& elemMapper = simulator_.model().elementMapper(); does not work.
//const auto& elemMapper = simulator_.model().elementMapper(); //does not work.
// Set it up manually
ElementMapper elemMapper(simulator_.vanguard().gridView(), Dune::mcmgElementLayout());
detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
@ -224,7 +224,8 @@ namespace Opm
useWellConn_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
if (numJacobiBlocks_ > 1) {
const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd();
detail::setWellConnections(simulator_.vanguard().grid(), wellsForConn, useWellConn_,
const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
detail::setWellConnections(simulator_.vanguard().grid(), cartMapper, wellsForConn, useWellConn_,
wellConnectionsGraph_, numJacobiBlocks_);
std::cout << "Create block-Jacobi pattern" << std::endl;
blockJacobiAdjacency();

10
opm/simulators/linalg/findOverlapRowsAndColumns.hpp
View File

@ -23,6 +23,11 @@
#include <vector>
#include <utility>
#include <opm/grid/common/WellConnections.hpp>
#include <opm/grid/common/CartesianIndexMapper.hpp>
namespace Dune {
template<class Grid> class CartesianIndexMapper;
}
namespace Opm
{
@ -38,12 +43,11 @@ namespace detail
/// \param wells List of wells contained in grid.
/// \param useWellConn Boolean that is true when UseWellContribusion is true
/// \param wellGraph Cell IDs of well cells stored in a graph.
template<class Grid, class W>
void setWellConnections(const Grid& grid, const W& wells, bool useWellConn, std::vector<std::set<int>>& wellGraph, int numJacobiBlocks)
template<class Grid, class CartMapper, class W>
void setWellConnections(const Grid& grid, const CartMapper& cartMapper, const W& wells, bool useWellConn, std::vector<std::set<int>>& wellGraph, int numJacobiBlocks)
{
if ( grid.comm().size() > 1 || numJacobiBlocks > 1)
{
Dune::CartesianIndexMapper< Grid > cartMapper( grid );
const int numCells = cartMapper.compressedSize(); // grid.numCells()
wellGraph.resize(numCells);

1
opm/simulators/utils/ParallelEclipseState.cpp
View File

@ -273,7 +273,6 @@ void ParallelEclipseState::switchToDistributedProps()
{
if (m_comm.size() == 1) // No need for the parallel frontend
return;
m_parProps = true;
}

1
opm/simulators/utils/ParallelRestart.cpp
View File

@ -36,6 +36,7 @@
#include <opm/output/data/Wells.hpp>
#include <opm/input/eclipse/EclipseState/Util/OrderedMap.hpp>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/input/eclipse/EclipseState/Util/OrderedMap.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/input/eclipse/Schedule/SummaryState.hpp>

13
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -97,6 +97,7 @@ namespace Opm {
typedef BlackoilModelParametersEbos<TypeTag> ModelParameters;
using Grid = GetPropType<TypeTag, Properties::Grid>;
using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
@ -224,14 +225,14 @@ namespace Opm {
{
initFromRestartFile(restartValues,
this->ebosSimulator_.vanguard().transferWTestState(),
UgGridHelpers::numCells(grid()),
grid().size(0),
param_.use_multisegment_well_);
}
data::Wells wellData() const
{
auto wsrpt = this->wellState()
.report(UgGridHelpers::globalCell(this->grid()),
.report(ebosSimulator_.vanguard().globalCell().data(),
[this](const int well_index) -> bool
{
return this->wasDynamicallyShutThisTimeStep(well_index);
@ -314,6 +315,9 @@ namespace Opm {
// a vector of all the wells.
std::vector<WellInterfacePtr > well_container_{};
// map from logically cartesian cell indices to compressed ones
std::vector<int> cartesian_to_compressed_;
std::vector<bool> is_cell_perforated_{};
@ -357,6 +361,9 @@ namespace Opm {
const Grid& grid() const
{ return ebosSimulator_.vanguard().grid(); }
const EquilGrid& equilGrid() const
{ return ebosSimulator_.vanguard().equilGrid(); }
const EclipseState& eclState() const
{ return ebosSimulator_.vanguard().eclState(); }
@ -378,6 +385,8 @@ namespace Opm {
// setting the well_solutions_ based on well_state.
void updatePrimaryVariables(DeferredLogger& deferred_logger);
void setupCartesianToCompressed_();
void updateAverageFormationFactor();
void computePotentials(const std::size_t widx,

2
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@ -639,7 +639,7 @@ initFromRestartFile(const RestartValue& restartValues,
this->local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
this->initializeWellProdIndCalculators();
this->initializeWellPerfData();
initializeWellPerfData();
if (! this->wells_ecl_.empty()) {
handle_ms_well &= anyMSWellOpenLocal();

2
opm/simulators/wells/BlackoilWellModelGeneric.hpp
View File

@ -439,6 +439,8 @@ protected:
std::vector<std::reference_wrapper<ParallelWellInfo>> local_parallel_well_info_;
std::vector<WellProdIndexCalculator> prod_index_calc_;
std::vector<int> cartesian_to_compressed_;
std::vector<int> pvt_region_idx_;

27
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -51,14 +51,17 @@ namespace Opm {
// Number of cells the global grid view
global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
// Set up parallel wells
auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
// Set up cartesian mapping.
{
// this->setupCartesianToCompressed_();
auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
this->parallel_well_info_.reserve(parallel_wells.size());
for( const auto& name_bool: parallel_wells)
{
this->parallel_well_info_.emplace_back(name_bool,
ebosSimulator_.gridView().comm());
grid().comm());
}
}
this->alternative_well_rate_init_ =
@ -125,7 +128,7 @@ namespace Opm {
for ( size_t c=0; c < connectionSet.size(); c++ )
{
const auto& connection = connectionSet.get(c);
int compressed_idx = compressedIndexForInterior(connection.global_index());
int compressed_idx = cartesian_to_compressed_.at(connection.global_index());
if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
wellCells.push_back(compressed_idx);
@ -1217,7 +1220,8 @@ namespace Opm {
{
std::vector<std::vector<int>> wells;
// Create cartesian to compressed mapping
const auto& globalCell = grid().globalCell();
const EquilGrid& equilGrid = ebosSimulator_.vanguard().equilGrid();
const auto& globalCell = equilGrid.globalCell();
auto cartMap = cartesianToCompressed(grid().size(0),
globalCell.data());
@ -1627,9 +1631,6 @@ namespace Opm {
updatePrimaryVariables(deferred_logger);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
@ -1727,12 +1728,10 @@ namespace Opm {
void
BlackoilWellModel<TypeTag>::extractLegacyDepth_()
{
const auto& grid = ebosSimulator_.vanguard().grid();
const unsigned numCells = grid.size(/*codim=*/0);
depth_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
depth_[cellIdx] = UgGridHelpers::cellCenterDepth( grid, cellIdx );
const auto& eclProblem = ebosSimulator_.problem();
depth_.resize(local_num_cells_);
for (unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
}
}

14
tests/test_equil.cc
View File

@ -503,7 +503,7 @@ BOOST_AUTO_TEST_CASE(DeckAllDead)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -582,7 +582,7 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillary)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
@ -622,7 +622,7 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillaryOverlap)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -683,7 +683,7 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveOil)
// Initialize the fluid system
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -760,7 +760,7 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveGas)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -840,7 +840,7 @@ BOOST_AUTO_TEST_CASE(DeckWithRSVDAndRVVD)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@ -940,7 +940,7 @@ BOOST_AUTO_TEST_CASE(DeckWithPBVDAndPDVD)
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
eclipseState, simulator->vanguard().gridView(),
eclipseState, simulator->vanguard(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);