Move invert3x3 out of general purpose utils.hxx header

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

include/bout/utils.hxx

@@ -410,59 +410,6 @@ bool operator==(const Tensor<T>& lhs, const Tensor<T>& 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 <typename T>
-int invert3x3(Matrix<T>& 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
  */

src/mesh/coordinates.cxx

@@ -18,6 +18,7 @@
 
 #include <bout/globals.hxx>
 
+#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;
     }
 

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 <http://www.gnu.org/licenses/>.
+ *
+ **************************************************************************/
+
+#pragma once
+
+#include <bout/utils.hxx>
+#include <optional>
+
+/// 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<BoutReal> invert3x3(Matrix<BoutReal>& 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<BoutReal>{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

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

tests/unit/mesh/test_invert3x3.cxx

@@ -0,0 +1,74 @@
+#include "../../src/mesh/invert3x3.hxx"
+
+#include "gtest/gtest.h"
+
+TEST(Invert3x3Test, Identity) {
+  Matrix<BoutReal> 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<BoutReal> rawDataMat = {0.05567105, 0.92458227, 0.19954631,
+                                      0.28581972, 0.54009039, 0.13234403,
+                                      0.8841194,  0.161224,   0.74853209};
+  std::vector<BoutReal> rawDataInv = {-2.48021781, 4.27410022,  -0.09449605,
+                                      0.6278449,   0.87275842,  -0.32168092,
+                                      2.79424897,  -5.23628123, 1.51684677};
+
+  Matrix<BoutReal> input(3, 3);
+  Matrix<BoutReal> 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<BoutReal> input(3, 3);
+  input = 0;
+  auto result = bout::invert3x3(input);
+  EXPECT_TRUE(result.has_value());
+}
+
+TEST(Invert3x3Test, BadCondition) {
+  Matrix<BoutReal> 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());
+}

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<BoutReal> 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<BoutReal> rawDataMat = {0.05567105, 0.92458227, 0.19954631,
-                                      0.28581972, 0.54009039, 0.13234403,
-                                      0.8841194,  0.161224,   0.74853209};
-  std::vector<BoutReal> rawDataInv = {-2.48021781, 4.27410022,  -0.09449605,
-                                      0.6278449,   0.87275842,  -0.32168092,
-                                      2.79424897,  -5.23628123, 1.51684677};
-
-  Matrix<BoutReal> input(3, 3);
-  Matrix<BoutReal> 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<BoutReal> input(3, 3);
-  input = 0;
-  EXPECT_THROW(invert3x3(input), BoutException);
-}
-
-TEST(Invert3x3Test, BadCondition) {
-  Matrix<BoutReal> 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));