Report where time goes inside the well solves

opm/models/discretization/common/fvbasediscretization.hh

@@ -67,8 +67,10 @@
 #include <algorithm>
 #include <array>
 #include <cstddef>
+#include <exception>
 #include <list>
 #include <memory>
+#include <mutex>
 #include <stdexcept>
 #include <sstream>
 #include <string>
@@ -726,49 +728,86 @@ public:
     {
         invalidateIntensiveQuantitiesCache(timeIdx);
 
+        // exceptions must not escape the parallel block below (that calls
+        // std::terminate()); tuck any exception away and rethrow it after the
+        // block, so that e.g. a failed flash in the property evaluation leads
+        // to a time step chop instead of an abort
+        std::mutex exceptionLock;
+        std::exception_ptr exceptionPtr = nullptr;
+
         // loop over all elements...
         ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
         {
-            ElementContext elemCtx(simulator_);
-            ElementIterator elemIt = threadedElemIt.beginParallel();
-            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
-                const Element& elem = *elemIt;
-                elemCtx.updatePrimaryStencil(elem);
-                elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
+            try {
+                ElementContext elemCtx(simulator_);
+                for (ElementIterator elemIt = threadedElemIt.beginParallel();
+                     !threadedElemIt.isFinished(elemIt);
+                     elemIt = threadedElemIt.increment())
+                {
+                    const Element& elem = *elemIt;
+                    elemCtx.updatePrimaryStencil(elem);
+                    elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
+                }
+            }
+            catch (...) {
+                std::lock_guard<std::mutex> take(exceptionLock);
+                exceptionPtr = std::current_exception();
+                threadedElemIt.setFinished();
             }
         }
+
+        if (exceptionPtr) {
+            std::rethrow_exception(exceptionPtr);
+        }
     }
 
     template <class GridViewType>
     void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
     {
+        // see the overload above for why exceptions are bridged out of the
+        // parallel block like this
+        std::mutex exceptionLock;
+        std::exception_ptr exceptionPtr = nullptr;
+
         // loop over all elements...
         ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
         {
-
-            ElementContext elemCtx(simulator_);
-            auto elemIt = threadedElemIt.beginParallel();
-            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
-                if (elemIt->partitionType() != Dune::InteriorEntity) {
-                    continue;
-                }
-                const Element& elem = *elemIt;
-                elemCtx.updatePrimaryStencil(elem);
-                // Mark cache for this element as invalid.
-                const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
-                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
-                    const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
-                    setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false);
+            try {
+                ElementContext elemCtx(simulator_);
+                for (auto elemIt = threadedElemIt.beginParallel();
+                     !threadedElemIt.isFinished(elemIt);
+                     elemIt = threadedElemIt.increment())
+                {
+                    if (elemIt->partitionType() != Dune::InteriorEntity) {
+                        continue;
+                    }
+                    const Element& elem = *elemIt;
+                    elemCtx.updatePrimaryStencil(elem);
+                    // Mark cache for this element as invalid.
+                    const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
+                    for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
+                        const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+                        setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false);
+                    }
+                    // Update for this element.
+                    elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
                 }
-                // Update for this element.
-                elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
             }
+            catch (...) {
+                std::lock_guard<std::mutex> take(exceptionLock);
+                exceptionPtr = std::current_exception();
+                threadedElemIt.setFinished();
+            }
+        }
+
+        if (exceptionPtr) {
+            std::rethrow_exception(exceptionPtr);
         }
     }
 

opm/simulators/flow/Banners.cpp

