Stress tensor output to restart and VTK files

CMakeLists_files.cmake

@@ -496,6 +496,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_interregflows.cpp
   tests/test_invert.cpp
   tests/test_keyword_validator.cpp
+  tests/test_linearleastsquares.cpp
   tests/test_LogOutputHelper.cpp
   tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
@@ -962,6 +963,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/models/utils/basicparameters.hh
   opm/models/utils/basicproperties.hh
   opm/models/utils/genericguard.hh
+  opm/models/utils/linearleastsquares.hpp
   opm/models/utils/parametersystem.hpp
   opm/models/utils/pffgridvector.hh
   opm/models/utils/prefetch.hh

opm/models/discretization/common/fvbaseprimaryvariables.hh

@@ -183,11 +183,11 @@ namespace Dune {
     };
 
     /** Specialization of FieldTraits for all PrimaryVariables derived from Opm::FvBasePrimaryVariables */
-    template<class TypeTag, template <class> class EwomsPrimaryVariable>
-    struct FieldTraits<EwomsPrimaryVariable<TypeTag>>
-      : public FieldTraitsImpl<TypeTag,
-                               std::is_base_of_v<Opm::FvBasePrimaryVariables<TypeTag>,
-                                                 EwomsPrimaryVariable<TypeTag>>>
+    template <class TypeTag, template <class> class EwomsPrimaryVariable>
+    requires std::is_base_of_v<Opm::FvBasePrimaryVariables<TypeTag>,
+                               EwomsPrimaryVariable<TypeTag> >
+    struct FieldTraits<EwomsPrimaryVariable<TypeTag> >
+        : public FieldTraitsImpl<TypeTag, true>
     {
     };
 

opm/models/discretization/common/tpsalinearizer.hpp

@@ -72,6 +72,8 @@ class TpsaLinearizer
     using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
     using VectorBlock = Dune::FieldVector<Scalar, numEq>;
 
+    using StressInfoVector = Dune::FieldVector<Scalar, 3>;
+
 public:
     // ///
     // Public functions
@@ -301,6 +303,16 @@ public:
     const LinearizationType& getLinearizationType() const
     { return linearizationType_; }
 
