Add Container for Dynamic WBPn Source Values

This commit adds a new container class,

    ParallelPAvgDynamicSourceData

which inherits from PAvgDynamicSourceData and provides a parallel
view of source contributions.  Member function

    collectLocalSources

will call the user-provided source term evaluation function for each
source location in its purview--typically those locations owned by
the current MPI rank.  Those values will be distributed to other MPI
ranks through member function synchroniseSources which will fill the
base class' 'src_' data member, and become available to clients
through read-only item spans.

CMakeLists.txt

@@ -363,6 +363,45 @@ opm_add_test(test_parallelwellinfo_mpi
   NO_COMPILE
     )
 
+opm_add_test(test_parallel_wbp_sourcevalues_np2
+  EXE_NAME
+    test_parallel_wbp_sourcevalues
+  CONDITION
+    MPI_FOUND AND Boost_UNIT_TEST_FRAMEWORK_FOUND
+  DRIVER_ARGS
+    -n 2
+    -b ${PROJECT_BINARY_DIR}
+  NO_COMPILE
+  PROCESSORS
+    2
+)
+
+opm_add_test(test_parallel_wbp_sourcevalues_np3
+  EXE_NAME
+    test_parallel_wbp_sourcevalues
+  CONDITION
+    MPI_FOUND AND Boost_UNIT_TEST_FRAMEWORK_FOUND
+  DRIVER_ARGS
+    -n 3
+    -b ${PROJECT_BINARY_DIR}
+  NO_COMPILE
+  PROCESSORS
+    3
+)
+
+opm_add_test(test_parallel_wbp_sourcevalues_np4
+  EXE_NAME
+    test_parallel_wbp_sourcevalues
+  CONDITION
+    MPI_FOUND AND Boost_UNIT_TEST_FRAMEWORK_FOUND
+  DRIVER_ARGS
+    -n 4
+    -b ${PROJECT_BINARY_DIR}
+  NO_COMPILE
+  PROCESSORS
+    4
+)
+
 opm_add_test(test_broadcast
   DEPENDS "opmsimulators"
   LIBRARIES opmsimulators ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}

CMakeLists_files.cmake

@@ -104,6 +104,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/MultisegmentWellGeneric.cpp
   opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
   opm/simulators/wells/MultisegmentWellSegments.cpp
+  opm/simulators/wells/ParallelPAvgDynamicSourceData.cpp
   opm/simulators/wells/ParallelWellInfo.cpp
   opm/simulators/wells/PerfData.cpp
   opm/simulators/wells/RateConverter.cpp
@@ -116,8 +117,8 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/StandardWellPrimaryVariables.cpp
   opm/simulators/wells/TargetCalculator.cpp
   opm/simulators/wells/VFPHelpers.cpp
-  opm/simulators/wells/VFPProdProperties.cpp
   opm/simulators/wells/VFPInjProperties.cpp
+  opm/simulators/wells/VFPProdProperties.cpp
   opm/simulators/wells/WellAssemble.cpp
   opm/simulators/wells/WellBhpThpCalculator.cpp
   opm/simulators/wells/WellConnectionAuxiliaryModule.cpp
