From be84969338ca1f1f7598a68029b2a62855fe4be2 Mon Sep 17 00:00:00 2001
From: Arne Morten Kvarving <arne.morten.kvarving@sintef.no>
Date: Mon, 11 Dec 2023 11:26:43 +0100
Subject: [PATCH] added: hook up RPTRST CONV output to restart file if
 requested

---
 ebos/collecttoiorank.hh                   |  3 ++
 ebos/eclgenericoutputblackoilmodule.cc    | 14 +++++++++
 ebos/eclgenericoutputblackoilmodule.hh    | 15 ++++++++--
 ebos/eclgenericwriter.hh                  |  5 ++++
 ebos/ecloutputblackoilmodule.hh           |  5 +---
 ebos/eclproblem.hh                        |  7 ++++-
 ebos/eclwriter.hh                         |  5 ++++
 opm/simulators/flow/BlackoilModelEbos.hpp | 36 +++++++++++++++++++++--
 8 files changed, 80 insertions(+), 10 deletions(-)

diff --git a/ebos/collecttoiorank.hh b/ebos/collecttoiorank.hh
index 5b5370ff7..4ae2b517d 100644
--- a/ebos/collecttoiorank.hh
+++ b/ebos/collecttoiorank.hh
@@ -91,6 +91,9 @@ public:
     const data::Solution& globalCellData() const
     { return globalCellData_; }
 
+    data::Solution& globalCellData()
+    { return globalCellData_; }
+
     const data::Wells& globalWellData() const
     { return globalWellData_; }
 
diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc
index 4f92597f6..2b98cd3c3 100644
--- a/ebos/eclgenericoutputblackoilmodule.cc
+++ b/ebos/eclgenericoutputblackoilmodule.cc
@@ -1499,6 +1499,20 @@ setupBlockData(std::function<bool(int)> isCartIdxOnThisRank)
     }
 }
 
+template<class FluidSystem, class Scalar>
+void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
+assignGlobalFieldsToSolution(data::Solution& sol)
+{
+    if (!this->cnvData_.empty()) {
+        constexpr std::array names = {"CNV_OIL", "CNV_GAS", "CNV_WAT"};
+        for (size_t i = 0; i < 3; ++i) {
+            if (!this->cnvData_[i].empty()) {
+                sol.insert(names[i], this->cnvData_[i], data::TargetType::RESTART_SOLUTION);
+            }
+        }
+    }
+}
+
 template class EclGenericOutputBlackoilModule<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,double>;
 
 } // namespace Opm
diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh
index 2a3151aa9..26a27262b 100644
--- a/ebos/eclgenericoutputblackoilmodule.hh
+++ b/ebos/eclgenericoutputblackoilmodule.hh
@@ -35,9 +35,7 @@
 #include <cstddef>
 #include <functional>
 #include <map>
-#include <numeric>
 #include <optional>
-#include <stdexcept>
 #include <unordered_map>
 #include <utility>
 #include <vector>
@@ -268,6 +266,11 @@ public:
         local_data_valid_ = true;
     }
 
+    void setCnvData(const std::vector<std::vector<int>>& data)
+    {
+        cnvData_ = data;
+    }
+
     // Virtual destructor for safer inheritance.
     virtual ~EclGenericOutputBlackoilModule();
 
@@ -277,6 +280,12 @@ public:
         serializer(initialInplace_);
     }
 
+    //! \brief Assign fields that are in global numbering to the solution.
+    //! \detail This is used to add fields that for some reason cannot be collected
+    //!         using the regular collect mechanism to the solution. In particular this
+    //!         is used with RPTRST CONV output.
+    void assignGlobalFieldsToSolution(data::Solution& sol);
+
 protected:
     using ScalarBuffer = std::vector<Scalar>;
     using StringBuffer = std::vector<std::string>;
@@ -490,6 +499,8 @@ protected:
     std::map<std::size_t, Scalar> gasConnectionSaturations_;
     std::map<std::pair<std::string, int>, double> blockData_;
 
+    std::vector<std::vector<int>> cnvData_; //!< Data for CNV_xxx arrays
+
     std::optional<Inplace> initialInplace_;
     bool local_data_valid_;
 };
diff --git a/ebos/eclgenericwriter.hh b/ebos/eclgenericwriter.hh
index f4d5c4696..73ff9fdc8 100644
--- a/ebos/eclgenericwriter.hh
+++ b/ebos/eclgenericwriter.hh
@@ -108,6 +108,11 @@ public:
         return outputNnc_;
     }
 
+    const CollectDataToIORankType& collectToIORank() const
+    {
+        return collectToIORank_;
+    }
+
 protected:
     const TransmissibilityType& globalTrans() const;
     unsigned int gridEquilIdxToGridIdx(unsigned int elemIndex) const;
diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh
index 35dc0ed6e..365d51538 100644
--- a/ebos/ecloutputblackoilmodule.hh
+++ b/ebos/ecloutputblackoilmodule.hh
@@ -29,6 +29,7 @@
 
 #include <dune/common/fvector.hh>
 
+#include <ebos/eclbasevanguard.hh>
 #include <ebos/eclgenericoutputblackoilmodule.hh>
 
 #include <opm/common/Exceptions.hpp>
