diff --git a/bluemira/builders/tf_coils.py b/bluemira/builders/tf_coils.py index 855d961842..04dde5a277 100644 --- a/bluemira/builders/tf_coils.py +++ b/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 diff --git a/bluemira/builders/tools.py b/bluemira/builders/tools.py index 14f82fe605..e19a887340 100644 --- a/bluemira/builders/tools.py +++ b/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]) diff --git a/bluemira/codes/_freecadapi.py b/bluemira/codes/_freecadapi.py index 0efe810eb4..fc08af31b3 100644 --- a/bluemira/codes/_freecadapi.py +++ b/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. diff --git a/bluemira/codes/_nlopt_api.py b/bluemira/codes/_nlopt_api.py index ac2835b4e8..f701f9ef12 100644 --- a/bluemira/codes/_nlopt_api.py +++ b/bluemira/codes/_nlopt_api.py @@ -14,6 +14,7 @@ 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, @@ -21,8 +22,6 @@ 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 diff --git a/bluemira/equilibria/positioner.py b/bluemira/equilibria/positioner.py index af14642376..e07e86f783 100644 --- a/bluemira/equilibria/positioner.py +++ b/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") diff --git a/bluemira/geometry/constants.py b/bluemira/geometry/constants.py index 046d8e690e..1c29eb6b26 100644 --- a/bluemira/geometry/constants.py +++ b/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] diff --git a/bluemira/geometry/coordinates.py b/bluemira/geometry/coordinates.py index c53605815d..f61c798966 100644 --- a/bluemira/geometry/coordinates.py +++ b/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): diff --git a/bluemira/geometry/parameterisations.py b/bluemira/geometry/parameterisations.py index e1ef6e9c03..df69dbb60f 100644 --- a/bluemira/geometry/parameterisations.py +++ b/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=}" diff --git a/bluemira/geometry/wire.py b/bluemira/geometry/wire.py index c3bcfb74e5..fe7940514e 100644 --- a/bluemira/geometry/wire.py +++ b/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. diff --git a/bluemira/magnetostatics/polyhedral_prism.py b/bluemira/magnetostatics/polyhedral_prism.py index 52c73ba8a2..e035be8d92 100644 --- a/bluemira/magnetostatics/polyhedral_prism.py +++ b/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) diff --git a/bluemira/optimisation/_nlopt/conditions.py b/bluemira/optimisation/_nlopt/conditions.py index 247f2faae5..d60952e01a 100644 --- a/bluemira/optimisation/_nlopt/conditions.py +++ b/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( diff --git a/bluemira/radiation_transport/advective_transport.py b/bluemira/radiation_transport/advective_transport.py index 53b2d20418..06b40b5e4e 100644 --- a/bluemira/radiation_transport/advective_transport.py +++ b/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) diff --git a/bluemira/utilities/positioning.py b/bluemira/utilities/positioning.py index 0b4856f396..87b144c02e 100644 --- a/bluemira/utilities/positioning.py +++ b/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( diff --git a/tests/builders/test_tools.py b/tests/builders/test_tools.py index 60689f2c81..b448e075fc 100644 --- a/tests/builders/test_tools.py +++ b/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: