From 5a3373f93f9a180ab59f2c857b3962ec3028130a Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Feb 2025 17:05:25 +0100 Subject: [PATCH 1/2] added: voigt array this is a container for vector data stored using voigt notation --- CMakeLists_files.cmake | 3 + opm/simulators/utils/VoigtArray.cpp | 85 +++++++++++++++++++++ opm/simulators/utils/VoigtArray.hpp | 112 ++++++++++++++++++++++++++++ tests/test_VoigtArray.cpp | 87 +++++++++++++++++++++ 4 files changed, 287 insertions(+) create mode 100644 opm/simulators/utils/VoigtArray.cpp create mode 100644 opm/simulators/utils/VoigtArray.hpp create mode 100644 tests/test_VoigtArray.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f366b92b1..3b7779d59 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -160,6 +160,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/utils/PressureAverage.cpp opm/simulators/utils/SerializationPackers.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 @@ -419,6 +420,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_stoppedwells.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 +987,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/ParallelSerialization.hpp opm/simulators/utils/readDeck.hpp opm/simulators/utils/satfunc/RelpermDiagnostics.hpp + opm/simulators/utils/VoigtArray.hpp opm/simulators/wells/ALQState.hpp opm/simulators/wells/BlackoilWellModel.hpp opm/simulators/wells/BlackoilWellModel_impl.hpp diff --git a/opm/simulators/utils/VoigtArray.cpp b/opm/simulators/utils/VoigtArray.cpp new file mode 100644 index 000000000..d930c555f --- /dev/null +++ b/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 . +*/ +#include +#include + +#include + +#include + +namespace Opm { + +template +template +VoigtContainer:: +VoigtContainer(const Array& array) +{ + std::copy(array.begin(), array.end(), data_.begin()); +} + +template +VoigtArray:: +VoigtArray(const std::size_t size) +{ + this->resize(size); +} + +template +void VoigtArray:: +resize(const std::size_t size) +{ + std::for_each(this->data_.begin(), this->data_.end(), + [size](auto& d) { d.resize(size); }); +} + +template +void VoigtArray:: +assign(const std::size_t i, const VoigtContainer& array) +{ + for (const auto idx : this->unique_indices) { + (*this)[idx][i] = array[idx]; + } +} + +template +Scalar VoigtArray:: +operator()(const VoigtIndex idx, const std::size_t i) const +{ + return (*this)[idx].at(i); +} + +template +Scalar& VoigtArray:: +operator()(const VoigtIndex idx, const std::size_t i) +{ + return (*this)[idx].at(i); +} + +#define INSTANTIATE_TYPE(T) \ + template class VoigtArray; \ + template VoigtContainer::VoigtContainer(const std::array&); \ + template VoigtContainer::VoigtContainer(const Dune::FieldVector&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/simulators/utils/VoigtArray.hpp b/opm/simulators/utils/VoigtArray.hpp new file mode 100644 index 000000000..b25912c76 --- /dev/null +++ b/opm/simulators/utils/VoigtArray.hpp @@ -0,0 +1,112 @@ +/* + 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 . +*/ + +#ifndef OPM_UTIL_VOIGT_ARRAY_HPP +#define OPM_UTIL_VOIGT_ARRAY_HPP + +#include +#include +#include +#include +#include +#include + +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 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 + VoigtContainer(const Array& array); + + VoigtContainer(std::initializer_list 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>(idx)]; } + + T& operator [](const VoigtIndex idx) + { return data_[static_cast>(idx)]; } + + constexpr std::size_t size() const { return data_.size(); } + +protected: + std::array data_{}; +}; + +template using SymmTensor = VoigtContainer; + +template +class VoigtArray : public VoigtContainer> +{ +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& array); +}; + +} // namespace Opm + +#endif // OPM_UTIL_VOIGT_ARRAY_HPP diff --git a/tests/test_VoigtArray.cpp b/tests/test_VoigtArray.cpp new file mode 100644 index 000000000..de524f9a0 --- /dev/null +++ b/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 . +*/ + +#include + +#include + +#define BOOST_TEST_MODULE VoigArrayTest +#include +#include + +#include + +namespace { + #if FLOW_INSTANTIATE_FLOAT + using Types = boost::mpl::list; + #else + using Types = boost::mpl::list; + #endif +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(DefaultConstructed, Scalar, Types) +{ + Opm::VoigtArray array; + static constexpr auto& indices = Opm::VoigtArray::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 array(4); + BOOST_CHECK_EQUAL(array.size(), 6); + static constexpr auto& indices = Opm::VoigtArray::indices; + static constexpr auto& uindices = Opm::VoigtArray::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(i)); + } + + for (const auto i: indices) { + BOOST_CHECK_EQUAL(array[i].size(), 4); + BOOST_CHECK_EQUAL(array[i][1], static_cast(i)); + BOOST_CHECK_EQUAL(array(i,1), static_cast(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)); +} From 41361baa443c3ca325b80b1dcf45a48d24737fa5 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 6 Feb 2025 09:43:00 +0100 Subject: [PATCH 2/2] added: SymmTensor this is a VoigtContainer with some algebraic operations --- CMakeLists_files.cmake | 3 + opm/simulators/utils/SymmTensor.cpp | 128 ++++++++++++++++++++++++++++ opm/simulators/utils/SymmTensor.hpp | 64 ++++++++++++++ opm/simulators/utils/VoigtArray.hpp | 2 - tests/test_SymmTensor.cpp | 87 +++++++++++++++++++ 5 files changed, 282 insertions(+), 2 deletions(-) create mode 100644 opm/simulators/utils/SymmTensor.cpp create mode 100644 opm/simulators/utils/SymmTensor.hpp create mode 100644 tests/test_SymmTensor.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 3b7779d59..6e6265b91 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -159,6 +159,7 @@ 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 @@ -418,6 +419,7 @@ 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 @@ -987,6 +989,7 @@ 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 diff --git a/opm/simulators/utils/SymmTensor.cpp b/opm/simulators/utils/SymmTensor.cpp new file mode 100644 index 000000000..6b0aef8fa --- /dev/null +++ b/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 . +*/ +#include +#include + +#include + +namespace Opm { + +template +void SymmTensor:: +operator+=(const T data) +{ + std::for_each(this->data_.begin(), this->data_.end(), + [&data](auto& v) { v += data; }); +} + +template +void SymmTensor:: +operator+=(const SymmTensor& data) +{ + auto it = data.data_.begin(); + std::for_each(this->data_.begin(), this->data_.end(), + [&it](auto& v) { v += *it++; }); +} + +template +void SymmTensor:: +operator*=(const T value) +{ + std::for_each(this->data_.begin(), this->data_.end(), + [&value](auto& v) { v *= value; }); +} + +template +void SymmTensor::reset() +{ + this->data_.fill(T{0}); +} + +template +T SymmTensor:: +trace() const +{ + return (*this)[VoigtIndex::XX] + + (*this)[VoigtIndex::YY] + + (*this)[VoigtIndex::ZZ]; +} + +template +T SymmTensor::traction(const Dune::FieldVector& normal) const +{ + T traction = 0.0; + constexpr auto& ind = Opm::SymmTensor::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 +SymmTensor& +SymmTensor::operator=(const T value) +{ + this->data_.fill(value); + return *this; +} + +template +SymmTensor operator*(const T2 value, SymmTensor t1) +{ + t1 *= value; + return t1; +} + +template +SymmTensor operator*(SymmTensor t1, const T2 value) +{ + t1 *= value; + return t1; +} + +template +SymmTensor operator+(SymmTensor t1, const SymmTensor& t2) +{ + t1 += t2; + return t1; +} + +#define INSTANTIATE_OPS(T1, T2) \ + template SymmTensor operator*(const T2, SymmTensor); \ + template SymmTensor operator*(SymmTensor, const T2); + +#define INSTANTIATE_TYPE(T) \ + template class SymmTensor; \ + INSTANTIATE_OPS(T, T) \ + INSTANTIATE_OPS(T, int) \ + INSTANTIATE_OPS(T, unsigned) \ + INSTANTIATE_OPS(T, std::size_t) \ + template SymmTensor operator+(SymmTensor, const SymmTensor&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/simulators/utils/SymmTensor.hpp b/opm/simulators/utils/SymmTensor.hpp new file mode 100644 index 000000000..50b0260ae --- /dev/null +++ b/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 . +*/ + +#ifndef OPM_UTIL_SYMM_TENSOR_HPP +#define OPM_UTIL_SYMM_TENSOR_HPP + +#include + +#include + +#include + +namespace Opm { + +template +class SymmTensor : public VoigtContainer +{ +public: + using field_type = T; + + SymmTensor() = default; + SymmTensor(std::initializer_list value) + : VoigtContainer(std::move(value)) + {} + + void operator+=(const T data); + void operator+=(const SymmTensor& data); + + void operator*=(const T data); + + SymmTensor& operator=(const T value); + + void reset(); + + T trace() const; + T traction(const Dune::FieldVector& normal) const; +}; + +template +SymmTensor operator*(const T2 value, SymmTensor t1); +template +SymmTensor operator*(SymmTensor t1, const T2 value); +template +SymmTensor operator+(SymmTensor t1, const SymmTensor& t2); + +} // namespace Opm + +#endif // OPM_UTIL_SYMM_TENSOR_HPP diff --git a/opm/simulators/utils/VoigtArray.hpp b/opm/simulators/utils/VoigtArray.hpp index b25912c76..e9917d57b 100644 --- a/opm/simulators/utils/VoigtArray.hpp +++ b/opm/simulators/utils/VoigtArray.hpp @@ -90,8 +90,6 @@ protected: std::array data_{}; }; -template using SymmTensor = VoigtContainer; - template class VoigtArray : public VoigtContainer> { diff --git a/tests/test_SymmTensor.cpp b/tests/test_SymmTensor.cpp new file mode 100644 index 000000000..1b2d4a6dd --- /dev/null +++ b/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 . +*/ + +#include + +#include + +#define BOOST_TEST_MODULE SymmTensorTest +#include +#include + +#include + +namespace { + #if FLOW_INSTANTIATE_FLOAT + using Types = boost::mpl::list; + #else + using Types = boost::mpl::list; + #endif +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(Basics, Scalar, Types) +{ + Opm::SymmTensor tensor; + constexpr auto& indices = Opm::SymmTensor::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 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 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::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); +}