@@ -227,7 +228,6 @@ endif()
 # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
 list (APPEND TEST_SOURCE_FILES
   tests/test_ALQState.cpp
-  tests/test_partitionCells.cpp
   tests/test_blackoil_amg.cpp
   tests/test_convergenceoutputconfiguration.cpp
   tests/test_convergencereport.cpp
@@ -244,7 +244,9 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
   tests/test_norne_pvt.cpp
+  tests/test_parallel_wbp_sourcevalues.cpp
   tests/test_parallelwellinfo.cpp
+  tests/test_partitionCells.cpp
   tests/test_preconditionerfactory.cpp
   tests/test_relpermdiagnostics.cpp
   tests/test_RestartSerialization.cpp
@@ -449,6 +451,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/utils/PropsCentroidsDataHandle.hpp
   opm/simulators/utils/SerializationPackers.hpp
   opm/simulators/utils/VectorVectorDataHandle.hpp
+  opm/simulators/utils/readDeck.hpp
   opm/simulators/wells/ALQState.hpp
   opm/simulators/wells/BlackoilWellModel.hpp
   opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -474,11 +477,11 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/MultisegmentWellGeneric.hpp
   opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
   opm/simulators/wells/MultisegmentWellSegments.hpp
+  opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp
   opm/simulators/wells/ParallelWellInfo.hpp
   opm/simulators/wells/PerfData.hpp
   opm/simulators/wells/PerforationData.hpp
   opm/simulators/wells/RateConverter.hpp
-  opm/simulators/utils/readDeck.hpp
   opm/simulators/wells/RegionAttributeHelpers.hpp
   opm/simulators/wells/RegionAverageCalculator.hpp
   opm/simulators/wells/SingleWellState.hpp

opm/simulators/wells/ParallelPAvgDynamicSourceData.cpp

@@ -0,0 +1,156 @@
+/*
+  Copyright 2023 Equinor ASA.
+
+  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 3 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/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <opm/grid/common/CommunicationUtils.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <functional>
+#include <iterator>
+#include <numeric>
+#include <utility>
+#include <vector>
+
+Opm::ParallelPAvgDynamicSourceData::
+ParallelPAvgDynamicSourceData(const Parallel::Communication&  comm,
+                              const std::vector<std::size_t>& sourceLocations,
+                              GlobalToLocal                   localCellIdx)
+    : PAvgDynamicSourceData { sourceLocations }
+    , comm_                 { comm }
+{
+    this->finaliseConstruction(sourceLocations, std::move(localCellIdx));
+}
+
+void Opm::ParallelPAvgDynamicSourceData::setToZero()
+{
+    std::fill_n(this->localSrc_.begin(), this->localSrc_.size(), 0.0);
+}
+
+void
+Opm::ParallelPAvgDynamicSourceData::
+reconstruct(const std::vector<std::size_t>& sourceLocations,
+            GlobalToLocal                   localCellIdx)
+{
+    PAvgDynamicSourceData::reconstruct(sourceLocations); // Reconstruct base
+    this->finaliseConstruction(sourceLocations, std::move(localCellIdx));
+}
+
+void Opm::ParallelPAvgDynamicSourceData::collectLocalSources(Evaluator eval)
+{
+    auto localIx = std::size_t{0};
+
+    for (const auto& location : this->locations_) {
+        eval(location.cell, this->localSourceTerm(localIx++));
+    }
+}
+
+void Opm::ParallelPAvgDynamicSourceData::synchroniseSources()
+{
+    this->comm_.get()
+        .allgatherv(this->localSrc_.data(),       // Input (from)
+                    static_cast<int>(this->localSrc_.size()),
+                    this->src_.data(),            // Output (to)
+                    this->allSizes_.data(),       // #elements per rank
+                    this->startPointers_.data()); // Output offsets
+}
+
+std::vector<double>::size_type
+Opm::ParallelPAvgDynamicSourceData::
+storageIndex(const std::vector<double>::size_type elemIndex) const
+{
+    return this->storageIndex_[elemIndex];
+}
+
+void
+Opm::ParallelPAvgDynamicSourceData::
+finaliseConstruction(const std::vector<std::size_t>& sourceLocations,
+                     GlobalToLocal                   localCellIdx)
+{
+    auto ix = std::size_t{0};
+
+    this->locations_.clear();
+
+    for (const auto& location : sourceLocations) {
+        if (const auto cell = localCellIdx(location); cell >= 0) {
+            this->locations_.push_back({ ix, cell });
+        }
+
+        ix += 1;
+    }
+
+    this->localSrc_.assign(numSpanItems() * this->locations_.size(), 0.0);
+
+    this->defineCommunication();
+}
+
+Opm::PAvgDynamicSourceData::SourceDataSpan<double>
+Opm::ParallelPAvgDynamicSourceData::localSourceTerm(const std::size_t localIx)
+{
+    return this->sourceTerm(localIx, this->localSrc_);
+}
+
+void Opm::ParallelPAvgDynamicSourceData::defineCommunication()
+{
+    // 1) Determine origins/owning ranks for all source terms.
+    auto ixVec = std::vector<std::size_t>(this->locations_.size());
+    std::transform(this->locations_.begin(), this->locations_.end(),
+                   ixVec.begin(),
+                   [](const auto& location) { return location.ix; });
+
+    constexpr auto numItems = numSpanItems();
+
+    const auto& [allIndices, allIxStart] = allGatherv(ixVec, this->comm_.get());
+
+    // -----------------------------------------------------------------------
+
+    // 2) Determine starting pointers/offsets/displacements for received
+    //    basic elements from each rank.  There are 'numItems' basic data
+    //    elements for each source term.
+    this->startPointers_.resize(allIxStart.size());
+    std::transform(allIxStart.begin(), allIxStart.end(),
+                   this->startPointers_.begin(),
+                   [](const int start)
+                   {
+                       return numItems * start;
+                   });
+
+    // -----------------------------------------------------------------------
+
+    // 3) Determine number of basic data elements to receive from each rank.
+    this->allSizes_.resize(allIxStart.size() - 1);
+    std::adjacent_difference(this->startPointers_.begin() + 1,
+                             this->startPointers_.end(),
+                             this->allSizes_.begin());
+
+    // -----------------------------------------------------------------------
+
+    // 4) Build translation mapping from source term element indices to
+    //    storage indices.
+    this->storageIndex_.resize(allIndices.size());
+    auto storageIx = std::vector<double>::size_type{0};
+    for (const auto& elemIndex : allIndices) {
+        this->storageIndex_[elemIndex] = storageIx++;
+    }
+}

opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp

@@ -0,0 +1,173 @@
+/*
+  Copyright 2023 Equinor ASA.
+
+  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 3 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/>.
+*/
+
+#ifndef PARALLEL_PAVG_DYNAMIC_SOURCE_DATA_HPP
+#define PARALLEL_PAVG_DYNAMIC_SOURCE_DATA_HPP
+
+#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <cstddef>
+#include <functional>
+#include <vector>
+
+namespace Opm {
+
+/// Dynamic source data for block-average pressure calculations.
+/// Specialisation for parallel runs.
+class ParallelPAvgDynamicSourceData : public PAvgDynamicSourceData
+{
+public:
+    /// Translate globally unique, linearised Cartesian cell indices to
+    /// local, on-rank, cell indices.  Assumed to return a negative value
+    /// result if the input cell index is not owned by the current rank.
+    using GlobalToLocal = std::function<int(const std::size_t)>;
+
+    /// Collect source term contributions from local, on-rank, cell.
+    ///
+    /// Called as
+    /// \code
+    ///    eval(cellIndex, sourceTerm)
+    /// \endcode
+    /// in which \c cellIndex is the local, on-rank, cell index in the range
+    /// 0 to #active cells on rank - 1.  Function \c eval is expected to
+    /// fill in/assign all \c sourceTerm items for this cell.
+    using Evaluator = std::function<void(int, SourceDataSpan<double>)>;
+
+    /// Constructor
+    ///
+    /// \param[in] comm MPI communication object.  Typically \code
+    ///   grid.comm() \endcode from the main simulation grid.
+    ///
+    /// \param[in] sourceLocations Known locations, typically linearised
+    ///   global call IDs, for which to enable collecting/reporting dynamic
+    ///   source data.  Typically \code allWBPCells() \endcode from a \c
+    ///   PAvgCalculatorCollection.
+    ///
+    /// \param[in] localCellIdx Translation from global, Cartesian cell
+    ///   indices to local, on-rank, cell indices.
+    ParallelPAvgDynamicSourceData(const Parallel::Communication&  comm,
+                                  const std::vector<std::size_t>& sourceLocations,
+                                  GlobalToLocal                   localCellIdx);
+
+    /// Clear contents of local source term contributions.
+    ///
+    /// Mostly useful when collecting source term contributions along the
+    /// well bore.
+    void setToZero();
+
+    /// Reconstruct Source Data backing storage and internal mapping tables
+    ///
+    /// Effectively replaces the original object formed by the constructor.
+    /// Mainly intended for updating objects as new wells and/or new
+    /// reservoir connections are introduced.
+    ///
+    /// \param[in] sourceLocations Known locations, typically linearised
+    ///   global call IDs, for which to enable collecting/reporting dynamic
+    ///   source data.  Typically \code allWBPCells() \endcode from a \c
+    ///   PAvgCalculatorCollection.
+    ///
+    /// \param[in] localCellIdx Translation from global, Cartesian cell
+    ///   indices to local, on-rank, cell indices.
+    void reconstruct(const std::vector<std::size_t>& sourceLocations,
+                     GlobalToLocal                   localCellIdx);
+
+    /// Compute local, on-rank, contributions to the collection of source
+    /// terms.
+    ///
+    /// \param[in] eval Source term evaluator object.
+    void collectLocalSources(Evaluator eval);
+
+    /// Exchange local contributions to build full, global view of all
+    /// source terms.
+    void synchroniseSources();
+
+private:
+    /// Identifier for local source term element.
+    struct LocalLocation
+    {
+        /// Source term element index.
+        std::size_t ix{};
+
+        /// Local cell index for this source term (0..#active on rank - 1).
+        int cell{};
+    };
+
+    /// MPI communication object.
+    std::reference_wrapper<const Parallel::Communication> comm_;
+
+    /// Subset of source locations owned by current rank.
+    std::vector<LocalLocation> locations_{};
+
+    /// Source data values owned by current rank.
+    std::vector<double> localSrc_{};
+
+    /// Translation map from element index to storage index in
+    /// PAvgDynamicSourceData::src_.
+    std::vector<std::vector<double>::size_type> storageIndex_{};
+
+    /// Receive size from all ranks (allgatherv()).
+    std::vector<int> allSizes_{}; // Type int to meet API requirements.
+
+    /// Receive displacements for all ranks (allgatherv()).
+    std::vector<int> startPointers_{}; // Type int to meet API requirements.
+
+    /// Translate element index into storage index.
+    ///
+    /// Customisation point.
+    ///
+    /// \param[in] elemIndex Source element index.
+    ///
+    /// \return Storage (starting) index in PAvgDynamicSourceData::src_.
+    [[nodiscard]] std::vector<double>::size_type
+    storageIndex(std::vector<double>::size_type elemIndex) const override;
+
+    /// Identify local source term elements on rank and build communication
+    /// pattern for all source terms.
+    ///
+    /// Assigns \c storageIndex_, \c allSizes_, and \c startPointers_.
+    ///
+    /// \param[in] sourceLocations Known locations, typically linearised
+    ///   global call IDs, for which to enable collecting/reporting dynamic
+    ///   source data.  Typically \code allWBPCells() \endcode from a \c
+    ///   PAvgCalculatorCollection.
+    ///
+    /// \param[in] localCellIdx Translation from global, Cartesian cell
+    ///   indices to local, on-rank, cell indices.
+    void finaliseConstruction(const std::vector<std::size_t>& sourceLocations,
+                              GlobalToLocal                   localCellIdx);
+
+    /// Form mutable data span into non-default backing store.
+    ///
+    /// \param[in] localIx Logical element index into \c localSrc_.
+    ///
+    /// \return Mutable view into \c localSrc_.
+    [[nodiscard]] SourceDataSpan<double>
+    localSourceTerm(const std::size_t localIx);
+
+    /// Build communication pattern for all source terms.
+    ///
+    /// Assigns \c storageIndex_, \c allSizes_, and \c startPointers_.
+    void defineCommunication();
+};
+
+} // namespace Opm
+
+#endif // PARALLEL_PAVG_DYNAMIC_SOURCE_DATA_HPP

tests/test_parallel_wbp_sourcevalues.cpp

@@ -0,0 +1,201 @@
+/*
+  Copyright 2023 Equinor.
+
+  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 3 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/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE Parallel_PAvg_Dynamic_Source_Data
+
+#define BOOST_TEST_NO_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <cstddef>
+#include <functional>
+#include <iostream>
+#include <numeric>
+#include <string>
+#include <string_view>
+#include <vector>
+
+namespace {
+
+#if HAVE_MPI
+struct MPIError
+{
+    MPIError(std::string_view errstr, const int ec)
+        : errorstring { errstr }
+        , errorcode   { ec }
+    {}
+
+    std::string errorstring;
+    int errorcode;
+};
+
+void MPI_err_handler(MPI_Comm*, int* err_code, ...)
+{
+    std::vector<char> err_string_vec(MPI_MAX_ERROR_STRING);
+    auto err_length = 0;
+
+    MPI_Error_string(*err_code, err_string_vec.data(), &err_length);
+
+    auto err_string = std::string_view {
+        err_string_vec.data(), static_cast<std::string_view::size_type>(err_length)
+    };
+
+    std::cerr << "An MPI Error ocurred:\n  -> " << err_string << '\n';
+
+    throw MPIError { err_string, *err_code };
+}
+#endif // HAVE_MPI
+
+bool init_unit_test_func()
+{
+    return true;
+}
+
+std::vector<std::size_t> sourceLocations(const std::size_t numLoc)
+{
+    auto srcLoc = std::vector<std::size_t>(numLoc);
+    std::iota(srcLoc.begin(), srcLoc.end(), std::size_t{0});
+    return srcLoc;
+}
+
+class LocalCellIndex
+{
+public:
+    explicit LocalCellIndex(const std::size_t rank,
+                            const std::size_t size)
+        : rank_ { rank }
+        , size_ { size }
+    {}
+
+    int operator()(const std::size_t i) const
+    {
+        return ((i % this->size_) == this->rank_)
+            ? static_cast<int>(i / this->size_)
+            : -1;
+    }
+
+private:
+    std::size_t rank_{};
+    std::size_t size_{};
+};
+
+class CalculateSourceTerm
+{
+public:
+    using SrcTerm = Opm::PAvgDynamicSourceData::SourceDataSpan<double>;
+
+    explicit CalculateSourceTerm(const std::size_t rank)
+        : rank_ { rank }
+    {}
+
+    void operator()(const int i, SrcTerm source_term) const
+    {
+        using Item = typename SrcTerm::Item;
+
+        source_term
+            .set(Item::Pressure      , this->rank_*314.15 -      i)
+            .set(Item::PoreVol       , this->rank_*172.9  + 10.0*i)
+            .set(Item::MixtureDensity, this->rank_*852.96 +      i);
+    }
+
+private:
+    std::size_t rank_{};
+};
+
+std::size_t
+sourceTermsAreCorrect(const std::size_t                 comm_size,
+                      const std::size_t                 num_src,
+                      const Opm::PAvgDynamicSourceData& source_data)
+{
+    using Item = Opm::PAvgDynamicSourceData::SourceDataSpan<const double>::Item;
+
+    auto num_correct = 0*num_src;
+
+    for (auto srcID = 0*num_src; srcID < num_src; ++srcID) {
+        const auto rank = srcID % comm_size;
+        const auto locI = srcID / comm_size;
+
+        const auto src = source_data[srcID];
+        const auto ok =
+            (src[Item::Pressure]       == rank*314.15 -    locI) &&
+            (src[Item::PoreVol]        == rank*172.9  + 10*locI) &&
+            (src[Item::MixtureDensity] == rank*852.96 +    locI);
+
+        num_correct += ok;
+    }
+
+    return num_correct == num_src;
+}
+
+} // Anonymous namespace
+
+BOOST_AUTO_TEST_CASE(Eval_and_collect)
+{
+    auto comm = Opm::Parallel::Communication {
+        Dune::MPIHelper::getCommunicator()
+    };
+
+    const auto comm_rank = static_cast<std::size_t>(comm.rank());
+    const auto comm_size = static_cast<std::size_t>(comm.size());
+
+    const auto num_src_loc = std::size_t{50};
+
+    auto source_data = Opm::ParallelPAvgDynamicSourceData {
+        comm, sourceLocations(num_src_loc),
+        LocalCellIndex { comm_rank, comm_size }
+    };
+
+    source_data.collectLocalSources(CalculateSourceTerm { comm_rank });
+    source_data.synchroniseSources();
+
+    const auto num_rank_correct = comm.sum
+        (sourceTermsAreCorrect(comm_size, num_src_loc, source_data));
+
+    if (comm_rank == 0) {
+        BOOST_CHECK_EQUAL(num_rank_correct, comm_size);
+    }
+}
+
+int main(int argc, char** argv)
+{
+    Dune::MPIHelper::instance(argc, argv);
+
+#if HAVE_MPI
+    // Register a throwing error handler to allow for debugging with
+    //
+    //   catch throw
+    //
+    // in GDB.
+    MPI_Errhandler handler{};
+    MPI_Comm_create_errhandler(MPI_err_handler, &handler);
+    MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler);
+#endif // HAVE_MPI
+
+    return boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv);
+}