+    /*!
+     * Container with information to compute stress output
+     *
+     * \returns Stress info. container
+     */
+    const auto& getStressInfo() const
+    {
+        return stressInfo_;
+    }
+
     /*!
     * \brief Get constraints map
     *
@@ -347,6 +359,11 @@ private:
         unsigned numCells = flowModel.numTotalDof();
         neighborInfo_.reserve(numCells, 6 * numCells);
         std::vector<NeighborInfo> loc_nbinfo;
+
+        stressInfo_.reserve(numCells, 6 * numCells);
+        std::vector<StressInfo> loc_stressinfo;
+        StressInfoVector nullTraction(0.0);
+        StressInfoVector nullFaceNormal(0.0);
         for (const auto& elem : elements(gridView_())) {
             // Loop over primary dofs in the element
             stencil.update(elem);
@@ -356,6 +373,10 @@ private:
                 const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
                 loc_nbinfo.resize(stencil.numDof() - 1);
 
+                // Build information for stress output
+                const int numFaces = stencil.numBoundaryFaces() + stencil.numInteriorFaces();
+                loc_stressinfo.resize(numFaces);
+
                 for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
                     // NOTE: NeighborInfo could/should be expanded with cell face parameters located in problem_()
                     // needed when computing face terms in LocalResidual
@@ -366,6 +387,10 @@ private:
                         const auto& scvf = stencil.interiorFace(scvfIdx);
                         const Scalar area = scvf.area();
                         loc_nbinfo[dofIdx - 1] = NeighborInfo{ neighborIdx, area, nullptr };
+
+                        int faceId = scvf.dirId();
+                        loc_stressinfo[dofIdx - 1] =
+                            StressInfo{faceId, nullTraction, nullFaceNormal, area};
                     }
                 }
 
@@ -379,6 +404,13 @@ private:
                         // Get boundary face direction
                         const auto& bf = stencil.boundaryFace(bfIndex);
                         const int dir_id = bf.dirId();
+                        const auto bfArea = bf.area();
+
+                        // Initialize boundary information for stress output container
+                        // OBS: NNCs are not implemented with TPSA, hence traction and face
+                        // normal vector will stay zero for NNC connections!
+                        loc_stressinfo[stencil.numInteriorFaces() + bfIndex] =
+                            StressInfo{dir_id, nullTraction, nullFaceNormal, bfArea};
 
                         // Skip NNCs
                         if (dir_id < 0) {
@@ -395,7 +427,7 @@ private:
                         }
 
                         // Insert boundary condition data in container
-                        BoundaryConditionData bcdata { type, displacement, bfIndex, bf.area() };
+                        BoundaryConditionData bcdata { type, displacement, bfIndex, bfArea };
                         boundaryInfo_.push_back( { myIdx, dir_id, bfIndex, bcdata } );
                         ++bfIndex;
                         continue;
@@ -405,6 +437,7 @@ private:
                         continue;
                     }
                 }
+                stressInfo_.appendRow(loc_stressinfo.begin(), loc_stressinfo.end());
             }
         }
 
@@ -462,6 +495,7 @@ private:
                 OPM_TIMEBLOCK_LOCAL(faceCalculationForEachCellTPSA, Subsystem::Assembly);
 
                 // Loop over neighboring cells
+                short loc = 0;
                 for (const auto& nbInfo : nbInfos) {
                     OPM_TIMEBLOCK_LOCAL(calculationForEachFaceTPSA, Subsystem::Assembly);
 
@@ -501,6 +535,16 @@ private:
                     // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
                     bMat *= -1.0;
                     *nbInfo.matBlockAddress += bMat;
+
+                    // Insert interior traction vector
+                    // OBS: Assume traction vector is the three first entries in residual!
+                    const auto& faceNormal = problem.cellFaceNormal(globI, globJ);
+                    stressInfo_[globI][loc].faceNormal = faceNormal;
+                    stressInfo_[globI][loc].faceArea = nbInfo.faceArea;
+                    for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
+                        stressInfo_[globI][loc].traction[tractionIdx] = res[tractionIdx];
+                    }
+                    ++loc;
                 }
             }
 
@@ -582,6 +626,17 @@ private:
             // Insert contribution to (globI, globI) sub-block
             // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
             *diagMatAddress_[globI] += bMat;
+
+            // Insert boundary traction vector
+            // OBS: Assume traction vector is the three first entries in residual!
+            const auto& nbInfos = neighborInfo_[globI];
+            short loc = nbInfos.size() + bdyInfo.bfIndex;
+            const auto& bndyNormal = problem.cellFaceNormalBoundary(globI, bdyInfo.bfIndex);
+            stressInfo_[globI][loc].faceNormal = bndyNormal;
+            stressInfo_[globI][loc].faceArea = bdyInfo.bcdata.faceArea;
+            for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
+                stressInfo_[globI][loc].traction[tractionIdx] = res[tractionIdx];
+            }
         }
     }
 
@@ -706,6 +761,19 @@ private:
     };
     SparseTable<NeighborInfo> neighborInfo_{};
 
+    /*!
+     * \brief Face stress information
+    */
+    struct StressInfo
+    {
+        int faceId;
+        StressInfoVector traction;
+        StressInfoVector faceNormal;
+        double faceArea;
+    };
+
+    SparseTable<StressInfo> stressInfo_{};
+
     /*!
     * \brief Information for boundary conditions
     */

opm/models/io/vtktpsamodule.hpp

@@ -25,6 +25,9 @@
 #ifndef VTK_TPSA_MODULE_HPP
 #define VTK_TPSA_MODULE_HPP
 
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/fvector.hh>
+
 #include <opm/material/densead/Math.hpp>
 
 #include <opm/models/io/baseoutputmodule.hh>
@@ -34,7 +37,6 @@
 #include <opm/models/utils/parametersystem.hpp>
 #include <opm/models/utils/propertysystem.hh>
 
