Merge pull request #3032 from blattms/test-spe9-distributed-wells

Support external loadbalancers in EclCpGridVanguard

CMakeLists.txt

@@ -299,6 +299,18 @@ opm_add_test(flow_poly
   $<TARGET_OBJECTS:moduleVersion>)
 target_compile_definitions(flow_poly PRIVATE USE_POLYHEDRALGRID)
 
+# the production oriented general-purpose ECL simulator
+opm_add_test(flow_distribute_z
+  ONLY_COMPILE
+  ALWAYS_ENABLE
+  DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF}
+  DEPENDS opmsimulators
+  LIBRARIES opmsimulators
+  SOURCES
+  flow/flow_distribute_z.cpp
+  ${FLOW_TGTS}
+  $<TARGET_OBJECTS:moduleVersion>
+  )
 
 if (NOT BUILD_FLOW_VARIANTS)
   set(FLOW_VARIANTS_DEFAULT_ENABLE_IF "FALSE")

compareECLFiles.cmake

@@ -904,6 +904,14 @@ if(MPI_FOUND)
                                        REL_TOL ${rel_tol_parallel}
                                        TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8)
 
+  # A test for distributed standard wells. We load distribute only along the z-axis
+  add_test_compare_parallel_simulation(CASENAME spe9
+                                       FILENAME SPE9_CP_SHORT
+                                       SIMULATOR flow_distribute_z
+                                       ABS_TOL ${abs_tol_parallel}
+                                       REL_TOL ${rel_tol_parallel}
+                                       TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-8)
+
   add_test_compare_parallel_simulation(CASENAME spe9group
                                        FILENAME SPE9_CP_GROUP
                                        SIMULATOR flow

ebos/eclbasevanguard.hh

@@ -50,7 +50,6 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQState.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/UDQ/UDQConfig.hpp>
 
-#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
 
 #include <opm/simulators/utils/readDeck.hpp>
 
@@ -175,6 +174,9 @@ struct AllowDistributedWells<TypeTag, TTag::EclBaseVanguard> {
     static constexpr bool value = false;
 };
 
+template<class T1, class T2>
+struct UseMultisegmentWell;
+
 // Same as in BlackoilModelParametersEbos.hpp but for here.
 template<class TypeTag>
 struct UseMultisegmentWell<TypeTag, TTag::EclBaseVanguard> {

ebos/eclcpgridvanguard.hh

@@ -42,6 +42,7 @@
 #include <dune/common/version.hh>
 
 #include <sstream>
+#include <functional>
 
 namespace Opm {
 template <class TypeTag>
@@ -183,26 +184,31 @@ public:
             // not work)
             const auto& gridView = grid_->leafGridView();
             unsigned numFaces = grid_->numFaces();
-            std::vector<double> faceTrans(numFaces, 0.0);
-            ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
-            auto elemIt = gridView.template begin</*codim=*/0>();
-            const auto& elemEndIt = gridView.template end</*codim=*/0>();
-            for (; elemIt != elemEndIt; ++ elemIt) {
-                const auto& elem = *elemIt;
-                auto isIt = gridView.ibegin(elem);
-                const auto& isEndIt = gridView.iend(elem);
-                for (; isIt != isEndIt; ++ isIt) {
-                    const auto& is = *isIt;
-                    if (!is.neighbor())
-                        continue;
+            std::vector<double> faceTrans;
+            int loadBalancerSet = externalLoadBalancer_.has_value();
+            grid_->comm().broadcast(&loadBalancerSet, 1, 0);
+            if (!loadBalancerSet){
+                faceTrans.resize(numFaces, 0.0);
+                ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
+                auto elemIt = gridView.template begin</*codim=*/0>();
+                const auto& elemEndIt = gridView.template end</*codim=*/0>();
+                for (; elemIt != elemEndIt; ++ elemIt) {
+                    const auto& elem = *elemIt;
+                    auto isIt = gridView.ibegin(elem);
+                    const auto& isEndIt = gridView.iend(elem);
+                    for (; isIt != isEndIt; ++ isIt) {
+                        const auto& is = *isIt;
+                        if (!is.neighbor())
+                            continue;
 
-                    unsigned I = elemMapper.index(is.inside());
-                    unsigned J = elemMapper.index(is.outside());
+                        unsigned I = elemMapper.index(is.inside());
+                        unsigned J = elemMapper.index(is.outside());
 
-                    // FIXME (?): this is not portable!
-                    unsigned faceIdx = is.id();
+                        // FIXME (?): this is not portable!
+                        unsigned faceIdx = is.id();
 
-                    faceTrans[faceIdx] = globalTrans_->transmissibility(I, J);
+                        faceTrans[faceIdx] = globalTrans_->transmissibility(I, J);
+                    }
                 }
             }
 
@@ -222,9 +228,22 @@ public:
 
                     PropsCentroidsDataHandle<Dune::CpGrid> handle(*grid_, eclState, eclGrid, this->centroids_,
                                                                   cartesianIndexMapper());
-                    this->parallelWells_ = std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, serialPartitioning,
-                                                                          faceTrans.data(), ownersFirst, false, 1, true, zoltanImbalanceTol,
-                                                                          enableDistributedWells));
+                    if (loadBalancerSet)
+                    {
+                        std::vector<int> parts;
+                        if (grid_->comm().rank() == 0)
+                        {
+                            parts =  (*externalLoadBalancer_)(*grid_);
+                        }
+                        this->parallelWells_ = std::get<1>(grid_->loadBalance(handle, parts, &wells, ownersFirst, false, 1));
+                    }
+                    else
+                    {
+                        this->parallelWells_ =
+                            std::get<1>(grid_->loadBalance(handle, edgeWeightsMethod, &wells, serialPartitioning,
+                                                           faceTrans.data(), ownersFirst, false, 1, true, zoltanImbalanceTol,
+                                                           enableDistributedWells));
+                    }
                 }
                 catch(const std::bad_cast& e)
                 {
@@ -311,6 +330,13 @@ public:
         globalTrans_.reset();
     }
 
+    /// \brief Sets a function that returns external load balancing information when passed the grid
+    ///
+    /// The information is a vector of integers indication the partition index for each cell id.
+    static void setExternalLoadBalancer(const std::function<std::vector<int> (const Grid&)>& loadBalancer)
+    {
+        externalLoadBalancer_ = loadBalancer;
+    }
 protected:
     void createGrids_()
     {
@@ -378,9 +404,18 @@ protected:
     std::unique_ptr<CartesianIndexMapper> equilCartesianIndexMapper_;
 
     std::unique_ptr<EclTransmissibility<TypeTag> > globalTrans_;
+
+    /// \brief optional functor returning external load balancing information
+    ///
+    /// If it is set then this will be used during loadbalance.
+    static std::optional<std::function<std::vector<int> (const Grid&)>> externalLoadBalancer_;
     int mpiRank;
 };
 
+template<class TypeTag>
+std::optional<std::function<std::vector<int>(const typename EclCpGridVanguard<TypeTag>::Grid&)>>
+Opm::EclCpGridVanguard<TypeTag>::externalLoadBalancer_;
+
 } // namespace Opm
 
 #endif

flow/flow_distribute_z.cpp

@@ -0,0 +1,67 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015, 2017 IRIS AS
+  Copyright 2021 OPM-OP AS
+
+  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/Main.hpp>
+#include <ebos/eclcpgridvanguard.hh>
+
+std::vector<int> loadBalanceInZOnly(const Dune::CpGrid& grid)
+{
+    auto cartMapper = Dune::CartesianIndexMapper<Dune::CpGrid>(grid);
+    auto dims = cartMapper.cartesianDimensions();
+    std::vector<int> parts(grid.leafGridView().size(0));
+    auto numCellsPerProc = dims[2]/grid.comm().size();
+
+    if (grid.size(0)>0 && numCellsPerProc == 0)
+    {
+        OPM_THROW(std::logic_error,
+                  "cartesian grid must have more cells in z direction than number of processes.");
+    }
+
+    using ElementMapper =
+        Dune::MultipleCodimMultipleGeomTypeMapper<typename Dune::CpGrid::LeafGridView>;
+    const auto& gridView = grid.leafGridView();
+    const auto& idSet = grid.localIdSet();
+    ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
+
+    for( const auto &element : elements(gridView) )
+    {
+        const auto& id = idSet.id(element);
+        unsigned elemIdx = elemMapper.index(element);
+        const auto& cartIndex = cartMapper.cartesianIndex(elemIdx);
+        std::array<int,3> cartCoord;
+        cartMapper.cartesianCoordinate(cartIndex, cartCoord);
+        using std::min;
+        auto rank = min(grid.comm().size() -1, cartCoord[2] / numCellsPerProc);
+        parts[id] = rank;
+    }
+    return parts;
+}
+
+int main(int argc, char** argv)
+{
+    auto mainObject = Opm::Main(argc, argv);
+    Opm::EclCpGridVanguard<Opm::Properties::TTag::EclFlowProblem>::setExternalLoadBalancer(loadBalanceInZOnly);
+    return mainObject.runDynamic();
+}
+

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -32,10 +32,6 @@ namespace TTag {
 struct FlowModelParameters {};
 }
 
-// forward declaration to make this header usable from eclbasevanguard.hh
-template<class TypeTag, class MyTypeTag>
-struct EclDeckFileName;
-
 template<class TypeTag, class MyTypeTag>
 struct DbhpMaxRel {
     using type = UndefinedProperty;