change[autograd]: set minimum discretization wavelength when computing autograd derivatives of cylinders and make adaptive spacing relative to minimum free space wavelength

Author: Gregory Roberts

CHANGELOG.md

@@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - Validate mode solver object for large number of grid points on the modal plane.
+- Adaptive minimum spacing for `PolySlab` integration is now wavelength relative and a minimum discretization is set for computing gradients for cylinders.
 
 ### Fixed
 - Fixed missing amplitude factor and handling of negative normal direction case when making adjoint sources from `DiffractionMonitor`.

tests/test_components/test_autograd.py

@@ -19,7 +19,11 @@
 
 import tidy3d as td
 import tidy3d.web as web
-from tidy3d.components.autograd.constants import MAX_NUM_TRACED_STRUCTURES
+from tidy3d.components.autograd.constants import (
+    MAX_NUM_TRACED_STRUCTURES,
+    MIN_WVL_FRACTION_CYLINDER_DISCRETIZE,
+    MINIMUM_SPACING_FRACTION,
+)
 from tidy3d.components.autograd.derivative_utils import DerivativeInfo
 from tidy3d.components.autograd.utils import is_tidy_box
 from tidy3d.components.data.data_array import DataArray
@@ -1638,6 +1642,70 @@ def f(eps_inf, poles):
             assert np.isclose(grads_computed[field_path], grad_poles[i][j])
 
 
+@pytest.mark.parametrize("eps_real", [1e6, -1e8])
+def test_adaptive_spacing(eps_real):
+    freq = 5e9
+
+    info = DerivativeInfo(
+        paths={},
+        E_der_map={},
+        D_der_map={},
+        E_fwd={},
+        D_fwd={},
+        E_adj={},
+        D_adj={},
+        eps_data={},
+        eps_in=eps_real,
+        eps_out=1.0,
+        frequencies=[freq],
+        bounds=((-1, -1, -1), (1, 1, 1)),
+        eps_no_structure={},
+        eps_inf_structure={},
+        bounds_intersect=((-1, -1, -1), (1, 1, 1)),
+    )
+
+    with AssertLogLevel("WARNING", contains_str="Based on the material, the adaptive spacing"):
+        expected_vjp_spacing = info.wavelength_min * MINIMUM_SPACING_FRACTION
+        vjp_spacing = info.adaptive_vjp_spacing()
+
+        assert np.isclose(expected_vjp_spacing, vjp_spacing), "Unexpected adaptive vjp spacing!"
+
+
+@pytest.mark.parametrize("eps_real", [1e6, -1e8])
+def test_cylinder_discretization(eps_real):
+    freq = 5e9
+
+    info = DerivativeInfo(
+        paths={},
+        E_der_map={},
+        D_der_map={},
+        E_fwd={},
+        D_fwd={},
+        E_adj={},
+        D_adj={},
+        eps_data={},
+        eps_in=eps_real,
+        eps_out=1.0,
+        frequencies=[freq],
+        bounds=((-1, -1, -1), (1, 1, 1)),
+        eps_no_structure={},
+        eps_inf_structure={},
+        bounds_intersect=((-1, -1, -1), (1, 1, 1)),
+    )
+
+    with AssertLogLevel(
+        "WARNING", contains_str="The minimum wavelength inside the cylinder material"
+    ):
+        cylinder = td.Cylinder(axis=2, length=info.wavelength_min, radius=2 * info.wavelength_min)
+
+        expected_wvl_mat = info.wavelength_min * MIN_WVL_FRACTION_CYLINDER_DISCRETIZE
+        wvl_mat = cylinder._discretization_wavelength(derivative_info=info)
+
+        assert np.isclose(expected_wvl_mat, wvl_mat), (
+            "Unexpected wavelength for discretizing cylinder!"
+        )
+
+
 def test_custom_pole_residue(monkeypatch):
     """Test that computed pole residue derivatives match."""
 

tidy3d/components/autograd/constants.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 
+# minimum fraction of minimum free space wavelength for discretizing cylinder in autograd derivative
+MIN_WVL_FRACTION_CYLINDER_DISCRETIZE = 5e-2
 # default number of points per wvl in material for discretizing cylinder in autograd derivative
 PTS_PER_WVL_MAT_CYLINDER_DISCRETIZE = 10
 
@@ -19,7 +21,8 @@
 AUTOGRAD_MONITOR_INTERVAL_SPACE_CUSTOM = (1, 1, 1)
 
 DEFAULT_WAVELENGTH_FRACTION = 0.1
-MINIMUM_SPACING = 1e-2
+# minimum fraction of minimum free space wavelength to be used when computing adaptive spacing
+MINIMUM_SPACING_FRACTION = 1e-2
 
 EDGE_CLIP_TOLERANCE = 1e-9
 

tidy3d/components/autograd/derivative_utils.py

@@ -11,12 +11,13 @@
 from tidy3d.components.data.data_array import FreqDataArray, ScalarFieldDataArray
 from tidy3d.components.types import ArrayLike, Bound, tidycomplex
 from tidy3d.constants import C_0, LARGE_NUMBER
+from tidy3d.log import log
 
 from .constants import (
     DEFAULT_WAVELENGTH_FRACTION,
     GRADIENT_DTYPE_COMPLEX,
     GRADIENT_DTYPE_FLOAT,
-    MINIMUM_SPACING,
+    MINIMUM_SPACING_FRACTION,
 )
 from .types import PathType
 from .utils import get_static
@@ -430,7 +431,7 @@ def _project_in_basis(
     def adaptive_vjp_spacing(
         self,
         wl_fraction: float = DEFAULT_WAVELENGTH_FRACTION,
-        min_allowed_spacing: float = MINIMUM_SPACING,
+        min_allowed_spacing_fraction: float = MINIMUM_SPACING_FRACTION,
     ) -> float:
         """Compute adaptive spacing for finite-difference gradient evaluation.
 
@@ -441,8 +442,9 @@ def adaptive_vjp_spacing(
         ----------
         wl_fraction : float = 0.1
             Fraction of wavelength/skin depth to use as spacing.
-        min_allowed_spacing : float = 1e-2
-            Minimum allowed spacing to prevent numerical issues.
+        min_allowed_spacing_fraction : float = 1e-2
+            Minimum allowed spacing fraction of free space wavelength to
+            prevent numerical issues.
 
         Returns
         -------
@@ -471,7 +473,17 @@ def adaptive_vjp_spacing(
             delta_min = C_0 / (omega * np.sqrt(np.abs(eps_neg).max()))
             dx_candidates.append(wl_fraction * delta_min)
 
-        return max(min(dx_candidates), min_allowed_spacing)
+        computed_spacing = min(dx_candidates)
+        min_allowed_spacing = self.wavelength_min * min_allowed_spacing_fraction
+
+        if computed_spacing < min_allowed_spacing:
+            log.warning(
+                f"Based on the material, the adaptive spacing for integrating the polyslab surface "
+                f"would be {computed_spacing:.3e} μm. The spacing has been clipped to {min_allowed_spacing:.3e} μm "
+                f"to prevent a performance degradation.",
+                log_once=True,
+            )
+        return max(computed_spacing, min_allowed_spacing)
 
     @property
     def wavelength_min(self) -> float:

tidy3d/components/geometry/primitives.py

@@ -11,12 +11,16 @@
 import shapely
 
 from tidy3d.components.autograd import AutogradFieldMap, TracedSize1D
-from tidy3d.components.autograd.constants import PTS_PER_WVL_MAT_CYLINDER_DISCRETIZE
+from tidy3d.components.autograd.constants import (
+    MIN_WVL_FRACTION_CYLINDER_DISCRETIZE,
+    PTS_PER_WVL_MAT_CYLINDER_DISCRETIZE,
+)
 from tidy3d.components.autograd.derivative_utils import DerivativeInfo
 from tidy3d.components.base import cached_property, skip_if_fields_missing
 from tidy3d.components.types import Axis, Bound, Coordinate, MatrixReal4x4, Shapely
 from tidy3d.constants import LARGE_NUMBER, MICROMETER
 from tidy3d.exceptions import SetupError, ValidationError
+from tidy3d.log import log
 from tidy3d.packaging import verify_packages_import
 
 from . import base
@@ -278,12 +282,29 @@ def _points_unit_circle(
         ys = np.sin(angles)
         return np.stack((xs, ys), axis=0)
 
+    def _discretization_wavelength(self, derivative_info: DerivativeInfo) -> float:
+        """Choose a reference wavelength for discretizing the cylinder into a `PolySlab`."""
+        wvl0_min = derivative_info.wavelength_min
+        wvl_mat = wvl0_min / np.max([1.0, np.max(np.sqrt(abs(derivative_info.eps_in)))])
+
+        min_wvl_mat = MIN_WVL_FRACTION_CYLINDER_DISCRETIZE * wvl0_min
+        if wvl_mat < min_wvl_mat:
+            log.warning(
+                f"The minimum wavelength inside the cylinder material is {wvl_mat:.3e} μm, which would "
+                f"create a large number of discretization points for computing the gradient. "
+                f"To prevent performance degradation, the discretization wavelength has "
+                f"been clipped to {min_wvl_mat:.3e} μm.",
+                log_once=True,
+            )
+        wvl_mat = max(wvl_mat, min_wvl_mat)
+
+        return wvl_mat
+
     def _compute_derivatives(self, derivative_info: DerivativeInfo) -> AutogradFieldMap:
         """Compute the adjoint derivatives for this object."""
 
         # compute circumference discretization
-        wvl0_min = derivative_info.wavelength_min
-        wvl_mat = wvl0_min / np.max([1.0, np.max(np.sqrt(abs(derivative_info.eps_in)))])
+        wvl_mat = self._discretization_wavelength(derivative_info=derivative_info)
 
         circumference = 2 * np.pi * self.radius
         wvls_in_circumference = circumference / wvl_mat