Add Facility for Collecting Connection Level Fracturing Statistics

In particular,

  * RunningStatistics calculates basic statics measures of a running
    sample
  * ConnFracStatistics aggregates those measures for connection
    level fracture pressure, fracture width, and flow rate

CMakeLists_files.cmake

@@ -1014,6 +1014,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/BlackoilWellModelWBP.hpp
   opm/simulators/wells/ConnectionIndexMap.hpp
   opm/simulators/wells/ConnFiltrateData.hpp
+  opm/simulators/wells/ConnFracStatistics.hpp
   opm/simulators/wells/FractionCalculator.hpp
   opm/simulators/wells/GasLiftCommon.hpp
   opm/simulators/wells/GasLiftGroupInfo.hpp
@@ -1045,6 +1046,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/RegionAttributeHelpers.hpp
   opm/simulators/wells/RegionAverageCalculator.hpp
   opm/simulators/wells/RuntimePerforation.hpp
+  opm/simulators/wells/RunningStatistics.hpp
   opm/simulators/wells/SingleWellState.hpp
   opm/simulators/wells/StandardWell.hpp
   opm/simulators/wells/StandardWell_impl.hpp

opm/simulators/wells/ConnFracStatistics.hpp

@@ -0,0 +1,139 @@
+/*
+  Copyright 2024 Equinor ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_CONNFRACSTATISTICS_HPP
+#define OPM_CONNFRACSTATISTICS_HPP
+
+#include <opm/simulators/wells/RunningStatistics.hpp>
+
+#include <array>
+#include <cstddef>
+#include <type_traits>
+
+namespace Opm {
+
+/// Collection of fracturing statistics measures at the connection level.
+///
+/// \tparam Scalar Statistics element type.  Typically a built-in arithmetic
+/// type like \c float or \c double.
+template <typename Scalar>
+class ConnFracStatistics
+{
+public:
+    /// Known quantities for which this collection provides statistics
+    /// measures.
+    enum class Quantity : std::size_t
+    {
+        /// Fracture pressure
+        Pressure,
+
+        /// Fracture flow rate
+        FlowRate,
+
+        /// Fracture width
+        Width,
+
+        // -------------------------------------------------------------
+        // Helper.  Must be last enumerator.
+        NumQuantities,
+    };
+
+    /// Sample point representation.
+    ///
+    /// Client code must populate an object of this type in order to collect
+    /// statistics.
+    using SamplePoint = std::array
+        <Scalar, static_cast<std::underlying_type_t<Quantity>>(Quantity::NumQuantities)>;
+
+    /// Convert between byte array and object representation.
+    ///
+    /// \tparam Serializer Byte array conversion protocol.
+    ///
+    /// \param[in,out] serializer Byte array conversion object.
+    template <class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(this->quantity_);
+    }
+
+    /// Reset internal counters to prepare for calculating a new set of
+    /// sample statistics.
+    void reset()
+    {
+        for (auto& q : this->quantity_) { q.reset(); }
+    }
+
+    /// Include new element into sample.
+    ///
+    /// Updates internal statistics counters.
+    ///
+    /// \param[in] samplePoint Collection of sample values for the
+    /// fracturing of the current well/reservoir connection.
+    void addSamplePoint(const SamplePoint& samplePoint)
+    {
+        for (auto qIdx = 0*samplePoint.size(); qIdx < samplePoint.size(); ++qIdx) {
+            this->quantity_[qIdx].addSamplePoint(samplePoint[qIdx]);
+        }
+    }
+
+    /// Retrieve collection of sample statistics for a single quantity.
+    ///
+    /// \param[in] q Quantity for which to retrieve sample statistics.
+    ///
+    /// \return Sample statistics for quantity \p q.
+    const RunningStatistics<Scalar>& statistics(const Quantity q) const
+    {
+        return this->quantity_[ static_cast<std::underlying_type_t<Quantity>>(q) ];
+    }
+
+    /// Create a serialisation test object.
+    static ConnFracStatistics serializationTestObject()
+    {
+        auto stat = ConnFracStatistics{};
+
+        stat.quantity_
+            .fill(RunningStatistics<Scalar>::serializationTestObject());
+
+        return stat;
+    }
+
+    /// Equality predicate.
+    ///
+    /// \param[in] that Object against which \code *this \endcode will
+    /// be tested for equality.
+    ///
+    /// \return Whether or not \code *this \endcode is the same as \p that.
+    bool operator==(const ConnFracStatistics& that) const
+    {
+        return this->quantity_ == that.quantity_;
+    }
+
+private:
+    using StatArray = std::array<
+    RunningStatistics<Scalar>,
+    static_cast<std::underlying_type_t<Quantity>>(Quantity::NumQuantities)
+    >;
+
+    /// Collection of connection level fracturing statistics.
+    StatArray quantity_{};
+};
+
+} // namespace Opm
+
+#endif // OPM_CONNFRACSTATISTICS_HPP

opm/simulators/wells/RunningStatistics.hpp

@@ -0,0 +1,163 @@
+/*
+  Copyright 2024 Equinor ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_RUNNING_STATISTICS_HPP
+#define OPM_RUNNING_STATISTICS_HPP
+
+#include <cmath>
+#include <limits>
+#include <optional>
+
+namespace Opm {
+
+/// Facility for calculating simple sample statistics without having full
+/// sample available.
+///
+/// \tparam Scalar Sample element type.  Typically a built-in arithmetic
+/// type like \c float or \c double.
+template <typename Scalar>
+class RunningStatistics
+{
+public:
+    /// Convert between byte array and object representation.
+    ///
+    /// \tparam Serializer Byte array conversion protocol.
+    ///
+    /// \param[in,out] serializer Byte array conversion object.
+    template <class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(this->sampleSize_);
+        serializer(this->min_);
+        serializer(this->max_);
+        serializer(this->mean_);
+        serializer(this->totalVariance_);
+    }
+
+    /// Create a serialisation test object.
+    static RunningStatistics serializationTestObject()
+    {
+        auto stat = RunningStatistics{};
+
+        stat.sampleSize_ = 12;
+        stat.min_ = -static_cast<Scalar>(1);
+        stat.max_ =  static_cast<Scalar>(2);
+        stat.mean_ = static_cast<Scalar>(0.03);
+        stat.totalVariance_ = static_cast<Scalar>(0.4);
+
+        return stat;
+    }
+
+    /// Equality predicate.
+    ///
+    /// \param[in] that Object against which \code *this \endcode will
+    /// be tested for equality.
+    ///
+    /// \return Whether or not \code *this \endcode is the same as \p that.
+    bool operator==(const RunningStatistics& that) const
+    {
+        return (this->sampleSize_ == that.sampleSize_)
+            && (this->min_ == that.min_)
+            && (this->max_ == that.max_)
+            && (this->mean_ == that.mean_)
+            && (this->totalVariance_ == that.totalVariance_)
+            ;
+    }
+
+    /// Reset internal counters to prepare for calculating a new set of
+    /// sample statistics.
+    void reset()
+    {
+        this->sampleSize_ = 0;
+        this->min_ = std::numeric_limits<Scalar>::max();
+        this->max_ = std::numeric_limits<Scalar>::lowest();
+        this->mean_ = Scalar{};
+        this->totalVariance_ = Scalar{};
+    }
+
+    /// Include new element into sample.
+    ///
+    /// Updates internal statistics counters.
+    ///
+    /// \param[in] x Sample point.
+    void addSamplePoint(const Scalar x)
+    {
+        if (x < this->min_) { this->min_ = x; }
+        if (x > this->max_) { this->max_ = x; }
+
+        // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
+
+        ++this->sampleSize_;
+
+        const auto d1 = x - this->mean();
+        this->mean_ += d1 / this->sampleSize_;
+
+        const auto d2 = x - this->mean();
+        this->totalVariance_ += d1 * d2;
+    }
+
+    /// Retrieve current sample size.
+    ///
+    /// Effectively returns the number of calls to addSamplePoint() since
+    /// object was constructed or since the previous call to reset().
+    std::size_t sampleSize() const { return this->sampleSize_; }
+
+    /// Retrieve smallest sample value seen so far.
+    Scalar min() const { return this->min_; }
+
+    /// Retrieve largest sample value seen so far.
+    Scalar max() const { return this->max_; }
+
+    /// Retrieve arithmetic average of all sample points seen so far.
+    Scalar mean() const { return this->mean_; }
+
+    /// Retrieve unbiased standard deviation of all sample points seen so
+    /// far.
+    ///
+    /// Returns nullopt if number of sample points is less than two.
+    std::optional<Scalar> stdev() const
+    {
+        if (this->sampleSize_ < 2) {
+            return {};
+        }
+
+        using std::sqrt;
+        return sqrt(this->totalVariance_ / (this->sampleSize_ - 1));
+    }
+
+private:
+    /// Current sample size.
+    std::size_t sampleSize_{};
+
+    /// Smallest sample value seen so far.
+    Scalar min_ { std::numeric_limits<Scalar>::max() };
+
+    /// Largest sample value seen so far.
+    Scalar max_ { std::numeric_limits<Scalar>::lowest() };
+
+    /// Arithmetic average of all sample points seen so far.
+    Scalar mean_{};
+
+    /// Variance measure.  In particular, N-1 * Var{x_i}.
+    Scalar totalVariance_{};
+};
+
+} // namespace Opm
+
+#endif // OPM_RUNNING_STATISTICS_HPP