@@ -53,12 +54,8 @@
 #include <array>
 #include <cassert>
 #include <cstddef>
-#include <initializer_list>
-#include <numeric>
-#include <optional>
 #include <stdexcept>
 #include <string>
-#include <tuple>
 #include <type_traits>
 #include <utility>
 #include <vector>
diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index 8bced897b..004e7adb1 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -760,7 +760,7 @@ public:
         if (enableDamarisOutput_) {
             damarisWriter_->writeOutput(localCellData, isSubStep) ;
         }
-#endif 
+#endif
         if (enableEclOutput_){
             eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep);
         }
@@ -1777,6 +1777,11 @@ public:
         return eclWriter_;
     }
 
+    void setConvData(const std::vector<std::vector<int>>& data)
+    {
+        eclWriter_->mutableEclOutputModule().setCnvData(data);
+    }
+
     template<class Serializer>
     void serializeOp(Serializer& serializer)
     {
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 548d74feb..cfb735168 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -431,6 +431,11 @@ public:
                                            /* interRegFlows = */ {},
                                            flowsn,
                                            floresn);
+            if (this->collectToIORank_.isIORank()) {
+                this->eclOutputModule_->assignGlobalFieldsToSolution(this->collectToIORank_.globalCellData());
+            }
+        } else {
+            this->eclOutputModule_->assignGlobalFieldsToSolution(localCellData);
         }
 
         if (this->collectToIORank_.isIORank()) {
diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp
index e305f7282..e514bd495 100644
--- a/opm/simulators/flow/BlackoilModelEbos.hpp
+++ b/opm/simulators/flow/BlackoilModelEbos.hpp
@@ -40,9 +40,10 @@
 #include <opm/simulators/aquifers/AquiferGridUtils.hpp>
 #include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
 #include <opm/simulators/flow/BlackoilModelEbosNldd.hpp>
+#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
 #include <opm/simulators/flow/countGlobalCells.hpp>
 #include <opm/simulators/flow/NonlinearSolverEbos.hpp>
-#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
+#include <opm/simulators/flow/RSTConv.hpp>
 #include <opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp>
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
@@ -233,6 +234,8 @@ namespace Opm {
         , phaseUsage_(phaseUsageFromDeck(eclState()))
         , param_( param )
         , well_model_ (well_model)
+        , rst_conv_(ebosSimulator_.problem().eclWriter()->collectToIORank().localIdxToGlobalIdxMapping(),
+                    grid_.comm())
         , terminal_output_ (terminal_output)
         , current_relaxation_(1.0)
         , dx_old_(ebosSimulator_.model().numGridDof())
@@ -303,6 +306,22 @@ namespace Opm {
 
             report.pre_post_time += perfTimer.stop();
 
+            auto getIdx = [](unsigned phaseIdx) -> int
+            {
+                if (FluidSystem::phaseIsActive(phaseIdx)) {
+                    const unsigned sIdx = FluidSystem::solventComponentIndex(phaseIdx);
+                    return Indices::canonicalToActiveComponentIndex(sIdx);
+                }
+
+                return -1;
+            };
+            const auto& schedule = ebosSimulator_.vanguard().schedule();
+            rst_conv_.init(ebosSimulator_.vanguard().globalNumCells(),
+                           schedule[timer.reportStepNum()].rst_config(),
+                           {getIdx(FluidSystem::oilPhaseIdx),
+                            getIdx(FluidSystem::gasPhaseIdx),
+                            getIdx(FluidSystem::waterPhaseIdx)});
+
             return report;
         }
 
@@ -375,14 +394,19 @@ namespace Opm {
                 convergence_reports_.back().report.reserve(11);
             }
 
+            SimulatorReportSingle result;
             if ((this->param_.nonlinear_solver_ != "nldd") ||
                 (iteration < this->param_.nldd_num_initial_newton_iter_))
             {
-                return this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
+                result = this->nonlinearIterationNewton(iteration, timer, nonlinear_solver);
             }
             else {
-                return this->nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
+                result = this->nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver);
             }
+
+            rst_conv_.update(ebosSimulator_.model().linearizer().residual());
+
+            return result;
         }
 
 
@@ -479,6 +503,7 @@ namespace Opm {
             Dune::Timer perfTimer;
             perfTimer.start();
             ebosSimulator_.problem().endTimeStep();
+            ebosSimulator_.problem().setConvData(rst_conv_.getData());
             report.pre_post_time += perfTimer.stop();
             return report;
         }
@@ -1056,6 +1081,9 @@ namespace Opm {
             }
         }
 
+        const std::vector<std::vector<int>>& getConvCells() const
+        { return rst_conv_.getData(); }
+
     protected:
         // ---------  Data members  ---------
 
@@ -1077,6 +1105,8 @@ namespace Opm {
         // Well Model
         BlackoilWellModel<TypeTag>& well_model_;
 
+        RSTConv rst_conv_; //!< Helper class for RPTRST CONV
+
         /// \brief Whether we print something to std::cout
         bool terminal_output_;
         /// \brief The number of cells of the global grid.