From 70ece6d25a26055b61ee48763dea4292fe712a5d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 May 2021 12:13:25 +0200 Subject: [PATCH] eclwriter: split in typetag dependent and typetag-independent parts --- CMakeLists_files.cmake | 1 + ebos/eclgenericwriter.cc | 538 +++++++++++++++++++++++++++++++++++ ebos/eclgenericwriter.hh | 129 +++++++++ ebos/eclproblem.hh | 3 +- ebos/eclwriter.hh | 595 +++++---------------------------------- tests/test_ecl_output.cc | 2 +- 6 files changed, 742 insertions(+), 526 deletions(-) create mode 100644 ebos/eclgenericwriter.cc create mode 100644 ebos/eclgenericwriter.hh diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 06699f610..e08965fd9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -29,6 +29,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/eclgenericthresholdpressure.cc ebos/eclgenerictracermodel.cc ebos/eclgenericvanguard.cc + ebos/eclgenericwriter.cc ebos/ecltransmissibility.cc opm/core/props/phaseUsageFromDeck.cpp opm/core/props/satfunc/RelpermDiagnostics.cpp diff --git a/ebos/eclgenericwriter.cc b/ebos/eclgenericwriter.cc new file mode 100644 index 000000000..99c404635 --- /dev/null +++ b/ebos/eclgenericwriter.cc @@ -0,0 +1,538 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#if HAVE_DUNE_FEM +#include +#include +#include +#endif + +#if HAVE_MPI +#include +#endif + +namespace { + +/*! + * \brief Detect whether two cells are direct vertical neighbours. + * + * I.e. have the same i and j index and all cartesian cells between them + * along the vertical column are inactive. + * + * \tparam CM The type of the cartesian index mapper. + * \param cartMapper The mapper onto cartesian indices. + * \param cartesianToActive The mapping of cartesian indices to active indices. + * \param smallGlobalIndex The cartesian cell index of the cell with smaller index + * \param largeGlobalIndex The cartesian cell index of the cell with larger index + * \return True if the cells have the same i and j indices and all cartesian cells + * between them are inactive. + */ +bool directVerticalNeighbors(const std::array& cartDims, + const std::unordered_map& cartesianToActive, + int smallGlobalIndex, int largeGlobalIndex) +{ + assert(smallGlobalIndex <= largeGlobalIndex); + std::array ijk1, ijk2; + auto globalToIjk = [cartDims](int gc) { + std::array ijk; + ijk[0] = gc % cartDims[0]; + gc /= cartDims[0]; + ijk[1] = gc % cartDims[1]; + ijk[2] = gc / cartDims[1]; + return ijk; + }; + ijk1 = globalToIjk(smallGlobalIndex); + ijk2 = globalToIjk(largeGlobalIndex); + assert(ijk2[2]>=ijk1[2]); + + if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) + { + assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); + for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; + gi += cartDims[0] * cartDims[1] ) + { + if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) + { + return false; + } + } + return true; + } else + return false; +} + + +struct EclWriteTasklet : public Opm::TaskletInterface +{ + Opm::Action::State actionState_; + Opm::SummaryState summaryState_; + Opm::UDQState udqState_; + Opm::EclipseIO& eclIO_; + int reportStepNum_; + bool isSubStep_; + double secondsElapsed_; + Opm::RestartValue restartValue_; + bool writeDoublePrecision_; + + explicit EclWriteTasklet(const Opm::Action::State& actionState, + const Opm::SummaryState& summaryState, + const Opm::UDQState& udqState, + Opm::EclipseIO& eclIO, + int reportStepNum, + bool isSubStep, + double secondsElapsed, + Opm::RestartValue restartValue, + bool writeDoublePrecision) + : actionState_(actionState) + , summaryState_(summaryState) + , udqState_(udqState) + , eclIO_(eclIO) + , reportStepNum_(reportStepNum) + , isSubStep_(isSubStep) + , secondsElapsed_(secondsElapsed) + , restartValue_(restartValue) + , writeDoublePrecision_(writeDoublePrecision) + { } + + // callback to eclIO serial writeTimeStep method + void run() + { + eclIO_.writeTimeStep(actionState_, + summaryState_, + udqState_, + reportStepNum_, + isSubStep_, + secondsElapsed_, + restartValue_, + writeDoublePrecision_); + } +}; + +} + +namespace Opm { + +template +EclGenericWriter:: +EclGenericWriter(const Schedule& schedule, + const EclipseState& eclState, + const SummaryConfig& summaryConfig, + const Grid& grid, + const EquilGrid* equilGrid, + const GridView& gridView, + const Dune::CartesianIndexMapper& cartMapper, + const Dune::CartesianIndexMapper* equilCartMapper, + const TransmissibilityType& globalTrans, + bool enableAsyncOutput) + : collectToIORank_(grid, + equilGrid, + gridView, + cartMapper, + equilCartMapper) + , grid_(grid) + , gridView_(gridView) + , schedule_(schedule) + , eclState_(eclState) + , summaryConfig_(summaryConfig) + , globalTrans_(globalTrans) + , cartMapper_(cartMapper) + , equilCartMapper_(equilCartMapper) + , equilGrid_(equilGrid) +{ + if (collectToIORank_.isIORank()) { + eclIO_.reset(new EclipseIO(eclState_, + UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()), + schedule_, + summaryConfig_)); + + const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 ); + wbp_index_list_ = wbp_calculators.index_list(); + } + if (collectToIORank_.isParallel()) { + const auto& comm = grid_.comm(); + unsigned long size = wbp_index_list_.size(); + comm.broadcast(&size, 1, collectToIORank_.ioRank); + if (!collectToIORank_.isIORank()) + wbp_index_list_.resize( size ); + comm.broadcast(wbp_index_list_.data(), size, collectToIORank_.ioRank); + } + // create output thread if enabled and rank is I/O rank + // async output is enabled by default if pthread are enabled + int numWorkerThreads = 0; + if (enableAsyncOutput && collectToIORank_.isIORank()) + numWorkerThreads = 1; + taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); +} + +template +const EclipseIO& EclGenericWriter:: +eclIO() const +{ + assert(eclIO_); + return *eclIO_; +} + +template +void EclGenericWriter:: +writeInit() +{ + if (collectToIORank_.isIORank()) { + std::map > integerVectors; + if (collectToIORank_.isParallel()) + integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); + auto cartMap = cartesianToCompressed(equilGrid_->size(0), + UgGridHelpers::globalCell(*equilGrid_)); + eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap)); + } +} + +template +data::Solution EclGenericWriter:: +computeTrans_(const std::unordered_map& cartesianToActive) const +{ + const auto& cartMapper = *equilCartMapper_; + const auto& cartDims = cartMapper.cartesianDimensions(); + const int globalSize = cartDims[0]*cartDims[1]*cartDims[2]; + + data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; + data::CellData trany = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; + data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; + + for (size_t i = 0; i < tranx.data.size(); ++i) { + tranx.data[0] = 0.0; + trany.data[0] = 0.0; + tranz.data[0] = 0.0; + } + + using GlobalGridView = typename EquilGrid::LeafGridView; + const GlobalGridView& globalGridView = equilGrid_->leafGridView(); + using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + // Ordering of compressed and uncompressed index should be the same + const int cartIdx1 = cartMapper.cartesianIndex( c1 ); + const int cartIdx2 = cartMapper.cartesianIndex( c2 ); + // Ordering of compressed and uncompressed index should be the same + assert(cartIdx1 <= cartIdx2); + int gc1 = std::min(cartIdx1, cartIdx2); + int gc2 = std::max(cartIdx1, cartIdx2); + + if (gc2 - gc1 == 1 && cartDims[0] > 1 ) { + tranx.data[gc1] = globalTrans_.transmissibility(c1, c2); + continue; // skip other if clauses as they are false, last one needs some computation + } + + if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { + trany.data[gc1] = globalTrans_.transmissibility(c1, c2); + continue; // skipt next if clause as it needs some computation + } + + if ( gc2 - gc1 == cartDims[0]*cartDims[1] || + directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) + tranz.data[gc1] = globalTrans_.transmissibility(c1, c2); + } + } + + return {{"TRANX", tranx}, + {"TRANY", trany}, + {"TRANZ", tranz}}; +} + +template +std::vector EclGenericWriter:: +exportNncStructure_(const std::unordered_map& cartesianToActive) const +{ + std::size_t nx = eclState_.getInputGrid().getNX(); + std::size_t ny = eclState_.getInputGrid().getNY(); + auto nncData = eclState_.getInputNNC().input(); + const auto& unitSystem = eclState_.getDeckUnitSystem(); + std::vector outputNnc; + std::size_t index = 0; + + for( const auto& entry : nncData ) { + // test whether NNC is not a neighboring connection + // cell2>=cell1 holds due to sortNncAndApplyEditnnc + assert( entry.cell2 >= entry.cell1 ); + auto cellDiff = entry.cell2 - entry.cell1; + + if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { + auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans); + // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 + if ( tt >= 1.0e-6 ) + outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans); + } + ++index; + } + + using GlobalGridView = typename EquilGrid::LeafGridView; + const GlobalGridView& globalGridView = equilGrid_->leafGridView(); + using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + + // Cartesian index mapper for the serial I/O grid + const auto& equilCartMapper = *equilCartMapper_; + const auto& cartDims = cartMapper_.cartesianDimensions(); + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1]; + std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2]; + + if ( cc2 < cc1 ) + std::swap(cc1, cc2); + + auto cellDiff = cc2 - cc1; + + if (cellDiff != 1 && + cellDiff != nx && + cellDiff != nx*ny && + !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { + // We need to check whether an NNC for this face was also specified + // via the NNC keyword in the deck (i.e. in the first origNncSize entries. + auto t = globalTrans_.transmissibility(c1, c2); + auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0)); + + while ( candidate != nncData.end() && candidate->cell1 == cc1 + && candidate->cell2 == cc2) { + t -= candidate->trans; + ++candidate; + } + // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC + // with corresponding EDITNNC above). In addition we do set small transmissibilties + // to zero when setting up the simulator. These will be ignored here, too. + auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t)); + if ( tt > 1e-12 ) + outputNnc.push_back({cc1, cc2, t}); + } + } + } + return outputNnc; +} + +template +void EclGenericWriter:: +doWriteOutput(const int reportStepNum, + const bool isSubStep, + data::Solution&& localCellData, + data::Wells&& localWellData, + data::GroupAndNetworkValues&& localGroupAndNetworkData, + const Action::State& actionState, + const UDQState& udqState, + const SummaryState& summaryState, + const std::vector& thresholdPressure, + Scalar curTime, + Scalar nextStepSize, + bool doublePrecision) +{ + const auto isParallel = this->collectToIORank_.isParallel(); + + RestartValue restartValue { + isParallel ? this->collectToIORank_.globalCellData() + : std::move(localCellData), + + isParallel ? this->collectToIORank_.globalWellData() + : std::move(localWellData), + + isParallel ? this->collectToIORank_.globalGroupAndNetworkData() + : std::move(localGroupAndNetworkData) + }; + + if (eclState_.getSimulationConfig().useThresholdPressure()) { + restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure, + thresholdPressure); + } + + // Add suggested next timestep to extra data. + if (! isSubStep) { + restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); + } + + // first, create a tasklet to write the data for the current time + // step to disk + auto eclWriteTasklet = std::make_shared( + actionState, summaryState, udqState, *this->eclIO_, + reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision); + + // then, make sure that the previous I/O request has been completed + // and the number of incomplete tasklets does not increase between + // time steps + this->taskletRunner_->barrier(); + + // finally, start a new output writing job + this->taskletRunner_->dispatch(std::move(eclWriteTasklet)); +} + +template +void EclGenericWriter:: +evalSummary(int reportStepNum, + Scalar curTime, + const std::map& wbpData, + const data::Wells& localWellData, + const data::GroupAndNetworkValues& localGroupAndNetworkData, + const std::map& localAquiferData, + const std::map, double>& blockData, + const std::map& miscSummaryData, + const std::map>& regionData, + SummaryState& summaryState, + UDQState& udqState, + const Inplace& inplace, + const Inplace& initialInPlace) +{ + std::vector buffer; + if (collectToIORank_.isIORank()) { + const auto& summary = eclIO_->summary(); + auto wbp_calculators = summary.wbp_calculators(reportStepNum); + + for (const auto& [global_index, pressure] : wbpData) + wbp_calculators.add_pressure( global_index, pressure ); + + const auto& wellData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalWellData() + : localWellData; + + const auto& groupAndNetworkData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalGroupAndNetworkData() + : localGroupAndNetworkData; + + const auto& aquiferData = this->collectToIORank_.isParallel() + ? this->collectToIORank_.globalAquiferData() + : localAquiferData; + + summary.eval(summaryState, + reportStepNum, + curTime, + wellData, + groupAndNetworkData, + miscSummaryData, + initialInPlace, + inplace, + wbp_calculators, + regionData, + blockData, + aquiferData); + + /* + Off-by-one-fun: The reportStepNum argument corresponds to the + report step these results will be written to, whereas the argument + to UDQ function evaluation corresponds to the report step we are + currently on. + */ + auto udq_step = reportStepNum - 1; + const auto& udq_config = schedule_.getUDQConfig(udq_step); + udq_config.eval( udq_step, schedule_.wellMatcher(udq_step), summaryState, udqState); + + buffer = summaryState.serialize(); + } + if (collectToIORank_.isParallel()) { +#ifdef HAVE_MPI + unsigned long buffer_size = buffer.size(); + MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD); + if (!collectToIORank_.isIORank()) + buffer.resize( buffer_size ); + + MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD); + if (!collectToIORank_.isIORank()) { + SummaryState& st = summaryState; + st.deserialize(buffer); + } +#endif + } +} + +#if HAVE_DUNE_FEM +template class EclGenericWriter>>, + Dune::MultipleCodimMultipleGeomTypeMapper>>>, + double>; +#else +template class EclGenericWriter>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double>; +#endif + +template class EclGenericWriter, + Dune::PolyhedralGrid<3,3,double>, + Dune::GridView>, + Dune::MultipleCodimMultipleGeomTypeMapper>, Dune::Impl::MCMGFailLayout>, + double>; + +} // namespace Opm diff --git a/ebos/eclgenericwriter.hh b/ebos/eclgenericwriter.hh new file mode 100644 index 000000000..3b6739ece --- /dev/null +++ b/ebos/eclgenericwriter.hh @@ -0,0 +1,129 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::EclWriter + */ +#ifndef EWOMS_ECL_GENERIC_WRITER_HH +#define EWOMS_ECL_GENERIC_WRITER_HH + +#include "collecttoiorank.hh" +#include + +#include + +#include +#include +#include +#include +#include + +namespace Opm { + +namespace Action { class State; } +class EclipseIO; +class EclipseState; +class Inplace; +struct NNCdata; +class Schedule; +class SummaryConfig; +class SummaryState; +class UDQState; + +template +class EclGenericWriter +{ + using CollectDataToIORankType = CollectDataToIORank; + using TransmissibilityType = EclTransmissibility; + +public: + // The Simulator object should preferably have been const - the + // only reason that is not the case is due to the SummaryState + // object owned deep down by the vanguard. + EclGenericWriter(const Schedule& schedule, + const EclipseState& eclState, + const SummaryConfig& summaryConfig, + const Grid& grid, + const EquilGrid* equilGrid, + const GridView& gridView, + const Dune::CartesianIndexMapper& cartMapper, + const Dune::CartesianIndexMapper* equilCartMapper, + const TransmissibilityType& globalTrans, + bool enableAsyncOutput); + + const EclipseIO& eclIO() const; + + void writeInit(); + +protected: + void doWriteOutput(const int reportStepNum, + const bool isSubStep, + data::Solution&& localCellData, + data::Wells&& localWellData, + data::GroupAndNetworkValues&& localGroupAndNetworkData, + const Action::State& actionState, + const UDQState& udqState, + const SummaryState& summaryState, + const std::vector& thresholdPressure, + Scalar curTime, + Scalar nextStepSize, + bool doublePrecision); + + void evalSummary(int reportStepNum, + Scalar curTime, + const std::map& wbpData, + const data::Wells& localWellData, + const data::GroupAndNetworkValues& localGroupAndNetworkData, + const std::map& localAquiferData, + const std::map, double>& blockData, + const std::map& miscSummaryData, + const std::map>& regionData, + SummaryState& summaryState, + UDQState& udqState, + const Inplace& inplace, + const Inplace& initialInPlace); + + CollectDataToIORankType collectToIORank_; + const Grid& grid_; + const GridView& gridView_; + const Schedule& schedule_; + const EclipseState& eclState_; + const SummaryConfig& summaryConfig_; + std::unique_ptr eclIO_; + std::unique_ptr taskletRunner_; + Scalar restartTimeStepSize_; + const TransmissibilityType& globalTrans_; + const Dune::CartesianIndexMapper& cartMapper_; + const Dune::CartesianIndexMapper* equilCartMapper_; + const EquilGrid* equilGrid_; + std::vector wbp_index_list_; + +private: + data::Solution computeTrans_(const std::unordered_map& cartesianToActive) const; + std::vector exportNncStructure_(const std::unordered_map& cartesianToActive) const; +}; + +} // namespace Opm + +#endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 6ae74ac8d..4ad200953 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -95,6 +95,7 @@ #include #include #include +#include #include #include #include @@ -838,7 +839,7 @@ public: ExtboModule::initFromState(vanguard.eclState()); // create the ECL writer - eclWriter_.reset(new EclWriterType(simulator)); + eclWriter_.reset(new EclWriterType(simulator, *this)); enableDriftCompensation_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableDriftCompensation); diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index d43e9831e..67f94984d 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -31,39 +31,13 @@ #include "collecttoiorank.hh" #include "ecloutputblackoilmodule.hh" -#include - -#include - -#include -#include -#include - -#include - -#include -#include #include -#include -#include -#include -#include #include -#include -#include -#include -#include +#include -#include -#include #include -#include - -#ifdef HAVE_MPI -#include -#endif namespace Opm::Properties { @@ -84,60 +58,9 @@ struct EclOutputDoublePrecision { namespace Opm { -template -class EclWriter; - -template -class EclOutputBlackOilModule; - -/*! - * \brief Detect whether two cells are direct vertical neighbours. - * - * I.e. have the same i and j index and all cartesian cells between them - * along the vertical column are inactive. - * - * \tparam CM The type of the cartesian index mapper. - * \param cartMapper The mapper onto cartesian indices. - * \param cartesianToActive The mapping of cartesian indices to active indices. - * \param smallGlobalIndex The cartesian cell index of the cell with smaller index - * \param largeGlobalIndex The cartesian cell index of the cell with larger index - * \return True if the cells have the same i and j indices and all cartesian cells - * between them are inactive. - */ -inline -bool directVerticalNeighbors(const std::array& cartDims, - const std::unordered_map& cartesianToActive, - int smallGlobalIndex, int largeGlobalIndex) -{ - assert(smallGlobalIndex <= largeGlobalIndex); - std::array ijk1, ijk2; - auto globalToIjk = [cartDims](int gc) { - std::array ijk; - ijk[0] = gc % cartDims[0]; - gc /= cartDims[0]; - ijk[1] = gc % cartDims[1]; - ijk[2] = gc / cartDims[1]; - return ijk; - }; - ijk1 = globalToIjk(smallGlobalIndex); - ijk2 = globalToIjk(largeGlobalIndex); - assert(ijk2[2]>=ijk1[2]); - - if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) - { - assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); - for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; - gi += cartDims[0] * cartDims[1] ) - { - if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) - { - return false; - } - } - return true; - } else - return false; -} +namespace Action { class State; } +class EclipseIO; +class UDQState; /*! * \ingroup EclBlackOilSimulator @@ -155,7 +78,11 @@ bool directVerticalNeighbors(const std::array& cartDims, * centered finite volume discretization. */ template -class EclWriter +class EclWriter : public EclGenericWriter, + GetPropType, + GetPropType, + GetPropType, + GetPropType> { using Simulator = GetPropType; using Vanguard = GetPropType; @@ -166,16 +93,13 @@ class EclWriter using ElementContext = GetPropType; using FluidSystem = GetPropType; using Element = typename GridView::template Codim<0>::Entity; + using ElementMapper = GetPropType; using ElementIterator = typename GridView::template Codim<0>::Iterator; - - using CollectDataToIORankType = CollectDataToIORank; - - typedef std::vector ScalarBuffer; + using BaseType = EclGenericWriter; enum { enableEnergy = getPropValue() }; enum { enableSolvent = getPropValue() }; - public: static void registerParameters() { @@ -188,73 +112,32 @@ public: // The Simulator object should preferably have been const - the // only reason that is not the case is due to the SummaryState // object owned deep down by the vanguard. - EclWriter(Simulator& simulator) - : simulator_(simulator) - , collectToIORank_(simulator_.vanguard().grid(), - simulator_.vanguard().grid().comm().rank() == 0 ? - &simulator_.vanguard().equilGrid() : nullptr, - simulator_.vanguard().gridView(), - simulator_.vanguard().cartesianIndexMapper(), - simulator_.vanguard().grid().comm().rank() == 0 ? - &simulator_.vanguard().equilCartesianIndexMapper() : nullptr) - + template + EclWriter(Simulator& simulator, const Problem& problem) + : BaseType(simulator.vanguard().schedule(), + simulator.vanguard().eclState(), + simulator.vanguard().summaryConfig(), + simulator.vanguard().grid(), + simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilGrid() : nullptr, + simulator.vanguard().gridView(), + simulator.vanguard().cartesianIndexMapper(), + simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilCartesianIndexMapper() : nullptr, + simulator.vanguard().grid().comm().size() > 1 ? simulator.vanguard().globalTransmissibility() : problem.eclTransmissibilities(), + EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput)) + , simulator_(simulator) { - std::vector wbp_index_list; - if (collectToIORank_.isIORank()) { - const auto& schedule = simulator_.vanguard().schedule(); - eclIO_.reset(new EclipseIO(simulator_.vanguard().eclState(), - UgGridHelpers::createEclipseGrid(globalGrid(), simulator_.vanguard().eclState().getInputGrid()), - schedule, - simulator_.vanguard().summaryConfig())); - - const auto& wbp_calculators = eclIO_->summary().wbp_calculators( schedule.size() - 1 ); - wbp_index_list = wbp_calculators.index_list(); - } - if (collectToIORank_.isParallel()) { - const auto& comm = simulator_.vanguard().grid().comm(); - unsigned long size = wbp_index_list.size(); - comm.broadcast(&size, 1, collectToIORank_.ioRank); - if (!collectToIORank_.isIORank()) - wbp_index_list.resize( size ); - comm.broadcast(wbp_index_list.data(), size, collectToIORank_.ioRank); - } - // create output thread if enabled and rank is I/O rank - // async output is enabled by default if pthread are enabled - bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput); - int numWorkerThreads = 0; - if (enableAsyncOutput && collectToIORank_.isIORank()) - numWorkerThreads = 1; - taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); - - this->eclOutputModule_ = std::make_unique>(simulator, wbp_index_list, this->collectToIORank_); + this->eclOutputModule_ = std::make_unique>(simulator, this->wbp_index_list_, this->collectToIORank_); + this->wbp_index_list_.clear(); } ~EclWriter() { } - const EclipseIO& eclIO() const - { - assert(eclIO_); - return *eclIO_; - } - const EquilGrid& globalGrid() const { return simulator_.vanguard().equilGrid(); } - void writeInit() - { - if (collectToIORank_.isIORank()) { - std::map > integerVectors; - if (collectToIORank_.isParallel()) - integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); - auto cartMap = cartesianToCompressed(globalGrid().size(0), - UgGridHelpers::globalCell(globalGrid())); - eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap)); - } - } - /*! * \brief collect and pass data and pass it to eclIO writer */ @@ -297,13 +180,13 @@ public: this->prepareLocalCellData(isSubStep, reportStepNum); - if (collectToIORank_.isParallel()) - collectToIORank_.collect({}, - eclOutputModule_->getBlockData(), - eclOutputModule_->getWBPData(), - localWellData, - localGroupAndNetworkData, - localAquiferData); + if (this->collectToIORank_.isParallel()) + this->collectToIORank_.collect({}, + eclOutputModule_->getBlockData(), + eclOutputModule_->getWBPData(), + localWellData, + localGroupAndNetworkData, + localAquiferData); std::map miscSummaryData; @@ -317,81 +200,25 @@ public: eclOutputModule_->outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput); eclOutputModule_->outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput); - - std::vector buffer; - if (this->collectToIORank_.isIORank()) { - const auto& summary = eclIO_->summary(); - auto wbp_calculators = summary.wbp_calculators(reportStepNum); - const auto& wbpData - = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalWBPData() - : this->eclOutputModule_->getWBPData(); - - for (const auto& [global_index, pressure] : wbpData) - wbp_calculators.add_pressure( global_index, pressure ); - - // Add TCPU - if (totalCpuTime != 0.0) { - miscSummaryData["TCPU"] = totalCpuTime; - } - - const auto& wellData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalWellData() - : localWellData; - - const auto& groupAndNetworkData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalGroupAndNetworkData() - : localGroupAndNetworkData; - - const auto& aquiferData = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalAquiferData() - : localAquiferData; - - const auto& blockData - = this->collectToIORank_.isParallel() - ? this->collectToIORank_.globalBlockData() - : this->eclOutputModule_->getBlockData(); - - summary.eval(summaryState(), - reportStepNum, - curTime, - wellData, - groupAndNetworkData, - miscSummaryData, - eclOutputModule_->initialInplace(), - inplace, - wbp_calculators, - regionData, - blockData, - aquiferData); - - /* - Off-by-one-fun: The reportStepNum argument corresponds to the - report step these results will be written to, whereas the argument - to UDQ function evaluation corresponds to the report step we are - currently on. - */ - auto udq_step = reportStepNum - 1; - const auto& udq_config = schedule().getUDQConfig(udq_step); - udq_config.eval( udq_step, schedule().wellMatcher(udq_step), summaryState(), udqState() ); - - buffer = summaryState().serialize(); + // Add TCPU + if (totalCpuTime != 0.0) { + miscSummaryData["TCPU"] = totalCpuTime; } - if (collectToIORank_.isParallel()) { -#ifdef HAVE_MPI - unsigned long buffer_size = buffer.size(); - MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD); - if (!collectToIORank_.isIORank()) - buffer.resize( buffer_size ); - - MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD); - if (!collectToIORank_.isIORank()) { - SummaryState& st = summaryState(); - st.deserialize(buffer); - } -#endif - } + this->evalSummary(reportStepNum, curTime, + this->collectToIORank_.isParallel() ? + this->collectToIORank_.globalWBPData() : + this->eclOutputModule_->getWBPData(), + localWellData, + localGroupAndNetworkData, + localAquiferData, + this->collectToIORank_.isParallel() ? + this->collectToIORank_.globalBlockData() : + this->eclOutputModule_->getBlockData(), + miscSummaryData, regionData, + summaryState(), udqState(), + inplace, + eclOutputModule_->initialInplace()); } @@ -416,19 +243,27 @@ public: } if (this->collectToIORank_.isParallel()) { - collectToIORank_.collect(localCellData, - eclOutputModule_->getBlockData(), - eclOutputModule_->getWBPData(), - localWellData, - localGroupAndNetworkData, - {}); + this->collectToIORank_.collect(localCellData, + eclOutputModule_->getBlockData(), + eclOutputModule_->getWBPData(), + localWellData, + localGroupAndNetworkData, + {}); } if (this->collectToIORank_.isIORank()) { - this->writeOutput(reportStepNum, isSubStep, - std::move(localCellData), - std::move(localWellData), - std::move(localGroupAndNetworkData)); + const Scalar curTime = simulator_.time() + simulator_.timeStepSize(); + const Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); + this->doWriteOutput(reportStepNum, isSubStep, + std::move(localCellData), + std::move(localWellData), + std::move(localGroupAndNetworkData), + this->actionState(), + this->udqState(), + this->summaryState(), + simulator_.problem().thresholdPressure().data(), + curTime, nextStepSize, + EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision)); } } @@ -470,10 +305,10 @@ public: { SummaryState& summaryState = simulator_.vanguard().summaryState(); Action::State& actionState = simulator_.vanguard().actionState(); - auto restartValues = loadParallelRestart(eclIO_.get(), actionState, summaryState, solutionKeys, extraKeys, + auto restartValues = loadParallelRestart(this->eclIO_.get(), actionState, summaryState, solutionKeys, extraKeys, gridView.grid().comm()); for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { - unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx); + unsigned globalIdx = this->collectToIORank_.localIdxToGlobalIdx(elemIdx); eclOutputModule_->setRestart(restartValues.solution, elemIdx, globalIdx); } @@ -502,247 +337,10 @@ public: Scalar restartTimeStepSize() const { return restartTimeStepSize_; } - private: static bool enableEclOutput_() { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); } - data::Solution computeTrans_(const std::unordered_map& cartesianToActive) const - { - const auto& cartMapper = simulator_.vanguard().equilCartesianIndexMapper(); - const auto& cartDims = cartMapper.cartesianDimensions(); - const int globalSize = cartDims[0]*cartDims[1]*cartDims[2]; - - data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; - data::CellData trany = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; - data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector(globalSize), data::TargetType::INIT}; - - for (size_t i = 0; i < tranx.data.size(); ++i) { - tranx.data[0] = 0.0; - trany.data[0] = 0.0; - tranz.data[0] = 0.0; - } - - typedef typename EquilGrid :: LeafGridView GlobalGridView; - const GlobalGridView& globalGridView = globalGrid().leafGridView(); - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - - using TransmissibilityType = typename Vanguard::TransmissibilityType; - - const TransmissibilityType* globalTrans; - - if (!collectToIORank_.isParallel()) - { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &simulator_.problem().eclTransmissibilities(); - } - else - { - globalTrans = &(simulator_.vanguard().globalTransmissibility()); - } - - auto elemIt = globalGridView.template begin(); - const auto& elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++ elemIt) { - const auto& elem = *elemIt; - - auto isIt = globalGridView.ibegin(elem); - const auto& isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - const auto& is = *isIt; - - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - // Ordering of compressed and uncompressed index should be the same - const int cartIdx1 = cartMapper.cartesianIndex( c1 ); - const int cartIdx2 = cartMapper.cartesianIndex( c2 ); - // Ordering of compressed and uncompressed index should be the same - assert(cartIdx1 <= cartIdx2); - int gc1 = std::min(cartIdx1, cartIdx2); - int gc2 = std::max(cartIdx1, cartIdx2); - - if (gc2 - gc1 == 1 && cartDims[0] > 1 ) { - tranx.data[gc1] = globalTrans->transmissibility(c1, c2); - continue; // skip other if clauses as they are false, last one needs some computation - } - - if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) { - trany.data[gc1] = globalTrans->transmissibility(c1, c2); - continue; // skipt next if clause as it needs some computation - } - - if ( gc2 - gc1 == cartDims[0]*cartDims[1] || - directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) - tranz.data[gc1] = globalTrans->transmissibility(c1, c2); - } - } - - return {{"TRANX", tranx}, - {"TRANY", trany}, - {"TRANZ", tranz}}; - } - - std::vector exportNncStructure_(const std::unordered_map& cartesianToActive) const - { - std::size_t nx = eclState().getInputGrid().getNX(); - std::size_t ny = eclState().getInputGrid().getNY(); - auto nncData = eclState().getInputNNC().input(); - const auto& unitSystem = simulator_.vanguard().eclState().getDeckUnitSystem(); - std::vector outputNnc; - std::size_t index = 0; - - for( const auto& entry : nncData ) { - // test whether NNC is not a neighboring connection - // cell2>=cell1 holds due to sortNncAndApplyEditnnc - assert( entry.cell2 >= entry.cell1 ); - auto cellDiff = entry.cell2 - entry.cell1; - - if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans); - // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 - if ( tt >= 1.0e-6 ) - outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans); - } - ++index; - } - - typedef typename EquilGrid :: LeafGridView GlobalGridView; - const GlobalGridView& globalGridView = globalGrid().leafGridView(); - typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); - - using TransmissibilityType = typename Vanguard::TransmissibilityType; - const TransmissibilityType* globalTrans; - if (!collectToIORank_.isParallel()) { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &simulator_.problem().eclTransmissibilities(); - } - else - { - globalTrans = &(simulator_.vanguard().globalTransmissibility()); - } - - // Cartesian index mapper for the serial I/O grid - const auto& equilCartMapper = simulator_.vanguard().equilCartesianIndexMapper(); - const auto& cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); - auto elemIt = globalGridView.template begin(); - const auto& elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++ elemIt) { - const auto& elem = *elemIt; - - auto isIt = globalGridView.ibegin(elem); - const auto& isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++ isIt) { - const auto& is = *isIt; - - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1]; - std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2]; - - if ( cc2 < cc1 ) - std::swap(cc1, cc2); - - auto cellDiff = cc2 - cc1; - - if (cellDiff != 1 && - cellDiff != nx && - cellDiff != nx*ny && - ! directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { - // We need to check whether an NNC for this face was also specified - // via the NNC keyword in the deck (i.e. in the first origNncSize entries. - auto t = globalTrans->transmissibility(c1, c2); - auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0)); - - while ( candidate != nncData.end() && candidate->cell1 == cc1 - && candidate->cell2 == cc2) { - t -= candidate->trans; - ++candidate; - } - // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC - // with corresponding EDITNNC above). In addition we do set small transmissibilties - // to zero when setting up the simulator. These will be ignored here, too. - auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t)); - if ( tt > 1e-12 ) - outputNnc.push_back({cc1, cc2, t}); - } - } - } - return outputNnc; - } - - struct EclWriteTasklet - : public TaskletInterface - { - Action::State actionState_; - SummaryState summaryState_; - UDQState udqState_; - EclipseIO& eclIO_; - int reportStepNum_; - bool isSubStep_; - double secondsElapsed_; - RestartValue restartValue_; - bool writeDoublePrecision_; - - explicit EclWriteTasklet(const Action::State& actionState, - const SummaryState& summaryState, - const UDQState& udqState, - EclipseIO& eclIO, - int reportStepNum, - bool isSubStep, - double secondsElapsed, - RestartValue restartValue, - bool writeDoublePrecision) - : actionState_(actionState) - , summaryState_(summaryState) - , udqState_(udqState) - , eclIO_(eclIO) - , reportStepNum_(reportStepNum) - , isSubStep_(isSubStep) - , secondsElapsed_(secondsElapsed) - , restartValue_(restartValue) - , writeDoublePrecision_(writeDoublePrecision) - { } - - // callback to eclIO serial writeTimeStep method - void run() - { - eclIO_.writeTimeStep(actionState_, - summaryState_, - udqState_, - reportStepNum_, - isSubStep_, - secondsElapsed_, - restartValue_, - writeDoublePrecision_); - } - }; - const EclipseState& eclState() const { return simulator_.vanguard().eclState(); } @@ -763,7 +361,7 @@ private: { const auto& gridView = simulator_.vanguard().gridView(); const int numElements = gridView.size(/*codim=*/0); - const bool log = collectToIORank_.isIORank(); + const bool log = this->collectToIORank_.isIORank(); eclOutputModule_->allocBuffers(numElements, reportStepNum, isSubStep, log, /*isRestart*/ false); @@ -782,59 +380,8 @@ private: } } - void writeOutput(const int reportStepNum, - const bool isSubStep, - ::Opm::data::Solution&& localCellData, - ::Opm::data::Wells&& localWellData, - ::Opm::data::GroupAndNetworkValues&& localGroupAndNetworkData) - { - const Scalar curTime = simulator_.time() + simulator_.timeStepSize(); - const Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); - const auto isParallel = this->collectToIORank_.isParallel(); - - RestartValue restartValue { - isParallel ? this->collectToIORank_.globalCellData() - : std::move(localCellData), - - isParallel ? this->collectToIORank_.globalWellData() - : std::move(localWellData), - - isParallel ? this->collectToIORank_.globalGroupAndNetworkData() - : std::move(localGroupAndNetworkData) - }; - - if (simulator_.vanguard().eclState().getSimulationConfig().useThresholdPressure()) { - restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure, - simulator_.problem().thresholdPressure().data()); - } - - // Add suggested next timestep to extra data. - if (! isSubStep) { - restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); - } - - // first, create a tasklet to write the data for the current time - // step to disk - auto eclWriteTasklet = std::make_shared( - this->actionState(), this->summaryState(), this->udqState(), *this->eclIO_, - reportStepNum, isSubStep, curTime, std::move(restartValue), - EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision) - ); - - // then, make sure that the previous I/O request has been completed - // and the number of incomplete tasklets does not increase between - // time steps - this->taskletRunner_->barrier(); - - // finally, start a new output writing job - this->taskletRunner_->dispatch(std::move(eclWriteTasklet)); - } - Simulator& simulator_; - CollectDataToIORankType collectToIORank_; std::unique_ptr> eclOutputModule_; - std::unique_ptr eclIO_; - std::unique_ptr taskletRunner_; Scalar restartTimeStepSize_; }; } // namespace Opm diff --git a/tests/test_ecl_output.cc b/tests/test_ecl_output.cc index 17a3cc474..3e9095be7 100644 --- a/tests/test_ecl_output.cc +++ b/tests/test_ecl_output.cc @@ -159,7 +159,7 @@ BOOST_AUTO_TEST_CASE(Summary) typedef Opm::EclWriter EclWriterType; // create the actual ECL writer - std::unique_ptr eclWriter = std::unique_ptr(new EclWriterType(*simulator)); + std::unique_ptr eclWriter = std::unique_ptr(new EclWriterType(*simulator, simulator->problem())); simulator->model().applyInitialSolution(); Opm::data::Wells dw;