Merge pull request #5959 from akva2/add_voigt_symmtensor

Add VoigtArray and SymmTensor classes

CMakeLists_files.cmake

@@ -159,7 +159,9 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/utils/PartiallySupportedFlowKeywords.cpp
   opm/simulators/utils/PressureAverage.cpp
   opm/simulators/utils/SerializationPackers.cpp
+  opm/simulators/utils/SymmTensor.cpp
   opm/simulators/utils/UnsupportedFlowKeywords.cpp
+  opm/simulators/utils/VoigtArray.cpp
   opm/simulators/utils/compressPartition.cpp
   opm/simulators/utils/gatherDeferredLogger.cpp
   opm/simulators/utils/phaseUsageFromDeck.cpp
@@ -417,8 +419,10 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_RestartSerialization.cpp
   tests/test_rstconv.cpp
   tests/test_stoppedwells.cpp
+  tests/test_SymmTensor.cpp
   tests/test_timer.cpp
   tests/test_vfpproperties.cpp
+  tests/test_VoigtArray.cpp
   tests/test_wellmodel.cpp
   tests/test_wellprodindexcalculator.cpp
   tests/test_wellstate.cpp
@@ -985,6 +989,8 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/utils/ParallelSerialization.hpp
   opm/simulators/utils/readDeck.hpp
   opm/simulators/utils/satfunc/RelpermDiagnostics.hpp
+  opm/simulators/utils/SymmTensor.hpp
+  opm/simulators/utils/VoigtArray.hpp
   opm/simulators/wells/ALQState.hpp
   opm/simulators/wells/BlackoilWellModel.hpp
   opm/simulators/wells/BlackoilWellModel_impl.hpp

opm/simulators/utils/SymmTensor.cpp

@@ -0,0 +1,128 @@
+/*
+  Copyright 2025 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/>.
+*/
+#include <config.h>
+#include <opm/simulators/utils/SymmTensor.hpp>
+
+#include <algorithm>
+
+namespace Opm {
+
+template<class T>
+void SymmTensor<T>::
+operator+=(const T data)
+{
+    std::for_each(this->data_.begin(), this->data_.end(),
+                  [&data](auto& v) { v += data; });
+}
+
+template<class T>
+void SymmTensor<T>::
+operator+=(const SymmTensor& data)
+{
+    auto it = data.data_.begin();
+    std::for_each(this->data_.begin(), this->data_.end(),
+                  [&it](auto& v) { v += *it++; });
+}
+
+template<class T>
+void SymmTensor<T>::
+operator*=(const T value)
+{
+    std::for_each(this->data_.begin(), this->data_.end(),
+                  [&value](auto& v) { v *= value; });
+}
+
+template<class T>
+void SymmTensor<T>::reset()
+{
+    this->data_.fill(T{0});
+}
+
+template<class T>
+T SymmTensor<T>::
+trace() const
+{
+    return (*this)[VoigtIndex::XX] +
+           (*this)[VoigtIndex::YY] +
+           (*this)[VoigtIndex::ZZ];
+}
+
+template<class T>
+T SymmTensor<T>::traction(const Dune::FieldVector<T,3>& normal) const
+{
+    T traction = 0.0;
+    constexpr auto& ind = Opm::SymmTensor<double>::diag_indices;
+    for (std::size_t i = 0; i < 3; ++i){
+        traction += (*this)[ind[i]] * normal[i] * normal[i];
+    }
+    traction += T{2} * (*this)[VoigtIndex::YZ] * normal[0] * normal[1]; // xy*nx*ny;
+    traction += T{2} * (*this)[VoigtIndex::XZ] * normal[0] * normal[2]; // xz*nx*nz
+    traction += T{2} * (*this)[VoigtIndex::XY] * normal[1] * normal[2]; // yz*ny*nz
+
+    return traction;
+}
+
+template<class T>
+SymmTensor<T>&
+SymmTensor<T>::operator=(const T value)
+{
+    this->data_.fill(value);
+    return *this;
+}
+
+template<class T1, class T2>
+SymmTensor<T1> operator*(const T2 value, SymmTensor<T1> t1)
+{
+    t1 *= value;
+    return t1;
+}
+
+template<class T1, class T2>
+SymmTensor<T1> operator*(SymmTensor<T1> t1, const T2 value)
+{
+    t1 *= value;
+    return t1;
+}
+
+template<class T>
+SymmTensor<T> operator+(SymmTensor<T> t1, const SymmTensor<T>& t2)
+{
+    t1 += t2;
+    return t1;
+}
+
+#define INSTANTIATE_OPS(T1, T2)                                  \
+    template SymmTensor<T1> operator*(const T2, SymmTensor<T1>); \
+    template SymmTensor<T1> operator*(SymmTensor<T1>, const T2);
+
+#define INSTANTIATE_TYPE(T)                                                \
+    template class SymmTensor<T>;                                          \
+    INSTANTIATE_OPS(T, T)                                                  \
+    INSTANTIATE_OPS(T, int)                                                \
+    INSTANTIATE_OPS(T, unsigned)                                           \
+    INSTANTIATE_OPS(T, std::size_t)                                        \
+    template SymmTensor<T> operator+(SymmTensor<T>, const SymmTensor<T>&);
+
+INSTANTIATE_TYPE(double)
+
+#if FLOW_INSTANTIATE_FLOAT
+INSTANTIATE_TYPE(float)
+#endif
+
+} // namespace Opm

