OutputCompositionalModule.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_OUTPUT_COMPOSITIONAL_MODULE_HPP
#define OPM_OUTPUT_COMPOSITIONAL_MODULE_HPP

#include <dune/grid/common/gridenums.hh>

#include <opm/simulators/utils/moduleVersion.hpp>

#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>

#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>

#include <opm/material/common/Valgrind.hpp>

#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/io/vtkblackoilparams.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>

#include <opm/simulators/flow/CompositionalContainer.hpp>
#include <opm/simulators/flow/FlowBaseVanguard.hpp>
#include <opm/simulators/flow/GenericOutputBlackoilModule.hpp>
#include <opm/simulators/flow/OutputExtractor.hpp>

#include <algorithm>
#include <cstddef>
#include <fstream>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>


namespace Opm {

// forward declaration
template <class TypeTag>
class EcfvDiscretization;

template <class TypeTag>
class OutputCompositionalModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
{
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Discretization = GetPropType<TypeTag, Properties::Discretization>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using BaseType = GenericOutputBlackoilModule<FluidSystem>;
    using Extractor = detail::Extractor<TypeTag>;

    enum { numPhases = FluidSystem::numPhases };
    enum { numComponents = FluidSystem::numComponents };
    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };

public:
    static void registerParameters()
    {
        BaseType::registerParameters();
        VtkBlackoilParams::registerParameters();
    }

    template <class CollectDataToIORankType>
    OutputCompositionalModule(const Simulator& simulator,
                              const SummaryConfig& smryCfg,
                              const CollectDataToIORankType& collectToIORank)
        : BaseType(simulator.vanguard().eclState(),
                   simulator.vanguard().schedule(),
                   smryCfg,
                   simulator.vanguard().summaryState(),
                   moduleVersionName(),
                   [this](const int idx)
                   { return simulator_.problem().eclWriter().collectOnIORank().localIdxToGlobalIdx(idx); },
                   [&collectToIORank](const int idx)
                   { return collectToIORank.isCartIdxOnThisRank(idx); },
                   simulator.vanguard().grid().comm(),
                   getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
                   getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::ConstantTemperature,
                   getPropValue<TypeTag, Properties::EnableMech>(),
                   getPropValue<TypeTag, Properties::EnableSolvent>(),
                   getPropValue<TypeTag, Properties::EnablePolymer>(),
                   getPropValue<TypeTag, Properties::EnableFoam>(),
                   getPropValue<TypeTag, Properties::EnableBrine>(),
                   getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
                   getPropValue<TypeTag, Properties::EnableExtbo>(),
                   getPropValue<TypeTag, Properties::EnableBioeffects>(),
                   getPropValue<TypeTag, Properties::EnableGeochemistry>())
        , simulator_(simulator)
    {
        for (auto& region_pair : this->regions_) {
            this->createLocalRegion_(region_pair.second);
        }

        auto isCartIdxOnThisRank = [&collectToIORank](const int idx) {
            return collectToIORank.isCartIdxOnThisRank(idx);
        };

        this->setupBlockData(isCartIdxOnThisRank);

        // Allocate LB* summary slots.  The compositional block-data fill
        // (processElementBlockData) is a stub, so no LB* values are produced
        // here and the parallel gather has nothing to collect; the allocation
        // is kept single-process and the ownership predicate is trivially true
        // (there are no per-rank values to keep disjoint).
        if (! collectToIORank.isParallel()) {
            // Empty name->level map: the compositional block-data fill is a stub,
            // so no LB* values are produced and nothing needs allocating.
            this->setupLgrBlockData({}, [](const int, const int) { return true; });
        }

        if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
            const std::string msg = "The output code does not support --owner-cells-first=false.";
            if (collectToIORank.isIORank()) {
                OpmLog::error(msg);
            }
            OPM_THROW_NOLOG(std::runtime_error, msg);
        }

        if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
            auto rset = this->eclState_.fieldProps().fip_regions();
            rset.push_back("PVTNUM");

