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

çonverted tabs to spaces

Browse Source
This commit is contained in:
goncalvesmachadoc 2019-11-01 12:56:20 +01:00
parent 9982ddf3fb
commit 8efae16087
2 changed files with 1047 additions and 1064 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
ebos/ecloutputblackoilmodule.hh
View File

@ -1277,7 +1277,6 @@ public:
case WType::WATER: return "Wat"; case WType::WATER: return "Wat";
case WType::GAS: return "Gas"; case WType::GAS: return "Gas";
case WType::MULTI: return "Multi"; case WType::MULTI: return "Multi";
default: default:
{ {
return ""; return "";
@ -1293,7 +1292,6 @@ public:
case CMode::THP: return "THP"; case CMode::THP: return "THP";
case CMode::BHP: return "BHP"; case CMode::BHP: return "BHP";
case CMode::GRUP: return "GRUP"; case CMode::GRUP: return "GRUP";
default: default:
{ {
return ""; return "";
@ -1405,7 +1403,6 @@ public:
case WType::WATER: return "Wat"; case WType::WATER: return "Wat";
case WType::GAS: return "Gas"; case WType::GAS: return "Gas";
case WType::MULTI: return "Multi"; case WType::MULTI: return "Multi";
default: default:
{ {
return ""; return "";
@ -1421,7 +1418,6 @@ public:
case CMode::THP: return "THP"; case CMode::THP: return "THP";
case CMode::BHP: return "BHP"; case CMode::BHP: return "BHP";
case CMode::GRUP: return "GRUP"; case CMode::GRUP: return "GRUP";
default: default:
{ {
return ""; return "";

499
ebos/eclwriter.hh
View File

@ -43,19 +43,19 @@
#include <opm/output/eclipse/Summary.hpp> #include <opm/output/eclipse/Summary.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp> #include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/simulators/utils/ParallelRestart.hpp>
#include <opm/grid/GridHelpers.hpp> #include <opm/grid/GridHelpers.hpp>
#include <opm/grid/utility/cartesianToCompressed.hpp> #include <opm/grid/utility/cartesianToCompressed.hpp>
#include <opm/simulators/utils/ParallelRestart.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/Valgrind.hpp> #include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/common/OpmLog/OpmLog.hpp> #include <opm/common/OpmLog/OpmLog.hpp>
#include <chrono>
#include <list> #include <list>
#include <string>
#include <utility> #include <utility>
#include <string>
#include <chrono>
#ifdef HAVE_MPI #ifdef HAVE_MPI
#include <mpi.h> #include <mpi.h>
@ -71,9 +71,11 @@ END_PROPERTIES
namespace Opm { namespace Opm {
template <class TypeTag> class EclWriter; template <class TypeTag>
class EclWriter;
template <class TypeTag> class EclOutputBlackOilModule; template <class TypeTag>
class EclOutputBlackOilModule;
/*! /*!
* \brief Detect whether two cells are direct vertical neighbours. * \brief Detect whether two cells are direct vertical neighbours.
@ -84,18 +86,16 @@ template <class TypeTag> class EclOutputBlackOilModule;
* \tparam CM The type of the cartesian index mapper. * \tparam CM The type of the cartesian index mapper.
* \param cartMapper The mapper onto cartesian indices. * \param cartMapper The mapper onto cartesian indices.
* \param cartesianToActive The mapping of cartesian indices to active indices. * \param cartesianToActive The mapping of cartesian indices to active indices.
* \param smallGlobalIndex The cartesian cell index of the cell with smaller * \param smallGlobalIndex The cartesian cell index of the cell with smaller index
* index * \param largeGlobalIndex The cartesian cell index of the cell with larger index
* \param largeGlobalIndex The cartesian cell index of the cell with larger * \return True if the cells have the same i and j indices and all cartesian cells
* index
* \return True if the cells have the same i and j indices and all cartesian
* cells
* between them are inactive. * between them are inactive.
*/ */
inline bool inline
directVerticalNeighbors(const std::array<int, 3> &cartDims, bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
const std::unordered_map<int, int> &cartesianToActive, const std::unordered_map<int,int>& cartesianToActive,
int smallGlobalIndex, int largeGlobalIndex) { int smallGlobalIndex, int largeGlobalIndex)
{
assert(smallGlobalIndex <= largeGlobalIndex); assert(smallGlobalIndex <= largeGlobalIndex);
std::array<int, 3> ijk1, ijk2; std::array<int, 3> ijk1, ijk2;
auto globalToIjk = [cartDims](int gc) { auto globalToIjk = [cartDims](int gc) {
@ -108,15 +108,16 @@ directVerticalNeighbors(const std::array<int, 3> &cartDims,
}; };
ijk1 = globalToIjk(smallGlobalIndex); ijk1 = globalToIjk(smallGlobalIndex);
ijk2 = globalToIjk(largeGlobalIndex); ijk2 = globalToIjk(largeGlobalIndex);
assert(ijk2[2] >= ijk1[2]); assert(ijk2[2]>=ijk1[2]);
if (ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) { if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
assert((largeGlobalIndex - smallGlobalIndex) % {
(cartDims[0] * cartDims[1]) == assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
0); for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
for (int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi += cartDims[0] * cartDims[1] )
gi < largeGlobalIndex; gi += cartDims[0] * cartDims[1]) { {
if (cartesianToActive.find(gi) != cartesianToActive.end()) { if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
{
return false; return false;
} }
} }
@ -140,7 +141,9 @@ directVerticalNeighbors(const std::array<int, 3> &cartDims,
* - This class requires to use the black oil model with the element * - This class requires to use the black oil model with the element
* centered finite volume discretization. * centered finite volume discretization.
*/ */
template <class TypeTag> class EclWriter { template <class TypeTag>
class EclWriter
{
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@ -158,58 +161,60 @@ template <class TypeTag> class EclWriter {
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
public: public:
static void registerParameters() { static void registerParameters()
{
EclOutputBlackOilModule<TypeTag>::registerParameters(); EclOutputBlackOilModule<TypeTag>::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncEclOutput, EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncEclOutput,
"Write the ECL-formated results in a non-blocking way " "Write the ECL-formated results in a non-blocking way (i.e., using a separate thread).");
"(i.e., using a separate thread).");
} }
// The Simulator object should preferably have been const - the // The Simulator object should preferably have been const - the
// only reason that is not the case is due to the SummaryState // only reason that is not the case is due to the SummaryState
// object owned deep down by the vanguard. // object owned deep down by the vanguard.
EclWriter(Simulator &simulator) EclWriter(Simulator& simulator)
: simulator_(simulator), collectToIORank_(simulator_.vanguard()), : simulator_(simulator)
eclOutputModule_(simulator, collectToIORank_) { , collectToIORank_(simulator_.vanguard())
, eclOutputModule_(simulator, collectToIORank_)
{
if (collectToIORank_.isIORank()) { if (collectToIORank_.isIORank()) {
globalGrid_ = simulator_.vanguard().grid(); globalGrid_ = simulator_.vanguard().grid();
globalGrid_.switchToGlobalView(); globalGrid_.switchToGlobalView();
eclIO_.reset(new Opm::EclipseIO( eclIO_.reset(new Opm::EclipseIO(simulator_.vanguard().eclState(),
simulator_.vanguard().eclState(), Opm::UgGridHelpers::createEclipseGrid(globalGrid_, simulator_.vanguard().eclState().getInputGrid()),
Opm::UgGridHelpers::createEclipseGrid(
globalGrid_, simulator_.vanguard().eclState().getInputGrid()),
simulator_.vanguard().schedule(), simulator_.vanguard().schedule(),
simulator_.vanguard().summaryConfig())); simulator_.vanguard().summaryConfig()));
} }
// create output thread if enabled and rank is I/O rank // create output thread if enabled and rank is I/O rank
// async output is enabled by default if pthread are enabled // async output is enabled by default if pthread are enabled
bool enableAsyncOutput = bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput);
EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput);
int numWorkerThreads = 0; int numWorkerThreads = 0;
if (enableAsyncOutput && collectToIORank_.isIORank()) if (enableAsyncOutput && collectToIORank_.isIORank())
numWorkerThreads = 1; numWorkerThreads = 1;
taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
} }
~EclWriter() {} ~EclWriter()
{ }
const Opm::EclipseIO &eclIO() const { const Opm::EclipseIO& eclIO() const
{
assert(eclIO_); assert(eclIO_);
return *eclIO_; return *eclIO_;
} }
void writeInit() { void writeInit()
{
if (collectToIORank_.isIORank()) { if (collectToIORank_.isIORank()) {
std::map<std::string, std::vector<int>> integerVectors; std::map<std::string, std::vector<int> > integerVectors;
if (collectToIORank_.isParallel()) if (collectToIORank_.isParallel())
integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
auto cartMap = Opm::cartesianToCompressed( auto cartMap = Opm::cartesianToCompressed(globalGrid_.size(0),
globalGrid_.size(0), Opm::UgGridHelpers::globalCell(globalGrid_)); Opm::UgGridHelpers::globalCell(globalGrid_));
eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap));
exportNncStructure_(cartMap));
} }
} }
@ -217,7 +222,8 @@ public:
* \brief collect and pass data and pass it to eclIO writer * \brief collect and pass data and pass it to eclIO writer
*/ */
void evalSummaryState(bool isSubStep) { void evalSummaryState(bool isSubStep)
{
int reportStepNum = simulator_.episodeIndex() + 1; int reportStepNum = simulator_.episodeIndex() + 1;
/* /*
The summary data is not evaluated for timestep 0, that is The summary data is not evaluated for timestep 0, that is
@ -240,32 +246,30 @@ public:
return; return;
Scalar curTime = simulator_.time() + simulator_.timeStepSize(); Scalar curTime = simulator_.time() + simulator_.timeStepSize();
Scalar totalCpuTime = simulator_.executionTimer().realTimeElapsed() + Scalar totalCpuTime =
simulator_.executionTimer().realTimeElapsed() +
simulator_.setupTimer().realTimeElapsed() + simulator_.setupTimer().realTimeElapsed() +
simulator_.vanguard().externalSetupTime(); simulator_.vanguard().externalSetupTime();
Opm::data::Wells localWellData = Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData();
simulator_.problem().wellModel().wellData();
const auto &gridView = simulator_.vanguard().gridView(); const auto& gridView = simulator_.vanguard().gridView();
int numElements = gridView.size(/*codim=*/0); int numElements = gridView.size(/*codim=*/0);
bool log = collectToIORank_.isIORank(); bool log = collectToIORank_.isIORank();
eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, /*isRestart*/ false);
/*isRestart*/ false);
ElementContext elemCtx(simulator_); ElementContext elemCtx(simulator_);
ElementIterator elemIt = gridView.template begin</*codim=*/0>(); ElementIterator elemIt = gridView.template begin</*codim=*/0>();
const ElementIterator &elemEndIt = gridView.template end</*codim=*/0>(); const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const Element &elem = *elemIt; const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_.processElement(elemCtx); eclOutputModule_.processElement(elemCtx);
} }
if (collectToIORank_.isParallel()) if (collectToIORank_.isParallel())
collectToIORank_.collect({}, eclOutputModule_.getBlockData(), collectToIORank_.collect({}, eclOutputModule_.getBlockData(), localWellData);
localWellData);
std::map<std::string, double> miscSummaryData; std::map<std::string, double> miscSummaryData;
std::map<std::string, std::vector<double>> regionData; std::map<std::string, std::vector<double>> regionData;
@ -280,64 +284,68 @@ public:
std::vector<char> buffer; std::vector<char> buffer;
if (collectToIORank_.isIORank()) { if (collectToIORank_.isIORank()) {
const auto &summary = eclIO_->summary(); const auto& summary = eclIO_->summary();
const auto &eclState = simulator_.vanguard().eclState(); const auto& eclState = simulator_.vanguard().eclState();
// Add TCPU // Add TCPU
if (totalCpuTime != 0.0) if (totalCpuTime != 0.0)
miscSummaryData["TCPU"] = totalCpuTime; miscSummaryData["TCPU"] = totalCpuTime;
const Opm::data::Wells &wellData = collectToIORank_.isParallel() const Opm::data::Wells& wellData = collectToIORank_.isParallel() ? collectToIORank_.globalWellData() : localWellData;
? collectToIORank_.globalWellData()
: localWellData;
const std::map<std::pair<std::string, int>, double> &blockData = const std::map<std::pair<std::string, int>, double>& blockData
collectToIORank_.isParallel() ? collectToIORank_.globalBlockData() = collectToIORank_.isParallel()
? collectToIORank_.globalBlockData()
: eclOutputModule_.getBlockData(); : eclOutputModule_.getBlockData();
summary.eval(summaryState(), reportStepNum, curTime, eclState, schedule(), summary.eval(summaryState(),
wellData, miscSummaryData, regionData, blockData); reportStepNum,
curTime,
eclState,
schedule(),
wellData,
miscSummaryData,
regionData,
blockData);
buffer = summaryState().serialize(); buffer = summaryState().serialize();
} }
if (collectToIORank_.isParallel()) { if (collectToIORank_.isParallel()) {
#ifdef HAVE_MPI #ifdef HAVE_MPI
unsigned long buffer_size = buffer.size(); unsigned long buffer_size = buffer.size();
MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD);
MPI_COMM_WORLD);
if (!collectToIORank_.isIORank()) if (!collectToIORank_.isIORank())
buffer.resize(buffer_size); buffer.resize( buffer_size );
MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD);
MPI_COMM_WORLD);
if (!collectToIORank_.isIORank()) { if (!collectToIORank_.isIORank()) {
Opm::SummaryState &st = summaryState(); Opm::SummaryState& st = summaryState();
st.deserialize(buffer); st.deserialize(buffer);
} }
#endif #endif
} }
} }
void writeOutput(bool isSubStep) {
void writeOutput(bool isSubStep)
{
Scalar curTime = simulator_.time() + simulator_.timeStepSize(); Scalar curTime = simulator_.time() + simulator_.timeStepSize();
Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
// output using eclWriter if enabled // output using eclWriter if enabled
Opm::data::Wells localWellData = Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData();
simulator_.problem().wellModel().wellData();
int reportStepNum = simulator_.episodeIndex() + 1; int reportStepNum = simulator_.episodeIndex() + 1;
const auto &gridView = simulator_.vanguard().gridView(); const auto& gridView = simulator_.vanguard().gridView();
int numElements = gridView.size(/*codim=*/0); int numElements = gridView.size(/*codim=*/0);
bool log = collectToIORank_.isIORank(); bool log = collectToIORank_.isIORank();
eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, /*isRestart*/ false);
/*isRestart*/ false);
ElementContext elemCtx(simulator_); ElementContext elemCtx(simulator_);
ElementIterator elemIt = gridView.template begin</*codim=*/0>(); ElementIterator elemIt = gridView.template begin</*codim=*/0>();
const ElementIterator &elemEndIt = gridView.template end</*codim=*/0>(); const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++elemIt) {
const Element &elem = *elemIt; const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem); elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_.processElement(elemCtx); eclOutputModule_.processElement(elemCtx);
@ -354,39 +362,35 @@ public:
eclOutputModule_.addRftDataToWells(localWellData, reportStepNum); eclOutputModule_.addRftDataToWells(localWellData, reportStepNum);
if (collectToIORank_.isParallel()) if (collectToIORank_.isParallel())
collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), localWellData);
localWellData);
if (collectToIORank_.isIORank()) { if (collectToIORank_.isIORank()) {
const auto &eclState = simulator_.vanguard().eclState(); const auto& eclState = simulator_.vanguard().eclState();
const auto &simConfig = eclState.getSimulationConfig(); const auto& simConfig = eclState.getSimulationConfig();
bool enableDoublePrecisionOutput = bool enableDoublePrecisionOutput = EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision);
EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision); const Opm::data::Solution& cellData = collectToIORank_.isParallel() ? collectToIORank_.globalCellData() : localCellData;
const Opm::data::Solution &cellData = const Opm::data::Wells& wellData = collectToIORank_.isParallel() ? collectToIORank_.globalWellData() : localWellData;
collectToIORank_.isParallel() ? collectToIORank_.globalCellData()
: localCellData;
const Opm::data::Wells &wellData = collectToIORank_.isParallel()
? collectToIORank_.globalWellData()
: localWellData;
Opm::RestartValue restartValue(cellData, wellData); Opm::RestartValue restartValue(cellData, wellData);
if (simConfig.useThresholdPressure()) if (simConfig.useThresholdPressure())
restartValue.addExtra("THRESHPR", Opm::UnitSystem::measure::pressure, restartValue.addExtra("THRESHPR", Opm::UnitSystem::measure::pressure, simulator_.problem().thresholdPressure().data());
simulator_.problem().thresholdPressure().data());
// Add suggested next timestep to extra data. // Add suggested next timestep to extra data.
if (!isSubStep) if (!isSubStep)
restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize)); restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
// first, create a tasklet to write the data for the current time step to // first, create a tasklet to write the data for the current time step to disk
// disk auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(summaryState(),
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>( *eclIO_,
summaryState(), *eclIO_, reportStepNum, isSubStep, curTime, reportStepNum,
restartValue, enableDoublePrecisionOutput); isSubStep,
curTime,
restartValue,
enableDoublePrecisionOutput);
// then, make sure that the previous I/O request has been completed and // then, make sure that the previous I/O request has been completed and the
// the
// number of incomplete tasklets does not increase between time steps // number of incomplete tasklets does not increase between time steps
taskletRunner_->barrier(); taskletRunner_->barrier();
@ -395,59 +399,45 @@ public:
} }
} }
void beginRestart() { void beginRestart()
bool enableHysteresis = {
simulator_.problem().materialLawManager()->enableHysteresis(); bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis();
bool enableSwatinit = simulator_.vanguard() bool enableSwatinit = simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
.eclState()
.get3DProperties()
.hasDeckDoubleGridProperty("SWATINIT");
std::vector<Opm::RestartKey> solutionKeys{ std::vector<Opm::RestartKey> solutionKeys{
{"PRESSURE", Opm::UnitSystem::measure::pressure}, {"PRESSURE", Opm::UnitSystem::measure::pressure},
{"SWAT", Opm::UnitSystem::measure::identity, {"SWAT", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
static_cast<bool>( {"SGAS", Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))},
FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))}, {"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy},
{"SGAS", Opm::UnitSystem::measure::identity, {"SSOLVENT" , Opm::UnitSystem::measure::identity, enableSolvent},
static_cast<bool>( {"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))}, {"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
{"TEMP", Opm::UnitSystem::measure::temperature, enableEnergy}, {"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()},
{"SSOLVENT", Opm::UnitSystem::measure::identity, enableSolvent},
{"RS", Opm::UnitSystem::measure::gas_oil_ratio,
FluidSystem::enableDissolvedGas()},
{"RV", Opm::UnitSystem::measure::oil_gas_ratio,
FluidSystem::enableVaporizedOil()},
{"SOMAX", Opm::UnitSystem::measure::identity,
simulator_.problem().vapparsActive()},
{"PCSWM_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, {"PCSWM_OW", Opm::UnitSystem::measure::identity, enableHysteresis},
{"KRNSW_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, {"KRNSW_OW", Opm::UnitSystem::measure::identity, enableHysteresis},
{"PCSWM_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, {"PCSWM_GO", Opm::UnitSystem::measure::identity, enableHysteresis},
{"KRNSW_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, {"KRNSW_GO", Opm::UnitSystem::measure::identity, enableHysteresis},
{"PPCW", Opm::UnitSystem::measure::pressure, enableSwatinit}}; {"PPCW", Opm::UnitSystem::measure::pressure, enableSwatinit}
};
const auto &inputThpres = const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
eclState().getSimulationConfig().getThresholdPressure(); std::vector<Opm::RestartKey> extraKeys = {{"OPMEXTRA", Opm::UnitSystem::measure::identity, false},
std::vector<Opm::RestartKey> extraKeys = {
{"OPMEXTRA", Opm::UnitSystem::measure::identity, false},
{"THRESHPR", Opm::UnitSystem::measure::pressure, inputThpres.active()}}; {"THRESHPR", Opm::UnitSystem::measure::pressure, inputThpres.active()}};
// The episodeIndex is rewined one back before beginRestart is called // The episodeIndex is rewined one back before beginRestart is called
// and can not be used here. // and can not be used here.
// We just ask the initconfig directly to be sure that we use the correct // We just ask the initconfig directly to be sure that we use the correct
// index. // index.
const auto &initconfig = simulator_.vanguard().eclState().getInitConfig(); const auto& initconfig = simulator_.vanguard().eclState().getInitConfig();
int restartStepIdx = initconfig.getRestartStep(); int restartStepIdx = initconfig.getRestartStep();
const auto &gridView = simulator_.vanguard().gridView(); const auto& gridView = simulator_.vanguard().gridView();
unsigned numElements = gridView.size(/*codim=*/0); unsigned numElements = gridView.size(/*codim=*/0);
eclOutputModule_.allocBuffers(numElements, restartStepIdx, eclOutputModule_.allocBuffers(numElements, restartStepIdx, /*isSubStep=*/false, /*log=*/false, /*isRestart*/ true);
/*isSubStep=*/false, /*log=*/false,
/*isRestart*/ true);
{ {
Opm::SummaryState &summaryState = simulator_.vanguard().summaryState(); Opm::SummaryState& summaryState = simulator_.vanguard().summaryState();
auto restartValues = auto restartValues = loadParallelRestart(eclIO_.get(), summaryState, solutionKeys, extraKeys,
loadParallelRestart(eclIO_.get(), summaryState, solutionKeys, gridView.grid().comm());
extraKeys, gridView.grid().comm());
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx); unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
@ -455,9 +445,9 @@ public:
} }
if (inputThpres.active()) { if (inputThpres.active()) {
Simulator &mutableSimulator = const_cast<Simulator &>(simulator_); Simulator& mutableSimulator = const_cast<Simulator&>(simulator_);
auto &thpres = mutableSimulator.problem().thresholdPressure(); auto& thpres = mutableSimulator.problem().thresholdPressure();
const auto &thpresValues = restartValues.getExtra("THRESHPR"); const auto& thpresValues = restartValues.getExtra("THRESHPR");
thpres.setFromRestart(thpresValues); thpres.setFromRestart(thpresValues);
} }
restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0]; restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0];
@ -467,34 +457,29 @@ public:
} }
} }
void endRestart() {} void endRestart()
{}
const EclOutputBlackOilModule<TypeTag> &eclOutputModule() const { const EclOutputBlackOilModule<TypeTag>& eclOutputModule() const
return eclOutputModule_; { return eclOutputModule_; }
}
Scalar restartTimeStepSize() const
{ return restartTimeStepSize_; }
Scalar restartTimeStepSize() const { return restartTimeStepSize_; }
private: private:
static bool enableEclOutput_() { static bool enableEclOutput_()
return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
}
Opm::data::Solution Opm::data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const
computeTrans_(const std::unordered_map<int, int> &cartesianToActive) const { {
const auto &cartMapper = simulator_.vanguard().cartesianIndexMapper(); const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
const auto &cartDims = cartMapper.cartesianDimensions(); const auto& cartDims = cartMapper.cartesianDimensions();
const int globalSize = cartDims[0] * cartDims[1] * cartDims[2]; const int globalSize = cartDims[0]*cartDims[1]*cartDims[2];
Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
std::vector<double>(globalSize), Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
Opm::data::TargetType::INIT}; Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, std::vector<double>(globalSize), Opm::data::TargetType::INIT};
Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility,
std::vector<double>(globalSize),
Opm::data::TargetType::INIT};
Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility,
std::vector<double>(globalSize),
Opm::data::TargetType::INIT};
for (size_t i = 0; i < tranx.data.size(); ++i) { for (size_t i = 0; i < tranx.data.size(); ++i) {
tranx.data[0] = 0.0; tranx.data[0] = 0.0;
@ -502,23 +487,21 @@ private:
tranz.data[0] = 0.0; tranz.data[0] = 0.0;
} }
typedef typename Grid::LeafGridView GlobalGridView; typedef typename Grid :: LeafGridView GlobalGridView;
const GlobalGridView &globalGridView = globalGrid_.leafGridView(); const GlobalGridView& globalGridView = globalGrid_.leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) #if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> ElementMapper;
ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else #else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView, typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView, Dune::MCMGElementLayout> ElementMapper;
Dune::MCMGElementLayout>
ElementMapper;
ElementMapper globalElemMapper(globalGridView); ElementMapper globalElemMapper(globalGridView);
#endif #endif
const auto &cartesianCellIdx = globalGrid_.globalCell(); const auto& cartesianCellIdx = globalGrid_.globalCell();
const EclTransmissibility<TypeTag> *globalTrans; const EclTransmissibility<TypeTag>* globalTrans;
if (!collectToIORank_.isParallel()) { if (!collectToIORank_.isParallel())
{
// in the sequential case we must use the transmissibilites defined by // in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does // the problem. (because in the sequential case, the grid manager does
// not compute "global" transmissibilities for performance reasons. in // not compute "global" transmissibilities for performance reasons. in
@ -526,19 +509,21 @@ private:
// because this object refers to the distributed grid and we need the // because this object refers to the distributed grid and we need the
// sequential version here.) // sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities(); globalTrans = &simulator_.problem().eclTransmissibilities();
} else { }
else
{
globalTrans = &(simulator_.vanguard().globalTransmissibility()); globalTrans = &(simulator_.vanguard().globalTransmissibility());
} }
auto elemIt = globalGridView.template begin</*codim=*/0>(); auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto &elemEndIt = globalGridView.template end</*codim=*/0>(); const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++ elemIt) {
const auto &elem = *elemIt; const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem); auto isIt = globalGridView.ibegin(elem);
const auto &isEndIt = globalGridView.iend(elem); const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++isIt) { for (; isIt != isEndIt; ++ isIt) {
const auto &is = *isIt; const auto& is = *isIt;
if (!is.neighbor()) if (!is.neighbor())
continue; // intersection is on the domain boundary continue; // intersection is on the domain boundary
@ -556,8 +541,7 @@ private:
if (gc2 - gc1 == 1) { if (gc2 - gc1 == 1) {
tranx.data[gc1] = globalTrans->transmissibility(c1, c2); tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
continue; // skip other if clauses as they are false, last one needs continue; // skip other if clauses as they are false, last one needs some computation
// some computation
} }
if (gc2 - gc1 == cartDims[0]) { if (gc2 - gc1 == cartDims[0]) {
@ -565,65 +549,61 @@ private:
continue; // skipt next if clause as it needs some computation continue; // skipt next if clause as it needs some computation
} }
if (gc2 - gc1 == cartDims[0] * cartDims[1] || if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
tranz.data[gc1] = globalTrans->transmissibility(c1, c2); tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
} }
} }
return {{"TRANX", tranx}, {"TRANY", trany}, {"TRANZ", tranz}}; return {{"TRANX", tranx},
{"TRANY", trany},
{"TRANZ", tranz}};
} }
Opm::NNC exportNncStructure_( Opm::NNC exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const
const std::unordered_map<int, int> &cartesianToActive) const { {
std::size_t nx = eclState().getInputGrid().getNX(); std::size_t nx = eclState().getInputGrid().getNX();
std::size_t ny = eclState().getInputGrid().getNY(); std::size_t ny = eclState().getInputGrid().getNY();
auto nncData = sortNncAndApplyEditnnc(eclState().getInputNNC().nncdata(), auto nncData = sortNncAndApplyEditnnc(eclState().getInputNNC().nncdata(),
eclState().getInputEDITNNC().data()); eclState().getInputEDITNNC().data());
const auto &unitSystem = simulator_.vanguard().deck().getActiveUnitSystem(); const auto& unitSystem = simulator_.vanguard().deck().getActiveUnitSystem();
std::vector<Opm::NNCdata> outputNnc; std::vector<Opm::NNCdata> outputNnc;
std::size_t index = 0; std::size_t index = 0;
for (const auto &entry : nncData) { for( const auto& entry : nncData ) {
// test whether NNC is not a neighboring connection // test whether NNC is not a neighboring connection
// cell2>=cell1 holds due to sortNncAndApplyEditnnc // cell2>=cell1 holds due to sortNncAndApplyEditnnc
assert(entry.cell2 >= entry.cell1); assert( entry.cell2 >= entry.cell1 );
auto cellDiff = entry.cell2 - entry.cell1; auto cellDiff = entry.cell2 - entry.cell1;
if (cellDiff != 1 && cellDiff != nx && cellDiff != nx * ny) { if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) {
auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, entry.trans);
entry.trans); // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6
// Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems if ( tt >= 1.0e-6 )
// like the threshold is 1.0e-6
if (tt >= 1.0e-6)
outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans); outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans);
} }
++index; ++index;
} }
auto nncCompare = [](const Opm::NNCdata &nnc1, const Opm::NNCdata &nnc2) { auto nncCompare = []( const Opm::NNCdata& nnc1, const Opm::NNCdata& nnc2){
return nnc1.cell1 < nnc2.cell1 || return nnc1.cell1 < nnc2.cell1 ||
(nnc1.cell1 == nnc2.cell1 && nnc1.cell2 < nnc2.cell2); ( nnc1.cell1 == nnc2.cell1 && nnc1.cell2 < nnc2.cell2);};
};
// Sort the nncData values from the deck as they need to be // Sort the nncData values from the deck as they need to be
// Checked when writing NNC transmissibilities from the simulation. // Checked when writing NNC transmissibilities from the simulation.
std::sort(nncData.begin(), nncData.end(), nncCompare); std::sort(nncData.begin(), nncData.end(), nncCompare);
typedef typename Grid::LeafGridView GlobalGridView; typedef typename Grid :: LeafGridView GlobalGridView;
const GlobalGridView &globalGridView = globalGrid_.leafGridView(); const GlobalGridView& globalGridView = globalGrid_.leafGridView();
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) #if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView> ElementMapper;
ElementMapper;
ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
#else #else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView, typedef Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView, Dune::MCMGElementLayout> ElementMapper;
Dune::MCMGElementLayout>
ElementMapper;
ElementMapper globalElemMapper(globalGridView); ElementMapper globalElemMapper(globalGridView);
#endif #endif
const EclTransmissibility<TypeTag> *globalTrans; const EclTransmissibility<TypeTag>* globalTrans;
if (!collectToIORank_.isParallel()) { if (!collectToIORank_.isParallel()) {
// in the sequential case we must use the transmissibilites defined by // in the sequential case we must use the transmissibilites defined by
// the problem. (because in the sequential case, the grid manager does // the problem. (because in the sequential case, the grid manager does
@ -632,21 +612,22 @@ private:
// because this object refers to the distributed grid and we need the // because this object refers to the distributed grid and we need the
// sequential version here.) // sequential version here.)
globalTrans = &simulator_.problem().eclTransmissibilities(); globalTrans = &simulator_.problem().eclTransmissibilities();
} else { }
else
{
globalTrans = &(simulator_.vanguard().globalTransmissibility()); globalTrans = &(simulator_.vanguard().globalTransmissibility());
} }
auto cartDims = auto cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions();
simulator_.vanguard().cartesianIndexMapper().cartesianDimensions();
auto elemIt = globalGridView.template begin</*codim=*/0>(); auto elemIt = globalGridView.template begin</*codim=*/0>();
const auto &elemEndIt = globalGridView.template end</*codim=*/0>(); const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) { for (; elemIt != elemEndIt; ++ elemIt) {
const auto &elem = *elemIt; const auto& elem = *elemIt;
auto isIt = globalGridView.ibegin(elem); auto isIt = globalGridView.ibegin(elem);
const auto &isEndIt = globalGridView.iend(elem); const auto& isEndIt = globalGridView.iend(elem);
for (; isIt != isEndIt; ++isIt) { for (; isIt != isEndIt; ++ isIt) {
const auto &is = *isIt; const auto& is = *isIt;
if (!is.neighbor()) if (!is.neighbor())
continue; // intersection is on the domain boundary continue; // intersection is on the domain boundary
@ -663,91 +644,97 @@ private:
std::size_t cc1 = globalGrid_.globalCell()[c1]; std::size_t cc1 = globalGrid_.globalCell()[c1];
std::size_t cc2 = globalGrid_.globalCell()[c2]; std::size_t cc2 = globalGrid_.globalCell()[c2];
if (cc2 < cc1) if ( cc2 < cc1 )
std::swap(cc1, cc2); std::swap(cc1, cc2);
auto cellDiff = cc2 - cc1; auto cellDiff = cc2 - cc1;
if (cellDiff != 1 && cellDiff != nx && cellDiff != nx * ny && if (cellDiff != 1 &&
!directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { cellDiff != nx &&
cellDiff != nx*ny &&
! directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) {
// We need to check whether an NNC for this face was also specified // We need to check whether an NNC for this face was also specified
// via the NNC keyword in the deck (i.e. in the first origNncSize // via the NNC keyword in the deck (i.e. in the first origNncSize entries.
// entries.
auto t = globalTrans->transmissibility(c1, c2); auto t = globalTrans->transmissibility(c1, c2);
auto candidate = auto candidate = std::lower_bound(nncData.begin(), nncData.end(), Opm::NNCdata(cc1, cc2, 0.0), nncCompare);
std::lower_bound(nncData.begin(), nncData.end(),
Opm::NNCdata(cc1, cc2, 0.0), nncCompare);
while (candidate != nncData.end() && candidate->cell1 == cc1 && while ( candidate != nncData.end() && candidate->cell1 == cc1
candidate->cell2 == cc2) { && candidate->cell2 == cc2) {
t -= candidate->trans; t -= candidate->trans;
++candidate; ++candidate;
} }
// eclipse ignores NNCs with zero transmissibility (different // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC
// threshold than for NNC // with corresponding EDITNNC above). In addition we do set small transmissibilties
// with corresponding EDITNNC above). In addition we do set small // to zero when setting up the simulator. These will be ignored here, too.
// transmissibilties auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, std::abs(t));
// to zero when setting up the simulator. These will be ignored here, if ( tt > 1e-12 )
// too.
auto tt = unitSystem.from_si(
Opm::UnitSystem::measure::transmissibility, std::abs(t));
if (tt > 1e-12)
outputNnc.push_back({cc1, cc2, t}); outputNnc.push_back({cc1, cc2, t});
} }
} }
} }
Opm::NNC ret; Opm::NNC ret;
for (const auto &nncItem : outputNnc) for(const auto& nncItem: outputNnc)
ret.addNNC(nncItem.cell1, nncItem.cell2, nncItem.trans); ret.addNNC(nncItem.cell1, nncItem.cell2, nncItem.trans);
return ret; return ret;
} }
struct EclWriteTasklet : public TaskletInterface { struct EclWriteTasklet
: public TaskletInterface
{
Opm::SummaryState summaryState_; Opm::SummaryState summaryState_;
Opm::EclipseIO &eclIO_; Opm::EclipseIO& eclIO_;
int reportStepNum_; int reportStepNum_;
bool isSubStep_; bool isSubStep_;
double secondsElapsed_; double secondsElapsed_;
Opm::RestartValue restartValue_; Opm::RestartValue restartValue_;
bool writeDoublePrecision_; bool writeDoublePrecision_;
explicit EclWriteTasklet(const Opm::SummaryState &summaryState, explicit EclWriteTasklet(const Opm::SummaryState& summaryState,
Opm::EclipseIO &eclIO, int reportStepNum, Opm::EclipseIO& eclIO,
bool isSubStep, double secondsElapsed, int reportStepNum,
bool isSubStep,
double secondsElapsed,
Opm::RestartValue restartValue, Opm::RestartValue restartValue,
bool writeDoublePrecision) bool writeDoublePrecision)
: summaryState_(summaryState), eclIO_(eclIO), : summaryState_(summaryState)
reportStepNum_(reportStepNum), isSubStep_(isSubStep), , eclIO_(eclIO)
secondsElapsed_(secondsElapsed), restartValue_(restartValue), , reportStepNum_(reportStepNum)
writeDoublePrecision_(writeDoublePrecision) {} , isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(restartValue)
, writeDoublePrecision_(writeDoublePrecision)
{ }
// callback to eclIO serial writeTimeStep method // callback to eclIO serial writeTimeStep method
void run() { void run()
eclIO_.writeTimeStep(summaryState_, reportStepNum_, isSubStep_, {
secondsElapsed_, restartValue_, eclIO_.writeTimeStep(summaryState_,
reportStepNum_,
isSubStep_,
secondsElapsed_,
restartValue_,
writeDoublePrecision_); writeDoublePrecision_);
} }
}; };
const Opm::EclipseState &eclState() const { const Opm::EclipseState& eclState() const
return simulator_.vanguard().eclState(); { return simulator_.vanguard().eclState(); }
}
Opm::SummaryState &summaryState() { Opm::SummaryState& summaryState()
return simulator_.vanguard().summaryState(); { return simulator_.vanguard().summaryState(); }
}
const Opm::Schedule &schedule() const { const Opm::Schedule& schedule() const
return simulator_.vanguard().schedule(); { return simulator_.vanguard().schedule(); }
}
Simulator &simulator_; Simulator& simulator_;
CollectDataToIORankType collectToIORank_; CollectDataToIORankType collectToIORank_;
EclOutputBlackOilModule<TypeTag> eclOutputModule_; EclOutputBlackOilModule<TypeTag> eclOutputModule_;
std::unique_ptr<Opm::EclipseIO> eclIO_; std::unique_ptr<Opm::EclipseIO> eclIO_;
Grid globalGrid_; Grid globalGrid_;
std::unique_ptr<TaskletRunner> taskletRunner_; std::unique_ptr<TaskletRunner> taskletRunner_;
Scalar restartTimeStepSize_; Scalar restartTimeStepSize_;
}; };
} // namespace Opm } // namespace Opm