-
 namespace Opm {
 
 /*!
@@ -50,6 +52,7 @@ class VtkTpsaModule : public BaseOutputModule<TypeTag>
     using GridView = GetPropType<TypeTag, Properties::GridView>;
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
     static constexpr auto vtkFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
     using VtkMultiWriter = Opm::VtkMultiWriter<GridView, vtkFormat>;
@@ -59,6 +62,10 @@ class VtkTpsaModule : public BaseOutputModule<TypeTag>
     using BufferType = typename ParentType::BufferType;
     using ScalarBuffer = typename ParentType::ScalarBuffer;
     using VectorBuffer = typename ParentType::VectorBuffer;
+    using TensorBuffer = typename ParentType::TensorBuffer;
+
+    using SymTensor = Dune::FieldVector<Scalar, 6>;
+    using Tensor = Dune::DynamicMatrix<Scalar>;
 
 public:
     /*!
@@ -111,6 +118,11 @@ public:
             if (params_.solidPressureOutput_) {
                 this->resizeScalarBuffer_(solidPres_, BufferType::Dof);
             }
+
+            // Stress
+            if (params_.stressOutput_) {
+                this->resizeTensorBuffer_(stress_, BufferType::Dof);
+            }
         }
     }
 
@@ -130,6 +142,7 @@ public:
             // Assign quantities from material state
             const auto& problem = elemCtx.problem();
             const auto& geoMechModel = problem.geoMechModel();
+            const auto& faceStressInfo = geoMechModel.linearizer().getStressInfo();
             for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
                 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
                 const auto& materialState = geoMechModel.materialState(globalDofIdx, /*timeIdx=*/0);
@@ -151,6 +164,13 @@ public:
                 if (params_.solidPressureOutput_) {
                     solidPres_[globalDofIdx] = scalarValue(materialState.solidPressure());
                 }
+
+                // Stress
+                if (params_.stressOutput_ && !faceStressInfo.empty()) {
+                    SymTensor stressSymTensor = geoMechModel.stress(globalDofIdx, false);
+                    Tensor& stressTensor = stress_[globalDofIdx];
+                    setTensorFromVoigt_(stressTensor, stressSymTensor);
+                }
             }
         }
     }
@@ -182,17 +202,42 @@ public:
             if (params_.solidPressureOutput_) {
                 this->commitScalarBuffer_(baseWriter, "solid_pressure", solidPres_, BufferType::Dof);
             }
+
+            if (params_.stressOutput_) {
+                this->commitTensorBuffer_(baseWriter, "stress", stress_, BufferType::Dof);
+            }
         }
     }
 
 private:
+    static void setTensorFromVoigt_(Tensor& tensor, const SymTensor& symTensor)
+    {
+        // Diagonal terms
+        for (std::size_t i = 0; i < 3; ++i) {
+            tensor[i][i] = symTensor[i];
+        }
+
+        // Off-diagonal terms
+        tensor[0][1] = symTensor[5];
+        tensor[0][2] = symTensor[4];
+        tensor[1][2] = symTensor[3];
+        for (std::size_t i = 0; i < 3; ++i) {
+            for (std::size_t j = 0; j < 3; ++j) {
+                if (i > j) {
+                    tensor[i][j] = tensor[j][i];
+                }
+            }
+        }
+    }
+
     VtkTpsaParams params_{};
 
     VectorBuffer displacement_{};
     VectorBuffer rotation_{};
     ScalarBuffer solidPres_{};
-};  // class VtkTpsaModule
+    TensorBuffer stress_{};
+}; // class VtkTpsaModule
 
-}  // namespace Opm
+} // namespace Opm
 
 #endif

opm/models/io/vtktpsaparams.cpp

@@ -39,6 +39,8 @@ void VtkTpsaParams::registerParameters()
         ("Include rotation in VTK output files");
     Parameters::Register<Parameters::VtkWriteSolidPressure>
         ("Include solid pressure in VTK output files");