            // Note: We explicitly use decltype(auto) here because the
            // default scheme (-> auto) will deduce an undesirable type.  We
            // need the "reference to vector" semantics in this instance.
            this->regionAvgDensity_
                .emplace(this->simulator_.gridView().comm(),
                         FluidSystem::numPhases, rset,
                         [fp = std::cref(this->eclState_.fieldProps())]
                         (const std::string& rsetName) -> decltype(auto)
                         { return fp.get().get_int(rsetName); });
        }
    }

    void
    allocBuffers(const unsigned bufferSize,
                 const unsigned reportStepNum,
                 const bool     substep,
                 const bool     log,
                 const bool     isRestart)
    {
        if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
            return;
        }

        auto rstKeywords = this->schedule_.rst_keywords(reportStepNum);
        this->compC_.allocate(bufferSize, rstKeywords);

        this->doAllocBuffers(bufferSize, reportStepNum, substep, log, isRestart,
                             /* hysteresisConfig = */ nullptr,
                             /* numOutputNnc =*/ 0,
                             std::move(rstKeywords));
    }

    void assignToSolution(data::Solution& sol)
    {
        this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]);
        BaseType::assignToSolution(sol);
    }

    void outputFipAndResvLog(const Inplace& inplace,
                             const std::size_t reportStepNum,
                             double elapsed,
                             boost::posix_time::ptime currentDate,
                             const bool substep,
                             const Parallel::Communication& comm)
    {
        if (comm.rank() != 0) {
            return;
        }

        std::unique_ptr<FIPConfig> fipSched;
        if (reportStepNum > 0) {
            const auto& rpt = this->schedule_[reportStepNum - 1].rpt_config.get();
            fipSched = std::make_unique<FIPConfig>(rpt);
        }

        const FIPConfig& fipc = reportStepNum == 0
            ? this->eclState_.getEclipseConfig().fip()
            : *fipSched;

        if (!substep && !this->forceDisableFipOutput_ && fipc.output(FIPConfig::OutputField::FIELD)) {
            this->logOutput_.timeStamp("BALANCE", elapsed, reportStepNum, currentDate);

            const auto& initial_inplace = *this->initialInplace();
            this->logOutput_.fip(inplace, initial_inplace, "");

            if (fipc.output(FIPConfig::OutputField::FIPNUM)) {
                this->logOutput_.fip(inplace, initial_inplace, "FIPNUM");

                if (fipc.output(FIPConfig::OutputField::RESV)) {
                    this->logOutput_.fipResv(inplace, "FIPNUM");
                }
            }

            if (fipc.output(FIPConfig::OutputField::FIP)) {
                for (const auto& reg : this->regions_) {
                    if (reg.first != "FIPNUM") {
                        std::ostringstream ss;
                        ss << "BAL" << reg.first.substr(3);
                        this->logOutput_.timeStamp(ss.str(), elapsed, reportStepNum, currentDate);
                        this->logOutput_.fip(inplace, initial_inplace, reg.first);

                        if (fipc.output(FIPConfig::OutputField::RESV)) {
                            this->logOutput_.fipResv(inplace, reg.first);
                        }
                    }
                }
            }
        }
    }

    void outputFipAndResvLogToCSV(const std::size_t reportStepNum,
                                  const bool substep,
                                  const Parallel::Communication& comm)
    {
        if (comm.rank() != 0) {
            return;
        }

        if ((reportStepNum == 0) && (!substep) &&
            (this->schedule_.initialReportConfiguration().has_value()) &&
            (this->schedule_.initialReportConfiguration()->contains("CSVFIP"))) {

            std::ostringstream csv_stream;

            this->logOutput_.csv_header(csv_stream);

            const auto& initial_inplace = *this->initialInplace();

            this->logOutput_.fip_csv(csv_stream, initial_inplace, "FIPNUM");

            for (const auto& reg : this->regions_) {
                if (reg.first != "FIPNUM") {
                    this->logOutput_.fip_csv(csv_stream, initial_inplace, reg.first);
                }
            }

            const IOConfig& io = this->eclState_.getIOConfig();
            auto csv_fname = io.getOutputDir() + "/" + io.getBaseName() + ".CSV";

            std::ofstream outputFile(csv_fname);
            outputFile << csv_stream.str();
        }
    }

    void setupExtractors(const bool         /*isSubStep*/,
                         const std::size_t  /*reportStepNum*/)
    {
        using Entry = typename Extractor::Entry;
        using ExtractContext = typename Extractor::Context;
        using ScalarEntry = typename Extractor::ScalarEntry;
        using PhaseEntry = typename Extractor::PhaseEntry;

        auto extractors = std::array{
            Entry{PhaseEntry{&this->saturation_,
                  [](const unsigned phase, const ExtractContext& ectx)
                  { return getValue(ectx.fs.saturation(phase)); }}
            },
            Entry{ScalarEntry{&this->fluidPressure_,
                  [](const ExtractContext& ectx)
                  {
                      if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
                          // Output oil pressure as default
                          return getValue(ectx.fs.pressure(oilPhaseIdx));
                      }
                      else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
                          // Output gas if oil is not present
                          return getValue(ectx.fs.pressure(gasPhaseIdx));
                      }
                      else {
                          // Output water if neither oil nor gas is present
                          return getValue(ectx.fs.pressure(waterPhaseIdx));
                      }
                  }}
            },
            Entry{ScalarEntry{&this->temperature_,
                  [](const ExtractContext& ectx)
                  { return getValue(ectx.fs.temperature(oilPhaseIdx)); }}
            },
            Entry{[&compC = this->compC_](const ExtractContext& ectx)
                  {
                      compC.assignMoleFractions(ectx.globalDofIdx,
                                                [&fs = ectx.fs](const unsigned compIdx)
                                                { return getValue(fs.moleFraction(compIdx)); });

                      if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
                          compC.assignGasFractions(ectx.globalDofIdx,
                                                   [&fs = ectx.fs](const unsigned compIdx)
                                                   { return getValue(fs.moleFraction(gasPhaseIdx, compIdx)); });
                      }

                      if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
                          compC.assignOilFractions(ectx.globalDofIdx,
                                                   [&fs = ectx.fs](const unsigned compIdx)
                                                   { return getValue(fs.moleFraction(oilPhaseIdx, compIdx)); });
                      }
                  }, this->compC_.allocated()
            },
        };

        this->extractors_ = Extractor::removeInactive(extractors);
    }

    void clearExtractors()
    { this->extractors_.clear(); }

    void processElement(const ElementContext& elemCtx)
    {
        OPM_TIMEBLOCK_LOCAL(processElement, Subsystem::Output);
        if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
            return;
        }

        typename Extractor::HysteresisParams hysterParams{};
        for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
            const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
            const auto& fs = intQuants.fluidState();

            const typename Extractor::Context ectx{
                elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0),
                0, // elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(),
                elemCtx.simulator().episodeIndex(),
                fs,
                intQuants,
                hysterParams
            };

            Extractor::process(ectx, extractors_);
        }
    }

    void processElementFlows(const ElementContext& /* elemCtx */)
    {
        OPM_TIMEBLOCK_LOCAL(processElementBlockData, Subsystem::Output);
        if (!std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>)
            return;
    }

    void processElementBlockData(const ElementContext& /* elemCtx */)
    {
        OPM_TIMEBLOCK_LOCAL(processElementBlockData, Subsystem::Output);
        if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
            return;
    }

    template <class ActiveIndex, class CartesianIndex>
    void processFluxes(const ElementContext& /* elemCtx */,
                       ActiveIndex&&         /* activeIndex*/,
                       CartesianIndex&&      /* cartesianIndex */)
    {
    }

    void initializeFluxData()
    {
        // Inter-region flow rates.  Note: ".clear()" prepares to accumulate
        // contributions per bulk connection between FIP regions.
        this->interRegionFlows_.clear();
    }

    void finalizeFluxData()
    {
        this->interRegionFlows_.compress();
    }

    const InterRegFlowMap& getInterRegFlows() const
    {
        return this->interRegionFlows_;
    }

    void updateFluidInPlace(const unsigned             /* globalDofIdx */,
                            const IntensiveQuantities& /* intQuants */,
                            const double               /* totVolume */)
    {
        // this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
    }

