OutputExtractor.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_EXTRACTORS_HPP
#define OPM_OUTPUT_EXTRACTORS_HPP
 
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/Visitor.hpp>
#include <opm/common/utility/VoigtArray.hpp>
 
#include <opm/material/common/Valgrind.hpp>
 
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/basicproperties.hh>
#include <opm/models/utils/propertysystem.hh>
 
#include <algorithm>
#include <array>
#include <optional>
#include <span>
#include <unordered_map>
#include <set>
#include <variant>
#include <vector>
 
#include <fmt/format.h>
 
namespace Opm::detail {
 
template<class TypeTag>
struct Extractor
{
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using FluidState = typename IntensiveQuantities::FluidState;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    static constexpr int numPhases = FluidSystem::numPhases;
 
    struct HysteresisParams
    {
        Scalar somax{}; 
        Scalar swmax{}; 
        Scalar swmin{}; 
        Scalar sgmax{}; 
        Scalar shmax{}; 
        Scalar somin{}; 
    };
 
    struct Context
    {
        unsigned globalDofIdx; 
        unsigned pvtRegionIdx; 
        int episodeIndex;      
        const FluidState& fs;  
        const IntensiveQuantities& intQuants; 
        const HysteresisParams& hParams; 
    };
 
    using AssignFunc = std::function<void(const Context&)>;
 
    using ScalarFunc = std::function<Scalar(const Context&)>;
 
    using PhaseFunc = std::function<Scalar(const unsigned /*phase*/, const Context&)>;
 
    using ScalarBuffer = std::vector<Scalar>; 
    using PhaseArray = std::array<ScalarBuffer,numPhases>; 
 
    struct ScalarEntry
    {
        ScalarBuffer* data; 
        ScalarFunc extract; 
    };
 
    struct PhaseEntry
    {
        PhaseArray* data; 
        PhaseFunc extract; 
    };
 
    struct Entry
    {
        std::variant<AssignFunc, ScalarEntry, PhaseEntry> data; 
        bool condition = true; 
    };
 
    template<std::size_t size>
    static std::vector<Entry> removeInactive(std::array<Entry,size>& input)
    {
        // Setup active extractors
        std::vector<Entry> filtered_extractors;
        filtered_extractors.reserve(input.size());
        std::copy_if(std::move_iterator(input.begin()),
                     std::move_iterator(input.end()),
                     std::back_inserter(filtered_extractors),
                     [](const Entry& e)
                     {
                         if (!e.condition) {
                            return false;
                         }
                         return std::visit(VisitorOverloadSet{
                                               [](const AssignFunc&)
                                               {
                                                   return true;
                                               },
                                               [](const ScalarEntry& v)
                                               {
                                                   return !v.data->empty();
                                               },
                                               [](const PhaseEntry& v)
                                               {
                                                   return std::ranges::any_of(*v.data,
                                                                              [](const auto& ve)
                                                                              { return !ve.empty(); });
                                               }
                                           }, e.data);
                     });
 
        return filtered_extractors;
    }
 
    static void process(const Context& ectx,
                        const std::vector<Entry>& extractors)
    {
        std::ranges::for_each(extractors,
                      [&ectx](const auto& entry)
                      {
                          std::visit(VisitorOverloadSet{
                              [&ectx](const ScalarEntry& v)
                              {
                                  auto& array = *v.data;
                                  array[ectx.globalDofIdx] = v.extract(ectx);
                                  Valgrind::CheckDefined(array[ectx.globalDofIdx]);
                              },
                              [&ectx](const PhaseEntry& v)
                              {
                                  std::ranges::for_each(*v.data,
                                                [phaseIdx = 0, &ectx, &v](auto& array) mutable
                                                {
                                                    if (!array.empty()) {
                                                        array[ectx.globalDofIdx] = v.extract(phaseIdx, ectx);
                                                        Valgrind::CheckDefined(array[ectx.globalDofIdx]);
                                                    }
                                                    ++phaseIdx;
                                                });
                              },
                              [&ectx](const typename Extractor::AssignFunc& extract)
                              { extract(ectx); },
                          }, entry.data);
                      });
    }
};
 
template<class TypeTag>
struct BlockExtractor
{
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using FluidState = typename IntensiveQuantities::FluidState;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    static constexpr int numPhases = FluidSystem::numPhases;
    static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
    struct Context
    {
        unsigned globalDofIdx; 
        unsigned dofIdx;
        const FluidState& fs;  
        const IntensiveQuantities& intQuants; 
        const ElementContext& elemCtx;
    };
 
