[WIP] [LGR] flow: parallel LGR INIT output from gathered simulator transmissibilities

CMakeLists_files.cmake

@@ -1036,6 +1036,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/HybridNewtonConfig.hpp
   opm/simulators/flow/InterRegFlows.hpp
   opm/simulators/flow/KeywordValidation.hpp
+  opm/simulators/flow/LgrOutputTransGather.hpp
   opm/simulators/flow/LogOutputHelper.hpp
   opm/simulators/flow/Main.hpp
   opm/simulators/flow/MechContainer.hpp

compareECLFiles.cmake

@@ -289,6 +289,7 @@ function(add_test_compare_parallel_simulation)
     DIR
     POSTFIX
     MPI_PROCS
+    COMPARE_MODE
   )
   set(multiValueArgs TEST_ARGS)
   cmake_parse_arguments(PARAM "$" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} )
@@ -313,6 +314,14 @@ function(add_test_compare_parallel_simulation)
     set(PARAM_SIMULATOR ${PARAM_DEV_SIMULATOR})
   endif()
 
+  # COMPARE_MODE "init": dry-run, compares EGRID+INIT only (parallel INIT-file regressions,
+  # e.g. parallel LGR transmissibility). Default (unset): compares SMRY+UNRST as before.
+  if(PARAM_COMPARE_MODE STREQUAL "init")
+    set(TEST_NAME_PREFIX compareParallelInitSim_)
+  else()
+    set(TEST_NAME_PREFIX compareParallelSim_)
+  endif()
+
   if(MPIEXEC_MAX_NUMPROCS GREATER_EQUAL MPI_PROCS)
     # Local computer system has at least ${MPI_PROCS} CPUs. Register test.
     set(RESULT_PATH ${BASE_RESULT_PATH}/parallel/${PARAM_SIMULATOR}+${PARAM_CASENAME})
@@ -325,17 +334,20 @@ function(add_test_compare_parallel_simulation)
                     -t ${PARAM_REL_TOL}
                     -c $<TARGET_FILE:compareECL>
                     -n ${MPI_PROCS})
+    if(PARAM_COMPARE_MODE STREQUAL "init")
+      list(APPEND DRIVER_ARGS -m init)
+    endif()
 
     # Add test that runs flow_mpi and outputs the results to file
-    opm_add_test(compareParallelSim_${PARAM_SIMULATOR}+${PARAM_FILENAME}${PARAM_POSTFIX}
+    opm_add_test(${TEST_NAME_PREFIX}${PARAM_SIMULATOR}+${PARAM_FILENAME}${PARAM_POSTFIX}
       EXE_TARGET
         ${PARAM_SIMULATOR}
       DRIVER_ARGS
         ${DRIVER_ARGS}
       TEST_ARGS
         ${TEST_ARGS}
     )
-    set_tests_properties(compareParallelSim_${PARAM_SIMULATOR}+${PARAM_FILENAME}${PARAM_POSTFIX}
+    set_tests_properties(${TEST_NAME_PREFIX}${PARAM_SIMULATOR}+${PARAM_FILENAME}${PARAM_POSTFIX}
                          PROPERTIES PROCESSORS ${MPI_PROCS})
   endif()
 endfunction()

opm/simulators/flow/EclGenericWriter.hpp

@@ -36,6 +36,7 @@
 #include <opm/simulators/flow/Transmissibility.hpp>
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
 
+#include <array>
 #include <map>
 #include <memory>
 #include <optional>
@@ -93,6 +94,24 @@ class EclGenericWriter
         globalTrans_ = globalTrans;
     }
 
+    // Parallel LGR INIT output: connection transmissibilities gathered from the
+    // ranks' own (distributed) simulator transmissibilities. When non-empty,
+    // computeTrans_ / exportNncStructure_ take values from these instead of
+    // querying a whole-grid transmissibility object. Complete on the I/O rank only.
+
+    // Same-level connections, keyed by (level, min/max level-Cartesian index).
+    void setGatheredLgrTrans(std::map<std::array<int,3>, double> gatheredTrans)
+    {
+        gatheredLgrTrans_ = std::move(gatheredTrans);
+    }
+
+    // Level-crossing connections (global<->LGR, LGR<->LGR NNCs), keyed by
+    // (smaller level, its level-Cartesian index, larger level, its level-Cartesian index).
+    void setGatheredLgrNncTrans(std::map<std::array<int,4>, double> gatheredNncTrans)
+    {
+        gatheredLgrNncTrans_ = std::move(gatheredNncTrans);
+    }
+
     void setSubStepReport(const SimulatorReportSingle& report)
     {
         sub_step_report_ = report;
@@ -118,6 +137,23 @@ class EclGenericWriter
     const TransmissibilityType& globalTrans() const;
     unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const;
 
+    // Output transmissibility value for a same-level connection: from the gathered
+    // per-rank simulator transmissibilities when available (parallel LGR runs),
+    // otherwise from the whole-grid transmissibility object (c1, c2).
+    double sameLevelOutputTrans_(int level,
+                                 int minLevelCartIdx,
+                                 int maxLevelCartIdx,
+                                 unsigned c1,
+                                 unsigned c2) const;
+
+    // Same, for a level-crossing (NNC) connection.
+    double crossLevelOutputTrans_(int smallerLevel,
+                                  int smallerLevelCartIdx,
+                                  int largerLevel,
+                                  int largerLevelCartIdx,
+                                  unsigned c1,
+                                  unsigned c2) const;
+
     void doWriteOutput(const int                          reportStepNum,
                        const std::optional<int>           timeStepNum,
                        const bool                         isSubStep,
@@ -165,6 +201,8 @@ class EclGenericWriter
     std::unique_ptr<TaskletRunner> taskletRunner_;
     Scalar restartTimeStepSize_;
     const TransmissibilityType* globalTrans_ = nullptr;
+    std::map<std::array<int,3>, double> gatheredLgrTrans_;
+    std::map<std::array<int,4>, double> gatheredLgrNncTrans_;
     const Dune::CartesianIndexMapper<Grid>& cartMapper_;
     const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper_;
     const EquilGrid* equilGrid_;

opm/simulators/flow/EclGenericWriter_impl.hpp

@@ -66,6 +66,7 @@
 #include <functional>
 #include <map>
 #include <memory>
+#include <stdexcept>
 #include <string>
 #include <unordered_map>
 #include <utility>
@@ -295,6 +296,8 @@ writeInit()
         }
         this->outputTrans_.reset();
     }
+    this->gatheredLgrTrans_.clear();
+    this->gatheredLgrNncTrans_.clear();
 }
 
 template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