@@ -121,15 +121,16 @@ void printFlowTrailer(int nprocs,
                       const double total_setup_time,
                       const double deck_read_time,
                       const SimulatorReport& report,
-                      const std::string_view extra_summary)
+                      const std::string_view extra_summary,
+                      const bool performance_details)
 {
     std::ostringstream ss;
     ss << "\n\n================    End of simulation     ===============\n\n";
     ss << fmt::format("Number of MPI processes: {:9}\n", nprocs);
     ss << fmt::format("Threads per MPI process: {:9}\n", nthreads);
     ss << fmt::format("Setup time:                 {:9.2f} s\n", total_setup_time);
     ss << fmt::format("  Deck input:               {:9.2f} s\n", deck_read_time);
-    report.reportFullyImplicit(ss);
+    report.reportFullyImplicit(ss, performance_details);
     if (!extra_summary.empty()) {
         if (extra_summary.front() != '\n') {
             ss << '\n';

opm/simulators/flow/Banners.hpp

@@ -46,7 +46,8 @@ void printFlowTrailer(int nprocs,
                       const double total_setup_time,
                       const double deck_read_time,
                       const SimulatorReport& report,
-                      std::string_view extra_summary);
+                      std::string_view extra_summary,
+                      const bool performance_details = false);
 
 } // namespace Opm
 

opm/simulators/flow/ConvergenceOutputConfiguration.cpp

@@ -57,15 +57,15 @@ namespace {
                 fmt::format("Unsupported convergence output "
                             "option value{}: {}\n"
                             "Supported values are \"none\", "
-                            "\"steps\", and \"iterations\"",
+                            "\"steps\", \"iterations\", and \"performance\"",
                             pl, fmt::join(unsupp.begin(), u, ", "))
             };
         }
 
         throw std::invalid_argument {
             fmt::format("Option {}:\n - Unsupported value{}: {}\n"
                         " - Supported values are \"none\", "
-                        "\"steps\", and \"iterations\"",
+                        "\"steps\", \"iterations\", and \"performance\"",
                         optionName, pl,
                         fmt::join(unsupp.begin(), u, ", "))
         };
@@ -79,11 +79,13 @@ namespace {
         auto opt = std::vector<Option>{};
 
         const auto values = std::unordered_map<std::string, Option> {
-            { "none"      , Option::None       },
-            { "step"      , Option::Steps      }, // Alias for 'steps' (plural)
-            { "steps"     , Option::Steps      },
-            { "iteration" , Option::Iterations }, // Alias for 'iterations' (plural)
-            { "iterations", Option::Iterations },
+            { "none"       , Option::None        },
+            { "step"       , Option::Steps       }, // Alias for 'steps' (plural)
+            { "steps"      , Option::Steps       },
+            { "iteration"  , Option::Iterations  }, // Alias for 'iterations' (plural)
+            { "iterations" , Option::Iterations  },
+            { "perf"       , Option::Performance }, // Alias for 'performance'
+            { "performance", Option::Performance },
         };
 
         auto unsupp = std::vector<std::string>{};

opm/simulators/flow/ConvergenceOutputConfiguration.hpp

@@ -39,6 +39,9 @@ namespace Opm {
 ///   * "iterations" -- Want additional convergence output pertaining to each
 ///                     non-linar ieration in each timestep.
 ///
+///   * "performance" -- Want additional performance details in the summary
+///                      at the end of the simulation.
+///
 /// Option value "none" overrides all other options.  In other words, if the
 /// user requests "none", then there will be no additional convergence
 /// output, even if there are other options in the option string.
@@ -54,6 +57,7 @@ class ConvergenceOutputConfiguration
         None = 0,
         Steps = 1 << 1,
         Iterations = 1 << 2,
+        Performance = 1 << 3,
     };
 
     /// Constructor

opm/simulators/flow/FlowGenericProblem_impl.hpp

@@ -81,9 +81,25 @@ FlowGenericProblem(const EclipseState& eclState,
     // 2. EQLDIMS item 2.  Default value from
     //    opm-common/opm/input/eclipse/share/keywords/000_Eclipse100/E/EQLDIMS
 
-    numPressurePointsEquil_ = Parameters::IsSet<Parameters::NumPressurePointsEquil>()
-        ? Parameters::Get<Parameters::NumPressurePointsEquil>()
-        : eclState.getTableManager().getEqldims().getNumDepthNodesP();
+    if (Parameters::IsSet<Parameters::NumPressurePointsEquil>()) {
+        numPressurePointsEquil_ = Parameters::Get<Parameters::NumPressurePointsEquil>();
+        if (numPressurePointsEquil_ < 1) {
+            throw std::invalid_argument {
+                fmt::format("--num-pressure-points-equil must be at least 1, "
+                            "but {} was given.", numPressurePointsEquil_)
+            };
+        }
+    }
+    else {
+        numPressurePointsEquil_ = eclState.getTableManager().getEqldims().getNumDepthNodesP();
+        if (numPressurePointsEquil_ < 1) {
+            throw std::invalid_argument {
+                fmt::format("EQLDIMS item 2, the number of depth nodes in the "
+                            "equilibration pressure tables, must be at least 1, "
+                            "but {} was given.", numPressurePointsEquil_)
+            };
+        }
+    }
 
     explicitRockCompaction_ = Parameters::Get<Parameters::ExplicitRockCompaction>();
 }

opm/simulators/flow/FlowMain.hpp

@@ -29,6 +29,7 @@
 #include <opm/models/utils/start.hh>
 
 #include <opm/simulators/flow/Banners.hpp>
+#include <opm/simulators/flow/ConvergenceOutputConfiguration.hpp>
 #include <opm/simulators/flow/FlowUtils.hpp>
 #include <opm/simulators/flow/NlddReporting.hpp>
 #include <opm/simulators/flow/SimulatorFullyImplicit.hpp>
@@ -435,8 +436,14 @@ namespace Opm {
                 = omp_get_max_threads();
 #endif
 
+            const auto extraConvOutput = ConvergenceOutputConfiguration {
+                Parameters::Get<Parameters::OutputExtraConvergenceInfo>(),
+                R"(OutputExtraConvergenceInfo (--output-extra-convergence-info))"
+            };
+
             printFlowTrailer(mpi_size_, threads, total_setup_time_, deck_read_time_, report,
-                             simulator_->model().simulator().problem().extraTrailerSummary());
+                             simulator_->model().simulator().problem().extraTrailerSummary(),
+                             extraConvOutput.want(ConvergenceOutputConfiguration::Option::Performance));
 
             detail::handleExtraConvergenceOutput(report,
                                                  Parameters::Get<Parameters::OutputExtraConvergenceInfo>(),

opm/simulators/flow/SimulatorFullyImplicit.cpp

@@ -40,7 +40,9 @@ void registerSimulatorParameters()
          "\"none\" gives no extra output and "
          "overrides all other options, "
          "\"steps\" generates an INFOSTEP file, "
-         "\"iterations\" generates an INFOITER file. "
+         "\"iterations\" generates an INFOITER file, "
+         "\"performance\" adds detailed timings to the "
+         "performance summary at the end of the simulation. "
          "Combine options with commas, e.g., "
          "\"steps,iterations\" for multiple outputs.");
     Parameters::Register<Parameters::SaveStep>

opm/simulators/flow/Transmissibility.hpp

@@ -100,6 +100,24 @@ class Transmissibility {
      */
     Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const;
 
+    /*!
+     * \brief One-sided (half) transmissibility of the face between two
+     *        cells, seen from the inside cell (NTG applied, before face
+     *        multipliers). Only available when setStoreHalfTrans(true)
+     *        was called before update(); used by the adjoint module for
+     *        the permeability chain rule
+     *        dT/dK_inside = (h_out/(h_in+h_out))^2 dh_in/dK_inside.
+     */
+    Scalar halfTransmissibility(unsigned insideElemIdx, unsigned outsideElemIdx) const;
+
+    //! \brief Enable storage of the one-sided half transmissibilities
+    //!        (call before update()).
+    void setStoreHalfTrans(bool yesno)
+    { storeHalfTrans_ = yesno; }
+
+    bool storeHalfTrans() const
+    { return storeHalfTrans_; }
+
     Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const;
 
     const std::map<std::pair<unsigned, unsigned>, Scalar>& getThermalHalfTransBoundary() const;
@@ -296,6 +314,8 @@ class Transmissibility {
     bool enableDispersivity_;
     bool warnEditNNC_ = true;
     std::unordered_map<std::uint64_t, Scalar> thermalHalfTrans_; //NB this is based on direction map size is ca 2*trans_ (diffusivity_)
+    std::unordered_map<std::uint64_t, Scalar> halfTrans_; // directional, only filled when storeHalfTrans_
+    bool storeHalfTrans_ = false;
     std::unordered_map<std::uint64_t, Scalar> diffusivity_;
     std::unordered_map<std::uint64_t, Scalar> dispersivity_;
 

opm/simulators/flow/Transmissibility_impl.hpp

@@ -138,6 +138,24 @@ thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
     return thermalHalfTrans_.at(details::directionalIsId(insideElemIdx, outsideElemIdx));
 }
 
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+Scalar Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
+halfTransmissibility(unsigned insideElemIdx, unsigned outsideElemIdx) const
+{
+    if (!storeHalfTrans_) {
+        // Without this the caller gets a bare std::out_of_range from the empty
+        // map.  One branch on a bool, and this is not a hot path: only the
+        // adjoint code asks for these today.
+        OPM_THROW(std::logic_error,
+                  "One-sided half transmissibilities were not stored. "
+                  "Call setStoreHalfTrans(true) before update() to have them "
+                  "computed; they are off by default because only the adjoint "
+                  "code needs them.");
+    }
+
+    return halfTrans_.at(details::directionalIsId(insideElemIdx, outsideElemIdx));
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 Scalar Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
@@ -212,6 +230,13 @@ update(bool global, const TransUpdateQuantities update_quantities,
 
     transBoundary_.clear();
 
+    if (storeHalfTrans_) {
+        halfTrans_.clear();
+        if (num_threads == 1) {
+            halfTrans_.reserve(numElements*6*1.05);
+        }
+    }
+
     // if energy is enabled, let's do the same for the "thermal half transmissibilities"
     if ( enableEnergy_ && !onlyTrans) {
         thermalHalfTrans_.clear();
@@ -312,6 +337,8 @@ update(bool global, const TransUpdateQuantities update_quantities,
                                                   MapBuilderInsertionMode::Insert_Or_Assign);
     ThreadSafeMapBuilder thermalHalfTrans(thermalHalfTrans_, num_threads,
                                           MapBuilderInsertionMode::Insert_Or_Assign);
+    ThreadSafeMapBuilder halfTransMap(halfTrans_, num_threads,
+                                      MapBuilderInsertionMode::Insert_Or_Assign);
     ThreadSafeMapBuilder diffusivity(diffusivity_, num_threads,
                                      MapBuilderInsertionMode::Insert_Or_Assign);
     ThreadSafeMapBuilder dispersivity(dispersivity_, num_threads,
@@ -457,6 +484,22 @@ update(bool global, const TransUpdateQuantities update_quantities,
 
                 Scalar trans = computeHalfMean(computeHalfTrans_, permeability_);
 
+                if (storeHalfTrans_) {
+                    // one-sided half transmissibilities (NTG applied),
+                    // for the adjoint permeability chain rule
+                    auto onesided = computeHalf(computeHalfTrans_,
+                                                permeability_[inside.elemIdx],
+                                                permeability_[outside.elemIdx]);
+                    applyNtg_(onesided[0], inside, ntg);
+                    applyNtg_(onesided[1], outside, ntg);
+                    halfTransMap.insert_or_assign(
+                        details::directionalIsId(inside.elemIdx, outside.elemIdx),
+                        onesided[0]);
+                    halfTransMap.insert_or_assign(
+                        details::directionalIsId(outside.elemIdx, inside.elemIdx),
+                        onesided[1]);
+                }
+
                 // apply the full face transmissibility multipliers
                 // for the inside ...
                 if (!pinchActive) {
@@ -559,6 +602,10 @@ update(bool global, const TransUpdateQuantities update_quantities,
 #pragma omp section
 #endif
         dispersivity.finalize();
+#ifdef _OPENMP
+#pragma omp section
+#endif
+        halfTransMap.finalize();
     }
 
     // Potentially overwrite and/or modify transmissibilities based on input from deck

opm/simulators/flow/equil/InitStateEquil_impl.hpp

@@ -323,6 +323,11 @@ RK4IVP<Scalar,RHS>::RK4IVP(const RHS& f,
     : N_(N)
     , span_(span)
 {
+    // stepsize() divides by N_ and operator() evaluates interval N_ - 1.  A
+    // non-positive sample count is rejected when the value is loaded in
+    // FlowGenericProblem, where its source is known.
+    assert(N >= 1);
+
     const Scalar h = stepsize();
     const Scalar h2 = h / 2;
     const Scalar h6 = h / 6;