opm/simulators/utils/SymmTensor.hpp

@@ -0,0 +1,64 @@
+/*
+  Copyright 2025 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_UTIL_SYMM_TENSOR_HPP
+#define OPM_UTIL_SYMM_TENSOR_HPP
+
+#include <dune/common/fvector.hh>
+
+#include <opm/simulators/utils/VoigtArray.hpp>
+
+#include <utility>
+
+namespace Opm {
+
+template<class T>
+class SymmTensor : public VoigtContainer<T>
+{
+public:
+    using field_type = T;
+
+    SymmTensor() = default;
+    SymmTensor(std::initializer_list<T> value)
+        : VoigtContainer<T>(std::move(value))
+    {}
+
+    void operator+=(const T data);
+    void operator+=(const SymmTensor<T>& data);
+
+    void operator*=(const T data);
+
+    SymmTensor<T>& operator=(const T value);
+
+    void reset();
+
+    T trace() const;
+    T traction(const Dune::FieldVector<T,3>& normal) const;
+};
+
+template<class T1, class T2>
+SymmTensor<T1> operator*(const T2 value, SymmTensor<T1> t1);
+template<class T1, class T2>
+SymmTensor<T1> operator*(SymmTensor<T1> t1, const T2 value);
+template<class T>
+SymmTensor<T> operator+(SymmTensor<T> t1, const SymmTensor<T>& t2);
+
+} // namespace Opm
+
+#endif // OPM_UTIL_SYMM_TENSOR_HPP

opm/simulators/utils/VoigtArray.cpp

@@ -0,0 +1,85 @@
+/*
+  Copyright 2025 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/>.
+*/
+#include <config.h>
+#include <opm/simulators/utils/VoigtArray.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <algorithm>
+
+namespace Opm {
+
+template<class T>
+template<class Array>
+VoigtContainer<T>::
+VoigtContainer(const Array& array)
+{
+    std::copy(array.begin(), array.end(), data_.begin());
+}
+
+template<class Scalar>
+VoigtArray<Scalar>::
+VoigtArray(const std::size_t size)
+{
+    this->resize(size);
+}
+
+template<class Scalar>
+void VoigtArray<Scalar>::
+resize(const std::size_t size)
+{
+    std::for_each(this->data_.begin(), this->data_.end(),
+                  [size](auto& d) { d.resize(size); });
+}
+
+template<class Scalar>
+void VoigtArray<Scalar>::
+assign(const std::size_t i, const VoigtContainer<Scalar>& array)
+{
+    for (const auto idx : this->unique_indices) {
+        (*this)[idx][i] = array[idx];
+    }
+}
+
+template<class Scalar>
+Scalar VoigtArray<Scalar>::
+operator()(const VoigtIndex idx, const std::size_t i) const
+{
+    return (*this)[idx].at(i);
+}
+
+template<class Scalar>
+Scalar& VoigtArray<Scalar>::
+operator()(const VoigtIndex idx, const std::size_t i)
+{
+    return (*this)[idx].at(i);
+}
+
+#define INSTANTIATE_TYPE(T)                                                    \
+    template class VoigtArray<T>;                                              \
+    template VoigtContainer<T>::VoigtContainer(const std::array<T,6>&);        \
+    template VoigtContainer<T>::VoigtContainer(const Dune::FieldVector<T,6>&);
+
+INSTANTIATE_TYPE(double)
+
+#if FLOW_INSTANTIATE_FLOAT
+INSTANTIATE_TYPE(float)
+#endif
+
+} // namespace Opm

opm/simulators/utils/VoigtArray.hpp

