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/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 <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 Entry = std::variant<ScalarEntry, PhaseEntry>;

    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>>;

    template<std::size_t size>
    static std::optional<ScalarFunc>
    makeExtractor(const std::string_view base_kw,
                  const std::array<Entry,size>& handlers)
    {
        using PhaseViewArray = std::array<std::string_view, numPhases>;
        using StringViewVec = std::vector<std::string_view>;

        unsigned phase{};
        const auto handler_info =
            std::ranges::find_if(
                handlers,
                [&base_kw, &phase](const auto& handler)
                {
                    // 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);
                                      }
                                  }, 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);
                                         };
                              }
                          }, *handler_info);
    }

    template<std::size_t size>
    static ExecMap setupExecMap(std::map<std::pair<std::string, int>, double>& blockData,
                                const std::array<Entry,size>& 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>>>;

    template<std::size_t size>
    static LgrExecMap setupLgrExecMap(std::map<std::tuple<std::string,int,int>,
                                               double>& lgrBlockData,
                                      const std::array<Entry,size>& 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