Add ecl output from opm-output

The new eclwriter output and restart using the eclIO from opm-output

All tests in opm-simulator pass and
MPI restart works.

Standard initialization is done prior to restart in order
to compute correct initial thpressure values. This is not
necessary if thpressures are written to the restart file

TODO: Some trivial fields are written out in order mimic legacy code,
this should be cleaned up

TODO: Output of wells, FIP, NNC is still done in opm-simulators.
This should be moved later.
This commit is contained in:
Tor Harald Sandve 2017-12-07 12:38:48 +01:00
parent 4957a301f1
commit 0132c3326e
4 changed files with 1036 additions and 527 deletions

View File

@ -23,6 +23,9 @@
#ifndef EWOMS_PARALLELSERIALOUTPUT_HH
#define EWOMS_PARALLELSERIALOUTPUT_HH
#include <opm/output/data/Cells.hpp>
#include <opm/output/data/Solution.hpp>
//#if HAVE_OPM_GRID
#include <dune/grid/common/p2pcommunicator.hh>
#include <dune/grid/utility/persistentcontainer.hh>
@ -90,9 +93,6 @@ namespace Ewoms
IndexMapType& localIndexMap_;
IndexMapStorageType& indexMaps_;
std::map< const int, const int > globalPosition_;
#ifndef NDEBUG
std::set< int > checkPosition_;
#endif
public:
DistributeIndexMapping( const std::vector<int>& globalIndex,
@ -111,7 +111,7 @@ namespace Ewoms
globalPosition_.insert( std::make_pair( globalIndex[ index ], index ) );
}
// on I/O rank we need to create a mapping from local to global
// we need to create a mapping from local to global
if( ! indexMaps_.empty() )
{
// for the ioRank create a localIndex to index in global state map
@ -122,10 +122,6 @@ namespace Ewoms
{
const int id = distributedGlobalIndex_[ localIndexMap_[ i ] ];
indexMap[ i ] = globalPosition_[ id ] ;
#ifndef NDEBUG
assert( checkPosition_.find( id ) == checkPosition_.end() );
checkPosition_.insert( id );
#endif
}
}
}
@ -164,10 +160,6 @@ namespace Ewoms
buffer.read( globalId );
assert( globalPosition_.find( globalId ) != globalPosition_.end() );
indexMap[ index ] = globalPosition_[ globalId ];
#ifndef NDEBUG
assert( checkPosition_.find( globalId ) == checkPosition_.end() );
checkPosition_.insert( globalId );
#endif
}
}
};
@ -178,12 +170,10 @@ namespace Ewoms
typename GridManager::Grid, typename GridManager::EquilGrid > :: value ;
CollectDataToIORank( const GridManager& gridManager )
: toIORankComm_( ),
isIORank_( gridManager.grid().comm().rank() == ioRank ),
isParallel_( gridManager.grid().comm().size() > 1 )
: toIORankComm_( )
{
// index maps only have to be build when reordering is needed
if( ! needsReordering && ! isParallel_ )
if( ! needsReordering && ! isParallel() )
{
return ;
}
@ -192,36 +182,35 @@ namespace Ewoms
{
std::set< int > send, recv;
typedef typename GridManager::EquilGrid::LeafGridView EquilGridView;
const EquilGridView equilGridView = gridManager.equilGrid().leafGridView() ;
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView, Dune::MCMGElementLayout> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView);
#endif
// We need a mapping from local to global grid, here we
// use equilGrid which represents a view on the global grid
const size_t globalSize = gridManager.equilGrid().leafGridView().size( 0 );
// reserve memory
globalCartesianIndex_.resize(globalSize, -1);
// loop over all elements (global grid) and store Cartesian index
auto elemIt = gridManager.equilGrid().leafGridView().template begin<0>();
const auto& elemEndIt = gridManager.equilGrid().leafGridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = equilElemMapper.index(*elemIt );
int cartElemIdx = gridManager.equilCartesianIndexMapper().cartesianIndex(elemIdx);
globalCartesianIndex_[elemIdx] = cartElemIdx;
}
// the I/O rank receives from all other ranks
if( isIORank() )
{
typedef typename GridManager::EquilGrid::LeafGridView EquilGridView;
const EquilGridView equilGridView = gridManager.equilGrid().leafGridView() ;
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
#else
typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView, Dune::MCMGElementLayout> EquilElementMapper;
EquilElementMapper equilElemMapper(equilGridView);
#endif
// the I/O rank needs a picture of the global grid, here we
// use equilGrid which represents a view on the global grid
const size_t globalSize = gridManager.equilGrid().leafGridView().size( 0 );
// reserve memory
globalCartesianIndex_.resize(globalSize, -1);
// loop over all elements (global grid) and store Cartesian index
auto elemIt = gridManager.equilGrid().leafGridView().template begin<0>();
const auto& elemEndIt = gridManager.equilGrid().leafGridView().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
int elemIdx = equilElemMapper.index(*elemIt );
int cartElemIdx = gridManager.equilCartesianIndexMapper().cartesianIndex(elemIdx);
globalCartesianIndex_[elemIdx] = cartElemIdx;
}
for(int i=0; i<comm.size(); ++i)
{
if( i != ioRank )
@ -254,15 +243,16 @@ namespace Ewoms
ElementMapper elemMapper(localGridView);
#endif
for( auto it = localGridView.template begin< 0, Dune::Interior_Partition >(),
end = localGridView.template end< 0, Dune::Interior_Partition >(); it != end; ++it )
// A mapping for the whole grid (including the ghosts) is needed for restarts
for( auto it = localGridView.template begin< 0 >(),
end = localGridView.template end< 0 >(); it != end; ++it )
{
const auto element = *it ;
int elemIdx = elemMapper.index( element );
distributedCartesianIndex[elemIdx] = gridManager.cartesianIndex( elemIdx );
// only store interior element for collection
assert( element.partitionType() == Dune :: InteriorEntity );
//assert( element.partitionType() == Dune :: InteriorEntity );
localIndexMap_.push_back( elemIdx );
}
@ -270,12 +260,9 @@ namespace Ewoms
// insert send and recv linkage to communicator
toIORankComm_.insertRequest( send, recv );
if( isIORank() )
{
// need an index map for each rank
indexMaps_.clear();
indexMaps_.resize( comm.size() );
}
// need an index map for each rank
indexMaps_.clear();
indexMaps_.resize( comm.size() );
// distribute global id's to io rank for later association of dof's
DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_ );
@ -283,33 +270,40 @@ namespace Ewoms
}
}
template <class BufferList>
class PackUnPackOutputBuffers : public P2PCommunicatorType::DataHandleInterface
class PackUnPack : public P2PCommunicatorType::DataHandleInterface
{
BufferList& bufferList_;
const Opm::data::Solution& localCellData_;
Opm::data::Solution& globalCellData_;
const IndexMapType& localIndexMap_;
const IndexMapStorageType& indexMaps_;
public:
PackUnPackOutputBuffers( BufferList& bufferList,
const IndexMapType& localIndexMap,
const IndexMapStorageType& indexMaps,
const size_t globalSize,
const bool isIORank )
: bufferList_( bufferList ),
PackUnPack( const Opm::data::Solution& localCellData,
Opm::data::Solution& globalCellData,
const IndexMapType& localIndexMap,
const IndexMapStorageType& indexMaps,
const size_t globalSize,
const bool isIORank )
: localCellData_( localCellData ),
globalCellData_( globalCellData ),
localIndexMap_( localIndexMap ),
indexMaps_( indexMaps )
{
if( isIORank )
{
// add missing data to global cell data
for (const auto& pair : localCellData_) {
const std::string& key = pair.first;
std::size_t container_size = globalSize;
auto OPM_OPTIM_UNUSED ret = globalCellData_.insert(key, pair.second.dim,
std::vector<double>(container_size),
pair.second.target);
assert(ret.second);
}
MessageBufferType buffer;
pack( 0, buffer );
// resize all buffers
for (auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
{
it->second->resize( globalSize );
}
// the last index map is the local one
doUnpack( indexMaps.back(), buffer );
@ -324,22 +318,26 @@ namespace Ewoms
OPM_THROW(std::logic_error,"link in method pack is not 0 as expected");
}
size_t buffers = bufferList_.size();
buffer.write( buffers );
for (auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
{
write( buffer, localIndexMap_, *(it->second) );
// write all cell data registered in local state
for (const auto& pair : localCellData_) {
const auto& data = pair.second.data;
// write all data from local data to buffer
write( buffer, localIndexMap_, data);
}
}
void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
{
size_t buffers = 0;
buffer.read( buffers );
assert( buffers == bufferList_.size() );
for( auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
{
read( buffer, indexMap, *(it->second) );
// we loop over the data as
// its order governs the order the data got received.
for (auto& pair : localCellData_) {
const std::string& key = pair.first;
auto& data = globalCellData_.data(key);
//write all data from local cell data to buffer
read( buffer, indexMap, data);
}
}
@ -351,71 +349,63 @@ namespace Ewoms
protected:
template <class Vector>
void write( MessageBufferType& buffer, const IndexMapType& localIndexMap, const Vector& data ) const
void write( MessageBufferType& buffer,
const IndexMapType& localIndexMap,
const Vector& vector,
const unsigned int offset = 0,
const unsigned int stride = 1 ) const
{
const size_t size = localIndexMap.size();
assert( size <= data.size() );
unsigned int size = localIndexMap.size();
buffer.write( size );
for( size_t i=0; i<size; ++i )
assert( vector.size() >= stride * size );
for( unsigned int i=0; i<size; ++i )
{
buffer.write( data[ localIndexMap[ i ] ] );
const unsigned int index = localIndexMap[ i ] * stride + offset;
assert( index < vector.size() );
buffer.write( vector[ index ] );
}
}
template <class Vector>
void read( MessageBufferType& buffer, const IndexMapType& indexMap, Vector& data ) const
void read( MessageBufferType& buffer,
const IndexMapType& indexMap,
Vector& vector,
const unsigned int offset = 0, const unsigned int stride = 1 ) const
{
size_t size = indexMap.size();
assert( size <= data.size() );
unsigned int size = 0;
buffer.read( size );
assert( size == indexMap.size() );
for( size_t i=0; i<size; ++i )
for( unsigned int i=0; i<size; ++i )
{
buffer.read( data[ indexMap[ i ] ] );
const unsigned int index = indexMap[ i ] * stride + offset;
assert( index < vector.size() );
buffer.read( vector[ index ] );
}
}
void writeString( MessageBufferType& buffer, const std::string& s) const
{
const int size = s.size();
buffer.write( size );
for( int i=0; i<size; ++i )
{
buffer.write( s[ i ] );
}
}
void readString( MessageBufferType& buffer, std::string& s) const
{
int size = -1;
buffer.read( size );
s.resize( size );
for( int i=0; i<size; ++i )
{
buffer.read( s[ i ] );
}
}
};
// gather solution to rank 0 for EclipseWriter
template <class BufferList>
void collect( BufferList& bufferList ) const
void collect( const Opm::data::Solution& localCellData )
{
globalCellData_ = {};
// index maps only have to be build when reordering is needed
if( ! needsReordering && ! isParallel_ )
if( ! needsReordering && ! isParallel() )
{
return ;
}
// this also packs and unpacks the local buffers one ioRank
PackUnPackOutputBuffers< BufferList >
packUnpack( bufferList,
PackUnPack
packUnpack( localCellData,
globalCellData_,
localIndexMap_,
indexMaps_,
numCells(),
isIORank() );
if ( ! isParallel_ )
if ( ! isParallel() )
{
// no need to collect anything.
return;
@ -430,9 +420,14 @@ namespace Ewoms
#endif
}
const Opm::data::Solution& globalCellData() const
{
return globalCellData_;
}
bool isIORank() const
{
return isIORank_;
return toIORankComm_.rank() == ioRank;
}
bool isParallel() const
@ -440,6 +435,23 @@ namespace Ewoms
return toIORankComm_.size() > 1;
}
int localIdxToGlobalIdx(const unsigned localIdx) {
if ( ! isParallel() )
{
return localIdx;
}
// the last indexMap is the local one
IndexMapType& indexMap = indexMaps_.back();
if( indexMap.empty() )
OPM_THROW(std::logic_error,"index map is not created on this rank");
if (localIdx > indexMap.size())
OPM_THROW(std::logic_error,"local index is outside map range");
return indexMap[localIdx];
}
size_t numCells () const { return globalCartesianIndex_.size(); }
protected:
@ -447,10 +459,7 @@ namespace Ewoms
IndexMapType globalCartesianIndex_;
IndexMapType localIndexMap_;
IndexMapStorageType indexMaps_;
// true if we are on I/O rank
bool isIORank_;
/// \brief True if there is more than one MPI process
bool isParallel_;
Opm::data::Solution globalCellData_;
};
} // end namespace Opm

File diff suppressed because it is too large Load Diff

View File

@ -54,13 +54,11 @@
#include "eclwellmanager.hh"
#include "eclequilinitializer.hh"
#include "eclwriter.hh"
#include "eclsummarywriter.hh"
#include "ecloutputblackoilmodule.hh"
#include "ecltransmissibility.hh"
#include "eclthresholdpressure.hh"
#include "ecldummygradientcalculator.hh"
#include "eclfluxmodule.hh"
#include "ecldeckunits.hh"
#include <ewoms/common/pffgridvector.hh>
#include <ewoms/models/blackoil/blackoilmodel.hh>
@ -80,6 +78,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
@ -87,6 +86,8 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <boost/date_time.hpp>
#include <vector>
@ -205,9 +206,6 @@ SET_BOOL_PROP(EclBaseProblem, EnableVtkOutput, false);
// ... but enable the ECL output by default
SET_BOOL_PROP(EclBaseProblem, EnableEclOutput, true);
// also enable the summary output.
SET_BOOL_PROP(EclBaseProblem, EnableEclSummaryOutput, true);
// the cache for intensive quantities can be used for ECL problems and also yields a
// decent speedup...
SET_BOOL_PROP(EclBaseProblem, EnableIntensiveQuantityCache, true);
@ -292,11 +290,13 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
/*enableTemperature=*/true> InitialFluidState;
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Ewoms::EclSummaryWriter<TypeTag> EclSummaryWriter;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef EclWriter<TypeTag> EclWriterType;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
struct RockParams {
Scalar referencePressure;
Scalar compressibility;
@ -310,8 +310,6 @@ public:
{
ParentType::registerParameters();
Ewoms::EclOutputBlackOilModule<TypeTag>::registerParameters();
EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWriteAllSolutions,
"Write all solutions to disk instead of only the ones for the "
"report steps");
@ -330,20 +328,17 @@ public:
, transmissibilities_(simulator.gridManager())
, thresholdPressures_(simulator)
, wellManager_(simulator)
, deckUnits_(simulator)
, eclWriter_( EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput)
? new EclWriterType(simulator) : nullptr )
, summaryWriter_(simulator)
, pffDofData_(simulator.gridView(), this->elementMapper())
{
// add the output module for the Ecl binary output
simulator.model().addOutputModule(new Ewoms::EclOutputBlackOilModule<TypeTag>(simulator));
// Tell the solvent module to initialize its internal data structures
// Tell the extra modules to initialize its internal data structures
const auto& gridManager = simulator.gridManager();
SolventModule::initFromDeck(gridManager.deck(), gridManager.eclState());
PolymerModule::initFromDeck(gridManager.deck(), gridManager.eclState());
// Hack to compute the initial thpressure values for restarts
restartApplied = false;
}
/*!
@ -502,14 +497,6 @@ public:
{
if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
wellManager_.beginTimeStep();
// this is a little hack to write the initial condition, which we need to do
// before the first time step has finished.
static bool initialWritten = false;
if (this->simulator().episodeIndex() == 0 && !initialWritten) {
summaryWriter_.write(wellManager_, /*isInitial=*/true);
initialWritten = true;
}
}
}
@ -548,9 +535,6 @@ public:
if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
wellManager_.endTimeStep();
// write the summary information after each time step
summaryWriter_.write(wellManager_);
}
// we no longer need the initial soluiton
@ -613,28 +597,41 @@ public:
*/
void writeOutput(bool verbose = true)
{
// calculate the time _after_ the time was updated
Scalar t = this->simulator().time() + this->simulator().timeStepSize();
// prepare the ECL and the VTK writers
if ( eclWriter_ )
eclWriter_->beginWrite(t);
Opm::data::Wells dw;
if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
using rt = Opm::data::Rates::opt;
for (unsigned wellIdx = 0; wellIdx < wellManager_.numWells(); ++wellIdx) {
const auto& well = wellManager_.well(wellIdx);
auto& wellOut = dw[ well->name() ];
wellOut.bhp = well->bottomHolePressure();
wellOut.thp = well->tubingHeadPressure();
wellOut.temperature = 0;
wellOut.rates.set( rt::wat, well->surfaceRate(waterPhaseIdx) );
wellOut.rates.set( rt::oil, well->surfaceRate(oilPhaseIdx) );
wellOut.rates.set( rt::gas, well->surfaceRate(gasPhaseIdx) );
}
}
Scalar totalSolverTime = 0.0;
Scalar nextstep = this->simulator().timeStepSize();
Opm::data::Solution fip;
writeOutput(dw, t, false, totalSolverTime, nextstep, fip, verbose);
}
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip, bool verbose = true)
{
// use the generic code to prepare the output fields and to
// write the desired VTK files.
ParentType::writeOutput(verbose);
// output using eclWriter if enabled
if ( eclWriter_ ) {
this->model().appendOutputFields(*eclWriter_);
eclWriter_->endWrite();
eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep, fip);
}
}
/*!
* \brief Returns the object which converts between SI and deck units.
*/
const EclDeckUnits<TypeTag>& deckUnits() const
{ return deckUnits_; }
}
/*!
* \copydoc FvBaseMultiPhaseProblem::intrinsicPermeability
@ -960,6 +957,7 @@ public:
*/
void initialSolutionApplied()
{
if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
// initialize the wells. Note that this needs to be done after initializing the
// intrinsic permeabilities and the after applying the initial solution because
@ -967,11 +965,24 @@ public:
wellManager_.init(this->simulator().gridManager().eclState(), this->simulator().gridManager().schedule());
}
// the initialSolutionApplied is called recursively by readEclRestartSolution_()
// in order to setup the inital threshold pressures correctly
if (restartApplied)
return;
// let the object for threshold pressures initialize itself. this is done only at
// this point, because determining the threshold pressures may require to access
// the initial solution.
thresholdPressures_.finishInit();
const auto& eclState = this->simulator().gridManager().eclState();
const auto& initconfig = eclState.getInitConfig();
if(initconfig.restartRequested()) {
restartApplied = true;
this->simulator().setEpisodeIndex(initconfig.getRestartStep());
readEclRestartSolution_();
}
// release the memory of the EQUIL grid since it's no longer needed after this point
this->simulator().gridManager().releaseEquilGrid();
}
@ -1013,6 +1024,13 @@ public:
return initialFluidStates_[globalDofIdx];
}
void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO) {
eclWriter_->setEclIO(std::move(eclIO));
}
const Opm::EclipseIO& eclIO() const
{return eclWriter_->eclIO();}
private:
Scalar cellCenterDepth( const Element& element ) const
{
@ -1205,14 +1223,15 @@ private:
void readInitialCondition_()
{
const auto& gridManager = this->simulator().gridManager();
const auto& deck = gridManager.deck();
const auto& deck = gridManager.deck();
if (!deck.hasKeyword("EQUIL"))
readExplicitInitialCondition_();
else
readEquilInitialCondition_();
readBlackoilExtentionsInitialConditions_();
}
@ -1235,6 +1254,35 @@ private:
}
}
void readEclRestartSolution_()
{
// since the EquilInitializer provides fluid states that are consistent with the
// black-oil model, we can use naive instead of mass conservative determination
// of the primary variables.
useMassConservativeInitialCondition_ = false;
eclWriter_->restartBegin();
size_t numElems = this->model().numGridDof();
initialFluidStates_.resize(numElems);
if (enableSolvent)
solventSaturation_.resize(numElems,0.0);
if (enablePolymer)
polymerConcentration_.resize(numElems,0.0);
for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
auto& elemFluidState = initialFluidStates_[elemIdx];
eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
if (enableSolvent)
solventSaturation_[elemIdx] = eclWriter_->eclOutputModule().getSolventSaturation(elemIdx);
if (enablePolymer)
polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
}
this->model().applyInitialSolution();
}
void readExplicitInitialCondition_()
{
const auto& gridManager = this->simulator().gridManager();
@ -1572,12 +1620,12 @@ private:
EclWellManager<TypeTag> wellManager_;
EclDeckUnits<TypeTag> deckUnits_;
std::unique_ptr< EclWriterType > eclWriter_;
EclSummaryWriter summaryWriter_;
PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
bool restartApplied;
};
} // namespace Ewoms

