GenericOutputModule.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_GENERIC_OUTPUT_MODULE_HPP
#define OPM_GENERIC_OUTPUT_MODULE_HPP
 
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
 
#include <opm/material/common/MathToolbox.hpp>
 
#include <opm/output/data/RegionVariableMapping.hpp>
#include <opm/output/data/Wells.hpp>
 
#include <opm/output/eclipse/Inplace.hpp>
#include <opm/output/eclipse/RegionVariableCollection.hpp>
 
#include <opm/simulators/flow/BioeffectsContainer.hpp>
#include <opm/simulators/flow/CO2H2Container.hpp>
#include <opm/simulators/flow/ExtboContainer.hpp>
#include <opm/simulators/flow/FIPContainer.hpp>
#include <opm/simulators/flow/FlowsContainer.hpp>
#include <opm/simulators/flow/InterRegFlows.hpp>
#include <opm/simulators/flow/LogOutputHelper.hpp>
#include <opm/simulators/flow/MechContainer.hpp>
#include <opm/simulators/flow/RegionPhasePVAverage.hpp>
#include <opm/simulators/flow/RFTContainer.hpp>
#include <opm/simulators/flow/RSTConv.hpp>
#include <opm/simulators/flow/GeochemistryContainer.hpp>
#include <opm/simulators/flow/TpsaContainer.hpp>
#include <opm/simulators/flow/TracerContainer.hpp>
 
#include <opm/simulators/utils/ParallelCommunication.hpp>
 
#include <array>
#include <cassert>
#include <cstddef>
#include <functional>
#include <map>
#include <memory>
#include <optional>
#include <string_view>
#include <tuple>
#include <unordered_map>
#include <utility>
#include <vector>
 
namespace Opm::Parameters {
 
struct ForceDisableFluidInPlaceOutput { static constexpr bool value = false; };
struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false; };
 
} // namespace Opm::Parameters
 
namespace Opm {
 
namespace data { class Solution; }
class EclHysteresisConfig;
class EclipseGrid;
class EclipseState;
class Schedule;
class SummaryConfig;
class SummaryConfigNode;
class SummaryState;
 
template<class FluidSystem>
class GenericOutputModule {
public:
    using Scalar = typename FluidSystem::Scalar;
 
    // Virtual destructor for safer inheritance.
    virtual ~GenericOutputModule();
 
    static void registerParameters();
 
    void outputTimeStamp(const std::string& lbl,
                         double elapsed,
                         int rstep,
                         boost::posix_time::ptime currentDate);
 
    void prepareDensityAccumulation();
 
    void accumulateDensityParallel();
 
    // write cumulative production and injection reports to output
    void outputCumLog(std::size_t reportStepNum,
                      const bool connData);
 
    // write production report to output
    void outputProdLog(std::size_t reportStepNum,
                       const bool connData);
 
    // write injection report to output
    void outputInjLog(std::size_t reportStepNum,
                      const bool connData);
 
    // write msw report to output
    void outputMSWLog(std::size_t reportStepNum);
 
    // calculate Initial Fluid In Place
    void calc_initial_inplace(const Parallel::Communication& comm);
 
    // calculate Fluid In Place
    Inplace calc_inplace(std::map<std::string, double>& miscSummaryData,
                         std::map<std::string, std::vector<double>>& regionData,
                         const Parallel::Communication& comm);
 
    void outputWellspecReport(const std::vector<std::string>& changedWells,
                              const bool changedWellLists,
                              const std::size_t reportStepNum,
                              const double elapsed,
                              boost::posix_time::ptime currentDate) const;
 
    void outputErrorLog(const Parallel::Communication& comm) const;
 
    void addRftDataToWells(data::Wells& wellDatas,
                           std::size_t reportStepNum,
                           const Parallel::Communication& comm)
    { this->rftC_.addToWells(wellDatas, reportStepNum, comm); }
 
    virtual void assignToSolution(data::Solution& sol);
 
