Merge pull request #5053 from akva2/add_rstconv_class

added: class for calculating RPTRST CONV

CMakeLists.txt

@@ -519,6 +519,19 @@ opm_add_test(test_HDF5Serializer_Parallel
     4
 )
 
+opm_add_test(test_rstconv_parallel
+  EXE_NAME
+    test_rstconv
+  CONDITION
+    MPI_FOUND AND Boost_UNIT_TEST_FRAMEWORK_FOUND
+  DRIVER_ARGS
+    -n 4
+    -b ${PROJECT_BINARY_DIR}
+  NO_COMPILE
+  PROCESSORS
+    4
+)
+
 include(OpmBashCompletion)
 
 if (NOT BUILD_FLOW)

CMakeLists_files.cmake

@@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/flow/LogOutputHelper.cpp
   opm/simulators/flow/Main.cpp
   opm/simulators/flow/NonlinearSolverEbos.cpp
+  opm/simulators/flow/RSTConv.cpp
   opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.cpp
   opm/simulators/flow/SimulatorSerializer.cpp
   opm/simulators/flow/ValidationFunctions.cpp
@@ -276,6 +277,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_privarspacking.cpp
   tests/test_relpermdiagnostics.cpp
   tests/test_RestartSerialization.cpp
+  tests/test_rstconv.cpp
   tests/test_stoppedwells.cpp
   tests/test_timer.cpp
   tests/test_vfpproperties.cpp
@@ -443,6 +445,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/FlowMainEbos.hpp
   opm/simulators/flow/Main.hpp
   opm/simulators/flow/NonlinearSolverEbos.hpp
+  opm/simulators/flow/RSTConv.hpp
   opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp
   opm/simulators/flow/SimulatorSerializer.hpp
   opm/simulators/flow/KeywordValidation.hpp

opm/simulators/flow/RSTConv.cpp

