Fix incorrect application of EPS (Papercut/ issue #1084) (#3000)

* Fixed all of the places where EPS may be incorrectly applied

* Fixed all places where EPS is used in atol instead of rtol incorrectly.

* Moved the EPS_FREECAD constant to geometry.constants; Identified and fixed the bug that caused the Interpolator's geometry to be considered infeasible (accidentally set the threshold too tight when the atol=0 overrode the default tolerance.).

* Increased test coverage

* Removed optimize.bisect and found an analytical solution using arctan instead.

* ruff check

* Corrected tolerance value to D_TOLERANCE, and added test coverage.

* deleted debugging print statement

* Minor edit to comments in geometry.constants

* Finished papercut.

bluemira/builders/tf_coils.py

@@ -24,10 +24,10 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import bluemira_debug_flush
 from bluemira.base.parameter_frame import Parameter, ParameterFrame, make_parameter_frame
 from bluemira.display import plot_2d
+from bluemira.geometry.constants import EPS_FREECAD
 from bluemira.geometry.coordinates import Coordinates
 from bluemira.geometry.face import BluemiraFace
 from bluemira.geometry.optimisation import GeomOptimisationProblem, KeepOutZone
@@ -323,7 +323,7 @@ def set_wire(self, wire: BluemiraWire):
         points = wire.discretize(byedges=True, ndiscr=200)
         arg_x_max = np.argmax(points.x)
         x_max_point = points[:, arg_x_max]
-        self._alpha_0 = wire.parameter_at(x_max_point, tolerance=10 * EPS)
+        self._alpha_0 = wire.parameter_at(x_max_point, tolerance=EPS_FREECAD)
 
     def make_ripple_constraint(
         self, parameterisation, solver, TF_ripple_limit, rip_con_tol

bluemira/builders/tools.py

@@ -27,9 +27,9 @@
 import bluemira.base.components as bm_comp
 import bluemira.geometry as bm_geo
 from bluemira.base.components import PhysicalComponent
-from bluemira.base.constants import EPS
 from bluemira.base.error import BuilderError, ComponentError
 from bluemira.builders._varied_offset import varied_offset
+from bluemira.geometry.constants import D_TOLERANCE
 from bluemira.geometry.face import BluemiraFace
 from bluemira.geometry.plane import BluemiraPlane
 from bluemira.geometry.tools import (
@@ -338,12 +338,18 @@ def make_circular_xy_ring(r_inner: float, r_outer: float) -> BluemiraFace:
     """
     centre = (0, 0, 0)
     axis = (0, 0, 1)
-    if np.isclose(r_inner, r_outer, rtol=0, atol=2 * EPS):
-        raise BuilderError(f"Cannot make an annulus where r_inner = r_outer = {r_inner}")
-
+    if r_inner <= 0 or r_outer <= 0:
+        raise BuilderError(f"Radii must be nonnegative {r_inner=}, {r_outer=}.")
     if r_inner > r_outer:
         r_inner, r_outer = r_outer, r_inner
 
+    # Make sure that the annulus is thick enough even when they're in float32
+    # (FreeCAD stores numbers as float32 if I understand correctly)
+    if np.isclose(r_inner, r_outer, rtol=0, atol=D_TOLERANCE):
+        raise BuilderError(
+            f"Cannot make an annulus so thin that {r_outer - r_inner = }mm"
+        )
+
     inner = make_circle(r_inner, center=centre, axis=axis)
     outer = make_circle(r_outer, center=centre, axis=axis)
     return BluemiraFace([outer, inner])

bluemira/codes/_freecadapi.py

@@ -54,7 +54,7 @@
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.codes._freecadconfig import _freecad_save_config
 from bluemira.codes.error import FreeCADError, InvalidCADInputsError
-from bluemira.geometry.constants import MINIMUM_LENGTH
+from bluemira.geometry.constants import EPS_FREECAD, MINIMUM_LENGTH
 from bluemira.utilities.tools import ColourDescriptor
 
 if TYPE_CHECKING:
@@ -349,7 +349,7 @@ def interpolate_bspline(
     if len(pntslist) < 2:  # noqa: PLR2004
         _err = "interpolate_bspline: not enough points"
         raise InvalidCADInputsError(_err + "\n")
-    if np.allclose(pntslist[0], pntslist[-1], rtol=0, atol=EPS):
+    if np.allclose(pntslist[0], pntslist[-1], rtol=EPS, atol=0):
         if len(pntslist) > 2:  # noqa: PLR2004
             if not closed:
                 bluemira_warn("interpolate_bspline: equal endpoints forced Closed")
@@ -959,7 +959,7 @@ def wire_value_at(wire: apiWire, distance: float) -> np.ndarray:
 
 
 def wire_parameter_at(
-    wire: apiWire, vertex: Iterable[float], tolerance: float = EPS * 10
+    wire: apiWire, vertex: Iterable[float], tolerance: float = EPS_FREECAD
 ) -> float:
     """
     Get the parameter value at a vertex along a wire.
@@ -1903,7 +1903,7 @@ def point_inside_shape(point: Iterable[float], shape: apiShape) -> bool:
     Whether or not the point is inside the shape
     """
     vector = apiVector(*point)
-    return shape.isInside(vector, EPS, True)
+    return shape.isInside(vector, EPS_FREECAD, True)
 
 
 # ======================================================================================
@@ -2026,7 +2026,9 @@ def _make_shapes_coaxis(shapes):
 # ======================================================================================
 
 
-def fix_wire(wire: apiWire, precision: float = EPS, min_length: float = MINIMUM_LENGTH):
+def fix_wire(
+    wire: apiWire, precision: float = EPS_FREECAD, min_length: float = MINIMUM_LENGTH
+):
     """
     Fix a wire by removing any small edges and joining the remaining edges.
 

bluemira/codes/_nlopt_api.py

@@ -14,15 +14,14 @@
 import nlopt
 import numpy as np
 
+from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.utilities.error import (
     ExternalOptError,
     OptUtilitiesError,
     OptVariablesError,
 )
 
-EPS = np.finfo(np.float64).eps
-
 NLOPT_ALG_MAPPING = {
     "SLSQP": nlopt.LD_SLSQP,
     "COBYLA": nlopt.LN_COBYLA,
@@ -86,8 +85,8 @@ def process_NLOPT_conditions(  # noqa: N802
         elif v > 0:
             if k in ["ftol_abs", "ftol_res", "xtol_abs", "xtol_res"] and v < EPS:
                 bluemira_warn(
-                    "You are setting an optimisation termination condition to below"
-                    " machine precision. Don't.."
+                    "You are setting an optimisation termination condition to a value"
+                    " too small given your machine's precision. Don't.."
                 )
 
             conditions[k] = v

bluemira/equilibria/positioner.py

@@ -923,6 +923,6 @@ def check_loop_feasibility(coords: Coordinates):
 
         """
         if not np.allclose(
-            ConvexHull(coords.xz.T).volume, get_area_2d(coords.x, coords.z), atol=EPS
+            ConvexHull(coords.xz.T).volume, get_area_2d(coords.x, coords.z), rtol=EPS
         ):
             raise GeometryError("Region must be a Convex Hull")

bluemira/geometry/constants.py

@@ -8,9 +8,16 @@
 Constants for the geometry module
 """
 
+import numpy as np
+
 # Absolute tolerance for equality in distances
 D_TOLERANCE = 1e-5  # [m]
 
+# FreeCAD's default tolerance is 1E-7, so correspondingly, we need a larger epsilon.
+# The following value happens to equal 1.1920929E-7, which satisifies the requirement.
+# We think that FreeCAD stores number as float 32 but we're not sure.
+EPS_FREECAD = np.finfo(np.float32).eps
+
 # Minimum length of a wire or sub-wire (edge)
 MINIMUM_LENGTH = 1e-5  # [m]
 

bluemira/geometry/coordinates.py

@@ -1251,7 +1251,7 @@ def closed(self) -> bool:
         Whether or not this is a closed set of Coordinates.
         """
         if len(self) > 2 and np.allclose(  # noqa: PLR2004
-            self[:, 0], self[:, -1], rtol=0, atol=EPS
+            self[:, 0], self[:, -1], rtol=EPS, atol=0
         ):
             return True
         return False
@@ -1428,7 +1428,7 @@ def __eq__(self, other: Coordinates) -> bool:
         counted as identical.
         """
         if isinstance(other, self.__class__):
-            return np.all(np.allclose(self._array, other._array, rtol=0, atol=EPS))
+            return np.all(np.allclose(self._array, other._array, rtol=EPS, atol=0))
         return False
 
     def __hash__(self):

bluemira/geometry/parameterisations.py

@@ -1699,6 +1699,7 @@ def _make_tapered_inner_leg(
             )
         theta = 2 * np.arctan2(x_mid - x_in, z_taper)
         r_taper = z_taper / np.sin(theta) - r_min
+
         if r_taper < r_min or theta >= (np.pi / 2):
             raise GeometryParameterisationError(
                 f"Cannot achieve radius of curvature <= {r_min=}"

bluemira/geometry/wire.py

@@ -13,10 +13,10 @@
 from typing import TYPE_CHECKING, Iterable, List, Optional, Tuple, Union
 
 import bluemira.codes._freecadapi as cadapi
-from bluemira.base.constants import EPS
 from bluemira.base.look_and_feel import LOGGER, bluemira_warn
 from bluemira.codes.error import FreeCADError
 from bluemira.geometry.base import BluemiraGeo, _Orientation
+from bluemira.geometry.constants import EPS_FREECAD
 from bluemira.geometry.coordinates import Coordinates
 from bluemira.geometry.error import (
     GeometryError,
@@ -191,7 +191,7 @@ def value_at(
         return cadapi.wire_value_at(self.shape, distance)
 
     def parameter_at(
-        self, vertex: Iterable[float], tolerance: float = EPS * 10
+        self, vertex: Iterable[float], tolerance: float = EPS_FREECAD
     ) -> float:
         """
         Get the parameter value at a vertex along a wire.

bluemira/magnetostatics/polyhedral_prism.py

@@ -28,7 +28,7 @@
 import numba as nb
 import numpy as np
 
-from bluemira.base.constants import MU_0_4PI
+from bluemira.base.constants import EPS, MU_0_4PI
 from bluemira.geometry.coordinates import Coordinates, get_area_2d
 from bluemira.magnetostatics.baseclass import (
     PolyhedralCrossSectionCurrentSource,
@@ -41,8 +41,6 @@
 # NOTE: Polyhedral kernels are not intended to be user-facing, but
 # it's useful for testing.
 
-ZERO_DIV_GUARD_EPS = 1e-14
-
 
 class PolyhedralKernel(abc.ABC):
     """
@@ -130,7 +128,7 @@ def _vector_norm_eps(r: np.ndarray) -> float:
     \t:math:`\\lvert \\mathbf{r}\\rvert = \\sqrt{\\lvert \\mathbf{r}\\rvert^2+\\epsilon^2}`
     """  # noqa: W505 E501
     r_norm = np.linalg.norm(r)
-    return np.sqrt(r_norm**2 + ZERO_DIV_GUARD_EPS)
+    return np.sqrt(r_norm**2 + EPS)  # guard against division by 0
 
 
 @nb.jit(nopython=True, cache=True)
@@ -181,7 +179,7 @@ def _omega_t(r: np.ndarray, r1: np.ndarray, r2: np.ndarray, r3: np.ndarray) -> f
     )
     a = np.dot(r1_r, np.cross(r2_r, r3_r))
     # Not sure this is an issue...
-    # if abs(a) < ZERO_GUARD_EPS and (-ZERO_GUARD_EPS < d < 0):
+    # if abs(a) < EPS and (-EPS < d < 0): # guard against division by 0
     #     return 0 # and not pi as per IEEE
     return 2 * np.arctan2(a, d)
 

bluemira/optimisation/_nlopt/conditions.py

@@ -40,8 +40,8 @@ def _validate(self) -> None:
         ]:
             if condition and condition < EPS:
                 bluemira_warn(
-                    "optimisation: Setting stopping condition to less than machine "
-                    "precision. This condition may never be met."
+                    "optimisation: Setting stopping condition is too small given this "
+                    "machine's precision. This condition may never be met."
                 )
         if self._no_stopping_condition_set():
             raise OptimisationConditionsError(

bluemira/radiation_transport/advective_transport.py

@@ -199,7 +199,7 @@ def _make_flux_surfaces_ob(self):
         self.flux_surfaces_ob_hfs = []
 
         x = self.x_sep_omp + self.dx_mp
-        while x < x_out_omp - EPS:
+        while x < x_out_omp:
             lfs, hfs = self._make_flux_surfaces(x, self._o_point.z)
             self.flux_surfaces_ob_lfs.append(lfs)
             self.flux_surfaces_ob_hfs.append(hfs)
@@ -214,7 +214,7 @@ def _make_flux_surfaces_ib(self):
         self.flux_surfaces_ib_lfs = []
         self.flux_surfaces_ib_hfs = []
         x = self.x_sep_imp - self.dx_mp
-        while x > x_out_imp + EPS:
+        while x > x_out_imp:
             lfs, hfs = self._make_flux_surfaces(x, self._o_point.z)
             self.flux_surfaces_ib_lfs.append(lfs)
             self.flux_surfaces_ib_hfs.append(hfs)

bluemira/utilities/positioning.py

@@ -20,8 +20,7 @@
 import numpy as np
 from scipy.spatial import ConvexHull
 
-from bluemira.base.constants import EPS
-from bluemira.geometry.constants import VERY_BIG
+from bluemira.geometry.constants import D_TOLERANCE, VERY_BIG
 from bluemira.geometry.plane import BluemiraPlane
 from bluemira.geometry.tools import slice_shape
 from bluemira.utilities.error import PositionerError
@@ -43,12 +42,12 @@ class XZGeometryInterpolator(abc.ABC):
     def __init__(self, geometry: BluemiraWire):
         self.geometry = geometry
 
-    def _get_xz_coordinates(self):
+    def _get_xz_coordinates(self, num_pts):
         """
         Get discretised x-z coordinates of the geometry.
         """
         coordinates = self.geometry.discretize(
-            byedges=True, dl=self.geometry.length / 1000
+            byedges=True, dl=self.geometry.length / num_pts
         )
         coordinates.set_ccw([0, 1, 0])
         return coordinates.xz
@@ -173,14 +172,17 @@ def _check_geometry_feasibility(self, geometry: BluemiraWire):
         if not self.geometry.is_closed:
             raise PositionerError("RegionInterpolator can only handle closed wires.")
 
-        xz_coordinates = self._get_xz_coordinates()
+        xz_coordinates = self._get_xz_coordinates(10000)
         hull = ConvexHull(xz_coordinates.T)
         # Yes, the "area" of a 2-D scipy ConvexHull is its perimeter...
-        if not np.allclose(hull.area, geometry.length, atol=EPS):
+        if not np.allclose(hull.area, geometry.length, rtol=1e-5, atol=D_TOLERANCE):
             raise PositionerError(
-                "RegionInterpolator can only handle convex geometries. Perimeter"
-                " difference between convex hull and geometry:"
-                f" {hull.area - geometry.length}"
+                "RegionInterpolator can only handle simple convex geometries. "
+                "Perimeter difference between convex hull and geometry:"
+                f"{hull.area} - {geometry.length} ="
+                f"{hull.area - geometry.length}\n"
+                "This suggests that the shape is"
+                "too complex to be discretized accurately."
             )
 
     def to_xz(

tests/builders/test_tools.py

@@ -294,9 +294,16 @@ def test_annulus_area_reversed_radii(self, r_in, r_out):
         face = make_circular_xy_ring(r_in, r_out)
         np.testing.assert_almost_equal(face.area, np.pi * (r_in**2 - r_out**2))
 
-    def test_raises_error_on_equal_radii(self):
+    @pytest.mark.parametrize(
+        ("r_in", "r_out"),
+        [
+            (1, 1),  # ring too thin
+            (-0.01, 1),  # negative radius
+        ],
+    )
+    def test_raises_error_on_equal_radii(self, r_in, r_out):
         with pytest.raises(BuilderError):
-            make_circular_xy_ring(1, 1)
+            make_circular_xy_ring(r_in, r_out)
 
 
 class TestBuildSectioned: