Collect Statistics on Connection Level Fracturing Process

This commit adds a new member function

    Fracture::assignGeomechWellState(ConnFracStatistics<Scalar>&)

which computes statistical measures (minimum value, maximum value,
average value/mean, and standard deviation) of fracture pressures,
injection flow rate, and fracture width for a single fracture.  We
also add a new member function

    FractureModel::assignGeomechWellState(WellState<Scalar>& wellState)

which calls the above for each fracture in the model.  Collectively,
these enable emitting the new summary vectors

    CFR<Q>{MAX,MIN,AVG,STD}

for the quantity <Q> as 'P' (pressure), 'IR' (injection rate), and
'WD' (width) through the regular mechanism.

opm/geomech/Fracture.cpp

@@ -11,6 +11,7 @@
 
 #include <opm/simulators/linalg/FlowLinearSolverParameters.hpp>
 #include <opm/simulators/linalg/setupPropertyTree.hpp>
+#include <opm/simulators/wells/ConnFracStatistics.hpp>
 
 #include <opm/geomech/DiscreteDisplacement.hpp>
 #include <opm/geomech/GridStretcher.hpp>
@@ -959,6 +960,34 @@ std::vector<std::tuple<int,double,double>> Fracture::wellIndices() const
     return wellIndices;
 }
 
+template <typename Scalar>
+void Fracture::assignGeomechWellState(ConnFracStatistics<Scalar>& stats) const
+{
+    using Quantity = typename ConnFracStatistics<Scalar>::Quantity;
+
+    constexpr auto pressIx = static_cast<std::underlying_type_t<Quantity>>(Quantity::Pressure);
+    constexpr auto rateIx  = static_cast<std::underlying_type_t<Quantity>>(Quantity::FlowRate);
+    constexpr auto widthIx = static_cast<std::underlying_type_t<Quantity>>(Quantity::Width);
+
+    const auto nCells = this->reservoir_cells_.size();
+
+    stats.reset();
+
+    for (auto cellIx = 0*nCells; cellIx < nCells; ++cellIx) {
+        auto samplePoint = typename ConnFracStatistics<Scalar>::SamplePoint{};
+
+        samplePoint[pressIx] = this->fracture_pressure_[cellIx][0];
+
+        samplePoint[rateIx] = this->leakof_[cellIx]
+            * (this->fracture_pressure_[cellIx][0] -
+               this->reservoir_pressure_[cellIx]);
+
+        samplePoint[widthIx] = this->fracture_width_[cellIx][0];
+
+        stats.addSamplePoint(samplePoint);
+    }
+}
+
 void Fracture::writePressureSystem() const
 {
     if(prm_.get<bool>("write_pressure_system")){
@@ -1255,4 +1284,7 @@ template void
 Fracture::updateReservoirCells(const external::cvf::ref<external::cvf::BoundingBoxTree>& cellSearchTree,
                                const Dune::PolyhedralGrid<3,3,double>& grid3D);
 
+template void Fracture::assignGeomechWellState(ConnFracStatistics<float>&) const;
+template void Fracture::assignGeomechWellState(ConnFracStatistics<double>&) const;
+
 } // namespace Opm

opm/geomech/Fracture.hpp

@@ -53,6 +53,11 @@
 
 #include <opm/geomech/GridStretcher.hpp>
 
+namespace Opm {
+    template <typename Scalar>
+    class ConnFracStatistics;
+}
+
 namespace Opm
 {
 struct WellInfo {
@@ -184,6 +189,9 @@ class Fracture
     Dune::FieldVector<double, 6> strain(const Dune::FieldVector<double, 3>& obs) const;
     Dune::FieldVector<double, 3> disp(const Dune::FieldVector<double, 3>& obs) const;
 
+    template <typename Scalar>
+    void assignGeomechWellState(ConnFracStatistics<Scalar>& stats) const;
+
 private:
      Dune::BlockVector<Dune::FieldVector<double, 3>> all_slips() const;
     void resetWriters();

opm/geomech/FractureModel.cpp

@@ -7,6 +7,9 @@
 
 #include <opm/simulators/linalg/PropertyTree.hpp>
 
+#include <opm/simulators/wells/SingleWellState.hpp>
+#include <opm/simulators/wells/WellState.hpp>
+
 #include <algorithm>
 #include <cassert>
 #include <iostream>
@@ -254,4 +257,39 @@ namespace Opm{
 
         return wellindices;
     }
+
+    template <typename Scalar>
+    void FractureModel::assignGeomechWellState(WellState<Scalar>& wellState) const
+    {
+        const auto nWells = this->wells_.size();
+        for (auto i = 0*nWells; i < nWells; ++i) {
+            const auto wsIx = wellState.index(this->wells_[i].name());
+            if (! wsIx.has_value()) { continue; }
+
+            auto& perfData = wellState[*wsIx].perf_data;
+
+            if (perfData.connFracStatistics.size() != perfData.cell_index.size()) {
+                perfData.connFracStatistics.resize(perfData.cell_index.size());
+            }
+
+            for (const auto& fracture : this->well_fractures_[i]) {
+                auto perfPos = std::find(perfData.cell_index.begin(),
+                                         perfData.cell_index.end(),
+                                         fracture.wellInfo().well_cell);
+                if (perfPos == perfData.cell_index.end()) { continue; }
+
+                // Possibly just "fracture.wellInfo().perf" instead.
+                const auto perfIx = std::distance(perfData.cell_index.begin(), perfPos);
+
+                fracture.assignGeomechWellState(perfData.connFracStatistics[perfIx]);
+            }
+        }
+    }
 }
+
+// ===========================================================================
+// Explicit specialisations.  No other code below separator.
+// ===========================================================================
+
+template void Opm::FractureModel::assignGeomechWellState(WellState<float>&) const;
+template void Opm::FractureModel::assignGeomechWellState(WellState<double>&) const;

opm/geomech/FractureModel.hpp

@@ -31,6 +31,11 @@ std::vector<T> as_vector(ptree const& pt, ptree::key_type const& key)
     return r;
 }
 
+namespace Opm {
+    template <typename Scalar>
+    class WellState;
+}
+
 // template <typename T>
 // std::vector<T> opm_as_vector(const Opm::PropertyTree& pt, const std::string& key)
 // {
@@ -143,6 +148,9 @@ class FractureModel
     std::vector<std::tuple<int,double,double>>
     getExtraWellIndices(const std::string& wellname) const;
 
+    template <typename Scalar>
+    void assignGeomechWellState(WellState<Scalar>& wellState) const;
+
     bool addPertsToSchedule(){return prm_.get<bool>("addperfs_to_schedule");}
 
     // probably this should be collected in one loop since all do full loop over fracture ++ well

opm/geomech/eclproblemgeomech.hh

@@ -263,6 +263,9 @@ namespace Opm{
                         assert(false);
                         this->addConnectionsToWell();
                     }
+
+                    this->geomechModel_.fractureModel()
+                        .assignGeomechWellState(this->wellModel_.wellState());
                 }
             }