@@ -0,0 +1,129 @@
+/*
+  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/flow/RSTConv.hpp>
+
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/input/eclipse/Schedule/RSTConfig.hpp>
+
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+
+namespace Opm {
+
+void RSTConv::init(const std::size_t numCells,
+                   const RSTConfig& rst_config,
+                   const std::array<int,3>& compIdx)
+{
+    const auto kw = rst_config.keywords.find("CONV");
+    if (kw == rst_config.keywords.end()) {
+        N_ = 0;
+        cnv_X_.clear();
+        return;
+    }
+
+    N_ = kw->second;
+    compIdx_ = compIdx;
+
+    cnv_X_.resize(3);
+    for (std::size_t i = 0; i < 3; ++i) {
+        if (compIdx_[i] != -1) {
+            cnv_X_[i].resize(numCells);
+        }
+    }
+}
+
+template<class ResidualVector>
+void RSTConv::update(const ResidualVector& resid)
+{
+    if (N_ == 0) {
+        return;
+    }
+
+    std::vector<int> cix(resid.size());
+    std::iota(cix.begin(), cix.end(), 0);
+    for (std::size_t i = 0; i < cnv_X_.size(); ++i) {
+        if (cnv_X_[i].empty()) {
+            continue;
+        }
+        std::partial_sort(cix.begin(), cix.begin() + N_, cix.end(),
+                          [&resid,c = compIdx_[i]](const int a, const int b)
+                          { return std::abs(resid[a][c]) > std::abs(resid[b][c]); });
+
+        this->gatherAndAccumulate(cix, resid, i);
+    }
+}
+
+template<class ResidualVector>
+void RSTConv::gatherAndAccumulate(const std::vector<int>& lIdx,
+                                  const ResidualVector& resid, int comp)
+{
+    if (comm_.size() == 1) {
+        for (int n = 0; n < N_; ++n) {
+            ++cnv_X_[comp][lIdx[n]];
+        }
+        return;
+    }
+
+    std::vector<double> values(comm_.rank() == 0 ? comm_.size() * N_ : N_);
+    std::vector<int> gIdx(comm_.rank() == 0 ? comm_.size() * N_ : N_);
+    for (int i = 0; i < N_; ++i) {
+        values[i] = std::abs(resid[lIdx[i]][compIdx_[comp]]);
+        gIdx[i] = globalCell_[lIdx[i]];
+    }
+
+    comm_.gather(gIdx.data(), gIdx.data(), N_, 0);
+    comm_.gather(values.data(), values.data(), N_, 0);
+    if (comm_.rank() != 0) {
+        return;
+    }
+
+    std::vector<int> valIdx(values.size());
+    std::iota(valIdx.begin(), valIdx.end(), 0);
+
+    std::partial_sort(valIdx.begin(), valIdx.begin() + N_, valIdx.end(),
+                      [&values](const int i1, const int i2)
+                      { return values[i1] > values[i2]; });
+
+    for (int n = 0; n < N_; ++n) {
+        ++cnv_X_[comp][gIdx[valIdx[n]]];
+    }
+}
+
+template<std::size_t Size>
+using BFV = Dune::BlockVector<Dune::FieldVector<double,Size>>;
+
+#define INSTANCE(SIZE) \
+    template void RSTConv::update<BFV<SIZE>>(const BFV<SIZE>&); \
+    template void RSTConv::gatherAndAccumulate<BFV<SIZE>>(const std::vector<int>&, \
+                                                          const BFV<SIZE>&, int);
+
+INSTANCE(1)
+INSTANCE(2)
+INSTANCE(3)
+INSTANCE(4)
+INSTANCE(5)
+INSTANCE(6)
+
+} // namespace Opm

opm/simulators/flow/RSTConv.hpp

@@ -0,0 +1,80 @@
+/*
+  Copyright 2023 SINTEF Digital
+
+  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 OPM_RST_CONV_HEADER_INCLUDED
+#define OPM_RST_CONV_HEADER_INCLUDED
+
+#include <array>
+#include <vector>
+
+#include <opm/simulators/utils/ParallelCommunication.hpp>
+
+namespace Opm {
+
+class RSTConfig;
+
+//! \brief Class computing RPTRST CONV output.
+class RSTConv
+{
+public:
+    //! \brief Constructor.
+    //! \param globalCell Mapping from local to global cell indices
+    //! \param comm Parallel communicator
+    RSTConv(const std::vector<int>& globalCell,
+            Parallel::Communication comm)
+        : globalCell_(globalCell)
+        , comm_(comm)
+    {}
+
+    //! \brief Init state at beginning of step.
+    //! \param numCells Global number of active cells in the model
+    //! \param rst_config RPTRST configuration
+    //! \param compIdx Component index for phases {OIL, GAS, WAT}, -1 if inactive
+    void init(const std::size_t numCells,
+              const RSTConfig& rst_config,
+              const std::array<int,3>& compIdx);
+
+    //! \brief Adds the CONV output for given residual vector.
+    template<class ResidualVector>
+    void update(const ResidualVector& resid);
+
+    //! \brief Obtain a const-ref to the accumulated data.
+    const std::vector<std::vector<int>>& getData() const
+    { return cnv_X_; }
+
+private:
+    //! \brief Gathers and accumulates results to the CONV arrays.
+    //! \param lIdx Vector of local indices (N first is used)
+    //! \param resid Residual vector
+    //! \param comp Component to consider
+    //! \details Handles parallel reduction
+    template<class ResidualVector>
+    void gatherAndAccumulate(const std::vector<int>& lIdx,
+                             const ResidualVector& resid, int comp);
+
+    const std::vector<int>& globalCell_; //!< Global cell indices
+    Parallel::Communication comm_; //!< Communicator
+    std::vector<std::vector<int>> cnv_X_; //!< Counts of worst cells for RPTRST CONV
+    std::array<int,3> compIdx_; //!< Component indices
+    int N_ = 0; //!< Number of cells to consider
+};
+
+} // namespace Opm
+
+#endif // OPM_RST_CONV_HEADER_INCLUDED

tests/test_rstconv.cpp

@@ -0,0 +1,194 @@
+/*
+  Copyright 2018 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2018 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 TestRstConv
+#define BOOST_TEST_NO_MAIN
+
+#include <boost/test/unit_test.hpp>
+#include <boost/version.hpp>
+#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 > 66
+#include <boost/test/data/test_case.hpp>
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+
+#include <opm/input/eclipse/Schedule/RSTConfig.hpp>
+#include <opm/simulators/flow/RSTConv.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <numeric>
+#include <random>
+
+#if HAVE_MPI
+struct MPIError
+{
+    MPIError(std::string s, int e) : errorstring(std::move(s)), errorcode(e){}
+    std::string errorstring;
+    int errorcode;
+};
+
+void MPI_err_handler(MPI_Comm*, int* err_code, ...)
+{
+    std::vector<char> err_string(MPI_MAX_ERROR_STRING);
+    int err_length;
+    MPI_Error_string(*err_code, err_string.data(), &err_length);
+    std::string s(err_string.data(), err_length);
+    std::cerr << "An MPI Error ocurred:" << std::endl << s << std::endl;
+    throw MPIError(s, *err_code);
+}
+#endif
+
+bool
+init_unit_test_func()
+{
+    return true;
+}
+
+struct TestCase {
+    std::size_t N;
+    std::array<int,3> phase;
+};
+
+std::ostream& operator<<(std::ostream& os, const TestCase& t)
+{
+    os << t.N;
+    for (int i : t.phase)
+        os << " " << i;
+    os << std::endl;
+
+    return os;
+}
+
+static const std::vector<TestCase> tests = {
+    {1, {0, 1, 2}},
+    {3, {0, 1, 2}},
+    {1, {0, -1, 1}},
+    {5, {0, -1, 1}},
+    {1, {-1, -1, 0}},
+    {4, {-1, -1, 0}},
+};
+
+#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 > 66
+BOOST_DATA_TEST_CASE(RstConvTest, tests)
+#else
+BOOST_AUTO_TEST_CASE(RstConvTest)
+#endif
+{
+#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 67
+    for (const auto& sample : tests) {
+#endif
+    const auto& cc = Dune::MPIHelper::getCommunication();
+
+    Opm::RSTConfig rst;
+    rst.keywords["CONV"] = sample.N;
+
+    std::vector<int> cellMapping(10);
+    std::iota(cellMapping.begin(), cellMapping.end(), cc.rank()*10);
+
+    Dune::BlockVector<Dune::FieldVector<double,3>> residual(10);
+
+    // generate data
+    std::vector<std::vector<int>> max(3);
+    if (cc.rank() == 0) {
+        std::random_device rng_device;
+        std::mt19937 mersenne_engine{rng_device()};
+        std::uniform_int_distribution<int> dist{0, 10*cc.size()-1};
+
+        for (int c = 0; c < 3; ++c) {
+            if (sample.phase[c] == -1) {
+                continue;
+            }
+            std::vector<int>& v = max[c];
+            while (v.size() < sample.N) {
+                int m = dist(mersenne_engine);
+                if (std::find(v.begin(), v.end(), m) == v.end()) {
+                    v.push_back(m);
+                }
+            }
+        }
+    }
+
+    for (int c = 0; c < 3; ++c) {
+        std::size_t size = max[c].size();
+        cc.broadcast(&size, 1, 0);
+        if (cc.rank() != 0) {
+            max[c].resize(size);
+        }
+        cc.broadcast(max[c].data(), max[c].size(), 0);
+    }
+
+    for (int i = 0; i < 10; ++i) {
+        for (int c = 0; c < 3; ++c) {
+            if (sample.phase[c] != -1) {
+                bool inMax = std::find(max[c].begin(),
+                                       max[c].end(),
+                                       cellMapping[i]) != max[c].end();
+                residual[i][sample.phase[c]] = inMax ? 1.0 : 1.0 / (i+2);
+            }
+        }
+    }
+
+    Opm::RSTConv cnv(cellMapping, cc);
+    cnv.init(10*cc.size(), rst, sample.phase);
+
+    cnv.update(residual);
+    cnv.update(residual);
+
+    if (cc.rank() == 0) {
+        BOOST_CHECK_EQUAL(cnv.getData().size(), 3);
+        BOOST_CHECK_EQUAL(cnv.getData()[0].size(),
+                          sample.phase[0] == -1 ? 0 : cc.size() * 10);
+        BOOST_CHECK_EQUAL(cnv.getData()[1].size(),
+                          sample.phase[1] == -1 ? 0 : cc.size() * 10);
+        BOOST_CHECK_EQUAL(cnv.getData()[2].size(),
+                          sample.phase[2] == -1 ? 0 : cc.size() * 10);
+
+        for (int i = 0; i < cc.size() * 10; ++i) {
+            for (int c = 0; c < 3; ++c) {
+                if (sample.phase[c] != -1) {
+                    bool inMax = std::find(max[c].begin(),
+                                           max[c].end(), i) != max[c].end();
+                    BOOST_CHECK_EQUAL(cnv.getData()[c][i], inMax ? 2 : 0);
+                }
+            }
+        }
+    }
+#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 67
+    }
+#endif
+}
+
+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
+    return boost::unit_test::unit_test_main(&init_unit_test_func, argc, argv);
+}