@@ -574,12 +577,14 @@ computeTrans_(const std::vector<std::unordered_map<int,int>>&  levelCartToLevelC
             }
 
             if (maxLevelCartIdx - minLevelCartIdx == 1 && levelCartDims[0] > 1 ) {
-                outputTrans_->at(level).at("TRANX").template data<double>()[minLevelCartIdx] = globalTrans().transmissibility(c1, c2);
+                outputTrans_->at(level).at("TRANX").template data<double>()[minLevelCartIdx] =
+                    sameLevelOutputTrans_(level, minLevelCartIdx, maxLevelCartIdx, c1, c2);
                 continue; // skip other if clauses as they are false, last one needs some computation
             }
 
             if (maxLevelCartIdx - minLevelCartIdx == levelCartDims[0] && levelCartDims[1] > 1) {
-                outputTrans_->at(level).at("TRANY").template data<double>()[minLevelCartIdx] = globalTrans().transmissibility(c1, c2);
+                outputTrans_->at(level).at("TRANY").template data<double>()[minLevelCartIdx] =
+                    sameLevelOutputTrans_(level, minLevelCartIdx, maxLevelCartIdx, c1, c2);
                 continue; // skipt next if clause as it needs some computation
             }
 
@@ -588,7 +593,8 @@ computeTrans_(const std::vector<std::unordered_map<int,int>>&  levelCartToLevelC
                                          levelCartToLevelCompressed[level],
                                          minLevelCartIdx,
                                          maxLevelCartIdx)) {
-                outputTrans_->at(level).at("TRANZ").template data<double>()[minLevelCartIdx] = globalTrans().transmissibility(c1, c2);
+                outputTrans_->at(level).at("TRANZ").template data<double>()[minLevelCartIdx] =
+                    sameLevelOutputTrans_(level, minLevelCartIdx, maxLevelCartIdx, c1, c2);
             }
         }
     }
@@ -720,7 +726,8 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
                 const auto& [smallerLevel, smallerLevelCartIdx] = smallerPair;
                 const auto& [largerLevel, largerLevelCartIdx] = largerPair;
 
-                auto t = this->globalTrans().transmissibility(c1, c2);
+                auto t = this->crossLevelOutputTrans_(smallerLevel, smallerLevelCartIdx,
+                                                      largerLevel, largerLevelCartIdx, c1, c2);
 
                 // ECLIPSE ignores NNCs with zero transmissibility
                 // (different threshold than for NNC with corresponding
@@ -781,7 +788,7 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
                                           levelCartIdxIn, levelCartIdxOut)) {
                     // We need to check whether an NNC for this face was also
                     // specified via the NNC keyword in the deck.
-                    auto t = this->globalTrans().transmissibility(c1, c2);
+                    auto t = this->sameLevelOutputTrans_(level, levelCartIdxIn, levelCartIdxOut, c1, c2);
 
                     if (level == 0) {
                         auto candidate = std::lower_bound(nncData.begin(), nncData.end(),
@@ -879,7 +886,8 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
                     continue;
                 }
 
-                trans = this->globalTrans().transmissibility(c1, c2);
+                trans = this->sameLevelOutputTrans_(0, static_cast<int>(entry.cell1),
+                                                    static_cast<int>(entry.cell2), c1, c2);
 
                 if (! generatedNnc.empty()) {
                     for (const auto& generated : generatedNnc) {
@@ -1101,6 +1109,59 @@ evalSummary(const int                                            reportStepNum,
 #endif
 }
 
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+double
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+sameLevelOutputTrans_(int level,
+                      int minLevelCartIdx,
+                      int maxLevelCartIdx,
+                      unsigned c1,
+                      unsigned c2) const
+{
+    if (!gatheredLgrTrans_.empty()) {
+        const auto candidate = gatheredLgrTrans_.find(
+            std::array{level, minLevelCartIdx, maxLevelCartIdx});
+        if (candidate == gatheredLgrTrans_.end()) {
+            throw std::logic_error {
+                "Gathered LGR transmissibilities: no value for connection (level " +
+                std::to_string(level) + ", " + std::to_string(minLevelCartIdx) +
+                " -> " + std::to_string(maxLevelCartIdx) + ')'
+            };
+        }
+        return candidate->second;
+    }
+    return this->globalTrans().transmissibility(c1, c2);
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+double
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+crossLevelOutputTrans_(int smallerLevel,
+                       int smallerLevelCartIdx,
+                       int largerLevel,
+                       int largerLevelCartIdx,
+                       unsigned c1,
+                       unsigned c2) const
+{
+    // Mode test on the same-level map: it is non-empty in every parallel LGR run,
+    // while this map may legitimately be empty (it only guards against a silent
+    // fallback for a connection class the gather does cover).
+    if (!gatheredLgrTrans_.empty()) {
+        const auto candidate = gatheredLgrNncTrans_.find(
+            std::array{smallerLevel, smallerLevelCartIdx, largerLevel, largerLevelCartIdx});
+        if (candidate == gatheredLgrNncTrans_.end()) {
+            throw std::logic_error {
+                "Gathered LGR NNC transmissibilities: no value for connection (level " +
+                std::to_string(smallerLevel) + ", " + std::to_string(smallerLevelCartIdx) +
+                " -> level " + std::to_string(largerLevel) + ", " +
+                std::to_string(largerLevelCartIdx) + ')'
+            };
+        }
+        return candidate->second;
+    }
+    return this->globalTrans().transmissibility(c1, c2);
+}
+
 template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
 const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
 EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -52,6 +52,7 @@
 #include <opm/simulators/flow/FlowProblem.hpp>
 #include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
 #include <opm/simulators/flow/FlowThresholdPressure.hpp>
+#include <opm/simulators/flow/LgrOutputTransGather.hpp>
 #include <opm/simulators/flow/MixingRateControls.hpp>
 #include <opm/simulators/flow/OutputBlackoilModule.hpp>
 #include <opm/simulators/flow/VtkTracerModule.hpp>
@@ -72,6 +73,7 @@
 #include <stdexcept>
 #include <string>
 #include <string_view>
+#include <type_traits>
 #include <vector>
 
 namespace Opm {
@@ -347,8 +349,35 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // we try to avoid for the parallel running, has both global trans_ and transmissibilities_ allocated at the same time
         if (enableEclOutput_) {
             if (simulator.vanguard().grid().comm().size() > 1) {
-                if (simulator.vanguard().grid().comm().rank() == 0)
+                bool wholeGridTransNeeded = simulator.vanguard().grid().comm().rank() == 0;
+                // Parallel LGR: reuse the simulator's own (distributed) transmissibilities for the
+                // INIT output -- each rank contributes its interior connections, gathered on the
+                // I/O rank and keyed by level-Cartesian indices so the output walk over the global
+                // (equil) grid can look them up directly. This reuses the values already computed
+                // in parallel for the simulation itself instead of recomputing a whole-grid
+                // transmissibility. finishTransmissibilities() is idempotent; calling it here only
+                // moves the (anyway required) local trans build before the INIT write.
+                // grid().maxLevel() is identical on every rank: CpGrid::maxLevel counts the level
+                // grids, one per global LGR, so it is not this rank's partition depth.
+                if constexpr (std::is_same_v<GetPropType<TypeTag, Properties::Grid>, Dune::CpGrid>) {
+                    if (simulator.vanguard().grid().maxLevel() > 0) {
+                        finishTransmissibilities();
+                        const auto& localTrans = simulator.problem().eclTransmissibilities();
+                        auto gathered =
+                            gatherLgrOutputTrans(simulator.vanguard().grid(),
+                                                 simulator.vanguard().gridView(),
+                                                 [&localTrans](unsigned c1, unsigned c2)
+                                                 { return static_cast<double>(localTrans.transmissibility(c1, c2)); });
+                        eclWriter_->setGatheredLgrTrans(std::move(gathered.sameLevel));
+                        eclWriter_->setGatheredLgrNncTrans(std::move(gathered.crossLevel));
+                        // All output values (TRANX/Y/Z and NNC) come from the gathered maps --
+                        // no whole-grid transmissibility object is needed on the I/O rank.
+                        wholeGridTransNeeded = false;
+                    }
+                }
+                if (wholeGridTransNeeded) {
                     eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility());
+                }
             } else {
                 finishTransmissibilities();
                 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());