+    Parameters::Register<Parameters::VtkWriteStress>
+        ("Include stress tensor in VTK output files");
 }
 
 /*!
@@ -49,6 +51,7 @@ void VtkTpsaParams::read()
     displacementOutput_ = Parameters::Get<Parameters::VtkWriteDisplacement>();
     rotationOutput_ = Parameters::Get<Parameters::VtkWriteRotation>();
     solidPressureOutput_ = Parameters::Get<Parameters::VtkWriteSolidPressure>();
+    stressOutput_ = Parameters::Get<Parameters::VtkWriteStress>();
 }
 
 }  // namespace Opm

opm/models/io/vtktpsaparams.hpp

@@ -32,6 +32,7 @@ namespace Opm::Parameters {
 struct VtkWriteDisplacement { static constexpr bool value = true; };
 struct VtkWriteRotation { static constexpr bool value = true; };
 struct VtkWriteSolidPressure { static constexpr bool value = true; };
+struct VtkWriteStress { static constexpr bool value = true; };
 
 }  // namespace Opm::Parameters
 
@@ -48,6 +49,7 @@ struct VtkTpsaParams
     bool displacementOutput_;
     bool rotationOutput_;
     bool solidPressureOutput_;
+    bool stressOutput_;
 };
 
 }  // namespace Opm

opm/models/tpsa/tpsamodel.hpp

@@ -27,16 +27,21 @@
 
 #include <dune/grid/common/gridenums.hh>
 
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
 #include <opm/grid/utility/ElementChunks.hpp>
 
 #include <opm/material/common/MathToolbox.hpp>
 
 #include <opm/models/parallel/threadmanager.hpp>
 #include <opm/models/tpsa/tpsabaseproperties.hpp>
+#include <opm/models/utils/linearleastsquares.hpp>
 #include <opm/models/utils/propertysystem.hh>
 
 #include <array>
 #include <memory>
+#include <unordered_map>
 
 
 namespace Opm {
@@ -466,12 +471,101 @@ public:
     * \param with_fracture Boolean to activate fracture output
     * \returns Stress tensor (Voigt notation) at grid cell
     *
-    * \note Needed in OutputBlackOilModule, but zero for now!
+    * \warning NNC are not implemented with TPSA, therefore cells with NNC faces will give
+    * incorrect stress output!
     */