    using AssignFunc = std::function<void(const Context&)>;
 
    using ScalarFunc = std::function<Scalar(const Context&)>;
 
    using PhaseFunc = std::function<Scalar(const unsigned /*phase*/, const Context&)>;
 
    struct ScalarEntry
    {
        std::variant<std::string_view, std::vector<std::string_view>> kw;
 
        ScalarFunc extract;
    };
 
    struct PhaseEntry
    {
        std::variant<std::array<std::string_view, numPhases>,
                     std::array<std::array<std::string_view, numPhases>, 2>> kw;
 
        PhaseFunc extract;
    };
 
    using TensorFunc = std::function<Scalar(const VoigtIndex, const Context&)>;
 
    struct TensorEntry
    {
        std::string_view kw;
 
        TensorFunc extract;
    };
 
    using Entry = std::variant<ScalarEntry, PhaseEntry, TensorEntry>;
 
    struct Exec
    {
        Exec(double* d, ScalarFunc&& e)
            : data(d), extract(std::move(e))
        {}
 
        double* data; 
        ScalarFunc extract; 
    };
 
    using ExecMap = std::unordered_map<int, std::vector<Exec>>;
 
    static std::optional<ScalarFunc>
    makeExtractor(const std::string_view base_kw,
                  const std::span<const Entry> handlers)
    {
        using PhaseViewArray = std::array<std::string_view, numPhases>;
        using StringViewVec = std::vector<std::string_view>;
 
        // Suffix table position must equal the VoigtIndex enum value
        // (XX = 0, YY = 1, ZZ = 2, YZ = 3, XZ = 4, XY = 5).
        static constexpr auto voigtSuffixes = std::array<std::string_view, 6>{
            "XX", "YY", "ZZ", "YZ", "XZ", "XY",
        };
 
        unsigned phase{};
        unsigned voigt{};
        const auto handler_info =
            std::ranges::find_if(
                handlers,
                [&base_kw, &phase, &voigt](const auto& handler)
                {
                    // Tensor keywords are matched as base name + Voigt
                    // component suffix, without building the names.
                    if (const auto* tensor = std::get_if<TensorEntry>(&handler)) {
                        if (base_kw.size() != tensor->kw.size() + 2 ||
                            !base_kw.starts_with(tensor->kw))
                        {
                            return false;
                        }
                        const auto pos =
                            std::ranges::find(voigtSuffixes,
                                              base_kw.substr(tensor->kw.size()));
                        if (pos == voigtSuffixes.end()) {
                            return false;
                        }
                        voigt = std::distance(voigtSuffixes.begin(), pos);
                        return true;
                    }
 
                    // Extract list of keyword names from handler
                    const auto gen_handlers =
                        std::visit(VisitorOverloadSet{
                                      [](const ScalarEntry& entry)
                                      {
                                          return std::visit(VisitorOverloadSet{
                                                                [](const std::string_view& kw) -> StringViewVec
                                                                {
                                                                    return {kw};
                                                                },
                                                                [](const StringViewVec& kws) -> StringViewVec
                                                                { return kws; }
                                                            }, entry.kw);
                                      },
                                      [](const PhaseEntry& entry)
                                      {
                                          return std::visit(VisitorOverloadSet{
                                                                [](const PhaseViewArray& data)
                                                                {
                                                                    return StringViewVec{data.begin(), data.end()};
                                                                },
                                                                [](const std::array<PhaseViewArray,2>& data)
                                                                {
                                                                    StringViewVec res;
                                                                    res.reserve(2*numPhases);
                                                                    res.insert(res.end(),
                                                                               data[0].begin(),
                                                                               data[0].end());
                                                                    res.insert(res.end(),
                                                                               data[1].begin(),
                                                                               data[1].end());
                                                                    return res;
                                                                }
                                                            }, entry.kw);
                                      },
                                      [](const TensorEntry&) -> StringViewVec
                                      {
                                          // handled by the prefix/suffix match above
                                          return {};
                                      }
                                  }, handler);
 
                    const auto found_handler =
                        std::ranges::find(gen_handlers, base_kw);
                    if (found_handler != gen_handlers.end()) {
                        phase = std::distance(gen_handlers.begin(), found_handler) % numPhases;
                    }
                    return found_handler != gen_handlers.end();
                }
            );
 
        if (handler_info == handlers.end()) {
            return std::nullopt;
        }
 
        return std::visit(VisitorOverloadSet{
                              [](const ScalarEntry& e) -> ScalarFunc
                              {
                                  return e.extract;
                              },
                              [phase](const PhaseEntry& e) -> ScalarFunc
                              {
                                  return [phase, extract = e.extract]
                                         (const Context& ectx)
                                         {
                                             static constexpr auto phaseMap = std::array{
                                                 waterPhaseIdx,
                                                 oilPhaseIdx,
                                                 gasPhaseIdx,
                                             };
                                             return extract(phaseMap[phase], ectx);
                                         };
                              },
                              [voigt](const TensorEntry& e) -> ScalarFunc
                              {
                                  return [voigt, extract = e.extract]
                                         (const Context& ectx)
                                         {
                                             return extract(static_cast<VoigtIndex>(voigt), ectx);
                                         };
                              }
                          }, *handler_info);
    }
 
