Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Move invert3x3 out of general purpose utils.hxx header #3018

Open
wants to merge 4 commits into
base: next
Choose a base branch
from

Conversation

ZedThree
Copy link
Member

Only used in Coordinates, so make private implementation detail

Also return bool instead of int

Copy link
Contributor

@dschwoerer dschwoerer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The negative small value is only used in tests.

As this is private, we can remove this usage.

Then the checks on the caller can be removed, making the code cleaner. Unfortunately I cannot suggest this, so I can push a commit if you agree, @ZedThree

Comment on lines 81 to 87

// 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));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// 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));

Comment on lines 37 to 38
/// 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.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
/// 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.

Comment on lines 52 to 55
if (small >= 0) {
throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small);
}
return false;
Copy link
Contributor

@dschwoerer dschwoerer Oct 30, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (small >= 0) {
throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small);
}
return false;
throw BoutException("Determinant of metric component {:e} < {:e} --> Poorly conditioned", std::abs(det), small);

/// 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) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
bool invert3x3(Matrix<T>& a, T small = 1.0e-15) {
void invert3x3(Matrix<T>& a, T small = 1.0e-15) {

Comment on lines 78 to 79

return true;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return true;

@ZedThree
Copy link
Member Author

Yes, I was tempted to do the same thing, although I realised the error message in Coordinates is actually more specific. We could rewrite the interface even further, we calculate all 9 elements and we only need 6!

So maybe we just remove the throw from invert3x3 but still return bool and do the error handling in Coordinates?

@dschwoerer
Copy link
Contributor

Right now the error is thrown from invert3x3 - so the more specific error is actually not used.

We can throw a more specific error, as it is only used to invert the metric tensor.

@ZedThree
Copy link
Member Author

Yes, exactly

@dschwoerer
Copy link
Contributor

We could rewrite the interface even further, we calculate all 9 elements and we only need 6!

Introduce a new SymmetricMatrix<T> type :-)

Allows throwing more specific error in coordinates
Copy link
Contributor

@dschwoerer dschwoerer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, thanks @ZedThree

Comment on lines +1247 to +1248
"\tERROR: metric tensor is singular at ({:d}, {:d}), determinant: {:d}\n",
i.x(), i.y(), det.value());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could be 3D compatible:

Suggested change
"\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());

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice, good catch, thanks!

Comment on lines +1304 to +1305
"\tERROR: metric tensor is singular at ({:d}, {:d}), determinant: {:d}\n",
i.x(), i.y(), det.value());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same:

Suggested change
"\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());

throw BoutException("Determinant of matrix < {:e} --> Poorly conditioned", small);
}
return false;
return std::optional<BoutReal>{det};
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the explicit constructor really needed? Could this not be:

Suggested change
return std::optional<BoutReal>{det};
return det;

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure, will check

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants