EclGenericWriter.hpp

Go to the documentation of this file.

// -*- 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.
*/
#ifndef OPM_ECL_GENERIC_WRITER_HPP
#define OPM_ECL_GENERIC_WRITER_HPP
 
#include <opm/models/parallel/tasklets.hpp>
 
#include <opm/output/data/Groups.hpp>
 
#include <opm/simulators/flow/CollectDataOnIORank.hpp>
#include <opm/simulators/flow/LgrOutputTransGather.hpp>
#include <opm/simulators/flow/Transmissibility.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
 
#include <array>
#include <map>
#include <memory>
#include <optional>
#include <string>
#include <utility>
#include <vector>
 
namespace Opm {
 
class EclipseIO;
class EclipseState;
class InterRegFlowMap;
class Inplace;
template <class Grid> class LevelCartesianIndexMapper;
struct NNCdata;
class RegionVariableCollection;
class Schedule;
class SummaryConfig;
class SummaryState;
class UDQState;
 
} // namespace Opm
 
namespace Opm::Action {
class State;
} // namespace Opm::Action
 
namespace Opm::data {
class RegionVariableMapping;
} // namespace Opm::data
 
namespace Opm {
 
template <class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
class EclGenericWriter
{
    using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
    using EquilCartesianIndexMapper = Dune::CartesianIndexMapper<EquilGrid>;
    using CollectDataOnIORankType = CollectDataOnIORank<Grid,EquilGrid,GridView>;
    using TransmissibilityType = Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>;
 
public:
    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,
                     bool enableAsyncOutput,
                     bool enableEsmry);
 
    const EclipseIO& eclIO() const;
 
    void writeInit();
 
    void setTransmissibilities(const TransmissibilityType* globalTrans)
    {
        globalTrans_ = globalTrans;
    }
 
    // Parallel LGR INIT output: connection transmissibilities gathered from the
    // ranks' own (distributed) simulator transmissibilities. When set,
    // computeTrans_ / exportNncStructure_ take values from these instead of
    // querying a whole-grid transmissibility object. Complete on the I/O rank only.
    void setGatheredLgrTrans(std::optional<GatheredLgrOutputTrans> gatheredTrans)
    {
        gatheredLgrTrans_ = std::move(gatheredTrans);
    }
 
    void setSubStepReport(const SimulatorReportSingle& report)
    {
        sub_step_report_ = report;
    }
    void setSimulationReport(const SimulatorReport& report)
    {
        simulation_report_ = report;
    }
 
    const std::vector<std::vector<NNCdata>>& getOutputNnc() const
    {
        return outputNnc_;
    }
 
    const CollectDataOnIORankType& collectOnIORank() const
    {
        return collectOnIORank_;
    }
 
    void extractOutputTransAndNNC(const std::function<unsigned int(unsigned int)>& map);
 
protected:
    const TransmissibilityType& globalTrans() const;
    unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const;
 
    // Value for `key` from the gathered records, or nullptr when the key is
    // absent or no gathered records are set. N = 3 selects the same-level
    // records, N = 4 the level-crossing (NNC) ones -- see GatheredLgrOutputTrans.
    template <std::size_t N>
    const double* findGatheredTrans_(const std::array<int,N>& key) const;
 
    // Output transmissibility value for a connection: from the gathered per-rank
    // simulator transmissibilities when set (parallel LGR runs; a missing key is
    // a hard error), otherwise from the whole-grid transmissibility object
    // (c1, c2). Key shape as in findGatheredTrans_.
    template <std::size_t N>
    double gatheredOrGlobalTrans_(const std::array<int,N>& key,
                                  unsigned c1,
                                  unsigned c2) const;
 
    void doWriteOutput(const int                          reportStepNum,
                       const std::optional<int>           timeStepNum,
                       const bool                         isSubStep,
                       const bool                         forcedSimulationFinished,
                       data::Solution&&                   localCellData,
                       data::Wells&&                      localWellData,
                       data::GroupAndNetworkValues&&      localGroupAndNetworkData,
                       data::Aquifers&&                   localAquiferData,
                       WellTestState&&                    localWTestState,
                       const Action::State&               actionState,
                       const UDQState&                    udqState,
                       const SummaryState&                summaryState,
                       const std::vector<Scalar>&         thresholdPressure,
                       Scalar                             curTime,
                       Scalar                             nextStepSize,
                       bool                               doublePrecision,
                       bool                               isFlowsn,
                       std::array<FlowsData<double>,3>&&  flowsn,
                       bool                               isFloresn,
                       std::array<FlowsData<double>, 3>&& floresn);
 
    void evalSummary(int                                                  reportStepNum,
                     Scalar                                               curTime,
                     const data::Wells&                                   localWellData,
                     const data::WellBlockAveragePressures&               localWBPData,
                     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::tuple<std::string, int, int>, double>& lgrBlockData,
                     const std::map<std::string, double>&                 miscSummaryData,
                     const std::map<std::string, std::vector<double>>&    regionData,
                     const data::RegionVariableMapping&                   regVarMap,
                     const RegionVariableCollection&                      regVars,
                     const Inplace&                                       inplace,
                     const Inplace*                                       initialInPlace,
                     const InterRegFlowMap&                               interRegFlows,
                     SummaryState&                                        summaryState,
                     UDQState&                                            udqState,
                     const data::ReservoirCouplingGroupRates*             rcGroupRates = nullptr);
 
    CollectDataOnIORankType collectOnIORank_;
    const Grid& grid_;
    const GridView& gridView_;
    const Schedule& schedule_;
    const EclipseState& eclState_;
    std::unique_ptr<EclipseIO> eclIO_;
    std::unique_ptr<TaskletRunner> taskletRunner_;
    Scalar restartTimeStepSize_;
    const TransmissibilityType* globalTrans_ = nullptr;
    std::optional<GatheredLgrOutputTrans> gatheredLgrTrans_;
    const Dune::CartesianIndexMapper<Grid>& cartMapper_;
    const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper_;
    const EquilGrid* equilGrid_;
    SimulatorReportSingle sub_step_report_;
    SimulatorReport simulation_report_;
 
    // Regular NNCs per grid: internal to a grid.
    // Both cells belong to the same level grid, either the main grid or a level/local grid.
    // nnc.cell1 (NNC1 in *.EGRID) Level/Local Cartesian index of cell1
    // nnc.cell2 (NNC2 in *.EGRID) Level/Local Cartesian index of cell2
    // Equivalent to TRANNNC in *.INIT
    mutable std::vector<std::vector<NNCdata>> outputNnc_;
 
    // NNCs between main (level zero) grid and local grids:
    // nnc.cell1 (NNCG in *.EGRID) Cartesian index of cell1 in the main grid where the cell belongs to.
    // nnc.cell2 (NNCL in *.EGRID) Level/Local Cartesian index of cell2 in the refined level grid where the cell belongs to.
    // Equivalent to TRANGL in *.INIT
    mutable std::vector<std::vector<NNCdata>> outputNncGlobalLocal_; // here GLOBAL refers to level 0 grid, local to any LGR (refined grid)
 
    // Amalgamated NNCs: nncs between different LGRs. For example, nested refinement or neighboring LGRs.
    // The cells belong to different refined level grids.
    // nnc.cell1 (NNA1 in *.EGRID) Level/Local Cartesian index of cell1 (in its level grid: level1)
    // nnc.cell2 (NNA2 in *.EGRID) Level/Local Cartesian index of cell2 (in its level grid: level2, with level2 > level1).
    // Equivalent to TRANLL in *.INIT
    mutable std::vector<std::vector<std::vector<NNCdata>>> outputAmalgamatedNnc_;
 
    mutable std::unique_ptr<std::vector<data::Solution>> outputTrans_;
 
private:
    template<typename LevelIndicesFunction, typename OriginIndicesFunction>
    void computeTrans_(const std::vector<std::unordered_map<int,int>>& levelCartToLevelCompressed,
                       const std::function<unsigned int(unsigned int)>& map,
                       const LevelIndicesFunction& computeLevelIndices,
                       const std::function<int(int, int)>& computeLevelCartIdx,
                       const std::function<std::array<int,3>(int)>& computeLevelCartDimensions,
                       const OriginIndicesFunction& computeOriginIndices) const;
 
    template<typename LevelIndicesFunction, typename OriginIndicesFunction>
    std::vector<std::vector<NNCdata>> exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToLevelCompressed,
                                                          const std::function<unsigned int(unsigned int)>& map,
                                                          const LevelIndicesFunction& computeLevelIndices,
                                                          const std::function<int(int, int)>& computeLevelCartIdx,
                                                          const std::function<std::array<int,3>(int)>& computeLevelCartDimensions,
                                                          const OriginIndicesFunction& computeOriginIndices) const;
 
    bool isCartesianNeighbour_(const std::array<int,3>& levelCartDims,
                               const std::size_t levelCartIdx1,
                               const std::size_t levelCartIdx2) const;
 
    bool isDirectNeighbours_(const std::unordered_map<int,int>& levelCartesianToActive,
                             const std::array<int,3>& levelCartDims,
                             const std::size_t levelCartIdx1,
                             const std::size_t levelCartIdx2) const;
 
    auto activeCell_(const std::unordered_map<int,int>& levelCartToLevelCompressed,
                     const std::size_t levelCartIdx) const;
 
 
 
    bool isNumAquCell_(const std::size_t cartIdx) const;
    bool isNumAquConn_(const std::size_t cartIdx1, const std::size_t cartIdx2) const;
 
    template<bool equilGridIsCpGrid>
    Opm::LevelCartesianIndexMapper<EquilGrid> createLevelCartMapp_() const;
 
    template<bool equilGridIsCpGrid>
    std::vector<std::unordered_map<int,int>> createCartesianToActiveMaps_(const Opm::LevelCartesianIndexMapper<EquilGrid>& levelCartMapp) const;
 
    template<bool equilGridIsCpGrid>
    std::function<std::array<int,3>(int)> computeLevelCartDimensions_(const Opm::LevelCartesianIndexMapper<EquilGrid>& levelCartMapp,
                                                                      const Dune::CartesianIndexMapper<EquilGrid>& equilCartMapp) const;
 
    template<bool equilGridIsCpGrid>
    std::function<int(int, int)> computeLevelCartIdx_(const Opm::LevelCartesianIndexMapper<EquilGrid>& levelCartMapp,
                                                      const Dune::CartesianIndexMapper<EquilGrid>& equilCartMapp) const;
 
    template <bool equilGridIsCpGrid>
    auto computeLevelIndices_() const;
 
    template <bool equilGridIsCpGrid>
    auto computeOriginIndices_() const;
 
    void allocateLevelTrans_(const std::array<int,3>& levelCartDims,
                             data::Solution& levelTrans) const;
    void allocateAllNncs_(int maxLevel) const;
};
 
} // namespace Opm
 
#endif // OPM_ECL_GENERIC_WRITER_HPP