-    SymTensor stress(const unsigned /*globalIdx*/, const bool /*with_fracture*/) const
+    SymTensor stress(const unsigned globalIdx, const bool /*with_fracture*/) const
     {
-        SymTensor val;
-        return val;
+        SymTensor stressOutput;
+        const auto& stressInfo = linearizer_->getStressInfo();
+        if (!stressInfo.empty()) {
+            const auto& stressInfoGlobI = stressInfo[globalIdx];
+
+            // Setup least-squares matrix of face normals and right-hand side of traction vectors
+            // per faces. Matrix shape = 3 rows per face x 6 columns for each (symmetric) stress
+            // tensor component.
+            std::unordered_multimap<int, std::size_t> faceIdMap;
+            std::size_t mapSize = 0;
+            for (std::size_t sInfoIdx = 0; sInfoIdx < stressInfoGlobI.size(); ++sInfoIdx) {
+                const auto faceId  = stressInfoGlobI[sInfoIdx].faceId;
+
+                // OBS: NNC faces are not added to least squares system, since TPSA does not handle
+                // NNC connections, like numerical aquifers, yet!
+                if (faceId < 0 || stressInfoGlobI[sInfoIdx].faceArea == 0.0) {
+                    continue;
+                }
+                if (!faceIdMap.contains(faceId)) {
+                    ++mapSize;
+                }
+                faceIdMap.insert({faceId, sInfoIdx});
+            }
+
+            // If no valid faces exist for cell, return zero SymTensor
+            if (mapSize == 0) {
+                return stressOutput;
+            }
+
+            Dune::DynamicMatrix<Scalar> mat(3 * mapSize, 6);
+            Dune::DynamicVector<Scalar> rhs(3 * mapSize);
+            std::size_t rowIdx = 0;
+            for (auto it = faceIdMap.begin(); it != faceIdMap.end();) {
+                auto [first, last] = faceIdMap.equal_range(it->first);
+
+                // Loop over possible multiple of the same face
+                Scalar sumFaceArea = 0.0;
+                for (auto inner = first; inner != last; ++inner) {
+                    const auto sInfoIdx = inner->second;
+                    const auto& faceStressInfo = stressInfoGlobI[sInfoIdx];
+
+                    // One normal per face
+                    if (inner == first) {
+                        const auto& fNormal = faceStressInfo.faceNormal;
+                        mat[3*rowIdx][0] = fNormal[0];
+                        mat[3*rowIdx][3] = fNormal[1];
+                        mat[3*rowIdx][4] = fNormal[2];
+
+                        mat[3*rowIdx + 1][1] = fNormal[1];
+                        mat[3*rowIdx + 1][3] = fNormal[0];
+                        mat[3*rowIdx + 1][5] = fNormal[2];
+
+                        mat[3*rowIdx + 2][2] = fNormal[2];
+                        mat[3*rowIdx + 2][4] = fNormal[0];
+                        mat[3*rowIdx + 2][5] = fNormal[1];
+                    }
+
+                    // Add traction vectors for same faces
+                    const auto& fTraction = faceStressInfo.traction;
+                    rhs[3*rowIdx] += fTraction[0];
+                    rhs[3*rowIdx + 1] += fTraction[1];
+                    rhs[3*rowIdx + 2] += fTraction[2];
+
+                    // Sum face area
+                    sumFaceArea += faceStressInfo.faceArea;
+                }
+                // Divide rhs for (unique) face by sum of face areas
+                rhs[3*rowIdx] /= sumFaceArea;
+                rhs[3*rowIdx + 1] /= sumFaceArea;
+                rhs[3*rowIdx + 2] /= sumFaceArea;
+
+                // Move to next unique face
+                ++rowIdx;
+                it = last;
+            }
+
+            // Use least squares solver for underdetermined system
+            LinearLeastSquares<Scalar> lsq(mat, rhs);
+            lsq.solve();
+
+            // Reconstructed stress tensor at cell center
+            const auto& stress = lsq.x();
+            stressOutput[0] = stress[0]; // XX
+            stressOutput[1] = stress[1]; // YY
+            stressOutput[2] = stress[2]; // ZZ
+            stressOutput[3] = stress[5]; // YZ
+            stressOutput[4] = stress[4]; // XZ
+            stressOutput[5] = stress[3]; // XY
+        }
+        return stressOutput;
     }
 
     /*!

opm/models/utils/linearleastsquares.hpp

@@ -0,0 +1,188 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+#ifndef LINEAR_LEAST_SQUARES_HPP
+#define LINEAR_LEAST_SQUARES_HPP
+
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
+
+#include <stdexcept>
+
+namespace Opm
+{
+
+/*!
+ * @brief Linear least squares calculations and properties
+ *
+ * @warning Do not use with large matrices since matrix inversion is LU decomposition
+ */
+template <class Scalar>
+class LinearLeastSquares
+{
+    using Matrix = Dune::DynamicMatrix<Scalar>;
+    using Vector = Dune::DynamicVector<Scalar>;
+
+public:
+    /*!
+     * @brief Constructor
+     *
+     * @param A Coefficient matrix
+     * @param b Right-hand side (data) vector
+     */
+    LinearLeastSquares(const Matrix& A, const Vector& b)
+        : A_(A), b_(b), x_(A.M())
+    {
+    }
+
+    /*!
+     * @brief Solve linear least squares system
+     *
+     * Done by solving the normal equations: (A^T*A)*x = A^T*b
+     */
+    void solve()
+    {
+        solveNormalEquations_();
+    }
+
+    /*!
+     * @brief Read-only vector of calculated coefficient vector
+     *
+     * @return Coeff. vector
+     */
+    const Vector& x() const
+    {
+        return x_;
+    }
+
+    /*!
+     * @brief Evaluate regression model at input x vector
+     *
+     * @param x Input vector
+     * @return Regression model value
+     */
+    Scalar evaluate(const Vector& x) const
+    {
+        return x_ * x;
+    }
+
+    /*!
+     * @brief Measure of the discrepancy between data and regression model
+     *
+     * @return Sum of the square of residual
+     */
+    Scalar residualSumOfSquares() const
+    {
+        Vector r(b_);
+        A_.mmv(x_, r);
+        return r.two_norm2();
+    }
+
+    /*!
+     * @brief Value for how well regression model represents the model
+     *
+     * @return Model sum of squares
+     */
+    Scalar explainedSumOfSquares() const
+    {
+        Scalar ymean = 0.0;
+        const auto ndata = b_.N();
+        for (size_t i = 0; i < ndata; ++i) {
+            ymean += b_[i];
+        }
+        ymean /= ndata;
+
+        Vector r(ndata, ymean);
+        A_.mmv(x_, r);
+
+        return r.two_norm2();
+    }
+
+    /*!
+     * @brief Sum of all squared differences
+     *
+     * Total sum of squares = explained sum of squares + residual sum of squares
+     *
+     * @return Total sum of squares
+     */
+    Scalar totalSumOfSquares() const
+    {
+        Scalar ymean = 0.0;
+        const auto ndata = b_.N();
+        for (size_t i = 0; i < ndata; ++i) {
+            ymean += b_[i];
+        }
+        ymean /= ndata;
+
+        Vector r(ndata, ymean);
+        r -= b_;
+
+        return r.two_norm2();
+    }
+
+    /*!
+     * @brief Coefficient of determination
+     *
+     * Typical value for how well a regression model fits the data
+     *
+     * @return R^2 value
+     */
+    Scalar RSquared() const
+    {
+        const auto tss = totalSumOfSquares();
+        if (tss < 1e-16) {
+            OPM_THROW(std::runtime_error,
+                      "Total sum of squares is close to zero, hence R^2 is undefined!");
+        }
+        return explainedSumOfSquares() / tss;
+    }
+
+private:
+    /*!
+     * @brief Solve the linear least squares system
+     */
+    void solveNormalEquations_()
+    {
+        // Right-hand side, bhat = A^T*b
+        Vector bhat(A_.M());
+        A_.mtv(b_, bhat);
+
+        // Normal matrix Ahat = A^T*A
+        Matrix Ahat(A_.M(), A_.M());
+        detail::multMatrixTransposed(A_, A_, Ahat);
+
+        // Solve normal equations x = Ahat^{-1}*bhat
+        Ahat.solve(x_, bhat, /*doPivoting=*/true);
+    }
+
+    Matrix A_{}; ///< Input matrix
+    Vector b_{}; ///< Data vector
+    Vector x_{}; ///< Coefficient vector
+};
+
+} // Opm
+
+#endif // LINEAR_LEAST_SQUARES_HPP