@@ -0,0 +1,110 @@
+/*
+  Copyright 2025 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_UTIL_VOIGT_ARRAY_HPP
+#define OPM_UTIL_VOIGT_ARRAY_HPP
+
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <initializer_list>
+#include <type_traits>
+#include <vector>
+
+namespace Opm {
+
+enum class VoigtIndex {
+    XX =  0, XY =  5, XZ = 4,
+    YX = XY, YY =  1, YZ = 3,
+    ZX = XZ, ZY = YZ, ZZ = 2,
+};
+
+template<class T>
+class VoigtContainer
+{
+public:
+    static constexpr auto indices = std::array{
+        Opm::VoigtIndex::XX,
+        Opm::VoigtIndex::XY,
+        Opm::VoigtIndex::XZ,
+        Opm::VoigtIndex::YX,
+        Opm::VoigtIndex::YY,
+        Opm::VoigtIndex::YZ,
+        Opm::VoigtIndex::ZX,
+        Opm::VoigtIndex::ZY,
+        Opm::VoigtIndex::ZZ,
+    };
+
+    static constexpr auto unique_indices = std::array{
+        Opm::VoigtIndex::XX,
+        Opm::VoigtIndex::YY,
+        Opm::VoigtIndex::ZZ,
+        Opm::VoigtIndex::YZ,
+        Opm::VoigtIndex::XZ,
+        Opm::VoigtIndex::XY
+    };
+
+    static constexpr auto diag_indices = std::array{
+        Opm::VoigtIndex::XX,
+        Opm::VoigtIndex::YY,
+        Opm::VoigtIndex::ZZ,
+    };
+
+    VoigtContainer() = default;
+
+    template<class Array>
+    VoigtContainer(const Array& array);
+
+    VoigtContainer(std::initializer_list<T> value)
+    {
+        std::copy_n(value.begin(),
+                    std::min(data_.size(), value.size()),
+                    data_.begin());
+    }
+
+    const T& operator[](const VoigtIndex idx) const
+    { return data_[static_cast<std::underlying_type_t<VoigtIndex>>(idx)]; }
+
+    T& operator [](const VoigtIndex idx)
+    { return data_[static_cast<std::underlying_type_t<VoigtIndex>>(idx)]; }
+
+    constexpr std::size_t size() const { return data_.size(); }
+
+protected:
+    std::array<T, 6> data_{};
+};
+
+template<class Scalar>
+class VoigtArray : public VoigtContainer<std::vector<Scalar>>
+{
+public:
+    VoigtArray() = default;
+    explicit VoigtArray(const std::size_t size);
+
+    void resize(const std::size_t size);
+
+    Scalar operator()(const VoigtIndex idx, const std::size_t i) const;
+    Scalar& operator()(const VoigtIndex idx, const std::size_t i);
+
+    void assign(const std::size_t i, const VoigtContainer<Scalar>& array);
+};
+
+} // namespace Opm
+
+#endif // OPM_UTIL_VOIGT_ARRAY_HPP

tests/test_SymmTensor.cpp

@@ -0,0 +1,87 @@
+/*
+  Copyright 2025 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/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/utils/SymmTensor.hpp>
+
+#define BOOST_TEST_MODULE SymmTensorTest
+#include <boost/mpl/list.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <stdexcept>
+
+namespace {
+    #if FLOW_INSTANTIATE_FLOAT
+    using Types = boost::mpl::list<float,double>;
+    #else
+    using Types = boost::mpl::list<double>;
+    #endif
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Basics, Scalar, Types)
+{
+    Opm::SymmTensor<Scalar> tensor;
+    constexpr auto& indices = Opm::SymmTensor<Scalar>::unique_indices;
+    for (const auto i : indices) {
+        BOOST_CHECK_EQUAL(tensor[i], Scalar{0});
+    }
+
+    tensor = Scalar{1};
+    for (const auto i : indices) {
+        BOOST_CHECK_EQUAL(tensor[i], Scalar{1});
+    }
+
+    const Opm::SymmTensor<Scalar> tensor2{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::XX], 1.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::YY], 2.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::ZZ], 3.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::YZ], 4.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::YZ], tensor2[Opm::VoigtIndex::ZY]);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::XZ], 5.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::XZ], tensor2[Opm::VoigtIndex::ZX]);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::XY], 6.0);
+    BOOST_CHECK_EQUAL(tensor2[Opm::VoigtIndex::XY], tensor2[Opm::VoigtIndex::YX]);
+
+    const Opm::SymmTensor<Scalar> tensor2s{1.0, 2.0, 3.0};
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::XX], 1.0);
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::YY], 2.0);
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::ZZ], 3.0);
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::YZ], 0.0);
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::XZ], 0.0);
+    BOOST_CHECK_EQUAL(tensor2s[Opm::VoigtIndex::XY], 0.0);
+
+    const auto tensor3 = Scalar{2} * tensor + tensor2;
+    auto tensor4 = tensor;
+    tensor4 *= Scalar{2};
+    tensor += tensor2;
+    constexpr auto& uindices = Opm::SymmTensor<Scalar>::unique_indices;
+    for (const auto i : uindices) {
+        BOOST_CHECK_EQUAL(tensor[i], tensor2[i] + 1.0);
+        BOOST_CHECK_EQUAL(tensor3[i], tensor2[i] + 2.0);
+        BOOST_CHECK_EQUAL(tensor4[i], 2.0);
+    }
+    BOOST_CHECK_EQUAL(tensor.trace(), Scalar{2} + Scalar{3} + Scalar{4});
+    BOOST_CHECK_EQUAL(tensor.traction({1.0, 0.0, 0.0}), 2.0);
+    BOOST_CHECK_EQUAL(tensor.traction({0.0, 1.0, 0.0}), 3.0);
+    BOOST_CHECK_EQUAL(tensor.traction({0.0, 0.0, 1.0}), 4.0);
+    BOOST_CHECK_EQUAL(tensor.traction({1.0, 1.0, 0.0}), 2.0 + 3.0 + 2 * 5.0);
+    BOOST_CHECK_EQUAL(tensor.traction({1.0, 0.0, 1.0}), 2.0 + 4.0 + 2 * 6.0);
+    BOOST_CHECK_EQUAL(tensor.traction({0.0, 1.0, 1.0}), 3.0 + 4.0 + 2 * 7.0);
+}

tests/test_VoigtArray.cpp

@@ -0,0 +1,87 @@
+/*
+  Copyright 2025 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/>.
+*/
+
+#include <config.h>
+
+#include <opm/simulators/utils/VoigtArray.hpp>
+
+#define BOOST_TEST_MODULE VoigArrayTest
+#include <boost/mpl/list.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <stdexcept>
+
+namespace {
+    #if FLOW_INSTANTIATE_FLOAT
+    using Types = boost::mpl::list<float,double>;
+    #else
+    using Types = boost::mpl::list<double>;
+    #endif
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(DefaultConstructed, Scalar, Types)
+{
+    Opm::VoigtArray<Scalar> array;
+    static constexpr auto& indices = Opm::VoigtArray<Scalar>::indices;
+    BOOST_CHECK_EQUAL(array.size(), 6);
+    for (const auto i : indices) {
+        BOOST_CHECK(array[i].empty());
+    }
+
+    for (const auto i : indices) {
+        BOOST_CHECK(array[i].empty());
+        BOOST_CHECK_THROW(array(i, 1), std::out_of_range);
+    }
+
+    array.resize(5);
+    for (const auto i : indices) {
+        BOOST_CHECK(!array[i].empty());
+        BOOST_CHECK_NO_THROW(array(i, 1));
+    }
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(SizeConstructed, Scalar, Types)
+{
+    Opm::VoigtArray<Scalar> array(4);
+    BOOST_CHECK_EQUAL(array.size(), 6);
+    static constexpr auto& indices = Opm::VoigtArray<Scalar>::indices;
+    static constexpr auto& uindices = Opm::VoigtArray<Scalar>::unique_indices;
+    for (const auto i : indices) {
+        BOOST_CHECK_EQUAL(array[i].size(), 4);
+    }
+
+    for (const auto i : uindices) {
+        std::fill(array[i].begin(), array[i].end(), static_cast<Scalar>(i));
+    }
+
+    for (const auto i: indices) {
+        BOOST_CHECK_EQUAL(array[i].size(), 4);
+        BOOST_CHECK_EQUAL(array[i][1], static_cast<Scalar>(i));
+        BOOST_CHECK_EQUAL(array(i,1), static_cast<Scalar>(i));
+    }
+
+    BOOST_CHECK_EQUAL(array(Opm::VoigtIndex::XY, 1),
+                      array(Opm::VoigtIndex::YX, 1));
+
+    BOOST_CHECK_EQUAL(array(Opm::VoigtIndex::XZ, 2),
+                      array(Opm::VoigtIndex::ZX, 2));
+
+    BOOST_CHECK_EQUAL(array(Opm::VoigtIndex::YZ, 3),
+                      array(Opm::VoigtIndex::ZY, 3));
+}