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/GenericOutputModule.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 GenericOutputModule<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 = GenericOutputModule<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) override
    {
        this->compC_.outputRestart(sol, this->saturation_[oilPhaseIdx]);
        BaseType::assignToSolution(sol);
 
        // Use the compositional restart names for phase densities and viscosities.
        this->assignPhaseProperties(sol, typename BaseType::PhasePropertyNames {
            .oilDensity     = "DENO",
            .gasDensity     = "DENG",
            .waterDensity   = "DENW",
            .oilViscosity   = "VOIL",
            .gasViscosity   = "VGAS",
            .waterViscosity = "VWAT",
        });
 
        // Report compositional relative permeabilities under their restart names.
        using M = UnitSystem::measure;
        this->assignBuffer(sol, "KRO", M::identity,
                           relativePermeability_[oilPhaseIdx], oilPhaseIdx);
        this->assignBuffer(sol, "KRG", M::identity,
                           relativePermeability_[gasPhaseIdx], gasPhaseIdx);
        if constexpr (numPhases > 2) {
            this->assignBuffer(sol, "KRW", M::identity,
                               relativePermeability_[waterPhaseIdx], waterPhaseIdx);
        }
    }
 
    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)); });
                  }, this->compC_.moleFractionsAllocated()
            },
 
            // A phase with zero saturation has no defined composition; report zero instead
            // of stale flash values. A positive, however small, saturation is meaningful.
            Entry{[&compC = this->compC_](const ExtractContext& ectx)
                  {
                      const bool hasGas =
                          getValue(ectx.fs.saturation(gasPhaseIdx)) > Scalar{0};
                      compC.assignGasFractions(ectx.globalDofIdx,
                                               [&fs = ectx.fs, hasGas](const unsigned compIdx)
                                               {
                                                   return hasGas
                                                       ? getValue(fs.moleFraction(gasPhaseIdx, compIdx))
                                                       : Scalar{0};
                                               });
                  }, FluidSystem::phaseIsActive(gasPhaseIdx) &&
                     this->compC_.gasFractionsAllocated()
            },
            Entry{[&compC = this->compC_](const ExtractContext& ectx)
                  {
                      const bool hasOil =
                          getValue(ectx.fs.saturation(oilPhaseIdx)) > Scalar{0};
                      compC.assignOilFractions(ectx.globalDofIdx,
                                               [&fs = ectx.fs, hasOil](const unsigned compIdx)
                                               {
                                                   return hasOil
                                                       ? getValue(fs.moleFraction(oilPhaseIdx, compIdx))
                                                       : Scalar{0};
                                               });
                  }, FluidSystem::phaseIsActive(oilPhaseIdx) &&
                     this->compC_.oilFractionsAllocated()
            },
            Entry{[&compC = this->compC_](const ExtractContext& ectx)
                  {
                      compC.assignPhasePressures(ectx.globalDofIdx,
                                                 getValue(ectx.fs.pressure(oilPhaseIdx)),
                                                 getValue(ectx.fs.pressure(gasPhaseIdx)));
                  }, this->compC_.phasePressuresAllocated()
            },
            // Vapour mole fraction of the total mixture from the flash.
            Entry{[&compC = this->compC_](const ExtractContext& ectx)
                  {
                      const Scalar liquidFraction = getValue(ectx.fs.L());
                      compC.assignVaporFraction(ectx.globalDofIdx,
                                                std::clamp(Scalar{1} - liquidFraction,
                                                           Scalar{0}, Scalar{1}));
                  }, this->compC_.vaporFractionAllocated()
            },
            // The phase densities and viscosities, reported where the phase is present.
            Entry{PhaseEntry{&this->relativePermeability_,
                  [](const unsigned phaseIdx, const ExtractContext& ectx)
                  { return getValue(ectx.intQuants.relativePermeability(phaseIdx)); }}
            },
            Entry{PhaseEntry{&this->density_,
                  [](const unsigned phaseIdx, const ExtractContext& ectx)
                  {
                      return getValue(ectx.fs.saturation(phaseIdx)) > 0.0
                          ? getValue(ectx.fs.density(phaseIdx))
                          : Scalar{0};
                  }}
            },
            Entry{PhaseEntry{&this->viscosity_,
                  [](const unsigned phaseIdx, const ExtractContext& ectx)
                  {
                      return getValue(ectx.fs.saturation(phaseIdx)) > 0.0
                          ? getValue(ectx.fs.viscosity(phaseIdx))
                          : Scalar{0};
                  }}
            },
        };
 
        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);
    }
 
protected:
    void allocFormulationBuffers(std::map<std::string, int>& rstKeywords,
                                 const unsigned bufferSize) override
    {
        // Name each phase explicitly: the compositional phase ordering is not
        // the one the black-oil keyword names assume.
        const auto named = std::array{
            std::pair{static_cast<unsigned>(oilPhaseIdx),   std::string_view{"KRO"}},
            std::pair{static_cast<unsigned>(gasPhaseIdx),   std::string_view{"KRG"}},
            std::pair{static_cast<unsigned>(waterPhaseIdx), std::string_view{"KRW"}},
        };
        for (const auto& [phase, kw] : named) {
            if (phase >= numPhases || !FluidSystem::phaseIsActive(phase)) {
                continue;
            }
            BaseType::allocBufferIfRequested(rstKeywords, bufferSize,
                                             relativePermeability_[phase], kw, true);
        }
    }
 
private:
    using ScalarBuffer = typename BaseType::ScalarBuffer;
    std::array<ScalarBuffer, numPhases> relativePermeability_;
 
    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