Return std::optional from invert3x3

Allows throwing more specific error in coordinates

Author: ZedThree

src/mesh/coordinates.cxx

@@ -1264,9 +1264,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;
     }
 
@@ -1320,9 +1321,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

@@ -29,42 +29,38 @@
 #pragma once
 
 #include <bout/utils.hxx>
+#include <optional>
 
 /// 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 false.
-/// This is ugly but can be used to support some use cases.
-template <typename T>
-bool invert3x3(Matrix<T>& a, T small = 1.0e-15) {
+/// 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
-  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);
+  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
-  T det = a(0, 0) * A + a(0, 1) * B + a(0, 2) * C;
-
+  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)) {
-    if (small >= 0) {
-      throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small);
-    }
-    return false;
+    return std::optional<BoutReal>{det};
   }
 
   // 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);
+  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
-  T detinv = 1.0 / det;
+  BoutReal detinv = 1.0 / det;
 
   a(0, 0) = A * detinv;
   a(0, 1) = D * detinv;
@@ -76,5 +72,6 @@ bool invert3x3(Matrix<T>& a, T small = 1.0e-15) {
   a(2, 1) = F * detinv;
   a(2, 2) = I * detinv;
 
-  return true;
+  return std::nullopt;
 }
+} // namespace bout

tests/unit/mesh/test_invert3x3.cxx

@@ -9,7 +9,7 @@ TEST(Invert3x3Test, Identity) {
     input(i, i) = 1.0;
   }
   auto expected = input;
-  invert3x3(input);
+  bout::invert3x3(input);
 
   for (int j = 0; j < 3; j++) {
     for (int i = 0; i < 3; i++) {
@@ -39,7 +39,7 @@ TEST(Invert3x3Test, InvertTwice) {
   }
 
   // Invert twice to check if we get back to where we started
-  invert3x3(input);
+  bout::invert3x3(input);
 
   for (int j = 0; j < 3; j++) {
     for (int i = 0; i < 3; i++) {
@@ -52,37 +52,23 @@ TEST(Invert3x3Test, InvertTwice) {
 TEST(Invert3x3Test, Singular) {
   Matrix<BoutReal> input(3, 3);
   input = 0;
-  EXPECT_THROW(invert3x3(input), BoutException);
+  auto result = bout::invert3x3(input);
+  EXPECT_TRUE(result.has_value());
 }
 
 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);
+  EXPECT_TRUE(bout::invert3x3(input).has_value());
 
-  // Default small -- not quite bad enough condition
+  // 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));
+  EXPECT_FALSE(bout::invert3x3(input).has_value());
 }