From 0132c3326ed5c320a7e91d3bcf542cbe6e17456b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 7 Dec 2017 12:38:48 +0100 Subject: [PATCH 1/4] 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. --- ebos/collecttoiorank.hh | 243 ++++----- ebos/ecloutputblackoilmodule.hh | 856 +++++++++++++++++++++++++++----- ebos/eclproblem.hh | 130 +++-- ebos/eclwriter.hh | 334 ++++--------- 4 files changed, 1036 insertions(+), 527 deletions(-) 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..282958f6f 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,158 @@ 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 log) { + std::map rstKeywords = restartConfig.getRestartKeywords(reportStepNum); + for (auto& keyValue : rstKeywords) { + keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); + } + if (!std::is_same >::value) return; - auto bufferType = ParentType::ElementBuffer; if (saturationsOutput_()) { for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - this->resizeScalarBuffer_(saturation_[phaseIdx], bufferType); + saturation_[phaseIdx].resize(bufferSize,0.0); } if (pressuresOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - this->resizeScalarBuffer_(pressure_[phaseIdx], bufferType); + oilPressure_.resize(bufferSize,0.0); } - if (gasDissolutionFactorOutput_()) - this->resizeScalarBuffer_(gasDissolutionFactor_, bufferType); - if (oilVaporizationFactorOutput_()) - this->resizeScalarBuffer_(oilVaporizationFactor_, bufferType); + if (temperatureOutput_()) { + temperature_.resize(bufferSize,0.0); + } + + if (gasDissolutionFactorOutput_() && rstKeywords["RSSAT"] > 0) { + rstKeywords["RSSAT"] = 0; + gasDissolutionFactor_.resize(bufferSize,0.0); + } + if (oilVaporizationFactorOutput_() && rstKeywords["RVSAT"] > 0) { + rstKeywords["RVSAT"] = 0; + oilVaporizationFactor_.resize(bufferSize,0.0); + } + if (gasFormationVolumeFactorOutput_()) - this->resizeScalarBuffer_(gasFormationVolumeFactor_, bufferType); + gasFormationVolumeFactor_.resize(bufferSize,0.0); + if (saturatedOilFormationVolumeFactorOutput_()) + saturatedOilFormationVolumeFactor_.resize(bufferSize,0.0); if (oilSaturationPressureOutput_()) - this->resizeScalarBuffer_(oilSaturationPressure_, bufferType); + oilSaturationPressure_.resize(bufferSize,0.0); + + if (rsOutput_()) + rs_.resize(bufferSize,0.0); + if (rvOutput_()) + rv_.resize(bufferSize,0.0); + if (invBOutput_()) { + 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 (densityOutput_() && rstKeywords["DEN"] > 0) { + rstKeywords["DEN"] = 0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + density_[phaseIdx].resize(bufferSize,0.0); + } + } + if (viscosityOutput_()) { + 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 (relativePermeabilityOutput_()) { + 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 (solventOutput_()) { + sSol_.resize(bufferSize,0.0); + } + if (polymerOutput_()) { + cPolymer_.resize(bufferSize,0.0); + } + + // TODO: Only needed if Vappars or hysteresis. + // Now: Output the same as legacy + 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); + + 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(); } /*! @@ -166,158 +239,669 @@ 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)) + if (saturation_[phaseIdx].size() == 0) continue; 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; + if (oilPressure_.size() > 0) { + oilPressure_[globalDofIdx] = Toolbox::value(fs.pressure(oilPhaseIdx)); + Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]); - pressure_[phaseIdx][globalDofIdx] = Toolbox::value(fs.pressure(phaseIdx)); - Opm::Valgrind::CheckDefined(pressure_[phaseIdx][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]); + } + + if (invBOutput_()) { + 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]); + } + } + + if (densityOutput_()) { + 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]); + } + } + + if (viscosityOutput_()) { + 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]); + } + } + + if (relativePermeabilityOutput_()) { + 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); + + if (hysteresisOutput_()) { + const auto& matLawManager = elemCtx.simulator().problem().materialLawManager(); + if (matLawManager->enableHysteresis()) { + matLawManager->oilWaterHysteresisParams( + pcSwMdcOw_[globalDofIdx], + krnSwMdcOw_[globalDofIdx], + globalDofIdx); + 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 + rv_[globalDofIdx] = fs_initial.Rv(); + 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 Add 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... - - typedef EclDeckUnits DeckUnits; - const DeckUnits& deckUnits = this->simulator_.problem().deckUnits(); - - typename ParentType::BufferType bufferType = ParentType::ElementBuffer; - if (pressuresOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - deckUnits.siToDeck(pressure_[phaseIdx], DeckUnits::pressure); - - this->commitScalarBuffer_(writer, "PRESSURE", pressure_[oilPhaseIdx], bufferType); - this->commitScalarBuffer_(writer, "PGAS", pressure_[gasPhaseIdx], bufferType); - this->commitScalarBuffer_(writer, "PWAT", pressure_[waterPhaseIdx], bufferType); + if ( oilPressure_.size() > 0 ) { + sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, oilPressure_, Opm::data::TargetType::RESTART_SOLUTION); } + + if ( temperature_.size() > 0 ) { + sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, temperature_, Opm::data::TargetType::RESTART_SOLUTION); + } + if (saturationsOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - deckUnits.siToDeck(saturation_[phaseIdx], DeckUnits::saturation); + if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size()>0 ) { + sol.insert( "SWAT", Opm::UnitSystem::measure::identity, saturation_[waterPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size()>0) { + sol.insert( "SGAS", Opm::UnitSystem::measure::identity, saturation_[gasPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); + } + } + if ( gasDissolutionFactor_.size() > 0 ) { + sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, 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, 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, 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, rv_ , Opm::data::TargetType::RESTART_SOLUTION ); } - if (oilSaturationPressureOutput_()) { - deckUnits.siToDeck(oilSaturationPressure_, DeckUnits::pressure); - this->commitScalarBuffer_(writer, "PSAT", oilSaturationPressure_); + if (invBOutput_()) { + if( FluidSystem::phaseIsActive(waterPhaseIdx) && invB_[waterPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, invB_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + if( FluidSystem::phaseIsActive(oilPhaseIdx) && invB_[oilPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, invB_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx) && invB_[gasPhaseIdx].size() > 0 ) { + sol.insert( "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, invB_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + } + if (densityOutput_()) { + if( FluidSystem::phaseIsActive(waterPhaseIdx) && density_[waterPhaseIdx].size() > 0 ) { + sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, density_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + if( FluidSystem::phaseIsActive(oilPhaseIdx) && density_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, density_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx) && density_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, density_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + } + if (viscosityOutput_()) { + if( FluidSystem::phaseIsActive(waterPhaseIdx) && viscosity_[waterPhaseIdx].size() > 0 ) { + sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + } + if( FluidSystem::phaseIsActive(oilPhaseIdx) && viscosity_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx) && viscosity_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + } + } + if (relativePermeabilityOutput_()) { + if( FluidSystem::phaseIsActive(waterPhaseIdx) && relativePermeability_[waterPhaseIdx].size() > 0) { + sol.insert( "WATKR", Opm::UnitSystem::measure::identity, relativePermeability_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + } + if( FluidSystem::phaseIsActive(oilPhaseIdx)&& relativePermeability_[oilPhaseIdx].size() > 0 ) { + sol.insert( "OILKR", Opm::UnitSystem::measure::identity, relativePermeability_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + } + if( FluidSystem::phaseIsActive(gasPhaseIdx)&& relativePermeability_[gasPhaseIdx].size() > 0 ) { + sol.insert( "GASKR", Opm::UnitSystem::measure::identity, relativePermeability_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + } + } + + if (hysteresisOutput_()) { + sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, pcSwMdcOw_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, krnSwMdcOw_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, pcSwMdcGo_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, krnSwMdcGo_, Opm::data::TargetType::RESTART_AUXILIARY); + } + + if (soMaxOutput_()) + sol.insert ("SOMAX", Opm::UnitSystem::measure::identity, soMax_, Opm::data::TargetType::RESTART_SOLUTION); + + if (solventOutput_()) + sol.insert ("SSOL", Opm::UnitSystem::measure::identity, sSol_, Opm::data::TargetType::RESTART_SOLUTION); + + if (polymerOutput_()) + sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, cPolymer_, Opm::data::TargetType::RESTART_SOLUTION); + + if (dewPointPressureOutput_() && dewPointPressure_.size() > 0) + sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, dewPointPressure_, Opm::data::TargetType::RESTART_AUXILIARY); + + if (bubbelPointPressureOutput_() && bubblePointPressure_.size() > 0) + sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, 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 + { + + if (saturationsOutput_()) { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + fs.setSaturation(phaseIdx, saturation_[phaseIdx][elemIdx]); + } + } + if (pressuresOutput_()) { + + // 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 (temperatureOutput_()) { + fs.setTemperature( temperature_[elemIdx]); + } + + if (rsOutput_()) { + fs.setRs(rs_[elemIdx]); + + } + if (rvOutput_()) { + fs.setRv(rv_[elemIdx]); + } + } + + void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const { + + if (soMaxOutput_()) + simulator.model().setMaxOilSaturation(soMax_[elemIdx], elemIdx); + + if (simulator.problem().materialLawManager()->enableHysteresis()) { + auto matLawManager = simulator.problem().materialLawManager(); + + matLawManager->setOilWaterHysteresisParams( + pcSwMdcOw_[elemIdx], + krnSwMdcOw_[elemIdx], + elemIdx); + matLawManager->setGasOilHysteresisParams( + pcSwMdcGo_[elemIdx], + krnSwMdcGo_[elemIdx], + elemIdx); + } + + } + + Scalar getSolventSaturation(unsigned elemIdx) const { + if(solventOutput_()) + return sSol_[elemIdx]; + + return 0; + } + + Scalar getPolymerConcentration(unsigned elemIdx) const { + if(polymerOutput_()) + return cPolymer_[elemIdx]; + + return 0; + } + + private: + // This should be cleaned up. + // For now output the same as legacy + // to make the tests pass static bool saturationsOutput_() - { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteSaturations); } + { return true; } static bool pressuresOutput_() - { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWritePressures); } + { return true; } + + static bool temperatureOutput_() + { + return true; + } + + static bool solventOutput_() + { + return GET_PROP_VALUE(TypeTag, EnableSolvent); + } + + static bool polymerOutput_() + { + return GET_PROP_VALUE(TypeTag, EnablePolymer); + } static bool gasDissolutionFactorOutput_() { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor); + return true; + //FluidSystem::enableDissolvedGas(); } static bool gasFormationVolumeFactorOutput_() { return FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor); + FluidSystem::phaseIsActive(gasPhaseIdx); } static bool oilVaporizationFactorOutput_() + { + return true; + //FluidSystem::enableVaporizedOil(); + } + + static bool saturatedOilFormationVolumeFactorOutput_() { return FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilVaporizationFactor); + FluidSystem::phaseIsActive(gasPhaseIdx); } static bool oilSaturationPressureOutput_() { return FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx) && - EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure); + FluidSystem::phaseIsActive(gasPhaseIdx); } + static bool rsOutput_() + { + // Output the same as legacy + return true; //FluidSystem::enableDissolvedGas(); + + } + + static bool rvOutput_() + { + // Output the same as legacy + return true;//FluidSystem::enableVaporizedOil(); + } + + + static bool invBOutput_() + { + return true; + } + + static bool densityOutput_() + { + return true; + } + + static bool viscosityOutput_() + { + return true; + } + + + static bool relativePermeabilityOutput_() + { + return true; + } + + static bool soMaxOutput_() + { + return true; + } + + static bool hysteresisOutput_() + { + return true; + } + + static bool dewPointPressureOutput_() + { + return true; + } + + static bool bubbelPointPressureOutput_() + { + return true; + } + + bool isIORank() const + { + const auto& comm = simulator_.gridView().comm(); + return comm.rank() == 0; + } + + 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..b457f032c 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,169 @@ 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(), 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); + 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(), 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 From 1969c3711da506752af447f6b098c40fa01721c2 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 9 Jan 2018 11:28:20 +0100 Subject: [PATCH 2/4] Move instead of copy data to solution vector --- ebos/ecloutputblackoilmodule.hh | 60 ++++++++++++++++----------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 282958f6f..222cd9d3b 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -509,7 +509,7 @@ public: } /*! - * \brief Add all buffers to data::Solution. + * \brief Move all buffers to data::Solution. */ void assignToSolution(Opm::data::Solution& sol) { @@ -517,101 +517,101 @@ public: return; if ( oilPressure_.size() > 0 ) { - sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, oilPressure_, Opm::data::TargetType::RESTART_SOLUTION); + sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, std::move(oilPressure_), Opm::data::TargetType::RESTART_SOLUTION); } if ( temperature_.size() > 0 ) { - sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, temperature_, Opm::data::TargetType::RESTART_SOLUTION); + sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, std::move(temperature_), Opm::data::TargetType::RESTART_SOLUTION); } if (saturationsOutput_()) { if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size()>0 ) { - sol.insert( "SWAT", Opm::UnitSystem::measure::identity, saturation_[waterPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); + 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, saturation_[gasPhaseIdx], Opm::data::TargetType::RESTART_SOLUTION ); + sol.insert( "SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION ); } } if ( gasDissolutionFactor_.size() > 0 ) { - sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, gasDissolutionFactor_, Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, std::move(gasDissolutionFactor_), Opm::data::TargetType::RESTART_AUXILIARY ); } if ( oilVaporizationFactor_.size() > 0 ) { - sol.insert( "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, oilVaporizationFactor_ , Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, std::move(oilVaporizationFactor_) , Opm::data::TargetType::RESTART_AUXILIARY ); } if (rs_.size() > 0 ) { - sol.insert( "RS", Opm::UnitSystem::measure::gas_oil_ratio, rs_, Opm::data::TargetType::RESTART_SOLUTION ); + sol.insert( "RS", Opm::UnitSystem::measure::gas_oil_ratio, std::move(rs_), Opm::data::TargetType::RESTART_SOLUTION ); } if (rv_.size() > 0 ) { - sol.insert( "RV", Opm::UnitSystem::measure::oil_gas_ratio, rv_ , Opm::data::TargetType::RESTART_SOLUTION ); + sol.insert( "RV", Opm::UnitSystem::measure::oil_gas_ratio, std::move(rv_) , Opm::data::TargetType::RESTART_SOLUTION ); } if (invBOutput_()) { if( FluidSystem::phaseIsActive(waterPhaseIdx) && invB_[waterPhaseIdx].size() > 0 ) { - sol.insert( "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, invB_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, std::move(invB_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } if( FluidSystem::phaseIsActive(oilPhaseIdx) && invB_[oilPhaseIdx].size() > 0 ) { - sol.insert( "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, invB_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, std::move(invB_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } if( FluidSystem::phaseIsActive(gasPhaseIdx) && invB_[gasPhaseIdx].size() > 0 ) { - sol.insert( "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, invB_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, std::move(invB_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } } if (densityOutput_()) { if( FluidSystem::phaseIsActive(waterPhaseIdx) && density_[waterPhaseIdx].size() > 0 ) { - sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, density_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, std::move(density_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } if( FluidSystem::phaseIsActive(oilPhaseIdx) && density_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, density_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, std::move(density_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } if( FluidSystem::phaseIsActive(gasPhaseIdx) && density_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, density_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, std::move(density_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } } if (viscosityOutput_()) { if( FluidSystem::phaseIsActive(waterPhaseIdx) && viscosity_[waterPhaseIdx].size() > 0 ) { - sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); } if( FluidSystem::phaseIsActive(oilPhaseIdx) && viscosity_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); } if( FluidSystem::phaseIsActive(gasPhaseIdx) && viscosity_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, viscosity_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY ); + sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); } } if (relativePermeabilityOutput_()) { if( FluidSystem::phaseIsActive(waterPhaseIdx) && relativePermeability_[waterPhaseIdx].size() > 0) { - sol.insert( "WATKR", Opm::UnitSystem::measure::identity, relativePermeability_[waterPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert( "WATKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); } if( FluidSystem::phaseIsActive(oilPhaseIdx)&& relativePermeability_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OILKR", Opm::UnitSystem::measure::identity, relativePermeability_[oilPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert( "OILKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); } if( FluidSystem::phaseIsActive(gasPhaseIdx)&& relativePermeability_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GASKR", Opm::UnitSystem::measure::identity, relativePermeability_[gasPhaseIdx], Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert( "GASKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); } } if (hysteresisOutput_()) { - sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, pcSwMdcOw_, Opm::data::TargetType::RESTART_AUXILIARY); - sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, krnSwMdcOw_, Opm::data::TargetType::RESTART_AUXILIARY); - sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, pcSwMdcGo_, Opm::data::TargetType::RESTART_AUXILIARY); - sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, krnSwMdcGo_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, std::move(pcSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, std::move(krnSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, std::move(pcSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, std::move(krnSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); } if (soMaxOutput_()) - sol.insert ("SOMAX", Opm::UnitSystem::measure::identity, soMax_, Opm::data::TargetType::RESTART_SOLUTION); + sol.insert ("SOMAX", Opm::UnitSystem::measure::identity, std::move(soMax_), Opm::data::TargetType::RESTART_SOLUTION); if (solventOutput_()) - sol.insert ("SSOL", Opm::UnitSystem::measure::identity, sSol_, Opm::data::TargetType::RESTART_SOLUTION); + sol.insert ("SSOL", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION); if (polymerOutput_()) - sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, cPolymer_, Opm::data::TargetType::RESTART_SOLUTION); + sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION); if (dewPointPressureOutput_() && dewPointPressure_.size() > 0) - sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, dewPointPressure_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, std::move(dewPointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); if (bubbelPointPressureOutput_() && bubblePointPressure_.size() > 0) - sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, bubblePointPressure_, Opm::data::TargetType::RESTART_AUXILIARY); + sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); } From abd3271f1ca60acf848a086996f044e358956d7d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 9 Jan 2018 13:08:05 +0100 Subject: [PATCH 3/4] Only collect to globalData for mpi --- ebos/eclwriter.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index b457f032c..1d5e15931 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -148,7 +148,8 @@ public: // collect all data to I/O rank and assign to sol Opm::data::Solution localCellData = fip; eclOutputModule_.assignToSolution(localCellData); - collectToIORank_.collect(localCellData); + if (collectToIORank_.isParallel()) + collectToIORank_.collect(localCellData); // write output on I/O rank if (collectToIORank_.isIORank()) { From f0ee6df1363b5251583024f3af8237f8dccba272 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 9 Jan 2018 14:01:30 +0100 Subject: [PATCH 4/4] Only compute auxiliary values if asked for Ask the restartConfig if restart files are written. Solution variables are always written since they may be needed by the summary files Remove output of saturated data to align with legacy code. --- ebos/ecloutputblackoilmodule.hh | 549 +++++++++++++------------------- ebos/eclwriter.hh | 4 +- 2 files changed, 223 insertions(+), 330 deletions(-) diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 222cd9d3b..ad96ece9b 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -90,66 +90,83 @@ public: * \brief Allocate memory for the scalar fields we would like to * write to ECL output files */ - void allocBuffers(unsigned bufferSize, unsigned reportStepNum, const Opm::RestartConfig& restartConfig, const bool log) + void allocBuffers(unsigned bufferSize, unsigned reportStepNum, const Opm::RestartConfig& restartConfig, const bool substep, const bool log) { + + if (!std::is_same >::value) + return; + std::map rstKeywords = restartConfig.getRestartKeywords(reportStepNum); for (auto& keyValue : rstKeywords) { keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum); } - if (!std::is_same >::value) + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + saturation_[phaseIdx].resize(bufferSize,0.0); + } + + 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; - if (saturationsOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) - saturation_[phaseIdx].resize(bufferSize,0.0); - } - if (pressuresOutput_()) { - oilPressure_.resize(bufferSize,0.0); - } - if (temperatureOutput_()) { - temperature_.resize(bufferSize,0.0); - } - - if (gasDissolutionFactorOutput_() && rstKeywords["RSSAT"] > 0) { - rstKeywords["RSSAT"] = 0; - gasDissolutionFactor_.resize(bufferSize,0.0); - } - if (oilVaporizationFactorOutput_() && rstKeywords["RVSAT"] > 0) { - rstKeywords["RVSAT"] = 0; - oilVaporizationFactor_.resize(bufferSize,0.0); - } - - if (gasFormationVolumeFactorOutput_()) - gasFormationVolumeFactor_.resize(bufferSize,0.0); - if (saturatedOilFormationVolumeFactorOutput_()) - saturatedOilFormationVolumeFactor_.resize(bufferSize,0.0); - if (oilSaturationPressureOutput_()) - oilSaturationPressure_.resize(bufferSize,0.0); - - if (rsOutput_()) - rs_.resize(bufferSize,0.0); - if (rvOutput_()) - rv_.resize(bufferSize,0.0); - if (invBOutput_()) { - if (FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["BW"] > 0) - { - rstKeywords["BW"] = 0; - invB_[waterPhaseIdx].resize(bufferSize,0.0); + // 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 (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["RVSAT"] > 0) { + rstKeywords["RVSAT"] = 0; + oilVaporizationFactor_.resize(bufferSize,0.0); } } - if (densityOutput_() && rstKeywords["DEN"] > 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)) @@ -157,59 +174,42 @@ public: density_[phaseIdx].resize(bufferSize,0.0); } } - if (viscosityOutput_()) { - 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; + 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) && hasVWAT) + { + rstKeywords["VWAT"] = 0; + viscosity_[waterPhaseIdx].resize(bufferSize,0.0); } - if (relativePermeabilityOutput_()) { - 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 (FluidSystem::phaseIsActive(oilPhaseIdx) && hasVOIL > 0) + { + rstKeywords["VOIL"] = 0; + viscosity_[oilPhaseIdx].resize(bufferSize,0.0); } - if (solventOutput_()) { - sSol_.resize(bufferSize,0.0); - } - if (polymerOutput_()) { - cPolymer_.resize(bufferSize,0.0); + if (FluidSystem::phaseIsActive(gasPhaseIdx) && hasVGAS > 0) + { + rstKeywords["VGAS"] = 0; + viscosity_[gasPhaseIdx].resize(bufferSize,0.0); } - // TODO: Only needed if Vappars or hysteresis. - // Now: Output the same as legacy - 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); + 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; @@ -231,6 +231,12 @@ public: failedCellsPb_.clear(); failedCellsPd_.clear(); + + // Not supported in flow legacy + if (false) + saturatedOilFormationVolumeFactor_.resize(bufferSize,0.0); + if (false) + oilSaturationPressure_.resize(bufferSize,0.0); } /*! @@ -251,15 +257,14 @@ public: 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 (saturation_[phaseIdx].size() == 0) - 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 (oilPressure_.size() > 0) { oilPressure_[globalDofIdx] = Toolbox::value(fs.pressure(oilPhaseIdx)); Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]); @@ -308,44 +313,36 @@ public: Opm::Valgrind::CheckDefined(rv_[globalDofIdx]); } - if (invBOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (invB_[phaseIdx].size() == 0) - continue; + 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]); - } + invB_[phaseIdx][globalDofIdx] = Toolbox::value(fs.invB(phaseIdx)); + Opm::Valgrind::CheckDefined(invB_[phaseIdx][globalDofIdx]); } - if (densityOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (density_[phaseIdx].size() == 0) - continue; + 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]); - } + density_[phaseIdx][globalDofIdx] = Toolbox::value(fs.density(phaseIdx)); + Opm::Valgrind::CheckDefined(density_[phaseIdx][globalDofIdx]); } - if (viscosityOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (viscosity_[phaseIdx].size() == 0) - continue; + 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]); - } + viscosity_[phaseIdx][globalDofIdx] = Toolbox::value(fs.viscosity(phaseIdx)); + Opm::Valgrind::CheckDefined(viscosity_[phaseIdx][globalDofIdx]); } - if (relativePermeabilityOutput_()) { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { - if (relativePermeability_[phaseIdx].size() == 0) - continue; + 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]); - } + relativePermeability_[phaseIdx][globalDofIdx] = Toolbox::value(intQuants.relativePermeability(phaseIdx)); + Opm::Valgrind::CheckDefined(relativePermeability_[phaseIdx][globalDofIdx]); } if (sSol_.size() > 0) { @@ -380,13 +377,15 @@ public: if (soMax_.size() > 0) soMax_[globalDofIdx] = elemCtx.simulator().model().maxOilSaturation(globalDofIdx); - if (hysteresisOutput_()) { - const auto& matLawManager = elemCtx.simulator().problem().materialLawManager(); - if (matLawManager->enableHysteresis()) { + 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], @@ -405,8 +404,11 @@ public: const auto& fs_initial = elemCtx.simulator().problem().initialFluidState(globalDofIdx); // use initial rs and rv values - rv_[globalDofIdx] = fs_initial.Rv(); - rs_[globalDofIdx] = fs_initial.Rs(); + 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) @@ -449,7 +451,7 @@ public: std::vector displPb, displPd, recvLenPb, recvLenPd; const auto& comm = simulator_.gridView().comm(); - if ( isIORank() ) + if ( isIORank_() ) { displPb.resize(comm.size()+1, 0); displPd.resize(comm.size()+1, 0); @@ -463,7 +465,7 @@ public: std::partial_sum(recvLenPd.begin(), recvLenPd.end(), displPd.begin()+1); std::vector globalFailedCellsPb, globalFailedCellsPd; - if ( isIORank() ) + if ( isIORank_() ) { globalFailedCellsPb.resize(displPb.back()); globalFailedCellsPd.resize(displPd.back()); @@ -524,14 +526,13 @@ public: sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, std::move(temperature_), Opm::data::TargetType::RESTART_SOLUTION); } - if (saturationsOutput_()) { - 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 ); - } + 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 ); + } + if ( gasDissolutionFactor_.size() > 0 ) { sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, std::move(gasDissolutionFactor_), Opm::data::TargetType::RESTART_AUXILIARY ); @@ -539,78 +540,78 @@ public: if ( oilVaporizationFactor_.size() > 0 ) { sol.insert( "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, std::move(oilVaporizationFactor_) , Opm::data::TargetType::RESTART_AUXILIARY ); } - if (rs_.size() > 0 ) { + if ( rs_.size() > 0 ) { sol.insert( "RS", Opm::UnitSystem::measure::gas_oil_ratio, std::move(rs_), Opm::data::TargetType::RESTART_SOLUTION ); } if (rv_.size() > 0 ) { sol.insert( "RV", Opm::UnitSystem::measure::oil_gas_ratio, std::move(rv_) , Opm::data::TargetType::RESTART_SOLUTION ); } - if (invBOutput_()) { - if( FluidSystem::phaseIsActive(waterPhaseIdx) && 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( FluidSystem::phaseIsActive(oilPhaseIdx) && 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( FluidSystem::phaseIsActive(gasPhaseIdx) && 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 (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 (densityOutput_()) { - if( FluidSystem::phaseIsActive(waterPhaseIdx) && density_[waterPhaseIdx].size() > 0 ) { - sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, std::move(density_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); - } - if( FluidSystem::phaseIsActive(oilPhaseIdx) && density_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, std::move(density_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); - } - if( FluidSystem::phaseIsActive(gasPhaseIdx) && density_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, std::move(density_[gasPhaseIdx]), 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 (viscosityOutput_()) { - if( FluidSystem::phaseIsActive(waterPhaseIdx) && viscosity_[waterPhaseIdx].size() > 0 ) { - sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); - } - if( FluidSystem::phaseIsActive(oilPhaseIdx) && viscosity_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); - } - if( FluidSystem::phaseIsActive(gasPhaseIdx) && viscosity_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY ); - } - } - if (relativePermeabilityOutput_()) { - if( FluidSystem::phaseIsActive(waterPhaseIdx) && relativePermeability_[waterPhaseIdx].size() > 0) { - sol.insert( "WATKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); - } - if( FluidSystem::phaseIsActive(oilPhaseIdx)&& relativePermeability_[oilPhaseIdx].size() > 0 ) { - sol.insert( "OILKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY); - } - if( FluidSystem::phaseIsActive(gasPhaseIdx)&& relativePermeability_[gasPhaseIdx].size() > 0 ) { - sol.insert( "GASKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[gasPhaseIdx]), 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 (hysteresisOutput_()) { + 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); - sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, std::move(krnSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY); - sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, std::move(pcSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); - sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, std::move(krnSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY); - } - if (soMaxOutput_()) + 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 (solventOutput_()) + if (sSol_.size() > 0) sol.insert ("SSOL", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION); - if (polymerOutput_()) + if (cPolymer_.size() > 0) sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION); - if (dewPointPressureOutput_() && dewPointPressure_.size() > 0) + if (dewPointPressure_.size() > 0) sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, std::move(dewPointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); - if (bubbelPointPressureOutput_() && bubblePointPressure_.size() > 0) + if (bubblePointPressure_.size() > 0) sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY); } @@ -681,15 +682,14 @@ public: void assignToFluidState(FluidState& fs, unsigned elemIdx) const { - 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; - fs.setSaturation(phaseIdx, saturation_[phaseIdx][elemIdx]); - } + fs.setSaturation(phaseIdx, saturation_[phaseIdx][elemIdx]); } - if (pressuresOutput_()) { + + 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.) @@ -707,48 +707,52 @@ public: } } - if (temperatureOutput_()) { + if (temperature_.size() > 0) { fs.setTemperature( temperature_[elemIdx]); } - if (rsOutput_()) { + if (rs_.size() > 0) { fs.setRs(rs_[elemIdx]); } - if (rvOutput_()) { + if (rv_.size() > 0) { fs.setRv(rv_[elemIdx]); } } void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const { - if (soMaxOutput_()) + 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); + pcSwMdcOw_[elemIdx], + krnSwMdcOw_[elemIdx], + elemIdx); + } + if (pcSwMdcGo_.size() > 0 && krnSwMdcGo_.size() > 0) { matLawManager->setGasOilHysteresisParams( - pcSwMdcGo_[elemIdx], - krnSwMdcGo_[elemIdx], - elemIdx); + pcSwMdcGo_[elemIdx], + krnSwMdcGo_[elemIdx], + elemIdx); + } } } Scalar getSolventSaturation(unsigned elemIdx) const { - if(solventOutput_()) + if(sSol_.size() > 0) return sSol_[elemIdx]; return 0; } Scalar getPolymerConcentration(unsigned elemIdx) const { - if(polymerOutput_()) + if(cPolymer_.size() > 0) return cPolymer_[elemIdx]; return 0; @@ -756,119 +760,8 @@ public: private: - // This should be cleaned up. - // For now output the same as legacy - // to make the tests pass - static bool saturationsOutput_() - { return true; } - static bool pressuresOutput_() - { return true; } - - static bool temperatureOutput_() - { - return true; - } - - static bool solventOutput_() - { - return GET_PROP_VALUE(TypeTag, EnableSolvent); - } - - static bool polymerOutput_() - { - return GET_PROP_VALUE(TypeTag, EnablePolymer); - } - - static bool gasDissolutionFactorOutput_() - { - return true; - //FluidSystem::enableDissolvedGas(); - } - - static bool gasFormationVolumeFactorOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx); - } - - static bool oilVaporizationFactorOutput_() - { - return true; - //FluidSystem::enableVaporizedOil(); - } - - static bool saturatedOilFormationVolumeFactorOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx); - } - - static bool oilSaturationPressureOutput_() - { - return - FluidSystem::phaseIsActive(oilPhaseIdx) && - FluidSystem::phaseIsActive(gasPhaseIdx); - } - - static bool rsOutput_() - { - // Output the same as legacy - return true; //FluidSystem::enableDissolvedGas(); - - } - - static bool rvOutput_() - { - // Output the same as legacy - return true;//FluidSystem::enableVaporizedOil(); - } - - - static bool invBOutput_() - { - return true; - } - - static bool densityOutput_() - { - return true; - } - - static bool viscosityOutput_() - { - return true; - } - - - static bool relativePermeabilityOutput_() - { - return true; - } - - static bool soMaxOutput_() - { - return true; - } - - static bool hysteresisOutput_() - { - return true; - } - - static bool dewPointPressureOutput_() - { - return true; - } - - static bool bubbelPointPressureOutput_() - { - return true; - } - - bool isIORank() const + bool isIORank_() const { const auto& comm = simulator_.gridView().comm(); return comm.rank() == 0; diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 1d5e15931..96ce8354e 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -132,7 +132,7 @@ public: 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); + eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), substep, log); ElementContext elemCtx(simulator_); ElementIterator elemIt = gridView.template begin(); @@ -199,7 +199,7 @@ public: 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); + 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) {