diff --git a/ebos/collecttoiorank.hh b/ebos/collecttoiorank.hh index 8c2f80496..6d4f94722 100644 --- a/ebos/collecttoiorank.hh +++ b/ebos/collecttoiorank.hh @@ -23,6 +23,9 @@ #ifndef EWOMS_PARALLELSERIALOUTPUT_HH #define EWOMS_PARALLELSERIALOUTPUT_HH +#include +#include + //#if HAVE_OPM_GRID #include #include @@ -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& 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 EquilElementMapper; + EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout()); +#else + typedef Dune::MultipleCodimMultipleGeomTypeMapper 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 EquilElementMapper; - EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout()); -#else - typedef Dune::MultipleCodimMultipleGeomTypeMapper 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(), - 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 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(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 - 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= stride * size ); + for( unsigned int i=0; i - 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 - 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 diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index e075aca01..ad96ece9b 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -28,13 +28,15 @@ #define EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH #include "eclwriter.hh" -#include "ecldeckunits.hh" -#include #include #include #include +#include +#include +#include + #include @@ -42,26 +44,9 @@ namespace Ewoms { namespace Properties { -// create new type tag for the VTK multi-phase output +// create new type tag for the Ecl-output NEW_TYPE_TAG(EclOutputBlackOil); -// create the property tags needed for the multi phase module -NEW_PROP_TAG(EclOutputWriteSaturations); -NEW_PROP_TAG(EclOutputWritePressures); -NEW_PROP_TAG(EclOutputWriteGasDissolutionFactor); -NEW_PROP_TAG(EclOutputWriteOilVaporizationFactor); -NEW_PROP_TAG(EclOutputWriteGasFormationVolumeFactor); -NEW_PROP_TAG(EclOutputWriteOilFormationVolumeFactor); -NEW_PROP_TAG(EclOutputWriteOilSaturationPressure); - -// set default values for what quantities to output -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteSaturations, true); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWritePressures, true); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasDissolutionFactor, true); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilVaporizationFactor, false); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasFormationVolumeFactor, true); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilFormationVolumeFactor, true); -SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilSaturationPressure, false); } // namespace Properties // forward declaration @@ -75,20 +60,18 @@ class EcfvDiscretization; * ECL binary format. */ template -class EclOutputBlackOilModule : public BaseOutputModule +class EclOutputBlackOilModule { - typedef BaseOutputModule ParentType; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef Ewoms::EclWriter EclWriter; - enum { numPhases = FluidSystem::numPhases }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; @@ -96,68 +79,164 @@ class EclOutputBlackOilModule : public BaseOutputModule enum { gasCompIdx = FluidSystem::gasCompIdx }; enum { oilCompIdx = FluidSystem::oilCompIdx }; - typedef typename ParentType::ScalarBuffer ScalarBuffer; + typedef std::vector ScalarBuffer; public: EclOutputBlackOilModule(const Simulator& simulator) - : ParentType(simulator) + : simulator_(simulator) { } - /*! - * \brief Register all run-time parameters for the multi-phase VTK output - * module. - */ - static void registerParameters() - { - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteSaturations, - "Include the saturations of all fluid phases in the " - "ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWritePressures, - "Include the absolute pressures of all fluid phases in the " - "ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor, - "Include the gas dissolution factor of saturated oil in the " - "ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilVaporizationFactor, - "Include the oil vaporization factor of saturated gas in the " - "ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor, - "Include the gas formation volume factor in the " - "ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilFormationVolumeFactor, - "Include the oil formation volume factor of saturated oil " - "in the ECL output files"); - EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure, - "Include the saturation pressure of oil in the " - "ECL output files"); - } - /*! * \brief Allocate memory for the scalar fields we would like to - * write to the VTK file. + * write to ECL output files */ - void allocBuffers() + void allocBuffers(unsigned bufferSize, unsigned reportStepNum, const Opm::RestartConfig& restartConfig, const bool substep, const bool log) { + if (!std::is_same >::value) return; - auto bufferType = ParentType::ElementBuffer; - if (saturationsOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - this->resizeScalarBuffer_(saturation_[phaseIdx], bufferType); + std::map rstKeywords = restartConfig.getRestartKeywords(reportStepNum); + for (auto& keyValue : rstKeywords) { + keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); } - if (pressuresOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - this->resizeScalarBuffer_(pressure_[phaseIdx], bufferType); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + saturation_[phaseIdx].resize(bufferSize,0.0); } - if (gasDissolutionFactorOutput_()) - this->resizeScalarBuffer_(gasDissolutionFactor_, bufferType); - if (oilVaporizationFactorOutput_()) - this->resizeScalarBuffer_(oilVaporizationFactor_, bufferType); - if (gasFormationVolumeFactorOutput_()) - this->resizeScalarBuffer_(gasFormationVolumeFactor_, bufferType); - if (oilSaturationPressureOutput_()) - this->resizeScalarBuffer_(oilSaturationPressure_, bufferType); + + oilPressure_.resize(bufferSize,0.0); + temperature_.resize(bufferSize,0.0); + + // Output the same as legacy + // TODO: Only needed if DISGAS or VAPOIL + if (true) + rs_.resize(bufferSize,0.0); + if (true) + rv_.resize(bufferSize,0.0); + + if (GET_PROP_VALUE(TypeTag, EnableSolvent)) { + sSol_.resize(bufferSize,0.0); + } + if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { + cPolymer_.resize(bufferSize,0.0); + } + + // Output the same as legacy + // TODO: Only needed if Vappars or hysteresis. + soMax_.resize(bufferSize,0.0); + pcSwMdcOw_.resize(bufferSize,0.0); + krnSwMdcOw_.resize(bufferSize,0.0); + pcSwMdcGo_.resize(bufferSize,0.0); + krnSwMdcGo_.resize(bufferSize,0.0); + + // Only provide RESTART_AUXILIARY if it is asked for by the user + if (!restartConfig.getWriteRestartFile(reportStepNum) || substep) + return; + + // Output the same as legacy + // TODO: Only needed if DISGAS or VAPOIL + if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) { + if (rstKeywords["RSSAT"] > 0) { + rstKeywords["RSSAT"] = 0; + gasDissolutionFactor_.resize(bufferSize,0.0); + } + if (rstKeywords["RVSAT"] > 0) { + rstKeywords["RVSAT"] = 0; + oilVaporizationFactor_.resize(bufferSize,0.0); + } + } + + if (FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["BW"] > 0) + { + rstKeywords["BW"] = 0; + invB_[waterPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(oilPhaseIdx) && rstKeywords["BO"] > 0) + { + rstKeywords["BO"] = 0; + invB_[oilPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx) && rstKeywords["BG"] > 0) + { + rstKeywords["BG"] = 0; + invB_[gasPhaseIdx].resize(bufferSize,0.0); + } + + if (rstKeywords["DEN"] > 0) { + rstKeywords["DEN"] = 0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + density_[phaseIdx].resize(bufferSize,0.0); + } + } + const bool hasVWAT = (rstKeywords["VISC"] > 0) || (rstKeywords["VWAT"] > 0); + const bool hasVOIL = (rstKeywords["VISC"] > 0) || (rstKeywords["VOIL"] > 0); + const bool hasVGAS = (rstKeywords["VISC"] > 0) || (rstKeywords["VGAS"] > 0); + rstKeywords["VISC"] = 0; + + if (FluidSystem::phaseIsActive(waterPhaseIdx) && hasVWAT) + { + rstKeywords["VWAT"] = 0; + viscosity_[waterPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(oilPhaseIdx) && hasVOIL > 0) + { + rstKeywords["VOIL"] = 0; + viscosity_[oilPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx) && hasVGAS > 0) + { + rstKeywords["VGAS"] = 0; + viscosity_[gasPhaseIdx].resize(bufferSize,0.0); + } + + if (FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["KRW"] > 0) + { + rstKeywords["KRW"] = 0; + relativePermeability_[waterPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(oilPhaseIdx) && rstKeywords["KRO"] > 0) + { + rstKeywords["KRO"] = 0; + relativePermeability_[oilPhaseIdx].resize(bufferSize,0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx) && rstKeywords["KRG"] > 0) + { + rstKeywords["KRG"] = 0; + relativePermeability_[gasPhaseIdx].resize(bufferSize,0.0); + } + + if (rstKeywords["PBPD"] > 0) { + rstKeywords["PBPD"] = 0; + bubblePointPressure_.resize(bufferSize,0.0); + dewPointPressure_.resize(bufferSize,0.0); + } + + //Warn for any unhandled keyword + if (log) { + for (auto& keyValue : rstKeywords) { + if (keyValue.second > 0) { + std::string logstring = "Keyword '"; + logstring.append(keyValue.first); + logstring.append("' is unhandled for output to file."); + Opm::OpmLog::warning("Unhandled output keyword", logstring); + } + } + } + + failedCellsPb_.clear(); + failedCellsPd_.clear(); + + // Not supported in flow legacy + if (false) + saturatedOilFormationVolumeFactor_.resize(bufferSize,0.0); + if (false) + oilSaturationPressure_.resize(bufferSize,0.0); } /*! @@ -166,158 +245,556 @@ public: */ void processElement(const ElementContext& elemCtx) { - if (!EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput)) - return; - typedef Opm::MathToolbox Toolbox; if (!std::is_same >::value) return; for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { - const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState(); + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); typedef typename std::remove_const::type>::type FluidState; unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(); - if (saturationsOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (saturation_[phaseIdx].size() == 0) + continue; - saturation_[phaseIdx][globalDofIdx] = Toolbox::value(fs.saturation(phaseIdx)); - Opm::Valgrind::CheckDefined(saturation_[phaseIdx][globalDofIdx]); - } + saturation_[phaseIdx][globalDofIdx] = Toolbox::value(fs.saturation(phaseIdx)); + Opm::Valgrind::CheckDefined(saturation_[phaseIdx][globalDofIdx]); } - if (pressuresOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - pressure_[phaseIdx][globalDofIdx] = Toolbox::value(fs.pressure(phaseIdx)); - Opm::Valgrind::CheckDefined(pressure_[phaseIdx][globalDofIdx]); - } + if (oilPressure_.size() > 0) { + oilPressure_[globalDofIdx] = Toolbox::value(fs.pressure(oilPhaseIdx)); + Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]); + } - if (gasDissolutionFactorOutput_()) { - Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx); - gasDissolutionFactor_[globalDofIdx] = - FluidSystem::template saturatedDissolutionFactor(fs, gasPhaseIdx, pvtRegionIdx, SoMax); - Opm::Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]); - } - if (oilVaporizationFactorOutput_()) { + if (gasDissolutionFactor_.size() > 0) { Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx); gasDissolutionFactor_[globalDofIdx] = FluidSystem::template saturatedDissolutionFactor(fs, oilPhaseIdx, pvtRegionIdx, SoMax); Opm::Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]); + } - if (gasFormationVolumeFactorOutput_()) { + if (oilVaporizationFactor_.size() > 0) { + Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx); + oilVaporizationFactor_[globalDofIdx] = + FluidSystem::template saturatedDissolutionFactor(fs, gasPhaseIdx, pvtRegionIdx, SoMax); + Opm::Valgrind::CheckDefined(oilVaporizationFactor_[globalDofIdx]); + + } + if (gasFormationVolumeFactor_.size() > 0) { gasFormationVolumeFactor_[globalDofIdx] = 1.0/FluidSystem::template inverseFormationVolumeFactor(fs, gasPhaseIdx, pvtRegionIdx); Opm::Valgrind::CheckDefined(gasFormationVolumeFactor_[globalDofIdx]); + } - if (oilSaturationPressureOutput_()) { + if (saturatedOilFormationVolumeFactor_.size() > 0) { + saturatedOilFormationVolumeFactor_[globalDofIdx] = + 1.0/FluidSystem::template saturatedInverseFormationVolumeFactor(fs, oilPhaseIdx, pvtRegionIdx); + Opm::Valgrind::CheckDefined(saturatedOilFormationVolumeFactor_[globalDofIdx]); + + } + if (oilSaturationPressure_.size() > 0) { oilSaturationPressure_[globalDofIdx] = FluidSystem::template saturationPressure(fs, oilPhaseIdx, pvtRegionIdx); Opm::Valgrind::CheckDefined(oilSaturationPressure_[globalDofIdx]); + } + + if (rs_.size()) { + rs_[globalDofIdx] = Toolbox::value(fs.Rs()); + Opm::Valgrind::CheckDefined(rs_[globalDofIdx]); + } + + if (rv_.size()) { + rv_[globalDofIdx] = Toolbox::value(fs.Rv()); + Opm::Valgrind::CheckDefined(rv_[globalDofIdx]); + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (invB_[phaseIdx].size() == 0) + continue; + + invB_[phaseIdx][globalDofIdx] = Toolbox::value(fs.invB(phaseIdx)); + Opm::Valgrind::CheckDefined(invB_[phaseIdx][globalDofIdx]); + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (density_[phaseIdx].size() == 0) + continue; + + density_[phaseIdx][globalDofIdx] = Toolbox::value(fs.density(phaseIdx)); + Opm::Valgrind::CheckDefined(density_[phaseIdx][globalDofIdx]); + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (viscosity_[phaseIdx].size() == 0) + continue; + + viscosity_[phaseIdx][globalDofIdx] = Toolbox::value(fs.viscosity(phaseIdx)); + Opm::Valgrind::CheckDefined(viscosity_[phaseIdx][globalDofIdx]); + } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (relativePermeability_[phaseIdx].size() == 0) + continue; + + relativePermeability_[phaseIdx][globalDofIdx] = Toolbox::value(intQuants.relativePermeability(phaseIdx)); + Opm::Valgrind::CheckDefined(relativePermeability_[phaseIdx][globalDofIdx]); + } + + if (sSol_.size() > 0) { + sSol_[globalDofIdx] = intQuants.solventSaturation().value(); + } + + if (cPolymer_.size() > 0) { + cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value(); + } + + if (bubblePointPressure_.size() > 0) + { + try { + bubblePointPressure_[globalDofIdx] = Toolbox::value(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex())); + } + catch (const Opm::NumericalProblem& e) { + const auto globalIdx = elemCtx.simulator().gridManager().grid().globalCell()[globalDofIdx]; + failedCellsPb_.push_back(globalIdx); + } + } + if (dewPointPressure_.size() > 0) + { + try { + dewPointPressure_[globalDofIdx] = Toolbox::value(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex())); + } + catch (const Opm::NumericalProblem& e) { + const auto globalIdx = elemCtx.simulator().gridManager().grid().globalCell()[globalDofIdx]; + failedCellsPd_.push_back(globalIdx); + } + } + + if (soMax_.size() > 0) + soMax_[globalDofIdx] = elemCtx.simulator().model().maxOilSaturation(globalDofIdx); + + const auto& matLawManager = elemCtx.simulator().problem().materialLawManager(); + if (matLawManager->enableHysteresis()) { + if (pcSwMdcOw_.size() > 0 && krnSwMdcOw_.size() > 0) { + matLawManager->oilWaterHysteresisParams( + pcSwMdcOw_[globalDofIdx], + krnSwMdcOw_[globalDofIdx], + globalDofIdx); + } + if (pcSwMdcGo_.size() > 0 && krnSwMdcGo_.size() > 0) { + matLawManager->gasOilHysteresisParams( + pcSwMdcGo_[globalDofIdx], + krnSwMdcGo_[globalDofIdx], + globalDofIdx); + } + } + + // hack to make the intial output of rs and rv Ecl compatible. + // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step + // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite + // rs and rv with the values computed in the initially. + // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values. + // This can be removed when ebos has 100% controll over output + if (elemCtx.simulator().episodeIndex() < 0 && FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx) ) { + + const auto& fs_initial = elemCtx.simulator().problem().initialFluidState(globalDofIdx); + + // use initial rs and rv values + if (rv_.size() > 0) + rv_[globalDofIdx] = fs_initial.Rv(); + + if (rs_.size() > 0) + rs_[globalDofIdx] = fs_initial.Rs(); + + // re-compute the volume factors, viscosities and densities if asked for + if (density_[oilPhaseIdx].size() > 0) + density_[oilPhaseIdx][globalDofIdx] = FluidSystem::density(fs_initial, + oilPhaseIdx, + intQuants.pvtRegionIndex()); + if (density_[gasPhaseIdx].size() > 0) + density_[gasPhaseIdx][globalDofIdx] = FluidSystem::density(fs_initial, + gasPhaseIdx, + intQuants.pvtRegionIndex()); + + if (invB_[oilPhaseIdx].size() > 0) + invB_[oilPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fs_initial, + oilPhaseIdx, + intQuants.pvtRegionIndex()); + if (invB_[gasPhaseIdx].size() > 0) + invB_[gasPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fs_initial, + gasPhaseIdx, + intQuants.pvtRegionIndex()); + if (viscosity_[oilPhaseIdx].size() > 0) + viscosity_[oilPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fs_initial, + oilPhaseIdx, + intQuants.pvtRegionIndex()); + if (viscosity_[gasPhaseIdx].size() > 0) + viscosity_[gasPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fs_initial, + gasPhaseIdx, + intQuants.pvtRegionIndex()); + } + + } + + } + + + void outputErrorLog() + { + const size_t maxNumCellsFaillog = 20; + + int pbSize = failedCellsPb_.size(), pd_size = failedCellsPd_.size(); + std::vector displPb, displPd, recvLenPb, recvLenPd; + const auto& comm = simulator_.gridView().comm(); + + if ( isIORank_() ) + { + displPb.resize(comm.size()+1, 0); + displPd.resize(comm.size()+1, 0); + recvLenPb.resize(comm.size()); + recvLenPd.resize(comm.size()); + } + + comm.gather(&pbSize, recvLenPb.data(), 1, 0); + comm.gather(&pd_size, recvLenPd.data(), 1, 0); + std::partial_sum(recvLenPb.begin(), recvLenPb.end(), displPb.begin()+1); + std::partial_sum(recvLenPd.begin(), recvLenPd.end(), displPd.begin()+1); + std::vector globalFailedCellsPb, globalFailedCellsPd; + + if ( isIORank_() ) + { + globalFailedCellsPb.resize(displPb.back()); + globalFailedCellsPd.resize(displPd.back()); + } + + comm.gatherv(failedCellsPb_.data(), static_cast(failedCellsPb_.size()), + globalFailedCellsPb.data(), recvLenPb.data(), + displPb.data(), 0); + comm.gatherv(failedCellsPd_.data(), static_cast(failedCellsPd_.size()), + globalFailedCellsPd.data(), recvLenPd.data(), + displPd.data(), 0); + std::sort(globalFailedCellsPb.begin(), globalFailedCellsPb.end()); + std::sort(globalFailedCellsPd.begin(), globalFailedCellsPd.end()); + + if (globalFailedCellsPb.size() > 0) { + std::stringstream errlog; + errlog << "Finding the bubble point pressure failed for " << globalFailedCellsPb.size() << " cells ["; + errlog << globalFailedCellsPb[0]; + const size_t max_elems = std::min(maxNumCellsFaillog, globalFailedCellsPb.size()); + for (size_t i = 1; i < max_elems; ++i) { + errlog << ", " << globalFailedCellsPb[i]; + } + if (globalFailedCellsPb.size() > maxNumCellsFaillog) { + errlog << ", ..."; + } + errlog << "]"; + Opm::OpmLog::warning("Bubble point numerical problem", errlog.str()); + } + if (globalFailedCellsPd.size() > 0) { + std::stringstream errlog; + errlog << "Finding the dew point pressure failed for " << globalFailedCellsPd.size() << " cells ["; + errlog << globalFailedCellsPd[0]; + const size_t max_elems = std::min(maxNumCellsFaillog, globalFailedCellsPd.size()); + for (size_t i = 1; i < max_elems; ++i) { + errlog << ", " << globalFailedCellsPd[i]; + } + if (globalFailedCellsPd.size() > maxNumCellsFaillog) { + errlog << ", ..."; + } + errlog << "]"; + Opm::OpmLog::warning("Dew point numerical problem", errlog.str()); } } /*! - * \brief Add all buffers to the VTK output writer. + * \brief Move all buffers to data::Solution. */ - void commitBuffers(BaseOutputWriter& writer) + void assignToSolution(Opm::data::Solution& sol) { if (!std::is_same >::value) return; - if (!dynamic_cast(&writer)) - return; // this module only consideres ecl writers... + if ( oilPressure_.size() > 0 ) { + sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, std::move(oilPressure_), Opm::data::TargetType::RESTART_SOLUTION); + } - typedef EclDeckUnits DeckUnits; - const DeckUnits& deckUnits = this->simulator_.problem().deckUnits(); + if ( temperature_.size() > 0 ) { + sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, std::move(temperature_), Opm::data::TargetType::RESTART_SOLUTION); + } - typename ParentType::BufferType bufferType = ParentType::ElementBuffer; - if (pressuresOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - deckUnits.siToDeck(pressure_[phaseIdx], DeckUnits::pressure); + if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size() > 0 ) { + sol.insert( "SWAT", Opm::UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION ); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size() > 0) { + sol.insert( "SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION ); + } - this->commitScalarBuffer_(writer, "PRESSURE", pressure_[oilPhaseIdx], bufferType); - this->commitScalarBuffer_(writer, "PGAS", pressure_[gasPhaseIdx], bufferType); - this->commitScalarBuffer_(writer, "PWAT", pressure_[waterPhaseIdx], bufferType); - } - if (saturationsOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - deckUnits.siToDeck(saturation_[phaseIdx], DeckUnits::saturation); + if ( gasDissolutionFactor_.size() > 0 ) { + sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, std::move(gasDissolutionFactor_), Opm::data::TargetType::RESTART_AUXILIARY ); - this->commitScalarBuffer_(writer, "SWAT", saturation_[waterPhaseIdx], bufferType); - this->commitScalarBuffer_(writer, "SGAS", saturation_[gasPhaseIdx], bufferType); - // the oil saturation is _NOT_ written to disk. Instead, it is calculated by - // the visualization tool. Wondering why is probably a waste of time... } - if (gasDissolutionFactorOutput_()) { - deckUnits.siToDeck(gasDissolutionFactor_, DeckUnits::gasDissolutionFactor); - this->commitScalarBuffer_(writer, "RS", gasDissolutionFactor_, bufferType); + if ( oilVaporizationFactor_.size() > 0 ) { + sol.insert( "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, std::move(oilVaporizationFactor_) , Opm::data::TargetType::RESTART_AUXILIARY ); } - if (oilVaporizationFactorOutput_()) { - deckUnits.siToDeck(oilVaporizationFactor_, DeckUnits::oilVaporizationFactor); - this->commitScalarBuffer_(writer, "RV", oilVaporizationFactor_, bufferType); + if ( rs_.size() > 0 ) { + sol.insert( "RS", Opm::UnitSystem::measure::gas_oil_ratio, std::move(rs_), Opm::data::TargetType::RESTART_SOLUTION ); + } - if (gasFormationVolumeFactorOutput_()) { - // no unit conversion required - this->commitScalarBuffer_(writer, "BG", gasFormationVolumeFactor_, bufferType); + if (rv_.size() > 0 ) { + sol.insert( "RV", Opm::UnitSystem::measure::oil_gas_ratio, std::move(rv_) , Opm::data::TargetType::RESTART_SOLUTION ); } - if (oilSaturationPressureOutput_()) { - deckUnits.siToDeck(oilSaturationPressure_, DeckUnits::pressure); - this->commitScalarBuffer_(writer, "PSAT", oilSaturationPressure_); + if (invB_[waterPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, std::move(invB_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + if (invB_[oilPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, std::move(invB_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + if (invB_[gasPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, std::move(invB_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + + if (density_[waterPhaseIdx].size() > 0 ) { + sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, std::move(density_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + if (density_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, std::move(density_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + if (density_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, std::move(density_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + + if (viscosity_[waterPhaseIdx].size() > 0 ) { + sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); + } + if (viscosity_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); + } + if (viscosity_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); + } + + if (relativePermeability_[waterPhaseIdx].size() > 0) { + sol.insert( "WATKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); + } + if (relativePermeability_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OILKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); + } + if (relativePermeability_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GASKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); + } + + if (pcSwMdcOw_.size() > 0 ) + sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, std::move(pcSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (krnSwMdcOw_.size() > 0) + sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, std::move(krnSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (pcSwMdcGo_.size() > 0) + sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, std::move(pcSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (krnSwMdcGo_.size() > 0) + sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, std::move(krnSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (soMax_.size() > 0) + sol.insert ("SOMAX", Opm::UnitSystem::measure::identity, std::move(soMax_), Opm::data::TargetType::RESTART_SOLUTION); + + if (sSol_.size() > 0) + sol.insert ("SSOL", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION); + + if (cPolymer_.size() > 0) + sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION); + + if (dewPointPressure_.size() > 0) + sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, std::move(dewPointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); + + if (bubblePointPressure_.size() > 0) + sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); + +} + + void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex) + { + + Scalar so = 1.0; + if( sol.has( "SWAT" ) ) { + saturation_[waterPhaseIdx][elemIdx] = sol.data("SWAT")[globalDofIndex]; + so -= sol.data("SWAT")[globalDofIndex]; + } + + if( sol.has( "SGAS" ) ) { + saturation_[gasPhaseIdx][elemIdx] = sol.data("SGAS")[globalDofIndex]; + so -= sol.data("SGAS")[globalDofIndex]; + } + + saturation_[oilPhaseIdx][elemIdx] = so; + + if( sol.has( "PRESSURE" ) ) { + oilPressure_[elemIdx] = sol.data( "PRESSURE" )[globalDofIndex]; + } + + if( sol.has( "TEMP" ) ) { + temperature_[elemIdx] = sol.data( "TEMP" )[globalDofIndex]; + } + + if( sol.has( "RS" ) ) { + rs_[elemIdx] = sol.data("RS")[globalDofIndex]; + } + + if( sol.has( "RV" ) ) { + rv_[elemIdx] = sol.data("RV")[globalDofIndex]; + } + + if ( sol.has( "SSOL" ) ) { + sSol_[elemIdx] = sol.data("SSOL")[globalDofIndex]; + } + + if ( sol.has("POLYMER" ) ) { + cPolymer_[elemIdx] = sol.data("POLYMER")[globalDofIndex]; + } + + if ( sol.has("SOMAX" ) ) { + soMax_[elemIdx] = sol.data("SOMAX")[globalDofIndex]; + } + + if ( sol.has("PCSWM_OW" ) ) { + pcSwMdcOw_[elemIdx] = sol.data("PCSWM_OW")[globalDofIndex]; + } + + if ( sol.has("KRNSW_OW" ) ) { + krnSwMdcOw_[elemIdx] = sol.data("KRNSW_OW")[globalDofIndex]; + } + + if ( sol.has("PCSWM_GO" ) ) { + pcSwMdcGo_[elemIdx] = sol.data("PCSWM_GO")[globalDofIndex]; + } + + if ( sol.has("KRNSW_GO" ) ) { + krnSwMdcGo_[elemIdx] = sol.data("KRNSW_GO")[globalDofIndex]; } } + + template + void assignToFluidState(FluidState& fs, unsigned elemIdx) const + { + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (saturation_[phaseIdx].size() == 0) + continue; + + fs.setSaturation(phaseIdx, saturation_[phaseIdx][elemIdx]); + } + + if (oilPressure_.size() > 0) { + + // this assumes that capillary pressures only depend on the phase saturations + // and possibly on temperature. (this is always the case for ECL problems.) + Dune::FieldVector< Scalar, numPhases > pc( 0 ); + const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx); + MaterialLaw::capillaryPressures(pc, matParams, fs); + Opm::Valgrind::CheckDefined(oilPressure_[elemIdx]); + Opm::Valgrind::CheckDefined(pc); + assert(FluidSystem::phaseIsActive(oilPhaseIdx)); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + fs.setPressure(phaseIdx, oilPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx])); + } + } + + if (temperature_.size() > 0) { + fs.setTemperature( temperature_[elemIdx]); + } + + if (rs_.size() > 0) { + fs.setRs(rs_[elemIdx]); + + } + if (rv_.size() > 0) { + fs.setRv(rv_[elemIdx]); + } + } + + void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const { + + if (soMax_.size() > 0) + simulator.model().setMaxOilSaturation(soMax_[elemIdx], elemIdx); + + if (simulator.problem().materialLawManager()->enableHysteresis()) { + auto matLawManager = simulator.problem().materialLawManager(); + + if (pcSwMdcOw_.size() > 0 && krnSwMdcOw_.size() > 0) { + matLawManager->setOilWaterHysteresisParams( + pcSwMdcOw_[elemIdx], + krnSwMdcOw_[elemIdx], + elemIdx); + } + if (pcSwMdcGo_.size() > 0 && krnSwMdcGo_.size() > 0) { + matLawManager->setGasOilHysteresisParams( + pcSwMdcGo_[elemIdx], + krnSwMdcGo_[elemIdx], + elemIdx); + } + } + + } + + Scalar getSolventSaturation(unsigned elemIdx) const { + if(sSol_.size() > 0) + return sSol_[elemIdx]; + + return 0; + } + + Scalar getPolymerConcentration(unsigned elemIdx) const { + if(cPolymer_.size() > 0) + return cPolymer_[elemIdx]; + + return 0; + } + + private: - static bool saturationsOutput_() - { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteSaturations); } - static bool pressuresOutput_() - { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWritePressures); } - - static bool gasDissolutionFactorOutput_() + bool isIORank_() const { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor); + const auto& comm = simulator_.gridView().comm(); + return comm.rank() == 0; } - static bool gasFormationVolumeFactorOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor); - } - - static bool oilVaporizationFactorOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilVaporizationFactor); - } - - static bool oilSaturationPressureOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure); - } + const Simulator& simulator_; ScalarBuffer saturation_[numPhases]; - ScalarBuffer pressure_[numPhases]; + ScalarBuffer oilPressure_; + ScalarBuffer temperature_; ScalarBuffer gasDissolutionFactor_; ScalarBuffer oilVaporizationFactor_; ScalarBuffer gasFormationVolumeFactor_; + ScalarBuffer saturatedOilFormationVolumeFactor_; ScalarBuffer oilSaturationPressure_; + ScalarBuffer rs_; + ScalarBuffer rv_; + ScalarBuffer invB_[numPhases]; + ScalarBuffer density_[numPhases]; + ScalarBuffer viscosity_[numPhases]; + ScalarBuffer relativePermeability_[numPhases]; + ScalarBuffer sSol_; + ScalarBuffer cPolymer_; + ScalarBuffer soMax_; + ScalarBuffer pcSwMdcOw_; + ScalarBuffer krnSwMdcOw_; + ScalarBuffer pcSwMdcGo_; + ScalarBuffer krnSwMdcGo_; + ScalarBuffer bubblePointPressure_; + ScalarBuffer dewPointPressure_; + std::vector failedCellsPb_; + std::vector failedCellsPd_; + }; } // namespace Ewoms diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 104e6b2a2..cac910c2f 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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 #include @@ -80,6 +78,7 @@ #include #include #include +#include #include #include @@ -87,6 +86,8 @@ #include #include +#include + #include #include @@ -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 Toolbox; - typedef Ewoms::EclSummaryWriter EclSummaryWriter; typedef Dune::FieldMatrix DimMatrix; typedef EclWriter EclWriterType; + typedef typename GridView::template Codim<0>::Iterator ElementIterator; + + struct RockParams { Scalar referencePressure; Scalar compressibility; @@ -310,8 +310,6 @@ public: { ParentType::registerParameters(); - Ewoms::EclOutputBlackOilModule::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(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& 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&& 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 wellManager_; - EclDeckUnits deckUnits_; - std::unique_ptr< EclWriterType > eclWriter_; - EclSummaryWriter summaryWriter_; PffGridVector pffDofData_; + + bool restartApplied; + }; } // namespace Ewoms diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 1955a0174..96ce8354e 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -30,11 +30,12 @@ #include -#include "ertwrappers.hh" #include "collecttoiorank.hh" +#include "ecloutputblackoilmodule.hh" #include #include +#include #include #include @@ -58,302 +59,170 @@ NEW_PROP_TAG(EnableEclOutput); template class EclWriter; -template -class EclWriterHelper -{ - friend class EclWriter; - - static void writeHeaders_(EclWriter& writer) - { - typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization; - if (!std::is_same >::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 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 ScalarBuffer; - friend class EclWriterHelper; 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&& 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::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(), substep, 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(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(); + const ElementIterator& elemEndIt = gridView.template end(); + 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); + if (collectToIORank_.isParallel()) + 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> extraRestartData; + std::map miscSummaryData; - std::shared_ptr> - bufKeyword(new ErtKeyword(name, buffer)); - solution.add(bufKeyword); + // Add suggested next timestep to extra data. + extraRestartData["OPMEXTRA"] = std::vector(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 - void serialize(Restarter& res) - { - res.serializeSectionBegin("EclWriter"); - res.serializeStream() << reportStepIdx_ << "\n"; - res.serializeSectionEnd(); + void restartBegin() { + std::map 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 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(), true, 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 - void deserialize(Restarter& res) - { - res.deserializeSectionBegin("EclWriter"); - res.deserializeStream() >> reportStepIdx_; - res.deserializeSectionEnd(); + + const EclOutputBlackOilModule& 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& 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::min()) { - b[i] = 0.0; - } - } - } - const Simulator& simulator_; - const GridView gridView_; - - ElementMapper elementMapper_; - VertexMapper vertexMapper_; + EclOutputBlackOilModule eclOutputModule_; CollectDataToIORankType collectToIORank_; + std::unique_ptr eclIO_; - double curTime_; - unsigned reportStepIdx_; - - std::list > attachedBuffers_; }; } // namespace Ewoms