    static ExecMap setupExecMap(std::map<std::pair<std::string, int>, double>& blockData,
                                const std::span<const Entry> handlers)
    {
        ExecMap extractors;
 
        std::ranges::for_each(
            blockData,
            [&handlers, &extractors](auto& bd_info)
            {
                const auto& [key, cell] = bd_info.first;
                if (auto extract = makeExtractor(std::string_view{key}, handlers)) {
                    extractors[cell - 1].emplace_back(&bd_info.second, std::move(*extract));
                }
                else {
                    OpmLog::warning("Unhandled output keyword",
                                    fmt::format("Keyword '{}' is unhandled for output "
                                                "to summary file.", key));
                }
            }
        );
 
        return extractors;
    }
 
    static void process(const std::vector<Exec>& blockExtractors,
                        const Context& ectx)
    {
        std::ranges::for_each(blockExtractors,
                             [&ectx](auto& bdata)
                             { *bdata.data = bdata.extract(ectx); });
    }
 
    using LgrExecMap = std::unordered_map<int,
                                          std::unordered_map<int,
                                                             std::vector<Exec>>>;
 
    static LgrExecMap setupLgrExecMap(std::map<std::tuple<std::string,int,int>,
                                               double>& lgrBlockData,
                                      const std::span<const Entry> handlers)
    {
        LgrExecMap extractors;
 
        std::ranges::for_each(
            lgrBlockData,
            [&handlers, &extractors](auto& bd_info)
            {
                const auto& [lgr_kw, level, localCart] = bd_info.first;
                // LB* keywords mirror the global B* names; drop the leading
                // 'L' so the shared handler lookup sees "BPR", "BOSAT", ...
                const auto base_kw = std::string_view{lgr_kw}.substr(1);
                if (auto extract = makeExtractor(base_kw, handlers)) {
                    extractors[level][localCart].emplace_back(&bd_info.second, std::move(*extract));
                }
                else {
                    OpmLog::warning("Unhandled LGR output keyword",
                                    fmt::format("Keyword '{}' (LGR mirror of '{}') is "
                                                "unhandled for output to summary file.",
                                                lgr_kw, base_kw));
                }
            }
        );
 
        return extractors;
    }
};
 
} // namespace Opm::detail
 
#endif // OPM_OUTPUT_EXTRACTORS_HPP