private:
    bool isDefunctParallelWell(const std::string& wname) const override
    {
        if (simulator_.gridView().comm().size() == 1)
            return false;
        const auto& parallelWells = simulator_.vanguard().parallelWells();
        std::pair<std::string, bool> value {wname, true};
        auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
        return candidate == parallelWells.end() || *candidate != value;
    }

    bool isOwnedByCurrentRank(const std::string& wname) const override
    {
        // Note: This statement is not correct for distributed wells and
        // will need additional logic once those are supported for
        // compositional flows.
        return ! this->isDefunctParallelWell(wname);
    }

    bool isOnCurrentRank(const std::string& wname) const override
    {
        // Note: This statement is not correct for distributed wells and
        // will need additional logic once those are supported for
        // compositional flows.
        return ! this->isDefunctParallelWell(wname);
    }

    void createLocalRegion_(std::vector<int>& region)
    {
        std::size_t elemIdx = 0;
        for (const auto& elem : elements(simulator_.gridView())) {
            if (elem.partitionType() != Dune::InteriorEntity) {
                region[elemIdx] = 0;
            }

            ++elemIdx;
        }
    }

    const Simulator& simulator_;
    CompositionalContainer<FluidSystem> compC_;
    std::vector<typename Extractor::Entry> extractors_;
};

} // namespace Opm

#endif // OPM_OUTPUT_COMPOSITIONAL_MODULE_HPP