diff --git a/CMakeLists.txt b/CMakeLists.txt index 5e1924ba1e..8b3b5256a6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -270,6 +270,7 @@ set(BOUT_SOURCES ./src/mesh/interpolation/interpolation_z.cxx ./src/mesh/interpolation/lagrange_4pt_xz.cxx ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx + ./src/mesh/invert3x3.hxx ./src/mesh/mesh.cxx ./src/mesh/parallel/fci.cxx ./src/mesh/parallel/fci.hxx diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx index b45152fbcc..d0281da5fc 100644 --- a/include/bout/utils.hxx +++ b/include/bout/utils.hxx @@ -410,59 +410,6 @@ bool operator==(const Tensor& lhs, const Tensor& rhs) { return std::equal(lhs.begin(), lhs.end(), rhs.begin()); } -/************************************************************************** - * Matrix routines - **************************************************************************/ -/// Explicit inversion of a 3x3 matrix \p a -/// -/// The input \p small determines how small the determinant must be for -/// us to throw due to the matrix being singular (ill conditioned); -/// If small is less than zero then instead of throwing we return 1. -/// This is ugly but can be used to support some use cases. -template -int invert3x3(Matrix& a, BoutReal small = 1.0e-15) { - TRACE("invert3x3"); - - // Calculate the first co-factors - T A = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1); - T B = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2); - T C = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0); - - // Calculate the determinant - T det = a(0, 0) * A + a(0, 1) * B + a(0, 2) * C; - - if (std::abs(det) < std::abs(small)) { - if (small >= 0) { - throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small); - } else { - return 1; - } - } - - // Calculate the rest of the co-factors - T D = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2); - T E = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0); - T F = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1); - T G = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1); - T H = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2); - T I = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0); - - // Now construct the output, overwrites input - T detinv = 1.0 / det; - - a(0, 0) = A * detinv; - a(0, 1) = D * detinv; - a(0, 2) = G * detinv; - a(1, 0) = B * detinv; - a(1, 1) = E * detinv; - a(1, 2) = H * detinv; - a(2, 0) = C * detinv; - a(2, 1) = F * detinv; - a(2, 2) = I * detinv; - - return 0; -} - /*! * Get Random number between 0 and 1 */ diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 4e515449ca..ce096cc8f9 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -18,6 +18,7 @@ #include +#include "invert3x3.hxx" #include "parallel/fci.hxx" #include "parallel/shiftedmetricinterp.hxx" @@ -1241,9 +1242,10 @@ int Coordinates::calcCovariant(const std::string& region) { a(1, 2) = a(2, 1) = g23[i]; a(0, 2) = a(2, 0) = g13[i]; - if (invert3x3(a)) { - output_error.write("\tERROR: metric tensor is singular at ({:d}, {:d})\n", i.x(), - i.y()); + if (const auto det = bout::invert3x3(a); det.has_value()) { + output_error.write( + "\tERROR: metric tensor is singular at ({:d}, {:d}), determinant: {:d}\n", + i.x(), i.y(), det.value()); return 1; } @@ -1297,9 +1299,10 @@ int Coordinates::calcContravariant(const std::string& region) { a(1, 2) = a(2, 1) = g_23[i]; a(0, 2) = a(2, 0) = g_13[i]; - if (invert3x3(a)) { - output_error.write("\tERROR: metric tensor is singular at ({:d}, {:d})\n", i.x(), - i.y()); + if (const auto det = bout::invert3x3(a); det.has_value()) { + output_error.write( + "\tERROR: metric tensor is singular at ({:d}, {:d}), determinant: {:d}\n", + i.x(), i.y(), det.value()); return 1; } diff --git a/src/mesh/invert3x3.hxx b/src/mesh/invert3x3.hxx new file mode 100644 index 0000000000..9e635d8150 --- /dev/null +++ b/src/mesh/invert3x3.hxx @@ -0,0 +1,77 @@ +/*!************************************************************************* + * \file invert3x3.hxx + * + * A mix of short utilities for memory management, strings, and some + * simple but common calculations + * + ************************************************************************** + * Copyright 2010-2024 B.D.Dudson, BOUT++ Team + * + * Contact: Ben Dudson, dudson2@llnl.gov + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#pragma once + +#include +#include + +/// Explicit inversion of a 3x3 matrix \p a +/// +/// If the matrix is singular (ill conditioned), the determinant is +/// return. Otherwise, an empty `std::optional` is return +namespace bout { +inline std::optional invert3x3(Matrix& a) { + TRACE("invert3x3"); + + // Calculate the first co-factors + BoutReal A = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1); + BoutReal B = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2); + BoutReal C = a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0); + + // Calculate the determinant + const BoutReal det = a(0, 0) * A + a(0, 1) * B + a(0, 2) * C; + constexpr BoutReal small = 1.0e-15; + if (std::abs(det) < std::abs(small)) { + return std::optional{det}; + } + + // Calculate the rest of the co-factors + BoutReal D = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2); + BoutReal E = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0); + BoutReal F = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1); + BoutReal G = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1); + BoutReal H = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2); + BoutReal I = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0); + + // Now construct the output, overwrites input + BoutReal detinv = 1.0 / det; + + a(0, 0) = A * detinv; + a(0, 1) = D * detinv; + a(0, 2) = G * detinv; + a(1, 0) = B * detinv; + a(1, 1) = E * detinv; + a(1, 2) = H * detinv; + a(2, 0) = C * detinv; + a(2, 1) = F * detinv; + a(2, 2) = I * detinv; + + return std::nullopt; +} +} // namespace bout diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index c4ffa4fa75..9d7a83e373 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -69,6 +69,7 @@ set(serial_tests_source ./mesh/test_coordinates.cxx ./mesh/test_coordinates_accessor.cxx ./mesh/test_interpolation.cxx + ./mesh/test_invert3x3.cxx ./mesh/test_mesh.cxx ./mesh/test_paralleltransform.cxx ./solver/test_fakesolver.cxx diff --git a/tests/unit/mesh/test_invert3x3.cxx b/tests/unit/mesh/test_invert3x3.cxx new file mode 100644 index 0000000000..3bc4ae69d8 --- /dev/null +++ b/tests/unit/mesh/test_invert3x3.cxx @@ -0,0 +1,74 @@ +#include "../../src/mesh/invert3x3.hxx" + +#include "gtest/gtest.h" + +TEST(Invert3x3Test, Identity) { + Matrix input(3, 3); + input = 0; + for (int i = 0; i < 3; i++) { + input(i, i) = 1.0; + } + auto expected = input; + bout::invert3x3(input); + + for (int j = 0; j < 3; j++) { + for (int i = 0; i < 3; i++) { + EXPECT_EQ(input(i, j), expected(i, j)); + } + } +} + +TEST(Invert3x3Test, InvertTwice) { + std::vector rawDataMat = {0.05567105, 0.92458227, 0.19954631, + 0.28581972, 0.54009039, 0.13234403, + 0.8841194, 0.161224, 0.74853209}; + std::vector rawDataInv = {-2.48021781, 4.27410022, -0.09449605, + 0.6278449, 0.87275842, -0.32168092, + 2.79424897, -5.23628123, 1.51684677}; + + Matrix input(3, 3); + Matrix expected(3, 3); + + int counter = 0; + for (int j = 0; j < 3; j++) { + for (int i = 0; i < 3; i++) { + input(i, j) = rawDataMat[counter]; + expected(i, j) = rawDataInv[counter]; + counter++; + } + } + + // Invert twice to check if we get back to where we started + bout::invert3x3(input); + + for (int j = 0; j < 3; j++) { + for (int i = 0; i < 3; i++) { + // Note we only check to single tolerance here + EXPECT_FLOAT_EQ(input(i, j), expected(i, j)); + } + } +} + +TEST(Invert3x3Test, Singular) { + Matrix input(3, 3); + input = 0; + auto result = bout::invert3x3(input); + EXPECT_TRUE(result.has_value()); +} + +TEST(Invert3x3Test, BadCondition) { + Matrix input(3, 3); + + input = 0.; + input(0, 0) = 1.0e-16; + input(1, 1) = 1.0; + input(2, 2) = 1.0; + EXPECT_TRUE(bout::invert3x3(input).has_value()); + + // not quite bad enough condition + input = 0.; + input(0, 0) = 1.0e-12; + input(1, 1) = 1.0; + input(2, 2) = 1.0; + EXPECT_FALSE(bout::invert3x3(input).has_value()); +} diff --git a/tests/unit/sys/test_utils.cxx b/tests/unit/sys/test_utils.cxx index 747257bafc..6d84813c48 100644 --- a/tests/unit/sys/test_utils.cxx +++ b/tests/unit/sys/test_utils.cxx @@ -386,91 +386,6 @@ TEST(TensorTest, ConstGetData) { std::all_of(std::begin(tensor), std::end(tensor), [](int a) { return a == 3; })); } -TEST(Invert3x3Test, Identity) { - Matrix input(3, 3); - input = 0; - for (int i = 0; i < 3; i++) { - input(i, i) = 1.0; - } - auto expected = input; - invert3x3(input); - - for (int j = 0; j < 3; j++) { - for (int i = 0; i < 3; i++) { - EXPECT_EQ(input(i, j), expected(i, j)); - } - } -} - -TEST(Invert3x3Test, InvertTwice) { - std::vector rawDataMat = {0.05567105, 0.92458227, 0.19954631, - 0.28581972, 0.54009039, 0.13234403, - 0.8841194, 0.161224, 0.74853209}; - std::vector rawDataInv = {-2.48021781, 4.27410022, -0.09449605, - 0.6278449, 0.87275842, -0.32168092, - 2.79424897, -5.23628123, 1.51684677}; - - Matrix input(3, 3); - Matrix expected(3, 3); - - int counter = 0; - for (int j = 0; j < 3; j++) { - for (int i = 0; i < 3; i++) { - input(i, j) = rawDataMat[counter]; - expected(i, j) = rawDataInv[counter]; - counter++; - } - } - - // Invert twice to check if we get back to where we started - invert3x3(input); - - for (int j = 0; j < 3; j++) { - for (int i = 0; i < 3; i++) { - // Note we only check to single tolerance here - EXPECT_FLOAT_EQ(input(i, j), expected(i, j)); - } - } -} - -TEST(Invert3x3Test, Singular) { - Matrix input(3, 3); - input = 0; - EXPECT_THROW(invert3x3(input), BoutException); -} - -TEST(Invert3x3Test, BadCondition) { - Matrix input(3, 3); - - // Default small - input = 0.; - input(0, 0) = 1.0e-16; - input(1, 1) = 1.0; - input(2, 2) = 1.0; - EXPECT_THROW(invert3x3(input), BoutException); - - // Default small -- not quite bad enough condition - input = 0.; - input(0, 0) = 1.0e-12; - input(1, 1) = 1.0; - input(2, 2) = 1.0; - EXPECT_NO_THROW(invert3x3(input)); - - // Non-default small - input = 0.; - input(0, 0) = 1.0e-12; - input(1, 1) = 1.0; - input(2, 2) = 1.0; - EXPECT_THROW(invert3x3(input, 1.0e-10), BoutException); - - // Non-default small - input = 0.; - input(0, 0) = 1.0e-12; - input(1, 1) = 1.0; - input(2, 2) = 1.0; - EXPECT_NO_THROW(invert3x3(input, -1.0e-10)); -} - TEST(NumberUtilitiesTest, SquareInt) { EXPECT_EQ(4, SQ(2)); EXPECT_EQ(4, SQ(-2));