Scale linear system of equations in plane_stress and skip test if incompatible system

documentation/installation/installation.md

@@ -14,6 +14,11 @@ PROCESS is a command-line Python program. As such, it should be supported by a w
     wsl --install
     ```
 
+!!! Warning "Differences between systems"
+    PROCESS is a complex numerical code that performs thousands of calculations and, on certain systems (e.g. Mac and some non-Ubuntu Linux systems), we have observed that floating-point rounding error can slightly alter the results of individual models--they are not wrong, just different. However, in a code like PROCESS, this error _could_ accumulate over many iterations of the optimiser to cause non-convergence/convergence to a very different design point. [Issue #3972](https://github.com/ukaea/PROCESS/issues/3972) shows that these differences and shows how these small errors can accumulate to slightly change a solution. 
+
+    In some sections of the code, these differences between systems are more pronounced (e.g. where operations occur on ill-conditioned matrices). Therefore, in our test suite, we skip a couple of tests if the system is not Linux or has an `ldd` version less than `2.31` because of this floating-point rounding difference: you can look for tests using the `skip_if_incompatible_system` fixture. While we encourage users of all systems to use PROCESS, we highly recommend you check your results on our officially support systems (listed above) and Python versions. 
+
 
 Start by downloading the PROCESS source code using `git`:
 

process/tf_coil.py

@@ -5066,12 +5066,25 @@ def plane_stress(nu, rad, ey, j, nlayers, n_radial_array):
     # out of numba, running the code, and then copying the result
     # back in. This means that the linear algebra solve is not compiled and runs
     # as if it were written in normal Python.
-    # This is necessary because numba compiles against the SciPy algebra library,
-    # not the Numpy one. There are differences in the Scipy solvers depending on
-    # operating system/architecture/Scipy version that cause tests to fail.
+    # This is done because numba compiles against the SciPy algebra library,
+    # not the Numpy one. We have observed some odd behaviour when using the SciPy library
+    # for this specific problem and so opt for using the numpy library instead.
     # https://github.com/ukaea/PROCESS/issues/3027
     # https://github.com/scipy/scipy/issues/23639
     with numba.objmode(cc="float64[:]"):
+        # These matrices can often end up being very ill-conditioned which can lead to numerical
+        # instability when solving below.
+        # Scaling the matrix can help reduce numerical instability by reducing the condition of matrix.
+        # 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.
+        # NOTE: this does not entirely solve the numerical instability and you can get above-floating point
+        # differences in the result of this function depending on system.
+        row_scale = np.max(np.abs(aa), axis=1)
+        # The transpose below ensures the scale is repeated along the row, not the column
+        aa /= np.broadcast_to(row_scale, aa.shape).T
+        bb /= row_scale
+
         cc = np.linalg.solve(aa, bb)
 
     #  Multiply c by (-1) (John Last, internal CCFE memorandum, 21/05/2013)

tests/conftest.py

@@ -92,29 +92,32 @@ def opt_params_only(request: SubRequest) -> bool:
     return request.config.getoption("--opt-params-only")
 
 
-@pytest.fixture(scope="session", autouse=True)
-def precondition(request):
-    """Check a user for an outdated system
-    Warn a user if their system is outdated and could have
-    regression test floating point issues.
+@pytest.fixture
+def skip_if_incompatible_system():
+    """Skip the test using this fixture if it is detcted that their system is incompatible
+    and may raise errors because of floating-point rounding error.
+    """
+    if not system_compatible():
+        pytest.skip(
+            "This test could fail on your system due to differences caused by "
+            "floating-point rounding error"
+        )
 
-    Exits the test suite if trying to overwrite tests
-    to stop inaccurate test assets being written.
 
-    e.g. "pytest --overwrite" returns True here.
-    :param request: request fixture to access CLI args
-    :type request: SubRequest
+@pytest.fixture(scope="session", autouse=True)
+def running_on_compatible_system_warning():
+    """Check for an outdated system
+
+    Warn the user if their system is outdated and could have test floating point issues.
     """
-    compatible = system_compatible()
-    if compatible:
+    if system_compatible():
         return
     warnings.warn(
         """
-        \u001b[33m\033[1mYou are running the PROCESS test suite on an outdated system.\033[0m
-        This can cause floating point rounding errors in regression tests.
+        \u001b[33m\033[1mYou are running the PROCESS test suite on an incompatible system.\033[0m
+        This can cause floating point rounding errors in tests.
 
-        Please see documentation for information on running PROCESS (and tests)
-        using a Docker/Singularity container.
+        Some unit tests may be skipped!
         """,
         UserWarning,
         stacklevel=2,

tests/integration/test_pfcoil_int.py

@@ -1662,7 +1662,7 @@ def test_mtrx(pfcoil: PFCoil):
     assert_array_almost_equal(bvec, bvec_exp)
 
 
-def test_solv(pfcoil: PFCoil):
+def test_solv(pfcoil: PFCoil, skip_if_incompatible_system):
     """Test solv() with simple arguments.
 
     Running baseline_2019 results in 2D array args with 740 elements: unfeasible

tests/unit/test_tfcoil.py

@@ -5128,7 +5128,7 @@ class PlaneStressParam(NamedTuple):
         ),
     ),
 )
-def test_plane_stress(planestressparam, monkeypatch):
+def test_plane_stress(planestressparam, skip_if_incompatible_system):
     """
     Automatically generated Regression Unit Test for plane_stress.