Scale linear system of equations in plane_stress and skip test if incompatible system#3884
Scale linear system of equations in plane_stress and skip test if incompatible system#3884timothy-nunn merged 4 commits intomainfrom
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #3884 +/- ##
==========================================
- Coverage 46.04% 46.04% -0.01%
==========================================
Files 123 123
Lines 28962 28965 +3
==========================================
Hits 13336 13336
- Misses 15626 15629 +3 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
@chris-ashe checked on WSL and no unit tests are skipped |
eaf954e to
9a920d6
Compare
|
Correctly skips on CSD3 (where #3027 was first observed) |
e0a8619 to
1394fab
Compare
jonmaddock
left a comment
There was a problem hiding this comment.
I approved of the changes: just a couple of comments for consideration. Well done of finding the root cause of these well-above-floating point errors from a single calculation between different systems!
| # change the solution provided each element of a given row is scaled by the same scalar | ||
| # and the corresponding entry in bb is also scaled the same amount. | ||
| row_scale = np.max(np.abs(aa), axis=1) | ||
| aa /= np.repeat(row_scale[:, np.newaxis], aa.shape[0], axis=1) |
There was a problem hiding this comment.
Wouldn't it be clearer to broadcast row_scale here?
| # Here, we scale aa such that the largest element on a given row is 1.0. This does not | ||
| # change the solution provided each element of a given row is scaled by the same scalar | ||
| # and the corresponding entry in bb is also scaled the same amount. | ||
| row_scale = np.max(np.abs(aa), axis=1) |
There was a problem hiding this comment.
This works nicely: thank you for helping me understand ill-posed problems and condition number. I can now see that this could cause floating-point errors to lead to (well-) above floating point errors in the solution for
It seems strange to me that there is apparently no convenience function for performing this in numpy, for example: (I see there's a cond() function). I believe this falls into the class of regularisation techniques with various other methods, but this custom solution is certainly the simplest and works fine. I don't see any reason to do any different.
Sadly, #3856 was short lived and enabling
numpy>=2re-introduced the issues observed in #3027.I have been unable to find an acceptably fast solution that will ensure parity on different systems when solving$Ax=b$ in the $10^{12}$ ).
plane_stressfunction. The issue arises (I believe) because the matrixaais ill-conditioned (on the order ofThe first thing this PR introduces is a row-wise scaling of the$10^{3}$ ). However, one of the
aamatrix andbbvector such that the largest value on a row ofaawill be1.0. This reduces the condition ofaasignificantly (to be on the order ofplane_stresstests still failed on Mac, the other two passed.The second thing I have done is skip this test if the system is deemed incompatible (we actually already had a function for this that I wrote back in 2021). This means that the
plane_stresstests are skipped on incompatible systems (e.g. Mac, Windows, CSD3) and will run when we think it should pass (on Ubuntu VMs and the CI). Maybe the reviewer could check that these tests don't get skipped on WSL/their Ubuntu laptop.I have also added a warning to the installation page to ensure this issue is abundantly obvious.