Merge pull request #3264 from akva2/eclwriter_typetag_split

eclwriter: split in typetag dependent and typetag-independent parts

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

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 <http://www.gnu.org/licenses/>.
+
+  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 <config.h>
+#include <ebos/eclgenericwriter.hh>
+
+#include <opm/grid/CpGrid.hpp>
+#include <opm/grid/cpgrid/GridHelpers.hpp>
+#include <opm/grid/polyhedralgrid.hh>
+#include <opm/grid/utility/cartesianToCompressed.hpp>
+
+#include <opm/output/eclipse/EclipseIO.hpp>
+#include <opm/output/eclipse/RestartValue.hpp>
+#include <opm/output/eclipse/Summary.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/SummaryState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQState.hpp>
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+#include <dune/grid/common/mcmgmapper.hh>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/gridpart/common/gridpart2gridview.hh>
+#include <ebos/femcpgridcompat.hh>
+#endif
+
+#if HAVE_MPI
+#include <mpi.h>
+#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<int, 3>& cartDims,
+                             const std::unordered_map<int,int>& cartesianToActive,
+                             int smallGlobalIndex, int largeGlobalIndex)
+{
+    assert(smallGlobalIndex <= largeGlobalIndex);
+    std::array<int, 3> ijk1, ijk2;
+    auto globalToIjk = [cartDims](int gc) {
+                           std::array<int, 3> 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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+EclGenericWriter(const Schedule& schedule,
+                 const EclipseState& eclState,
+                 const SummaryConfig& summaryConfig,
+                 const Grid& grid,
+                 const EquilGrid* equilGrid,
+                 const GridView& gridView,
+                 const Dune::CartesianIndexMapper<Grid>& cartMapper,
+                 const Dune::CartesianIndexMapper<EquilGrid>* 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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+eclIO() const
+{
+    assert(eclIO_);
+    return *eclIO_;
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+writeInit()
+{
+    if (collectToIORank_.isIORank()) {
+        std::map<std::string, std::vector<int> > 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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+computeTrans_(const std::unordered_map<int,int>& 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<double>(globalSize), data::TargetType::INIT};
+    data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
+    data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>(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<GlobalGridView>;
+    GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
+
+    auto elemIt = globalGridView.template begin</*codim=*/0>();
+    const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
+    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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+exportNncStructure_(const std::unordered_map<int,int>& 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<NNCdata> 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<GlobalGridView>;
+    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</*codim=*/0>();
+    const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
+    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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+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<Scalar>& 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<double>(1, nextStepSize));
+    }
+
+    // first, create a tasklet to write the data for the current time
+    // step to disk
+    auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
+        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<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+evalSummary(int reportStepNum,
+            Scalar curTime,
+            const std::map<std::size_t, double>& wbpData,
+            const data::Wells& localWellData,
+            const data::GroupAndNetworkValues& localGroupAndNetworkData,
+            const std::map<int,data::AquiferData>& localAquiferData,
+            const std::map<std::pair<std::string, int>, double>& blockData,
+            const std::map<std::string, double>& miscSummaryData,
+            const std::map<std::string, std::vector<double>>& regionData,
+            SummaryState& summaryState,
+            UDQState& udqState,
+            const Inplace& inplace,
+            const Inplace& initialInPlace)
+{
+    std::vector<char> 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::CpGrid,
+                                Dune::CpGrid,
+                                Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
+                                Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
+                                double>;
+#else
+template class EclGenericWriter<Dune::CpGrid,
+                                Dune::CpGrid,
+                                Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
+                                Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>, Dune::Impl::MCMGFailLayout>,
+                                double>;
+#endif
+
+template class EclGenericWriter<Dune::PolyhedralGrid<3,3,double>,
+                                Dune::PolyhedralGrid<3,3,double>,
+                                Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>,
+                                Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>, Dune::Impl::MCMGFailLayout>,
+                                double>;
+
+} // namespace Opm

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 <http://www.gnu.org/licenses/>.
+
+  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 <ebos/ecltransmissibility.hh>
+
+#include <opm/models/parallel/tasklets.hh>
+
+#include <cstddef>
+#include <map>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Opm {
+
+namespace Action { class State; }
+class EclipseIO;
+class EclipseState;
+class Inplace;
+struct NNCdata;
+class Schedule;
+class SummaryConfig;
+class SummaryState;
+class UDQState;
+
+template <class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+class EclGenericWriter
+{
+    using CollectDataToIORankType = CollectDataToIORank<Grid,EquilGrid,GridView>;
+    using TransmissibilityType = EclTransmissibility<Grid,GridView,ElementMapper,Scalar>;
+
+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<Grid>& cartMapper,
+                     const Dune::CartesianIndexMapper<EquilGrid>* 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<Scalar>& thresholdPressure,
+                       Scalar curTime,
+                       Scalar nextStepSize,
+                       bool doublePrecision);
+
+    void evalSummary(int reportStepNum,
+                     Scalar curTime,
+                     const std::map<std::size_t, double>& wbpData,
+                     const data::Wells& localWellData,
+                     const data::GroupAndNetworkValues& localGroupAndNetworkData,
+                     const std::map<int,data::AquiferData>& localAquiferData,
+                     const std::map<std::pair<std::string, int>, double>& blockData,
+                     const std::map<std::string, double>& miscSummaryData,
+                     const std::map<std::string, std::vector<double>>& 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<EclipseIO> eclIO_;
+    std::unique_ptr<TaskletRunner> taskletRunner_;
+    Scalar restartTimeStepSize_;
+    const TransmissibilityType& globalTrans_;
+    const Dune::CartesianIndexMapper<Grid>& cartMapper_;
+    const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper_;
+    const EquilGrid* equilGrid_;
+    std::vector<std::size_t> wbp_index_list_;
+
+private:
+    data::Solution computeTrans_(const std::unordered_map<int,int>& cartesianToActive) const;
+    std::vector<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive) const;
+};
+
+} // namespace Opm
+
+#endif

ebos/eclproblem.hh

@@ -95,6 +95,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionContext.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Action/ActionX.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
 #include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RockwnodTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/OverburdTable.hpp>
@@ -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);
 

ebos/eclwriter.hh

@@ -31,39 +31,13 @@
 #include "collecttoiorank.hh"
 #include "ecloutputblackoilmodule.hh"
 
-#include <opm/models/blackoil/blackoilmodel.hh>
-
-#include <opm/simulators/wells/BlackoilWellModel.hpp>
-
-#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
-#include <opm/models/io/baseoutputwriter.hh>
-#include <opm/models/parallel/tasklets.hh>
-
-#include <opm/output/eclipse/EclipseIO.hpp>
-
-#include <opm/output/eclipse/RestartValue.hpp>
-#include <opm/output/eclipse/Summary.hpp>
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Action/State.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQConfig.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQState.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Well/PAvgCalculatorCollection.hpp>
 
 #include <opm/simulators/utils/ParallelRestart.hpp>
-#include <opm/grid/GridHelpers.hpp>
-#include <opm/grid/utility/cartesianToCompressed.hpp>
 
-#include <opm/material/common/Valgrind.hpp>
+#include <ebos/eclgenericwriter.hh>
-#include <opm/common/OpmLog/OpmLog.hpp>
 
-#include <list>
-#include <utility>
 #include <string>
-#include <chrono>
-
-#ifdef HAVE_MPI
-#include <mpi.h>
-#endif
 
 namespace Opm::Properties {
 
@@ -84,60 +58,9 @@ struct EclOutputDoublePrecision {
 
 namespace Opm {
 
-template <class TypeTag>
+namespace Action { class State; }
-class EclWriter;
+class EclipseIO;
-
+class UDQState;
-template <class TypeTag>
-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<int, 3>& cartDims,
-                             const std::unordered_map<int,int>& cartesianToActive,
-                             int smallGlobalIndex, int largeGlobalIndex)
-{
-    assert(smallGlobalIndex <= largeGlobalIndex);
-    std::array<int, 3> ijk1, ijk2;
-    auto globalToIjk = [cartDims](int gc) {
-                           std::array<int, 3> 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;
-}
 
 /*!
  * \ingroup EclBlackOilSimulator
@@ -155,7 +78,11 @@ bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
  *   centered finite volume discretization.
  */
 template <class TypeTag>
-class EclWriter
+class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>,
+                                          GetPropType<TypeTag, Properties::EquilGrid>,
+                                          GetPropType<TypeTag, Properties::GridView>,
+                                          GetPropType<TypeTag, Properties::ElementMapper>,
+                                          GetPropType<TypeTag, Properties::Scalar>>
 {
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
     using Vanguard = GetPropType<TypeTag, Properties::Vanguard>;
@@ -166,16 +93,13 @@ class EclWriter
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Element = typename GridView::template Codim<0>::Entity;
+    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
     using ElementIterator = typename GridView::template Codim<0>::Iterator;
-
+    using BaseType = EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>;
-    using CollectDataToIORankType = CollectDataToIORank<Grid,EquilGrid,GridView>;
-
-    typedef std::vector<Scalar> ScalarBuffer;
 
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
 
-
 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)
+    template<class Problem>
-        : simulator_(simulator)
+    EclWriter(Simulator& simulator, const Problem& problem)
-        , collectToIORank_(simulator_.vanguard().grid(),
+        : BaseType(simulator.vanguard().schedule(),
-                           simulator_.vanguard().grid().comm().rank() == 0 ?
+                   simulator.vanguard().eclState(),
-                           &simulator_.vanguard().equilGrid() : nullptr,
+                   simulator.vanguard().summaryConfig(),
-                           simulator_.vanguard().gridView(),
+                   simulator.vanguard().grid(),
-                           simulator_.vanguard().cartesianIndexMapper(),
+                   simulator.vanguard().grid().comm().rank() == 0 ? &simulator.vanguard().equilGrid() : nullptr,
-                           simulator_.vanguard().grid().comm().rank() == 0 ?
+                   simulator.vanguard().gridView(),
-                           &simulator_.vanguard().equilCartesianIndexMapper() : nullptr)
+                   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<std::size_t> wbp_index_list;
+        this->eclOutputModule_ = std::make_unique<EclOutputBlackOilModule<TypeTag>>(simulator, this->wbp_index_list_, this->collectToIORank_);
-        if (collectToIORank_.isIORank()) {
+        this->wbp_index_list_.clear();
-            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<EclOutputBlackOilModule<TypeTag>>(simulator, wbp_index_list, this->collectToIORank_);
     }
 
     ~EclWriter()
     { }
 
-    const EclipseIO& eclIO() const
-    {
-        assert(eclIO_);
-        return *eclIO_;
-    }
-
     const EquilGrid& globalGrid() const
     {
         return simulator_.vanguard().equilGrid();
     }
 
-    void writeInit()
-    {
-        if (collectToIORank_.isIORank()) {
-            std::map<std::string, std::vector<int> > 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,8 +180,8 @@ public:
 
         this->prepareLocalCellData(isSubStep, reportStepNum);
 
-        if (collectToIORank_.isParallel())
+        if (this->collectToIORank_.isParallel())
-            collectToIORank_.collect({},
+            this->collectToIORank_.collect({},
                                            eclOutputModule_->getBlockData(),
                                            eclOutputModule_->getWBPData(),
                                            localWellData,
@@ -317,81 +200,25 @@ public:
         eclOutputModule_->outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput);
         eclOutputModule_->outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput);
 
-
-        std::vector<char> 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->evalSummary(reportStepNum, curTime,
-                ? this->collectToIORank_.globalWellData()
+                          this->collectToIORank_.isParallel() ?
-                : localWellData;
+                            this->collectToIORank_.globalWBPData() :
-
+                            this->eclOutputModule_->getWBPData(),
-            const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
+                          localWellData,
-                ? this->collectToIORank_.globalGroupAndNetworkData()
+                          localGroupAndNetworkData,
-                : localGroupAndNetworkData;
+                          localAquiferData,
-
+                          this->collectToIORank_.isParallel() ?
-            const auto& aquiferData = this->collectToIORank_.isParallel()
+                            this->collectToIORank_.globalBlockData() :
-                ? this->collectToIORank_.globalAquiferData()
+                            this->eclOutputModule_->getBlockData(),
-                : localAquiferData;
+                          miscSummaryData, regionData,
-
+                          summaryState(), udqState(),
-            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,
+                          eclOutputModule_->initialInplace());
-                         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
-        }
     }
 
 
@@ -416,7 +243,7 @@ public:
         }
 
         if (this->collectToIORank_.isParallel()) {
-            collectToIORank_.collect(localCellData,
+            this->collectToIORank_.collect(localCellData,
                                            eclOutputModule_->getBlockData(),
                                            eclOutputModule_->getWBPData(),
                                            localWellData,
@@ -425,10 +252,18 @@ public:
         }
 
         if (this->collectToIORank_.isIORank()) {
-            this->writeOutput(reportStepNum, isSubStep,
+            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));
+                                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<int,int>& 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<double>(globalSize), data::TargetType::INIT};
-        data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>(globalSize), data::TargetType::INIT};
-        data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>(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<GlobalGridView> 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</*codim=*/0>();
-        const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
-        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<NNCdata> exportNncStructure_(const std::unordered_map<int,int>& 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<NNCdata> 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<GlobalGridView> 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</*codim=*/0>();
-        const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
-        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<double>(1, nextStepSize));
-        }
-
-        // first, create a tasklet to write the data for the current time
-        // step to disk
-        auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
-            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<EclOutputBlackOilModule<TypeTag>> eclOutputModule_;
-    std::unique_ptr<EclipseIO> eclIO_;
-    std::unique_ptr<TaskletRunner> taskletRunner_;
     Scalar restartTimeStepSize_;
 };
 } // namespace Opm

tests/test_ecl_output.cc

@@ -159,7 +159,7 @@ BOOST_AUTO_TEST_CASE(Summary)
 
     typedef Opm::EclWriter<TypeTag> EclWriterType;
     // create the actual ECL writer
-    std::unique_ptr<EclWriterType> eclWriter = std::unique_ptr<EclWriterType>(new EclWriterType(*simulator));
+    std::unique_ptr<EclWriterType> eclWriter = std::unique_ptr<EclWriterType>(new EclWriterType(*simulator, simulator->problem()));
 
     simulator->model().applyInitialSolution();
     Opm::data::Wells dw;