opm/simulators/flow/MechContainer.cpp

@@ -167,7 +167,7 @@ outputRestart(data::Solution& sol)
         DataEntry{"MECHPOTF", UnitSystem::measure::pressure, &potentialForce_},
         DataEntry{"PRESPOTF", UnitSystem::measure::pressure, &potentialPressForce_},
         DataEntry{"STRAIN",   UnitSystem::measure::identity, &strain_},
-        DataEntry{"STRESS",   UnitSystem::measure::length,   &stress_},
+        DataEntry{"STRESS",   UnitSystem::measure::pressure, &stress_},
         DataEntry{"TEMPPOTF", UnitSystem::measure::pressure, &potentialTempForce_},
     };
 

tests/test_linearleastsquares.cpp

@@ -0,0 +1,127 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  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 2 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/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+#include "config.h"
+
+#define BOOST_TEST_MODULE TestLinearLeastSquares
+
+#include <boost/mpl/list.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <dune/common/dynmatrix.hh>
+#include <dune/common/dynvector.hh>
+
+#include <opm/models/utils/linearleastsquares.hpp>
+
+using namespace Opm;
+
+using ScalarTypes = boost::mpl::list<double, float>;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(SimpleOverdeterminedSystem, Scalar, ScalarTypes)
+{
+    // Solve an overdetermined system Ax=b using linear least squares (i.e. normal equations)
+    // Example from: https://en.wikipedia.org/wiki/Linear_least_squares
+    Dune::DynamicMatrix<Scalar> A(3, 2);
+    A[0][0] = 1.0;
+    A[0][1] = 0.0;
+    A[1][0] = 0.0;
+    A[1][1] = 1.0;
+    A[2][0] = 1.0;
+    A[2][1] = 1.0;
+
+    Dune::DynamicVector<Scalar> b(3);
+    b[0] = 1.0;
+    b[1] = 1.0;
+    b[2] = 0.0;
+
+    // Instantiate LinearLeastSquares and solve normal equations
+    LinearLeastSquares<Scalar> lsq(A, b);
+    lsq.solve();
+
+    // Check results
+    Dune::DynamicVector<Scalar> xExpected(2);
+    xExpected[0] = 1.0 / 3.0;
+    xExpected[1] = 1.0 / 3.0;
+
+    const auto& x = lsq.x();
+    BOOST_CHECK_EQUAL(x.N(), xExpected.N());
+    for (std::size_t i = 0; i < x.N(); ++i) {
+        BOOST_CHECK_CLOSE(x[i], xExpected[i], 1e-4);
+    }
+
+    // Check sum of squared residuals
+    BOOST_CHECK_CLOSE(lsq.residualSumOfSquares(), 4.0 / 3.0, 1e-4);
+
+    // Check total sum of squares
+    BOOST_CHECK_CLOSE(lsq.totalSumOfSquares(), 2.0 / 3.0, 1e-4);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(LinearRegression, Scalar, ScalarTypes)
+{
+    // Linear regression for line: y = beta_0 + beta_1*x, with some y-data
+    // Example from: https://en.wikipedia.org/wiki/Linear_least_squares#Example
+    Dune::DynamicMatrix<Scalar> X(4, 2);
+    X[0][0] = 1.0;
+    X[0][1] = 1.0;
+    X[1][0] = 1.0;
+    X[1][1] = 2.0;
+    X[2][0] = 1.0;
+    X[2][1] = 3.0;
+    X[3][0] = 1.0;
+    X[3][1] = 4.0;
+
+    Dune::DynamicVector<Scalar> y(4);
+    y[0] = 6.0;
+    y[1] = 5.0;
+    y[2] = 7.0;
+    y[3] = 10.0;
+
+    // Solve linear least squares
+    LinearLeastSquares<Scalar> lsq(X, y);
+    lsq.solve();
+
+    // Check against solution
+    Dune::DynamicVector<Scalar> betaExpected(2);
+    betaExpected[0] = 3.5;
+    betaExpected[1] = 1.4;
+
+    const auto& beta = lsq.x();
+    BOOST_CHECK_EQUAL(beta.N(), betaExpected.N());
+    for (std::size_t i = 0; i < beta.N(); ++i) {
+        BOOST_CHECK_CLOSE(beta[i], betaExpected[i], 1e-4);
+    }
+
+    // Check sum of squared residuals
+    BOOST_CHECK_CLOSE(lsq.residualSumOfSquares(), 4.2, 1e-4);
+
+    // Check total sum of squares
+    BOOST_CHECK_CLOSE(lsq.totalSumOfSquares(), 14.0, 1e-4);
+
+    // Check R-squared, goodness-of-fit
+    BOOST_CHECK_CLOSE(lsq.RSquared(), 0.7, 1e-4);
+
+    // Check interpolation, i.e. y = beta^T * x for an arbitrary x
+    Dune::DynamicVector<Scalar> xPred(2);
+    xPred[0] = 1.0;
+    xPred[1] = 5.0;
+    BOOST_CHECK_CLOSE(lsq.evaluate(xPred), 10.5, 1e-4);
+}