    virtual void setRestart(const data::Solution& sol,
                    unsigned elemIdx,
                    unsigned globalDofIndex);
 
    Scalar getSolventSaturation(unsigned elemIdx) const
    {
        if (sSol_.size() > elemIdx)
            return sSol_[elemIdx];
 
        return 0;
    }
 
    Scalar getSolventRsw(unsigned elemIdx) const
    {
        if (rswSol_.size() > elemIdx)
            return rswSol_[elemIdx];
 
        return 0;
    }
 
    Scalar getPolymerConcentration(unsigned elemIdx) const
    {
        if (cPolymer_.size() > elemIdx)
            return cPolymer_[elemIdx];
 
        return 0;
    }
 
    Scalar getFoamConcentration(unsigned elemIdx) const
    {
        if (cFoam_.size() > elemIdx)
            return cFoam_[elemIdx];
 
        return 0;
    }
 
    Scalar getSaltConcentration(unsigned elemIdx) const
    {
        if (cSalt_.size() > elemIdx)
            return cSalt_[elemIdx];
 
        return 0;
    }
 
    Scalar getSaltSaturation(unsigned elemIdx) const
    {
        if (pSalt_.size() > elemIdx)
            return pSalt_[elemIdx];
 
        return 0;
    }
 
    Scalar getPermFactor(unsigned elemIdx) const
    {
        if (permFact_.size() > elemIdx)
            return permFact_[elemIdx];
 
        return 0;
    }
 
    const std::vector<Scalar>& getFluidPressure() const
    { return fluidPressure_; }
 
    const BioeffectsContainer<Scalar>& getBioeffects() const
    { return this->bioeffectsC_; }
 
    const CO2H2Container<Scalar>& getCO2H2() const
    { return this->CO2H2C_; }
 
    const FlowsContainer<FluidSystem>& getFlows() const
    { return this->flowsC_; }
 
    bool needInterfaceFluxes([[maybe_unused]] const bool isSubStep) const
    {
        return this->interRegionFlows_.wantInterRegflowSummary();
    }
 
    const std::map<std::pair<std::string, int>, double>& getBlockData()
    {
        return blockData_;
    }
 
    std::map<std::pair<std::string, int>, double>& getExtraBlockData()
    {
        return extraBlockData_;
    }
 
    // Per-LGR-cell block-summary results.  Cell identity is the
    // (keyword, grid level, level-local linearised Cartesian cell index)
    // tuple — the same identity data::Connection carries as
    // (lgr_grid, index) for LC* connection vectors.  Empty for runs
    // without LB* summary requests.
    const std::map<std::tuple<std::string, int, int>, double>& getLgrBlockData() const
    {
        return lgrBlockData_;
    }
 
    const Inplace* initialInplace() const
    {
        return this->initialInplace_.has_value()
            ? &*this->initialInplace_
            : nullptr;
    }
 
    bool localDataValid() const{
        return local_data_valid_;
    }
 
    void invalidateLocalData(){
        local_data_valid_ = false;
    }
 
    virtual void validateLocalData()
    {
        local_data_valid_ = true;
    }
 
    template<class Serializer>
    void serializeOp(Serializer& serializer)
    {
        serializer(initialInplace_);
    }
 
    RSTConv& getConv()
    { return this->rst_conv_; }
 
    const RSTConv& getConv() const
    { return this->rst_conv_; }
 
    const data::RegionVariableMapping& regVarMapping() const
    {
        return this->regVarMap_;
    }
 
    RegionVariableCollection& regionVariables()
    {
        return this->regionVars_;
    }
 
    void assignGlobalFieldsToSolution(data::Solution& sol);
 
protected:
    using ScalarBuffer = std::vector<Scalar>;
    using StringBuffer = std::vector<std::string>;
    enum { numPhases = FluidSystem::numPhases };
    enum { numComponents = FluidSystem::numComponents };
    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
    static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
    enum { gasCompIdx = FluidSystem::gasCompIdx };
    enum { oilCompIdx = FluidSystem::oilCompIdx };
    enum { waterCompIdx = FluidSystem::waterCompIdx };
    using Dir = FaceDir::DirEnum;
 
    template <class FluidState>
    Scalar hydroCarbonFraction(const FluidState& fs) const
    {
        if (this->eclState_.runspec().co2Storage()) {
            // CO2 storage uses the full pore volume.
            return 1.0;
        }
 
        // Sum the saturations of the active hydrocarbon phases.
        auto hydrocarbon = Scalar {0};
        if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
            hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
        }
 
        if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
            hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
        }
 
        return hydrocarbon;
    }
 
    template <class IntensiveQuantities>
    void updateTotalVolumesAndPressures_(const unsigned             globalDofIdx,
                                         const IntensiveQuantities& intQuants,
                                         const double               totVolume,
                                         const Scalar               referencePorosity)
    {
        const auto& fs = intQuants.fluidState();
 
        const double pv = totVolume * getValue(intQuants.porosity());
        const auto hydrocarbon = this->hydroCarbonFraction(fs);
 
        this->fipC_.assignPoreVolume(globalDofIdx,
                                    totVolume * referencePorosity,
                                    pv);
        if (! this->hydrocarbonPoreVolume_.empty()) {
            this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
        }
 
        if (!this->pressureTimesHydrocarbonVolume_.empty() &&
            !this->pressureTimesPoreVolume_.empty())
        {
            assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
            assert(this->fipC_.get(Inplace::Phase::PoreVolume).size() == this->pressureTimesPoreVolume_.size());
 
            if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
                this->pressureTimesPoreVolume_[globalDofIdx] =
                    getValue(fs.pressure(oilPhaseIdx)) * pv;
 
                this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
                    this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
            }
            else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
                this->pressureTimesPoreVolume_[globalDofIdx] =
                    getValue(fs.pressure(gasPhaseIdx)) * pv;
 
                this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
                    this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
            }
            else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
                this->pressureTimesPoreVolume_[globalDofIdx] =
                    getValue(fs.pressure(waterPhaseIdx)) * pv;
            }
        }
    }
 
    struct PhasePropertyNames
    {
        std::string_view oilDensity{};
        std::string_view gasDensity{};
        std::string_view waterDensity{};
        std::string_view oilViscosity{};
        std::string_view gasViscosity{};
        std::string_view waterViscosity{};
    };
 
    void assignPhaseProperties(data::Solution& sol,
                               const PhasePropertyNames& names);
 
    static bool allocBufferIfRequested(std::map<std::string, int>& rstKeywords,
                                       unsigned bufferSize,
                                       std::vector<Scalar>& buffer,
                                       std::string_view kw,
                                       bool supported,
                                       bool required = false);
 
    void assignBuffer(data::Solution& sol,
                      std::string_view name,
                      UnitSystem::measure measure,
                      std::vector<Scalar>& buffer);
 
    void assignBuffer(data::Solution& sol,
                      std::string_view name,
                      UnitSystem::measure measure,
                      std::vector<Scalar>& buffer,
                      int index);
 
    GenericOutputModule(const EclipseState& eclState,
                        const Schedule& schedule,
                        const SummaryConfig& summaryConfig,
                        const SummaryState& summaryState,
                        const std::string& moduleVersionName,
                        RSTConv::LocalToGlobalCellFunc globalCell,
                        std::function<bool(const unsigned)> isInterior,
                        const Parallel::Communication& comm,
                        bool enableEnergy,
                        bool constantTemperature,
                        bool enableMech,
                        bool enableSolvent,
                        bool enablePolymer,
                        bool enableFoam,
                        bool enableBrine,
                        bool enableSaltPrecipitation,
                        bool enableExtbo,
                        bool enableBioeffects,
                        bool enableGeochemistry);
 
    void doAllocBuffers(unsigned bufferSize,
                        unsigned reportStepNum,
                        const bool substep,
                        const bool log,
                        const bool forceRestartFieldAllocation,
                        const EclHysteresisConfig* hysteresisConfig,
                        unsigned numOutputNnc = 0,
                        std::map<std::string, int> rstKeywords = {});
 
    virtual void allocFormulationBuffers(std::map<std::string, int>& rstKeywords,
                                         unsigned bufferSize) = 0;
 
    void makeRegionSum(Inplace& inplace,
                       const std::string& region_name,
                       const Parallel::Communication& comm) const;
 
    Inplace accumulateRegionSums(const Parallel::Communication& comm);
 
    void updateSummaryRegionValues(const Inplace& inplace,
                                   std::map<std::string, double>& miscSummaryData,
                                   std::map<std::string, std::vector<double>>& regionData) const;
 
    static bool isOutputCreationDirective_(const std::string& keyword);
 
    // Sum Fip values over regions.
    static ScalarBuffer regionSum(const ScalarBuffer& property,
                                  const std::vector<int>& regionId,
                                  const std::size_t maxNumberOfRegions,
                                  const Parallel::Communication& comm);
 
    static int regionMax(const std::vector<int>& region,
                         const Parallel::Communication& comm);
 
    static void update(Inplace& inplace,
                       const std::string& region_name,
                       const Inplace::Phase phase,
                       const std::size_t ntFip,
                       const ScalarBuffer& values);
 
    static Scalar sum(const ScalarBuffer& v);
 
    void setupBlockData(std::function<bool(int)> isCartIdxOnThisRank);
    void setupExtraBlockData(const std::size_t        reportStepNum,
                             std::function<bool(int)> isCartIdxOnThisRank);
 
    // Allocate lgrBlockData_ slots for every Block-category summary node
    // that names an LGR cell.  The level component of the key is resolved
    // through EclipseGrid::get_lgr_cell_index (the same call data::Connection
    // uses for lgr_grid on LC* connection vectors); the level-local
    // linearised Cartesian index is recovered from node.number - 1.
    //
    // No rank check here — ownership for LGR leaf cells is decided in the
    // walk via Dune::InteriorEntity on the leaf entity; the existing
    // isCartIdxOnThisRank predicate operates in the global Cartesian space,
    // which LGR cells do not inhabit.
    //
    // Caller-supplied via parameter rather than read from a member to keep
    // this module decoupled from EclipseState ownership (matching how the
    // existing block path receives its rank predicate as an argument).
    void setupLgrBlockData(const std::map<std::string, int>& lgrNameToLevel,
                           const std::function<bool(int, int)>& isLgrCellOnThisRank);
 
    virtual bool isDefunctParallelWell(const std::string& wname) const = 0;
    virtual bool isOwnedByCurrentRank(const std::string& wname) const = 0;
    virtual bool isOnCurrentRank(const std::string& wname) const = 0;
 
    const EclipseState& eclState_;
    const Schedule& schedule_;
    const SummaryState& summaryState_;
 
    SummaryConfig summaryConfig_;
 
    InterRegFlowMap interRegionFlows_;
    LogOutputHelper<Scalar> logOutput_;
 
    bool enableEnergy_{false};
    bool constantTemperature_{false};
    bool enableMech_{false};
 
    bool enableSolvent_{false};
    bool enablePolymer_{false};
    bool enableFoam_{false};
    bool enableBrine_{false};
    bool enableSaltPrecipitation_{false};
    bool enableExtbo_{false};
    bool enableBioeffects_{false};
    bool enableGeochemistry_{false};
 
    bool forceDisableFipOutput_{false};
    bool forceDisableFipresvOutput_{false};
    bool computeFip_{false};
 
    RegionVariableCollection regionVars_;
 
    FIPContainer<FluidSystem> fipC_;
    std::unordered_map<std::string, std::vector<int>> regions_;
    std::unordered_map<Inplace::Phase, std::vector<SummaryConfigNode>> regionNodes_;
 
    std::vector<SummaryConfigNode> RPRNodes_;
    std::vector<SummaryConfigNode> RPRPNodes_;
 
    data::RegionVariableMapping regVarMap_{};
 
    std::vector<int> failedCellsPb_;
    std::vector<int> failedCellsPd_;
 
    ScalarBuffer gasFormationVolumeFactor_;
    ScalarBuffer hydrocarbonPoreVolume_;
    ScalarBuffer pressureTimesPoreVolume_;
    ScalarBuffer pressureTimesHydrocarbonVolume_;
    ScalarBuffer rPorV_;
    ScalarBuffer fluidPressure_;
    ScalarBuffer temperature_;
    ScalarBuffer rsw_;
    ScalarBuffer rvw_;
    ScalarBuffer overburdenPressure_;
    ScalarBuffer oilSaturationPressure_;
    ScalarBuffer drsdtcon_;
    ScalarBuffer sSol_;
    ScalarBuffer rswSol_;
    ScalarBuffer cPolymer_;
    ScalarBuffer cFoam_;
    ScalarBuffer cSalt_;
    ScalarBuffer pSalt_;
    ScalarBuffer permFact_;
    ExtboContainer<Scalar> extboC_;
    ScalarBuffer soMax_;
    ScalarBuffer swMax_;
    ScalarBuffer sgmax_;
    ScalarBuffer shmax_;
    ScalarBuffer somin_;
    ScalarBuffer swmin_;
    ScalarBuffer ppcw_;
    ScalarBuffer gasDissolutionFactor_;
    ScalarBuffer oilVaporizationFactor_;
    ScalarBuffer gasDissolutionFactorInWater_;
    ScalarBuffer waterVaporizationFactor_;
    ScalarBuffer rockCompPorvMultiplier_;
    ScalarBuffer minimumOilPressure_;
    ScalarBuffer saturatedOilFormationVolumeFactor_;
    ScalarBuffer rockCompTransMultiplier_;
    BioeffectsContainer<Scalar> bioeffectsC_;
    CO2H2Container<Scalar> CO2H2C_;
    ScalarBuffer pcgw_;
    ScalarBuffer pcow_;
    ScalarBuffer pcog_;
 
    // buffers for mechanical output
    MechContainer<Scalar> mech_;
    TpsaContainer<Scalar> tpsaC_;
 
    std::array<ScalarBuffer, numPhases> saturation_;
    std::array<ScalarBuffer, numPhases> density_;
    std::array<ScalarBuffer, numPhases> viscosity_;
 
    GeochemistryContainer<Scalar> geochemC_;
 
    TracerContainer<FluidSystem> tracerC_;
 
    std::array<ScalarBuffer, numPhases> residual_;
 
    FlowsContainer<FluidSystem> flowsC_;
 
    RFTContainer<FluidSystem> rftC_;
    RSTConv rst_conv_; 
 
    std::map<std::pair<std::string, int>, double> blockData_;
    // Extra block data required for non-summary output reasons
    // Example is the block pressures for RPTSCHED WELLS=2
    std::map<std::pair<std::string, int>, double> extraBlockData_;
 
    // Per-LGR-cell block-summary results.  Keyed by (keyword, grid level,
    // level-local linearised Cartesian cell index).  The map type matches
    // Summary::DynamicSimulatorState::LgrBlockValues so the writer can hand
    // a pointer to this member directly to the summary engine.  Empty for
    // runs that request no LB* summary vectors (zero non-LGR cost).
    std::map<std::tuple<std::string, int, int>, double> lgrBlockData_;
 
    std::optional<Inplace> initialInplace_;
    bool local_data_valid_{false};
 
    std::optional<RegionPhasePoreVolAverage> regionAvgDensity_;
 
private:
    void initialiseRegionVariableSupport();
};
 
} // namespace Opm
 
#endif // OPM_GENERIC_OUTPUT_MODULE_HPP