From 77adc6ea3b606fb1faa792bf3b7e79c8d7ab3b88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 6 Jun 2023 18:45:10 +0200 Subject: [PATCH] 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 | 39 ++++ CMakeLists_files.cmake | 9 +- .../wells/ParallelPAvgDynamicSourceData.cpp | 156 ++++++++++++++ .../wells/ParallelPAvgDynamicSourceData.hpp | 173 +++++++++++++++ tests/test_parallel_wbp_sourcevalues.cpp | 201 ++++++++++++++++++ 5 files changed, 575 insertions(+), 3 deletions(-) create mode 100644 opm/simulators/wells/ParallelPAvgDynamicSourceData.cpp create mode 100644 opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp create mode 100644 tests/test_parallel_wbp_sourcevalues.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5d3a208b4..2d54c0270 100644 --- a/CMakeLists.txt +++ b/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} diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index fff1ba507..b8070f8ed 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/simulators/wells/ParallelPAvgDynamicSourceData.cpp b/opm/simulators/wells/ParallelPAvgDynamicSourceData.cpp new file mode 100644 index 000000000..1edf6da91 --- /dev/null +++ b/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 . +*/ + +#include + +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +Opm::ParallelPAvgDynamicSourceData:: +ParallelPAvgDynamicSourceData(const Parallel::Communication& comm, + const std::vector& 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& 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(this->localSrc_.size()), + this->src_.data(), // Output (to) + this->allSizes_.data(), // #elements per rank + this->startPointers_.data()); // Output offsets +} + +std::vector::size_type +Opm::ParallelPAvgDynamicSourceData:: +storageIndex(const std::vector::size_type elemIndex) const +{ + return this->storageIndex_[elemIndex]; +} + +void +Opm::ParallelPAvgDynamicSourceData:: +finaliseConstruction(const std::vector& 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 +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(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::size_type{0}; + for (const auto& elemIndex : allIndices) { + this->storageIndex_[elemIndex] = storageIx++; + } +} diff --git a/opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp b/opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp new file mode 100644 index 000000000..f5c8a73f3 --- /dev/null +++ b/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 . +*/ + +#ifndef PARALLEL_PAVG_DYNAMIC_SOURCE_DATA_HPP +#define PARALLEL_PAVG_DYNAMIC_SOURCE_DATA_HPP + +#include + +#include + +#include +#include +#include + +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; + + /// 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)>; + + /// 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& 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& 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 comm_; + + /// Subset of source locations owned by current rank. + std::vector locations_{}; + + /// Source data values owned by current rank. + std::vector localSrc_{}; + + /// Translation map from element index to storage index in + /// PAvgDynamicSourceData::src_. + std::vector::size_type> storageIndex_{}; + + /// Receive size from all ranks (allgatherv()). + std::vector allSizes_{}; // Type int to meet API requirements. + + /// Receive displacements for all ranks (allgatherv()). + std::vector 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::size_type + storageIndex(std::vector::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& 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 + 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 diff --git a/tests/test_parallel_wbp_sourcevalues.cpp b/tests/test_parallel_wbp_sourcevalues.cpp new file mode 100644 index 000000000..cbb1186ac --- /dev/null +++ b/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 . +*/ + +#include + +#define BOOST_TEST_MODULE Parallel_PAvg_Dynamic_Source_Data + +#define BOOST_TEST_NO_MAIN + +#include + +#include + +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +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 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(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 sourceLocations(const std::size_t numLoc) +{ + auto srcLoc = std::vector(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(i / this->size_) + : -1; + } + +private: + std::size_t rank_{}; + std::size_t size_{}; +}; + +class CalculateSourceTerm +{ +public: + using SrcTerm = Opm::PAvgDynamicSourceData::SourceDataSpan; + + 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::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(comm.rank()); + const auto comm_size = static_cast(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); +}