From e594a890306a547d5ffa8297ff9c7a67c1fb9638 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 30 Oct 2024 13:58:31 +0000 Subject: [PATCH 01/10] Move `invert3x3` out of general purpose `utils.hxx` header Only used in `Coordinates`, so make private implementation detail --- CMakeLists.txt | 1 + include/bout/utils.hxx | 53 ------------------ src/mesh/coordinates.cxx | 1 + src/mesh/invert3x3.hxx | 81 +++++++++++++++++++++++++++ tests/unit/CMakeLists.txt | 1 + tests/unit/mesh/test_invert3x3.cxx | 89 ++++++++++++++++++++++++++++++ tests/unit/sys/test_utils.cxx | 85 ---------------------------- 7 files changed, 173 insertions(+), 138 deletions(-) create mode 100644 src/mesh/invert3x3.hxx create mode 100644 tests/unit/mesh/test_invert3x3.cxx 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..e636c8cd1a 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -20,6 +20,7 @@ #include "parallel/fci.hxx" #include "parallel/shiftedmetricinterp.hxx" +#include "invert3x3.hxx" // use anonymous namespace so this utility function is not available outside this file namespace { diff --git a/src/mesh/invert3x3.hxx b/src/mesh/invert3x3.hxx new file mode 100644 index 0000000000..dce208338d --- /dev/null +++ b/src/mesh/invert3x3.hxx @@ -0,0 +1,81 @@ +/*!************************************************************************* + * \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 + +/// 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; +} 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..02beeec644 --- /dev/null +++ b/tests/unit/mesh/test_invert3x3.cxx @@ -0,0 +1,89 @@ +#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; + 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)); +} + 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)); From 6e839d6514fe9ade6662a683dbe22c469a275912 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Wed, 30 Oct 2024 14:20:14 +0000 Subject: [PATCH 02/10] Return `bool` instead of `int` from `invert3x3` --- src/mesh/coordinates.cxx | 4 ++-- src/mesh/invert3x3.hxx | 9 ++++----- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index e636c8cd1a..b57ad055da 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1242,7 +1242,7 @@ 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)) { + if (!invert3x3(a)) { output_error.write("\tERROR: metric tensor is singular at ({:d}, {:d})\n", i.x(), i.y()); return 1; @@ -1298,7 +1298,7 @@ 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)) { + if (!invert3x3(a)) { output_error.write("\tERROR: metric tensor is singular at ({:d}, {:d})\n", i.x(), i.y()); return 1; diff --git a/src/mesh/invert3x3.hxx b/src/mesh/invert3x3.hxx index dce208338d..84278e2e43 100644 --- a/src/mesh/invert3x3.hxx +++ b/src/mesh/invert3x3.hxx @@ -34,10 +34,10 @@ /// /// 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. +/// 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 -int invert3x3(Matrix& a, BoutReal small = 1.0e-15) { +bool invert3x3(Matrix& a, T small = 1.0e-15) { TRACE("invert3x3"); // Calculate the first co-factors @@ -51,9 +51,8 @@ int invert3x3(Matrix& a, BoutReal small = 1.0e-15) { if (std::abs(det) < std::abs(small)) { if (small >= 0) { throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small); - } else { - return 1; } + return false; } // Calculate the rest of the co-factors @@ -77,5 +76,5 @@ int invert3x3(Matrix& a, BoutReal small = 1.0e-15) { a(2, 1) = F * detinv; a(2, 2) = I * detinv; - return 0; + return true; } From a7114c89c930438a2debdde4c0e0f6c3f28eadbf Mon Sep 17 00:00:00 2001 From: ZedThree Date: Wed, 30 Oct 2024 14:21:46 +0000 Subject: [PATCH 03/10] Apply clang-format changes --- src/mesh/coordinates.cxx | 2 +- tests/unit/mesh/test_invert3x3.cxx | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index b57ad055da..3694b3d1c1 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -18,9 +18,9 @@ #include +#include "invert3x3.hxx" #include "parallel/fci.hxx" #include "parallel/shiftedmetricinterp.hxx" -#include "invert3x3.hxx" // use anonymous namespace so this utility function is not available outside this file namespace { diff --git a/tests/unit/mesh/test_invert3x3.cxx b/tests/unit/mesh/test_invert3x3.cxx index 02beeec644..77b08354cc 100644 --- a/tests/unit/mesh/test_invert3x3.cxx +++ b/tests/unit/mesh/test_invert3x3.cxx @@ -86,4 +86,3 @@ TEST(Invert3x3Test, BadCondition) { input(2, 2) = 1.0; EXPECT_NO_THROW(invert3x3(input, -1.0e-10)); } - From 1aa0e2387212f7b24eb61b54b03520ae344e1a0d Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 7 Nov 2024 16:15:28 +0000 Subject: [PATCH 04/10] Return `std::optional` from `invert3x3` Allows throwing more specific error in coordinates --- src/mesh/coordinates.cxx | 14 +++++----- src/mesh/invert3x3.hxx | 43 ++++++++++++++---------------- tests/unit/mesh/test_invert3x3.cxx | 28 +++++-------------- 3 files changed, 35 insertions(+), 50 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 3694b3d1c1..ce096cc8f9 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1242,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; } @@ -1298,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 index 84278e2e43..9e635d8150 100644 --- a/src/mesh/invert3x3.hxx +++ b/src/mesh/invert3x3.hxx @@ -29,42 +29,38 @@ #pragma once #include +#include /// 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 -bool invert3x3(Matrix& 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 invert3x3(Matrix& 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{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& a, T small = 1.0e-15) { a(2, 1) = F * detinv; a(2, 2) = I * detinv; - return true; + return std::nullopt; } +} // namespace bout diff --git a/tests/unit/mesh/test_invert3x3.cxx b/tests/unit/mesh/test_invert3x3.cxx index 77b08354cc..3bc4ae69d8 100644 --- a/tests/unit/mesh/test_invert3x3.cxx +++ b/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 input(3, 3); input = 0; - EXPECT_THROW(invert3x3(input), BoutException); + auto result = bout::invert3x3(input); + EXPECT_TRUE(result.has_value()); } 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); + 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()); } From d6aa81f2177a19c32074428c328cd5e238fb91bb Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 26 Nov 2024 11:10:28 +0100 Subject: [PATCH 05/10] simplify return statement --- src/mesh/invert3x3.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/invert3x3.hxx b/src/mesh/invert3x3.hxx index 9e635d8150..9c6b614168 100644 --- a/src/mesh/invert3x3.hxx +++ b/src/mesh/invert3x3.hxx @@ -48,7 +48,7 @@ inline std::optional invert3x3(Matrix& a) { 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}; + return det; } // Calculate the rest of the co-factors From e70fb7a2573be5b13c111cf68cee010430767f9a Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 26 Nov 2024 11:11:26 +0100 Subject: [PATCH 06/10] Use formatter for SpecificInd This works for 2D and 3D fields (and is also shorter code) --- src/mesh/coordinates.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index ce096cc8f9..e59e48042c 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1244,8 +1244,8 @@ int Coordinates::calcCovariant(const std::string& region) { 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()); + "\tERROR: metric tensor is singular at {}, determinant: {:d}\n", + i, det.value()); return 1; } @@ -1301,8 +1301,8 @@ int Coordinates::calcContravariant(const std::string& region) { 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()); + "\tERROR: metric tensor is singular at {}, determinant: {:d}\n", + i, det.value()); return 1; } From 12292e861f9a5edc3e563f52ff72e4a4260ae7d1 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 26 Nov 2024 10:35:39 +0000 Subject: [PATCH 07/10] Apply clang-format changes --- externalpackages/PVODE/include/pvode/band.h | 12 +++--------- externalpackages/PVODE/precon/band.h | 12 +++--------- externalpackages/PVODE/precon/pvbbdpre.cpp | 5 ++--- src/mesh/coordinates.cxx | 10 ++++------ 4 files changed, 12 insertions(+), 27 deletions(-) diff --git a/externalpackages/PVODE/include/pvode/band.h b/externalpackages/PVODE/include/pvode/band.h index 1fd04a2057..7ba083f91d 100644 --- a/externalpackages/PVODE/include/pvode/band.h +++ b/externalpackages/PVODE/include/pvode/band.h @@ -57,7 +57,6 @@ namespace pvode { - /****************************************************************** * * * Type: BandMat * @@ -118,7 +117,6 @@ namespace pvode { * * ******************************************************************/ - typedef struct bandmat_type { integer size; integer mu, ml, smu; @@ -128,7 +126,6 @@ typedef struct bandmat_type { /* BandMat accessor macros */ - /****************************************************************** * * * Macro : BAND_ELEM * @@ -141,8 +138,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_ELEM(A,i,j) ((A->data)[j][i-j+(A->smu)]) - +#define BAND_ELEM(A, i, j) ((A->data)[j][i - j + (A->smu)]) /****************************************************************** * * @@ -157,8 +153,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_COL(A,j) (((A->data)[j])+(A->smu)) - +#define BAND_COL(A, j) (((A->data)[j]) + (A->smu)) /****************************************************************** * * @@ -173,8 +168,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_COL_ELEM(col_j,i,j) (col_j[i-j]) - +#define BAND_COL_ELEM(col_j, i, j) (col_j[i - j]) /* Functions that use the BandMat representation for a band matrix */ diff --git a/externalpackages/PVODE/precon/band.h b/externalpackages/PVODE/precon/band.h index 1fd04a2057..7ba083f91d 100644 --- a/externalpackages/PVODE/precon/band.h +++ b/externalpackages/PVODE/precon/band.h @@ -57,7 +57,6 @@ namespace pvode { - /****************************************************************** * * * Type: BandMat * @@ -118,7 +117,6 @@ namespace pvode { * * ******************************************************************/ - typedef struct bandmat_type { integer size; integer mu, ml, smu; @@ -128,7 +126,6 @@ typedef struct bandmat_type { /* BandMat accessor macros */ - /****************************************************************** * * * Macro : BAND_ELEM * @@ -141,8 +138,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_ELEM(A,i,j) ((A->data)[j][i-j+(A->smu)]) - +#define BAND_ELEM(A, i, j) ((A->data)[j][i - j + (A->smu)]) /****************************************************************** * * @@ -157,8 +153,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_COL(A,j) (((A->data)[j])+(A->smu)) - +#define BAND_COL(A, j) (((A->data)[j]) + (A->smu)) /****************************************************************** * * @@ -173,8 +168,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define BAND_COL_ELEM(col_j,i,j) (col_j[i-j]) - +#define BAND_COL_ELEM(col_j, i, j) (col_j[i - j]) /* Functions that use the BandMat representation for a band matrix */ diff --git a/externalpackages/PVODE/precon/pvbbdpre.cpp b/externalpackages/PVODE/precon/pvbbdpre.cpp index 3a1181dcf1..93493ccfa8 100644 --- a/externalpackages/PVODE/precon/pvbbdpre.cpp +++ b/externalpackages/PVODE/precon/pvbbdpre.cpp @@ -364,14 +364,13 @@ static void PVBBDDQJac(integer Nlocal, integer mudq, integer mldq, /* Restore ytemp, then form and load difference quotients */ for (j=group-1; j < Nlocal; j+=width) { ytemp_data[j] = y_data[j]; - col_j = BAND_COL(J,j); + col_j = BAND_COL(J, j); inc = MAX(rely*ABS(y_data[j]), minInc/ewt_data[j]); inc_inv = ONE/inc; i1 = MAX(0, j-mukeep); i2 = MIN(j+mlkeep, Nlocal-1); for (i=i1; i <= i2; i++) - BAND_COL_ELEM(col_j,i,j) = - inc_inv * (gtemp_data[i] - gy_data[i]); + BAND_COL_ELEM(col_j, i, j) = inc_inv * (gtemp_data[i] - gy_data[i]); } } } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index e59e48042c..6355f246de 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1243,9 +1243,8 @@ int Coordinates::calcCovariant(const std::string& region) { a(0, 2) = a(2, 0) = g13[i]; if (const auto det = bout::invert3x3(a); det.has_value()) { - output_error.write( - "\tERROR: metric tensor is singular at {}, determinant: {:d}\n", - i, det.value()); + output_error.write("\tERROR: metric tensor is singular at {}, determinant: {:d}\n", + i, det.value()); return 1; } @@ -1300,9 +1299,8 @@ int Coordinates::calcContravariant(const std::string& region) { a(0, 2) = a(2, 0) = g_13[i]; if (const auto det = bout::invert3x3(a); det.has_value()) { - output_error.write( - "\tERROR: metric tensor is singular at {}, determinant: {:d}\n", - i, det.value()); + output_error.write("\tERROR: metric tensor is singular at {}, determinant: {:d}\n", + i, det.value()); return 1; } From 00a156d7e02b4a888e5fc0c335eb118599e0464c Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 26 Nov 2024 13:23:10 +0100 Subject: [PATCH 08/10] Add missing header to format SpecificInd --- src/mesh/coordinates.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 6355f246de..34c524d1e7 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -15,6 +15,7 @@ #include #include #include +#include #include From f6f78f6d257aa575008f3b3011d01716eb7ad5c7 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 26 Nov 2024 12:26:43 +0000 Subject: [PATCH 09/10] Apply clang-format changes --- externalpackages/PVODE/include/pvode/band.h | 2 +- externalpackages/PVODE/precon/band.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/externalpackages/PVODE/include/pvode/band.h b/externalpackages/PVODE/include/pvode/band.h index 8c881f2f3d..1c2d21f7ef 100644 --- a/externalpackages/PVODE/include/pvode/band.h +++ b/externalpackages/PVODE/include/pvode/band.h @@ -168,7 +168,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define PVODE_BAND_COL_ELEM(col_j,i,j) (col_j[i-j]) +#define PVODE_BAND_COL_ELEM(col_j, i, j) (col_j[i - j]) /* Functions that use the BandMat representation for a band matrix */ diff --git a/externalpackages/PVODE/precon/band.h b/externalpackages/PVODE/precon/band.h index ae5a16af49..0817e3cdc0 100644 --- a/externalpackages/PVODE/precon/band.h +++ b/externalpackages/PVODE/precon/band.h @@ -168,7 +168,7 @@ typedef struct bandmat_type { * * ******************************************************************/ -#define PVODE_BAND_COL_ELEM(col_j,i,j) (col_j[i-j]) +#define PVODE_BAND_COL_ELEM(col_j, i, j) (col_j[i - j]) /* Functions that use the BandMat representation for a band matrix */ From 5a60a5265a8c99b5ce5cfec00d22f5d80d7c6006 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 26 Nov 2024 14:10:43 +0100 Subject: [PATCH 10/10] Prefere const --- src/mesh/invert3x3.hxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh/invert3x3.hxx b/src/mesh/invert3x3.hxx index 9c6b614168..c011f55bf7 100644 --- a/src/mesh/invert3x3.hxx +++ b/src/mesh/invert3x3.hxx @@ -40,9 +40,9 @@ 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); + const BoutReal A = a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1); + const BoutReal B = a(1, 2) * a(2, 0) - a(1, 0) * a(2, 2); + const 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; @@ -52,15 +52,15 @@ inline std::optional invert3x3(Matrix& a) { } // 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); + const BoutReal D = a(0, 2) * a(2, 1) - a(0, 1) * a(2, 2); + const BoutReal E = a(0, 0) * a(2, 2) - a(0, 2) * a(2, 0); + const BoutReal F = a(0, 1) * a(2, 0) - a(0, 0) * a(2, 1); + const BoutReal G = a(0, 1) * a(1, 2) - a(0, 2) * a(1, 1); + const BoutReal H = a(0, 2) * a(1, 0) - a(0, 0) * a(1, 2); + const 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; + const BoutReal detinv = 1.0 / det; a(0, 0) = A * detinv; a(0, 1) = D * detinv;