View File

@ -30,11 +30,12 @@
#include <opm/material/densead/Evaluation.hpp>
#include "ertwrappers.hh"
#include "collecttoiorank.hh"
#include "ecloutputblackoilmodule.hh"
#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
#include <ewoms/io/baseoutputwriter.hh>
#include <opm/output/eclipse/EclipseIO.hpp>
#include <opm/common/Valgrind.hpp>
#include <opm/common/ErrorMacros.hpp>
@ -58,302 +59,169 @@ NEW_PROP_TAG(EnableEclOutput);
template <class TypeTag>
class EclWriter;
template <class TypeTag, class GridManagerType>
class EclWriterHelper
{
friend class EclWriter<TypeTag>;
static void writeHeaders_(EclWriter<TypeTag>& writer)
{
typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
OPM_THROW(std::logic_error,
"Ecl binary output only works for the element centered "
"finite volume discretization.");
#if ! HAVE_ERT
OPM_THROW(std::logic_error,
"Ecl binary output requires the ERT libraries");
#else
// set the index of the first time step written to 0...
writer.reportStepIdx_ = 0;
char* egridRawFileName = ecl_util_alloc_filename(/*outputDir=*/"./",
writer.caseName().c_str(),
ECL_EGRID_FILE,
/*formatted=*/false, // -> write binary output
writer.reportStepIdx_);
std::string egridFileName(egridRawFileName);
std::free(egridRawFileName);
ErtGrid ertGrid(writer.simulator_.gridManager().eclState().getInputGrid(),
writer.simulator_.gridManager().grid(),
writer.simulator_.gridManager().cartesianIndexMapper(),
writer.simulator_.problem().deckUnits());
ertGrid.write(egridFileName, writer.reportStepIdx_);
#endif
}
};
/*!
* \ingroup EclBlackOilSimulator
*
* \brief Implements writing Ecl binary output files.
* \brief Collects necessary output values and pass it to opm-output.
*
* Caveats:
* - For this class to do do anything meaningful, you need to have the
* ERT libraries with development headers installed and the ERT
* build system test must pass sucessfully.
* - For this class to do do anything meaningful, you will have to
* have the OPM module opm-output.
* - The only DUNE grid which is currently supported is Dune::CpGrid
* from the OPM module "opm-core". Using another grid won't
* from the OPM module "opm-grid". Using another grid won't
* fail at compile time but you will provoke a fatal exception as
* soon as you try to write an ECL output file.
* - This class requires to use the black oil model with the element
* centered finite volume discretization.
* - MPI-parallel computations are not (yet?) supported.
*/
template <class TypeTag>
class EclWriter : public BaseOutputWriter
class EclWriter
{
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, GridManager) GridManager;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef CollectDataToIORank< GridManager > CollectDataToIORankType;
typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
typedef BaseOutputWriter::VectorBuffer VectorBuffer;
typedef BaseOutputWriter::TensorBuffer TensorBuffer;
typedef std::vector<Scalar> ScalarBuffer;
friend class EclWriterHelper<TypeTag, GridManager>;
public:
EclWriter(const Simulator& simulator)
: simulator_(simulator)
, gridView_(simulator_.gridView())
#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
, elementMapper_(gridView_, Dune::mcmgElementLayout())
, vertexMapper_(gridView_, Dune::mcmgVertexLayout())
#else
, elementMapper_(gridView_)
, vertexMapper_(gridView_)
#endif
, eclOutputModule_(simulator)
, collectToIORank_( simulator_.gridManager() )
{
reportStepIdx_ = 0;
Grid globalGrid = simulator_.gridManager().grid();
globalGrid.switchToGlobalView();
eclIO_.reset(new Opm::EclipseIO(simulator_.gridManager().eclState(),
Opm::UgGridHelpers::createEclipseGrid( globalGrid , simulator_.gridManager().eclState().getInputGrid() ),
simulator_.gridManager().schedule(),
simulator_.gridManager().summaryConfig()));
}
~EclWriter()
{ }
/*!
* \brief Returns the name of the simulation.
*
* This is the prefix of the files written to disk.
*/
std::string caseName() const
{ return boost::to_upper_copy(simulator_.problem().name()); }
void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO) {
eclIO_ = std::move(eclIO);
}
const Opm::EclipseIO& eclIO() const
{return *eclIO_;}
/*!
* \brief Updates the internal data structures after mesh
* refinement.
*
* If the grid changes between two calls of beginWrite(), this
* method _must_ be called before the second beginWrite()!
* \brief collect and pass data and pass it to eclIO writer
*/
void gridChanged()
void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip)
{
elementMapper_.update();
vertexMapper_.update();
}
/*!
* \brief Called whenever a new time step must be written.
*/
void beginWrite(double t OPM_UNUSED)
{
if (enableEclOutput_() && reportStepIdx_ == 0 && collectToIORank_.isIORank() )
EclWriterHelper<TypeTag, GridManager>::writeHeaders_(*this);
}
#if !HAVE_OPM_OUTPUT
OPM_THROW(std::runtime_error,
"Opm-output must be available to write ECL output!");
#else
/*
* \brief Add a vertex-centered scalar field to the output.
*
* For the EclWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachScalarVertexData(ScalarBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
{
OPM_THROW(std::logic_error,
"The EclWriter can only write element based quantities!");
}
int episodeIdx = simulator_.episodeIndex() + 1;
const auto& gridView = simulator_.gridManager().gridView();
int numElements = gridView.size(/*codim=*/0);
bool log = collectToIORank_.isIORank();
eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), log);
/*
* \brief Add a vertex-centered vector field to the output.
*
* For the EclWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachVectorVertexData(VectorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
{
OPM_THROW(std::logic_error,
"The EclWriter can only write element based quantities!");
}
/*
* \brief Add a vertex-centered tensor field to the output.
*/
void attachTensorVertexData(TensorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
{
OPM_THROW(std::logic_error,
"The EclWriter can only write element based quantities!");
}
/*!
* \brief Add a scalar quantity to the output.
*
* The buffer must exist at least until the call to endWrite()
* finishes. Modifying the buffer between the call to this method
* and endWrite() results in _undefined behavior_.
*/
void attachScalarElementData(ScalarBuffer& buf, std::string name)
{
attachedBuffers_.push_back(std::pair<std::string, ScalarBuffer*>(name, &buf));
}
/*
* \brief Add a element-centered vector field to the output.
*
* For the EclWriter, this method is a no-op which throws a
* std::logic_error exception
*/
void attachVectorElementData(VectorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
{
OPM_THROW(std::logic_error,
"Currently, the EclWriter can only write scalar quantities!");
}
/*
* \brief Add a element-centered tensor field to the output.
*/
void attachTensorElementData(TensorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
{
OPM_THROW(std::logic_error,
"Currently, the EclWriter can only write scalar quantities!");
}
/*!
* \brief Finalizes the current writer.
*
* This means that everything will be written to disk, except if
* the onlyDiscard argument is true. In this case, no output is
* written but the 'janitorial' jobs at the end of a time step are
* still done.
*/
void endWrite(bool onlyDiscard = false)
{
if (onlyDiscard || !enableEclOutput_() || !simulator_.episodeWillBeOver()) {
// detach all buffers
attachedBuffers_.clear();
return;
ElementContext elemCtx(simulator_);
ElementIterator elemIt = gridView.template begin</*codim=*/0>();
const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
eclOutputModule_.processElement(elemCtx);
}
eclOutputModule_.outputErrorLog();
#if !HAVE_ERT
OPM_THROW(std::runtime_error,
"The ERT libraries must be available to write ECL output!");
#else
// collect all data to I/O rank and store in attachedBuffers_
// this also reorders the data such that it fits the underlying eclGrid
collectToIORank_.collect( attachedBuffers_ );
// collect all data to I/O rank and assign to sol
Opm::data::Solution localCellData = fip;
eclOutputModule_.assignToSolution(localCellData);
collectToIORank_.collect(localCellData);
// write output on I/O rank
if (collectToIORank_.isIORank()) {
ErtRestartFile restartFile(simulator_, reportStepIdx_);
restartFile.writeHeader(simulator_, reportStepIdx_);
ErtSolution solution(restartFile);
auto bufIt = attachedBuffers_.begin();
const auto& bufEndIt = attachedBuffers_.end();
for (; bufIt != bufEndIt; ++ bufIt) {
const std::string& name = bufIt->first;
const ScalarBuffer& buffer = *bufIt->second;
std::map<std::string, std::vector<double>> extraRestartData;
std::map<std::string, double> miscSummaryData;
std::shared_ptr<const ErtKeyword<float>>
bufKeyword(new ErtKeyword<float>(name, buffer));
solution.add(bufKeyword);
// Add suggested next timestep to extra data.
extraRestartData["OPMEXTRA"] = std::vector<double>(1, nextstep);
// Add TCPU if simulatorReport is not defaulted.
if (totalSolverTime != 0.0) {
miscSummaryData["TCPU"] = totalSolverTime;
}
const Opm::data::Solution& cellData = collectToIORank_.isParallel() ? collectToIORank_.globalCellData() : localCellData;
eclIO_->writeTimeStep(episodeIdx,
substep,
t,
cellData,
dw,
miscSummaryData,
extraRestartData,
false);
}
// detach all buffers
attachedBuffers_.clear();
// next time we take the next report step
++ reportStepIdx_;
#endif
}
/*!
* \brief Write the multi-writer's state to a restart file.
*/
template <class Restarter>
void serialize(Restarter& res)
{
res.serializeSectionBegin("EclWriter");
res.serializeStream() << reportStepIdx_ << "\n";
res.serializeSectionEnd();
void restartBegin() {
std::map<std::string, Opm::RestartKey> solution_keys {{"PRESSURE" , Opm::RestartKey(Opm::UnitSystem::measure::pressure)},
{"SWAT" , Opm::RestartKey(Opm::UnitSystem::measure::identity)},
{"SGAS" , Opm::RestartKey(Opm::UnitSystem::measure::identity)},
{"TEMP" , Opm::RestartKey(Opm::UnitSystem::measure::temperature)},
{"RS" , Opm::RestartKey(Opm::UnitSystem::measure::gas_oil_ratio)},
{"RV" , Opm::RestartKey(Opm::UnitSystem::measure::oil_gas_ratio)},
{"SOMAX", {Opm::UnitSystem::measure::identity, false}},
{"PCSWM_OW", {Opm::UnitSystem::measure::identity, false}},
{"KRNSW_OW", {Opm::UnitSystem::measure::identity, false}},
{"PCSWM_GO", {Opm::UnitSystem::measure::identity, false}},
{"KRNSW_GO", {Opm::UnitSystem::measure::identity, false}}};
std::map<std::string, bool> extra_keys {
{"OPMEXTRA" , false}
};
unsigned episodeIdx = simulator_.episodeIndex();
const auto& gridView = simulator_.gridManager().gridView();
unsigned numElements = gridView.size(/*codim=*/0);
eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), false);
auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
eclOutputModule_.setRestart(restart_values.solution, elemIdx, globalIdx);
}
}
/*!
* \brief Read the multi-writer's state from a restart file.
*/
template <class Restarter>
void deserialize(Restarter& res)
{
res.deserializeSectionBegin("EclWriter");
res.deserializeStream() >> reportStepIdx_;
res.deserializeSectionEnd();
const EclOutputBlackOilModule<TypeTag>& eclOutputModule() const {
return eclOutputModule_;
}
private:
static bool enableEclOutput_()
{ return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
// make sure the field is well defined if running under valgrind
// and make sure that all values can be displayed by paraview
void sanitizeBuffer_(std::vector<float>& b)
{
static bool warningPrinted = false;
for (size_t i = 0; i < b.size(); ++i) {
Opm::Valgrind::CheckDefined(b[i]);
if (!warningPrinted && !std::isfinite(b[i])) {
std::cerr << "WARNING: data field written to disk contains non-finite entries!\n";
warningPrinted = true;
}
// set values which are too small to 0 to avoid possible
// problems
if (std::abs(b[i]) < std::numeric_limits<float>::min()) {
b[i] = 0.0;
}
}
}
const Simulator& simulator_;
const GridView gridView_;
ElementMapper elementMapper_;
VertexMapper vertexMapper_;
EclOutputBlackOilModule<TypeTag> eclOutputModule_;
CollectDataToIORankType collectToIORank_;
std::unique_ptr<Opm::EclipseIO> eclIO_;
double curTime_;
unsigned reportStepIdx_;
std::list<std::pair<std::string, ScalarBuffer*> > attachedBuffers_;
};
} // namespace Ewoms